Advices to simulate turbulent channel flow

[Alejandro_Claro]

Hello everyone,rnrnI have been working with OpenLB since a few months. I am trying to implement a turbulent inclined channel flow in a 2D case. Unfortunately, my simulation always diverges, the density and energy give a “”nan”” values for the first hundred time-steps 🙁 . The configuration of my case is:rn

rnReynolds number: 5600rnFroude number: 0.01rnTurbulence model: ShearSmagorinskyForcedBGKdynamicsrnSmagorinsky parameter: 0.16rnGravity acceleration = 9.81 m/srnLatticeL = 0.001 // In order to be at y+=1-2rnLatticeU = LatticeL // In order to be deltaT -> deltaX²rnCharU is function of the Froude numberrnCharL = is defined as the hydraulic diameter (4A/P). Its value is 0.4 m in this case.rnCharRho = 1000 kg/m³rnCharNu = 0.000001 m²/s // water at 20°CrnThe channel is a 0.2×0.4×2.0 m (height x width x long).rn

rnrnI used the forced-poiseuille2D example as my reference to implement my case.rnrnThe Geometry have been defined as:rn

Code:

rnvoid prepareGeometry(LBconverter<T> const& converter,rn SuperGeometry2D<T>& superGeometry) {rnrn OstreamManager clout(std::cout,””prepareGeometry””);rn clout << “”Prepare Geometry …”” << std::endl;rnrn superGeometry.rename(0,2);rnrn std::vector<T> extend(2,T());rn extend[0] = lx; // 2.0rn extend[1] = ly – 1.8*converter.getLatticeL(); //0.2 – 1.8*0.001rn std::vector<T> origin(2,T());rn origin[1] = 0.9*converter.getLatticeL();rn IndicatorCuboid2D<T> cuboid2(extend, origin);rnrn superGeometry.rename(2,1,cuboid2);rn rn origin[1] = ly – 0.9*converter.getLatticeL();rn extend[1] = converter.getLatticeL();rn IndicatorCuboid2D<T> TopWall(extend, origin);rn rn superGeometry.rename(2,3,TopWall);rn rn origin[1] = 0.9*converter.getLatticeL();rn extend[0] = 0.9*converter.getLatticeL();rn extend[1] = ly – 1.8*converter.getLatticeL();rn IndicatorCuboid2D<T> Inlet(extend, origin);rn rn superGeometry.rename(1,4,Inlet);rn rn origin[0] = lx – 0.9*converter.getLatticeL();rn extend[0] = converter.getLatticeL();rn IndicatorCuboid2D<T> Outlet(extend, origin);rn rn superGeometry.rename(1,5,Outlet);rn rn /// Removes all not needed boundary voxels outside the surfacern superGeometry.clean();rn /// Removes all not needed boundary voxels inside the surfacern superGeometry.innerClean();rn superGeometry.checkForErrors();rnrn superGeometry.print();rnrn clout << “”Prepare Geometry … OK”” << std::endl;rn

rnrnI want to simulated a inclined channel with a horizontal numerical channel. Thus, I need to define the gravity acceleration in the “”x”” and “”y”” direction by means of a slope. This slope will be calculated using the Manning equation and the Froude number. I know that for low Froude number this approach does not hold but it is the simplest one to obtain this value. rnrnThe lattice is prepared as:rn

Code:

rnvoid prepareLattice(LBconverter<T> const& converter,rn SuperLattice2D<T, DESCRIPTOR>& sLattice,rn Dynamics<T, DESCRIPTOR>& bulkDynamics,rn sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& sBoundaryCondition,rn SuperGeometry2D<T>& superGeometry ) {rnrn OstreamManager clout(std::cout,””prepareLattice””);rn clout << “”Prepare Lattice …”” << std::endl;rnrn T omega = converter.getOmega();rnrn /// Material=0 –>do nothingrn sLattice.defineDynamics(superGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>());rnrn /// Material=1 –>bulk dynamicsrn sLattice.defineDynamics(superGeometry, 1, &bulkDynamics);rnrn /// Material=2 –>bulk dynamicsrn sLattice.defineDynamics(superGeometry, 2, &bulkDynamics);rnrn /// Setting of the boundary conditionsrn sBoundaryCondition.addVelocityBoundary(superGeometry, 2, omega);rn rn /// Material=3 –> bulk dynamicsrn sLattice.defineDynamics(superGeometry, 3, &bulkDynamics);rn /// Setting of the boundary conditionsrn sBoundaryCondition.addSlipBoundary(superGeometry, 3);rn rn /// Material=4 –> bulk dynamicsrn sLattice.defineDynamics(superGeometry, 4, &bulkDynamics);rn /// Setting of the boundary conditionsrn sBoundaryCondition.addVelocityBoundary(superGeometry, 4, omega);rn rn /// Material=5 –> bulk dynamicsrn sLattice.defineDynamics(superGeometry, 5, &bulkDynamics);rn /// Setting of the boundary conditionsrn //sBoundaryCondition.addConvectionBoundary(superGeometry, 5, omega);rn sBoundaryCondition.addPressureBoundary(superGeometry, 5, omega); rnrn /// Initial conditionsrn /// Initial avShear parameterrn std::vector<T> avShearVec(1,T(0.1));rn AnalyticalConst2D<T,T> avShear(avShearVec);rn sLattice.defineExternalField(superGeometry, 1,rn DESCRIPTOR<T>::ExternalField::avShearIsAt,rn DESCRIPTOR<T>::ExternalField::sizeOfAvShear, avShear );rn sLattice.defineExternalField(superGeometry, 2,rn DESCRIPTOR<T>::ExternalField::avShearIsAt,rn DESCRIPTOR<T>::ExternalField::sizeOfAvShear, avShear );rn sLattice.defineExternalField(superGeometry, 3,rn DESCRIPTOR<T>::ExternalField::avShearIsAt,rn DESCRIPTOR<T>::ExternalField::sizeOfAvShear, avShear );rn sLattice.defineExternalField(superGeometry, 4,rn DESCRIPTOR<T>::ExternalField::avShearIsAt,rn DESCRIPTOR<T>::ExternalField::sizeOfAvShear, avShear );rn sLattice.defineExternalField(superGeometry, 5,rn DESCRIPTOR<T>::ExternalField::avShearIsAt,rn DESCRIPTOR<T>::ExternalField::sizeOfAvShear, avShear );rnrn rn std::vector<T> extend(2,T()), origin(2,T());rn origin[0] = 0.9*converter.getLatticeL();rn origin[1] = 0.9*converter.getLatticeL();rn extend[0] = lx – 1.8*converter.getLatticeL();rn extend[1] = ly – 1.8*converter.getLatticeL();rn IndicatorCuboid2D<T> bulk(extend, origin);rn rn origin[0] = (T)0.;rn //origin[1] = 0.9*converter.getLatticeL();rn extend[0] = 0.9*converter.getLatticeL();rn //extend[1] = ly – 1.8*converter.getLatticeL();rn IndicatorCuboid2D<T> Inlet(extend, origin);rn rn origin[0] = lx – 0.9*converter.getLatticeL();rn //origin[1] = 0.9*converter.getLatticeL();rn extend[0] = converter.getLatticeL();rn //extend[1] = ly – 1.8*converter.getLatticeL();rn IndicatorCuboid2D<T> Outlet(extend, origin);rn rn origin[0] = (T)0.;rn origin[1] = -0.9*converter.getLatticeL();rn extend[0] = lx;rn extend[1] = converter.getLatticeL();rn IndicatorCuboid2D<T> Bottom_Wall(extend, origin);rn rn //origin[0] = (T)0.;rn origin[1] = ly – 0.9*converter.getLatticeL();rn //extend[0] = lx;rn //extend[1] = converter.getLatticeL();rn IndicatorCuboid2D<T> Top_Wall(extend, origin);rn rn T GravityBodyForceScalar = (T)((converter.getLatticeU()*converter.getLatticeU())/(converter.getCharU()*converter.getCharU()))rn *converter.getDeltaX()*converter.getCharL()*gravityconst; rn // Body force in lattice unitsrn std::vector<T> GravityBodyForce(2,T());rn rn T ks = (T)100.; // Manning rugosity coefficientrn T A = width*ly; // 0.4 x 0.2rn T P = width+(2.*ly); // 0.4 + 2*0.2rn T Rh = A/P; // 0.4rnrn T Slope_1 = pow(converter.getCharU()/(ks*pow(Rh,(2/3))),2);rn rn GravityBodyForce[0] = sin(atan(Slope_1))*GravityBodyForceScalar;rn GravityBodyForce[1] = cos(atan(Slope_1))*GravityBodyForceScalar; rn rn AnalyticalLinear2D<T,T> force_y(0., GravityBodyForce[1], GravityBodyForce[1]*(-ly));rn AnalyticalLinear2D<T,T> force_x(0., GravityBodyForce[0], GravityBodyForce[0]*(-ly));rn AnalyticalComposed2D<T,T> force(force_x, force_y); // in order to account the force at each direction.rn rn // Initialize forcern sLattice.addExternalField(superGeometry, bulk,rn DESCRIPTOR<T>::ExternalField::forceBeginsAt,rn DESCRIPTOR<T>::ExternalField::sizeOfForce, force );rn sLattice.addExternalField(superGeometry, Bottom_Wall,rn DESCRIPTOR<T>::ExternalField::forceBeginsAt,rn DESCRIPTOR<T>::ExternalField::sizeOfForce, force );rn sLattice.addExternalField(superGeometry, Top_Wall,rn DESCRIPTOR<T>::ExternalField::forceBeginsAt,rn DESCRIPTOR<T>::ExternalField::sizeOfForce, force );rn sLattice.addExternalField(superGeometry, Inlet,rn DESCRIPTOR<T>::ExternalField::forceBeginsAt,rn DESCRIPTOR<T>::ExternalField::sizeOfForce, force );rn sLattice.addExternalField(superGeometry, Outlet,rn DESCRIPTOR<T>::ExternalField::forceBeginsAt,rn DESCRIPTOR<T>::ExternalField::sizeOfForce, force );rnrn /// Make the lattice ready for simulationrn sLattice.initialize();rnrn clout << “”Prepare Lattice … OK”” << std::endl;rn}rn

rnrnThe inlet of the domain (material number = 4) is defined as a velocity boundary condition. The outlet (material number = 5) is defined as a pressure boundary. The Top-Wall is defined as a rigid lid condition using the “”freeslip”” boundary condition, this condition is based in specular reflection. The Bottom-Wall boundary condition is defined as a velocity boundary condition with fixed velocity equal to zero. All the boundaries conditions are implemented with the interpolate approach (Skodors 1993 with regularized formula). The fluid, as I said at the beginning of the post, is treated with the shear Smagorinsky turbulence model (ShearSmagorinskyForcedBGKdynamics).rnrn

Code:

rnvoid setBoundaryValues(SuperLattice2D<T, DESCRIPTOR>& sLattice,rn LBconverter<T> const& converter, int iT,rn SuperGeometry2D<T>& superGeometry)rn{rn rn OstreamManager clout(std::cout,””setBoundaryValues””);rnrn if(iT==0) {rn ////////////////////////////////////////////////////////////////////rn /// Initial conditionsrn AnalyticalConst2D<T,T> rho(1000);rnrn //T maxVelocity = converter.getLatticeU();rn //AnalyticalConst2D<T,T> u(maxVelocity);rn AnalyticalConst2D<T,T> u(0.);rn AnalyticalConst2D<T,T> u_inlet(converter.getLatticeU());rnrn // Initialize all values of distribution functions to their local equilibriumrn sLattice.defineRhoU(superGeometry, 1, rho, u);rn sLattice.defineRhoU(superGeometry, 3, rho, u);rn sLattice.defineRhoU(superGeometry, 4, rho, u_inlet);rn sLattice.defineRhoU(superGeometry, 5, rho, u);rn rn sLattice.iniEquilibrium(superGeometry, 1, rho, u);rn sLattice.iniEquilibrium(superGeometry, 2, rho, u);rn sLattice.iniEquilibrium(superGeometry, 3, rho, u);rn sLattice.iniEquilibrium(superGeometry, 4, rho, u_inlet);rn sLattice.iniEquilibrium(superGeometry, 5, rho, u);rn rn ////////////////////////////////////////////////////////////////////rnrn /// Make the lattice ready for simulationrn sLattice.initialize();rn }rn}rn

rnrnI have the next questions:rnrn1. This initial conditions are well implemented? This could be the cause to get a diverged simulation?rnrn2. I want to simulated the outlet as an outflow condition, I already tested the “”convection-boundary”” option but I continue to obtain “”nan”” values.rnrn3. The body force implementation is incorrect? I know that this force should be in lattice units and I preformed a dimensional analysis to well define this force in lattice units.rnrn4. The spatial resolution of the model is enough to implement the shear Smagorinsky model, rigth?rnrn5. Is it the BGK model the principal source of error for this kind of implementation? I have checked the doxygen and this turbulent model have been implemented only to the “”Forced-BGK-Collision””. Thus, I will need to create a code to implement the same turbulence model but for the case of MRT 😕 rnrnrnI know that it is a long post, but I prefers to give the maximum of elements that I consider as a possible source of error. I really will appreciate anyone who could give me some advises to try to implement this case.rnrnI can shear the all the code if anyone is interested to check what I have tried to implement.rnrnBest regards,rnrnAlejandro rn

[mathias]

Hi,rnrn1/ Reduce to about Re=500. Make your case work with ForcedBGK.rn2/ Then, swich using a turbulent model, like Smago with a high SmagoConst. rn3/ Decrease SmagoConst.rnrnBestrnMathias

[Alejandro_Claro]

Hi Mathias,rnrnThank you for quick reply. rnrnAccording to the turbulence model, I should use the Shear-Smagorinky model with the “”normal”” parameter (Cs = 0.16-0.13) in order to take into account the presence of the wall, where the Smagorinsky parameter tends to decrease.rnrnThe body force term (gravity in this case) have been decomposed in its x/y directions according to the slope of the channel. The simulation diverges when I take in account both directions (in my case “”y””). So, I use just the horizontal component. Do you know which could be the problem? Do I not compute the gravity body force term right? Should I just take into account the horizontal direction of this force because I do not use a “”Volume-of-Fluid”” approach in order to simulate the water level variation ?rnrnPreviously, I used a smooth start curve (PolynomialStartScale<T,T>) in order to begin from a Re = 300 to Re = 3000. All along the simulation I used the shear-smagorinsky turbulence model. However, the simulation diverged. Is it necessary to begin the simulation at low Re without a turbulence model and then to switch to moderate Re with a turbulence model for low spatial resolution?rnrnBest Regards,rnrnAlejandrornrn

[mathias]

Dear Alejandro,rnrn There are several issues to be discussed for that topic. Further, we will have a spring school in March 2017 for LBM and OpenLB, where we can help to get started. Please contact me to proceed.rnrnBestrnMathias

[mathias]

Dear Alejandro,rnrn There are several issues to be discussed for that topic. Further, we will have a spring school in March 2017 for LBM and OpenLB, where we can help to get started. Please contact me to proceed.rnrnBestrnMathias

[Author]

Posts