Thermal expansion with beam bending

[felixR]

Hello everyone,

I am trying to simulate the interaction of a beam under bending while the whole body (every node) is subjected to a (unidirectional/orthotropic) thermal expansion:
A cantilever with a fixed support on the left side is subjected to a force perpendicular to the beam axis at the free end on the right side. The maximum displacement occurs at the free end.

Additionally, a thermal expansion is applied, which makes the beam longer (only in longitudinal direction so the cross-section remains constant) and therefore the deflection should be higher compared to a purely mechanical simulation.

The purely mechanical simulation without thermal expansion works fine and the maximum displacement corresponds to the analytical solution

wmax = 1/3*F*L³/(EI) = 1/3*(1050000N)*(2.00m)³/(210000e6N/m²*1/7500 m⁴) 
     = 0.1 m
with a rectangular cross-section b/h = 0.2/0.2 m

If I apply a thermal expansion and the length of the beam becomes Lth = 2.8m, the maximum deflection should be

wmax = 0.1m*(2.8m)³/(2.0m)³
     = 0.2744 m

but I still receive the same deflection of 0.1m.

I tried different simulation types (*Static, *Uncoupled Temperature-Displacement, *Coupled Temperature-Displacement) and different ways to increase the temperature (*Initial Conditions=Temperature with *Temperature or *Boundary DOF 11). Sometimes I received error messages, that either the initial temperature or the final temperature in a node is not specified.

Here is one input file of a working simulation, even though it took 3 hours to finish.

************************** Mesh ***************************
** 
*INCLUDE, INPUT= ./all.msh
**
**
****************** Node and Element Sets ******************
**
**In ./all.msh
**
**
************************* Materials ************************
**
*Material, Name=DummyMat
*Elastic
210000e6, 0.1
*Expansion, Zero=273.0,TYPE=ORTHO
0.005,0.0,0.0
*DENSITY
7850.0
*SPECIFIC HEAT
500.0
**
**
************************* Sections *************************
**
*Solid section, Elset=Eall, Material=DummyMat
**
**
*************************** Step ***************************
**
**Time Points
*TIME POINTS,NAME=T1,GENERATE 
0.0,10.0,1.0
**
**
**
** Supports
*Boundary,fixed
Nsurfleft,1,3
**
** Initial Temperature
*INITIAL CONDITIONS,TYPE=TEMPERATURE
Nall,273.
**
**
*Step,NLGEOM=YES,INC=10000000000
*COUPLED TEMPERATURE-DISPLACEMENT
1.0e-2,10.0
**
**
*Cload
Nsurfright,3,-16150
**
**
**Temperature
**Nall,353.0
**
*BOUNDARY
Nall,11,11,353.
**
** Output .frd
*Node file,TIMEPOINTS=T1
RF, U, NT
*El file,TIMEPOINTS=T1
S, E, NOE, ENER
**
** Outputs .dat
*NODE PRINT, NSET=Nsurfright, TOTALS=ONLY,TIMEPOINTS=T1
RF
*NODE PRINT,NSET=Nall,TIMEPOINTS=T1
U
**
*End step
**

I used C3D20R elements. The cload is applied to 65 nodes at the free end.

Any idea how I can adjust the simulation for increased deflection due to thermal expansion?

Thanks,

Felix

[Calc_em]

Try with 2 subsequent steps. Not necessarily coupled temperature-displacement, it might be sufficient to use static steps with *TEMPERATURE (prescribed temperature field).

[felixR]

Hey Calc_em,

thank you for your reply.
I changed it to a static simulation. Now it only takes a few seconds to finish.

a) First I tried it to keep the thermal expansion and the cload in one step, modified the final temperature in several simulations and observed the following behavior:

While the displacement and the mechanical strain in x-direction are the same in all simulations and only correct if L does not change (no thermal expansion), The stress in the longitudinal direction changes depending on the temperature expansion. It also corresponds to the analytical solution: σ_xx = M/W = F*L_th/(0.2^3/6) with L_th as the length of the stretched beam.

(I set the Poisson number to 0.0 for a clean stress output at the fixed support)

It seems that the stresses are calculated correctly, but the strain and displacement are not.

Then I tried to split the thermal expansion and the cload into two steps.
b) With the *Temperature in the first step the simulation does not converge. I just removed the *cload card from that step.

c) If I switch the order and apply the force first, the simulation runs successfully to the end. However, the displacement remains at 0.1m.

My three simulations (Temperature and cload in one step and in two steps with different order) are uploaded here: bwSync&Share

It would be interesting to know why the simulation does not converge if the *cload card is removed and why the stresses are correct but the displacement field is not.

Thank you very much

Felix

[Calc_em]

Try adding some plasticity to make this deformation from the first step (due to thermal expansion) permanent. Then the second step should take it into account. Nlgeom is enabled all the time, right ?

[JuanP74]

@felixR

*Expansion, Zero=273.0,TYPE=ORTHO
0.005,0.0,0.0

Is this correct? If I look at the example in the manual:

Example:
*EXPANSION,ZERO=273.
12.E-6,373.
20.E-6,573.

[Calc_em]

This example is isotropic temperature-dependent thermal expansion. Here it’s an orthotropic one.

[JuanP74]

right, so where is the temperature in the input deck above? that’s what I meant.

[Calc_em]

There’s no temperature dependency of the thermal expansion coefficient(s) in it. Only the reference (zero strain) temperature. Temperature dependency is optional.

[vicmw]

This does look worrying. I reproduced it starting from scratch and get the same effect. Here the length doubles due to thermal expansion but lateral displacement is the same.

Correct stress but wrong displacement is a sign of incorrect stiffness. It seems like the stiffness matrix somehow isn’t being updated correctly with the length change.

I wonder if this is isolated to thermal expansion. I would like to increase the length some other way. Just applying a tensile force would bury the effect in increased bending stiffness due to that load.

[vicmw]

It would be interesting to know why the simulation does not converge if the *cload card is removed

That’s a thing that tends to happen, and my model does that too. You can put a tiny dummy load anywhere and it gets it going. I don’t think it’s related to the main problem. My guess is something to do with a force convergence criterion when there are only zero forces everywhere.

[JuanP74]

@vicmw

Correct stress but wrong displacement is a sign of incorrect stiffness.

but stresses are calculated from displacements in FEM

[vicmw]

Yea but displacements are also calculated from stiffness and that cancels out scaling of stiffness. If you set the wrong Young’s modulus in linear analysis, stress still comes out right under force load.

[Calc_em]

In Abaqus, when using 2 steps and *TEMPERATURE, I get the following deflection:

• without NLGEOM: 101 mm (the length of the beam is 2800 mm, as predicted analytically)
• with NLGEOM: 332 mm (the length of the beam reaches 2985 mm)

Submitting the same nonlinear input file in CalculiX results in non-convergence so I can’t compare them directly.

[Disla]

Shouldn’t the material definition be orthotropic too when Ortho is used for the expansion properties.?

@vicmw

Correct stress but wrong displacement

If I use isotropic material my displacements are as expected but the stresses are suspicious. I would expect a beam free to expand to end up with zero stresses and it ends up with some VM Stress.

[vicmw]

If I use isotropic material my displacements are as expected but the stresses are suspicious.

I get zero stress (~1e-5) with this material, even at temperature that doubles the length:

*ELASTIC,TYPE=ISOTROPIC
200e9,0
*EXPANSION,TYPE=ORTHO
0.01,0,0

[Calc_em]

@Disla

Shouldn’t the material definition be orthotropic too when Ortho is used for the expansion properties.?

I don’t think so, it should be independent. At least theoretically. Here’s an input deck from Abaqus example adapted for ccx:

*Node
1, 0, 0, 0
2, 1, 0, 0
3, 1, 1, 0
4, 0, 1, 0
5, 0, 0, 1
6, 1, 0, 1
7, 1, 1, 1
8, 0, 1, 1
*Element, type=C3D8R
1, 1, 2, 3, 4, 5, 6, 7, 8
*Nset, nset=IC3D8R_NALL, generate
1,  8,  1
*Elset, elset=IC3D8R_EALL
1,
*Elset, elset=IC3D8R_Set-3
1,
*Orientation, name=IC3D8R-Ori-1
1,0,0,0,1,0
3,0
*Solid Section, elset=IC3D8R_Set-3, orientation=IC3D8R-Ori-1, material=MELAPLA
,
*Elset, elset=EALL
1,
*Nset, nset=NALL, generate
1,  8,  1
*Nset, nset=NFIX13, generate
1,  8,  1
*Material, name=MELAPLA
*Density
8032,
*Elastic
1.931e+11, 0.3
*Expansion, type=ORTHO
1e-05,   5e-06, 2.5e-06
*Plastic
206893,0
413786,1
*Boundary
NFIX13, 1, 1
NFIX13, 2, 2
NFIX13, 3, 3
*Initial Conditions, type=TEMPERATURE
NALL, 0
*Step, nlgeom=YES, inc=5000
*Static
0.1, 1, 1e-05, 1
*Temperature, op=NEW
NALL, 100
*Node file
U,NT
*El file
S,E,PE,THE
*End Step

I do get different stresses than in Abaqus.

[Disla]

@vicmw

I get zero stress (~1e-5) with this material

I wouldn’t want to add more confusion to the post. I think that can be a different post if you both don’t mind.
The reason I suspect , could be that we are considering different types of analysis and I think it would be better to deal with this separately.

[Disla]

Looks like ORTHO deactivates NLGEOM for both , mechanical and thermal geometry update if they are computed simultaneously. With isotropic expansion works fine but it’s a slightly different problem. If that’s acceptable, results agree with expected result.

[felixR]

Hello everyone,

I am reading your replies while trying out the ideas.
With a hyperelastic material (Neo-Hooke) or with a UMAT (Ciarlet) I observe the same behavior of a constant deflection and a varying stress, depending on the thermal expansion.

@Disla

Looks like ORTHO deactivates NLGEOM for both , mechanical and thermal geometry update if they are computed simultaneously

According to the terminal output, my simulations were always with nlgeom. When I start two simulations with a) isotropic th. exp. and b) orthotropic th. exp. but with the same value in all three directions, the results are equal. So basically

*Expansion, Zero=273.0,TYPE=ORTHO
0.005,0.005,0.005

Is it the same in your calculation?
Maybe my choice of parameters is a bit inconvenient since an increased length with a scaled cross-section always provide a deflection of approx. 100 mm.

Right now I am trying to get set up a simulation with some plasticity as recommended here

@Calc_em

Try adding some plasticity to make this deformation from the first step (due to thermal expansion) permanent.

however, until now there is no convergence.

If I go to equation 10 in the manual and look at the calculation of the stresses for the Neo-Hooke material in rubber.f, then the mechanical strain is simply used for the stress calculation. A thermal Jacobi determinant is also defined there and this is also used in some equations, but it is set to 1.0.

I am still trying to find out why stresses and displacements are not consistent with each other in my nonlinear simulations.

[Disla]

This is my approach with isotropic expansion. Comparision with the expected value is above.

Bending_And_Thermal_Expansion.inp

Note that just with thermal expansion (if you remove load) , displacement results match with the Abaqus linear solution. That’s why we suspect nonlinear is not active for thermal expansion.
Abaqus reports .49m of displacement to each direction with 0.98m overall length increase for the nonlinear process but that’s not a good reference to compare with ccx as they use different formulation.(SEE NOTE)

Ccx Strains without load shows an homogeneous S11=0.48 across the bar without load and zero Stresses.That means it’s using a Lagrange (Green) Strain formulation.
Why do you want to include plasticity in the model and add more uncertainties before surfing all doubts?

NOTE: Don Guido commented about discrepancies between Abaqus and ccx with expansion an NLGEOM. Take a look at this :

Comparison CalculiX-Abaqus with NLGEOM

Dear CalculiX users, in the past differences were reported in the results between CalculiX and Abaqus in the presence of NLGEOM on the *STEP card. Further analysis has revealed at least 3 sources for these differences: In Abaqus the logarithmic strain is used (unless a hyperelastic material is selected). This implies that also thermal strains, if present, are logarithmic, i.e. the thermal expansion coefficient must have been determined based on the actual length of the specimen, not the initi…