| piotr.mecht |
September 27, 2019 03:25 |
Highly unstable buoyantPimpleFoam simulation
4 Attachment(s)
Simulation background
========================
I'm trying to perform a buoyantPimpleFoam simulation, but it quickly goes unstable in the first timesteps. The goal is to estimate inlet water jet deviation in domestic hot water tank possibly with different initial conditions and dimensions.
Flow conditions:
initial water temperature (internalField) 42°C
inlet water temperature 20°C
inlet/outlet connector diameter 0.023 m
tank diameter 0.45 m
connectors distance 1.0 m
water flow 3 l/min -> ~0.12 m/s on average (~0.15 m/s max with distributed profile)
Problem solution history
========================
I've tried do reduce timestep, but my Courant number didn't get any smaller.
I also thought my mesh might be bad, because of some concave cells or some single cells with low determinant (checkMesh results), but concave cells are often hard to avoid with snappyHexMesh and they seem to be acceptable in many cases.
I tried to run the same case with fine tetmesh + layers addition + several nOrthogonalCorrectors, but the tendency was the same, so i think now that some other settings might be wrong or they require more tuning.
I blindly played with some fvSolution and fvSchemes, but with similiar results.
Residuals
========================
Enthaply residual in each timestep starts with 1, when static pressure isn't initialized so i helped it with funkySetFields, its initial value reduced significantly but it starts with ~0.7 initial value in second iteration
p_rgh initial residual is >0.4 it seems to be high but it goes serveral orders of magnitude lower in piso and pimple loop (I don't think I understand this output clearly)
Second iteration starts with low max Courant number e.g. ~0.2 , but third or fourth instantly goes to ~100
Omega goes wild in the second iteration with minimum negative value and maximum value of >4.5E+4
Solver crashes whith being unable to converge correct
Output example:
Code:
Starting time loop
Courant Number mean: 1.49349e-08 max: 0.02674
Time = 0.0001
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for Ux, Initial residual = 1, Final residual = 5.48763e-08, No Iterations 3
DILUPBiCG: Solving for Uy, Initial residual = 1, Final residual = 1.80726e-08, No Iterations 3
DILUPBiCG: Solving for Uz, Initial residual = 1, Final residual = 2.18823e-08, No Iterations 3
DILUPBiCG: Solving for h, Initial residual = 0.000321462, Final residual = 1.67432e-07, No Iterations 1
DICPCG: Solving for p_rgh, Initial residual = 0.996965, Final residual = 0.00952103, No Iterations 479
DICPCG: Solving for p_rgh, Initial residual = 0.000555127, Final residual = 5.38467e-06, No Iterations 86
DICPCG: Solving for p_rgh, Initial residual = 5.38737e-06, Final residual = 5.28314e-08, No Iterations 167
DICPCG: Solving for p_rgh, Initial residual = 5.28319e-08, Final residual = 9.66085e-09, No Iterations 392
[...char limit]
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 1.32918e-13, global = -9.31681e-17, cumulative = 3.77663e-16
DICPCG: Solving for p_rgh, Initial residual = 0.0396398, Final residual = 0.000376473, No Iterations 479
DICPCG: Solving for p_rgh, Initial residual = 0.000235303, Final residual = 2.23339e-06, No Iterations 85
DICPCG: Solving for p_rgh, Initial residual = 2.23373e-06, Final residual = 2.03148e-08, No Iterations 170
DICPCG: Solving for p_rgh, Initial residual = 2.0315e-08, Final residual = 9.81896e-09, No Iterations 378
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 1.35542e-13, global = -1.49899e-17, cumulative = 3.62673e-16
DILUPBiCG: Solving for omega, Initial residual = 1.47021e-05, Final residual = 4.75973e-09, No Iterations 1
bounding omega, min: -52.5363 max: 95779.9 average: 224.474
DILUPBiCG: Solving for k, Initial residual = 1, Final residual = 2.35803e-08, No Iterations 3
ExecutionTime = 215.67 s ClockTime = 216 s
Courant Number mean: 1.90016e-05 max: 0.239056
Time = 0.0002
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for Ux, Initial residual = 0.0965107, Final residual = 1.40458e-09, No Iterations 2
DILUPBiCG: Solving for Uy, Initial residual = 0.233133, Final residual = 3.71724e-09, No Iterations 2
DILUPBiCG: Solving for Uz, Initial residual = 0.238017, Final residual = 1.34334e-09, No Iterations 2
DILUPBiCG: Solving for h, Initial residual = 0.608475, Final residual = 9.47492e-10, No Iterations 2
DICPCG: Solving for p_rgh, Initial residual = 0.488246, Final residual = 0.00484784, No Iterations 461
DICPCG: Solving for p_rgh, Initial residual = 2.29527e-05, Final residual = 2.23269e-07, No Iterations 87
DICPCG: Solving for p_rgh, Initial residual = 2.23269e-07, Final residual = 9.84747e-09, No Iterations 85
DICPCG: Solving for p_rgh, Initial residual = 9.84746e-09, Final residual = 9.84746e-09, No Iterations 0
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 3.94984e-11, global = -1.08582e-13, cumulative = -1.08219e-13
[...skipped some for character limit]]
DILUPBiCG: Solving for omega, Initial residual = 9.8383e-06, Final residual = 2.19925e-07, No Iterations 1
bounding omega, min: -328.374 max: 95813.5 average: 225.052
DILUPBiCG: Solving for k, Initial residual = 0.0024148, Final residual = 3.30216e-07, No Iterations 4
ExecutionTime = 335.97 s ClockTime = 336 s
Courant Number mean: 0.00278018 max: 77.708
Time = 0.0003
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for Ux, Initial residual = 0.00294713, Final residual = 5.09134e-07, No Iterations 5
DILUPBiCG: Solving for Uy, Initial residual = 0.0498396, Final residual = 1.30056e-07, No Iterations 9
DILUPBiCG: Solving for Uz, Initial residual = 0.245217, Final residual = 8.18823e-07, No Iterations 8
DILUPBiCG: Solving for h, Initial residual = 0.999938, Final residual = 5.25371e-07, No Iterations 3
DICPCG: Solving for p_rgh, Initial residual = 0.869755, Final residual = 0.00846313, No Iterations 487
DICPCG: Solving for p_rgh, Initial residual = 5.09343e-05, Final residual = 4.87965e-07, No Iterations 88
DICPCG: Solving for p_rgh, Initial residual = 4.87965e-07, Final residual = 9.50119e-09, No Iterations 78
DICPCG: Solving for p_rgh, Initial residual = 9.50298e-09, Final residual = 9.50298e-09, No Iterations 0
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
[...]
time step continuity errors : sum local = 1.18932e-09, global = 4.01695e-11, cumulative = 4.48692e-11
DICPCG: Solving for p_rgh, Initial residual = 0.0262778, Final residual = 0.000244534, No Iterations 513
DICPCG: Solving for p_rgh, Initial residual = 0.000155778, Final residual = 1.50242e-06, No Iterations 79
DICPCG: Solving for p_rgh, Initial residual = 1.50243e-06, Final residual = 1.48662e-08, No Iterations 215
DICPCG: Solving for p_rgh, Initial residual = 1.48865e-08, Final residual = 9.91452e-09, No Iterations 4
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 1.14702e-09, global = 1.97249e-11, cumulative = 6.45941e-11
DILUPBiCG: Solving for omega, Initial residual = 0.999996, Final residual = 4.51237e-08, No Iterations 8
bounding omega, min: -154.316 max: 1.68053e+13 average: 4.75822e+08
DILUPBiCG: Solving for k, Initial residual = 0.965169, Final residual = 5.95112e-07, No Iterations 7
ExecutionTime = 508.91 s ClockTime = 510 s
Courant Number mean: 0.3024 max: 2550.14
Time = 0.0004
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for Ux, Initial residual = 0.0150419, Final residual = 5.99711e-07, No Iterations 11
DILUPBiCG: Solving for Uy, Initial residual = 0.236398, Final residual = 4.35313e-07, No Iterations 15
DILUPBiCG: Solving for Uz, Initial residual = 0.312253, Final residual = 2.32705e-07, No Iterations 15
DILUPBiCG: Solving for h, Initial residual = 0.99989, Final residual = 6.77965e-07, No Iterations 6
From function Foam::scalar Foam::species::thermo<Thermo, Type>::T(Foam::scalar, Foam::scalar, Foam::scalar, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar, Foam::scalar) const, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar, Foam::scalar) const, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar) const, bool) const [with Thermo = Foam::hPolynomialThermo<Foam::icoPolynomial<Foam::specie> >; Type = Foam::sensibleEnthalpy; Foam::scalar = double; Foam::species::thermo<Thermo, Type> = Foam::species::thermo<Foam::hPolynomialThermo<Foam::icoPolynomial<Foam::specie> >, Foam::sensibleEnthalpy>]
in file /home/ubuntu/OpenFOAM/OpenFOAM-dev/src/thermophysicalModels/specie/lnInclude/thermoI.H at line 70
Energy -> temperature conversion failed to converge:
[0] iter Test e/h Cv/p Tnew
[0] 0 316.267 2.73901e+06 4200 -341.927
[0] 1 -341.927 -3.12693e+06 4200 396.53
[0] 2 396.53 3.20771e+06 4200 -373.26
[0] 3 -373.26 -3.19648e+06 4200 381.757
[0] 4 381.757 3.11801e+06 4200 -366.676
[0] 5 -366.676 -3.18078e+06 4200 384.601
[0] 6 384.601 3.13528e+06 4200 -367.945
[0] 7 -367.945 -3.18376e+06 4200 384.043
[char limit saving]
[0] 100 384.135 3.13245e+06 4200 -367.737
[0] 101 -367.737 -3.18327e+06 4200 384.135
[1] iter Test e/h Cv/p Tnew
[1] 0 317.791 1.11566e+07 4200 -1245.64
[1] 1 -1245.64 -6.61605e+06 4200 1422.51
[1] 2 1422.51 5.20461e+06 4200 1276.22
[1] 3 1276.22 4.85664e+06 4200 1212.77
[...]
[1] 98 1104.47 4.7527e+06 4200 1065.77
[1] 99 1065.77 4.81153e+06 4200 1013.06
[1] 100 1013.06 4.96851e+06 4200 922.981
[1] 101 922.981 5.52505e+06 4200 700.39
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 0 in communicator MPI COMMUNICATOR 3 SPLIT FROM 0
with errorcode 1.
Problem
========================
No clue on how to make anything work. Are my BCs ok? Is my thermophysicalProperties entry any wrong? Should i change anything in fv* dictionaries? Is my mesh to coarse maybe? It's already several milion cells!
========================
Thermophysical
Code:
thermoType
{
type heRhoThermo;
mixture pureMixture;
transport polynomial;
thermo hPolynomial;
equationOfState icoPolynomial;
specie specie;
energy sensibleEnthalpy;
}
pRef 100000;
//dpdt off;
mixture
{
specie
{
molWeight 18;
}
equationOfState
{
rhoCoeffs<8> (117.045534 7.913913 -0.02254739 1.97987E-05 0 0 0 0);
}
thermodynamics
{
Hf 0;
Sf 0;
CpCoeffs<8> (4200 0 0 0 0 0 0 0);
}
transport
{
muCoeffs<8> (0.1368992 -0.0012040579 3.5597468E-06 -3.526231E-09 0 0 0 0);
kappaCoeffs<8> (-2.4283206 0.0233779 -5.9953665E-05 5.2602452E-08 0 0 0 0);
}
}
// ************************************************************************* //
p_rgh
Code:
internalField uniform 1E+5;
boundaryField
{
Inlet
{
type fixedFluxPressure;
gradient uniform 0;
value $internalField;
}
Outlet
{
type fixedValue;
value $internalField;
}
Mantle
{
type fixedFluxPressure;
gradient uniform 0;
value $internalField;
}
Connectors
{
type fixedFluxPressure;
gradient uniform 0;
value $internalField;
}
SymmetryY
{
type symmetry;
}
}
p
Code:
internalField uniform 1E+5; // initialize with funkySetFields
boundaryField
{
Inlet
{
type calculated;
value $internalField;
}
Outlet
{
type calculated;
value $internalField;
}
Mantle
{
type calculated;
value $internalField;
}
Connectors
{
type calculated;
value $internalField;
}
SymmetryY
{
type symmetry;
}
}
velocity
Code:
internalField uniform (0 0 0);
boundaryField
{
Inlet
{
type groovyBC;
//valueExpression "vector(0.015*pow((1-(sqrt(pow(pos().z,2)+pow(pos().y,2))/0.0115)),0.14),0,0)";
valueExpression "vector(0.15*pow((1-(sqrt(pow(pos().z,2)+pow(pos().y,2))/0.0115)),0.14),0,0)";
value uniform (0.12 0 0);
}
Outlet
{
type zeroGradient;
}
Mantle
{
type noSlip;
}
Connectors
{
type noSlip;
}
SymmetryY
{
type symmetry;
}
}
omega
Code:
internalField uniform 0.937;
boundaryField
{
Inlet
{
type fixedValue;
value uniform 0.937;
}
Outlet
{
type zeroGradient;
}
Mantle
{
type omegaWallFunction;
value uniform 0.937;
}
Connectors
{
type omegaWallFunction;
value uniform 0.937;
}
SymmetryY
{
type symmetry;
}
}
k
Code:
internalField uniform 8.44E-5;
boundaryField
{
Inlet
{
type fixedValue;
value $internalField;
}
Outlet
{
type zeroGradient;
}
Mantle
{
type kLowReWallFunction;// kqRWallFunction;
value $internalField;
}
Connectors
{
type kLowReWallFunction;// kqRWallFunction;
value $internalField;
}
SymmetryY
{
type symmetry;
}
}
nut
Code:
internalField uniform 0;
boundaryField
{
Inlet
{
type calculated;
value $internalField;
}
Outlet
{
type calculated;
value $internalField;
}
Mantle
{
type nutLowReWallFunction; // nutkWallFunction; // nutLowReWallFunction
value uniform 0;
}
Connectors
{
type nutLowReWallFunction; // nutkWallFunction;
value uniform 0;
}
SymmetryY
{
type symmetry;
}
}
T
Code:
internalField uniform 315.15;
boundaryField
{
Inlet
{
type fixedValue;
value uniform 293.15;
}
Outlet
{
type zeroGradient; //inletOutlet;
//inletValue 315.15 ;//$internalField;
//value $internalField;
}
Mantle
{
type zeroGradient;
}
Connectors
{
type zeroGradient;
}
SymmetryY
{
type symmetry;
}
}
alphat
Code:
internalField uniform 0;
boundaryField
{
Inlet
{
type calculated;
value uniform 0;
}
Outlet
{
type calculated;
value uniform 0;
}
Mantle
{
type compressible::alphatJayatillekeWallFunction;
Prt 0.85;
value uniform 0;
}
Connectors
{
type compressible::alphatJayatillekeWallFunction;
Prt 0.85;
value uniform 0;
}
SymmetryY
{
type symmetry;
}
}
fvSolution
Code:
solvers
{
cellDisplacement // OpenFOAM+ only
{
solver GAMG;
smoother GaussSeidel;
tolerance 1E-07;
relTol 0.01;
}
"rho.*"
{
solver PCG;
preconditioner DIC;
tolerance 0;
relTol 0;
}
p_rgh
{
solver PCG;
preconditioner DIC;
tolerance 1e-8;
relTol 0.01;
maxIter 2000; // non-default
}
p_rghFinal
{
$p_rgh;
relTol 0;
}
omega
{
solver PBiCGStab; // PBiCGStab;
preconditioner DILU;
tolerance 1e-7;
relTol 0.1;
}
h
{
solver PBiCGStab; // PBiCGStab;
preconditioner DILU;
tolerance 1e-7;
relTol 0.01;
}
"(U|h|e|k|R)"
{
solver PBiCG; // PBiCGStab;
preconditioner DILU;
tolerance 1e-6;
relTol 0.1;
maxIter 1500;
}
"(U|h|e|k|omega|R)Final"
{
$U;
relTol 0;
}
}
PIMPLE
{
momentumPredictor yes;
nOuterCorrectors 1; //3;
nCorrectors 4; //3
nNonOrthogonalCorrectors 1; // 0;
pRefCell 0;
pRefValue 1e5;
}
relaxationFactors
{
fields
{
p_rgh 0.6;
p 0.6;
rho 0.6;
}
equations
{
"U.*" 0.7;
"rho.*" 0.7;
"p_rgh.*" 0.7;
"epsilon" 0.3;
"omega" 0.3;
"h" 0.5;
"k" 0.6;
}
}
fvSchemes
Code:
ddtSchemes
{
default Euler;
}
gradSchemes
{
default Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss linearUpwindV grad(U);
div(phi,h) Gauss linearUpwind grad(U);
div(phi,e) Gauss linearUpwind grad(U);
div(phi,k) Gauss linearUpwind grad(omega);
div(phi,omega) Gauss linearUpwind grad(omega);
div(phi,R) Gauss linearUpwind grad(U);
div(phi,K) Gauss linearUpwind grad(U);
div(phi,Ekp) Gauss linearUpwind grad(U);
div(R) Gauss linearUpwind grad(U);
div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}
laplacianSchemes
{
default Gauss linear uncorrected;
laplacian(1,p) Gauss linear corrected;
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default corrected;
}
wallDist
{
method meshWave;
}
fluxRequired
{
default no;
p_rgh;
}
// ************************************************************************* //
|
