Coupled plasticity hydrogen diffusion problem

Hellos,

Has anyone tried to use Calculix for solving a coupled plasticity hydrogen diffusion problem using a UMAT+UMATHT?
I am trying to port the one I have from Abaqus to Calculix. I was able to port the UMAT and produce comparable results. But the moment I use a *Coupled Temperature-displacement with UMATHT, the solver crashes. I get abnormally high strain increments in the first increment in to UMAT, even when I set the flux in UMATHT to zero, without any thermal boundary condition. I also tried using the inbuilt plasticity model with my UMATHT with the same results.

Has anyone tried to use the analogy between the heat equation and the mass diffusion equations successfully, as done in Abaqus?

Thanks!

1 Like

Did you try running only the thermal (or in fact diffusion mass transfer) analysis in CalculiX ? Does it work with your subroutine when there’s no coupling with solid mechanics ? You could also compare the results with built-in diffusion mass transfer procedure to see if it’s working correctly.

True! I will give it a try tomorrow.
But, the UMATHT works fine in Abaqus anyway. So the difference should arise from the way the coupled solver is implemented in Calculix or I am missing some key definition in the input file.
I do understand that Calculix uses a ‘staggered’ technique. I suppose it is thermal followed by mechanical. It is strange though that the strain increment is so high after the first thermal calculation. No thermal expansion is defined, I wonder what part leads to a ‘thermal’ strain, and that high.

Which type of thermo-mechanical analysis step are you using - *COUPLED TEMPERATURE-DISPLACEMENT or *UNCOUPLED TEMPERATURE-DISPLACEMENT ?

I use *COUPLED TEMPERATURE DISPLACEMENT.
It was probably confusing, but I meant staggered scheme in each iteration. Isn’t that how the coupled solver in Calculix work?

Coupled step means that thermal and mechanical analyses are solved simultaneously (two-way coupling) while uncoupled step means that they are solved sequentially (thermal followed by mechanical in each increment and thus no feedback from mechanical to thermal within the same increment).

Yes you are right, but the manual says ‘the coupling is not done at the matrix level and only by updating the boundary condition after each iteration’. Which could be the reason why one doesn’t have to define the Internal thermal energy and its variations in UMATHT.
This coupling is a bit confusing for me, and I don’t seem to get enough information in the manual.

So I decided to test the coupled solver with in-built model and compare it to Abaqus.
I consider a simple cube with 8 C3D20R(T) elements, subjected to uniaxial displacement with 0 initial temperature and no other thermal loading → Purely mechanical loading. So temperature should remain zero.
Below is a comparison of displacements in the model between calculix *STATIC and *COUPLED Temp-Disp to the Abaqus *Coupled Temp-Disp. One notices a strange displacement behaviour in the coupled calculix case. As expected the temperatures however remain 0 in both solvers.
It could be a bug, or I understood the concept of coupled solver in calculix totally wrong.


Would be great to know if anyone else noticed such behaviour in ccx2.20.
The forum doesn’t allow me to attach an input file. Below is an excerpt:

** 
*NODE, NSET=NALL
      1,           1.,           1.,           1.
      2,           1.,          0.5,           1.
      3,           1.,           0.,           1.
      4,           1.,           1.,          0.5
      5,           1.,          0.5,          0.5
      6,           1.,           0.,          0.5
      7,           1.,           1.,           0.
      8,           1.,          0.5,           0.
      9,           1.,           0.,           0.
     10,          0.5,           1.,           1.
     11,          0.5,          0.5,           1.
     12,          0.5,           0.,           1.
     13,          0.5,           1.,          0.5
     14,          0.5,          0.5,          0.5
     15,          0.5,           0.,          0.5
     16,          0.5,           1.,           0.
     17,          0.5,          0.5,           0.
     18,          0.5,           0.,           0.
     19,           0.,           1.,           1.
     20,           0.,          0.5,           1.
     21,           0.,           0.,           1.
     22,           0.,           1.,          0.5
     23,           0.,          0.5,          0.5
     24,           0.,           0.,          0.5
     25,           0.,           1.,           0.
     26,           0.,          0.5,           0.
     27,           0.,           0.,           0.
     28,          0.5,           1.,         0.75
     29,          0.5,         0.75,          0.5
     30,          0.5,          0.5,         0.75
     31,          0.5,         0.75,           1.
     32,           1.,         0.75,           1.
     33,           1.,          0.5,         0.75
     34,           1.,         0.75,          0.5
     35,           1.,           1.,         0.75
     36,         0.75,          0.5,           1.
     37,         0.75,           1.,           1.
     38,         0.75,          0.5,          0.5
     39,         0.75,           1.,          0.5
     40,          0.5,         0.25,          0.5
     41,          0.5,           0.,         0.75
     42,          0.5,         0.25,           1.
     43,           1.,         0.25,           1.
     44,           1.,           0.,         0.75
     45,           1.,         0.25,          0.5
     46,         0.75,           0.,           1.
     47,         0.75,           0.,          0.5
     48,          0.5,           1.,         0.25
     49,          0.5,         0.75,           0.
     50,          0.5,          0.5,         0.25
     51,           1.,          0.5,         0.25
     52,           1.,         0.75,           0.
     53,           1.,           1.,         0.25
     54,         0.75,          0.5,           0.
     55,         0.75,           1.,           0.
     56,          0.5,         0.25,           0.
     57,          0.5,           0.,         0.25
     58,           1.,           0.,         0.25
     59,           1.,         0.25,           0.
     60,         0.75,           0.,           0.
     61,           0.,           1.,         0.75
     62,           0.,         0.75,          0.5
     63,           0.,          0.5,         0.75
     64,           0.,         0.75,           1.
     65,         0.25,          0.5,           1.
     66,         0.25,           1.,           1.
     67,         0.25,          0.5,          0.5
     68,         0.25,           1.,          0.5
     69,           0.,         0.25,          0.5
     70,           0.,           0.,         0.75
     71,           0.,         0.25,           1.
     72,         0.25,           0.,           1.
     73,         0.25,           0.,          0.5
     74,           0.,           1.,         0.25
     75,           0.,         0.75,           0.
     76,           0.,          0.5,         0.25
     77,         0.25,          0.5,           0.
     78,         0.25,           1.,           0.
     79,           0.,         0.25,           0.
     80,           0.,           0.,         0.25
     81,         0.25,           0.,           0.
*ELEMENT, TYPE=C3D20R, ELSET=EALL
1, 10, 11, 14, 13,  1,  2,  5,  4, 31, 30, 29, 28, 32, 33, 34,
   35, 37, 36, 38, 39
2, 11, 12, 15, 14,  2,  3,  6,  5, 42, 41, 40, 30, 43, 44, 45,
   33, 36, 46, 47, 38
3, 13, 14, 17, 16,  4,  5,  8,  7, 29, 50, 49, 48, 34, 51, 52,
   53, 39, 38, 54, 55
4, 14, 15, 18, 17,  5,  6,  9,  8, 40, 57, 56, 50, 45, 58, 59,
   51, 38, 47, 60, 54
5, 19, 20, 23, 22, 10, 11, 14, 13, 64, 63, 62, 61, 31, 30, 29,
   28, 66, 65, 67, 68
6, 20, 21, 24, 23, 11, 12, 15, 14, 71, 70, 69, 63, 42, 41, 40,
   30, 65, 72, 73, 67
7, 22, 23, 26, 25, 13, 14, 17, 16, 62, 76, 75, 74, 29, 50, 49,
   48, 68, 67, 77, 78
8, 23, 24, 27, 26, 14, 15, 18, 17, 69, 80, 79, 76, 40, 57, 56,
   50, 67, 73, 81, 77
*NSET, NSET=N-XX1
 19, 20, 21, 22, 23, 24, 25, 26, 27, 61, 62, 63, 64, 69, 70, 71
 74, 75, 76, 79, 80
*NSET, NSET=N-XX2
  1,  2,  3,  4,  5,  6,  7,  8,  9, 32, 33, 34, 35, 43, 44, 45
 51, 52, 53, 58, 59
*NSET, NSET=N-YY1
  3,  6,  9, 12, 15, 18, 21, 24, 27, 41, 44, 46, 47, 57, 58, 60
 70, 72, 73, 80, 81
*NSET, NSET=N-YY2
  1,  4,  7, 10, 13, 16, 19, 22, 25, 28, 35, 37, 39, 48, 53, 55
 61, 66, 68, 74, 78
*NSET, NSET=N-ZZ1
  7,  8,  9, 16, 17, 18, 25, 26, 27, 49, 52, 54, 55, 56, 59, 60
 75, 77, 78, 79, 81
*NSET, NSET=N-ZZ2
  1,  2,  3, 10, 11, 12, 19, 20, 21, 31, 32, 36, 37, 42, 43, 46
 64, 65, 66, 71, 72
**
*MATERIAL, NAME=vonMiese1
*DENSITY
1.0000e+00 
*SPECIFIC HEAT
1.0000e+00
*ELASTIC,
210000.0,0.3
*PLASTIC
400.0,0.0
800.0,1.0
*CONDUCTIVITY
1.0
**
*SOLID SECTION, ELSET=EALL, MATERIAL=vonMiese1 
*TIME POINTS,NAME=T1,GENERATE
0.,1.0,0.1
*Initial Conditions, type=TEMPERATURE
NALL, 0.0
** STEP: Load
*Step, nlgeom=YES, inc=100000
*COUPLED TEMPERATURE-DISPLACEMENT
**STATIC
0.001, 1.0, 1e-09, 1.0
**
*BOUNDARY
N-XX1, 1, 1
*BOUNDARY
N-YY1, 2, 2
*BOUNDARY
N-ZZ1, 3, 3
*BOUNDARY
N-YY2, 2, 2, 0.1
*Node file, TIME POINTS=T1
RF, U, NT
*El File, TIME POINTS=T1
E, S, SDV,HFL
*El print, ELSET=EALL, TIME POINTS=T1
E, S, SDV, EVOL,HFL
*End Step

Full input files can be shared by adding a link to some hosting website like WeTransfer, Dropbox or Google Drive. When the mesh contains just a few elements (like here), it might be also sufficient to use a preformatted text option like you did - then the post is not too long because there’s an additional scroll bar for the code.

I wouldn’t be very surprised if there was a bug leading to unusual behavior when coupled tem-disp step is used for purely mechanical case. I’ve encountered something like that but for opposite case (pure heat transfer in coupled temp-disp analysis): Wrong results if *Coupled temp-disp is used instead of *Heat transfer

Okay perfect. I just updated the mesh. Thanks for the link to the post. I will look in to it.
Now I have to just wait for a newer realease, or play around with different versions to see if one them is clean.

I notice the same behaviour in version 2.18 and 2.19, and for *UNCOUPLED TEMPERATURE-DISPLACEMENT solver

This one works

**
*MATERIAL, NAME=vonMiese1
*Density
7.8E-09
*Elastic
210000, 0.28
*PLASTIC
400.0,0.0
800.0,1.0
*Expansion, Zero=275
1.1E-05,275

*Conductivity
14,275

*Specific heat
440000000,275

*Amplitude, Name=Linear
0, 0, 0.1, 0,1,1

**
*SOLID SECTION, ELSET=EALL, MATERIAL=vonMiese1 
*TIME POINTS,NAME=T1,GENERATE
0.,1.0,0.1

*Initial Conditions, type=TEMPERATURE
NALL, 275.0

** STEP: Load

*Step, nlgeom=YES, inc=100000
*COUPLED TEMPERATURE-DISPLACEMENT

**STATIC
0.001, 1.0, 1e-09, 1.0

*BOUNDARY
N-XX1, 1, 1
*BOUNDARY
N-YY1, 2, 2
*BOUNDARY
N-ZZ1, 3, 3
*BOUNDARY,AMPLITUDE=Linear
N-YY2, 2, 2, 0.1

*Node file, TIME POINTS=T1
RF, U, NT

*El File, TIME POINTS=T1
E, S, SDV,HFL,PEEQ

*El print, ELSET=EALL, TIME POINTS=T1
E, S, SDV, EVOL,HFL

*End Step

@Disla thanks for the input. I investigated a bit with my material parameters.
Looks like using the *Amplitude for the mechanical displacement alone resolves the issue. Which is quite strange.