Hi, I am trying to program a
 

Hi,
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

(

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

==
m

);

IEqn.relax();

IEqn.solve();

thermo->correct();

}

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

Julienh

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
(
IOobject
(
"I",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);

Yours should look exactly the same.

Hrv

Hi, Thank you for this quick
 

Hi,
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

julienh

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
Tommaso

Hi, I found the mistake, I h
 

Hi,
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

julienh

Usually hEqn.H
 

Usually hEqn.H

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.

Hello, Thank you for this
 

Hello,

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

Julienh

Try volScalarField L = (2.
 

Try

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.

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

Hello, So, I can calculate th
 

Hello,
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;

boundayField
{
walls
{
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

Hello, Finally, I think I do
 

Hello,
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;
}

I.write();

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

julienh