Case 22 solved ok and is stable but epsilon is very large at 45 degrees or so: 38518.
Maybe this is ok. It implies a small length scale. The cells ARE very small in the
region of high epsilon. It implies a very low turbulence viscosity at this point.
-----------------------------------------------------
FoamFile{version 2.0; format ascii; class dictionary;
object transportProperties;}
transportModel Newtonian;
nu nu [0 2 -1 0 0 0 0] 1.5E-5;
// kinematic viscosity. We are incompressible so no density required.
// The kinematic pressures generated must be multiplied by density to
// get actual pressures.
-----------------------------------------------------
FoamFile{version 2.0; format ascii; class dictionary;
object RASProperties;}
// Use erroneous RASModel to get list of available models.
RASModel LaunderGibsonRSTM;
// RASModel kOmegaSST;
printCoeffs on;
turbulence on;
// Can define non-standard coeffs.
-----------------------------------------------------
FoamFile{version 2.0; format ascii; class volScalarField; object epsilon;}
internalField uniform 1;
boundaryField{
in { type freestream; freestreamValue uniform 0.05; value uniform 0.05; }
out { type freestream; freestreamValue uniform 0.05; value uniform 0.05; }
sym { type symmetryPlane; }
top { type symmetryPlane; }
wall { type zeroGradient; }
emptyf { type empty; }
emptyb { type empty; }
}
--------------------------------------------------------
FoamFile{version 2.0; format ascii; class volScalarField; object k;}
boundaryField {
in{type freestream; freestreamValue uniform 0.05;}
out{type freestream; freestreamValue uniform 0.05;}
top{type symmetryPlane;}
wall{type zeroGradient;} // zeroGradient? calculated? fixedValue; value uniform 0? 1E-10?
sym{type symmetryPlane;}
emptyf{type empty;}
emptyb{type empty;}
}
internalField uniform 0;
dimensions [0 2 -2 0 0 0 0];
--------------------------------------------------------
FoamFile{version 2.0; format ascii; class volScalarField; object p;}
boundaryField { in{type freestreamPressure;} // zeroGradient
out{type freestreamPressure;} // fixedValue; value uniform 0
top{type symmetryPlane;}
wall{type zeroGradient;}
sym{type symmetryPlane;}
emptyf{type empty;}
emptyb{type empty;}
}
internalField uniform 0;
dimensions [0 2 -2 0 0 0 0];
-----------------------------------------------------
FoamFile{version 2.0; format ascii; class volVectorField; object U;}
boundaryField{
wall{ type fixedValue; value uniform (0 0 0);}
top{ type symmetryPlane;}
sym{ type symmetryPlane;}
emptyf{ type empty;}
emptyb{ type empty;}
in{ type freestream; freestreamValue uniform (10 0 0); } // fixedValue; value uniform (10 0 0);}
out{ type freestream; freestreamValue uniform (10 0 0);} // zeroGradient;}
}
internalField uniform (10 0 0);
dimensions [0 1 -1 0 0 0 0];
-----------------------------------------------------
FoamFile{version 2.0; format ascii; class volSymmTensorField; object R;}
boundaryField{
wall{ type zeroGradient; } // fixedValue; value uniform (0 0 0 0 0 0);}
top{ type symmetryPlane;}
sym{ type symmetryPlane;}
emptyf{ type empty;}
emptyb{ type empty;}
in{ type freestream; freestreamValue uniform (0 0 0 0 0 0); } // fixedValue; value uniform (10 0 0);}
out{ type freestream; freestreamValue uniform (0 0 0 0 0 0);} // zeroGradient;}
}
internalField uniform (0 0 0 0 0 0);
dimensions [0 2 -2 0 0 0 0];
------------------------------------------------------
FoamFile{version 2.0; format ascii; class dictionary; location "system"; object fvSolution;}
SIMPLE {
nNonOrthogonalCorrectors 4; // Unpredictable. 0/1/2 blows. 3,4,5,6 ok with 2D.
convergence 1e-5; // convergence criteria steady state
pRefCell 0; // ref cell is 0 in the airfoil tutorial
pRefValue 0; // ref value is 0 in the airfoil case.
residualControl{p 1E-5; U 1E-5; R 1E-7; epsilon 1E-7;} // added due to airfoil case.
}
relaxationFactors {
// Critical section. Steady solutions will diverge with the default
// 1 setting for a value to make sure they're all in here.
// p is usually just under 1/2 of U, k, omega. Not sure about
// R and epsilon. p 0.3, U 0.7 for good case. p 0.1, U 0.3 considered
// a bad setting by some but slow and steady beats fast and blown...
p 0.1; // 0.3 is stable, decrease for bad mesh
U 0.3; // 0.7 is stable, decrease for bad mesh
R 0.1;
epsilon 0.1;
}
solvers {
U { // linear equation system solver for U
solver smoothSolver; // solver type. smoothSolver is more
// robust than some.
// smoother GaussSeidel; // smoother type. OK serial. Bad parallel.
smoother DILUGaussSeidel; // Good parallel.
tolerance 1e-06; //
relTol 0.01; // Some say this should = 0
nSweeps 1; // setting for smoothSolver
maxIter 100; // limitation of iterations number
}
p // linear equation system solver for p
{
solver GAMG; // very efficient multigrid solver
tolerance 1e-07; // solver finishes if this
relTol 0.001; // or this is reached
// some suggest this should be zero. Others say 0.1 or 0.01 for
// poor convergence cases.
minIter 3; // minimum number of iterations
maxIter 100; // limitation of iterions number
smoother DIC; // setting for GAMG.
// smoother DICGaussSeidel; // Good for parallel runs.
nPreSweeps 1; // 1 for p, set to 0 for all other!
nPostSweeps 2; // 2 is fine
nFinestSweeps 2; // 2 is fine
scaleCorrection true; // true is fine
directSolveCoarsestLevel false; // false is fine
cacheAgglomeration on; // on is fine; set to off, if dynamic
// mesh refinement is used!
nCellsInCoarsestLevel 500; // 500 is fine,
// otherwise sqrt(number of cells)
agglomerator faceAreaPair; // faceAreaPair is fine
mergeLevels 1; // 1 is fine
}
//epsilon{solver PBiCG; preconditioner DILU; tolerance 1E-9; relTol 0;}
//R{solver PBiCG; preconditioner DILU; tolerance 1E-9; relTol 0;}
epsilon{solver smoothSolver; smoother DILUGaussSeidel; tolerance 1E-9;
relTol 0; nSweeps 1; maxIter 100;}
R{solver smoothSolver; smoother DILUGaussSeidel; tolerance 1E-9;
relTol 0; nSweeps 1; maxIter 100;}
}
-------------------------------------------------------------
FoamFile{version 2.0; format ascii; class dictionary;
location "system"; object fvSchemes;}
snGradSchemes { // surface normal gradient on faces calculation methods.
default limited 0.7; // limited 0.7; for bad meshes;
// corrected; for good meshes;
// must fit to laplacianSchemes
}
laplacianSchemes { // Diffusion schemes.
default none;
// --------------- The following is my current favourite.
// default Gauss linear corrected; // low non-orthogonality (<30 deg?)
// default Gauss linear limited 0.5; // medium non-orthogonality (30 < 45 deg)
// default Gauss linear limited 0.33; // low non-orthogonality (>45 deg)
// --------------- or use the following.
// laplacian(nuEff,U) Gauss linear corrected; // for good meshes
laplacian(nuEff,U) Gauss linear limited 0.7; // for bad meshes
//laplacian((1|A(U)),p) Gauss linear corrected; // for good meshes
laplacian((1|A(U)),p) Gauss linear limited 0.7; // for bad meshes
// for turbulent flow only:
// laplacian(DkEff,k) Gauss linear corrected;
laplacian(DepsilonEff,epsilon) Gauss linear corrected;
laplacian(DREff,R) Gauss linear corrected;
// laplacian(DnuTildaEff,nuTilda) Gauss linear corrected; // SA only
}
fluxRequired {
default no;
p ;
phi ;
}
gradSchemes {
default Gauss linear; // Gauss linear or leastSquares
grad(p) faceLimited leastSquares 0 1;
grad(U) Gauss linear;
// grad(U) cellLimited Gauss linear 1; // Poor quality mesh.
}
interpolationSchemes { // interpolation used to calculate values on cell faces.
default linear;
U linear;
// interpolate(U) reconCentral phi cellLimited leastSquare 1.0; // new tet mesh interpolator. Use with reconCentral in divSchemes
}
divSchemes {
default none;
div(phi,U) Gauss upwind; // stable: Gauss upwind
// div(phi,U) Gauss linearUpwind Gauss; // faster linearUpwind Gauss
// div(phi,U) Gauss linearUpwindV Gauss linear; // *
// div(phi,U) Gauss reconCentral cellLimited leastSquares 1.0; // new, tets mesh setting
// div(phi,k) Gauss upwind; // for turbulent flow only
// div(phi,omega) 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 turbulent flow only
// div(phi,nuTilda) Gauss upwind; // for turbulent SA flow only
div((nuEff*dev(grad(U).T()))) Gauss linear; // this scheme is not used for
// steady and laminar flow, but
// it is read by the solver, so
// it must be defined.
div((nuEff*dev(T(grad(U))))) Gauss linear; // not advised, taken from motorbike case
}
ddtSchemes {default steadyState;} // no time dependence here
-------------------------------------------------------
In case 23 I calculated k and epsilon for 1% ti. It blew.
In case 24 I increased input epsilon from 0.05 to 0.06. Finish epsilon changed from 38800 to 38500.
It decreased as I increased the input value.
I calculated a new input epsilon using the cylinder diameter instead of the tunnel as
0.017. 24 was updated and restarted at 0.04 s. The out region is not forcing the
value after all; it is marked as,
out
{
type freestream;
freestreamValue uniform 0.017;
value nonuniform List
70
(
17.1646
16.6899
16.1215
15.448
... more items
) ;
)
varying from 18 to 5E-15; so it IS being calculated in spite of me specifying a freeStream value for it.
That is the point of the freeStream, perhaps. Outflow is calculated, inflow is specified.
The maximum epsilon increased as the input was decreased.
--------------------------------------------------------------
I think this is ok. However, starting with other values of k and epsilon
resulted in the solution blowing up. Perhaps different solvers would help.
Also, what about the k value used at the wall? Zero blows the solution;
1E-6 as used here results in ok results sometimes. First of all, try
dropping it to 1E-10 in 25.
I'm using smoothSolver for U but PBiCG for epsilon and R. Why not smoothsolver for all?
Also, I believe FLUENT uses cellLimited approaches; this might be interesting.
Try later.
case 25 is as case 24 but with 1E-10 for k on the wall. Inspection of
k on the wall shows it is NOT being set; it is being calculated. So that's
another experiment gone.
Take then, the unstable solution and try with different solvers; looking
for stability.
25: -----------------------------------------
Copied from 23 which blew by 0.0005. Solver changes:
epsilon{solver PBiCG; preconditioner DILU; tolerance 1E-9; relTol 0;}
->
epsilon{solver smoothSolver; smoother GaussSeidel; tolerance 1e-6; relTol 0.001; nSweeps 1; maxIter 100;}
This blew by 0.00367. Use same tolerances:
epsilon{solver smoothSolver; smoother GaussSeidel; tolerance 1e-9; relTol 0; nSweeps 1; maxIter 100;}
This blew at 0.000223. Earlier than the relaxed tolerance case. Do the same for R:
R{solver PBiCG; preconditioner DILU; tolerance 1E-9; relTol 0;}
->
R{solver smoothSolver; smoother GaussSeidel; tolerance 1e-9; relTol 0; nSweeps 1; maxIter 100;}
This blew at 0.000233. Time-step continuity errors. Epsilon spent 100 iterations at time=0.000222
and the time-step continuity became poor. I'll try changing nSweeps to 2 for all smooth solved
parameters. It ran for longer, to 0.000292. Still blew though.
Perhaps it's the p that needs attention......
fvSolution:SIMPLE had residualControl{p 1E-5;} Changed to 1E-7. Blew at the same point.
So we still have: high epsilon, low k start and cellLimited solution.
The smoother was changed from GaussSeidel to DILUGaussSeidel after an attempt to use the `better with poor meshes, DICGaussSeidel'. This blew, but slowly. Good!
cellLimited: This is a change to fvSchemes. interpolationSchemes{default linear; U linear;}
was changed to {default linear;} It was not possible to limit this.
Note: grad(p) is faceLimited leastSquares 0 1; why?
This simulation worked ok with specific turbulence parameters:
epsilon in = 0.05; k in = 0.0015;
It blew with epsilon in = 0.0007; k in = 0.005, That's progress. Drop k in, increase epsilon in to stabilise.
