CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM (
-   -   Radiation Model (

julienh May 18, 2006 07:25

Hi, I am trying to program a
I am trying to program a radiation model using the Gibb's method (spherical harmonics)

I used the "reactingFoam" solver, where I had a file IEqn.H with the following radiation equation :


fvScalarMatrix IEqn


- (I/(4*boltzmannCoeff)
+ pow(T,4))*(12*pow(knu,2 *boltzmannCoeff)







where I is the radiativ intensity I wanted to get from this equation

After the compilation, when I execute reactingFoam, I get this message :

--> FOAM FATAL ERROR : gradientInternalCoeffs cannot be called for a calculatedFvPatchField.

You are probably trying to solve for a field with a calculated or default boundary conditions.

From function calculatedFvPatchField<type>::gradientInternalCoef fs() const in file fields/fvPatchFields/basicFvPatchFields/calculated/calculatedFvPatchField.H at line 174.

In order to solve the I equation, I need of course some boundary conditions. I though the file I in the "0" directory of the case with the initial and boundary conditions was enough but it is apparently not the case.

I try to see in the code an equation file like UEqn.H or hEqn.H but I did not see how the boundary equations are coded

Could somebody give me some advices please ?

Thank you in advance


hjasak May 18, 2006 08:26

Looks like one of the b.c.-s t
Looks like one of the b.c.-s that you specify in 0/I is not correct or that the field I has been created by calculation rather then read in. Have a look at the place where I is created and pay attention to the MUST_READ and only a mesh argument after the IOobject. For example:

volScalarField I

Yours should look exactly the same.


julienh May 18, 2006 10:13

Hi, Thank you for this quick
Thank you for this quick answer. I have exactly the same thing for I in the creatField file.
I tried to change the boundary conditions but nothing work, I still have the same error about boundary conditions

Best regards


lucchini May 18, 2006 16:55

Hi Julien, it can be a troubl
Hi Julien,
it can be a trouble of field definition (or in the createFields or the 0 directory), or a trouble when you try to solve the equation.
Could you paste here how is defined I in the 0 directory? And how did you define it in the createFields?
What is m?
thanks. bye

julienh May 19, 2006 07:24

Hi, I found the mistake, I h
I found the mistake, I had a problem in the calculation of the emissivity, my code works quite well now and calculates the radiativ energy field.

I will now try to couple this equation with the basic thermal equation because the convective part of the energy flow emitted by a flame has a direct influence on the radiativ part.

But I have difficulties to find where the thermal equation is coded. Does somebody know ?

Best regards


mattijs May 23, 2006 04:28

Usually hEqn.H
Usually hEqn.H

mprinkey June 22, 2006 06:51

As you have it coded, I believ
As you have it coded, I believe the (I/(4*boltzmannCoeff)) term is being treated as an explicit source term and ending up on the RHS. I think what you want to do is use fvm::Sp() to linearize that term and place it on the center coefficient instead. That should improve solution stability.

julienh June 23, 2006 03:57

Hello, Thank you for this

Thank you for this quick answer.

I tried to do it with fvm::Sp() but it did not work. But I found one of my mistakes. The maxDeltaT was to large to satisfy the convergence criteria. Now it works but my model uses a fixed value of the optical thickness L and now I want my model to calculate L for each cell with the formula L = (2/3)*(volume of the cell)^(1/3)

I coded : volScalarField L = (2/3)*pow (mesh.V(), (1/3)) but it did not work.

Could somebody help me please ?

Thank you in advance

Best regards


eugene June 23, 2006 04:09

Try volScalarField L = (2.

volScalarField L = (2.0/3.0)*Foam::pow(mesh.V(), (1.0/3.0));

If you do this
a = (1/3)
a = 0

Because you are doing integer division.

julienh June 23, 2006 08:46

Thank you very much ! It works
Thank you very much ! It works like this.

So my very simplified radiation model enables me to calculate the radiation flux I inside my combustion chamber. But to do it, I put boundary conditions for I in the 0 file like fixedValue or zeroGradient. However, the goal of my work is to calculate the field I on the wall of a furnace. The idea is to modify the wallHeatFlux file to add the radiatif flux to the convectif flux. But to do that, instead a BC like zeroGradient for I, I need to put a BC (I0) which is the solution of a new equation :

(emissivityWall/4)*I0 - boltzmannCoeff*emissivityWall*Twall^4=(1/(3*a))*grad(I(0))

Am I obliged to create a new BC class to do that or is there a other way to do it ?

Can somebody give me some advice please ?

Thank you in advance

Best regards


julienh June 27, 2006 03:36

Hello, So, I can calculate th
So, I can calculate the radiation flux within the furnace if I put existing boundary conditions like zeroGradient. But I need to put a calculated value for the boundary condition Iwall of field I. The code calculates Iwall but now I want to put Iwall as a boundary condition for I. So, I think I must code a new BC class. This class must just enable me to fill the /0 field for I like this :

internalField uniform 25000;

type boundaryConditionsforI;
value Iwall;

It is perhaps easy to do that butI am not a very good programmer.

Could somebody give me some advise to program my own bounday conditions class please ?

Thank you in advance

Julien Heintz

julienh June 28, 2006 10:16

Hello, Finally, I think I do
Finally, I think I do not need to do my own boundary condition class for the field I so I tried to do it differently : here the code I had to my radiation solver in order to calculate the boundary field Iwall :

Info<< "\nCreate new boundary condition\n" << endl;

label inletPatchID = mesh.boundaryMesh().findPatchID("walls");

fvPatchScalarField& Iwall = I.boundaryField()[inletPatchID];

forAll(Iwall, faceI)

surfaceScalarField Iwall((0.6/epsilonWall)*fvc::snGrad(I) + 4*5.67e-8*pow(Tparoi,4));

I.boundaryField()[faceI] = Iwall;


Info<< "\n ExecutionTime = " << runTime.elapsedCpuTime() << " s\n" << endl;


but it doesn't work. The field I is calculated but the boundary values of I remain constant, equal to the value I put in the 0/I file.

What am I doing wrong ?

Thank you in advance


All times are GMT -4. The time now is 19:45.