CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM (
-   -   forced convection of air inside a heated rectangular channel (

Ank June 3, 2013 08:11

forced convection of air inside a heated rectangular channel

I am simulating the forced convection of air inside a heated rectangular channel. The air velocity is 1 m/s at the inlet. The temperature of the all four walls is 400 K and air is at 300 K. Initially I have taken the case of pitzdaily in simpleFoam solver and included the energy equation to incorporate the heat transfer. Then I used buoyantBoussinesqPimpleFoam solver but I am not able to get the convergence in pressure. I have given the velocity in x direction, but it does not show any velocity in that direction, instead it is showing the effect of buoyancy only, and show the circulating cells. I have tried the totalPressure, fixedPressure boundary conditions at the inlet and zeroGradient at the outlet, but no improvement.
I am using kEpsilon model with first order upwind. For testing purpose the grid is a bit coarser. I am able to run the case without energy equation successfully with the simpleFoam solver, but problem comes whenever I include the energy equation.

Please help me, it seems an easy problem but I am really stuck...


Ank June 4, 2013 03:23

1 Attachment(s)

I just checked my results. The velocity vectors are showing that the velocity reaches to zero value near the outlet. It is non physical as velocity should be 1 m/s everywhere.
I am attaching the snapshot. I dont know what is happening.

Please have a look.

Ahmed Khattab June 5, 2013 05:31

Hi Ankur,

could you please upload case files.


Ank June 5, 2013 06:38

3 Attachment(s)

I found the solution to my earlier problem, Now the problem is different and the velocity which is 1 m/s at the inlet is going above to 1.14 m/s at the outlet. I am posting the case files and snapshot.

I have made a silly mistake in the calculation of epsilon. I will get back to you once fix it. To get the pure forced convection my gravity and expansion coefficients are both zero.

Ank June 5, 2013 06:50

1 Attachment(s)
The problem remains same with correct calculations also. I am attaching more files now. The velocity is going beyond its correct value.

deji June 5, 2013 07:19

Perhaps you should try buoyantPressure BC for the pressure at the inlet and outlet. This is the "perturbation" pressure that you are referring to, correct?

Ank June 5, 2013 07:42

Hey Deji

I do not understand your reply. I have taken zero gravity so do you still think It will work?
secondly I tried it but it is diverging.
Now pressure is also converging but the velocity is going beyond its value as I posted above.


Tushar@cfd June 5, 2013 08:10


Originally Posted by Ank (Post 431618)

For testing purpose the grid is a bit coarser. I am able to run the case without energy equation successfully with the simpleFoam solver, but problem comes whenever I include the energy equation.


What do you mean to say by the mentioned Quote? Are you running your own solver "simpleFoam" with included energy equation in to it?

Ank June 5, 2013 08:13

Hey Tushar

Good to see your reply again..

No I am using my own solver. Sorry for the language. It meant when I shifted from simpleFoam to buoyantBoussinesqSimpleFoam. It seems to be a bit easy problem but I am getting strange velocity increment at the outlet.

And the case I have taken is with very coarse grids, dimensions 0.1 X 0.05 X 0.5 with 17000 hexa cells.

deji June 5, 2013 08:16

My mistake. Since g is set to zero, that BC would be invalid.

deji June 5, 2013 08:19

It seems like your flow is accelerating inside the channel due to a negative pressure gradient, -dP/dx. What kind of pressure boundaries do you have at the inlet and outlet? And is this pressure the "perturbation pressure"?

Ank June 5, 2013 08:22

Hey Deji

Thanks for the reply..

I am using zeroGradient at the inlet and fixedvalue (uniform 0) at the outlet. I am sorry but I dont know the meaning of perturbation pressure.

Normally the simulation runs fine with these BCs when heating is not there.

deji June 5, 2013 08:26

Do me a favor, look in solver and paste the pressure equation and how the pressure is solved. You are probably solving a pressure equation, I would think.

Ank June 5, 2013 08:33


I think I got the perturbation pressure, is it the pressure correction term in the SIMPLE algorithm?

Here is the solver

00002 volScalarField rUA("rUA", 1.0/UEqn().A());
00003 surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
00005 U = rUA*UEqn().H();
00006 UEqn.clear();
00008 phi = fvc::interpolate(U) & mesh.Sf();
00009 adjustPhi(phi, U, p_rgh);
00011 surfaceScalarField buoyancyPhi = rUAf*ghf*fvc::snGrad(rhok)*mesh.magSf();
00012 phi -= buoyancyPhi;
00014 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
00015 {
00016 fvScalarMatrix p_rghEqn
00017 (
00018 fvm::laplacian(rUAf, p_rgh) == fvc::div(phi)
00019 );
00021 p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
00023 // retain the residual from the first iteration
00024 if (nonOrth == 0)
00025 {
00026 eqnResidual = p_rghEqn.solve().initialResidual();
00027 maxResidual = max(eqnResidual, maxResidual);
00028 }
00029 else
00030 {
00031 p_rghEqn.solve();
00032 }
00034 if (nonOrth == nNonOrthCorr)
00035 {
00036 // Calculate the conservative fluxes
00037 phi -= p_rghEqn.flux();
00039 // Explicitly relax pressure for momentum corrector
00040 p_rgh.relax();
00042 // Correct the momentum source with the pressure gradient flux
00043 // calculated from the relaxed pressure
00044 U -= rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rUAf);
00045 U.correctBoundaryConditions();
00046 }
00047 }
00049 #include <finiteVolume/continuityErrs.H>
00051 p = p_rgh + rhok*gh;
00053 if (p_rgh.needReference())
00054 {
00055 p += dimensionedScalar
00056 (
00057 "p",
00058 p.dimensions(),
00059 pRefValue - getRefCellValue(p, pRefCell)
00060 );
00061 p_rgh = p - rhok*gh;
00062 }
00063 }

Tushar@cfd June 6, 2013 00:02

Hello Ank,

Since, Boussinesq approximation is used for the pressure, as you are using buoyantBoussinesqSimpleFoam. So, you need to provide the BC for p_rgh something like this....

type buoyantPressure;
rho rhok;
value uniform 0;

see the tutorial on heatTransfer/buoyantBoussinesqSimpleFoam/ you will get an idea.

Ank June 6, 2013 01:42

Hey Tushar

I tried that also, but it didnt work. Another thing is that its totally forced convection and my gravity is also zero so buoyant does not mean anything here as I have understood. I do not know what else can go wrong in this simple problem.

I will try more things now, lets see how it goes. You can also think what else can be wrong.


Tushar@cfd June 6, 2013 01:54


I mean to say by putting the mentioned BC for "p_rgh" and setting g=0, your case will be of forced convection only (i.e., independent of buoyancy).
I have a doubt with pressure "p" BC. Try setting the BC for Outlet & inlet.
If you are not getting solution.
Also check:
Are you getting convergence?

Ank June 6, 2013 02:19

I am posting my p_rgh and U files, kindly suggest me what I should change,

dimensions [0 2 -2 0 0 0 0];

internalField uniform 0;

type zeroGradient;

type zeroGradient;


type fixedValue;

value uniform 0;


type zeroGradient;


dimensions [0 1 -1 0 0 0 0];

internalField uniform (0 0 0);


type fixedValue;

value uniform (1 0 0);

type fixedValue;
value uniform (0 0 0);

type zeroGradient;


type fixedValue;
value uniform (0 0 0);

I have tried changinp_rgh to buoyantPressure it was not converging.

Tushar@cfd June 6, 2013 02:25

Check my earlier post for "p" BC...

Ank June 6, 2013 02:45

I changed the p_rgh to buoyantPressure and p file to zeroGradient at the inlet and fixedValue at the outlet.

It is showing this error

Continuity error cannot be removed by adjusting the outflow.
Please check the velocity boundary conditions and/or run potentialFoam to initialise the outflow.
Total flux : 0.000771162
Specified mass inflow : 0.005
Specified mass outflow : 0
Adjustable mass outflow : 2.34371e-110

From function adjustPhi(surfaceScalarField& phi, const volVectorField& U,const volScalarField& p
in file cfdTools/general/adjustPhi/adjustPhi.C at line 118.

FOAM exiting

All times are GMT -4. The time now is 18:49.