Hello everyone!
I'm trying to replicate a simulation described in: https://www.sciencedirect.com/science/article/abs/pii/S0020740317323512 .
As a precursor to this, I'm running a number of smaller simulations to get familiarized with LS-Dyna (I'm a simulation engineer, but haven't used LS-Dyna so far). They're all solving a very simple case: a spherical indentor with a diameter of 10 mm (modelled with 1 mm shells, rigid) is being pressed into a 120 mm by 120 mm DC04 sheet with a thickness of 1 mm (modelled with 1x1 shell elements, element type 25 as described in the paper linked above).
From reading the paper it's clear that the authors used Chaboche-Roussilier kinematic hardening. Having read recently that MAT036 with Chaboche-Roussilier was apparently broken for some time, I chose to use MAT133 (Barlat yield 2000). As the isotropic hardening properties of the DC04 material used by the authors were not disclosed, I opted for a simplified approach assuming no in-plane anisotropy (R00=R45=R90,SY00=SY45=SY90) and isotropic hardening according to the BaoSteel curve with a yield stress of 163 MPa.
So far everything worked magnificently, but with the addition of kinematic hardening I am running into a problem that leaves me baffled.
The authors of the paper pointed out that they made use of a single backstress term of the C-R model, with A=80 MPa, C=25500 MPa. Now, I realize that by assuming certain values for the isotropic hardening I will certainly not be able to fully replicate their results - however, I'd still like to make use of the kinematic hardening parameters they used. I did a number of test runs with (A, C) going from (0, 0) to (80, 2500) - all of which were successful. However, setting C to 25500 (or even 5000 or 7500 for that matter) invariably leads to a plasticity convergence failure very early in the simulation, at very benign deformations and strains. I've been banging my head against this problem for the last 2 weeks, playing around with contacts, control parameters, el. formulations and material parameters, but I cannot resolve it.
I would have thought nothing of it if the model failed towards the end of the aforementioned indentation case at 40%, 50% strain, but I find it very odd indeed that it should fail basically as soon as plasticity is activated.
I've attached to this post the last working version of the indentation case (A=80, C=2500 instead of the desired 25500). Perhaps somebody could help me understand:
1) Why the plasticity model fails right away at very low plastic strains or successfully completes (the usual case being either a fundamental error breaking it at the beginning or a difficult model failing at high strains... thus very different from my problem)?
2) Is there perhaps a way to make use of the other C-R terms to "ease" the model into kinematic hardening? Some combination of parameters that would lead to an "equivalent" value of (A, C) closer to (80, 2500) at low strains and closer to (80, 25500) at high strains?
3) Is there any other way to resolve this issue and get robust combined kinematic and isotropic hardening?
Looking forward to all replies!
Best regards,
Leibhussar
Let me help you troubleshoot this plasticity convergence issue with MAT133 and kinematic hardening. This is a complex problem that can arise from several factors:
1. Numerical Integration:
*CONTROL_SOLUTION NLSOLVR = 1 ! Try non-linear solver type 1 *CONTROL_IMPLICIT_SOLUTION NSOLVR = 12 ! Try Pardiso solver ILIMIT = 30 ! Increase iteration limit DCTOL = 1.0E-3 ! Adjust divergence control tolerance
2. Time Step Control:
*CONTROL_IMPLICIT_AUTO IAUTO = 1 ! Automatic time stepping ITEOPT = 25 ! Target iterations per step DTMIN = 1.0E-6 ! Minimum time step size
3. Material Parameter Scaling:
Try scaling down your parameters while maintaining their ratio:
- Instead of C=25500, A=80
- Try C=2550, A=8 (factor of 10)
- Or implement gradual ramping of parameters
4. Alternative Approaches:
a) Consider using MAT125 (Yoshida_Uemori):
*MAT_YOSHIDA_UEMORI C1 = [scaled value] C2 = [scaled value]
b) Or MAT103 (Anisotropic_Plastic):
*MAT_ANISOTROPIC_PLASTIC CP = [adjusted hardening parameter]
5. Debugging Steps:
- Monitor plastic strain increments
- Check for element distortion
- Verify contact behavior
- Review energy balance
Hi @leibhussar let me try to help answer your questions:
- The immediate failure at low strains typically occurs because:
- The high value of C (25500 MPa) creates very steep stress gradients in the yield surface evolution
- This leads to numerical instability in the return mapping algorithm
- At very small plastic strains, the algorithm struggles to converge because:
- The backstress evolution rate is extremely high
- The stress state jumps dramatically between iterations
- The plastic corrector step becomes unstable
- Yes, you can use multiple backstress terms to create a more stable evolution. Here's a suggested approach:
Term 1 (dominant at low strains): - C₁ = 2500 MPa - A₁ = 40 MPa Term 2 (activates at medium strains): - C₂ = 10000 MPa - A₂ = 25 MPa Term 3 (activates at higher strains): - C₃ = 25500 MPa - A₃ = 15 MPa
This creates an "effective" response that:
- Starts softer for numerical stability
- Gradually builds up to your target stiffness
- Maintains similar overall hardening behavior
- Alternative approaches for robust hardening:
- Use implicit time integration
- Increase number of integration points
- Add artificial bulk viscosity
- Use adaptive mesh refinement in high-strain regions