CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Bugs

(Heattransfer) Temperature boundary condition problem

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   November 23, 2009, 05:35
Default (Heattransfer) Temperature boundary condition problem
  #1
Member
 
Join Date: Jun 2009
Location: Germany
Posts: 38
Rep Power: 8
suitup is on a distinguished road
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.
Attached Files
File Type: gz zHeatsink.tar.gz (11.7 KB, 85 views)
suitup is offline   Reply With Quote

Old   January 15, 2010, 09:38
Default
  #2
kpm
New Member
 
kpm
Join Date: Jan 2010
Location: Germany
Posts: 9
Rep Power: 7
kpm is on a distinguished road
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

Code:
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:
Code:
    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:

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

Code:
            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
kpm is offline   Reply With Quote

Old   April 3, 2010, 05:13
Default
  #3
Super Moderator
 
Mattijs Janssens
Join Date: Mar 2009
Posts: 1,416
Rep Power: 16
mattijs is on a distinguished road
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
mattijs is offline   Reply With Quote

Old   April 12, 2010, 14:36
Default Getting FOAM FATAL IO ERROR with new the update
  #4
Member
 
Join Date: Jun 2009
Location: Germany
Posts: 38
Rep Power: 8
suitup is on a distinguished road
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.:

Code:
    
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:

Code:
--> 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.
suitup is offline   Reply With Quote

Old   April 13, 2010, 12:15
Default
  #5
Super Moderator
 
Mattijs Janssens
Join Date: Mar 2009
Posts: 1,416
Rep Power: 16
mattijs is on a distinguished road
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;
mattijs is offline   Reply With Quote

Old   June 24, 2010, 08:06
Default
  #6
mik
Member
 
Michea Ferrari
Join Date: Mar 2010
Location: Switzerland
Posts: 30
Rep Power: 7
mik is on a distinguished road
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:
Code:
    SOLID_to_FLUID
    {
        type            solidWallHeatFluxTemperature;
        K                 K;
        q                 uniform 25e04; // [W m⁻2]
        value          uniform 293; // initial Temp [K]
    }
I have to insert:
Code:
    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:
Code:
    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:
Code:
    SOLID_to_FLUID
    {     
        type            solidWallHeatFluxTemperature;
        K               K;
        q            uniform 0; // [W m⁻2]
        value           uniform 293; // initial Temp [K]
    }
or it is sufficent:
Code:
    type            zeroGradient;
(and maybe instead of 'q', 'gradient')

Many thanks

Best Regards

Michea
mik is offline   Reply With Quote

Old   June 25, 2010, 00:40
Default
  #7
Super Moderator
 
Mattijs Janssens
Join Date: Mar 2009
Posts: 1,416
Rep Power: 16
mattijs is on a distinguished road
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.
mattijs is offline   Reply With Quote

Old   June 25, 2010, 03:11
Default
  #8
mik
Member
 
Michea Ferrari
Join Date: Mar 2010
Location: Switzerland
Posts: 30
Rep Power: 7
mik is on a distinguished road
Hi Mattijs,

so the right formulation is the one:
Code:
    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
mik is offline   Reply With Quote

Old   June 25, 2010, 08:30
Default some more questions
  #9
Senior Member
 
maddalena's Avatar
 
maddalena
Join Date: Mar 2009
Posts: 436
Rep Power: 12
maddalena is on a distinguished road
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: solidWallHeatFluxTemperature at the solid solid interface in chtMultiRegionSimpleFoam
Thanks for your time.
Cheers,

maddalena
maddalena is offline   Reply With Quote

Old   June 28, 2010, 13:08
Default
  #10
mik
Member
 
Michea Ferrari
Join Date: Mar 2010
Location: Switzerland
Posts: 30
Rep Power: 7
mik is on a distinguished road
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
mik is offline   Reply With Quote

Old   July 7, 2010, 09:44
Default
  #11
Senior Member
 
maddalena's Avatar
 
maddalena
Join Date: Mar 2009
Posts: 436
Rep Power: 12
maddalena is on a distinguished road
Quote:
Originally Posted by mattijs View Post
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
Code:
    internalField   uniform (0 0 0);
and
Code:
 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.
Code:
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
maddalena is offline   Reply With Quote

Old   July 8, 2010, 05:28
Default
  #12
Super Moderator
 
Mattijs Janssens
Join Date: Mar 2009
Posts: 1,416
Rep Power: 16
mattijs is on a distinguished road
You are right, the time step calculation was dividing by zero (Courant number). I've pushed a fix into 17x.

Thanks,

Mattijs
mattijs is offline   Reply With Quote

Old   July 8, 2010, 05:37
Default
  #13
Senior Member
 
maddalena's Avatar
 
maddalena
Join Date: Mar 2009
Posts: 436
Rep Power: 12
maddalena is on a distinguished road
Quote:
Originally Posted by mattijs View Post
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
maddalena is offline   Reply With Quote

Old   September 17, 2010, 18:15
Default Same problem
  #14
New Member
 
Charles McCreary
Join Date: Jun 2010
Posts: 9
Rep Power: 7
crmccreary is on a distinguished road
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
crmccreary is offline   Reply With Quote

Old   September 24, 2010, 03:08
Default solidWallMixedTemperatureCoupled in chtMultiRegionFoam
  #15
Member
 
Richard Kenny
Join Date: Mar 2009
Posts: 59
Rep Power: 9
richpaj is on a distinguished road
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 is offline   Reply With Quote

Old   October 14, 2010, 22:18
Default
  #16
Member
 
Richard Kenny
Join Date: Mar 2009
Posts: 59
Rep Power: 9
richpaj is on a distinguished road
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.
richpaj is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Wind turbine simulation Saturn CFX 34 October 16, 2014 05:27
Free Stream Temperature wall boundary condition emanuele FLUENT 0 March 19, 2007 11:45
Dynamic Temperature Boundary Condition Kshitij FLUENT 0 October 12, 2005 13:40
2D simulation as boundary condition for 3D problem Jan De smet FLUENT 1 October 11, 2002 02:02
a problem with Boundary condition M Rad Main CFD Forum 12 November 27, 1998 13:49


All times are GMT -4. The time now is 04:42.