ComputeStress Function – Error ?
OpenLB – Open Source Lattice Boltzmann Code › Forums › on OpenLB › General Topics › ComputeStress Function – Error ?
 This topic has 5 replies, 2 voices, and was last updated 6 years, 2 months ago by mathias.

AuthorPosts

June 5, 2016 at 6:42 pm #1834Alejandro_ClaroMember
Hello everyone,rnrnI am interested to simulate with the LES (Large Eddy Simulation) model. I want to use the Smagorinsky model with the BGK collision operator. The class for this option is olb::SmagorinskyBGKdynamics< T, Lattice >. In order to implement the collision step it is used the function momenta.computeAllMomenta(cell, rho, u, pi) which for the fluid nodes I supposed it is used olb::BulkMomenta< T, Lattice > Struct Template Reference. In this struct the computeAllMomenta is computed with this function (defined at line 271 of file lbHelpers.h.):rn
Code:rntemplate<typename T, class Descriptor>rnstatic void olb::lbDynamicsHelpers< T, Descriptor >::computeStress (rn CellBase< T, Descriptor > const & cell,rn T rho,rn const T u[Descriptor::d],rn T pi[util::TensorVal< Descriptor >::n] rn ) rnrnrn
Code:rnstatic void computeStress(rn CellBase<T,Descriptor> const& cell, rn T rho, rn const T u[Descriptor::d],rn T pi[util::TensorVal<Descriptor>::n] ) rn{rn int iPi = 0;rn for (int iAlpha=0; iAlpha < Descriptor::d; ++iAlpha) {rn for (int iBeta=iAlpha; iBeta < Descriptor::d; ++iBeta) {rn pi[iPi] = T();rn for (int iPop=0; iPop < Descriptor::q; ++iPop) {rn pi[iPi] += Descriptor::c[iPop][iAlpha]*rn Descriptor::c[iPop][iBeta] * cell[iPop];rn }rn // stripe off equilibrium contributionrn pi[iPi] = rho*u[iAlpha]*u[iBeta];rn if (iAlpha==iBeta) {rn pi[iPi] = 1./Descriptor::invCs2*(rho(T)1);rn }rn ++iPi;rn }rn }rn}rnrnrnThe pi vector is the stress tensor for the nonequilibrium function. The stripe off equilibrium contribution part of the code seems to be incorrect. When iAlpha==iBeta the code should be:rn
Code:rnif (iAlpha==iBeta) {rn //pi[iPi] = 1./Descriptor::invCs2*(rho(T)1);rn pi[iPi] = 1./Descriptor::invCs2*(rho);rn }rnrnrnCould someone confirms me that this “”correction”” is correct or not?rnrnBest Regards,rnrnAlejandro
June 6, 2016 at 5:50 am #2349mathiasKeymasterHi Alejandro,rnrnwhy? It is the pressure that is taken. We checked with an analytical solution of NS equation. I am quite sure that this is correct.rnrnBestrnMathias
June 6, 2016 at 6:50 am #2350Alejandro_ClaroMemberHi Mathias,rnrnThe ComputeStress seems to give the secondorder moment of the nonequilibrium distribution fi(1). The secondorder moment of the distribution function is Sigma(Ci_A*Ci_B*[fi – fi(0)]). Ci is the discrete velocity; fi and fi(0) are the actual and the the equilibrium distribution functions in the lattice, respectively.rnrnSupposing that BGK collision model is used with a rectangular method for the collision operator. Then, the kinematic shear viscosity is
Code:nu = (tau – 0.5)*pow(Cs,2)Where Cs is the velocity of the sound in the lattice system and tau is the relaxation time. The pressure is defined by
Code:P = pow(Cs,2)*rhoThe secondorder moment of the equilibrium distribution fi(0) is rn
Code:rnSecond_Order_Moment_fi0 = rho*u_A*u_B + P* KroneckerDelta(A,B);rnSecond_Order_Moment_fi0 = rho*u_A*u_B + (pow(Cs,2)*rho)* KroneckerDelta(A,B);rnrnrnSo, the pressure is already is taken by the second rightside term of the previous expression, or I am wrong ? Thus, I do not understand why the extra term 1.*pow(Cs,2) is added to the previous expression? rnrnIn the user guide is indicated that the Smagorinsky model is based in the paper Nathen et al. (2013) An extension of the Lattice Boltzmann Method for simulating turbulent flows around rotating geometries of arbitrary shape. Unfortunately, I do not have access rights to read this document. The pressure term in this paper is obtained differently from the previous “”explanation”” ?rnrnBest Regards,rnrnAlejandro
June 6, 2016 at 5:59 pm #2351mathiasKeymasterIn NS equation, it is the pressure gradient that is to be found and not the pressure. Therefore, pressure is definded up to a constant, so that “”rho1″” or “”rho”” should result in the same solution of NS equation. Using “”rho1″” brings the pressure constant to zero, since usualy rho is set to be 1 + h^2. Do you see any differences in the result? Best Mathias
June 8, 2016 at 7:23 am #2352Alejandro_ClaroMemberHi Mathias,rnrnThank you for your response. When you said that “”rho”” is usually set to be 1 + h², this h is the grid spacing?rnrnI am in the process to implement a case of turbulent flow. Thus, as soon as I get some results I will share with you.rnrnBest Regards, rnAlejandro
June 8, 2016 at 7:53 am #2362mathiasKeymasterHi Alejandro,rnrnyes, h=latticeL (grid spacing) and h=latticeU (Mach no.). This is diffusive scaling. I am happy to see your results.rnrnBestrnMathias

AuthorPosts
 You must be logged in to reply to this topic.