ComputeStress Function – Error ?

[Alejandro_Claro]

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 ) rn

rnrn

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}rn

rnrnThe pi vector is the stress tensor for the non-equilibrium 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 }rn

rnrnCould someone confirms me that this “”correction”” is correct or not?rnrnBest Regards,rnrnAlejandro

[mathias]

Hi 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

[Alejandro_Claro]

Hi Mathias,rnrnThe ComputeStress seems to give the second-order moment of the non-equilibrium distribution fi(1). The second-order 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)*rho

The second-order 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);rn

rnrnSo, the pressure is already is taken by the second right-side 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

[mathias]

In 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 “”rho-1″” or “”rho”” should result in the same solution of NS equation. Using “”rho-1″” 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

[Alejandro_Claro]

Hi 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

[mathias]

Hi Alejandro,rnrnyes, h=latticeL (grid spacing) and h=latticeU (Mach no.). This is diffusive scaling. I am happy to see your results.rnrnBestrnMathias

[Author]

Posts