|
[Sponsors] |
(Heattransfer) Temperature boundary condition problem |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
November 23, 2009, 05:35 |
(Heattransfer) Temperature boundary condition problem
|
#1 |
Member
Join Date: Jun 2009
Location: Germany
Posts: 38
Rep Power: 17 |
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. |
|
January 15, 2010, 09:38 |
|
#2 |
New Member
kpm
Join Date: Jan 2010
Location: Germany
Posts: 9
Rep Power: 16 |
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; } 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); $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" } 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 } |
|
April 3, 2010, 06:13 |
|
#3 |
Senior Member
Mattijs Janssens
Join Date: Mar 2009
Posts: 1,419
Rep Power: 26 |
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 |
|
April 12, 2010, 15:36 |
Getting FOAM FATAL IO ERROR with new the update
|
#4 |
Member
Join Date: Jun 2009
Location: Germany
Posts: 38
Rep Power: 17 |
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; } 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 Best regards. |
|
April 13, 2010, 13:15 |
|
#5 |
Senior Member
Mattijs Janssens
Join Date: Mar 2009
Posts: 1,419
Rep Power: 26 |
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; |
|
June 24, 2010, 09:06 |
|
#6 |
Member
Michea Ferrari
Join Date: Mar 2010
Location: Switzerland
Posts: 30
Rep Power: 16 |
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] } Code:
SOLID_to_FLUID { type solidWallHeatFluxTemperature; K K; q uniform 25e04; // [W m⁻2] value uniform 293; // initial Temp [K] gradient uniform 0; } Code:
SOLID_to_FLUID { type solidWallHeatFluxTemperature; K K; gradient uniform 25e04; // [W m⁻2] value uniform 293; // initial Temp [K] } 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] } Code:
type zeroGradient; Many thanks Best Regards Michea |
|
June 25, 2010, 01:40 |
|
#7 |
Senior Member
Mattijs Janssens
Join Date: Mar 2009
Posts: 1,419
Rep Power: 26 |
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.
|
|
June 25, 2010, 04:11 |
|
#8 |
Member
Michea Ferrari
Join Date: Mar 2010
Location: Switzerland
Posts: 30
Rep Power: 16 |
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 |
|
June 25, 2010, 09:30 |
some more questions
|
#9 |
Senior Member
maddalena
Join Date: Mar 2009
Posts: 436
Rep Power: 23 |
Hi Mattijs,
one more question regarding the use of solidWallHeatFluxTemperatureCoupled BC: will the solver take into account the contribution of different fluxes? I mean:
Thanks for your time. Cheers, maddalena |
|
June 28, 2010, 14:08 |
|
#10 |
Member
Michea Ferrari
Join Date: Mar 2010
Location: Switzerland
Posts: 30
Rep Power: 16 |
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 |
|
July 7, 2010, 10:44 |
|
#11 | |
Senior Member
maddalena
Join Date: Mar 2009
Posts: 436
Rep Power: 23 |
Quote:
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); Code:
minX { type fixedValue; value uniform ( 0 0 0 ); } 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" Could you check that? Regards, maddalena |
||
July 8, 2010, 06:28 |
|
#12 |
Senior Member
Mattijs Janssens
Join Date: Mar 2009
Posts: 1,419
Rep Power: 26 |
You are right, the time step calculation was dividing by zero (Courant number). I've pushed a fix into 17x.
Thanks, Mattijs |
|
July 8, 2010, 06:37 |
|
#13 |
Senior Member
maddalena
Join Date: Mar 2009
Posts: 436
Rep Power: 23 |
||
September 17, 2010, 19:15 |
Same problem
|
#14 |
New Member
Charles McCreary
Join Date: Jun 2010
Posts: 12
Rep Power: 16 |
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 |
|
September 24, 2010, 04:08 |
solidWallMixedTemperatureCoupled in chtMultiRegionFoam
|
#15 |
Member
Richard Kenny
Join Date: Mar 2009
Posts: 64
Rep Power: 18 |
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. |
|
October 14, 2010, 23:18 |
|
#16 |
Member
Richard Kenny
Join Date: Mar 2009
Posts: 64
Rep Power: 18 |
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. |
|
Thread Tools | Search this Thread |
Display Modes | |
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
Wind turbine simulation | Saturn | CFX | 60 | July 17, 2024 06:45 |
Free Stream Temperature wall boundary condition | emanuele | FLUENT | 0 | March 19, 2007 11:45 |
Dynamic Temperature Boundary Condition | Kshitij | FLUENT | 0 | October 12, 2005 14:40 |
2D simulation as boundary condition for 3D problem | Jan De smet | FLUENT | 1 | October 11, 2002 03:02 |
a problem with Boundary condition | M Rad | Main CFD Forum | 12 | November 27, 1998 13:49 |