Drag Force (Coefficient)
OpenLB – Open Source Lattice Boltzmann Code › Forums › on OpenLB › General Topics › Drag Force (Coefficient)
- This topic has 9 replies, 3 voices, and was last updated 8 years, 5 months ago by aliiab.
-
AuthorPosts
-
August 30, 2016 at 2:37 pm #1856aliiabMember
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
August 30, 2016 at 3:10 pm #2437mathiasKeymasterDear 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
August 31, 2016 at 6:47 am #2439aliiabMemberThanks, 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.
August 31, 2016 at 8:40 am #2441aliiabMemberOne 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
August 31, 2016 at 7:30 pm #2443mathiasKeymasterThere is the CFL condition. But there is also a lot of experiance due to boundary conditions, the LBM modell, the actual flow domain, …
September 8, 2016 at 1:53 pm #2450aliiabMemberMathias, 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…
September 8, 2016 at 3:21 pm #2451robin.trunkKeymasterHi 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
September 8, 2016 at 6:05 pm #2452aliiabMemberRobin, 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.
September 9, 2016 at 9:19 am #2448robin.trunkKeymasterHi 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
September 9, 2016 at 9:31 am #2447aliiabMemberOk, thank you. I hoped that the dependence is not linear:). Thank you, Robin.
-
AuthorPosts
- You must be logged in to reply to this topic.