Multiplicative split of the deformation gradient for thermal expansion

[felixR96]

I have noticed that in simulations with thermal expansion, the deformation is not separated multiplicatively. In the ‘calcmechstrain.f’ file, the thermal strain is subtracted from the total strain to obtain the mechanical strain. One example problem has been discussed here: https://calculix.discourse.group/t/thermal-expansion-with-beam-bending/2359

I have implemented the multiplicative split and orientated myself on these papers:
http://dx.doi.org/10.2298/TAM0229379V
http://dx.doi.org/10.1016/j.ijsolstr.2014.04.010
https://doi.org/10.1016/j.finel.2018.03.001

My solution so far (see attached files):

• In resultsmech.f I calculate the deformation gradient F in case of thermal simulations (change in if-statement in line 756) and calculate thermal deformation gradient, thermal jacobian and the mechanical deformation gradient Fme = F*Fth^-1 (line 949 ff.).
• In calcmechstrain.f I pass the mechanical deformation gradient and compute the mechanical strain by Eme = 0.5*(Fme^T*Fme - I).
• Back in resultsmech.f, the stress and the tangent modulus are calculated as usual (no change in mechmodel.f). These quantities are living in the intermediate configuration now.
• After that, I perform a pull-back to the reference configuration for the stress by S = det(Fth)Fth^-1Sme*Fth^-T and a similar operation for the tangent modulus (line 1046 ff.).
• The commented-out code later, if used, (1137 ff.) does not affect the simulation. With the pull-back of the Lagrange-strain (E_mech,ref = Fth^TEmech,intFth), the thermal strain is calculated correctly since E_total = E_th + Fth^TE_mech,intFth. However, I don't use that because the mechanical strain w.r.t. the intermediate configuration is more meaningful for me than the thermal strain. An idea would be to leave E_mech untransformed and calculate the thermal strain by E_th = 0.5*(Fth^T*Fth - I).

My examples files attached:

• A cantilever (St. Venant-Kirchhoff material, Poisson ratio = 0.0) under a single load at the free end deflects with wmax = 1/3FL³/(EI) = 0.1m. With orthotropic thermal expansion (alphaTx > 0, alphaTx,y = 0) the beam lengthens by a factor of 1.4 while the cross-section remains the same. The deflection now is (1.4)^3*0.1m = 0.2744m (more information here https://calculix.discourse.group/t/thermal-expansion-with-beam-bending/2359/23). Even though I used an isotropic material with *ELASTIC, I needed the orthotropic (engineering constants) attribute, otherwise there is no convergence. I assume this is related to the fact that in linel.f for isotropy only the first two elements of elas(21) are calculated and therefore the pull-back of the tangent modulus does not work here.
• Same cantilever with Neo-Hooke material (Poisson ratio = 0.475 here). Just to check if the deflection also is increased by the factor (1.4)^3
• A cube (nearly incompressible) under thermal expansion which is constrained in y-direction (displacement = 0) and free to move in x- and z-direction with different thermal expansion coefficients. This works more or less fine, check my attached PDF file with the analytical solution. strains, stresses, displacements, these are correct. I explain minor deviations by the fact that my material is not perfectly incompressible.

What I did not do:

• I have only analysed *STATIC simulations with prescribed temperature field. All other simulation types has not been analysed. Especially not those ones where time derivatives play a role.
• I did not modify other code files where it should be necessary. For example, I added the mechanical deformation gradient as an argument to calcmechstrain, but only adjusted the code in resultsmech.f where the subroutine "calcmechstrain" is called, not in other files like objective_shapeener_dx.f.
• I have not modified the calculation of the thermal strain as an output for the frd-file or dat-file (here E_thermal could be calculated directly from F_thermal)
• I have rearranged something in the file for the calculation of the mechanical strain. I don't know whether this has consequences for certain types of simulation.

I am providing my solution approach, even if it is certainly not complete:
2024-09-16_Github_4.tar.gz.zip

Best regards,

Felix

[edit] I solved an issue in my cube example. I used my own UMAT for a nearly incompressible Neo-Hooke material. Some problems arose in the numerics so I decided to replace it with a simpler example. The implemented Neo-Hooke material is used and to make the example easier I used simpler expansion coefficients and just one hexaeder element.

[felixR96]

I have replaced the appendix because I have solved the problems with the cube example. The source code is the same. It was due to my model or rather the UMAT. The example is now simplified without Umat. Please refer to the current appendix and ignore the previous version.