Drag Force (Coefficient)

[aliiab]

Hello!rnrnI’m interested to simulate the steady fluid flow in a rectangular channel with sphere in its center. So this task is almost the same as example cylinder3d (but I don’t use stl – file). I want to get drag force acting on the sphere.rnrnInput parameters: liquid density = 1000 kg/m^3 (water), so charRho in converter = 1000.0; fluid velocity = 0.01 m/s (so charU = 0.01); dynamic viscosity = 8.9e-04 Pa*s (so charNu = 8.9e-07 m^2/s (dynamic_visvocity/rho)).rnrnAnd everything works fine, but in one moment drag force (or coefficient) starts increasing and becomes equal to the “”-nan””. I want to understand the reason of such behaviour…rnrnSo my code is presented below.rnrnMain function:rnrn

Code:

typedef double T;rn#define DESCRIPTOR D3Q19Descriptorrnrnrn// Parameters for the simulation setuprnconst int N = 1; // resolution of the modelrnconst int M = 1; // time discretization refinementrnconst T Re = 11*1000.0; // Reynolds numberrnconst T maxPhysT = 32.; // max. simulation time in s, SI unitrnconst T L = 0.01/N;rnrnconst T lengthX = 2.0;rnconst T lengthY = 1.0;rnconst T lengthZ = 1.0;rnrnconst T centerSphereX = 0.5*lengthX;rnconst T centerSphereY = 0.5*lengthY;rnconst T centerSphereZ = 0.5*lengthZ;rnrnconst T sphereDiam = 0.2;rnrnrnint main( int argc, char* argv[] )rn{rnrn /// === 1st Step: Initialization ===rn olbInit( &argc, &argv );rn singleton::directories().setOutputDir( “”./tmp/”” );rn OstreamManager clout( std::cout,””main”” );rn // display messages from every single mpi processrn //clout.setMultiOutput(true);rnrn LBconverter<T> converter(rn ( int ) 3, // dimrn ( T ) L, // latticeL_rn ( T ) 1.0e-02, // latticeU_rn ( T ) 8.9e-07, // charNu_rn ( T ) sphereDiam, // charL_rn ( T ) 0.01, // charU_ rn ( T ) 1000.0 // charRho_rn );rn converter.print();rnrnrn////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////rn /// === 2rd Step: Prepare Geometry ===rn Vector<T, 3> extend(lengthX, lengthY, lengthZ);rn Vector<T, 3> origin(0.0, 0.0, 0.0);rn IndicatorCuboid3D<T> cuboid(extend, origin);rnrnrn /// Instantiation of a cuboidGeometry with weightsrn const int noOfCuboids = 16;rnrn CuboidGeometry3D<T> cuboidGeometry(cuboid, L, noOfCuboids);rnrn Vector<T, 3> center(centerSphereX, centerSphereY, centerSphereZ);rn IndicatorSphere3D<T> sphere(center, 0.5*sphereDiam);rnrn /// Instantiation of a loadBalancerrn HeuristicLoadBalancer<T> loadBalancer(cuboidGeometry);rnrn /// Instantiation of a superGeometryrn SuperGeometry3D<T> superGeometry(cuboidGeometry, loadBalancer, 2);rnrn prepareGeometry(converter, superGeometry, sphere);rnrnrn /// === 3rd Step: Prepare Lattice ===rn SuperLattice3D<T, DESCRIPTOR> sLattice( superGeometry );rn BGKdynamics<T, DESCRIPTOR> bulkDynamics( converter.getOmega(), instances::getBulkMomenta<T, DESCRIPTOR>() );rnrn // choose between local and non-local boundary conditionrn sOnLatticeBoundaryCondition3D<T,DESCRIPTOR> sBoundaryCondition( sLattice );rn //createInterpBoundaryCondition3D<T,DESCRIPTOR>( sBoundaryCondition );rn createLocalBoundaryCondition3D<T,DESCRIPTOR>(sBoundaryCondition);rnrn sOffLatticeBoundaryCondition3D<T, DESCRIPTOR> sOffBoundaryCondition( sLattice );rn createBouzidiBoundaryCondition3D<T, DESCRIPTOR> ( sOffBoundaryCondition );rnrn prepareLattice( sLattice, converter, bulkDynamics, sBoundaryCondition, sOffBoundaryCondition, superGeometry, sphere );rnrn /// === 4th Step: Main Loop with Timer ===rn clout << “”starting simulation…”” << endl;rn Timer<T> timer( converter.numTimeSteps( maxPhysT ), superGeometry.getStatistics().getNvoxel() );rn timer.start();rnrn for ( int iT = 0; iT < converter.numTimeSteps( maxPhysT ); ++iT ) {rnrn /// === 5th Step: Definition of Initial and Boundary Conditions ===rn setSteadyBoundaryValues( sLattice, converter, superGeometry, iT );rnrn /// === 6th Step: Collide and Stream Execution ===rn sLattice.collideAndStream();rnrn /// === 7th Step: Computation and Output of the Results ===rn getResults( sLattice, converter, iT, superGeometry, sphere, timer );rn }rnrn timer.stop();rn timer.printSummary();rn}

rnrnGeometry preparation (I think, there are not any problems in this function, I checked the geometry in Paraview):rnrn

Code:

void prepareGeometry( LBconverter<T> const& converter, SuperGeometry3D<T>& superGeometry, IndicatorSphere3D<T>& sphere )rn{rnrn OstreamManager clout( std::cout,””prepareGeometry”” );rn clout << “”Prepare Geometry …”” << std::endl;rnrn superGeometry.rename(0, 2);rnrn superGeometry.rename(2, 1, 1, 1, 1);rnrn Vector<T,3> minPhysR1 = superGeometry.getStatistics().getMinPhysR( 1 );rn Vector<T,3> maxPhysR1 = superGeometry.getStatistics().getMaxPhysR( 1 );rnrn Vector<T,3> minPhysR2 = superGeometry.getStatistics().getMinPhysR( 2 );rn Vector<T,3> maxPhysR2 = superGeometry.getStatistics().getMaxPhysR( 2 );rnrn Vector<T,3> origin(minPhysR2[0] – L, minPhysR1[1], minPhysR1[2]);rn Vector<T,3> extend(minPhysR1[0] – origin[0], maxPhysR1[1] – origin[1], maxPhysR1[2] – origin[2]);rnrnrn /// Set material number for inflowrn IndicatorCuboid3D<T> inflow(extend, origin);rn superGeometry.rename(2, 3, 1, inflow);rnrn /// Set material number for outflowrn origin[0] = maxPhysR2[0] – L;rn origin[1] = minPhysR1[1];rn origin[2] = minPhysR1[2];rnrn extend[0] = maxPhysR2[0] + L – origin[0];rn extend[1] = maxPhysR1[1] – origin[1];rn extend[2] = maxPhysR1[2] – origin[2];rnrn IndicatorCuboid3D<T> outflow(extend, origin);rn superGeometry.rename(2, 4, 1, outflow);rnrn /// Set material number for spherern superGeometry.rename(1, 5, sphere);rnrn /// Removes all not needed boundary voxels outside the surfacern superGeometry.clean();rn superGeometry.checkForErrors();rnrn superGeometry.print();rnrn clout << “”Prepare Geometry … OK”” << std::endl;rn}

rnrnLattice preparation:rnrn

Code:

void prepareLattice( SuperLattice3D<T,DESCRIPTOR>& sLattice,rn LBconverter<T> const& converter,rn Dynamics<T, DESCRIPTOR>& bulkDynamics,rn sOnLatticeBoundaryCondition3D<T,DESCRIPTOR>& bc,rn sOffLatticeBoundaryCondition3D<T,DESCRIPTOR>& offBc,rn SuperGeometry3D<T>& superGeometry,rn IndicatorSphere3D<T>& sphere )rn{rnrn OstreamManager clout( std::cout,””prepareLattice”” );rn clout << “”Prepare Lattice …”” << std::endl;rnrn const 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 –>bounce backrn sLattice.defineDynamics( superGeometry, 2, &instances::getBounceBack<T, DESCRIPTOR>() );rnrn /// Material=3 –>bulk dynamics (inflow)rn sLattice.defineDynamics( superGeometry, 3, &bulkDynamics );rnrn /// Material=4 –>bulk dynamics (outflow)rn sLattice.defineDynamics( superGeometry, 4, &bulkDynamics );rnrn /// Setting of the boundary conditionsrn bc.addVelocityBoundary( superGeometry, 3, omega );rn bc.addPressureBoundary( superGeometry, 4, omega );rnrn /// Material=5 –>bouzidirn sLattice.defineDynamics( superGeometry, 5, &instances::getNoDynamics<T,DESCRIPTOR>() );rn offBc.addZeroVelocityBoundary( superGeometry, 5, sphere );rnrn /// Initial conditionsrn AnalyticalConst3D<T,T> rhoF( 1 );rn std::vector<T> velocity( 3,T() );rn AnalyticalConst3D<T,T> uF( velocity );rnrn // Initialize all values of distribution functions to their local equilibriumrn sLattice.defineRhoU( superGeometry, 1, rhoF, uF );rn sLattice.iniEquilibrium( superGeometry, 1, rhoF, uF );rn sLattice.defineRhoU( superGeometry, 3, rhoF, uF );rn sLattice.iniEquilibrium( superGeometry, 3, rhoF, uF );rn sLattice.defineRhoU( superGeometry, 4, rhoF, uF );rn sLattice.iniEquilibrium( superGeometry, 4, rhoF, uF ); rnrn /// Make the lattice ready for simulationrn sLattice.initialize();rnrn clout << “”Prepare Lattice … OK”” << std::endl;rn}

rnrnrnBoundary condition (for the first step only):rn

Code:

rn/// Generates a steady inflow for the first steprnvoid setSteadyBoundaryValues( SuperLattice3D<T, DESCRIPTOR>& sLattice, LBconverter<T> const& converter, SuperGeometry3D<T>& superGeometry, int iT )rn{rnrn OstreamManager clout( std::cout,””setSteadyBoundaryValues”” );rn rn if (iT==0) {rn std::vector<T> velocity(3,T(0));rn velocity[0] = 1.0*converter.getLatticeU();rn AnalyticalConst3D<T,T> uF(velocity);rn sLattice.defineU(superGeometry, 3, uF);rn }rn rn}

rnrnTo get drag force I use functor:rn

Code:

SuperLatticePhysDrag3D<T,DESCRIPTOR> drag( sLattice, superGeometry, 5, converter)

rnrnBut here (http://optilb.org/openlb/forum?mingleforumaction=viewtopic&t=183) I found that this functor provides not drag force, but drag coefficient. Is it right?rnrnAt the start of simulation I got such result:rn[LBconverter] LBconverter informationrn[LBconverter] characteristical valuesrn[LBconverter] Dimension(d): dim=3rn[LBconverter] Characteristical length(m): charL=0.2rn[LBconverter] Characteristical speed(m/s): charU=0.01rn[LBconverter] Characteristical time(s): charT=20rn[LBconverter] Density factor(kg/m^d): charRho=1000rn[LBconverter] Characterestical mass(kg): charMass=8rn[LBconverter] Characterestical force(N): charForce=0.004rn[LBconverter] Characterestical pressure(Pa): charPressure=0.1rn[LBconverter] Pressure level(Pa): pressureLevel=0rn[LBconverter] Phys. kinematic viscosity(m^2/s): charNu=8.9e-07rn[LBconverter] lattice valuesrn[LBconverter] DeltaX: deltaX=0.05rn[LBconverter] Lattice velocity: latticeU=0.01rn[LBconverter] DeltaT: deltaT=0.0005rn[LBconverter] Reynolds number: Re=2247.19rn[LBconverter] DimlessNu: dNu=0.000445rn[LBconverter] Viscosity for computation: latticeNu=8.9e-05rn[LBconverter] Relaxation time: tau=0.500267rn[LBconverter] Relaxation frequency: omega=1.99893rn[LBconverter] conversion factorsrn[LBconverter] latticeL(m): latticeL=0.01rn[LBconverter] Time step (s): physTime=0.01rn[LBconverter] Velocity factor(m/s): physVelocity=1rn[LBconverter] FlowRate factor(m^d/s): physFlowRate=0.0001rn[LBconverter] Mass factor(kg): physMass=0.001rn[LBconverter] Force factor(N): physForce=0.1rn[LBconverter] Force factor massless(N/kg): physMasslessForce=100rn[LBconverter] Pressure factor(Pa): physPressure=4000rn[LBconverter] latticePressure: latticeP=0.001rnrnMaybe there are some problems in input numerical parameters (density, viscosity, …)? I see that Reynolds number is equal to the 2247.19, but I excepted that Re will be about 10^4… How is it calculated?rnrnThanks,rnAlina.rn

[mathias]

Dear Alina,rnrnyou have chosen CharL=0.2m, CharU=0.01m/s and CharNu=8.9e-07 m^2/s in the Converter. Thus, you get a Re=2.247.19. Your Re=10.000 does not enter the code anywhere.rnrnTo get stable results you need to choose much smaller grids and time steps or you can use a turbulent model.rnrnBestrnMathias

[aliiab]

Thanks, Mathias. Yes, I know than my Re is not used in the code. But I expected that calculated Re = CharRo*CharU*CharL/CharNu. And now I see that OpenLB calculates Re correctly using this formula. Sorry, it was my mistake in calculations (probably I used another value of sphereDiam (and charL) when calculate Re value manually).rnrnThank you for you advice about grid and turbulent model. I will experiment.rnrnThanks,rnAlina.

[aliiab]

One more question. Are there any criterias for best space and time resolutions?rnrnI know that one can use CharU and LatticeU to compute the time step size in SI units (seconds) as:rn

Quote:

physTime = latticeU/charU * latticeL

(http://optilb.org/openlb/forum?mingleforumaction=viewtopic&t=78)rnrnI suppose that it’s possible to choose appropriate latticeU and latticeL values to get good time step. But what value time step has to be equal to in concrete physical conditions (Re, Rho, charU)?rn

[mathias]

There is the CFL condition. But there is also a lot of experiance due to boundary conditions, the LBM modell, the actual flow domain, …

[aliiab]

Mathias, thanks.rnOne more question, please:) I have detected, that for low kinematic viscosity values (~10^(-7) – 10^(-9)) one needs to make very-very small grids than for biger viscosity (~10^(-3) – 10^(-5)). I want to know is there a exact relation between viscosity value and latticeL value. My suggestion is that If one reduces viscosity value by 3 orders one has to reduce latticeL value by 3 orders too…

[robin.trunk]

Hi Alina,rnrnhaving a look at the relaxation time and viscosity one can see that the system gets unstable if the relaxation time tau gets close to 0.5. You can have a look how tau is computed in the LBconverter:rntau = (U_lattice * char_nu) / (L_lattice * U_char)rnSo your suggestion is correct to keep the value of tau constant.rnrnBestrnRobin

[aliiab]

Robin, thanks. The last question. Can you tell me the approximation dependence OpenLB simulation time on the L_lattice value. In fact, I just want to keep values of the tau, u_lattice, u_char constant. And I want to reduce char_nu by 3 orders. So I have to reduce L_lattice by 3 orders too. So I want to estimate the simulation time increase after theese changes.

[robin.trunk]

Hi Alina,rnrnif you are doing a 3D simulation, you can estimate that by dividing the L_lattice by 3, your number of cells will increase by 3^3=27. Depending on what you simulate (multiphase or single phase) the rquired simulation time will increase by this factor.rnrnBestrnRobin

[aliiab]

Ok, thank you. I hoped that the dependence is not linear:). Thank you, Robin.

[Author]

Posts