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())
 This topic has 7 replies, 2 voices, and was last updated 9 years, 10 months ago by Randy Roberts.

AuthorPosts

June 28, 2013 at 4:04 pm #1736Randy RobertsMember
I notice that the components of the distribution are stored with an overall effect of representing rho1 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 = (rho1)/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,rnRandyJune 28, 2013 at 4:10 pm #2105Randy RobertsMemberCorrection: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 = (rho1)*c^2rn
July 2, 2013 at 3:21 pm #2106mathiasKeymasterDear 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 ff^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 = (rho1)/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
July 3, 2013 at 3:20 pm #2107Randy RobertsMemberDear 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
July 30, 2013 at 9:48 am #2108mathiasKeymasterDear 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
July 30, 2013 at 3:22 pm #2112Randy RobertsMemberMathias,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
July 30, 2013 at 3:47 pm #2113mathiasKeymasterHi,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
July 30, 2013 at 3:55 pm #2114Randy RobertsMemberMathias,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 liquidgas and binary fluid systems, Phys. Rev. E. 54(5) (1996)rnrnRegards,rnRandy

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