rho vs. rho +/- 1 (esp. for computePressureField())
- This topic has 7 replies, 2 voices, and was last updated 11 years, 2 months ago by .
-
Posts
-
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,rnRandyDear 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 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
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
- You must be logged in to reply to this topic.