 |
 |
 |
|
|
|
[Sponsors] | ||||
September 27, 2019, 03:25
|
Highly unstable buoyantPimpleFoam simulation
|
#1 |
|
Member
Piotr Ładyński
Join Date: Apr 2017
Posts: 55
Rep Power: 9  |
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.
======================== 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);
}
}
// ************************************************************************* //
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;
}
}
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;
}
}
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;
}
}
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;
}
}
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;
}
}
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;
}
}
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;
}
}
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;
}
}
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;
}
}
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;
}
// ************************************************************************* //
|
|
|
|
|
|
May 20, 2020, 07:38
|
|
#2 |
|
New Member
Khaled Yassin
Join Date: Aug 2019
Location: Jülich-Germany
Posts: 12
Rep Power: 6 |
Hi piotr.mecht,
Have you found the soultion of this problem??
__________________
-- Best Regards, Khaled Yassin, Research Assistant Institute for Energy and Environmental Research (IEK-14) Forschungszentrum Jülich |
|
|
|
|
|
|
May 23, 2020, 09:15
|
|
#3 |
|
Member
Piotr Ładyński
Join Date: Apr 2017
Posts: 55
Rep Power: 9 |
Yes!
As the fluid i tried to simulate (water) was changing its density only due to temperature changes, the solver was having hard time trying to calculate pressure-work term. From OpenFOAM-v1906 source code energy equation looks like: Code:
fvScalarMatrix EEqn
(
fvm::ddt(rho, he) + fvm::div(phi, he)
+ fvc::ddt(rho, K) + fvc::div(phi, K)
+ (
he.name() == "e"
? fvc::div
(
fvc::absolute(phi/fvc::interpolate(rho), U),
p,
"div(phiv,p)"
)
: -dpdt
)
- fvm::laplacian(turbulence->alphaEff(), he)
==
rho*(U&g)
+ radiation->Sh(thermo, he)
+ fvOptions(rho, he)
);
Code:
dpdt off; I remember, that buoyantPimpleFoam worked somehow differently in openfoam 7 and OpenFOAM-v1906, and that particular case was crashing in OF7. My thermophysicalProperties dictionary looked like this: Code:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1906 |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object thermophysicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
thermoType
{
type heRhoThermo;
mixture pureMixture;
transport polynomial;
thermo hPolynomial;
equationOfState icoPolynomial;
specie specie;
energy sensibleEnthalpy;
}
pRef 1E+5;
dpdt off;
mixture
{
specie
{
nMoles 1;
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);
}
}
|
|
|
|
|
|
|
May 23, 2020, 09:33
|
|
#4 |
|
New Member
Khaled Yassin
Join Date: Aug 2019
Location: Jülich-Germany
Posts: 12
Rep Power: 6 |
Thank you very much for your answer....I will try that
__________________
-- Best Regards, Khaled Yassin, Research Assistant Institute for Energy and Environmental Research (IEK-14) Forschungszentrum Jülich |
|
|
|
|
|
 |
|
|
Similar Threads
|
||||
| Thread | Thread Starter | Forum | Replies | Last Post |
| Some questions about flow boiling simulation in Fluent | beastieboys6 | FLUENT | 8 | November 20, 2017 23:47 |
| Simulation FPEs - turbulence for transient and steady-state? | DaveR | OpenFOAM Running, Solving & CFD | 5 | March 5, 2017 15:06 |
| Unstable LES simulation - Courant going higher throughout time - Flow around a sphere | beluiz93 | OpenFOAM Running, Solving & CFD | 11 | November 3, 2016 06:05 |
| AMI simulation unstable within first few timesteps | craig0a | OpenFOAM Running, Solving & CFD | 0 | February 24, 2014 12:09 |
| 3-D Contaminant Dispersal Simulation | Apple L S Chan | Main CFD Forum | 1 | December 23, 1998 10:06 |
1Likes
Linear Mode
