[suitup]

Hello,

sorry that I ve to start a new bug thread, but I have really no idea why the boundary condition "solidWallHeatFluxTemperature" is not working properly. The normal temperature boundary seems to be working correct...

I'm working with the chtMultiReagion-Solver, the case example is a little heat sink in an air domain.

If I'm using solidWallHeatFluxTemperature as boundary condition, the heatsource is on the bottom side of the heat sink, the surface is about 0.003713 m², so q is 4040 W/m² (converted from 15 Watt). Strangely enough the temperaturevalue starts always at the starttemperaturevalue and rise upvery slowly. In about 240 seconds max 7Kelvin, a normal value should be around 70-90Kelvin.

Maybe there is a connection problem with the scaling related to the ProE-Bug, but actually I ve rescaled the mesh into millimeter.

The boundary of bc is declared as a normal wall, so maybe there is the problem, and the boundary shouldn't be a wall?

solidWallHeatFluxTemperature after 240sec @ q=4040W/m²



type solidWallHeatFluxTemperature;
K K;
q uniform 4040;
value 293.15;

vixed temperature value= 398.15K after 40sec



type fixedValue;
value uniform 398.15;


Or is there another bc which I could use as heat source, which can declared with heat power in the chtmultireagion?

I've additionally attached the case-file.

Thank you and best regards.

[kpm]

I don't know wether there's still anybody interested in this, but I just happened to run into the same issue with a similar test case.

There actually is another BC that You can use as heatsource, fixedGradient with gradient = q/K

But that won't solve the problem.

The reason for this "inaccuracy" is a lack of automatic time step control for heat conduction in the solid regions.
The solver determines the time step only by considering the Courant number in the fluid regions, so the time step could become too large for solid regions with unsteady temperature fields, and the transient results might be inaccurate to a certain extent (nominated for understatement of the year).
This will occur during the beginning of the heating process, when Your air is not moving.
Once the heatsink becomes hot, the air will start moving and Your automatic time step
control will reduce the timestep and improve the accuracy of Your transient
temperature calculations in the solid regions.
Well, in Your case, the heatsink never became hot within the given time period.

So it's not actually a bug, but merely an inconvenience.

To solve this problem, You might want to limit the maximum time step
manually by adding a "maxDeltaT" entry to Your controlDict.
A possible choice for maxDeltaT would be
"maxDeltaT = 0.1 * rho * cp * (deltaX_of_Your_mesh_in_meters)^2 / K",
using the properties of Your solid region, for an aluminium heatsink maybe something like
maxDeltaT = 0.1 * 2700 * 897 * (deltaX_of_Your_mesh_in_meters)^2 / 235 = 1000 * (deltaX_of_Your_mesh_in_meters)^2 s/m²,
e.g. "maxDeltaT = 0.001" for a mesh with an equidistant resolution of 1mm, so the entry in Your controlDict could be
"maxDeltaT 0.001;"
Just don't forget to adjust Your "maxDeltaT" in case You are refining Your mesh.

The resulting transient temperature field might still be inaccurate, but definitely closer to reality.
The result should improve if You further reduce the time step, e.g. by using a factor of 0.01 instead of the 0.1 used above.
You might also want to use a higher initial temperature in Your solid region to speed up the process if You're just interested in the equilibrium temperature field.

The disadvantage of this approach is that You also have to calculate more time steps in the fluid regions, if Your chosen maxDeltaT is smaller than the time step required for the fluid regions.
Solving the fluid regions makes up a large portion of the computational effort with such a heatsink test case, so You end up wasting CPU time where You don't need it.

I am currently experimenting with some ugly hack in chtMultiRegionFoam.
I use subcycles with smaller time steps in the solid regions.
In my case, the solid regions solve pretty fast, so it doesn't really hurt
performance if I perform 10, 100, or even 1000 subcycles on them.

If You want to try this subcycle experiment, do the following:

--------------------------------------------------------------

1) Add a new entry in the fvSolution for each solid region

SolidRegionSubCycleControl
{
    maxDeltaTForSolidRegion 1E-4;
}

2) Choose Your maxDeltaTForSolidRegion (not necessarily 1E-4)

3) Add a new file in

$WM_PROJECT_DIR/applications/solvers/heatTransfer/chtMultiRegionFoam/solid

named "readSolidRegionSubCycleControls.H"

with content:

    fvMesh& mesh = solidRegions[i];
    const dictionary& SolidRegionSubCycleControl = mesh.solutionDict().subDict("SolidRegionSubCycleControl");

    double maxDeltaTForSolidRegion =
        SolidRegionSubCycleControl.lookupOrDefault<double>("maxDeltaTForSolidRegion", 1.0);

4) Change a part of the code in file

$WM_PROJECT_DIR/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C

from:

   forAll(solidRegions, i)
            {
                Info<< "\nSolving for solid region "
                    << solidRegions[i].name() << endl;
                #include "setRegionSolidFields.H"
                #include "readSolidMultiRegionPIMPLEControls.H"
                #include "solveSolid.H"
            }

to:

            forAll(solidRegions, i)
            {
                        // this part was original in the solver
                        // begin original
                        //Info<< "\nSolving for solid region "
                        //    << solidRegions[i].name() << endl;
                        //#include "setRegionSolidFields.H"
                        //#include "readSolidMultiRegionPIMPLEControls.H"
                        //#include "solveSolid.H"
                        // end original
                // hack: we do some subcycles for solid regions
                // begin hack
                #include "readSolidRegionSubCycleControls.H"
                    Info<< "\nSolving for solid region "
                        << solidRegions[i].name() << endl;
                    Info<< "\nrunTime.deltaT().value() "
                        << runTime.deltaT().value() << endl;   
                    Info<< "\nmaxDeltaTForSolidRegion "
                        << maxDeltaTForSolidRegion << endl;                   
                if ( runTime.deltaT().value() > maxDeltaTForSolidRegion )
                    {
                        int numberOfSubCycles = ceil(runTime.deltaT().value()/maxDeltaTForSolidRegion);
                        Info<< "\nNeed subcycles "
                            << numberOfSubCycles << endl;
                        Foam::TimeState MyNewSubCycleTimeState = runTime.subCycle(numberOfSubCycles);
                        for (int subCycleCounter=1; subCycleCounter<=numberOfSubCycles; subCycleCounter++)
                            {
                                Info<< "\nSubCycle "
                                << subCycleCounter << endl;
                                #include "setRegionSolidFields.H"
                                #include "readSolidMultiRegionPIMPLEControls.H"
                                #include "solveSolid.H"
                                runTime++;
                            }
                        runTime.endSubCycle(MyNewSubCycleTimeState);
                    }
                else // business as usual
                    {
                        #include "setRegionSolidFields.H"
                        #include "readSolidMultiRegionPIMPLEControls.H"
                        #include "solveSolid.H"
                    }
                // end hack
            }

5) Recompile

[mattijs]

Hi kpm,suitup

I've added the solid time step control (based on diffusion number) to 16x. The system/controlDict now has an optional 'maxDi' entry. I've updated the chtMultiRegionFoam tutorials accordingly.

The heat-flux bc is now a gradient bc and should behave much better.

Thanks,

Mattijs

[suitup]

Hey, thank you for updating.

I have tested the new version, but the usage of solidWallHeatFluxTemperature seems to be changed, I always getting FOAM FATAL IO ERROR with the new update.

My old case with e.g.:

    
platine_to_s9
{
        type            solidWallHeatFluxTemperature;
        K               K;
        q               uniform 2000;
        value           uniform 298.15;
    }

Don't work any longer, with $internalField as well. Anything missed?

I'm always getting:

--> FOAM FATAL IO ERROR: 
keyword gradient is undefined in dictionary "/home/.../test/0.1/platine/T::boundaryField::platine_to_s9"

file: /home/caelinux/.../test/0.1/platine/T::boundaryField::platine_to_s9 from line 26 to line 20.

    From function dictionary::lookupEntry(const word&, bool, bool) const
    in file db/dictionary/dictionary.C at line 396.

FOAM exiting

The usage-introduction wasn't changed under "applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperature".


Best regards.

[mattijs]

Just add a 'gradient' entry which is your starting gradient. If you have no idea what it is you can use 0.

gradient uniform 0;

[mik]

Hi all,

I have some questions related to the modification of the sofware:

So you are saying I have to add a gradient entry, so instead of:

    SOLID_to_FLUID
    {
        type            solidWallHeatFluxTemperature;
        K                 K;
        q                 uniform 25e04; // [W m⁻2]
        value          uniform 293; // initial Temp [K]
    }

I have to insert:

    SOLID_to_FLUID
    {
        type            solidWallHeatFluxTemperature;
        K                  K;
        q                  uniform 25e04; // [W m⁻2]
        value           uniform 293; // initial Temp [K]
        gradient      uniform 0;
    }

or I have to replace the entry 'q' with 'gradient', also:

    SOLID_to_FLUID
     {
         type            solidWallHeatFluxTemperature;
         K                  K;
         gradient       uniform 25e04; // [W m⁻2]
         value           uniform 293; // initial Temp [K]
     }

Because I do not inserted the entry 'gradient' using the modified solver (chtMultiRegionFoam) and the simulation run however (not sure about results)!

Then, if I modify the solver but I don't add the maxDeltaT in the controlDict could create some problems?

Finally...
I have a cubic cell (which represents the fluid region) and a cylinder in the middle (which represents the solid region), I have inserted an heat flux at the fluid bottom as well as the solid bottom (the same heat flux), yet I have to insert a properly entry for the interface (lateral surface of the cylinder (mantel)), what is the correct one?
Because there is surely an heat flux but I don't know the magnitude.. can I insert something like:

    SOLID_to_FLUID
    {     
        type            solidWallHeatFluxTemperature;
        K               K;
        q            uniform 0; // [W m⁻2]
        value           uniform 293; // initial Temp [K]
    }

or it is sufficent:

    type            zeroGradient;

(and maybe instead of 'q', 'gradient')

Many thanks

Best Regards

Michea

[mattijs]

The 'gradient' entry in the solidWallHeatFluxTemperature is needed only for startup and gets recalculated using q/K. Using q=0 should be equivalent to zeroGradient.

[mik]

Hi Mattijs,

so the right formulation is the one:

    FLUID_BOTTOM
    {
        type            solidWallHeatFluxTemperature;
        K               K;
        q               uniform 25e04; // [W m⁻2]
    value           uniform 293; // initial Temp [K]
    gradient    uniform 0;
    }

?
where FLUID_BOTTOM is the name of the patch.

Then, for the maxDeltaT:
It is not obligatory to insert also the maxDeltaT in the controlDict file if I change the solver as proposed above ?
(also beacuse if I give a maxDeltaT it is not respected by the solver...)

Many thanks

Best Regards

Michea

[maddalena]

Hi Mattijs,
one more question regarding the use of solidWallHeatFluxTemperatureCoupled BC: will the solver take into account the contribution of different fluxes? I mean:

• what the BC will do on the interface between two regions that both have a heat generation?
• what the BC will be if an adding heat flux (that is not known a priori) will appear in a interface between two regions that have a solidWallHeatFluxTemperatureCoupled applied?

Some other information are here: http://www.cfd-online.com/Forums/ope...implefoam.html
Thanks for your time.
Cheers,

maddalena

[mik]

Hi all,

I have tried to use both chtMultiRegionFoam "version", the original one and the one modified as described above (using "sub loops" for the solid part) but I don't see any difference!

Can someone please help me to understand what is the point of do this ?

It is only related to the solid region ?

Many thanks

Best Regards

Michea

[maddalena]

Originally Posted by mattijs

I've added the solid time step control (based on diffusion number) to 16x. The system/controlDict now has an optional 'maxDi' entry. I've updated the chtMultiRegionFoam tutorials accordingly.

Hi Mattijs,
I compared 1.6 and 1.6.x results on the same multiRegionHeater tutorial. While 1.6 allows you to specify a uniform ( 0 0 0) value everywhere, the 1.6.x does not. This can be verified
setting U as

    internalField   uniform (0 0 0);

and

 minX
            {
                type            fixedValue;
                value           uniform ( 0 0 0 );
            }

in the bottomAir and topAir regions of the multiRegionHeater tutorial. While the 1.6 version runs without problem, the 1.6.x version gives sigFpe.

Region: bottomAir Courant Number mean: 0 max: 0
Region: topAir Courant Number mean: 0 max: 0
Region: heater Diffusion Number mean: 2.792023e-06 max: 3.333334e-06
Region: leftSolid Diffusion Number mean: 2.78481e-06 max: 3.333334e-06
Region: rightSolid Diffusion Number mean: 2.78481e-06 max: 3.333334e-06
#0  Foam::error::printStack(Foam::Ostream&) in "/root/OpenFOAM/OpenFOAM-1.6.x/lib/linux64GccDPOpt/libOpenFOAM.so"
#1  Foam::sigFpe::sigFpeHandler(int) in "/root/OpenFOAM/OpenFOAM-1.6.x/lib/linux64GccDPOpt/libOpenFOAM.so"
#2  ?? in "/lib/libc.so.6"
#3  main in "/root/OpenFOAM/OpenFOAM-1.6.x/applications/bin/linux64GccDPOpt/chtMultiRegionFoam"
#4  __libc_start_main in "/lib/libc.so.6"
#5  _start in "/root/OpenFOAM/OpenFOAM-1.6.x/applications/bin/linux64GccDPOpt/chtMultiRegionFoam"

What I noticed is that the calculation of the first time step is changed in the two versions: while 1.6 uses setInitialDeltaT.H, the 1.6.x uses setInitialMultiRegionDeltaT.H, thus I am thinking that the problem may be there.

Could you check that?

Regards,

maddalena

[mattijs]

You are right, the time step calculation was dividing by zero (Courant number). I've pushed a fix into 17x.

Thanks,

Mattijs

[maddalena]

Originally Posted by mattijs

You are right, the time step calculation was dividing by zero (Courant number). I've pushed a fix into 17x.

Thanks,

Mattijs

Perfect! Thank you for your fast support.

mad

[crmccreary]

I am having the same problem. Using 1.7.0 with the sub-cycle patch applied. A transient finite element analysis predicts SS in ~3000s. At ~120s, the solution blows up. At blowup, the max temp in the solid is ~1degK above initial conditions whereas in the transient thermal FEA, the max temp should be ~10degK above initial.
I created a tarball of the problem and placed it on my dropbox public folder

[richpaj]

This doesn't directly follow from the above, sorry, but the title of the
original thread appears to cover the observation made below.



"solidWallMixedTemperatureCoupled"

In this boundary condition the enthalpy flux across a given
solid/fluid interface is balanced to yield an estimate for the temperature on the wall.
Currently, this is implemented with the following expression for enthalpy flux
K * snGrad(T) ......(1)
which is fine, for example in many solids, if Cp is a constant for the temperatures
under consideration. For fluids, however, it might be more appropriate to use
alphaEff * snGrad( hs ) .......(2)
where hs is the sensible enthalpy (or alternatively implemented with h, the total enthalpy).

It's true that only a slight modification in T on the interface will probably result in
most cases but there are exceptions when steep thermal gradients near to the walls
might be present, for example.

In system/solid/changeSolidDict it should be possible to effect something along the
lines of the following:


T
{
internalField uniform 300;

boundaryField
{
....

solid_to_fluid
{
type solidWallMixedTemperatureCoupled;
neighbourFieldName T;// ==> hs
K;// ==>alphat or alphaSGS or whatever the turbulence model provides
$internalField;
}

.......

}

}

I haven't tested the precise implications of the above but at first sight it seems a relatively painless way
to provide a more accurate fluid flux estimate.

Regards,

Richard K.

[richpaj]

After some testing and re-coding it was in the event possible to balance
alphaEff * snGrad( hs )
either side of the interface. In the particular case I considered the integrated enthalpy flux over the common interface yielded values that were approximately
10% different from earlier results that employed K snGrad(T). This discrepancy
though can increase dramatically if patchInternalField values of alphaEff are chosen instead of those located on the patch itself. The formulation then
appears to be highly problem/mesh-dependent.


Regards,

Richard K.