[maddalena]

Dear FOAMers,
in order to simulate a conduction + convection problem with openFoam, I added a heat source to the standard chtMultiRegionFoam release. This was straightforward: in solveSolid.H I added the heat source H [W m-3] contribution to the general conduction equation:

{
    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
    {
        tmp<fvScalarMatrix> TEqn
        (
            fvm::ddt(rho*cp, T)
          - fvm::laplacian(K, T)
          - H
        );
        TEqn().relax();
        TEqn().solve();
    }

    Info<< "Min/max T:" << min(T) << ' ' << max(T) << endl;
}

This formulation runs fine if I use a single region, but fails when using two or more regions. In that case, the temperature inside every region raises to unacceptable values, and the theoretical solution known for the easiest problems is not reached (see also here)
After some more investigations, I noticed that the problem may be due to the coupling between interfaces. Indeed, the solidWallMixedTemperatureCoupled formulation calculates Twall as:

        scalarField Twall
        (
            (myKDelta()*intFld() + nbrKDelta*nbrIntFld)
          / (myKDelta() + nbrKDelta)
        );

I am guessing that this equation implements formula 3.15, pag 32 of this book: given three grid points W, P, and E and being dxw and dxe the distance between WP and PE respectively, than a_P*T_P = a_E*T_E+a_W*T_W + H*deltaX , where a_P = a_E+a_W, a_E = k_E/dxe, a_W = k_W/dxw. (T temperature, k conducibility). Therefore, in the case of heat source I should have something like

        scalarField Twall
        (
            (myKDelta()*intFld() + nbrKDelta*nbrIntFld + myHDelta())
          / (myKDelta() + nbrKDelta)
        );

The last term represents the heat source*DeltaX contribution. However, this is valid for a Temperature inside the grid and not at the boundary. Indeed, the heat source contribution in this case is applied to the neighbor cell as well.
I am thinking that I need a more general review of solidWallMixedTemperatureCoupled.C, but I am asking for suggestions and hints before going further:

1. Has anyone have any kind of experience in adding heat source to chtMultiRegionFoam? I do not think I am the first to try it, and I know for sure that somewhere else had similar problem with OF (see here) ... How did you solved that?
2. In the case I am the pioneer in this sense , is my approach correct? Is there anyone that can comment on that?
3. Beside this book, is there any reference that gives a numerical formulation for chtt?

Thanks for any help I may get.

Regards,

mad

[maddalena]

Following a presentation of the coupling made by HRV, I ended up with the following.
In the case of coupling between two regions A and B, both of them with a heat source H, the coupling condition should be something like:

• T_A = T_B
• K_A*grad(T_A) +H_A*x = K_B*grad(T_B) +H_B*x

This will help to keep the continuity between the two regions & take into consideration the changing slopes of the two "x vs T" curves.
I was thinking to implement this as:

tmp<scalarField> myHDelta = 2*H()/patch().deltaCoeffs();    
this->refValue() = nbrIntFld;
this->refGrad() = - myHDelta() / K();
this->valueFraction() = nbrKDelta / (nbrKDelta + myKDelta());

but without success.
Is there anyone that can comment on that? Is this the right approach or am I on the wrong path?
thank you

mad.

[maddalena]

Ok, the ideas presented above were unsuccessful, probably due to my poor C++ and / or something on the math formulation of the problem that I am missing.
Someone else agrees with me that the solidWallMixedTemperatureCoupled is not general enough to allow a heat source inside one or more regions. Unfortunately, it seems like that everybody facing this problem has given up and turned to a fixedValue approach. That is what I probably will do as well, at least for the moment.
On the other hand, I am really interested on a new coupling formulation, and I'd like to spend some of my time on it. However, I need some general ideas on how to approach the problem, since I am still not completely into it. The first step, for example, can be some comments on the ideas described above...

regards,

mad

[andrea]

Hi,
first I would suggest you to check if the problem is really the boundary condition, did you try to compute the continuity error?
For "engineering" porpoise the coupling mixed condition gives good result, but maybe the source term is giving troubles, did you try to solve implicitly, using SuSp()?

If the power source is generated also along the boundary, which maybe is unphysical, you should add it in the balancing eq. otherwise the original should do the trick.

Andrea

[maddalena]

Hi Andrea and thank you for replying!

Originally Posted by andrea

first I would suggest you to check if the problem is really the boundary condition, did you try to compute the continuity error?

What I have done up to now is to create a sort of matrix (see attached picture) with one or two regions, and applying energy using a temperature fixed value or a heat source. As you can see, the only simulation that fails is the 2 region + heat source, thus I concluded that:

Originally Posted by andrea

For "engineering" porpoise the coupling mixed condition gives good result, but maybe the source term is giving troubles

1. How can I check continuity? If I use two solids as I did, it does not appear on the standard output.
2. How can I switch to implicitly? where using SuSp()?

Originally Posted by andrea

If the power source is generated also along the boundary, which maybe is unphysical, you should add it in the balancing eq. otherwise the original should do the trick.

The power is generated on every cell in the domain, thus on the boundary cell as well. The solidWallMixedTemperatureCoupled does not take into account the T gradient between the last cell center and the boundary, as you have if applying a heat source, and the Twall does not include the heat source as well. This is what I was trying to do on the previous posts: adding the heat source contribution on the boundary balancing equation. Is this what you mean?

cheers,

mad

[andrea]

Hi,
about source term read this:
http://www.cfd-online.com/Forums/ope...r-details.html

If you don't want source at boundary create a vol Field with 0 fixed value boundary condition and add it to the equation to be solved, maybe this is worth a try.

To compute errors you should write your own little code... not too difficult.

This is what I was trying to do on the previous posts: adding the heat source contribution on the boundary balancing equation. Is this what you mean?

yes I meant that, but I think it's not physical, my guess.
Hope it helps
Andrea

[maddalena]

Hello,

Originally Posted by andrea

If you don't want source at boundary create a vol Field with 0 fixed value boundary condition and add it to the equation to be solved, maybe this is worth a try.

this was easy: I defined H in my domain in the same way I declare K or cp, thus I simply had to vary the boundary field as follow:

internalField   uniform 500;

boundaryField
{
    SyCube
    {
        type            symmetryPlane;
        value           uniform 500;
    }
    FB
    {
        type            empty;
        value           uniform 0;
    }
    cube1_to_cube2
    {
        type            fixedValue;
        value           uniform 0;
    }
    CyCube1
    {
        type            cyclic;
    }
}

well, it did not work.
I will try the SuSp() as soon as I understand where I should change the code..

cheers,

mad

[andrea]

Hi,
ok then, I'm curious, how much this heat source is large?
Have you tried to lower H's value and see if there's any hope to get
convergence?

I will try the SuSp() as soon as I understand where I should change the code..

if I remember correctly it should be as simple as

SuSp(H)

or

H.SuSp()

into the solve function or the fvMatrix inizialization.

[maddalena]

Hi Andrea,

Originally Posted by andrea

how much this heat source is large? Have you tried to lower H's value and see if there's any hope to get convergence?

The heat source is 500 W/m^3, but I tried to lower till 0.5 W/m^3. The most weird thing is that the time vs T tendency is linear, and I have this even if I let the simulation run longer or decrease H or change the material properties.

Some thoughts about that... Let us for a moment drop the heat source contribution and consider the (working) case if I impose a fixedValue temperature. From what I can understand, the solidWallMixedTemperatureValue copies the temperature value from the neighbour region, thus from the "cold" region, on my cell in the "warm" region. On the "cold" region, a zeroGradient BC is applied. In different (and easier) words the temperature of the "cold" region is diffused on the "warm" region and viceversa. This works fine if I do not have any heat source.
When I add H, the contribution of the "cold" region on the "warm" region is not enough, since the temperature gradient obtained on the boundary cell by coping "cold" cell value onto the "warm" cell value is lower than the heat source contribution. In other words, the temperature diffusion in the "cold" region works fine, but not in the "warm" region. Therefore, the "warm" region is somewhat isolated from the outside, and its temperature increases until the simulation ends. This is, at least, the explanation I am giving to the weird results I am getting.

Originally Posted by andrea

if I remember correctly it should be as simple as

SuSp(H)

or

H.SuSp()

into the solve function or the fvMatrix inizialization.

Indeed, I should have written:

        tmp<fvScalarMatrix> TEqn
        (
            fvm::ddt(rho*cp, T)
          - fvm::laplacian(K, T)
          - fvm::Su(H)
        );

on my solveSolid.H, to let the solver consider only the constant part of the source term linearization. In that case, I get a "Unknow matrix combination" error before starting the solver, and the LU or something operator complaining. (sorry, I am not at the office now, I will report the complete error message as soon as I will go back there). As for using a SuSp(), should not I have a sort of linear equation for the heat source term to use it? Something like H = H0 + H1*T? At the moment I only know the H0, that is why I wrote the solveSolid.H equation as in the first post. I am going to google a bit to see if I can find a sort of relation that is applicable in my case, and than update solveSolid.H as you suggested.

In the meanwhile, another question come into my mind. chtMultiRegionFoam is a compressible solver, while I want to consider my simulation as incompressible. May I get problems if I use it in my case? I would say no, but maybe there is something that I am missing.

Thanks for your time,
regards,

mad

[andrea]

Hi again,
maybe you know about that but there was a sort of "bug" in the coupling patch in previous releases of OF, in 1.5 was fixed.

In the meanwhile, another question come into my mind. chtMultiRegionFoam is a compressible solver, while I want to consider my simulation as incompressible. May I get problems if I use it in my case? I would say no, but maybe there is something that I am missing.

No, the coupling works fine since you are trying to force T and Heat Flux matching between two regions what ever the eq. to be solved on both sides . Actually, we used this thing for an incompressible case.

I still would like to see the coupling error between the two patches, can you try something as

 Info << "error: "<<Ts[i].boundaryField()[patchNumber]-Tf[i].boundaryField()[patchNumber]<< endl;

this should gives you how bad is the coupling working.

you're welcome,
Andrea

[maddalena]

Hi Andrea,

Originally Posted by andrea

maybe you know about that but there was a sort of "bug" in the coupling patch in previous releases of OF, in 1.5 was fixed.

I am using 1.6.x coupling, and modified version of setInitialMultiRegionDeltaT.H to take into account the time step selection on the base of maxDi or maxCo.

Originally Posted by andrea

No, the coupling works fine since you are trying to force T and Heat Flux matching between two regions what ever the eq. to be solved on both sides . Actually, we used this thing for an incompressible case.

That sounds nice. This is something I expected.

Originally Posted by andrea

I still would like to see the coupling error between the two patches, can you try something as

 Info << "error: "<<Ts[i].boundaryField()[patchNumber]-Tf[i].boundaryField()[patchNumber]<< endl;

this should gives you how bad is the coupling working.

In place of changing the code, I looked at the lastTimeStep/../T file for both regions. In my simulation, H is 500 W/m^3 and placed inside cube1. The coupling is based on the 1.6.x unmodified version of solidWallMixedTemperatureCoupled.H. So... At time 0 the temperature is 300 K everywhere. This is 100000/cube1/T:

cube1_to_cube2
{
type solidWallMixedTemperatureCoupled;
refValue nonuniform List<scalar> 10(520.5101 520.5101 520.5101 520.5101 520.5101 520.5101 520.5101 520.5101 520.5101 520.5101);
refGradient uniform 0;
valueFraction nonuniform List<scalar> 10(0.002617133 0.002617133 0.002617133 0.002617133 0.002617133 0.002617133 0.002617133 0.002617133 0.002617133 0.002617133);
value nonuniform List<scalar> 10(546.0636 546.0636 546.0636 546.0636 546.0636 546.0636 546.0636 546.0636 546.0636 546.0636);
neighbourFieldName T;
K K;
}

and this is 100000/cube2/T

cube2_to_cube1
{
type solidWallMixedTemperatureCoupled;
refValue nonuniform List<scalar> 10(546.1307 546.1307 546.1307 546.1307 546.1307 546.1307 546.1307 546.1307 546.1307 546.1307);
refGradient uniform 0;
valueFraction uniform 0.9973829;
value nonuniform List<scalar> 10(546.0637 546.0638 546.0638 546.0638 546.0638 546.0638 546.0638 546.0638 546.0638 546.0638);
neighbourFieldName T;
K K;
}

As you can see, cube1 has a valueFraction close to 0: a fixedValue BC is applied, and the 520.5101 value is the innerField value of the cube2 first cell. cube2 has a valueFraction close to 1: a fixedGradient value is applied, and the refValue is the boundary value. Comparing the temperature value at the interface (value entry), I can see the temperatures are (almost) the same, thus the coupling is working fine.

Still working on the heat source definition. This is the error I am getting when using fvm::Su(H, T):

--> FOAM FATAL ERROR:
Unknown matrix type combination

From function lduMatrix:: operator-=(const lduMatrix& A)
in file matrices/lduMatrix/lduMatrix/lduMatrixOperations.C at line 274.

FOAM aborting

Stay tuned for update.
regards,

mad

[andrea]

Hi,
to check the coupling properly you should also check the heat flux, and for this you have to add two code lines , don't be too scared... If the heat is not flowing to the next regions of course the temperature will rise whatever the source implementation, I think you had a feeling about this too.
I wait for some new updates, good luck!
Andrea

[maddalena]

Originally Posted by andrea

to check the coupling properly you should also check the heat flux, and for this you have to add two code lines , don't be too scared... I

I am not scared... I only do not want to loose time in implementing things that I already knows... For example, the heat flux can be obtained activating the debug flag. This is the last simulated time:

Solving for solid region cube1
cube1:cube1_to_cube2:T <- cube2:cube2_to_cube1:T : heat[W]:-1.338517 walltemperature min:546.0153 max:546.0153 avg:546.0153
[...]
Solving for solid region cube2
cube2:cube2_to_cube1:T <- cube1:cube1_to_cube2:T : heat[W]:1.34105 walltemperature min:546.0154 max:546.0154 avg:546.0154

The two flux are (almost) the same, so the flow seems to be modeled accurately.
In any case, where should I add the lines you suggested? In solidWallTemperatureCoupled.H or in chtMultiRegionFoam.H? Is Ts the solid temperature at the boundary and Tf the fluid temperature? If so, what does change if I add the lines or check the T value at the coupling region?

mad

[andrea]

Ok then,
from what I see it seems fine except for a little heat flux error, maybe not enough for a totally wrong simulation, is this error so little in the first time steps?
The problem is that you get unrealistic too high values inside the source domain, right?
I noticed that the coupling gives slower temperature rising than it should, maybe you "accumulate" too much heat in the beginning (I'm just thinking aloud now)

Andrea

[maddalena]

Hi andrea, and thank you to keep replying.

Originally Posted by andrea

from what I see it seems fine except for a little heat flux error, maybe not enough for a totally wrong simulation, is this error so little in the first time steps?

yes, that small error is always there. After 5000 s simulating time, the output is as follows:

Region: cube1 Diffusion Number mean: 0.009960159 max: 0.009960159
Region: cube2 Diffusion Number mean: 0.004417433 max: 0.004417433
deltaT = 19.92032

Time = 5000
Solving for solid region cube1
cube1:cube1_to_cube2:T <- cube2:cube2_to_cube1:T : heat[W]:-0.1132085 walltemperature min:312.3995 max:312.3995 avg:312.3995
DICPCG: Solving for T, Initial residual = 0.8195623, Final residual = 1.043602e-07, No Iterations 4
cube1:cube1_to_cube2:T <- cube2:cube2_to_cube1:T : heat[W]:-0.1158047 walltemperature min:312.4491 max:312.4491 avg:312.4491
DICPCG: Solving for T, Initial residual = 2.881033e-07, Final residual = 2.881033e-07, No Iterations 0
Min/max T:min(T) [0 0 0 1 0 0 0] 312.4491 max(T) [0 0 0 1 0 0 0] 312.4942

Solving for solid region cube2
cube2:cube2_to_cube1:T <- cube1:cube1_to_cube2:T : heat[W]:0.1158047 walltemperature min:312.3996 max:312.3996 avg:312.3996
DICPCG: Solving for T, Initial residual = 0.00406011, Final residual = 2.108958e-07, No Iterations 2
cube2:cube2_to_cube1:T <- cube1:cube1_to_cube2:T : heat[W]:0.1134764 walltemperature min:312.4492 max:312.4492 avg:312.4492
DICPCG: Solving for T, Initial residual = 2.102993e-07, Final residual = 2.102993e-07, No Iterations 0
Min/max T:min(T) [0 0 0 1 0 0 0] 300 max(T) [0 0 0 1 0 0 0] 312.4492
ExecutionTime = 0.2 s ClockTime = 0 s

At 5000 s the steady state is not reached: for cube1, the time constant should be around 40000s. Note that the delta T is chosen by OF automatically, but it matches the max admissible delta T, calculated as suggested in other threads, which is 100 s.

Originally Posted by andrea

The problem is that you get unrealistic too high values inside the source domain, right?

This is one of the problems... summarizing:

1. the temperature is too high inside the source region. As a consequence, the temperature is too high inside the "cold" region as well.
2. The time vs temperature has not the theoretical tendency: in my simulation, I have a linear time vs temperature!

Originally Posted by andrea

I noticed that the coupling gives slower temperature rising than it should

How did you noticed that? The final 100000 seconds simulation should be already on a steady state condition....

Originally Posted by andrea

maybe you "accumulate" too much heat in the beginning (I'm just thinking aloud now)

What do you mean? Should I increase the heat source little by little, as one usually does with turbulence on simpleFoam simulation? I have already done that, using three steps:
1. from 0 to 25000 s with H = 5 W/m^3
2. from 25000 to 50000 s with H = 50 W/m^3
3. from 50000 to 100000 s with H = 500 W/m^3
no success. The only things that changes is the time vs temperature line angular coefficient. The theoretical curve rising is not matched.

mad

[andrea]

Hi,
so the coupling seems to work, even if with the "small" enough error in heat flow that unlikely it is the problem.
So what about the others boundary conditions? It seems that you possibly have a zeroGradient somewhere in the colder domain. Let me know.

Bye
Andrea

[maddalena]

Hello Andrea,

find attached a simple figure of my domain. I have zeroGradient BC on K, cp and H (heat source) both at the interface between cube1 and cube2 and on the other vertical patch. As you can see, the grid is rough: should this be a problem? I guess that the solution may be wrong to the poor grid quality, but not SO wrong...

In the meanwhile, I found a working case. If I add the heat source to the multiRegionHeater tutorial, everything works as it should: time vs temperature is not linear but it has the well known curved tendency. Note that the grid is really really rough on the tutorial.
I do not now if I should be happy with that (OF does not betray us) or not happy (there is something wrong with my model). In any case... In order to find what is wrong with it, I cross checked what changes between the two models (tutorial and my own model):

• time step: my time step is much higher than the tutorial time step. Decrease it bring no improvements.
• mesh size: a rough mesh cannot be capable of describing what is happening. I am going to increase the mesh resolution close to the interface.
• cyclic: I am using cyclic, the tutorial is not. May this be cause of the problem?
• Material properties: I "created" a new material for cube1, with some properties I like in order to reduce the computational time, while using air properties for the cube2 (which is calculated as solid). I am going to change them. However, this will not explain the linear time vs temperature...

What do you think?
cheers,

mad

[andrea]

Well,
I agree with your analysis, and then the obviously thing to do is drop the cyclic patch type... then we can re-think on "why you cyclic?"

fingers crossed
Andrea

[maddalena]

So...

• mesh size: Increased mesh resolution at the interface --> No changes on the solution.
• cyclic: run with zeroGradient in place of cyclic --> No changes on the solution.
• material properties: adjusted with tutorial values --> solution changes, but not the time vs temperature.
• time step: tried to reduce time step, well below the prescribed limit of Dt < rho*cp*Dx^2/2/K --> no changes.
• 2D vs 3D: tutorial uses a 3D mesh, while I have a 2D case, maybe patch defined as empty is the problem? -> no, the time vs temperature keeps to be linear.

what should I do now? I am really out of ideas! It seems unbelievable that the only running case is the modified tutorial!

cheers,

mad

[andrea]

ok,
don't loose your temper, the thing is weird but maybe we're still missing something. Can you post the case, along with the solution? I will take a look at it

cheers
Andrea