CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Running, Solving & CFD (
-   -   can t reach steady state chtMultiregionFaom (

mabinty May 7, 2009 14:21

can t reach steady state chtMultiregionFaom
Dear all!

I m about to simulate half of a symmetric 2D 1x1m2 domain, with hot air entering from an inlet in the floor. The sides are open (ambient air can enter). My interest lays in the heat transfer to the ceiling at steady state. Unfortunately I m not able to reach a steady state as the velocities in the domain constantly increase over simulation time (total time 400s). Another strange fact is a min temperature in the air region slightly below ambient temperature (see log data). I observed that the increase in velocity is smaller with finer mesh (hex mesh) but still there. Below some specifications of my set-up are posted.

I would greatly appreciate your comments. Thanks in advance!

--- discretization schemes and solver settings ---

taken from default settings

--- BC for the air region ---

inlet: fixedValue uniform (0 0.8 0);
floor: fixedValue uniform (0 0 0);
air_to_ceiling: fixedValue uniform (0 0 0);
minX: symmetryPlane;
maxX: pressureDirectedInletOutletVelocity;
inletDirection uniform (1 0 0);
value uniform (0 0 0);

internalField: uniform 293;
inlet: fixedValue uniform 473;
floor: zeroGradient
air_to_ceiling: solidWallMixedTemperatureCoupled;
neighbourRegionName ceiling;
neighbourPatchName ceiling_to_air;
neighbourFieldName T;
K K;
value uniform 293;
minX: symmetryPlane;
maxX: inletOutlet;
inletValue uniform 293;
value uniform 293;

_epsilon (for standard k-eps): initial value 0.002;
inlet: fixedValue uniform 0.002;
zeroGradient for all other boundaries

_k (for standard k-eps): initial value 0.01;
inlet: fixedValue uniform 0.01;
zeroGradient for all other boundaries

maxX: zeroGradient;
fixedFluxBuoyantPressure uniform 0; for all other boundaries

initial value 100000
inlet: zeroGradient
floor and air_to_ceiling: calculated
minX: symmetryPlane;
maxX: fixedValue uniform 100000;

--- log data from the last time step ---

Solving for fluid region air
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for h, Initial residual = 1.271221e-06, Final residual = 6.479811e-08, No Iterations 1
Min/max T:min(T) [0 0 0 1 0 0 0] 292.4249 max(T) [0 0 0 1 0 0 0] 473
GAMG: Solving for pd, Initial residual = 0.3610631, Final residual = 0.0001952987, No Iterations 1
GAMG: Solving for pd, Initial residual = 0.0002406985, Final residual = 9.19176e-06, No Iterations 2
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors (air): sum local = 4.091744e-10, global = 1.526135e-10, cumulative = 7.550246e-05
GAMG: Solving for pd, Initial residual = 0.3604305, Final residual = 0.0001962409, No Iterations 1
GAMG: Solving for pd, Initial residual = 0.0002416951, Final residual = 5.678556e-07, No Iterations 6
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors (air): sum local = 2.535464e-11, global = -8.27864e-12, cumulative = 7.550245e-05
DILUPBiCG: Solving for epsilon, Initial residual = 5.23003e-06, Final residual = 1.120886e-08, No Iterations 2
DILUPBiCG: Solving for k, Initial residual = 6.230139e-06, Final residual = 2.278377e-08, No Iterations 2

Solving for solid region ceiling
DICPCG: Solving for T, Initial residual = 5.053041e-07, Final residual = 5.053041e-07, No Iterations 0
DICPCG: Solving for T, Initial residual = 5.05303e-07, Final residual = 5.05303e-07, No Iterations 0
Min/max T:min(T) [0 0 0 1 0 0 0] 293 max(T) [0 0 0 1 0 0 0] 295.2498
ExecutionTime = 13231.56 s ClockTime = 13290 s


mabinty May 11, 2009 12:38


I tried to change different parameters like BC for the open side (maxX) as well as parameters in the fvSchemes/fvSolutions. What I can observe most of the times is a temperature drop below ambient which seems not to be realistic to me. Hmm, I couldn t figure out yet whats the problem but keep on investigating and appreciate your comments!!

Thx in advance,

psosnows May 12, 2009 06:40

hello Aram,

I am working on a bit similar problem using a bit modified chtMultiRegionFoam.

When I was using the boundaries proposed in chtMultiRegionFoam I also got some strange temperatures which were lower than everything around. I thought that this was due to some numerical errors (since those drops were quite small).

Right now I am not so sure about it and in fact have some doubt in the solidWallTemperatureCoupled and solidWallHeatFluxTemperatureCoupled boundary conditions.
The problem may be with the K parameter which dimension does not mach for both boundaries.

try to add line:

Info << "K dim in swhftc " << K.dimensions() << endl;
inside updateCoeffs() function of solidWallHeatFluxTemperatureCoupled

and line:

Info << "Kw dim in swtc " << Kw.dimensions() << endl;
inside flux() function of solidWallTemperatureCoupled

in the tutorial case the dimensions do not mach...

I asked about it in thread:
unfortunately did not get any response.

in the end I couple the boundaries by putting on the boundary of the region a fixed value of the temperature inside the cell lying in the neighbour region (it took some effort but was worth it).
I do it for both sides of the multi region wall, get no strange temperature drops etc and the simulation seems to work fine.


mabinty May 12, 2009 15:05

Hi Pawel!

Thanks a lot for your remark on the cht-BC, in special the dimension of conductivity K. I checked the code again for the solidWallTemperatureCoupled/solidWallHeatFluxTemperatureCoupled BC as well as for the dimensions of K (I use OF 1.5.x, updated somewhere in February 09), and figured out the results below.

Fluid (see chtMultiRegionFoam/fluid/createFluidFields.H): Kf=thermof[i].Cp()*thermof[i].alpha().
The basicThermo.H gives us the dimenions for Cp and alpha, resulting in the following dimenion of Kf: [K]=J/(kgK)*kg/(ms)=W/(mK).

Solid (see 0.001/solidRegion/K): dimensions [ 1 1 -3 -1 0 0 0 ]; thus [Ks]=kgm/(s^3K)=W/(mK).

It seems that the problem you mentioned was fixed for the version I m currrently using.

Back to the my set-up, I observed the unrealistic low temperatures in the region of the adiabatic floor where no cht is activated. Hence I assume the problem does not lay in the way how cht is treated but in the open boundary (maxX) where hot air exits the domain at the top and ambient air enters at the bottom. You mentioned that you are working on a similar problem. Would you mind to tell me the BCs you use for velocity, pressure diffrence and pressure for this type of open boundary?

Thanks again for your comments!

mabinty May 12, 2009 15:13

Hi again!

I forgot to mention that the unrealistic temperature drop is even greater for my fines grid (uniform grid spacing of 0.005 m).

Best regards,

psosnows May 12, 2009 16:34

Hi Aram!

Regarding the dimensions misfit-
my version (OF-1.5x from around Oct 08 ) had in fluid/createFluidFields.H :
thermof[ i ].rho()*thermof[ i ].Cp()*thermof[ i ].alpha()
so as I thought- a rho too much ;)

and to the similarity:
I work with closed volume- natural convection inside light bulbs (filament, gas, glass, surrounding, etc.)
so unfortunately I do not have experience in open boundary conditions.

In the end I allowed myself to take closer look at the boundaries you described in the first post.
I understand that the case looks in rought description: no walls, roof and a hole in the floor which blows hot air.
here are some questions/suggestions about the boundaries:

the _U field:
since the inlet value is (0 0.8 0), why its direcion points in +X?
inletDirection uniform (1 0 0);

the _T field:
you wrote that air_to_ceiling has type:
air_to_ceiling: solidWallMixedTemperatureCoupled;
is this "Mixed" a misspeling? I understand from your last replay that there were changes made in chtMultiReigion but I was quite sure there were two boundaries for coupling: solidWallHeatFluxTemperatureCoupled and solidWallTemperatureCoupled

and the second thing (similar to velocity)
the inlet has:
inlet: fixedValue uniform 473;
but the inletValue is:
inletValue uniform 293;

regarding _pd:
I never used fixedFluxBuoyantPressure but from what I see in the code:
the value you introduce does nothing: in the beggining the value at the boundary is set the same as the internal field next to the patch and in the next calls for values only the gradient is fixed and dependent on rho.
Is this what you want?

and last thing about whole minX patch
it is set as symmetryPlane. From what I understood you want the air to go in and out freely on the horizontal walls; maybe it would be better to set the condition the same as for maxX?

Those are just some free thoughts and I hope you will find something useful and that it will help a bit.

best regards,

mabinty May 13, 2009 05:07

Hi Pawel!

Again, thanks a lot for your comments!!

Yes, it seems that the last version of the cht-BC had a mistake concerning dimensions but is fixed now. And additionally with the up-date a new type of cht-BC was introduced, the mentioned solidWallMixedTemperatureCoupled BC I m using (solidWallTemperatureCoupled/solidWallHeatFluxTemperatureCoupled still available).

Concerning your remarks on the BCs I use:

the _U field:
the keyword inletDirection uniform (1 0 0); belongs to the pressureDirectedInletOutletVelocity BC prescribed at maxX (the open boundary), and sets the direction of the flux from which U is determined. A brief description on this type of BC can be found on page U-129 of the user manual.

the _T field:
the value for the inlet in the floor is 473. the keyword inletVlue uniform 293;
belongs to the inletOutlet BC set at the open boundary maxX (if flux > 0, set T to zeroGrad; for flux < 0, set T to inletValue).

regarding _pd:
exact, fixedFluxBuoyantPressure fixes the gradient dependent on rho; I already tried it on the open boundary maxX, but it gave me strange results, as hot air was not exiting the domain. Henc, I set pd to zeroGradient which produced a plausible U-field.

about whole minX patch:
symmetryPlane is set there because the investigated set-up is symmetrical. Therefore minX is not an open boundary.

What I m currently trying is to study the isothermal case (temperature of air entering from the inlet in the floor = to ambient temperature). The flow field and temperatures (292.99 to 293.0) look plausible, but cumulative time step continuity errors are always around 0.0001179611. Arn t they too big??

Best regrads,

mabinty May 18, 2009 10:44

Dear all!!

After playing around with several parameters, from discretization schemes, to solver settings and BCs I was finally able to obtain a stable calculation. Like most of the times, the problem was to find the appropriate BCs. First, I changed the pressure BC at the open boundary (maxX) to fixedMeanValue with a meanValue=100000 ( This change did not stopped the solution to blow up but allowed to get lower continuity errors. Finally the BCs for k and epsilon turned out to be the critical point. On the same patch I prescribed inletOutlet for k and epsilon, where the inletValues were set to 1/10 of the values I took for the inlet. This settings then allowed to achiev a steady state after at least 100 s.

All the best,

mabinty May 26, 2009 08:20


Due to the fact that I managed to set up a stable simulation I was able to analyse the described configuration in detail. What is strange for me are the results for the temperature field. Even though the temperature of incoming air is 473 K (200 C) and the distance between inlet and ceiling is 1 m, the max predicted ceiling temperatures are only about 10 K above ambient (293 K). I really expected higher values for temperatures along the ceiling. I also refined the mesh constantly to 0.005 m and got same results. Does anybody have a comment on that??

Thanks in advance,

Khelian973 July 29, 2009 03:58

Hi Aram,

First thanks for you help, i could solve my problem and it runs now :)

I looked the previous threads. You are using solidWallMixedTemperatureCoupled bc but i dont have this one in my version (OF 1.5). Only the two bc like Pawel : solidWallHeatFluxTemperatureCoupled and solidWallTemperatureCoupled.

I have two questions.
- What is the difference between Mixed bc and the two others?
- Is there a package to add this bc or do i need to install the last version (which i guess it's 1.5-dev) ?


mabinty July 29, 2009 04:26

Hi Kyian!

The solidWallMixedTemperatureCoupled BC came with an update for OF 1.5.x, and I would assume that it is also included in the new version 1.6. So I think the best is to update your insatllation and check the codes for the BCs which you can find in solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields. In my cases I tried both types of BC for coupled heat transfer and always got same results.

Hope that helps! All the best,

Khelian973 July 29, 2009 05:34

Hi Aram

Thank you for the information. I will probably do the update this weekend.

I explain you my problem. In fact Im not sure how to define a BC. I have a fluid and a solid above that fluid. The top of the solid is warmed (a fixed value of temperature) and i want to study how long it takes to warm the fluid and how will be the temperature contours.

So the patch between the solid and the fluid is a simple solid wall coupled BC. What do you suggest me to use for this BC?

Thanks in advance

mabinty July 30, 2009 04:46

Dear Kyian!

As you will see in the new release 1.6, OF does only provide the mixed type BC, i.e. solidWallMixedTemperatureCoupled, and solidWallHeatFluxTemperature (for prescription of heat flux) for the solid-fluid coupling with chtMultiRegionFoam.

All the best,

All times are GMT -4. The time now is 14:54.