CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

Thermal stress prediction with elasticThermalSolidFoam - No contribution to UEqn?

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

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   November 10, 2014, 12:28
Default Thermal stress prediction with elasticThermalSolidFoam - No contribution to UEqn?
  #1
Member
 
A. Bernath
Join Date: Jun 2011
Location: Karlsruhe, Germany
Posts: 39
Rep Power: 10
derkermit is on a distinguished road
Hello all,
I'm trying to predict thermal stresses in a polymer that cures in a mold. At the same time I model the same problem with FEM (Abaqus) in order to compare the results.

In order to simplify my problem, I started with a cube with 20x20x20 elements and 0.1m edge length. Only thermal expansion is taken into account at the moment.

In the source code of the solver one can find the following lines:

Code:
    divSigmaExp = fvc::div
      (
       mu*gradU.T() + lambda*(I*tr(gradU)) - (mu + lambda)*gradU,
       "div(sigma)"
       )
      - gradThreeKalphaDeltaT;
This is part of the source term of the momentum equation an the thing I don't understand is the behaviour of gradThreeKalphaDeltaT in the presence of homogenous temperature fields. As an example, consider the cubes temperature to be 323.15 K and the specified reference temperatur T0 = 373.15 K. This will give a non zero threeKalphaDeltaT but as soon as the gradient of this is calculated and put in the momentum equation, the thermal contribution is just zero, at least for my understanding. As a consquence, the cube doesn't deform as it should and is only showing very small values of epsilon. From my calculations, it should reach epsilon_xx = -0.00275 of contraction (which is confirmed by abaqus) but elasticThermalSolidFoam only gives around 2e-5.
Does anyone have an idea what is wrong here?
I've attached the case to this post.

Thanks in advance for any help!
Alex
Attached Files
File Type: zip elasticThermalSolidFoam_dT-50K.zip (7.2 KB, 15 views)
derkermit is offline   Reply With Quote

Old   November 15, 2014, 07:06
Default
  #2
Super Moderator
 
bigphil's Avatar
 
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 782
Rep Power: 26
bigphil will become famous soon enoughbigphil will become famous soon enough
Hi,

I found a bug that was introduced in $FOAM_SRC/solidModels/constitutiveModel/tractionBoundaryGradient/tractionBoundaryGradient.C in foam-extend-3.1:

Change this line:
Quote:
if (!incremental)
{
const fvPatchScalarField& DT =
patch.lookupPatchField<volScalarField, scalar>("DT");

gradient += n*threeKalpha*DT;
}
to this line:
Quote:
if (incremental) // bugfix philipc Nov-14 remove not
{
const fvPatchScalarField& DT =
patch.lookupPatchField<volScalarField, scalar>("DT");

gradient += n*threeKalpha*DT;
}
This bug meant that the traction boundary conditions were incorrect with thermal stresses.

I tried your case and it seems to work now.

One comment on your case: in general I have found that the solution tolerance for U should typically be at least 1e-6, in the solidMechanics dict in system/fvSolution:
Quote:
solidMechanics
{
nCorrectors 10000;
U 1e-6;
//U 1e-3;
}
Let me know how it goes,
Philip


EDIT:
Ah I see you already found this problem and reported it:
http://sourceforge.net/p/openfoam-ex...ndrelease/258/

Last edited by bigphil; November 15, 2014 at 07:48. Reason: Just found bug report
bigphil is offline   Reply With Quote

Old   November 15, 2014, 11:35
Default
  #3
Member
 
A. Bernath
Join Date: Jun 2011
Location: Karlsruhe, Germany
Posts: 39
Rep Power: 10
derkermit is on a distinguished road
Hi Philip,
thanks for your answer. Your right, it works now, at least for the isotropic case =)

I also added TEqn to your orthotropic solver (elasticOrthoSolidFoam) and modified divSigmaExp like this:

Code:
if(divSigmaExpMethod == "standard")
  {
    //- calculating the full gradient has good convergence and no high freq oscillations
    divSigmaExp =
      fvc::div(
           (C && (symm(gradU) - alphaOrtho*(T-T0)))
           - (K & gradU),
           "div(sigma)"
           );
  }
tractionBoudnaryGradient.C has also been modified (http://www.cfd-online.com/Forums/ope...ent-class.html) in order to account for the thermal effects.
However, the solver doesn't converge. I attached the case to this post. Do you have an idea what could be the problem here?

Quote:
Originally Posted by bigphil View Post
One comment on your case: in general I have found that the solution tolerance for U should typically be at least 1e-6, in the solidMechanics dict in system/fvSolution:
I tried different values for this and didn't got that much of a difference besides a much faster simulation. I also searched for the criteria abaqus uses and this seems to be about 1e-3. But I have to mention that this was solely on simple meshes/geometries. Maybe on more complex meshes this affects the solution accuracy more. How did you figure out this specific value of 1e-6?

Have a nice weekend!
Alex
Attached Files
File Type: zip elasticOrthoThermalSolidFoam_dT-50K.zip (11.0 KB, 11 views)
derkermit is offline   Reply With Quote

Old   November 16, 2014, 07:14
Default
  #4
Super Moderator
 
bigphil's Avatar
 
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 782
Rep Power: 26
bigphil will become famous soon enoughbigphil will become famous soon enough
I will have a look at the addition of an orthotropic thermal term.

Quote:
Originally Posted by derkermit View Post
I tried different values for this and didn't got that much of a difference besides a much faster simulation. I also searched for the criteria abaqus uses and this seems to be about 1e-3. But I have to mention that this was solely on simple meshes/geometries. Maybe on more complex meshes this affects the solution accuracy more. How did you figure out this specific value of 1e-6?
This is just based on what I found works in general; but yep for a given case you just have to check that the tolerance is low enough that the results are independent of it.

Philip
bigphil is offline   Reply With Quote

Old   November 26, 2014, 07:07
Default
  #5
Member
 
A. Bernath
Join Date: Jun 2011
Location: Karlsruhe, Germany
Posts: 39
Rep Power: 10
derkermit is on a distinguished road
Hi Philip,
I'm facing another problem with the thermal part of the solver.
If there is only very small change in temperature within the domain, the TEqn doesn't converge even after 10000 iterations. However, if temperature changes occur, it converges really fast.
Did you also observe this behavior?
Is there a reason why TEqn is first solved and then T (not TEqn) is relaxed? I looked into the fluid solver codes and they always relax TEqn first and then solve it. However, modifying the solver that way doesn't help so it has to be something else.

Thanks,
Alex
derkermit is offline   Reply With Quote

Old   December 11, 2014, 06:53
Default
  #6
Super Moderator
 
bigphil's Avatar
 
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 782
Rep Power: 26
bigphil will become famous soon enoughbigphil will become famous soon enough
Yes I've noticed this too.
A very small amount of equation under-relaxation will fix this.

Put this line after the TEqn is created (and before it is solved):
Code:
TEqn.relax(TEqnRelaxFac);
And you read in the relaxation factor somewhere at the start of the solver, like this:
Code:
scalar TEqnRelaxFac = -1;
if (mesh.solutionDict().relax("TEqn"))
{
    TEqnRelaxFac = mesh.solutionDict().relaxationFactor("TEqn");
}
Then in your case system/fvSolution, add this:
Code:
relaxationFactors
{
    TEqn 0.999;
}
Note: the equation convergence is very sensitive to this factor, I have found a big difference between 0.99 and 0.999.

Philip
bigphil is offline   Reply With Quote

Old   December 15, 2014, 18:12
Default
  #7
Member
 
A. Bernath
Join Date: Jun 2011
Location: Karlsruhe, Germany
Posts: 39
Rep Power: 10
derkermit is on a distinguished road
Hi Philip,
thanks, will give it a try
Is it possible to under-relaxate UEqn/DUEqn too or would that introduce some error?
In CFD it's applied to both p and U so maybe this could be of advantage here too.

Regards, Alex
derkermit is offline   Reply With Quote

Old   December 16, 2014, 12:21
Default
  #8
Super Moderator
 
bigphil's Avatar
 
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 782
Rep Power: 26
bigphil will become famous soon enoughbigphil will become famous soon enough
Quote:
Originally Posted by derkermit View Post
Is it possible to under-relaxate UEqn/DUEqn too or would that introduce some error?
Yes, this is allowed as outer iterations are performed within each time-step, and this can be very helpful in some situations.
Typically only a very small amount is needed, between 0.99 and 0.999 often works well.

Philip
bigphil is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
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 Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Thermal and stress modelling of a dam quimperval OpenFOAM Running, Solving & CFD 0 August 22, 2014 12:08
CFX-Structural thermal stress analysis Ravikiran CFX 1 February 22, 2010 09:57
Transient Thermal Stress Distribution Ijaz CFX 0 February 19, 2009 17:27
thermal stress analysis? Jinfeng Main CFD Forum 1 January 21, 2009 09:26
Info: Short Course On Thermal Design of Electronic Equipment Arnold Free Main CFD Forum 0 August 10, 1999 11:18


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