Mechanical Models

Small strain mechanical constitutive models including elasticity, plasticity, viscoelasticity, and damage.

void abaqus2smart_M_light(const double *stress, const double *ddsdde, const int &nstatev, double *statev, const int &ndi, const int &nshr, arma::vec &sigma, arma::mat &Lt, arma::vec &Wm, arma::vec &statev_smart)

Light transfer from Abaqus output arrays to simcoon format (mechanical)

Converts stress, tangent, and state variables from Abaqus format. Use this after calling an external Abaqus UMAT to retrieve results.

Parameters:

  • stress – Abaqus stress array (ntens)

  • ddsdde – Abaqus tangent operator array (ntens×ntens, column-major)

  • nstatev – Number of state variables

  • statev – Abaqus state variables array

  • ndi – Number of direct stress components

  • nshr – Number of shear stress components

  • sigma – [out] simcoon stress vector (6)

  • Lt – [out] simcoon tangent matrix (6×6)

  • Wm – [out] simcoon work quantities vector (4): Wm, Wm_r, Wm_ir, Wm_d

  • statev_smart – [out] simcoon state variables vector

void abaqus2smart_M(const double *stress, const double *ddsdde, const double *stran, const double *dstran, const double *time, const double &dtime, const double &temperature, const double &Dtemperature, const int &nprops, const double *props, const int &nstatev, double *statev, const int &ndi, const int &nshr, const double *drot, arma::vec &sigma, arma::mat &Lt, arma::vec &Etot, arma::vec &DEtot, double &T, double &DT, double &Time, double &DTime, arma::vec &props_smart, arma::vec &Wm, arma::vec &statev_smart, arma::mat &DR, bool &start)

Full transfer from Abaqus input arrays to simcoon format (mechanical)

Converts all UMAT inputs from Abaqus format to simcoon Armadillo format. Use this before calling a simcoon constitutive model from an Abaqus UMAT wrapper.

Parameters:

  • stress – Abaqus stress array

  • ddsdde – Abaqus tangent operator array

  • stran – Abaqus total strain array

  • dstran – Abaqus strain increment array

  • time – Abaqus time array [step_time, total_time]

  • dtime – Time increment

  • temperature – Temperature at start of increment

  • Dtemperature – Temperature increment

  • nprops – Number of material properties

  • props – Abaqus material properties array

  • nstatev – Number of state variables

  • statev – Abaqus state variables array

  • ndi – Number of direct stress components

  • nshr – Number of shear stress components

  • drot – Abaqus rotation increment matrix (3×3, column-major)

  • sigma – [out] simcoon stress vector

  • Lt – [out] simcoon tangent matrix

  • Etot – [out] simcoon total strain vector

  • DEtot – [out] simcoon strain increment vector

  • T – [out] simcoon temperature

  • DT – [out] simcoon temperature increment

  • Time – [out] simcoon total time

  • DTime – [out] simcoon time increment

  • props_smart – [out] simcoon properties vector

  • Wm – [out] simcoon work quantities vector

  • statev_smart – [out] simcoon state variables vector

  • DR – [out] simcoon rotation matrix

  • start – [out] true if this is the first increment

void abaqus2smart_T(const double *stress, const double *ddsdde, const double *ddsddt, const double *drplde, const double &drpldt, const double *stran, const double *dstran, const double *time, const double &dtime, const double &temperature, const double &Dtemperature, const int &nprops, const double *props, const int &nstatev, double *statev, const int &ndi, const int &nshr, const double *drot, arma::vec &sigma, arma::mat &dSdE, arma::mat &dSdT, arma::mat &drpldE, arma::mat &drpldT, arma::vec &Etot, arma::vec &DEtot, double &T, double &DT, double &Time, double &DTime, arma::vec &props_smart, arma::vec &Wm, arma::vec &Wt, arma::vec &statev_smart, arma::mat &DR, bool &start)

Full transfer from Abaqus input arrays to simcoon format (thermomechanical)

Converts all thermomechanical UMAT inputs from Abaqus format. Includes thermal tangent operators and heat generation terms.

Parameters:

  • stress – Abaqus stress array

  • ddsdde – Abaqus mechanical tangent array

  • ddsddt – Abaqus stress-temperature tangent array

  • drplde – Abaqus heat generation-strain tangent array

  • drpldt – Abaqus heat generation-temperature tangent

  • stran – Abaqus total strain array

  • dstran – Abaqus strain increment array

  • time – Abaqus time array

  • dtime – Time increment

  • temperature – Temperature at start of increment

  • Dtemperature – Temperature increment

  • nprops – Number of material properties

  • props – Abaqus material properties array

  • nstatev – Number of state variables

  • statev – Abaqus state variables array

  • ndi – Number of direct stress components

  • nshr – Number of shear stress components

  • drot – Abaqus rotation increment matrix

  • sigma – [out] simcoon stress vector

  • dSdE – [out] simcoon mechanical tangent matrix

  • dSdT – [out] simcoon stress-temperature tangent matrix

  • drpldE – [out] simcoon heat generation-strain tangent matrix

  • drpldT – [out] simcoon heat generation-temperature tangent matrix

  • Etot – [out] simcoon total strain vector

  • DEtot – [out] simcoon strain increment vector

  • T – [out] simcoon temperature

  • DT – [out] simcoon temperature increment

  • Time – [out] simcoon total time

  • DTime – [out] simcoon time increment

  • props_smart – [out] simcoon properties vector

  • Wm – [out] simcoon mechanical work quantities vector

  • Wt – [out] simcoon thermal work quantities vector

  • statev_smart – [out] simcoon state variables vector

  • DR – [out] simcoon rotation matrix

  • start – [out] true if this is the first increment

void smart2abaqus_M(double *stress, double *ddsdde, double *statev, const int &ndi, const int &nshr, const arma::vec &sigma, const arma::vec &statev_smart, const arma::vec &Wm, const arma::mat &Lt)

Transfer simcoon output to Abaqus format (mechanical, simple)

Converts stress, tangent, and state variables to Abaqus format. Use this to return results from a simcoon model to an Abaqus UMAT.

Parameters:

  • stress – [out] Abaqus stress array

  • ddsdde – [out] Abaqus tangent operator array

  • statev – [out] Abaqus state variables array

  • ndi – Number of direct stress components

  • nshr – Number of shear stress components

  • sigma – simcoon stress vector

  • statev_smart – simcoon state variables vector

  • Wm – simcoon work quantities vector

  • Lt – simcoon tangent matrix

void smart2abaqus_M_full(double *stress, double *ddsdde, double *stran, double *dstran, double *time, double &dtime, double &temperature, double &Dtemperature, int &nprops, double *props, int &nstatev, double *statev, const int &ndi, const int &nshr, double *drot, const arma::vec &sigma, const arma::mat &Lt, const arma::vec &Etot, const arma::vec &DEtot, const double &T, const double &DT, const double &Time, const double &DTime, const arma::vec &props_smart, const arma::vec &Wm, const arma::vec &statev_smart, const arma::mat &DR, bool &start)

Full transfer from simcoon to Abaqus format (mechanical)

Converts all simcoon outputs to Abaqus UMAT format. Use this for complete state transfer in external UMAT wrappers.

Parameters:

  • stress – [out] Abaqus stress array

  • ddsdde – [out] Abaqus tangent operator array

  • stran – [out] Abaqus total strain array

  • dstran – [out] Abaqus strain increment array

  • time – [out] Abaqus time array

  • dtime – [out] Time increment

  • temperature – [out] Temperature

  • Dtemperature – [out] Temperature increment

  • nprops – Number of material properties

  • props – [out] Abaqus material properties array

  • nstatev – Number of state variables

  • statev – [out] Abaqus state variables array

  • ndi – Number of direct stress components

  • nshr – Number of shear stress components

  • drot – [out] Abaqus rotation increment matrix

  • sigma – simcoon stress vector

  • Lt – simcoon tangent matrix

  • Etot – simcoon total strain vector

  • DEtot – simcoon strain increment vector

  • T – simcoon temperature

  • DT – simcoon temperature increment

  • Time – simcoon total time

  • DTime – simcoon time increment

  • props_smart – simcoon properties vector

  • Wm – simcoon work quantities vector

  • statev_smart – simcoon state variables vector

  • DR – simcoon rotation matrix

  • start – simcoon start flag

void smart2abaqus_T(double *stress, double *ddsdde, double *ddsddt, double *drplde, double &drpldt, double &rpl, double *statev, const int &ndi, const int &nshr, const arma::vec &sigma, const arma::vec &statev_smart, const double &r, const arma::vec &Wm, const arma::vec &Wt, const arma::mat &dSdE, const arma::mat &dSdT, const arma::mat &drpldE, const arma::mat &drpldT)

Transfer simcoon output to Abaqus format (thermomechanical)

Converts thermomechanical outputs to Abaqus UMAT format. Includes thermal tangent operators and heat generation.

Parameters:

  • stress – [out] Abaqus stress array

  • ddsdde – [out] Abaqus mechanical tangent array

  • ddsddt – [out] Abaqus stress-temperature tangent array

  • drplde – [out] Abaqus heat generation-strain tangent array

  • drpldt – [out] Abaqus heat generation-temperature tangent

  • rpl – [out] Abaqus heat generation rate

  • statev – [out] Abaqus state variables array

  • ndi – Number of direct stress components

  • nshr – Number of shear stress components

  • sigma – simcoon stress vector

  • statev_smart – simcoon state variables vector

  • r – simcoon heat generation rate

  • Wm – simcoon mechanical work quantities vector

  • Wt – simcoon thermal work quantities vector

  • dSdE – simcoon mechanical tangent matrix

  • dSdT – simcoon stress-temperature tangent matrix

  • drpldE – simcoon heat generation-strain tangent matrix

  • drpldT – simcoon heat generation-temperature tangent matrix

void ansys2smart_M(const double *stress, const double *dstran, const double &sedEl, const double &sedPl, const double &epseq, const double *statev, const double *props, const double &Time, const double &DTime, const double &temperature, const double &Dtemperature, const int &ncomp, const int &nprops, const int &nstatev, arma::vec &sigma, arma::vec &DEtot, arma::vec &Wm, arma::vec &statev_smart, arma::vec &props_smart, double &T, double &DT)

Transfer Ansys USERMAT input arrays to simcoon format (mechanical)

Converts all USERMAT inputs from Ansys format to simcoon format. Handles Ansys Voigt notation: (11,22,33,12,23,13) by swapping components 4 and 5.

Parameters:

  • stress – Ansys stress array (ncomp)

  • dstran – Ansys strain increment array (ncomp)

  • sedEl – Specific elastic strain energy

  • sedPl – Specific plastic strain energy

  • epseq – Equivalent plastic strain

  • statev – Ansys state variables array

  • props – Ansys material properties array

  • Time – Total time

  • DTime – Time increment

  • temperature – Temperature

  • Dtemperature – Temperature increment

  • ncomp – Number of stress/strain components

  • nprops – Number of material properties

  • nstatev – Number of state variables

  • sigma – [out] simcoon stress vector

  • DEtot – [out] simcoon strain increment vector

  • Wm – [out] simcoon work quantities vector

  • statev_smart – [out] simcoon state variables vector

  • props_smart – [out] simcoon properties vector

  • T – [out] simcoon temperature

  • DT – [out] simcoon temperature increment

void smart2ansys_M(double *stress, double *ddsdde, double &sedEl, double &sedPl, double *statev, const int &ncomp, const arma::vec &sigma, const arma::mat &Lt, const arma::vec &Wm, const arma::vec &statev_smart)

Transfer simcoon output to Ansys USERMAT format (mechanical)

Converts stress, tangent, and state variables to Ansys format. Handles Ansys Voigt notation by swapping components 4 and 5.

Parameters:

  • stress – [out] Ansys stress array

  • ddsdde – [out] Ansys tangent operator array

  • sedEl – [out] Specific elastic strain energy

  • sedPl – [out] Specific plastic strain energy

  • statev – [out] Ansys state variables array

  • ncomp – Number of stress/strain components

  • sigma – simcoon stress vector

  • Lt – simcoon tangent matrix

  • Wm – simcoon work quantities vector

  • statev_smart – simcoon state variables vector

void umat_damage_LLD_0(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, arma::vec &sigma_in, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt)

Ladevèze-Le Dantec (LLD) anisotropic damage model with coupled plasticity for composite materials.

This function implements the Ladevèze-Le Dantec damage model specifically designed for unidirectional fiber-reinforced composite materials. The model features:

  • Transversely isotropic elastic behavior

  • Anisotropic damage in matrix-dominated modes (transverse tension, in-plane shear)

  • Coupled transverse-shear plasticity

  • Thermodynamically consistent damage evolution based on energy release rate

  • Separate damage variables for different failure modes

  • No fiber failure (matrix-dominated damage only)

Material Symmetry:

The model assumes transverse isotropy with the fiber direction as the axis of symmetry:

  • Direction 1: Fiber direction (longitudinal, no damage)

  • Directions 2-3: Transverse plane (isotropic, subject to damage)

Damage Variables:

Two independent scalar damage variables characterize the material state:

  • \( d_{22} \in [0,1] \): Transverse damage (matrix cracking perpendicular to fibers)

  • \( d_{12} \in [0,1] \): Shear damage (matrix/interface damage in fiber-matrix plane)

where \( d = 0 \) is undamaged and \( d = 1 \) is fully damaged.

Effective Stress Concept:

The effective stress acting on the undamaged material configuration is:

\[\tilde{\boldsymbol{\sigma}} = \mathbf{M}^{-1} : \boldsymbol{\sigma} \]
where \( \mathbf{M} \) is the damage effect tensor:
\[\mathbf{M} = \text{diag}(1, 1-d_{22}, 1-d_{22}, 1-d_{12}, 1-d_{12}, 1) \]

Damaged Elastic Stiffness:

The elastic stiffness degrades with damage:

\[\mathbf{L}(d_{22}, d_{12}) = \mathbf{L}_0 : \mathbf{M} \]
where \( \mathbf{L}_0 \) is the undamaged transversely isotropic stiffness.

Specific moduli degradation:

  • \( E_2 = E_{2,0} (1 - d_{22}) \) (transverse Young’s modulus)

  • \( E_3 = E_{3,0} (1 - d_{22}) \) (out-of-plane Young’s modulus)

  • \( G_{12} = G_{12,0} (1 - d_{12}) \) (in-plane shear modulus)

  • \( G_{13} = G_{13,0} (1 - d_{12}) \) (out-of-plane shear modulus)

  • \( E_1 \) (fiber direction) remains constant (no fiber damage)

Thermodynamic Forces:

The energy release rates driving damage evolution are:

\[Y_{22} = \frac{1}{2} \boldsymbol{\varepsilon}^e : \frac{\partial \mathbf{L}}{\partial d_{22}} : \boldsymbol{\varepsilon}^e \]
\[Y_{12} = \frac{1}{2} \boldsymbol{\varepsilon}^e : \frac{\partial \mathbf{L}}{\partial d_{12}} : \boldsymbol{\varepsilon}^e \]

Damage Evolution Laws:

Transverse Damage (d₂₂):

\[\begin{split}\dot{d}_{22} = \begin{cases} 0 & \text{if } Y_{22} < Y_{22,0} \\ \left( \frac{Y_{22} - Y_{22,c}}{Y_{22,u} - Y_{22,c}} \right)^b & \text{if } Y_{22,0} \leq Y_{22} < Y_{22,u} \\ \infty & \text{if } Y_{22} \geq Y_{22,u} \end{cases} \end{split}\]
where:

  • \( Y_{22,0} \) is the damage initiation threshold

  • \( Y_{22,c} \) is the characteristic energy release rate

  • \( Y_{22,u} \) is the ultimate energy release rate (failure)

  • \( b \) is the damage evolution exponent

Shear Damage (d₁₂):

\[\begin{split}\dot{d}_{12} = \begin{cases} 0 & \text{if } Y_{12} < Y_{12,0} \\ \frac{Y_{12} - Y_{12,c}}{Y_{12,c}} & \text{if } Y_{12} \geq Y_{12,0} \end{cases} \end{split}\]
where:

  • \( Y_{12,0} \) is the shear damage initiation threshold

  • \( Y_{12,c} \) is the characteristic shear energy release rate

Coupled Transverse-Shear Plasticity:

A Hill-type yield criterion couples transverse and shear stresses:

\[\Phi = \sqrt{\left( \frac{\sigma_{22}}{A_{ts}} \right)^2 + \sigma_{12}^2} - \sigma_{ts,0} - \alpha_{ts} p_{ts} - \beta_{ts} p_{ts}^2 \leq 0 \]
where:

  • \( A_{ts} \) is the transverse-shear coupling parameter

  • \( \sigma_{ts,0} \) is the initial yield stress

  • \( \alpha_{ts}, \beta_{ts} \) are hardening parameters

  • \( p_{ts} \) is the accumulated plastic strain

Plastic Flow Rule:

Associative flow in the transverse-shear plane:

\[\dot{\boldsymbol{\varepsilon}}^p = \dot{p}_{ts} \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} \]

Material Parameters (props):

Index

Symbol

Description

Units

Typical Range (CFRP)

props[0]

axis

Fiber orientation axis (1, 2, or 3)

-

1

props[1]

\( E_L \)

Longitudinal Young’s modulus

Stress

100-200 GPa

props[2]

\( E_T \)

Transverse Young’s modulus

Stress

5-15 GPa

props[3]

\( \nu_{TL} \)

Major Poisson’s ratio

-

0.25-0.35

props[4]

\( \nu_{TT} \)

Transverse Poisson’s ratio

-

0.35-0.45

props[5]

\( G_{LT} \)

In-plane shear modulus

Stress

3-8 GPa

props[6]

\( \alpha_L \)

Longitudinal CTE

1/Temperature

-0.5 to 0 e-6 /K

props[7]

\( \alpha_T \)

Transverse CTE

1/Temperature

20-40 e-6 /K

props[8]

\( Y_{12,0} \)

Shear damage initiation threshold

Energy/Volume

0.05-0.2 MPa

props[9]

\( Y_{12,c} \)

Characteristic shear energy release

Energy/Volume

0.1-0.5 MPa

props[10]

\( Y_{22,0} \)

Transverse damage initiation

Energy/Volume

0.1-0.3 MPa

props[11]

\( Y_{22,c} \)

Characteristic transverse energy

Energy/Volume

0.2-0.6 MPa

props[12]

\( Y_{22,u} \)

Ultimate transverse energy

Energy/Volume

0.5-2.0 MPa

props[13]

\( b \)

Transverse damage exponent

-

1-5

props[14]

\( A_{ts} \)

Transverse-shear coupling

-

1-3

props[15]

\( \sigma_{ts,0} \)

Initial yield stress

Stress

30-80 MPa

props[16]

\( \alpha_{ts} \)

Linear hardening parameter

Stress

0-500 MPa

props[17]

\( \beta_{ts} \)

Quadratic hardening parameter

Stress

0-5000 MPa

State Variables (statev):

Total state variables required: \( n_{statev} = 10 \)

Index

Symbol

Description

Units

statev[0]

\( T_{init} \)

Initial/reference temperature

Temperature

statev[1]

\( d_{22} \)

Transverse damage variable

-

statev[2]

\( d_{12} \)

Shear damage variable

-

statev[3]

\( p_{ts} \)

Accumulated plastic strain

Strain

statev[4]

\( \varepsilon^p_{11} \)

Plastic strain component 11

Strain

statev[5]

\( \varepsilon^p_{22} \)

Plastic strain component 22

Strain

statev[6]

\( \varepsilon^p_{33} \)

Plastic strain component 33

Strain

statev[7]

\( \varepsilon^p_{12} \)

Plastic strain component 12 (engineering)

Strain

statev[8]

\( \varepsilon^p_{13} \)

Plastic strain component 13 (engineering)

Strain

statev[9]

\( \varepsilon^p_{23} \)

Plastic strain component 23 (engineering)

Strain 

// Example usage: Carbon/Epoxy unidirectional composite (T300/914)
vec props(18);
props(0) = 1;           // axis = 1 (fibers along direction 1)
props(1) = 138000;      // EL = 138 GPa (fiber-dominated)
props(2) = 11000;       // ET = 11 GPa (matrix-dominated)
props(3) = 0.28;        // nuTL = 0.28
props(4) = 0.40;        // nuTT = 0.40
props(5) = 5500;        // GLT = 5.5 GPa
props(6) = -0.3e-6;     // alphaL = -0.3e-6 /K (negative for carbon fibers)
props(7) = 30e-6;       // alphaT = 30e-6 /K (resin dominated)
props(8) = 0.10;        // Y_12_0 = 0.10 MPa (shear damage threshold)
props(9) = 0.30;        // Y_12_c = 0.30 MPa
props(10) = 0.15;       // Y_22_0 = 0.15 MPa (transverse damage threshold)
props(11) = 0.40;       // Y_22_c = 0.40 MPa
props(12) = 1.20;       // Y_22_u = 1.20 MPa (ultimate failure)
props(13) = 2.5;        // b = 2.5
props(14) = 2.0;        // A_ts = 2.0
props(15) = 50;         // sigma_ts_0 = 50 MPa
props(16) = 300;        // alpha_ts = 300 MPa
props(17) = 1000;       // beta_ts = 1000 MPa

vec statev = zeros(10);
statev(0) = 20.0;  // Reference temperature 20°C

vec Etot = {0.0, 0.005, 0.0, 0.0, 0.0, 0.0};  // 0.5% transverse strain (induces damage)
vec DEtot = {0.0, 0.0001, 0.0, 0.0, 0.0, 0.0};
vec sigma = zeros(6);
mat Lt = zeros(6,6);
mat L = zeros(6,6);
vec sigma_in = zeros(6);
mat DR = eye(3,3);

umat_damage_LLD_0(Etot, DEtot, sigma, Lt, L, sigma_in, DR,
                  18, props, 10, statev, 20.0, 0.0, 0.0, 1.0,
                  Wm, Wm_r, Wm_ir, Wm_d, 3, 3, false, 0, tnew_dt);

// Check damage state
double d22 = statev(1);  // Transverse damage
double d12 = statev(2);  // Shear damage
cout << "Transverse damage: " << d22 << ", Shear damage: " << d12 << endl;

See also

L_isotrans() for transversely isotropic stiffness construction

See also

Lagrange_exp() for penalty function (Damage Functions bounds enforcement)

See also

rotate_strain() for strain tensor rotation

See also

rotate_stress() for Stress Functions tensor rotation

References:

  • Ladevèze, P., & Le Dantec, E. (1992). “Damage modelling of the elementary ply for laminated composites.” Composites Science and Technology, 43(3), 257-267.

  • Allix, O., & Ladevèze, P. (1992). “Interlaminar interface modelling for the prediction of delamination.” Composite Structures, 22(4), 235-242.

  • Ladevèze, P. (1992). “A damage computational method for composite structures.” Computers & Structures, 44(1-2), 79-87.

  • Matzenmiller, A., Lubliner, J., & Taylor, R. L. (1995). “A constitutive model for anisotropic damage in fiber-composites.” Mechanics of Materials, 20(2), 125-152.

  • Pinho, S. T., et al. (2012). “Material and structural response of polymer-matrix fibre-reinforced composites.” Journal of Composite Materials, 46(19-20), 2313-2341.

Note

The LLD model is specifically designed for unidirectional fiber composites

Note

Damage is irreversible and monotonically increasing

Note

The model does NOT account for fiber failure (compression/tension in direction 1)

Note

Suitable for matrix-dominated failure modes: transverse cracking, delamination

Note

Parameter identification requires multiple test configurations:

Note

- Transverse tension for d₂₂ parameters

Note

- In-plane shear for d₁₂ parameters

Note

- Off-axis tests for plasticity coupling

Note

Convergence requires small load steps once damage initiates

Note

Material axes must be properly oriented relative to global coordinates

Parameters:

  • Etot – Total strain tensor at beginning of increment (Voigt notation: 6×1 vector)

  • DEtot – Strain increment tensor (Voigt notation: 6×1 vector)

  • sigma – Stress tensor (Voigt notation: 6×1 vector) [output]

  • Lt – Consistent tangent modulus (6×6 matrix) [output]

  • L – Damaged elastic stiffness tensor (6×6 matrix) [output]

  • sigma_in – Internal stress contribution for explicit solvers (6×1 vector) [output]

  • DR – Rotation increment matrix (3×3) for objective integration

  • nprops – Number of material properties

  • props – Material properties vector (see table above)

  • nstatev – Number of state variables

  • statev – State variables vector (see table above) [input/output]

  • T – Temperature at beginning of increment

  • DT – Temperature increment

  • Time – Time at beginning of increment

  • DTime – Time increment

  • Wm – Total mechanical work [output]

  • Wm_r – Recoverable (elastic) work [output]

  • Wm_ir – Irrecoverable work stored in damage [output]

  • Wm_d – Dissipated work (damage + plasticity) [output]

  • ndi – Number of direct stress components (typically 3)

  • nshr – Number of shear stress components (typically 3)

  • start – Flag indicating first increment (true) or continuation (false)

  • solver_type – Solver type: 0=implicit, 1=explicit, 2=dynamic implicit

  • tnew_dt – Suggested new time step size for adaptive time stepping [output]

void umat_damage_weibull(const arma::vec&, const arma::vec&, arma::vec&, arma::mat&, arma::mat&, arma::vec&, const arma::mat&, const int&, const arma::vec&, const int&, arma::vec&, const double&, const double&, const double&, const double&, double&, double&, double&, double&, const int&, const int&, const bool&, const int&, double&)
void umat_elasticity_iso(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, arma::vec &sigma_in, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt)

Linear elastic isotropic constitutive model with thermal expansion.

This function implements the fundamental linear elastic isotropic material model based on Hooke’s law, suitable for small strain analysis of materials exhibiting reversible elastic deformation.

Constitutive Equations:

The stress-strain relationship follows the generalized Hooke’s law:

\[\boldsymbol{\sigma} = \mathbf{L} : (\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^{th}) \]

where:

  • \( \boldsymbol{\sigma} \) is the stress tensor (Voigt notation)

  • \( \mathbf{L} \) is the fourth-order elastic stiffness tensor

  • \( \boldsymbol{\varepsilon} \) is the total strain tensor

  • \( \boldsymbol{\varepsilon}^{th} = \alpha (T - T_0) \mathbf{I} \) is the thermal strain

  • \( \alpha \) is the coefficient of thermal expansion

  • \( T \) is the current temperature

  • \( T_0 \) is the reference temperature

Elastic Stiffness Tensor:

For an isotropic material, the stiffness tensor is constructed from two independent constants:

\[\mathbf{L} = 3K \mathbf{I}_{vol} + 2\mu \mathbf{I}_{dev} \]

where:

  • \( K = \frac{E}{3(1-2\nu)} \) is the bulk modulus

  • \( \mu = \frac{E}{2(1+\nu)} \) is the shear modulus

  • \( \mathbf{I}_{vol} \) is the volumetric projection tensor

  • \( \mathbf{I}_{dev} \) is the deviatoric projection tensor

In component form (Voigt notation):

\[\begin{split}\mathbf{L} = \frac{E}{(1+\nu)(1-2\nu)} \begin{bmatrix} 1-\nu & \nu & \nu & 0 & 0 & 0 \\ \nu & 1-\nu & \nu & 0 & 0 & 0 \\ \nu & \nu & 1-\nu & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1-2\nu}{2} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1-2\nu}{2} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1-2\nu}{2} \end{bmatrix} \end{split}\]

Stress Update:

For an increment of strain \( \Delta \boldsymbol{\varepsilon} \):

\[\boldsymbol{\sigma}_{n+1} = \boldsymbol{\sigma}_n + \mathbf{L} : \left( \Delta \boldsymbol{\varepsilon} - \alpha \Delta T \mathbf{I} \right) \]

Tangent Modulus:

For linear elasticity, the tangent modulus is constant and equal to the elastic stiffness:

\[\mathbf{L}_t = \mathbf{L} \]

This ensures optimal convergence in Newton-Raphson iterations.

Material Parameters (props):

Index

Symbol

Description

Units

Typical Range

props[0]

\( E \)

Young’s modulus

Stress

1-500 GPa

props[1]

\( \nu \)

Poisson’s ratio

-

0.0-0.5

props[2]

\( \alpha \)

Coefficient of thermal expansion

1/Temperature

1e-6 to 1e-4 /K

Constraint on Poisson’s ratio:

  • For 3D: \( -1 < \nu < 0.5 \)

  • \( \nu = 0.5 \) implies incompressibility (infinite bulk modulus)

  • \( \nu = 0 \) implies no lateral strain under uniaxial loading

State Variables (statev):

Index

Symbol

Description

Units

statev[0]

\( T_{init} \)

Initial/reference temperature

Temperature

Total state variables required: \( n_{statev} = 1 \)

// Example usage: Steel
vec props = {210000, 0.3, 1.2e-5};  // E=210 GPa, nu=0.3, alpha=12e-6 /K
vec statev = zeros(1);
statev(0) = 20.0;  // Reference temperature 20°C

vec Etot = {0.001, 0.0, 0.0, 0.0, 0.0, 0.0};  // 0.1% strain in direction 1
vec DEtot = {0.0001, 0.0, 0.0, 0.0, 0.0, 0.0};
vec sigma = zeros(6);
mat Lt = zeros(6,6);
mat L = zeros(6,6);
vec sigma_in = zeros(6);
mat DR = eye(3,3);

umat_elasticity_iso(Etot, DEtot, sigma, Lt, L, sigma_in, DR,
                    3, props, 1, statev, 25.0, 5.0, 0.0, 1.0,
                    Wm, Wm_r, Wm_ir, Wm_d, 3, 3, false, 0, tnew_dt);

// Expected stress: sigma(0) ≈ E*Etot(0)*(1-nu)/((1+nu)(1-2nu)) - E*alpha*DT/(1-2nu)

See also

L_iso() for elastic stiffness tensor construction

See also

Ivol() for volumetric projection tensor

See also

Idev() for deviatoric projection tensor

See also

Ith() for thermal expansion vector

References:

  • Timoshenko, S. P., & Goodier, J. N. (1970). Theory of Elasticity (3rd ed.). McGraw-Hill.

  • Bower, A. F. (2009). Applied Mechanics of Solids. CRC Press.

  • Gurtin, M. E. (1981). An Introduction to Continuum Mechanics. Academic Press.

Note

This model is restricted to small strains (< 1%)

Note

For large strains, use hyperelastic models (Neo-Hookean, Mooney-Rivlin)

Note

The model is path-independent (no history dependence)

Note

No energy dissipation occurs (Wm_ir = Wm_d = 0, Wm = Wm_r)

Note

For plane stress/strain, set ndi=2 or ndi=1 accordingly

Parameters:

  • Etot – Total strain tensor at beginning of increment (Voigt notation: 6×1 vector)

  • DEtot – Strain increment tensor (Voigt notation: 6×1 vector)

  • sigma – Stress tensor (Voigt notation: 6×1 vector) [output]

  • Lt – Tangent modulus \( \mathbf{L}_t = \mathbf{L} \) (6×6 matrix) [output]

  • L – Elastic stiffness tensor (6×6 matrix) [output]

  • sigma_in – Internal stress contribution for explicit solvers (6×1 vector) [output]

  • DR – Rotation increment matrix (3×3) for objective integration

  • nprops – Number of material properties

  • props – Material properties vector (see table above)

  • nstatev – Number of state variables

  • statev – State variables vector (see table above) [input/output]

  • T – Temperature at beginning of increment

  • DT – Temperature increment

  • Time – Time at beginning of increment

  • DTime – Time increment

  • Wm – Total mechanical work [output]

  • Wm_r – Recoverable (elastic) work [output]

  • Wm_ir – Irrecoverable work (0 for elastic model) [output]

  • Wm_d – Dissipated work (0 for elastic model) [output]

  • ndi – Number of direct stress components (typically 3)

  • nshr – Number of shear stress components (typically 3)

  • start – Flag indicating first increment (true) or continuation (false)

  • solver_type – Solver type: 0=implicit, 1=explicit, 2=dynamic implicit

  • tnew_dt – Suggested new time step size for adaptive time stepping [output]

void umat_elasticity_ortho(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, arma::vec &sigma_in, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt)

Linear elastic orthotropic constitutive model with thermal expansion.

This function implements a linear elastic orthotropic material model suitable for materials with three perpendicular planes of symmetry, such as:

  • Wood and timber

  • Rolled metals

  • Unidirectional composites (approximation)

  • Crystals with orthorhombic symmetry

Constitutive Equations:

The stress-strain relationship in the material principal axes:

\[\begin{split}\begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ \gamma_{12} \\ \gamma_{13} \\ \gamma_{23} \end{bmatrix} = \begin{bmatrix} 1/E_1 & -\nu_{12}/E_1 & -\nu_{13}/E_1 & 0 & 0 & 0 \\ -\nu_{12}/E_1 & 1/E_2 & -\nu_{23}/E_2 & 0 & 0 & 0 \\ -\nu_{13}/E_1 & -\nu_{23}/E_2 & 1/E_3 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1/G_{12} & 0 & 0 \\ 0 & 0 & 0 & 0 & 1/G_{13} & 0 \\ 0 & 0 & 0 & 0 & 0 & 1/G_{23} \end{bmatrix} \begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \tau_{12} \\ \tau_{13} \\ \tau_{23} \end{bmatrix} + \begin{bmatrix} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ 0 \\ 0 \\ 0 \end{bmatrix} \Delta T \end{split}\]

Symmetry Requirements:

The compliance matrix must be symmetric, which imposes:

\[\frac{\nu_{ij}}{E_i} = \frac{\nu_{ji}}{E_j} \]

Material Parameters (props):

Index

Symbol

Description

Units

props[0]

\( E_1 \)

Young’s modulus in direction 1

Stress

props[1]

\( E_2 \)

Young’s modulus in direction 2

Stress

props[2]

\( E_3 \)

Young’s modulus in direction 3

Stress

props[3]

\( \nu_{12} \)

Poisson’s ratio (strain in 2 due to stress in 1)

-

props[4]

\( \nu_{13} \)

Poisson’s ratio (strain in 3 due to stress in 1)

-

props[5]

\( \nu_{23} \)

Poisson’s ratio (strain in 3 due to stress in 2)

-

props[6]

\( G_{12} \)

Shear modulus in 1-2 plane

Stress

props[7]

\( G_{13} \)

Shear modulus in 1-3 plane

Stress

props[8]

\( G_{23} \)

Shear modulus in 2-3 plane

Stress

props[9]

\( \alpha_1 \)

CTE in direction 1

1/Temperature

props[10]

\( \alpha_2 \)

CTE in direction 2

1/Temperature

props[11]

\( \alpha_3 \)

CTE in direction 3

1/Temperature

State Variables (statev):

Index

Symbol

Description

Units

statev[0]

\( T_{init} \)

Initial/reference temperature

Temperature

See also

L_ortho() for orthotropic stiffness tensor construction

Note

Material axes must be aligned with global axes. For rotated materials, use local orientation definitions in the FE software.

Note

The stiffness matrix must be positive definite for physical stability

Parameters:

  • Etot – Total strain tensor at beginning of increment (Voigt notation: 6×1)

  • DEtot – Strain increment tensor (Voigt notation: 6×1)

  • sigma – Stress tensor [output] (Voigt notation: 6×1)

  • Lt – Tangent modulus (= L for linear elasticity) [output] (6×6)

  • L – Elastic stiffness tensor [output] (6×6)

  • sigma_in – Internal stress for explicit solvers [output] (6×1)

  • DR – Rotation increment matrix (3×3)

  • nprops – Number of material properties (12)

  • props – Material properties vector

  • nstatev – Number of state variables (1)

  • statev – State variables vector [input/output]

  • T – Temperature at beginning of increment

  • DT – Temperature increment

  • Time – Time at beginning of increment

  • DTime – Time increment

  • Wm – Total mechanical work [output]

  • Wm_r – Recoverable (elastic) work [output]

  • Wm_ir – Irrecoverable work (0 for elastic) [output]

  • Wm_d – Dissipated work (0 for elastic) [output]

  • ndi – Number of direct stress components

  • nshr – Number of shear stress components

  • start – Flag indicating first increment

  • solver_type – Solver type (0=implicit, 1=explicit)

  • tnew_dt – Suggested new time step [output]

void umat_elasticity_trans_iso(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, arma::vec &sigma_in, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt)

Linear elastic transversely isotropic constitutive model with thermal expansion.

This function implements a linear elastic transversely isotropic material model suitable for materials with one preferred direction, such as:

  • Unidirectional fiber composites

  • Hexagonal crystals

  • Biological tissues (muscles, tendons)

  • Extruded or drawn materials

Material Symmetry:

The material has rotational symmetry about axis 1 (fiber direction). Properties are isotropic in the 2-3 plane (transverse plane).

Constitutive Equations:

The compliance matrix in the material principal axes:

\[\begin{split}\begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ \gamma_{12} \\ \gamma_{13} \\ \gamma_{23} \end{bmatrix} = \begin{bmatrix} 1/E_L & -\nu_{LT}/E_L & -\nu_{LT}/E_L & 0 & 0 & 0 \\ -\nu_{LT}/E_L & 1/E_T & -\nu_{TT}/E_T & 0 & 0 & 0 \\ -\nu_{LT}/E_L & -\nu_{TT}/E_T & 1/E_T & 0 & 0 & 0 \\ 0 & 0 & 0 & 1/G_{LT} & 0 & 0 \\ 0 & 0 & 0 & 0 & 1/G_{LT} & 0 \\ 0 & 0 & 0 & 0 & 0 & 2(1+\nu_{TT})/E_T \end{bmatrix} \begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \tau_{12} \\ \tau_{13} \\ \tau_{23} \end{bmatrix} \end{split}\]

Note: The transverse shear modulus is not independent:

\[G_{TT} = \frac{E_T}{2(1 + \nu_{TT})} \]

Material Parameters (props):

Index

Symbol

Description

Units

props[0]

\( E_L \)

Longitudinal Young’s modulus (fiber direction)

Stress

props[1]

\( E_T \)

Transverse Young’s modulus

Stress

props[2]

\( \nu_{LT} \)

Poisson’s ratio (transverse strain due to longitudinal stress)

-

props[3]

\( \nu_{TT} \)

Transverse Poisson’s ratio (in the isotropic plane)

-

props[4]

\( G_{LT} \)

Longitudinal shear modulus

Stress

props[5]

\( \alpha_L \)

Longitudinal CTE

1/Temperature

props[6]

\( \alpha_T \)

Transverse CTE

1/Temperature

Typical Values (UD Carbon/Epoxy):

  • \( E_L \) ≈ 140 GPa

  • \( E_T \) ≈ 10 GPa

  • \( \nu_{LT} \) ≈ 0.3

  • \( \nu_{TT} \) ≈ 0.4

  • \( G_{LT} \) ≈ 5 GPa

State Variables (statev):

Index

Symbol

Description

Units

statev[0]

\( T_{init} \)

Initial/reference temperature

Temperature

See also

L_isotrans() for transversely isotropic stiffness tensor construction

Note

The fiber direction (axis 1) must be aligned with the global x-axis. For rotated fibers, use local orientation definitions.

Note

The 2-3 plane is the plane of isotropy (transverse plane)

Parameters:

  • Etot – Total strain tensor at beginning of increment (Voigt notation: 6×1)

  • DEtot – Strain increment tensor (Voigt notation: 6×1)

  • sigma – Stress tensor [output] (Voigt notation: 6×1)

  • Lt – Tangent modulus (= L for linear elasticity) [output] (6×6)

  • L – Elastic stiffness tensor [output] (6×6)

  • sigma_in – Internal stress for explicit solvers [output] (6×1)

  • DR – Rotation increment matrix (3×3)

  • nprops – Number of material properties (7)

  • props – Material properties vector

  • nstatev – Number of state variables (1)

  • statev – State variables vector [input/output]

  • T – Temperature at beginning of increment

  • DT – Temperature increment

  • Time – Time at beginning of increment

  • DTime – Time increment

  • Wm – Total mechanical work [output]

  • Wm_r – Recoverable (elastic) work [output]

  • Wm_ir – Irrecoverable work (0 for elastic) [output]

  • Wm_d – Dissipated work (0 for elastic) [output]

  • ndi – Number of direct stress components

  • nshr – Number of shear stress components

  • start – Flag indicating first increment

  • solver_type – Solver type (0=implicit, 1=explicit)

  • tnew_dt – Suggested new time step [output]

void umat_ani_chaboche_CCP(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, arma::vec &sigma_in, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt)

Elastic-plastic model with general anisotropic yield criterion and Chaboche kinematic hardening.

This function implements an advanced anisotropic plasticity model combining a general quadratic anisotropic yield criterion with Chaboche nonlinear kinematic hardening. This model provides maximum flexibility for capturing complex anisotropic cyclic plasticity behavior.

Key Features:

  • General quadratic anisotropic yield criterion (21 independent parameters)

  • Multiple Armstrong-Frederick backstresses for nonlinear kinematic hardening

  • Optional Voce-type isotropic hardening

  • Convex Cutting Plane algorithm for return mapping

  • Consistent tangent modulus for implicit FE analysis

General Anisotropic Yield Function:

The yield function is defined as:

\[\Phi(\boldsymbol{\sigma}, \mathbf{X}, p) = \sqrt{(\boldsymbol{\sigma} - \mathbf{X}) : \mathbf{H} : (\boldsymbol{\sigma} - \mathbf{X})} - R(p) - \sigma_Y \leq 0 \]
where \( \mathbf{H} \) is a general fourth-order anisotropy tensor with up to 21 independent components for a fully anisotropic material.

**Anisotropy Tensor \( \mathbf{H} \):**

In Voigt notation, the anisotropy tensor is represented as a 6×6 symmetric matrix:

\[\begin{split}\mathbf{H} = \begin{bmatrix} H_{11} & H_{12} & H_{13} & H_{14} & H_{15} & H_{16} \\ H_{12} & H_{22} & H_{23} & H_{24} & H_{25} & H_{26} \\ H_{13} & H_{23} & H_{33} & H_{34} & H_{35} & H_{36} \\ H_{14} & H_{24} & H_{34} & H_{44} & H_{45} & H_{46} \\ H_{15} & H_{25} & H_{35} & H_{45} & H_{55} & H_{56} \\ H_{16} & H_{26} & H_{36} & H_{46} & H_{56} & H_{66} \end{bmatrix} \end{split}\]

Special Cases:

  • von Mises (isotropic): \( \mathbf{H} = \mathbf{I}_{dev} \) (deviatoric identity)

  • Hill orthotropic: \( H_{ij} = 0 \) for \( i \leq 3 < j \) or \( j \leq 3 < i \)

  • Transversely isotropic: 5 independent parameters

Armstrong-Frederick Backstress Evolution:

Each backstress component evolves according to:

\[\dot{\mathbf{X}}_i = \frac{2}{3} C_i \dot{\boldsymbol{\varepsilon}}^p - D_i \mathbf{X}_i \dot{p} \]

Isotropic Hardening (Voce Law):

\[R(p) = \sum_{j=1}^{N_{iso}} Q_j (1 - e^{-b_j p}) \]

Applications:

This model is suited for:

  • Single crystals with complex slip system interactions

  • Highly textured polycrystals requiring full anisotropy

  • Materials with non-orthotropic symmetry (monoclinic, triclinic)

  • Advanced composite materials with complex anisotropy

State Variables (statev):

Total: \( n_{statev} = 1 + 1 + 6 + 6 \times N_{kin} \)

Index

Symbol

Description

Units

statev[0]

\( T_{init} \)

Initial temperature

Temperature

statev[1]

\( p \)

Accumulated plastic strain

Strain

statev[2:7]

\( \boldsymbol{\varepsilon}^p \)

Plastic strain tensor (Voigt)

Strain

statev[8+6i:13+6i]

\( \mathbf{X}_i \)

Backstress i (Voigt)

Stress

References:

  • Chaboche, J. L. (1986). “Time-independent constitutive theories for cyclic plasticity.” Int. J. Plasticity, 2(2), 149-188.

  • Barlat, F., et al. (2005). “Linear transformation-based anisotropic yield functions.” Int. J. Plasticity, 21(5), 1009-1039.

See also

umat_hill_chaboche_CCP() for orthotropic Hill + Chaboche

See also

umat_plasticity_chaboche_CCP() for isotropic von Mises + Chaboche

Note

The anisotropy tensor must be positive semi-definite for convexity

Note

For orthotropic materials, consider using Hill_chaboche_ccp instead

Note

Parameter identification requires extensive multiaxial testing

Parameters:

  • Etot – Total strain tensor at beginning of increment (Voigt notation: 6×1 vector)

  • DEtot – Strain increment tensor (Voigt notation: 6×1 vector)

  • sigma – Stress tensor (Voigt notation: 6×1 vector) [output]

  • Lt – Consistent tangent modulus (6×6 matrix) [output]

  • L – Elastic stiffness tensor (6×6 matrix) [output]

  • sigma_in – Internal stress for explicit solvers (6×1 vector) [output]

  • DR – Rotation increment matrix (3×3) for objective integration

  • nprops – Number of material properties

  • props – Material properties vector

  • nstatev – Number of state variables

  • statev – State variables vector [input/output]

  • T – Temperature at beginning of increment

  • DT – Temperature increment

  • Time – Time at beginning of increment

  • DTime – Time increment

  • Wm – Total mechanical work [output]

  • Wm_r – Recoverable (elastic) work [output]

  • Wm_ir – Irrecoverable work stored in hardening [output]

  • Wm_d – Dissipated (plastic) work [output]

  • ndi – Number of direct stress components (typically 3)

  • nshr – Number of shear stress components (typically 3)

  • start – Flag indicating first increment (true) or continuation (false)

  • solver_type – Solver type: 0=implicit, 1=explicit, 2=dynamic implicit

  • tnew_dt – Suggested new time step size for adaptive time stepping [output]

void umat_dfa_chaboche_CCP(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, arma::vec &sigma_in, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt)

Elastic-plastic model with Distortional/Directional Forming Anisotropy (DFA) and Chaboche hardening.

This function implements an advanced plasticity model that captures the evolution of plastic anisotropy during deformation (distortional hardening). Unlike classical models where the yield surface only translates (kinematic) or expands (isotropic), the DFA model allows the yield surface shape to evolve with plastic straining.

Key Features:

  • Distortional hardening: yield surface shape evolves with deformation

  • Multiple Armstrong-Frederick backstresses for kinematic hardening

  • Optional Voce-type isotropic hardening

  • Captures cross-hardening and directional effects

  • Convex Cutting Plane algorithm for return mapping

Physical Motivation:

Classical hardening models (isotropic + kinematic) cannot capture:

  • Cross-hardening: Yield stress in one direction affected by prior loading in another

  • Permanent softening: Reduced yield in tension after prior compression (or vice versa)

  • Yield surface distortion: Changes in yield locus shape observed in metals after pre-straining

The DFA model addresses these limitations by allowing the anisotropy tensor \( \mathbf{H} \) to evolve with plastic deformation.

Yield Function with Evolving Anisotropy:

\[\Phi = \sqrt{(\boldsymbol{\sigma} - \mathbf{X}) : \mathbf{H}(p, \boldsymbol{\varepsilon}^p) : (\boldsymbol{\sigma} - \mathbf{X})} - R(p) - \sigma_Y \leq 0 \]

where the anisotropy tensor \( \mathbf{H} \) evolves according to:

\[\dot{\mathbf{H}} = f(\mathbf{H}, \boldsymbol{\sigma}, \dot{\boldsymbol{\varepsilon}}^p, \dot{p}) \]

Directional Hardening:

The model tracks directional hardening through internal variables that remember the loading history orientation. This allows capturing:

  • Latent hardening: Increased yield stress in inactive slip systems

  • Texture evolution effects: Macroscopic manifestation of microstructural changes

  • Anisotropic damage accumulation: Direction-dependent material degradation

Applications:

  • Sheet metal forming with complex strain paths

  • Materials exhibiting strong cross-hardening (aluminum alloys, TWIP steels)

  • Multi-pass forming operations with changing strain directions

  • Accurate springback prediction in automotive stamping

State Variables (statev):

Index

Symbol

Description

Units

statev[0]

\( T_{init} \)

Initial temperature

Temperature

statev[1]

\( p \)

Accumulated plastic strain

Strain

statev[2:7]

\( \boldsymbol{\varepsilon}^p \)

Plastic strain tensor (Voigt)

Strain

statev[8+6i:13+6i]

\( \mathbf{X}_i \)

Backstress i (Voigt)

Stress

statev[…]

\( H_{ij} \)

Evolving anisotropy components

-

References:

  • Barlat, F., et al. (2011). “An alternative to kinematic hardening in classical plasticity.” Int. J. Plasticity, 27(9), 1309-1327.

  • Teodosiu, C., & Hu, Z. (1998). “Microstructure in the continuum modeling of plastic anisotropy.” Proc. 19th Riso Int. Symp., 149-168.

  • Holmedal, B. (2019). “Bauschinger effect modelled by yield surface distortions.” Int. J. Plasticity, 123, 86-100.

See also

umat_plasticity_chaboche_CCP() for standard Chaboche without distortional hardening

See also

umat_ani_chaboche_CCP() for fixed anisotropy Chaboche

Note

Requires careful parameter identification from cross-hardening tests

Note

Computationally more expensive than standard Chaboche due to evolving anisotropy

Note

Useful for complex forming simulations with strain path changes

Parameters:

  • Etot – Total strain tensor at beginning of increment (Voigt notation: 6×1 vector)

  • DEtot – Strain increment tensor (Voigt notation: 6×1 vector)

  • sigma – Stress tensor (Voigt notation: 6×1 vector) [output]

  • Lt – Consistent tangent modulus (6×6 matrix) [output]

  • L – Elastic stiffness tensor (6×6 matrix) [output]

  • sigma_in – Internal stress for explicit solvers (6×1 vector) [output]

  • DR – Rotation increment matrix (3×3) for objective integration

  • nprops – Number of material properties

  • props – Material properties vector

  • nstatev – Number of state variables

  • statev – State variables vector [input/output]

  • T – Temperature at beginning of increment

  • DT – Temperature increment

  • Time – Time at beginning of increment

  • DTime – Time increment

  • Wm – Total mechanical work [output]

  • Wm_r – Recoverable (elastic) work [output]

  • Wm_ir – Irrecoverable work stored in hardening [output]

  • Wm_d – Dissipated (plastic) work [output]

  • ndi – Number of direct stress components (typically 3)

  • nshr – Number of shear stress components (typically 3)

  • start – Flag indicating first increment (true) or continuation (false)

  • solver_type – Solver type: 0=implicit, 1=explicit, 2=dynamic implicit

  • tnew_dt – Suggested new time step size for adaptive time stepping [output]

void umat_generic_chaboche_CCP(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, arma::vec &sigma_in, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt)

Generic elastic-plastic framework with flexible yield criterion and Chaboche kinematic hardening.

This function implements a generic plasticity framework that allows for flexible configuration of the yield criterion while using Chaboche-type hardening. It serves as a template for implementing various plasticity models with different yield functions.

Key Features:

  • Configurable yield criterion through function pointers or template parameters

  • Multiple Armstrong-Frederick backstresses for nonlinear kinematic hardening

  • Optional Voce-type isotropic hardening

  • Convex Cutting Plane algorithm for return mapping

  • Consistent tangent modulus for implicit FE analysis

Generic Yield Function:

The yield function takes the general form:

\[\Phi(\boldsymbol{\sigma}, \mathbf{X}, p, \boldsymbol{\alpha}) = f(\boldsymbol{\sigma} - \mathbf{X}) - R(p) - \sigma_Y \leq 0 \]
where:

  • \( f(\boldsymbol{\eta}) \) is a general equivalent stress function (configurable)

  • \( \boldsymbol{\alpha} \) represents additional internal variables

Supported Yield Criteria:

The generic framework can accommodate:

  • von Mises: \( f = \sqrt{\frac{3}{2} \mathbf{s}:\mathbf{s}} \)

  • Hill (1948): Quadratic orthotropic

  • Barlat Yld2000-2d: Advanced anisotropic for sheet metals

  • Drucker-Prager: Pressure-dependent for geomaterials

  • Hosford: Non-quadratic isotropic

  • Custom user-defined: Through function interface

Armstrong-Frederick Backstress Evolution:

Each backstress evolves according to:

\[\dot{\mathbf{X}}_i = \frac{2}{3} C_i \dot{\boldsymbol{\varepsilon}}^p - D_i \mathbf{X}_i \dot{p} \]

Isotropic Hardening (Voce Law):

\[R(p) = \sum_{j=1}^{N_{iso}} Q_j (1 - e^{-b_j p}) \]

Framework Architecture:

The generic model separates:

  1. Yield criterion evaluation: \( f(\boldsymbol{\eta}) \) and \( \partial f / \partial \boldsymbol{\eta} \)

  2. Hardening laws: Evolution of \( \mathbf{X}_i \) and \( R \)

  3. Return mapping: CCP algorithm for stress update

  4. Consistent tangent: Algorithmic modulus computation

This allows easy extension to new yield criteria while reusing the hardening and integration framework.

State Variables (statev):

Index

Symbol

Description

Units

statev[0]

\( T_{init} \)

Initial temperature

Temperature

statev[1]

\( p \)

Accumulated plastic strain

Strain

statev[2:7]

\( \boldsymbol{\varepsilon}^p \)

Plastic strain tensor (Voigt)

Strain

statev[8+6i:13+6i]

\( \mathbf{X}_i \)

Backstress i (Voigt)

Stress

statev[…]

\( \boldsymbol{\alpha} \)

Additional criterion-specific variables

-

References:

  • Chaboche, J. L. (2008). “A review of some plasticity and viscoplasticity constitutive theories.” Int. J. Plasticity, 24(10), 1642-1693.

  • Simo, J. C., & Hughes, T. J. R. (1998). Computational Inelasticity. Springer.

See also

umat_plasticity_chaboche_CCP() for standard von Mises Chaboche

See also

umat_hill_chaboche_CCP() for Hill anisotropic Chaboche

See also

umat_ani_chaboche_CCP() for general anisotropic Chaboche

Note

Use this as a base for implementing new yield criteria

Note

Requires implementation of yield function and its derivatives

Note

See specialized models (Hill, Ani) for specific applications

Parameters:

  • Etot – Total strain tensor at beginning of increment (Voigt notation: 6×1 vector)

  • DEtot – Strain increment tensor (Voigt notation: 6×1 vector)

  • sigma – Stress tensor (Voigt notation: 6×1 vector) [output]

  • Lt – Consistent tangent modulus (6×6 matrix) [output]

  • L – Elastic stiffness tensor (6×6 matrix) [output]

  • sigma_in – Internal stress for explicit solvers (6×1 vector) [output]

  • DR – Rotation increment matrix (3×3) for objective integration

  • nprops – Number of material properties

  • props – Material properties vector (criterion-dependent)

  • nstatev – Number of state variables

  • statev – State variables vector [input/output]

  • T – Temperature at beginning of increment

  • DT – Temperature increment

  • Time – Time at beginning of increment

  • DTime – Time increment

  • Wm – Total mechanical work [output]

  • Wm_r – Recoverable (elastic) work [output]

  • Wm_ir – Irrecoverable work stored in hardening [output]

  • Wm_d – Dissipated (plastic) work [output]

  • ndi – Number of direct stress components (typically 3)

  • nshr – Number of shear stress components (typically 3)

  • start – Flag indicating first increment (true) or continuation (false)

  • solver_type – Solver type: 0=implicit, 1=explicit, 2=dynamic implicit

  • tnew_dt – Suggested new time step size for adaptive time stepping [output]

void umat_hill_chaboche_CCP(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, arma::vec &sigma_in, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt)

Elastic-plastic model combining Hill anisotropic yield criterion with Chaboche kinematic hardening.

This function implements an advanced plasticity model combining the Hill (1948) anisotropic yield criterion with Chaboche nonlinear kinematic hardening. This allows modeling of anisotropic materials under cyclic loading conditions. The model features:

  • Hill (1948) quadratic anisotropic yield criterion

  • Multiple Armstrong-Frederick backstresses for nonlinear kinematic hardening

  • Optional Voce-type isotropic hardening

  • Convex Cutting Plane algorithm for return mapping

  • Consistent tangent modulus for implicit FE analysis

Constitutive Equations:

The anisotropic yield function with kinematic hardening is defined as:

\[\Phi(\boldsymbol{\sigma}, \mathbf{X}, p) = \sigma_{eq}^{Hill}(\boldsymbol{\sigma} - \mathbf{X}) - R(p) - \sigma_Y \leq 0 \]
where the Hill equivalent stress of the shifted stress tensor is:
\[\sigma_{eq}^{Hill}(\boldsymbol{\eta}) = \sqrt{\boldsymbol{\eta} : \mathbf{P}^{Hill} : \boldsymbol{\eta}} \]
with \( \boldsymbol{\eta} = \boldsymbol{\sigma} - \mathbf{X} \) being the effective (shifted) stress.

Hill Anisotropic Yield Surface:

The Hill projection tensor \( \mathbf{P}^{Hill} \) is defined by six parameters \( (F, G, H, L, M, N) \):

\[(\sigma_{eq}^{Hill})^2 = F(\eta_{22} - \eta_{33})^2 + G(\eta_{33} - \eta_{11})^2 + H(\eta_{11} - \eta_{22})^2 + 2L\eta_{23}^2 + 2M\eta_{13}^2 + 2N\eta_{12}^2 \]

Armstrong-Frederick Backstress Evolution:

Each backstress component evolves with dynamic recovery:

\[\dot{\mathbf{X}}_i = \frac{2}{3} C_i \dot{\boldsymbol{\varepsilon}}^p - D_i \mathbf{X}_i \dot{p} \]
where the total backstress is: \( \mathbf{X} = \sum_{i=1}^{N_{kin}} \mathbf{X}_i \)

Flow Rule:

The plastic strain rate follows the associative flow rule:

\[\dot{\boldsymbol{\varepsilon}}^p = \dot{p} \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} = \dot{p} \frac{\mathbf{P}^{Hill} : \boldsymbol{\eta}}{\sigma_{eq}^{Hill}(\boldsymbol{\eta})} \]

Applications:

This model is suited for:

  • Rolled sheet metals under cyclic loading

  • Textured polycrystalline materials with Bauschinger effect

  • Drawn tubes and wires with anisotropic yield and kinematic hardening

  • Any orthotropic material exhibiting cyclic plasticity

Material Parameters (props):

Index

Symbol

Description

Units

props[0]

\( E \)

Young’s modulus

Stress

props[1]

\( \nu \)

Poisson’s ratio

-

props[2]

\( \alpha \)

Thermal expansion coefficient

1/Temperature

props[3]

\( \sigma_Y \)

Initial yield stress

Stress

props[4]

\( N_{iso} \)

Number of isotropic hardening terms

-

props[5]

\( N_{kin} \)

Number of kinematic hardening terms

-

props[6]

\( F \)

Hill parameter F

1/Stress²

props[7]

\( G \)

Hill parameter G

1/Stress²

props[8]

\( H \)

Hill parameter H

1/Stress²

props[9]

\( L \)

Hill parameter L

1/Stress²

props[10]

\( M \)

Hill parameter M

1/Stress²

props[11]

\( N \)

Hill parameter N

1/Stress²

props[12+2j]

\( Q_j \)

Isotropic saturation stress (j-th term)

Stress

props[13+2j]

\( b_j \)

Isotropic hardening rate (j-th term)

1/Strain

props[12+2N_{iso}+2i]

\( C_i \)

Kinematic modulus (i-th backstress)

Stress

props[13+2N_{iso}+2i]

\( D_i \)

Dynamic recovery (i-th backstress)

-

State Variables (statev):

Total: \( n_{statev} = 1 + 1 + 6 + 6 \times N_{kin} \)

Index

Symbol

Description

Units

statev[0]

\( T_{init} \)

Initial temperature

Temperature

statev[1]

\( p \)

Accumulated plastic strain

Strain

statev[2:7]

\( \boldsymbol{\varepsilon}^p \)

Plastic strain tensor (Voigt)

Strain

statev[8+6i:13+6i]

\( \mathbf{X}_i \)

Backstress i (Voigt)

Stress

References:

  • Hill, R. (1948). “A theory of the yielding and plastic flow of anisotropic metals.” Proc. Roy. Soc. A, 193(1033), 281-297.

  • Chaboche, J. L. (1986). “Time-independent constitutive theories for cyclic plasticity.” Int. J. Plasticity, 2(2), 149-188.

  • Barlat, F., et al. (2005). “Linear transformation-based anisotropic yield functions.” Int. J. Plasticity, 21(5), 1009-1039.

See also

umat_plasticity_hill_isoh_CCP() for Hill with isotropic hardening only

See also

umat_plasticity_chaboche_CCP() for isotropic (von Mises) Chaboche model

Note

Material axes must align with global coordinates or use rotation tensors

Note

For isotropic yield with Chaboche hardening, use umat_plasticity_chaboche_CCP

Note

Captures both anisotropic yielding and Bauschinger effect

Parameters:

  • Etot – Total strain tensor at beginning of increment (Voigt notation: 6×1 vector)

  • DEtot – Strain increment tensor (Voigt notation: 6×1 vector)

  • sigma – Stress tensor (Voigt notation: 6×1 vector) [output]

  • Lt – Consistent tangent modulus (6×6 matrix) [output]

  • L – Elastic stiffness tensor (6×6 matrix) [output]

  • sigma_in – Internal stress for explicit solvers (6×1 vector) [output]

  • DR – Rotation increment matrix (3×3) for objective integration

  • nprops – Number of material properties

  • props – Material properties vector (see table above)

  • nstatev – Number of state variables

  • statev – State variables vector (see table above) [input/output]

  • T – Temperature at beginning of increment

  • DT – Temperature increment

  • Time – Time at beginning of increment

  • DTime – Time increment

  • Wm – Total mechanical work [output]

  • Wm_r – Recoverable (elastic) work [output]

  • Wm_ir – Irrecoverable work stored in hardening [output]

  • Wm_d – Dissipated (plastic) work [output]

  • ndi – Number of direct stress components (typically 3)

  • nshr – Number of shear stress components (typically 3)

  • start – Flag indicating first increment (true) or continuation (false)

  • solver_type – Solver type: 0=implicit, 1=explicit, 2=dynamic implicit

  • tnew_dt – Suggested new time step size for adaptive time stepping [output]

void umat_plasticity_hill_isoh_CCP(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt)

Elastic-plastic constitutive model with Hill anisotropic yield criterion and isotropic hardening solved by the Convex Cutting Plane (CCP) algorithm.

This function implements the Hill anisotropic plasticity model for small and finite strain analysis of orthotropic materials. The model features:

  • Hill (1948) quadratic anisotropic yield criterion

  • Isotropic hardening with power law

  • Associative flow rule

  • Convex Cutting Plane algorithm for return mapping

  • Thermal expansion effects

  • Consistent tangent modulus for implicit FE analysis

Constitutive Equations:

The Hill anisotropic yield function is defined as:

\[\Phi(\boldsymbol{\sigma}, p) = \sigma_{eq}^{Hill} - H_p(p) - \sigma_Y \leq 0 \]
where:
\[\sigma_{eq}^{Hill} = \sqrt{F(\sigma_{22} - \sigma_{33})^2 + G(\sigma_{33} - \sigma_{11})^2 + H(\sigma_{11} - \sigma_{22})^2 + 2L\sigma_{23}^2 + 2M\sigma_{13}^2 + 2N\sigma_{12}^2} \]

Hill Anisotropy Parameters:

The six Hill parameters \( (F, G, H, L, M, N) \) define the material’s plastic anisotropy. For orthotropic symmetry with principal axes aligned with material axes:

  • \( F, G, H \) control yielding in normal stress components

  • \( L, M, N \) control yielding in shear stress components

Relationship to Yield Stress Ratios:

The Hill parameters can be related to yield stresses in different directions:

\[F = \frac{1}{2}\left(\frac{1}{\sigma_{Y22}^2} + \frac{1}{\sigma_{Y33}^2} - \frac{1}{\sigma_{Y11}^2}\right) \]
\[G = \frac{1}{2}\left(\frac{1}{\sigma_{Y33}^2} + \frac{1}{\sigma_{Y11}^2} - \frac{1}{\sigma_{Y22}^2}\right) \]
\[H = \frac{1}{2}\left(\frac{1}{\sigma_{Y11}^2} + \frac{1}{\sigma_{Y22}^2} - \frac{1}{\sigma_{Y33}^2}\right) \]
\[L = \frac{3}{2\sigma_{Y23}^2}, \quad M = \frac{3}{2\sigma_{Y13}^2}, \quad N = \frac{3}{2\sigma_{Y12}^2} \]
where \( \sigma_{Yii} \) is the yield stress in uniaxial loading along direction \( i \), and \( \sigma_{Yij} \) is the yield stress in pure shear in the \( ij \) plane.

Special Case - von Mises (Isotropic):

For isotropic materials, the Hill criterion reduces to von Mises when:

\[F = G = H = \frac{1}{2}, \quad L = M = N = \frac{3}{2} \]

Isotropic Hardening (Power Law):

The isotropic hardening follows a power law:

\[H_p(p) = k \cdot p^m \]
where:

  • \( k \) is the hardening coefficient

  • \( m \) is the hardening exponent

  • \( p = \int_0^t \dot{p} \, dt \) is the accumulated plastic strain

Plastic Flow Rule:

The plastic strain rate follows the associative flow rule:

\[\dot{\boldsymbol{\varepsilon}}^p = \dot{p} \mathbf{n} \]
where:
\[\mathbf{n} = \frac{\partial \sigma_{eq}^{Hill}}{\partial \boldsymbol{\sigma}} = \frac{1}{\sigma_{eq}^{Hill}} \mathbf{P}^{Hill} : \boldsymbol{\sigma} \]
is the flow direction normal to the Hill yield surface, and \( \mathbf{P}^{Hill} \) is the Hill projection tensor.

Hill Projection Tensor (Voigt Notation):

\[\begin{split}\mathbf{P}^{Hill} = \begin{bmatrix} G+H & -H & -G & 0 & 0 & 0 \\ -H & H+F & -F & 0 & 0 & 0 \\ -G & -F & F+G & 0 & 0 & 0 \\ 0 & 0 & 0 & 2N & 0 & 0 \\ 0 & 0 & 0 & 0 & 2M & 0 \\ 0 & 0 & 0 & 0 & 0 & 2L \end{bmatrix} \end{split}\]

Incremental Form (CCP Algorithm):

For an increment \( \Delta p \), the stress is updated as:

\[\boldsymbol{\sigma}_{n+1} = \boldsymbol{\sigma}^{trial} - \Delta p \left( \mathbf{L} : \mathbf{n} + \frac{\partial \mathbf{n}}{\partial \boldsymbol{\sigma}} : \boldsymbol{\sigma} \right) \]

The CCP algorithm solves for \( \Delta p \) using the Fischer-Burmeister complementarity condition:

\[FB(\Delta p, \Phi) = \sqrt{(\Delta p)^2 + \Phi^2} - \Delta p - \Phi = 0 \]

Consistent Tangent Modulus:

The algorithmic tangent modulus is:

\[\mathbf{L}_t = \mathbf{L} - \frac{(\mathbf{L}:\mathbf{n}) \otimes (\mathbf{n}:\mathbf{L})}{\mathbf{n}:\mathbf{L}:\mathbf{n} + H'} \]
where \( H' = \frac{dH_p}{dp} = k \cdot m \cdot p^{m-1} \) is the hardening modulus.

Material Parameters (props):

Index

Symbol

Description

Units

Typical Range

props[0]

\( E \)

Young’s modulus

Stress

50-500 GPa

props[1]

\( \nu \)

Poisson’s ratio

-

0.2-0.45

props[2]

\( \alpha \)

Coefficient of thermal expansion

1/Temperature

1e-6 to 1e-4 /K

props[3]

\( \sigma_Y \)

Initial yield stress

Stress

100-1000 MPa

props[4]

\( k \)

Hardening coefficient

Stress

100-2000 MPa

props[5]

\( m \)

Hardening exponent

-

0.05-0.5

props[6]

\( F \)

Hill parameter F

1/Stress²

0-10

props[7]

\( G \)

Hill parameter G

1/Stress²

0-10

props[8]

\( H \)

Hill parameter H

1/Stress²

0-10

props[9]

\( L \)

Hill parameter L

1/Stress²

0-10

props[10]

\( M \)

Hill parameter M

1/Stress²

0-10

props[11]

\( N \)

Hill parameter N

1/Stress²

0-10

Constraint on Hill Parameters:

  • For physical consistency: \( F + G + H > 0 \)

  • For von Mises isotropy: \( F = G = H = 0.5, L = M = N = 1.5 \)

  • For plane stress: Ensure appropriate ratios for the stress state

State Variables (statev):

Total state variables required: \( n_{statev} = 8 \)

Index

Symbol

Description

Units

statev[0]

\( T_{init} \)

Initial/reference temperature

Temperature

statev[1]

\( p \)

Accumulated plastic strain

Strain

statev[2]

\( \varepsilon^p_{11} \)

Plastic strain component 11

Strain

statev[3]

\( \varepsilon^p_{22} \)

Plastic strain component 22

Strain

statev[4]

\( \varepsilon^p_{33} \)

Plastic strain component 33

Strain

statev[5]

\( \varepsilon^p_{12} \)

Plastic strain component 12 (engineering)

Strain

statev[6]

\( \varepsilon^p_{13} \)

Plastic strain component 13 (engineering)

Strain

statev[7]

\( \varepsilon^p_{23} \)

Plastic strain component 23 (engineering)

Strain 

// Example usage: Aluminum rolled sheet (transverse isotropy in 1-2 plane)
vec props(12);
props(0) = 70000;       // E = 70 GPa
props(1) = 0.33;        // nu = 0.33
props(2) = 2.3e-5;      // alpha = 23e-6 /K
props(3) = 200;         // sigma_Y = 200 MPa
props(4) = 400;         // k = 400 MPa
props(5) = 0.15;        // m = 0.15

// Hill parameters for rolled sheet (stronger in rolling direction)
// Assume: sigma_Y11 = 200 MPa (rolling), sigma_Y22 = 220 MPa (transverse), sigma_Y33 = 210 MPa
// sigma_Y12 = sigma_Y13 = sigma_Y23 = 115.5 MPa (shear, sqrt(3)*sigma_Y/3)
props(6) = 0.00115;     // F
props(7) = 0.00116;     // G
props(8) = 0.00126;     // H
props(9) = 0.0130;      // L = 1.5 for isotropy
props(10) = 0.0130;     // M = 1.5 for isotropy
props(11) = 0.0130;     // N = 1.5 for isotropy

vec statev = zeros(8);
statev(0) = 20.0;  // Reference temperature 20°C

vec Etot = {0.003, 0.0, 0.0, 0.0, 0.0, 0.0};  // 0.3% uniaxial strain
vec DEtot = {0.0001, 0.0, 0.0, 0.0, 0.0, 0.0};
vec sigma = zeros(6);
mat Lt = zeros(6,6);
mat DR = eye(3,3);

umat_plasticity_hill_isoh_CCP(Etot, DEtot, sigma, Lt, DR,
                              12, props, 8, statev, 20.0, 0.0, 0.0, 1.0,
                              Wm, Wm_r, Wm_ir, Wm_d, 3, 3, false, tnew_dt);

// Expected: Different yield stress depending on loading direction

See also

L_iso() for isotropic elastic stiffness construction

See also

Ireal() for real identity tensor

See also

Fischer_Burmeister() for complementarity solver

See also

Eq_stress_Hill() for Hill equivalent Stress Functions computation

See also

denom_FB_Hill() for CCP denominator with Hill criterion

References:

  • Hill, R. (1948). “A theory of the yielding and plastic flow of anisotropic metals.” Proceedings of the Royal Society of London A, 193(1033), 281-297.

  • Hill, R. (1950). The Mathematical Theory of Plasticity. Oxford University Press.

  • Barlat, F., & Lian, J. (1989). “Plastic behavior and stretchability of sheet metals. Part I: A yield function for orthotropic sheets under plane stress conditions.” International Journal of Plasticity, 5(1), 51-66.

  • Banabic, D. (2010). Sheet Metal Forming Processes: Constitutive Modelling and Numerical Simulation. Springer.

  • Ortiz, M., & Simo, J. C. (1986). “An analysis of a new class of integration algorithms for elastoplastic constitutive relations.” International Journal for Numerical Methods in Engineering, 23(3), 353-366.

Note

The Hill criterion is widely used for metals with crystallographic texture

Note

Common applications: rolled sheets, drawn wires, forged components

Note

Parameter identification requires yield stress measurements in multiple directions

Note

For fiber composites, consider more advanced criteria (Tsai-Wu, Hoffman)

Note

The model assumes plastic incompressibility (plastic volume change = 0)

Note

Material axes must align with the global coordinate system or use rotation

Parameters:

  • Etot – Total strain tensor at beginning of increment (Voigt notation: 6×1 vector)

  • DEtot – Strain increment tensor (Voigt notation: 6×1 vector)

  • sigma – Stress tensor (Voigt notation: 6×1 vector) [output]

  • Lt – Consistent tangent modulus \( \mathbf{L}_t \) (6×6 matrix) [output]

  • DR – Rotation increment matrix (3×3) for objective integration

  • nprops – Number of material properties

  • props – Material properties vector (see table above)

  • nstatev – Number of state variables

  • statev – State variables vector (see table above) [input/output]

  • T – Temperature at beginning of increment

  • DT – Temperature increment

  • Time – Time at beginning of increment

  • DTime – Time increment

  • Wm – Total mechanical work [output]

  • Wm_r – Recoverable (elastic) work [output]

  • Wm_ir – Irrecoverable work stored in plastic deformation [output]

  • Wm_d – Dissipated (plastic) work [output]

  • ndi – Number of direct stress components (typically 3)

  • nshr – Number of shear stress components (typically 3)

  • start – Flag indicating first increment (true) or continuation (false)

  • tnew_dt – Suggested new time step size for adaptive time stepping [output]

void umat_plasticity_hill_isoh_CCP_N(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt)

Elastic-plastic model with Hill anisotropic criterion and N isotropic hardening terms.

This function implements the Hill anisotropic plasticity model with multiple isotropic hardening terms for enhanced flexibility in capturing complex hardening behavior. The model features:

  • Hill (1948) quadratic anisotropic yield criterion

  • Multiple Voce-type isotropic hardening terms: \( R(p) = \sum_{i=1}^{N} Q_i (1 - e^{-b_i p}) \)

  • Associative flow rule

  • Convex Cutting Plane algorithm for return mapping

  • Consistent tangent modulus for implicit FE analysis

Constitutive Equations:

The Hill anisotropic yield function is defined as:

\[\Phi(\boldsymbol{\sigma}, p) = \sigma_{eq}^{Hill} - R(p) - \sigma_Y \leq 0 \]
where:
\[\sigma_{eq}^{Hill} = \sqrt{F(\sigma_{22} - \sigma_{33})^2 + G(\sigma_{33} - \sigma_{11})^2 + H(\sigma_{11} - \sigma_{22})^2 + 2L\sigma_{23}^2 + 2M\sigma_{13}^2 + 2N\sigma_{12}^2} \]

Multiple Isotropic Hardening (Voce Law):

The isotropic hardening is the sum of N exponential terms:

\[R(p) = \sum_{i=1}^{N_{iso}} Q_i \left( 1 - e^{-b_i p} \right) \]
where:

  • \( Q_i \) is the saturation stress of the i-th hardening term

  • \( b_i \) is the hardening rate parameter of the i-th term

  • \( p \) is the accumulated plastic strain

Hardening Modulus:

The derivative of the isotropic hardening function:

\[H_{iso} = \frac{dR}{dp} = \sum_{i=1}^{N_{iso}} Q_i b_i e^{-b_i p} \]

Advantages of Multiple Hardening Terms:

  • Better fit to experimental stress-strain curves

  • Captures both rapid initial hardening and gradual saturation

  • Different terms can represent different physical mechanisms:

    • Fast term (high b): dislocation pile-up, initial work hardening

    • Slow term (low b): long-range backstress development, texture evolution

Material Parameters (props):

Index

Symbol

Description

Units

props[0]

\( E \)

Young’s modulus

Stress

props[1]

\( \nu \)

Poisson’s ratio

-

props[2]

\( \alpha \)

Thermal expansion coefficient

1/Temperature

props[3]

\( \sigma_Y \)

Initial yield stress

Stress

props[4]

\( N_{iso} \)

Number of isotropic hardening terms

-

props[5]

\( F \)

Hill parameter F

1/Stress²

props[6]

\( G \)

Hill parameter G

1/Stress²

props[7]

\( H \)

Hill parameter H

1/Stress²

props[8]

\( L \)

Hill parameter L

1/Stress²

props[9]

\( M \)

Hill parameter M

1/Stress²

props[10]

\( N \)

Hill parameter N

1/Stress²

props[11+2j]

\( Q_j \)

Saturation stress of j-th term

Stress

props[12+2j]

\( b_j \)

Hardening rate of j-th term

1/Strain

State Variables (statev):

Total state variables required: \( n_{statev} = 8 \)

Index

Symbol

Description

Units

statev[0]

\( T_{init} \)

Initial temperature

Temperature

statev[1]

\( p \)

Accumulated plastic strain

Strain

statev[2]

\( \varepsilon^p_{11} \)

Plastic strain component 11

Strain

statev[3]

\( \varepsilon^p_{22} \)

Plastic strain component 22

Strain

statev[4]

\( \varepsilon^p_{33} \)

Plastic strain component 33

Strain

statev[5]

\( \varepsilon^p_{12} \)

Plastic strain component 12

Strain

statev[6]

\( \varepsilon^p_{13} \)

Plastic strain component 13

Strain

statev[7]

\( \varepsilon^p_{23} \)

Plastic strain component 23

Strain

References:

  • Hill, R. (1948). “A theory of the yielding and plastic flow of anisotropic metals.” Proc. Roy. Soc. A, 193(1033), 281-297.

  • Voce, E. (1948). “The relationship between stress and strain for homogeneous deformation.” J. Inst. Met., 74, 537-562.

  • Barlat, F., et al. (2003). “Plane stress yield function for aluminum alloy sheets.” Int. J. Plasticity, 19(9), 1297-1319.

See also

umat_plasticity_hill_isoh_CCP() for single power-law hardening version

See also

umat_hill_chaboche_CCP() for combined kinematic and isotropic hardening

Note

Typically 2-3 hardening terms are sufficient for most metals

Note

For single power-law hardening, use umat_plasticity_hill_isoh_CCP instead

Note

The Hill parameters must satisfy physical constraints (positive definite)

Parameters:

  • Etot – Total strain tensor at beginning of increment (Voigt notation: 6×1 vector)

  • DEtot – Strain increment tensor (Voigt notation: 6×1 vector)

  • sigma – Stress tensor (Voigt notation: 6×1 vector) [output]

  • Lt – Consistent tangent modulus (6×6 matrix) [output]

  • DR – Rotation increment matrix (3×3) for objective integration

  • nprops – Number of material properties

  • props – Material properties vector (see table above)

  • nstatev – Number of state variables

  • statev – State variables vector (see table above) [input/output]

  • T – Temperature at beginning of increment

  • DT – Temperature increment

  • Time – Time at beginning of increment

  • DTime – Time increment

  • Wm – Total mechanical work [output]

  • Wm_r – Recoverable (elastic) work [output]

  • Wm_ir – Irrecoverable work stored in plastic deformation [output]

  • Wm_d – Dissipated (plastic) work [output]

  • ndi – Number of direct stress components (typically 3)

  • nshr – Number of shear stress components (typically 3)

  • start – Flag indicating first increment (true) or continuation (false)

  • tnew_dt – Suggested new time step size for adaptive time stepping [output]

void umat_plasticity_chaboche_CCP(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, arma::vec &sigma_in, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt)

Elastic-plastic constitutive model with Chaboche kinematic hardening solved by the Convex Cutting Plane (CCP) algorithm.

This function implements the Chaboche unified viscoplasticity model for small and finite strain analysis. The model features:

  • J2 (von Mises) plasticity with associative flow rule

  • Multiple Armstrong-Frederick backstresses for nonlinear kinematic hardening

  • Optional isotropic hardening with Voce law

  • Dynamic recovery (recall terms) in backstress evolution

  • Convex Cutting Plane algorithm for return mapping

  • Thermal expansion effects

  • Consistent tangent modulus for implicit FE analysis

Constitutive Equations:

The yield function is defined as:

\[\Phi(\boldsymbol{\sigma}, \mathbf{X}, R) = \sigma_{eq}(\boldsymbol{\sigma} - \mathbf{X}) - R - \sigma_Y \leq 0 \]
where:

  • \( \sigma_{eq}(\boldsymbol{\eta}) = \sqrt{\frac{3}{2} \boldsymbol{\eta}_{dev} : \boldsymbol{\eta}_{dev}} \) is the von Mises equivalent stress

  • \( \boldsymbol{\eta} = \boldsymbol{\sigma} - \mathbf{X} \) is the shifted (effective) stress tensor

  • \( \mathbf{X} = \sum_{i=1}^{N_{kin}} \mathbf{X}_i \) is the total backstress (sum of individual backstresses)

  • \( R \) is the isotropic hardening stress

  • \( \sigma_Y \) is the initial yield stress

Armstrong-Frederick Backstress Evolution:

Each backstress component evolves according to the Armstrong-Frederick law with dynamic recovery:

\[\dot{\mathbf{X}}_i = \frac{2}{3} C_i \dot{\boldsymbol{\varepsilon}}^p - D_i \mathbf{X}_i \dot{p} \]
where:

  • \( C_i \) is the kinematic hardening modulus of the i-th backstress

  • \( D_i \) is the dynamic recovery (recall) parameter of the i-th backstress

  • \( \dot{\boldsymbol{\varepsilon}}^p \) is the plastic strain rate tensor

  • \( \dot{p} = \sqrt{\frac{2}{3} \dot{\boldsymbol{\varepsilon}}^p : \dot{\boldsymbol{\varepsilon}}^p} \) is the accumulated plastic strain rate

Physical Interpretation:

  • The \( C_i \) term represents strain hardening (backstress growth)

  • The \( D_i \mathbf{X}_i \dot{p} \) term represents dynamic recovery (backstress fading)

  • Multiple backstresses capture different scales of cyclic hardening behavior

  • Typically 2-3 backstress components are used

Isotropic Hardening (Optional):

The isotropic hardening follows the Voce law:

\[R(p) = \sum_{j=1}^{N_{iso}} Q_j \left( 1 - e^{-b_j p} \right) \]
where:

  • \( Q_j \) is the saturation value of the j-th isotropic hardening component

  • \( b_j \) is the hardening rate parameter

  • \( p \) is the accumulated plastic strain

Plastic Flow Rule:

The plastic strain rate follows the associative flow rule:

\[\dot{\boldsymbol{\varepsilon}}^p = \dot{p} \mathbf{n} \]
where:
\[\mathbf{n} = \frac{3}{2} \frac{\boldsymbol{\eta}_{dev}}{\sigma_{eq}(\boldsymbol{\eta})} \]
is the flow direction (normal to the yield surface).

Incremental Form (CCP Algorithm):

For an increment \( \Delta p \), the stress and backstresses are updated as:

\[\boldsymbol{\sigma}_{n+1} = \boldsymbol{\sigma}^{trial} - \Delta p \left( \mathbf{L} : \mathbf{n} + \frac{\partial \mathbf{n}}{\partial \boldsymbol{\sigma}} : \boldsymbol{\sigma} \right) \]
\[\mathbf{X}_{i,n+1} = \mathbf{X}_{i,n} + \frac{2}{3} C_i \Delta \boldsymbol{\varepsilon}^p - D_i \mathbf{X}_{i,n} \Delta p \]

The CCP algorithm solves for \( \Delta p \) using the Fischer-Burmeister complementarity condition:

\[FB(\Delta p, \Phi) = \sqrt{(\Delta p)^2 + \Phi^2} - \Delta p - \Phi = 0 \]

Consistent Tangent Modulus:

The algorithmic tangent modulus accounts for both plastic flow and backstress evolution:

\[\mathbf{L}_t = \mathbf{L} - \frac{(\mathbf{L}:\mathbf{n}) \otimes (\mathbf{n}:\mathbf{L})}{\mathbf{n}:\mathbf{L}:\mathbf{n} + H_{tot}} \]
where:
\[H_{tot} = \sum_{i=1}^{N_{kin}} \frac{C_i}{1 + D_i \Delta p} + \sum_{j=1}^{N_{iso}} Q_j b_j e^{-b_j p} \]
is the total hardening modulus combining kinematic and isotropic contributions.

Material Parameters (props):

For \( N_{iso} \) isotropic hardening terms and \( N_{kin} \) kinematic hardening terms:

Index

Symbol

Description

Units

Typical Range

props[0]

\( E \)

Young’s modulus

Stress

50-500 GPa

props[1]

\( \nu \)

Poisson’s ratio

-

0.2-0.45

props[2]

\( \alpha \)

Coefficient of thermal expansion

1/Temperature

1e-6 to 1e-4 /K

props[3]

\( \sigma_Y \)

Initial yield stress

Stress

100-1000 MPa

props[4]

\( N_{iso} \)

Number of isotropic hardening terms

-

0-3

props[5]

\( N_{kin} \)

Number of kinematic hardening terms

-

1-3

props[6+2j]

\( Q_j \)

Saturation stress of j-th isotropic term

Stress

0-500 MPa

props[7+2j]

\( b_j \)

Hardening rate of j-th isotropic term

1/Strain

1-100

props[6+2N_{iso}+2i]

\( C_i \)

Kinematic modulus of i-th backstress

Stress

10-500 GPa

props[7+2N_{iso}+2i]

\( D_i \)

Dynamic recovery parameter of i-th backstress

-

0-100

Notes on Parameter Selection:

  • For pure kinematic hardening, set \( N_{iso} = 0 \)

  • For pure isotropic hardening, set \( N_{kin} = 0 \) (not typical for Chaboche)

  • First backstress (i=0) typically has high \( C_1 \), low \( D_1 \) (long-range backstress)

  • Second backstress (i=1) typically has moderate \( C_2 \), moderate \( D_2 \)

  • Third backstress (i=2) typically has low \( C_3 \), high \( D_3 \) (short-range backstress)

State Variables (statev):

Total state variables: \( n_{statev} = 1 + 1 + 6 + 6 \times N_{kin} \)

Index

Symbol

Description

Units

statev[0]

\( T_{init} \)

Initial/reference temperature

Temperature

statev[1]

\( p \)

Accumulated plastic strain

Strain

statev[2]

\( \varepsilon^p_{11} \)

Plastic strain component 11

Strain

statev[3]

\( \varepsilon^p_{22} \)

Plastic strain component 22

Strain

statev[4]

\( \varepsilon^p_{33} \)

Plastic strain component 33

Strain

statev[5]

\( \varepsilon^p_{12} \)

Plastic strain component 12 (engineering)

Strain

statev[6]

\( \varepsilon^p_{13} \)

Plastic strain component 13 (engineering)

Strain

statev[7]

\( \varepsilon^p_{23} \)

Plastic strain component 23 (engineering)

Strain

statev[8+6i]

\( X_{i,11} \)

Backstress i, component 11

Stress

statev[9+6i]

\( X_{i,22} \)

Backstress i, component 22

Stress

statev[10+6i]

\( X_{i,33} \)

Backstress i, component 33

Stress

statev[11+6i]

\( X_{i,12} \)

Backstress i, component 12

Stress

statev[12+6i]

\( X_{i,13} \)

Backstress i, component 13

Stress

statev[13+6i]

\( X_{i,23} \)

Backstress i, component 23

Stress

For \( N_{kin} = 2 \) backstresses: \( n_{statev} = 1 + 1 + 6 + 12 = 20 \)

// Example usage: 316 stainless steel with 2 backstresses, 1 isotropic hardening
int N_iso = 1;
int N_kin = 2;
vec props(6 + 2*N_iso + 2*N_kin);
props(0) = 200000;      // E = 200 GPa
props(1) = 0.3;         // nu = 0.3
props(2) = 1.7e-5;      // alpha = 17e-6 /K
props(3) = 150;         // sigma_Y = 150 MPa
props(4) = N_iso;       // 1 isotropic hardening term
props(5) = N_kin;       // 2 kinematic hardening terms

// Isotropic hardening (Voce)
props(6) = 100;         // Q1 = 100 MPa (saturation stress)
props(7) = 10;          // b1 = 10 (hardening rate)

// First backstress (long-range)
props(8) = 300000;      // C1 = 300 GPa
props(9) = 1000;        // D1 = 1000 (moderate recovery)

// Second backstress (short-range)
props(10) = 50000;      // C2 = 50 GPa
props(11) = 100;        // D2 = 100 (strong recovery)

vec statev = zeros(1 + 1 + 6 + 6*N_kin);  // 20 state variables
statev(0) = 20.0;  // Reference temperature 20°C

vec Etot = {0.002, -0.0006, -0.0006, 0.0, 0.0, 0.0};  // 0.2% axial strain
vec DEtot = {0.0001, -0.00003, -0.00003, 0.0, 0.0, 0.0};
vec sigma = zeros(6);
mat Lt = zeros(6,6);
mat L = zeros(6,6);
vec sigma_in = zeros(6);
mat DR = eye(3,3);

umat_plasticity_chaboche_CCP(Etot, DEtot, sigma, Lt, L, sigma_in, DR,
                             12, props, 20, statev, 20.0, 0.0, 0.0, 1.0,
                             Wm, Wm_r, Wm_ir, Wm_d, 3, 3, false, 0, tnew_dt);

// Check backstress components
vec X1 = statev.subvec(8, 13);   // First backstress
vec X2 = statev.subvec(14, 19);  // Second backstress

See also

Ireal() for real identity tensor (Voigt 6×6)

See also

Idev() for deviatoric projection tensor

See also

Fischer_Burmeister() for complementarity solver

See also

eta_stress() for shifted Stress Functions computation

See also

denom_FB_N_Mises() for CCP denominator

References:

  • Chaboche, J. L. (1986). “Time-independent constitutive theories for cyclic plasticity.” International Journal of Plasticity, 2(2), 149-188.

  • Chaboche, J. L. (1989). “Constitutive equations for cyclic plasticity and cyclic viscoplasticity.” International Journal of Plasticity, 5(3), 247-302.

  • Chaboche, J. L. (2008). “A review of some plasticity and viscoplasticity constitutive theories.” International Journal of Plasticity, 24(10), 1642-1693.

  • Armstrong, P. J., & Frederick, C. O. (1966). “A mathematical representation of the multiaxial Bauschinger effect.” CEGB Report RD/B/N731.

  • Lemaitre, J., & Chaboche, J. L. (1990). Mechanics of Solid Materials. Cambridge University Press.

  • Ortiz, M., & Simo, J. C. (1986). “An analysis of a new class of integration algorithms for elastoplastic constitutive relations.” International Journal for Numerical Methods in Engineering, 23(3), 353-366.

Note

The Chaboche model excels at capturing ratcheting and cyclic plasticity

Note

Multiple backstresses are essential for accurate multiaxial loading predictions

Note

Parameter identification typically requires cyclic test data (tension-compression, torsion)

Note

For monotonic loading, simpler isotropic hardening may suffice

Note

The model assumes small strains (< 10%); use updated Lagrangian for larger strains

Note

Dynamic recovery (D > 0) prevents unbounded backstress growth under cycling

Parameters:

  • Etot – Total strain tensor at beginning of increment (Voigt notation: 6×1 vector)

  • DEtot – Strain increment tensor (Voigt notation: 6×1 vector)

  • sigma – Stress tensor (Voigt notation: 6×1 vector) [output]

  • Lt – Consistent tangent modulus \( \mathbf{L}_t \) (6×6 matrix) [output]

  • L – Elastic stiffness tensor (6×6 matrix) [output]

  • sigma_in – Internal stress contribution for explicit solvers (6×1 vector) [output]

  • DR – Rotation increment matrix (3×3) for objective integration

  • nprops – Number of material properties

  • props – Material properties vector (see table above)

  • nstatev – Number of state variables

  • statev – State variables vector (see table above) [input/output]

  • T – Temperature at beginning of increment

  • DT – Temperature increment

  • Time – Time at beginning of increment

  • DTime – Time increment

  • Wm – Total mechanical work [output]

  • Wm_r – Recoverable (elastic) work [output]

  • Wm_ir – Irrecoverable work stored in plastic deformation [output]

  • Wm_d – Dissipated (plastic) work [output]

  • ndi – Number of direct stress components (typically 3)

  • nshr – Number of shear stress components (typically 3)

  • start – Flag indicating first increment (true) or continuation (false)

  • solver_type – Solver type: 0=implicit, 1=explicit, 2=dynamic implicit

  • tnew_dt – Suggested new time step size for adaptive time stepping [output]

void umat_plasticity_iso_CCP(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, arma::vec &sigma_in, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt)

Elastic-plastic constitutive model with isotropic hardening solved by the Convex Cutting Plane (CCP) algorithm.

This function implements an elastic-plastic material model for small and finite strain analysis. The model features:

  • J2 (von Mises) plasticity with associative flow rule

  • Isotropic hardening following a power law: \( H_p = k \cdot p^m \)

  • Convex Cutting Plane algorithm for return mapping

  • Thermal expansion effects

  • Consistent tangent modulus for implicit FE analysis

Constitutive Equations:

The yield function is defined as:

\[\Phi(\boldsymbol{\sigma}, p) = \sigma_{eq} - H_p(p) - \sigma_Y \leq 0 \]
where:

  • \( \sigma_{eq} = \sqrt{\frac{3}{2} \mathbf{s} : \mathbf{s}} \) is the von Mises equivalent stress

  • \( \mathbf{s} \) is the deviatoric stress tensor

  • \( H_p(p) = k \cdot p^m \) is the isotropic hardening function

  • \( p = \int \sqrt{\frac{2}{3} \dot{\boldsymbol{\varepsilon}}^p : \dot{\boldsymbol{\varepsilon}}^p} \, dt \) is the accumulated plastic strain

  • \( \sigma_Y \) is the initial yield stress

The plastic flow rule (associative plasticity):

\[\dot{\boldsymbol{\varepsilon}}^p = \dot{\lambda} \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} = \dot{\lambda} \frac{3}{2} \frac{\mathbf{s}}{\sigma_{eq}} \]

Evolution of accumulated plastic strain:

\[\dot{p} = \sqrt{\frac{2}{3} \dot{\boldsymbol{\varepsilon}}^p : \dot{\boldsymbol{\varepsilon}}^p} = \dot{\lambda} \]

Convex Cutting Plane Algorithm:

The CCP algorithm solves the return mapping problem by reformulating it as a complementarity problem:

  • Find \( \Delta p \geq 0 \) such that \( \Phi(\boldsymbol{\sigma}, p) \leq 0 \) and \( \Delta p \cdot \Phi = 0 \)

  • This is solved using the Fischer-Burmeister function for robust convergence

  • The method provides a consistent tangent modulus for implicit finite element analysis

Material Parameters (props):

The material properties vector must contain 6 constants:

Index

Symbol

Description

Units

props[0]

\( E \)

Young’s modulus

Stress

props[1]

\( \nu \)

Poisson’s ratio

-

props[2]

\( \alpha \)

Isotropic thermal expansion coefficient

1/Temperature

props[3]

\( \sigma_Y \)

Initial yield stress

Stress

props[4]

\( k \)

Hardening parameter

Stress

props[5]

\( m \)

Hardening exponent

-

State Variables (statev):

The state variables vector contains 8 internal variables:

Index

Symbol

Description

Units

statev[0]

\( T_{init} \)

Initial temperature

Temperature

statev[1]

\( p \)

Accumulated plastic strain

-

statev[2]

\( \varepsilon^p_{11} \)

Plastic strain component 11

-

statev[3]

\( \varepsilon^p_{22} \)

Plastic strain component 22

-

statev[4]

\( \varepsilon^p_{33} \)

Plastic strain component 33

-

statev[5]

\( \varepsilon^p_{12} \)

Plastic strain component 12 (×2 in Voigt)

-

statev[6]

\( \varepsilon^p_{13} \)

Plastic strain component 13 (×2 in Voigt)

-

statev[7]

\( \varepsilon^p_{23} \)

Plastic strain component 23 (×2 in Voigt)

- 

// Example usage:
vec Etot = {0.001, 0.0, 0.0, 0.0, 0.0, 0.0};
vec DEtot = {0.0001, 0.0, 0.0, 0.0, 0.0, 0.0};
vec sigma = zeros(6);
mat Lt = zeros(6,6);
mat L = zeros(6,6);
vec sigma_in = zeros(6);
mat DR = eye(3,3);
vec props = {70000, 0.3, 1e-5, 200, 500, 0.2};
vec statev = zeros(8);

umat_plasticity_iso_CCP(Etot, DEtot, sigma, Lt, L, sigma_in, DR, 6, props, 8, statev,
                        20.0, 0.0, 0.0, 1.0, Wm, Wm_r, Wm_ir, Wm_d,
                        3, 3, true, 0, tnew_dt);

See also

Fischer_Burmeister_m() for the complementarity solver

See also

L_iso() for isotropic elastic stiffness construction

See also

eta_stress() for plastic flow direction computation

Note

Voigt notation convention: [11, 22, 33, 12, 13, 23] with engineering shear strains (γ = 2ε)

Note

The consistent tangent modulus Lt ensures quadratic convergence in implicit Newton-Raphson schemes

Note

For explicit solvers (solver_type=1), use sigma_in instead of Lt

Parameters:

  • Etot – Total strain tensor at beginning of increment (Voigt notation: 6×1 vector)

  • DEtot – Strain increment tensor (Voigt notation: 6×1 vector)

  • sigma – Stress tensor (Voigt notation: 6×1 vector) [output]

  • Lt – Consistent tangent modulus \( \mathbf{L}_t = \frac{\partial \boldsymbol{\sigma}}{\partial \boldsymbol{\varepsilon}} \) (6×6 matrix) [output]

  • L – Elastic stiffness tensor (6×6 matrix) [output]

  • sigma_in – Internal stress contribution for explicit solvers (6×1 vector) [output]

  • DR – Rotation increment matrix (3×3) for objective integration

  • nprops – Number of material properties

  • props – Material properties vector (see table above)

  • nstatev – Number of state variables

  • statev – State variables vector (see table above) [input/output]

  • T – Temperature at beginning of increment

  • DT – Temperature increment

  • Time – Time at beginning of increment

  • DTime – Time increment

  • Wm – Total mechanical work [output]

  • Wm_r – Recoverable (elastic) work [output]

  • Wm_ir – Irrecoverable work stored in hardening [output]

  • Wm_d – Dissipated (plastic) work [output]

  • ndi – Number of direct stress components (typically 3)

  • nshr – Number of shear stress components (typically 3)

  • start – Flag indicating first increment (true) or continuation (false)

  • solver_type – Solver type: 0=implicit, 1=explicit, 2=dynamic implicit

  • tnew_dt – Suggested new time step size for adaptive time stepping [output]

void umat_plasticity_kin_iso_CCP(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, arma::vec &sigma_in, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt)

Elastic-plastic constitutive model with combined isotropic and linear kinematic hardening.

This function implements an elastic-plastic material model combining both isotropic and linear kinematic hardening mechanisms. The model features:

  • J2 (von Mises) plasticity with associative flow rule

  • Power law isotropic hardening: \( R(p) = k \cdot p^m \)

  • Linear (Prager) kinematic hardening: \( \dot{\mathbf{X}} = \frac{2}{3} H_{kin} \dot{\boldsymbol{\varepsilon}}^p \)

  • Convex Cutting Plane algorithm for return mapping

  • Consistent tangent modulus for implicit FE analysis

Constitutive Equations:

The yield function is defined as:

\[\Phi(\boldsymbol{\sigma}, \mathbf{X}, p) = \sigma_{eq}(\boldsymbol{\sigma} - \mathbf{X}) - R(p) - \sigma_Y \leq 0 \]
where:

  • \( \sigma_{eq}(\boldsymbol{\eta}) = \sqrt{\frac{3}{2} \boldsymbol{\eta}_{dev} : \boldsymbol{\eta}_{dev}} \) is the von Mises equivalent stress

  • \( \boldsymbol{\eta} = \boldsymbol{\sigma} - \mathbf{X} \) is the shifted (effective) stress

  • \( \mathbf{X} \) is the backstress (kinematic hardening variable)

  • \( R(p) = k \cdot p^m \) is the isotropic hardening function

  • \( \sigma_Y \) is the initial yield stress

Kinematic Hardening (Prager Linear):

The backstress evolves according to linear (Prager) kinematic hardening:

\[\dot{\mathbf{X}} = \frac{2}{3} H_{kin} \dot{\boldsymbol{\varepsilon}}^p = \frac{2}{3} H_{kin} \dot{p} \mathbf{n} \]
where:

  • \( H_{kin} \) is the kinematic hardening modulus

  • \( \mathbf{n} = \frac{3}{2} \frac{\boldsymbol{\eta}_{dev}}{\sigma_{eq}} \) is the flow direction

Combined Hardening Effect:

The combined hardening provides:

  • Isotropic hardening: Expands the yield surface uniformly (captures strain hardening)

  • Kinematic hardening: Translates the yield surface center (captures Bauschinger effect)

Plastic Flow Rule:

The plastic strain rate follows the associative flow rule:

\[\dot{\boldsymbol{\varepsilon}}^p = \dot{p} \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} = \dot{p} \mathbf{n} \]

Incremental Update:

For an increment \( \Delta p \):

\[\boldsymbol{\sigma}_{n+1} = \boldsymbol{\sigma}^{trial} - \Delta p \cdot \mathbf{L} : \mathbf{n} \]
\[\mathbf{X}_{n+1} = \mathbf{X}_n + \frac{2}{3} H_{kin} \Delta p \cdot \mathbf{n} \]

Consistent Tangent Modulus:

The algorithmic tangent accounts for both hardening mechanisms:

\[\mathbf{L}_t = \mathbf{L} - \frac{(\mathbf{L}:\mathbf{n}) \otimes (\mathbf{n}:\mathbf{L})}{\mathbf{n}:\mathbf{L}:\mathbf{n} + H_{iso} + H_{kin}} \]
where \( H_{iso} = \frac{dR}{dp} = k \cdot m \cdot p^{m-1} \) is the isotropic hardening modulus.

Material Parameters (props):

Index

Symbol

Description

Units

props[0]

\( E \)

Young’s modulus

Stress

props[1]

\( \nu \)

Poisson’s ratio

-

props[2]

\( \alpha \)

Thermal expansion coefficient

1/Temperature

props[3]

\( \sigma_Y \)

Initial yield stress

Stress

props[4]

\( k \)

Isotropic hardening coefficient

Stress

props[5]

\( m \)

Isotropic hardening exponent

-

props[6]

\( H_{kin} \)

Kinematic hardening modulus

Stress

State Variables (statev):

Total state variables required: \( n_{statev} = 14 \)

Index

Symbol

Description

Units

statev[0]

\( T_{init} \)

Initial temperature

Temperature

statev[1]

\( p \)

Accumulated plastic strain

Strain

statev[2]

\( \varepsilon^p_{11} \)

Plastic strain component 11

Strain

statev[3]

\( \varepsilon^p_{22} \)

Plastic strain component 22

Strain

statev[4]

\( \varepsilon^p_{33} \)

Plastic strain component 33

Strain

statev[5]

\( \varepsilon^p_{12} \)

Plastic strain component 12

Strain

statev[6]

\( \varepsilon^p_{13} \)

Plastic strain component 13

Strain

statev[7]

\( \varepsilon^p_{23} \)

Plastic strain component 23

Strain

statev[8]

\( X_{11} \)

Backstress component 11

Stress

statev[9]

\( X_{22} \)

Backstress component 22

Stress

statev[10]

\( X_{33} \)

Backstress component 33

Stress

statev[11]

\( X_{12} \)

Backstress component 12

Stress

statev[12]

\( X_{13} \)

Backstress component 13

Stress

statev[13]

\( X_{23} \)

Backstress component 23

Stress

References:

  • Prager, W. (1956). “A new method of analyzing stresses and strains in work-hardening plastic solids.” J. Appl. Mech., 23, 493-496.

  • Ziegler, H. (1959). “A modification of Prager’s hardening rule.” Q. Appl. Math., 17, 55-65.

  • Lemaitre, J., & Chaboche, J. L. (1990). Mechanics of Solid Materials. Cambridge University Press.

  • Simo, J. C., & Hughes, T. J. R. (1998). Computational Inelasticity. Springer.

See also

umat_plasticity_iso_CCP() for pure isotropic hardening

See also

umat_plasticity_chaboche_CCP() for nonlinear kinematic hardening

Note

Linear kinematic hardening cannot capture cyclic softening/hardening saturation

Note

For saturating kinematic hardening, use the Chaboche model (Armstrong-Frederick)

Note

The Bauschinger effect is proportional to the kinematic hardening modulus

Note

Setting \( H_{kin} = 0 \) recovers pure isotropic hardening

Note

Setting \( k = 0 \) (or \( m = 0 \)) recovers pure kinematic hardening

Parameters:

  • Etot – Total strain tensor at beginning of increment (Voigt notation: 6×1 vector)

  • DEtot – Strain increment tensor (Voigt notation: 6×1 vector)

  • sigma – Stress tensor (Voigt notation: 6×1 vector) [output]

  • Lt – Consistent tangent modulus (6×6 matrix) [output]

  • L – Elastic stiffness tensor (6×6 matrix) [output]

  • sigma_in – Internal stress contribution for explicit solvers (6×1 vector) [output]

  • DR – Rotation increment matrix (3×3) for objective integration

  • nprops – Number of material properties

  • props – Material properties vector (see table above)

  • nstatev – Number of state variables

  • statev – State variables vector (see table above) [input/output]

  • T – Temperature at beginning of increment

  • DT – Temperature increment

  • Time – Time at beginning of increment

  • DTime – Time increment

  • Wm – Total mechanical work [output]

  • Wm_r – Recoverable (elastic) work [output]

  • Wm_ir – Irrecoverable work stored in hardening [output]

  • Wm_d – Dissipated (plastic) work [output]

  • ndi – Number of direct stress components (typically 3)

  • nshr – Number of shear stress components (typically 3)

  • start – Flag indicating first increment (true) or continuation (false)

  • solver_type – Solver type: 0=implicit, 1=explicit, 2=dynamic implicit

  • tnew_dt – Suggested new time step size for adaptive time stepping [output]

void umat_sma_mono(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt)

Micromechanical monocrystal model for SMA based on Patoor et al. (1996)

This function implements the micromechanical monocrystal model of Patoor et al. (1996) for shape memory alloys. The model describes the martensitic phase transformation in a single crystal by tracking N individual martensite variants, where N is typically 12 or 24 depending on the crystallographic system.

Key Features:

  • Explicit tracking of N martensite variant volume fractions \( f_n \) (n = 1, …, N)

  • Crystallographic transformation strains for each variant

  • Variant selection driven by resolved thermodynamic driving force

  • Inter-variant interactions through hardening matrix

  • Applicable to single crystal SMA behavior

Physical Background:

In single crystal SMAs, the austenite-to-martensite transformation occurs through the formation of distinct crystallographic variants. For cubic-to-monoclinic transformations (e.g., NiTi), there are typically 12 or 24 habit plane variants, each characterized by:

  • A specific transformation strain tensor \( \boldsymbol{\varepsilon}^{tr}_n \)

  • A habit plane normal and transformation direction

Variant Volume Fractions:

The microstructural state is characterized by N variant volume fractions:

\[f_n \geq 0, \quad \sum_{n=1}^{N} f_n \leq 1, \quad f_A = 1 - \sum_{n=1}^{N} f_n \]
where:

  • \( f_n \) is the volume fraction of martensite variant n

  • \( f_A \) is the austenite volume fraction

  • The total martensite fraction is \( \xi = \sum_{n=1}^{N} f_n \)

Transformation Strain:

The macroscopic transformation strain is the volume-weighted sum over all variants:

\[\boldsymbol{\varepsilon}^{tr} = \sum_{n=1}^{N} f_n \boldsymbol{\varepsilon}^{tr}_n \]
where \( \boldsymbol{\varepsilon}^{tr}_n \) is the crystallographic transformation strain tensor of variant n, computed from the lattice correspondence.

Thermodynamic Driving Force:

For each variant n, the driving force for transformation is:

\[F_n = \boldsymbol{\sigma} : \boldsymbol{\varepsilon}^{tr}_n - \Delta G^{chem}(T) - \sum_{m=1}^{N} H_{nm} f_m \]
where:

  • \( \boldsymbol{\sigma} : \boldsymbol{\varepsilon}^{tr}_n \) is the mechanical driving force

  • \( \Delta G^{chem}(T) \) is the chemical free energy difference (temperature-dependent)

  • \( H_{nm} \) is the interaction matrix describing variant-variant hardening

Transformation Criteria:

Forward Transformation (A → M_n):

\[\Phi_n^{fwd} = F_n - F_c^{fwd} \leq 0, \quad \dot{f}_n \geq 0 \]

Reverse Transformation (M_n → A):

\[\Phi_n^{rev} = -F_n - F_c^{rev} \leq 0, \quad \dot{f}_n \leq 0 \]

where \( F_c^{fwd} \) and \( F_c^{rev} \) are critical driving forces for forward and reverse transformations.

Interaction Matrix:

The hardening matrix \( H_{nm} \) captures:

  • Self-hardening ( \( H_{nn} \)): resistance to growth of variant n

  • Latent hardening ( \( H_{nm}, n \neq m \)): interaction between different variants

Variant Selection:

Under applied stress, variants with favorable orientation (high resolved stress on transformation system) are preferentially activated. This leads to:

  • Single variant formation under uniaxial loading along specific orientations

  • Multi-variant microstructures under complex loading

  • Texture-dependent macroscopic response

Number of Variants:

Common crystallographic systems:

  • Cubic → Orthorhombic: N = 6 variants

  • Cubic → Monoclinic (NiTi): N = 12 or 24 variants

  • Cubic → Tetragonal: N = 3 variants

State Variables (statev):

Total state variables required: \( n_{statev} = 1 + N + 6 \) (for N variants)

Index

Symbol

Description

Units

statev[0]

\( T_{init} \)

Initial/reference temperature

Temperature

statev[1]

\( f_1 \)

Volume fraction of martensite variant 1

-

statev[2]

\( f_2 \)

Volume fraction of martensite variant 2

-

…

…

…

…

statev[N]

\( f_N \)

Volume fraction of martensite variant N

-

statev[N+1]

\( \varepsilon^{tr}_{11} \)

Macroscopic transformation strain component 11

Strain

statev[N+2]

\( \varepsilon^{tr}_{22} \)

Macroscopic transformation strain component 22

Strain

statev[N+3]

\( \varepsilon^{tr}_{33} \)

Macroscopic transformation strain component 33

Strain

statev[N+4]

\( \varepsilon^{tr}_{12} \)

Macroscopic transformation strain component 12

Strain

statev[N+5]

\( \varepsilon^{tr}_{13} \)

Macroscopic transformation strain component 13

Strain

statev[N+6]

\( \varepsilon^{tr}_{23} \)

Macroscopic transformation strain component 23

Strain

The macroscopic transformation strain is computed as: \( \boldsymbol{\varepsilon}^{tr} = \sum_{n=1}^{N} f_n \boldsymbol{\varepsilon}^{tr}_n \)

References:

  • Patoor, E., Eberhardt, A., & Berveiller, M. (1996). “Micromechanical modelling of

    superelasticity in shape memory alloys.”

    Journal de Physique IV, 6(C1), 277-292.

  • Patoor, E., Lagoudas, D. C., Entchev, P. B., Brinson, L. C., & Gao, X. (2006). “Shape memory alloys, Part I: General properties and modeling of single crystals.” Mechanics of Materials, 38(5-6), 391-429.

  • Gall, K., & Sehitoglu, H. (1999). “The role of texture in tension-compression

    asymmetry in polycrystalline NiTi.”

    International Journal of Plasticity, 15(1), 69-92.

Note

This model is specifically for single crystal SMA behavior

Note

For polycrystalline SMAs, use this model within a homogenization scheme

Note

The number of variants N depends on the crystallographic transformation system

Note

Variant transformation strains must be provided based on crystallographic data

Parameters:

  • Etot – Total strain tensor at beginning of increment (Voigt notation: 6×1 vector)

  • DEtot – Strain increment tensor (Voigt notation: 6×1 vector)

  • sigma – Cauchy stress tensor (Voigt notation: 6×1 vector) [output]

  • Lt – Consistent tangent modulus (6×6 matrix) [output]

  • L – Elastic stiffness tensor (6×6 matrix) [output]

  • DR – Rotation increment matrix (3×3) for objective integration

  • nprops – Number of material properties

  • props – Material properties vector

  • nstatev – Number of state variables (includes N variant volume fractions)

  • statev – State variables vector containing variant fractions [input/output]

  • T – Temperature at beginning of increment

  • DT – Temperature increment

  • Time – Time at beginning of increment

  • DTime – Time increment

  • Wm – Total mechanical work [output]

  • Wm_r – Recoverable (elastic) work [output]

  • Wm_ir – Irrecoverable work stored in transformation [output]

  • Wm_d – Dissipated work (hysteresis) [output]

  • ndi – Number of direct stress components (typically 3)

  • nshr – Number of shear stress components (typically 3)

  • start – Flag indicating first increment (true) or continuation (false)

  • tnew_dt – Suggested new time step size for adaptive time stepping [output]

void umat_prony_Nfast(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt)

Linear viscoelastic constitutive model with Prony series representation (Generalized Maxwell model)

This function implements a linear viscoelastic material model using the Prony series decomposition, also known as the Generalized Maxwell model. The model consists of multiple Maxwell elements in parallel, each characterized by a relaxation time and modulus.

Constitutive Equations:

The relaxation modulus is expressed as a Prony series:

\[E(t) = E_\infty + \sum_{i=1}^N E_i e^{-t/\tau_i} \]
where:

  • \( E_\infty \) is the long-term (equilibrium) modulus

  • \( E_i \) is the modulus of the i-th Maxwell element

  • \( \tau_i \) is the relaxation time of the i-th element

  • \( N \) is the number of Maxwell elements

Incremental Form:

The stress at time \( t + \Delta t \) is computed as:

\[\boldsymbol{\sigma}(t + \Delta t) = \mathbf{L}_\infty : \boldsymbol{\varepsilon}(t + \Delta t) + \sum_{i=1}^N \mathbf{q}_i(t + \Delta t) \]

where \( \mathbf{q}_i \) are the internal state variables (stress-like quantities) that evolve according to:

\[\mathbf{q}_i(t + \Delta t) = e^{-\Delta t/\tau_i} \mathbf{q}_i(t) + \frac{E_i}{E_\infty} \mathbf{L}_\infty : \left( e^{-\Delta t/\tau_i} - 1 \right) \Delta \boldsymbol{\varepsilon} \]

Bulk and Shear Decomposition:

The model handles volumetric and deviatoric responses independently:

  • Volumetric: \( K(t) = K_\infty + \sum_{i=1}^N K_i e^{-t/\tau_i^K} \)

  • Deviatoric: \( G(t) = G_\infty + \sum_{i=1}^N G_i e^{-t/\tau_i^G} \)

Consistent Tangent Modulus:

For implicit finite element analysis:

\[\mathbf{L}_t = \mathbf{L}_\infty + \sum_{i=1}^N \left( 1 - e^{-\Delta t/\tau_i} \right) \mathbf{L}_i \]

Material Parameters (props):

For N Maxwell elements, the props vector contains:

Index

Symbol

Description

Units

props[0]

\( N \)

Number of Maxwell elements

-

props[1]

\( E_\infty \)

Equilibrium Young’s modulus

Stress

props[2]

\( \nu_\infty \)

Equilibrium Poisson’s ratio

-

props[3]

\( \alpha \)

Thermal expansion coefficient

1/Temperature

props[4]

\( E_1 \)

Modulus of 1st Maxwell element

Stress

props[5]

\( \tau_1 \)

Relaxation time of 1st element

Time

…

…

…

…

props[4+2(i-1)]

\( E_i \)

Modulus of i-th element

Stress

props[5+2(i-1)]

\( \tau_i \)

Relaxation time of i-th element

Time

State Variables (statev):

For each Maxwell element, 6 internal state variables store the stress-like quantities \( \mathbf{q}_i \):

Index

Symbol

Description

Units

statev[0]

\( T_{init} \)

Initial temperature

Temperature

statev[1:6]

\( \mathbf{q}_1 \)

Internal stresses of 1st element (Voigt)

Stress

statev[7:12]

\( \mathbf{q}_2 \)

Internal stresses of 2nd element (Voigt)

Stress

…

…

…

…

statev[1+6(i-1):6i]

\( \mathbf{q}_i \)

Internal stresses of i-th element (Voigt)

Stress

Total state variables required: \( n_{statev} = 1 + 6N \)

// Example usage: 2 Maxwell elements
int N = 2;
vec props(8);
props(0) = N;          // 2 Maxwell elements
props(1) = 1000;       // E_inf = 1000 MPa
props(2) = 0.3;        // nu_inf = 0.3
props(3) = 1e-5;       // alpha = 1e-5 /K
props(4) = 500;        // E_1 = 500 MPa
props(5) = 0.1;        // tau_1 = 0.1 s
props(6) = 200;        // E_2 = 200 MPa
props(7) = 1.0;        // tau_2 = 1.0 s

vec statev = zeros(1 + 6*N);  // 1 + 12 state variables
vec Etot = {0.001, 0, 0, 0, 0, 0};
vec DEtot = {0.0001, 0, 0, 0, 0, 0};
vec sigma = zeros(6);
mat Lt = zeros(6,6);
mat DR = eye(3,3);

umat_prony_Nfast(Etot, DEtot, sigma, Lt, DR, 8, props, 13, statev,
                 20.0, 0.0, 0.0, 0.01, Wm, Wm_r, Wm_ir, Wm_d,
                 3, 3, false, tnew_dt);

See also

L_iso() for isotropic elastic stiffness construction

See also

Ivol() for volumetric projection tensor

See also

Idev() for deviatoric projection tensor

References:

  • Ferry, J. D. (1980). Viscoelastic Properties of Polymers (3rd ed.). Wiley.

  • Simo, J. C., & Hughes, T. J. R. (1998). Computational Inelasticity. Springer.

  • Park, S. W., & Schapery, R. A. (1999). “Methods of interconversion between linear viscoelastic material functions. Part I—A numerical method based on Prony series.” International Journal of Solids and Structures, 36(11), 1653-1675.

Note

The Prony series provides an efficient representation for viscoelastic relaxation

Note

Relaxation times should span the expected loading time scales

Note

For frequency-domain data, use fitting algorithms to extract Prony parameters

Note

The model assumes small strains and linear viscoelasticity

Note

Time step should be small relative to shortest relaxation time for accuracy

Parameters:

  • Etot – Total strain tensor at beginning of increment (Voigt notation: 6×1 vector)

  • DEtot – Strain increment tensor (Voigt notation: 6×1 vector)

  • sigma – Stress tensor (Voigt notation: 6×1 vector) [output]

  • Lt – Consistent tangent modulus (6×6 matrix) [output]

  • DR – Rotation increment matrix (3×3) for objective integration

  • nprops – Number of material properties

  • props – Material properties vector (see table above)

  • nstatev – Number of state variables

  • statev – State variables vector (see table above) [input/output]

  • T – Temperature at beginning of increment

  • DT – Temperature increment

  • Time – Time at beginning of increment

  • DTime – Time increment

  • Wm – Total mechanical work [output]

  • Wm_r – Recoverable (elastic) work [output]

  • Wm_ir – Irrecoverable work stored in viscous elements [output]

  • Wm_d – Dissipated (viscous) work [output]

  • ndi – Number of direct stress components (typically 3)

  • nshr – Number of shear stress components (typically 3)

  • start – Flag indicating first increment (true) or continuation (false)

  • tnew_dt – Suggested new time step size for adaptive time stepping [output]

void umat_zener_fast(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt)

Linear viscoelastic constitutive model using the Zener (Standard Linear Solid) model.

This function implements a linear viscoelastic material model using the Zener model, also known as the Standard Linear Solid (SLS). The Zener model consists of a spring in series with a Kelvin-Voigt element (spring and dashpot in parallel), providing both instantaneous elastic response and time-dependent relaxation.

Rheological Representation:

  E_0 (instantaneous)
─────────┬─────────
         │
  ┌──────┴──────┐
  │   E_1      η_1   │  (Kelvin-Voigt element)
  │ ──/\/\──┬──[=]──│
  │         │        │
  └─────────┴────────┘

Constitutive Equations:

The Zener model stress-strain relationship is given by:

\[ \boldsymbol{\sigma} + \tau \dot{\boldsymbol{\sigma}} = E_0 \boldsymbol{\varepsilon} + (E_0 + E_1) \tau \dot{\boldsymbol{\varepsilon}} \]
where:

  • \( E_0 \) is the instantaneous (glassy) modulus

  • \( E_1 \) is the modulus of the Kelvin-Voigt element

  • \( \tau = \eta_1 / E_1 \) is the relaxation time

  • \( \eta_1 \) is the viscosity of the dashpot

Relaxation Modulus:

The relaxation modulus for the Zener model is:

\[ E(t) = E_\infty + (E_0 - E_\infty) e^{-t/\tau} \]
where:

  • \( E_\infty = \frac{E_0 E_1}{E_0 + E_1} \) is the long-term (rubbery) modulus

  • \( E_0 \) is the instantaneous modulus

  • \( \tau \) is the relaxation time

Creep Compliance:

The creep compliance is:

\[ J(t) = \frac{1}{E_0} + \frac{1}{E_1} \left( 1 - e^{-t/\tau_c} \right) \]
where \( \tau_c = \eta_1 / E_1 \) is the retardation time.

Incremental Form:

For numerical integration over a time step \( \Delta t \):

\[ \boldsymbol{\sigma}_{n+1} = \mathbf{L}_\infty : \boldsymbol{\varepsilon}_{n+1} + \mathbf{q}_{n+1} \]
where the internal variable \( \mathbf{q} \) evolves as:
\[ \mathbf{q}_{n+1} = e^{-\Delta t/\tau} \mathbf{q}_n + \frac{E_1 - E_\infty}{E_\infty} \mathbf{L}_\infty : \left( e^{-\Delta t/\tau} - 1 \right) \Delta \boldsymbol{\varepsilon} \]

Material Parameters (props):

Index

Symbol

Description

Units

props[0]

\( E_\infty \)

Equilibrium (long-term) Young’s modulus

Stress

props[1]

\( \nu \)

Poisson’s ratio (assumed constant)

-

props[2]

\( \alpha \)

Thermal expansion coefficient

1/Temperature

props[3]

\( E_1 \)

Modulus of Kelvin-Voigt element

Stress

props[4]

\( \tau \)

Relaxation time

Time

State Variables (statev):

Total state variables required: \( n_{statev} = 7 \)

Index

Symbol

Description

Units

statev[0]

\( T_{init} \)

Initial temperature

Temperature

statev[1:6]

\( \mathbf{q} \)

Internal stress-like variable (Voigt 6×1)

Stress

References:

  • Zener, C. (1948). Elasticity and Anelasticity of Metals. University of Chicago Press.

  • Ferry, J. D. (1980). Viscoelastic Properties of Polymers (3rd ed.). Wiley.

  • Simo, J. C., & Hughes, T. J. R. (1998). Computational Inelasticity. Springer.

Note

The Zener model is the simplest model capturing both creep and relaxation

Note

For multiple relaxation times, use the Prony series model (umat_prony_Nfast)

Note

Time step should be small relative to relaxation time for accuracy

Note

The model assumes small strains and linear viscoelasticity

Parameters:

  • Etot – Total strain tensor at beginning of increment (Voigt notation: 6×1 vector)

  • DEtot – Strain increment tensor (Voigt notation: 6×1 vector)

  • sigma – Stress tensor (Voigt notation: 6×1 vector) [output]

  • Lt – Consistent tangent modulus (6×6 matrix) [output]

  • DR – Rotation increment matrix (3×3) for objective integration

  • nprops – Number of material properties

  • props – Material properties vector (see table above)

  • nstatev – Number of state variables

  • statev – State variables vector (see table above) [input/output]

  • T – Temperature at beginning of increment

  • DT – Temperature increment

  • Time – Time at beginning of increment

  • DTime – Time increment

  • Wm – Total mechanical work [output]

  • Wm_r – Recoverable (elastic) work [output]

  • Wm_ir – Irrecoverable work stored in viscous elements [output]

  • Wm_d – Dissipated (viscous) work [output]

  • ndi – Number of direct stress components (typically 3)

  • nshr – Number of shear stress components (typically 3)

  • start – Flag indicating first increment (true) or continuation (false)

  • tnew_dt – Suggested new time step size for adaptive time stepping [output]

void umat_zener_Nfast(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt)

Generalized Zener (Standard Linear Solid) viscoelastic model with N parallel branches.

This function implements a generalized linear viscoelastic material model using N Zener (Maxwell) elements in parallel with a long-term equilibrium spring. This is equivalent to a Generalized Maxwell model and provides a discrete approximation of the continuous relaxation spectrum.

Rheological Representation:

                 E_∞ (equilibrium spring)
────────────────────/\/\/\────────────────
         │                       │
         ├──── E_1 ─[η_1]────────┤   Branch 1
         │                       │
         ├──── E_2 ─[η_2]────────┤   Branch 2
         │                       │
         │         ...           │
         │                       │
         └──── E_N ─[η_N]────────┘   Branch N

Constitutive Equations:

The relaxation modulus is expressed as a sum of exponentials (Prony series):

\[ E(t) = E_\infty + \sum_{i=1}^N E_i e^{-t/\tau_i} \]
where:

  • \( E_\infty \) is the long-term (equilibrium) modulus

  • \( E_i \) is the modulus of the i-th Maxwell element

  • \( \tau_i = \eta_i / E_i \) is the relaxation time of the i-th element

  • \( N \) is the number of Maxwell elements

Instantaneous Modulus:

At \( t = 0 \):

\[ E_0 = E(0) = E_\infty + \sum_{i=1}^N E_i \]

Stress Decomposition:

The total stress is the sum of the equilibrium stress and internal stresses:

\[ \boldsymbol{\sigma} = \mathbf{L}_\infty : \boldsymbol{\varepsilon} + \sum_{i=1}^N \mathbf{q}_i \]
where \( \mathbf{q}_i \) are the internal stress-like variables for each Maxwell element.

Internal Variable Evolution:

Each internal variable evolves according to:

\[ \dot{\mathbf{q}}_i + \frac{1}{\tau_i} \mathbf{q}_i = \frac{E_i}{E_\infty} \mathbf{L}_\infty : \dot{\boldsymbol{\varepsilon}} \]

Incremental Form:

For numerical integration over a time step \( \Delta t \):

\[ \mathbf{q}_{i,n+1} = e^{-\Delta t/\tau_i} \mathbf{q}_{i,n} + \frac{E_i}{E_\infty} \mathbf{L}_\infty : \frac{\tau_i}{\Delta t} \left( 1 - e^{-\Delta t/\tau_i} \right) \Delta \boldsymbol{\varepsilon} \]

Consistent Tangent Modulus:

For implicit finite element analysis:

\[ \mathbf{L}_t = \mathbf{L}_\infty \left( 1 + \sum_{i=1}^N \frac{E_i}{E_\infty} \frac{\tau_i}{\Delta t} \left( 1 - e^{-\Delta t/\tau_i} \right) \right) \]

Material Parameters (props):

For N Maxwell elements, the props vector contains:

Index

Symbol

Description

Units

props[0]

\( N \)

Number of Maxwell elements

-

props[1]

\( E_\infty \)

Equilibrium Young’s modulus

Stress

props[2]

\( \nu \)

Poisson’s ratio (assumed constant)

-

props[3]

\( \alpha \)

Thermal expansion coefficient

1/Temperature

props[4]

\( E_1 \)

Modulus of 1st Maxwell element

Stress

props[5]

\( \tau_1 \)

Relaxation time of 1st element

Time

…

…

…

…

props[4+2(i-1)]

\( E_i \)

Modulus of i-th element

Stress

props[5+2(i-1)]

\( \tau_i \)

Relaxation time of i-th element

Time

State Variables (statev):

For each Maxwell element, 6 internal state variables store the stress-like quantities:

Index

Symbol

Description

Units

statev[0]

\( T_{init} \)

Initial temperature

Temperature

statev[1:6]

\( \mathbf{q}_1 \)

Internal stresses of 1st element (Voigt)

Stress

statev[7:12]

\( \mathbf{q}_2 \)

Internal stresses of 2nd element (Voigt)

Stress

…

…

…

…

statev[1+6(i-1):6i]

\( \mathbf{q}_i \)

Internal stresses of i-th element (Voigt)

Stress

Total state variables required: \( n_{statev} = 1 + 6N \)

References:

  • Zener, C. (1948). Elasticity and Anelasticity of Metals. University of Chicago Press.

  • Ferry, J. D. (1980). Viscoelastic Properties of Polymers (3rd ed.). Wiley.

  • Simo, J. C., & Hughes, T. J. R. (1998). Computational Inelasticity. Springer.

  • Park, S. W., & Schapery, R. A. (1999). “Methods of interconversion between linear viscoelastic material functions.” Int. J. Solids Struct., 36(11), 1653-1675.

See also

umat_zener_fast() for single Maxwell element version

See also

umat_prony_Nfast() for equivalent Prony series formulation

Note

Relaxation times should span the expected loading time scales (decades in log-time)

Note

Typically 5-10 Maxwell elements are sufficient for most polymers

Note

For frequency-domain data, use fitting algorithms to extract Prony parameters

Note

The model assumes small strains and linear viscoelasticity

Note

Time step should be small relative to shortest relaxation time for accuracy

Parameters:

  • Etot – Total strain tensor at beginning of increment (Voigt notation: 6×1 vector)

  • DEtot – Strain increment tensor (Voigt notation: 6×1 vector)

  • sigma – Stress tensor (Voigt notation: 6×1 vector) [output]

  • Lt – Consistent tangent modulus (6×6 matrix) [output]

  • DR – Rotation increment matrix (3×3) for objective integration

  • nprops – Number of material properties

  • props – Material properties vector (see table above)

  • nstatev – Number of state variables

  • statev – State variables vector (see table above) [input/output]

  • T – Temperature at beginning of increment

  • DT – Temperature increment

  • Time – Time at beginning of increment

  • DTime – Time increment

  • Wm – Total mechanical work [output]

  • Wm_r – Recoverable (elastic) work [output]

  • Wm_ir – Irrecoverable work stored in viscous elements [output]

  • Wm_d – Dissipated (viscous) work [output]

  • ndi – Number of direct stress components (typically 3)

  • nshr – Number of shear stress components (typically 3)

  • start – Flag indicating first increment (true) or continuation (false)

  • tnew_dt – Suggested new time step size for adaptive time stepping [output]

void get_L_elastic(phase_characteristics&)