The process here is:
1) Specify the mesh.
2) Set flow / material parameters.
3) Solution control:
3a) Set turbulence modelling approach.
3b) Specify solution schemes / types.
3c) Specify solution parameters.
3d) Specify software settings (save rate, start time etc)
4) Create the mesh.
5) Solve.
The turbulent analysis benefits from a laminar pre-solution to set the initial field for the turbulent analysis and progressively complex turbulence models benefit from progressively more accurate initial conditions. For this reason, a Reynolds Stress analysis may benefit from an initial laminar solution, then say a k-epsilon solution to generate the initial field. That's another post, anyway. This post runs from scratch with turbulence.
1) Create the mesh: Convergence studies suggested little change was occurring with further refinement to the boundary layer mesh. The file,
constant/polyMesh/blockMeshDict
FoamFile{version 2.0; format ascii; class dictionary; object blockMeshDict;}
convertToMeters 0.05;
// Vertex dimensions are in diameters!
// space between `vertices' and `(' is required.
vertices (
(0.5 0 0) // cylinder wall
(1 0 0) // b/l
(20 0 0) // exit
(20 0.707107 0)
(0.707107 0.707107 0)
(0.353553 0.353553 0)
(20 10 0)
(0.707107 10 0)
(0 10 0)
(0 1 0)
(0 0.5 0)
(0.5 0 0.5)
(1 0 0.5)
(20 0 0.5) // downstream exit at 2
(20 0.707107 0.5)
(0.707107 0.707107 0.5)
(0.353553 0.353553 0.5)
(20 10 0.5)
(0.707107 10 0.5)
(0 10 0.5)
(0 1 0.5)
(0 0.5 0.5)
(-0.5 0 0)
(-1 0 0)
(-10 0 0)
(-10 0.707107 0)
(-0.707107 0.707107 0)
(-0.353553 0.353553 0)
(-10 10 0)
(-0.707107 10 0)
(-0.5 0 0.5)
(-1 0 0.5)
(-10 0 0.5)
(-10 0.707107 0.5)
(-0.707107 0.707107 0.5)
(-0.353553 0.353553 0.5)
(-10 10 0.5)
(-0.707107 10 0.5)
);
// Defined anti-clockwise on front face, then back face.
blocks(
// cylinder wall, 0 to 45 deg. radial, circumferential.
hex(5 4 9 10 16 15 20 21) (80 10 1) simpleGrading (10 1 1)
// cylinder wall, 45 to 90 deg. radial, circumferential.
hex(0 1 4 5 11 12 15 16) (80 10 1) simpleGrading (10 1 1)
// downstream, behind cylinder.
hex(1 2 3 4 12 13 14 15) (160 10 1) simpleGrading (15 1 1)
// downstream, upper 1/4
hex(4 3 6 7 15 14 17 18) (160 60 1) simpleGrading (15 14 1)
// above cylinder, 45 to 90 deg.
hex(9 4 7 8 20 15 18 19) (10 60 1) simpleGrading (1 14 1)
// cylinder wall, 90 to 135 deg. radial, circumferential
hex(10 9 26 27 21 20 34 35) (80 10 1) simpleGrading (10 1 1)
// cylinder wall, 135 to 180 deg. radial, circ.
hex(27 26 23 22 35 34 31 30) (80 10 1) simpleGrading (10 1 1)
// above cylinder, 90 to 135 deg.
hex(23 26 25 24 31 34 33 32) (10 60 1) simpleGrading (1 14 1)
// upstream, upper 1/4
hex(26 29 28 25 34 37 36 33) (60 60 1) simpleGrading (14 14 1)
// directly upstream of cylinder.
hex(26 9 8 29 34 20 19 37) (10 60 1) simpleGrading (1 14 1)
);
edges(
arc 0 5 (0.469846 0.17101 0)
arc 5 10 (0.17101 0.469846 0)
arc 1 4 (0.939693 0.34202 0)
arc 4 9 (0.34202 0.939693 0)
arc 11 16 (0.469846 0.17101 0.5)
arc 16 21 (0.17101 0.469846 0.5)
arc 12 15 (0.939693 0.34202 0.5)
arc 15 20 (0.34202 0.939693 0.5)
arc 22 27 (-0.469846 0.17101 0)
arc 27 10 (-0.17101 0.469846 0)
arc 23 26 (-0.939693 0.34202 0)
arc 26 9 (-0.34202 0.939693 0)
arc 30 35 (-0.469846 0.17101 0.5)
arc 35 21 (-0.17101 0.469846 0.5)
arc 31 34 (-0.939693 0.34202 0.5)
arc 34 20 (-0.34202 0.939693 0.5) );
// patches ( type name (nodes)
// These are GEOMETRY effect, not flow etc.
// So flow inlet is a patch in this file,
// while physical walls are wall.
patches(
patch inlet(( 28 25 33 36 )
( 25 24 32 33 ))
patch outlet((6 3 14 17)
(3 2 13 14))
symmetry symmetry((24 23 31 32)
(23 22 30 31)
(0 1 12 11)
(1 2 13 12))
wall top((28 29 37 36)
(29 8 19 37)
(8 7 18 19)
(7 6 17 18))
wall wall((22 27 35 30)
(27 10 21 35)
(10 5 16 21)
(5 0 11 16))
empty emptyf((28 29 26 25)
(29 8 9 26)
(26 9 10 27)
(26 27 22 23)
(25 26 23 24)
(8 7 4 9)
(7 6 3 4)
(9 4 5 10)
(4 3 2 1)
(4 1 0 5))
empty emptyb((36 37 34 33)
(37 19 20 34)
(34 20 21 35)
(34 35 30 31)
(33 34 31 32)
(19 18 15 20)
(18 17 14 15)
(20 15 16 21)
(15 14 13 12)
(15 12 11 16)));
mergePatchPairs();
This file creates some skewed cells; these will require
nNonOrthogonalCorrectors = 4or greater. This was found experimentally. This parameter may effect the number of iterations of the pressure solver, per `time-step'.
2) Specify the fluid properties as kinematic viscosity. Specify nu where nu = mu / rho with mu, the viscosity and rho the density. mu = 1.78E-5 kg/(m s) = 1.78E-5 Pa s at 15 degrees centigrade. Air density is 1.225 kg/m3 at this temperature. nu therefore is, 1.78E-5 kg/(ms) / 1.225 kg/m3 = 1.45E-5 m2 s-1. (recall that viscosity, mu sets the proportionality between the force required to move two parallel, massless plates with fluid between them; and the plate area, their relative velocity and separation. F = mu A u/y (well, du/dy). Its units are prescribed by 1/mu = A u / (y F) -> 1/(m2 m s-1 m-1 N-1) = m-2 s N = Pa s). Of similar interest for this study, the Reynolds number Re = u D rho / mu = u D / nu with velocity, u; and diameter, D. For a 0.05m cylinder in 10m/s flow, the Reynolds number is 10m/s 0.05m 1.225 kg m-3 / (1.78E-5 kg m-1 s-1) = 34400. This cylinder will display transition in shear layer behaviour (Zdravkovitch) so the boundary layer will be laminar when it separates and separation will occur around 80 degrees.
constant/transportProperties
FoamFile{version 2.0; format ascii; class dictionary;object transportProperties;}
transportModel Newtonian;
nu nu [0 2 -1 0 0 0 0] 1.45E-5;
3) Solution control. 3a) Set the turbulence treatment for RAS. Create
constant/RASPropertiesSet as laminar for now (turbulence off):
FoamFile{version 2.0; format ascii; class dictionary; object RASProperties;}
RASModel LaunderGibsonRSTM;
turbulence off;
3b) Solution types / schemes are specified in
system/fvSchemes
FoamFile{version 2.0; format ascii; class dictionary; location "system"; object fvSchemes;}
snGradSchemes{
// limited 0.7; for bad meshes, corrected for good meshes.
// must fit to laplacianSchemes
default limited 0.7;
}
laplacianSchemes{
default none; // If something is missing, tell me.
// for good meshes
//laplacian(nuEff,U) Gauss linear corrected;
// for bad meshes
laplacian(nuEff,U) Gauss linear limited 0.7;
// for good meshes
//laplacian((1|A(U)),p) Gauss linear corrected;
// for bad meshes
laplacian((1|A(U)),p) Gauss linear limited 0.7;
// for turbulent flow involving k only:
// laplacian(DkEff,k) Gauss linear corrected;
// for turbulent flow involving epsilon only:
laplacian(DepsilonEff,epsilon) Gauss linear corrected;
// for turbulent flow including Reynolds stress only:
// Needed for LRR. corrected or limited.
laplacian(DREff,R) Gauss linear corrected;
// for turbulent flow requiring nuTilda only, ie SA
// laplacian(DnuTildaEff,nuTilda) Gauss linear corrected;
// These are from $FOAM_TUTORIALS/incompressible/simpleFoam/airFoil2D
// SA not kE
// laplacian(nuEff,U) Gauss linear corrected; // got it
// laplacian((1|A(U)),p) Gauss linear corrected; // got it
// Why not limited? nuTilda for SA only.
// laplacian(DnuTildaEff,nuTilda) Gauss linear corrected;
// didnt have it. limited?
// laplacian(1,p) Gauss linear corrected;
}
fluxRequired{
default no;
p;
//phi;
}
gradSchemes{
// Gauss linear or leastSquares
default Gauss linear;
grad(p) faceLimited leastSquares 0 1;
grad(U) Gauss linear;
}
interpolationSchemes{
default linear; // linear is fine
U linear;
// new tet mesh interpolator. Use with reconCentral in divSchemes.
// This is not available in vanilla openFoam yet. There is an extended
// version. I'll try compile it in another article.
// I've read articles where people stated it gave little benefit,
// but none that stated it helped.
// interpolate(U) reconCentral phi cellLimited leastSquare 1.0;
}
divSchemes{
default none;
// stable: Gauss upwind
div(phi,U) Gauss upwind;
// faster linearUpwind Gauss
// div(phi,U) Gauss linearUpwind Gauss;
// new, tets mesh setting:
// div(phi,U) Gauss reconCentral cellLimited leastSquares 1.0;
// for turbulent flow only
// div(phi,k) Gauss upwind;
// for turbulent flow only
div(phi,epsilon) Gauss upwind;
// for turbulent flow only
div(phi,R) Gauss upwind;
// for turbulent flow only
div(R) Gauss linear;
// for SA turbulent flow only.
// div(phi,nuTilda) Gauss upwind;
// this scheme is not used for steady and laminar flow, but
// it is read by the solver. (I haven't tested this)
div((nuEff*dev(grad(U).T()))) Gauss linear;
// not advised, taken from motorbike case
div((nuEff*dev(T(grad(U))))) Gauss linear;
// from airFoil2D:
// default none; // same
// div(phi,U) Gauss linearUpwind grad(U);
// set above but this I think is the fast version. SA only
// div(phi,nuTilda) Gauss linearUpwind grad(nuTilda);
// I think this is the same as the term above.
// div((nuEff*dev(T(grad(U))))) Gauss linear;
}
ddtSchemes{
default steadyState;
// ddt(epsilon)
}
3c) Specify the solution parameters. Create
system/fvSolution
FoamFile{version 2.0; format ascii; class dictionary; location "system"; object fvSolution;}
SIMPLE{
// 0 is fine for perfect mesh; increase to 1 or 2 for non orthogonal mesh
nNonOrthogonalCorrectors 2;
// convergence criteria steady state
convergence 1e-7;
// ref cell is 0 in the airfoil tutorial
pRefCell 0;
// ref value is 0 in the airfoil case.
pRefValue 0;
residualControl{
p 1E-5;
U 1E-5;
k 1E-5;
epsilon 1E-5;
}
}
relaxationFactors{
// 0.3 is stable, decrease for bad mesh
p 0.1;
// 0.7 is stable, decrease for bad mesh
U 0.3;
// 0.7 is stable, decrease for bad mesh
//k 0.3;
// Not sure yet about R relaxation.
//R 0.1;
// Not sure yet about epsilon relaxation.
//epsilon 0.3;
}
solvers{
// linear equation system solver for U
U{
// solver type
solver smoothSolver;
// smoother type
smoother GaussSeidel;
// solver finishes if either absolute or relative
// tolerance is reached. airFoil2D uses 1E-8
tolerance 1e-06;
relTol 0.01;
// smoothSolver setting.
//: airFoil uses 2 sweeps with 0.1 relTol.
nSweeps 1;
maxIter 100;
}
// linear equation system solver for p
p{
// very efficient multigrid solver
solver GAMG;
// solver finishes if either absolute or relative
// tolerance is reached.
// airFoil2D uses relTol 0.1, tol 1E-6
tolerance 1e-07;
relTol 0.001;
// a minimum number of iterations
minIter 3;
// limitation of iterions number
maxIter 100;
// setting for GAMG. airFoil2D uses GaussSeidel.
smoother DIC;
// 1 for p, set to 0 for all other! airfoil2D uses 0
nPreSweeps 1;
// 2 is fine. airfoil2D uses 2
nPostSweeps 2;
// 2 is fine
nFinestSweeps 2;
scaleCorrection true;
directSolveCoarsestLevel false;
// on unless dynamic mesh refinement is used.
cacheAgglomeration on;
// 500 or sqrt(cell count)
nCellsInCoarsestLevel 1000;
agglomerator faceAreaPair;
mergeLevels 1;
}
//epsilon{solver PBiCG; preconditioner DILU; tolerance 1E-9; relTol 0;}
//R{solver PBiCG; preconditioner DILU; tolerance 1E-9; relTol 0;}
}
3d) Create the control dictionary: file system/controlDict
FoamFile{version 2.0; format ascii; class dictionary; object controlDict;}
application simpleFoam;
startFrom latestTime; // latestTime / startTime
startTime 0;
stopAt endTime; endTime 100000;
deltaT 1;
writeControl timeStep; writeInterval 300; purgeWrite 0;
writeFormat ascii; writePrecision 6; writeCompression uncompressed;
timeFormat general;timePrecision 6;
runTimeModifiable yes;
4) Create the mesh by typing,
$ blockMesh5) Solve with a single processor by typing,
$ nohup simpleFoam > log &or in parallel with
$ mpirun -np 4 simpleFoam -parallel > log &
No comments:
Post a Comment