Skip to content

rho vs. rho +/- 1 (esp. for computePressureField())

OpenLB – Open Source Lattice Boltzmann Code Forums on OpenLB General Topics rho vs. rho +/- 1 (esp. for computePressureField())

Viewing 8 posts - 1 through 8 (of 8 total)
  • Author
  • #1736

    I notice that the components of the distribution are stored with an overall effect of representing rho-1 instead of rho. From lbHelpers.h… (Notice the extra (- Descriptor::t[iPop]) term)rnrn/// All helper functions are inside this structurerntemplate<typename T, class Descriptor>rnstruct lbDynamicsHelpers {rn /// Computation of equilibrium distributionrn static T equilibrium(int iPop, T rho, const T u[Descriptor::d], const T uSqr) {rn T c_u = T();rn for (int iD=0; iD < Descriptor::d; ++iD) {rn c_u += Descriptor::c[iPop][iD]*u[iD];rn }rn return rho * Descriptor::t[iPop] * (rn (T)1 + Descriptor::invCs2 * c_u +rn Descriptor::invCs2 * Descriptor::invCs2/(T)2 * c_u*c_u -rn Descriptor::invCs2/(T)2 * uSqrrn ) – Descriptor::t[iPop]; // <<

    extra term, not in usual LBM equationsrn }rnrnThis leads to calculations of rho like this member in lbDynamicsHelpersrn (Notice the rho += (T)1; line)…rnrn /// Computation of densityrn static T computeRho(T const* cell) {rn T rho = T();rn for (int iPop=0; iPop < Descriptor::q; ++iPop) {rn rho += cell[iPop];rn }rn rho += (T)1; // <<

    compensates for extra term shown abovern return rho;rn }rnrnI have two questions..rnrn1) Why do you do this? Does this yield a more stable or accurate simulation?rnrn2) Why is pressure calculated the way it is (shown below)? Isn’t the equation of statern p = rho/c^2,rnnotrn p = (rho-1)/c^2rn rnIn DataAnalysis2D.hh (computeRho() should be computing rho not rho+1)rnrntemplate<typename T, template<typename U> class Lattice>rnvoid DataAnalysis2D<T,Lattice>::computePressureField() const {rn if (flags.pressureFieldComputed) return;rn fields.pressureField.construct();rn for (int iX=0; iX<fields.pressureField.getNx(); ++iX) {rn for (int iY=0; iY<fields.pressureField.getNy(); ++iY) {rn fields.pressureField.get(iX,iY) = block.get(iX,iY).computeRho();rn fields.pressureField.get(iX,iY) -= (T)1; // <<

    why is this needed?rn fields.pressureField.get(iX,iY) /= Lattice<T>::invCs2;rn }rn }rn flags.pressureFieldComputed = true;rn}rnrnadvTHANKSance,rnRandy


    Correction:rn I meant for the c^2 term to be multiplied, not divided. rnrnCorrected question 2)…rnrn2) Why is pressure calculated the way it is (shown below)? Isn’t the equation of staternp = rho*c^2,rnnotrnp = (rho-1)*c^2rn


    Dear Randy,rnrn1) Why do you do this? Does this yield a more stable or accurate simulation?rnA: Yes indeed, we win one or more order of magnitude for the representation of the distribution function f which leads to more stability and accuracy. This is because f-f^eq is small and will be distributed around zero while f is distributes around f^eq which is not zero.rn2) Why is pressure calculated the way it is (shown below)? Isn’t the equation of staternp = rho/c^2,rnnotrnp = (rho-1)/c^2rnA: Have a look at some standard text books on LBM or “Lattice Boltzmann Model for the Incompressible Navier–Stokes“. In the later paper of He/Lou it is explained very nicely in my opinion. rnrnMathias


    Dear Mathias,rnrnThank you for your reply.rnrnI need to sit down and calculate some values for f^eq and f to get a better feel for question 1).rnrnI’m still a bit uncertain as to question 2). In He/Lou page 942, right after eq. (A13b) it saysrn P = c_s^2 rho / rho_0rnrnSince you are calculating the actual rho (not rho+1) inrn block.get(iX,iY).computeRho(),rnshould you not need thern fields.pressureField.get(iX,iY) -= (T)1rnexpression in order to calculate the pressure?rnrnRandy


    Dear Randy,rnrnSorry for the late reply. Dividing by rho0 is just doing scaling. You need to do that also for the results you get in OpenLB for which you can use the class UnitConverter. The term -1 comes from the followed implemetation strategy in OpenLB (cf. answer question 1).rnrnMathias


    Mathias,rnrnThank you again for your patience.rnrnI think you and I are missing an important point…rnrnblock.get(iX,iY).computeRho() computes the fluid’s actual RHO not (RHO + 1).rnrnUsing the actual RHO for the fluid, your computed pressure is (I think incorrectly)…rnrnP = (RHO – 1) * (Cs*Cs)rnrnI don’t think this is what you want.rnrnRegards,rnRandyrn


    Hi,rnrnI understand. Sorry for confusing you. By doing it the way “”P = (RHO – 1) * (Cs*Cs)”” you obtain “”P=0″” for “”RHO=1″”. If you dont do that scaling, you will get different “”pressure levels”” for different discretisations. Then, with the help of the unit converter you can add any constant “”pressure level”” you want. rnrnMathias


    Mathias,rnrnThank you! I need to look into the unit converter code.rnrnBTW, the reason I am so interested in this is that I am trying to implement Swift and Osborn’s model:rnrnLattice Boltzmann Simulation of Nonideal Fluids, Phys. Rev. Lett. 75, 4031 (1995)rnLattice Boltzmann Simulations of liquid-gas and binary fluid systems, Phys. Rev. E. 54(5) (1996)rnrnRegards,rnRandy

Viewing 8 posts - 1 through 8 (of 8 total)
  • You must be logged in to reply to this topic.