[Thomas Baumann]

Dear all,

I make simulations with chtmultiRegionSimpleFoam and with chtmultiRegionFoam (OF 5.0). It is a turbulent channel flow with a row of heated rods (heat source in fvOptions in the rods). The inlet temperature of the air has already ca. 1000 K. The temperature of the fluid by the rods is increased by ca. 30K

Around the channel is additional insulation with low kappa values (0.1) and a lower ambient temperature ca. 300K).
So in total I have 3 regions (fluid channel, Solid_Heater_Rods, Solid_Insulation).

Everything works fine without radiation, but the rods are of course getting too hot.

Adding radiation (FVDOM) effects the simulation crashes after some iterations due to negative temperatures (always directly after calculating the ILambda values) in the fluid and in the insulation solid. It starts that the temperatures are getting below the defined bc values and then are getting negative.


DILUPBiCGStab: Solving for ILambda_223_0, Initial residual = 0.0111591837, Final residual = 3.29815797e-005, No Iterations 2
Radiation solver iter: 1
DILUPBiCGStab: Solving for ILambda_56_0, Initial residual = 0.00921928238, Final residual = 1.80269034e-005, No Iterations 2
DILUPBiCGStab: Solving for ILambda_57_0, Initial residual = 0.00933455961, Final residual = 1.65074399e-005, No Iterations 2
DILUPBiCGStab: Solving for ILambda_82_0, Initial residual = 0.00924826846, Final residual = 1.77962475e-005, No Iterations 2
DILUPBiCGStab: Solving for ILambda_83_0, Initial residual = 0.0091733533, Final residual = 1.86632476e-005, No Iterations 2
DILUPBiCGStab: Solving for ILambda_84_0, Initial residual = 0.00868638911, Final residual = 2.44328267e-005, No Iterations 2
DILUPBiCGStab: Solving for ILambda_85_0, Initial residual = 0.00869436019, Final residual = 1.89730532e-005, No Iterations 2
DILUPBiCGStab: Solving for ILambda_86_0, Initial residual = 0.00822859442, Final residual = 1.65260782e-005, No Iterations 2
DILUPBiCGStab: Solving for ILambda_87_0, Initial residual = 0.00818240314, Final residual = 1.54362665e-005, No Iterations 2
DILUPBiCGStab: Solving for ILambda_109_0, Initial residual = 0.00806299195, Final residual = 1.76092394e-005, No Iterations 2
DILUPBiCGStab: Solving for ILambda_110_0, Initial residual = 0.00856531701, Final residual = 2.0917397e-005, No Iterations 2
DILUPBiCGStab: Solving for ILambda_111_0, Initial residual = 0.0086226024, Final residual = 2.41867149e-005, No Iterations 2
DILUPBiCGStab: Solving for ILambda_112_0, Initial residual = 0.00926088609, Final residual = 2.32284362e-005, No Iterations 2
DILUPBiCGStab: Solving for ILambda_113_0, Initial residual = 0.00875198209, Final residual = 2.10845819e-005, No Iterations 2
DILUPBiCGStab: Solving for ILambda_114_0, Initial residual = 0.00829104384, Final residual = 1.78349556e-005, No Iterations 2
DILUPBiCGStab: Solving for ILambda_137_0, Initial residual = 0.00811726417, Final residual = 7.63348464e-006, No Iterations 2
DILUPBiCGStab: Solving for ILambda_138_0, Initial residual = 0.00866982572, Final residual = 1.87350822e-005, No Iterations 2
DILUPBiCGStab: Solving for ILambda_139_0, Initial residual = 0.00868832268, Final residual = 1.98727924e-005, No Iterations 2
Min/max T:-1188.42124 2500.04226
GAMG: Solving for p_rgh, Initial residual = 0.00495497568, Final residual = 2.8333931e-005, No Iterations 3
time step continuity errors : sum local = 0.000158236471, global = -6.61631456e-006, cumulative = 0.000945905867
Min/max rho:0.1 1.55286807
DILUPBiCGStab: Solving for epsilon, Initial residual = 0.000381979405, Final residual = 9.9774055e-008, No Iterations 2
DILUPBiCGStab: Solving for k, Initial residual = 0.000950370334, Final residual = 3.44609781e-007, No Iterations 2

Solving for solid region Solid_Insulation
GAMG: Solving for h, Initial residual = 0.0101350153, Final residual = 9.61943473e-009, No Iterations 2
GAMG: Solving for h, Initial residual = 0.00914097844, Final residual = 8.0324505e-009, No Iterations 2
GAMG: Solving for h, Initial residual = 0.00848300327, Final residual = 6.80844891e-009, No Iterations 2
Min/max T:-1368.33306 2347.67783

Solving for solid region Solid_HeaterZone
GAMG: Solving for h, Initial residual = 0.000345150206, Final residual = 8.77786335e-011, No Iterations 2
GAMG: Solving for h, Initial residual = 0.000332840555, Final residual = 3.80969835e-011, No Iterations 2
GAMG: Solving for h, Initial residual = 0.00032599808, Final residual = 2.37158585e-011, No Iterations 2
Min/max T:2407.78594 2515.43428
ExecutionTime = 440.985 s ClockTime = 441 s




I added temperature limiters in both zones, I tried a lot with the relaxation parameters without any success.

When I increase the kappa of the insulation to e.g. 6 it works again but then of course the heat transfer is much too high in the insulation region.

It also works when I reduce the inlet temperature of the air flow or reduce the heat source, or reduce the paramters of the radiation properties.

So I guess the bigger the radiation effects and the less thermal diffusion/capacity in the solid zones the more unstable is the simulation.


Here the radiationPorperties of the insulation:
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object radiationProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //


radiationModel opaqueSolid;

absorptionEmissionModel constantAbsorptionEmission;

constantAbsorptionEmissionCoeffs
{
absorptivity absorptivity [ 0 -1 0 0 0 0 0 ] 0.7;
emissivity emissivity [ 0 -1 0 0 0 0 0 ] 0.7;
E E [ 1 -1 -3 0 0 0 0 ] 0;
}

scatterModel none;





Here the thermalproperties of the insulation:
*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.1.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object solidThermophysicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

thermoType
{
type heSolidThermo;
mixture pureMixture;
transport constIso;
thermo hConst;
equationOfState rhoConst;
specie specie;
energy sensibleEnthalpy;
}

mixture
{
specie
{
nMoles 1;
molWeight 60;
}

transport
{
kappa 0.2;// 6;//120;
}

thermodynamics
{
Hf 0;
Cp 800;
}

equationOfState
{
rho 200;
}
}


Here are the radiationProperties of the Fluid Region:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.1.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object radiationProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

radiation on;

radiationModel fvDOM;

fvDOMCoeffs
{
nPhi 7;// 3; // polar angles in PI (from Z to X-Y plane)
nTheta 8;//5; // azimuthal angles in PI/2 on X-Y.(from Y to X)
tolerance 1e-2; // 1e-2convergence criteria for radiation iteration
maxIter 6; // maximum number of iterations
}

// Number of flow iterations per radiation iteration
solverFreq 10;

radiationModel fvDOM;


// Number of flow iterations per radiation iteration
solverFreq 10;

absorptionEmissionModel constantAbsorptionEmission; // constantAbsorptionEmission;

constantAbsorptionEmissionCoeffs
{
absorptivity absorptivity [ m^-1 ] 0.001;
emissivity emissivity [ m^-1 ] 0.001;
E E [ kg m^-1 s^-3 ] 0;
}

scatterModel none;//constantScatter;

constantScatterCoeffs
{
sigma sigma [ 0 -1 0 0 0 0 0 ] 0;
C C [ 0 0 0 0 0 0 0 ] 0;
}

sootModel none;

transmissivityModel none;




fvSchemes Fluid Region:
*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.2.2 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
default steadyState;
}

gradSchemes
{
default cellLimited Gauss linear 1;//Gauss linear;
grad(U) cellLimited Gauss linear 1;
}

divSchemes
{
default bounded Gauss upwind;
div(phi,U) bounded Gauss upwind;
div(phi,K) bounded Gauss upwind;
div(phi,h) bounded Gauss upwind;
div(phi,k) bounded Gauss upwind;
div(phi,K) bounded Gauss upwind;
div(phi,epsilon) bounded Gauss upwind;
div(phi,R) bounded Gauss upwind;
div(R) Gauss linear;
div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
default Gauss linear corrected;//Gauss linear corrected 0.33;
}

interpolationSchemes
{
default linear;
}

snGradSchemes
{
default corrected;
}
/*
fluxRequired
{
default no;
p_rgh;
}*/




Has anybody an idea what is not correct?
Are the values to st the absorbtitivy and the emissivity constant for the solids correct tor a typical grey radiator?

I think it could be effected by local too high emissivity effects where small regions have a much too high heat flux by radiation and thus by the small kappa value of the insulation this regions are getting negative.

Thanks a lot, Thomas

[Alczem]

Hey!


I also have had trouble with chtMultiRegionSimpleFoam and radiation at the same time, with the same issue of negative temperatures.


I would suggest to try to solve radiation every flow iteration (solverFreq 1). It is more costly, but I feel like it helps the enery equation to remain more stable. When you solve the radiation every 10 iterations, I think it computes a radiative heat loss that remains constant for 10 iterations before being computed again. During those 10 iterations, the heat loss remains constant and unaffected by the other heat transfers by convection and conduction, so that could lead to really low temperatures. My understanding might be flawed though!


The other thing is using relaxation either for the energy and radiation equations/fields or both. If you are using the BC externalWallHeatFluxTemperature, you can also relax temperature specifically for this BC.


I would like to know if you manage to find a stable setup as I am still struggling with this too

[Thomas Baumann]

Hi Alczem,

I tried your advices. Changing the radiation frequency to 1 does not help. The simulation crushes again.

I changed also the relaxation paramters to low values or I changed so that the solid h relaxation is much lower than the fluid relaxation and also in opposite direction.

I also increased the fvDom Coeffs to much higher values, so that you don't have these peaks in the radiation fluxes. It seems that more numbers in the Coeffs can help and additional more cells around the rods help a bit. But the simulations are still crushing.

I changed also radiation model to viewFactor. Here the radiation directions are only calculated at the beginning and so there are no changes of the direction. This simulation also crushes, too.

I am using the externalWallHeatfluxBoundary condition. I relax here the qr. But I have not seen how to relax the temperature directly at the bc condition here. How can I do this?

But in general no externalWallHeatfluxBoundary.
I use from solid to bc-conditions the externalWallHeatFluxTemperature. Here no radiation should have influence:
"Wall.*|wall.*"
{
type externalWallHeatFluxTemperature;

mode coefficient;

Ta uniform 300;
h uniform 5.0;
thicknessLayers ();
kappaLayers ();

kappaMethod solidThermo;

value $internalField;
// qr qr;
// qrRelaxation 0.05;
// value uniform 600;
}


For bc conditions of the Fluid-part I relax the qr. But in this case I don't have such a bc-condition. because the fluid is connected only with solid regions.

Wall_channel
{
type externalWallHeatFluxTemperature;

mode coefficient;

Ta uniform 298.15;
h uniform 5.0;
thicknessLayers (0.3);
kappaLayers (0.1);

kappaMethod fluidThermo;

value $internalField;
qr qr;
qrRelaxation 0.05;
/*value uniform 600;*/
}

[Alczem]

Hey!


Too bad it is not helping... You can relax the temperature for externalWallHeatFluxTemperature simply by using the field relaxation (might not be available for the foundation version, I usually work with the ESI version of OpenFOAM).


But I remember that externalWallHeatFluxTemperature gave me headaches when used with radiation Sorry I could not be more helpful! Do you mind sharing the fvSolution and fvSchemes just in case?

[Thomas Baumann]

Hello,

sorry for the late response. I have no externalwallheatflux bc direct at the fluid fuild, only at the bc condition for the solid field to the ambient and here are not the problems. It would be very helpful to relax the coupling bc

Here is my fvSolution from the fluid field:


/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.1.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
/* rho
{
solver PCG
preconditioner DIC;
tolerance 1e-7;
relTol 0;
}

p_rgh
{
solver GAMG;
tolerance 1e-7;
relTol 0.01;

smoother GaussSeidel;

cacheAgglomeration true;
nCellsInCoarsestLevel 10;
agglomerator faceAreaPair;
mergeLevels 1;

maxIter 10;
}

"(U|h|k|epsilon)"
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-7;
relTol 0.0001;
}*/

"p_rgh.*"
{
solver GAMG;
smoother symGaussSeidel;
tolerance 1e-7;
relTol 0.01;
}

"(U).*"
{
solver PBiCGStab;
preconditioner DILU;
tolerance 1e-8;
relTol 0.0001;
}
"(h|k|epsilon).*"
{
solver PBiCGStab;
preconditioner DILU;
tolerance 1e-7;
relTol 0.001;
}

G
{
solver PCG;
preconditioner DIC;
tolerance 1e-7;
relTol 0.005;
}

GFinal
{
$G;
relTol 0;
}


Ii
{
solver PBiCGStab;
preconditioner DILU;
tolerance 1e-7;
relTol 0.005;
minIter 1;
}
}

/*
SIMPLE
{
momentumPredictor on;
nNonOrthogonalCorrectors 3;

pMinFactor 0.4;
pMaxFactor 3.0;
pRefCell 0;
pRefValue 0;
rhoMin rhoMin [1 -3 0 0 0] 0.1;
rhoMax rhoMax [1 -3 0 0 0] 3;
residualControl
{
p 1e-7;
U 1e-7;
e 1e-7;

// Possibly check turbulence fields
"(k|epsilon|omega)" 1e-7;
}
}*/


SIMPLE
{
momentumPredictor on;
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
rhoMin rhoMin [1 -3 0 0 0] 0.1;
rhoMax rhoMax [1 -3 0 0 0] 3;
}



relaxationFactors
{
fields
{
rho 1.0;
p_rgh 0.35; //start with reduced values
}
equations
{
U 0.35;////start wth reduced values (like 0.35)
h 0.755; //start wth reduced values (like 0.5)
e 0.5;
"(k|epsilon|omega)" 0.4;
// G 0.37;
"ILambda.*" 0.57;
// qr 0.57;
}
}



// ************************************************** *********************** //





for the solid:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.1.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
"e.*"
{
solver GAMG;
smoother symGaussSeidel;
tolerance 1e-7;
relTol 0.001;
}
"h.*"
{
solver GAMG;
smoother symGaussSeidel;
tolerance 1e-7;
relTol 0.001;
minIter 2;
}
}

SIMPLE
{
nNonOrthogonalCorrectors 1;
}

relaxationFactors
{
fields
{
}
equations
{
h 0.5;
// e 0.5;
// ".*" 0.6;
}
}
I also tested different combinations of the h field relaxations. (e.g. fluid field 0.7 and for the solid field 0.3.




If found also an iteresting post:
Heat transfer convergence problem

I will increase the number of cells near the rods (maybe use also less rods) to check the effect.
Maybe I will also switch to a newer version of OF to see if it is getting better hereby.

Regards Thomas

[dlahaye]

@Thomas:

what is your solverFreq in constant/radiationProperties ?

Typically it is set to 10. This implies that radiation is solved every 10 flow steps. For harder problem (larger radiative heart losses) it might (emphasis on "might") make sense to update the radiative field more often, i.e., reduce solverFreq to e.g. 5, 2 or 1.

Good luck.

[Alczem]

Originally Posted by dlahaye

@Thomas:

what is your solverFreq in constant/radiationProperties ?

Typically it is set to 10. This implies that radiation is solved every 10 flow steps. For harder problem (larger radiative heart losses) it might (emphasis on "might") make sense to update the radiative field more often, i.e., reduce solverFreq to e.g. 5, 2 or 1.

I had the same reasoning and suggested this fix too, but it had no effect apparently.


Thomas, if you have not already, can you try setting the absorptivity and emissivity to something like 0.5, reduce the angular discretization to 2 and 1 for phi and theta, set the tolerance to 1e-4 and maxIter to 1 (with solverFreq set to 1)? I am just sharing the most generic and robust settings I found during some tests setting the absorptivity and emissivity to a really small value or to 0 when using constantAbsorptionEmission gave me a lot of issues with temperature.


Maybe try the P1 model to see if it has to do with fvDOM?

[dlahaye]

I imagine that comparing with P1 might indeed be valuable here.

In case that radiative heat transfer is important, I would suggest using many rays (i.e. large values for phi and theta), solve the rays accurately (maxit =1000, atol = 1e-8, rtol = 1e-8), and solve the rays frequently (solverFreq = 1, 2 or 5).

To have grip on the convergence of the solver, I suggest to monitor min(T) and max(T) during the convergence of the solver (see below). In case that radiative heat transfer is important, min(T) and max(T) should increase/decrease fast and faster in case that solverFreq is low.

Good luck.

[Thomas Baumann]

I tested in the last days several modifications. I can tell you especially what does not help

I changed the frequency of the radiation to 1. The simulation runs more iterations, but then crashes, too. The numerical effort is enorm with such a change.

The simulations work if I reduce the absorbtivity and emissivity too lower values (like 0.1) or if I increase the kappa value of the insulation.

When I change nphi, theta ez. to lower values. Here the temperatures are getting very high and not realistic (1900K of instead 1600K). I think this is related due to the very coarse direction of of the radiation.
I changed also to much higher values nhi & theta to 10, but here the simuations crashes always, too.

It seems that with a much finer grid the simulation is more stable, but it still crashes with low kappa values and higher emissivity & absorbtivity. I test now cases with crashed with low grid discretication with ca 300 K cells with a much finer grid 4.3 mio cells to see, if it really helps. But the simulations are still running.

I changed all bc conditions to have no externalWallHeat bc's anymore, but this didn't help.

I reduced the relaxation parameters in the fvSoultion, which helpes for stabilization only a bit.

I also tested the simulation unsteady with the solver chtMultiRegionFoam, but it also crashed.


Here are the results of the temperture field for the carse grid with different scaling.
Simulation_with_Coarse_Mesh.jpg
Simulation_with_Coarse_Mesh_other_Scaling.jpg
You cn see that the temperature increases by ca. 7 Kelvin by flowing through the rods.
But there is also a low temperature field at the insulation dirct before the rods.
The same effect exists with the much finer mesh, but with much more cold spots of the rod.Simulation_with_Rods_Fine_Mesh.jpg

Could it be effected by the inlet conditions in the radiation field?


Next planned steps:
I will also test the P1 model, maybe with additional scattering.
I will test also a newer OF-Version.
I will change also the thermophyiscalProperties to another setup.