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 \(\times\) 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 \times 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. STATEV layout: Wm (4), then, when ndi < 3, the 6 - ntens stress-free strain components the material point solves (condensed_state.hpp, simcoon order, rotated with DROT and returned in Etot), then the law’s state variables (nstatev - 4 - (6 - ntens)). 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 \times 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. STATEV layout: Wm (4), Wt (3), then the stress-free strain components when ndi < 3 (as abaqus2smart_M), then the law’s state variables.

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, const bool &nlgeom = true, const arma::vec &Etot_end = arma::vec())

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 (dsigma/deps of the corotated kernel), condensed on the element’s components when ndi < 3 (select_umat_M): ddsdde is its block on them
nlgeom – whether the step is geometrically nonlinear (abaqus_nlgeom(dfgrd1)): ddsdde then gets abaqus_jacobian of Lt for the 3D, plane-strain and axisymmetric cases; plane stress and 1D take the condensed tangent as it is
Etot_end – total strain at the end of the increment (Etot + DEtot): when ndi < 3 its stress-free components are stored in STATEV after Wm (layout in abaqus2smart_M); empty leaves those slots untouched

arma::mat abaqus_jacobian(const arma::mat &Lt, const arma::vec &sigma)

The Abaqus material Jacobian of a corotated kernel tangent.

Abaqus defines DDSDDE through the Jaumann rate of the Kirchhoff stress divided by J: \( \frac{1}{J}\frac{\partial (J\boldsymbol{\sigma})}{\partial \boldsymbol{\varepsilon}} = \frac{\partial \boldsymbol{\sigma}}{\partial \boldsymbol{\varepsilon}} + \boldsymbol{\sigma} \otimes \mathbf{I} \). The kernels return the first term in Abaqus’s corotated frame; this adds the symmetric part of the second (the default symmetric solver keeps that part). Used by smart2abaqus_M and smart2abaqus_T for NLGEOM steps. Without NLGEOM the Jacobian is the kernel tangent (J = 1 by assumption), which also keeps linear-perturbation procedures (*FREQUENCY, *BUCKLE) on the right stiffness.

Parameters:

Lt – kernel tangent (6x6, Voigt, engineering shears)
sigma – Cauchy stress at the end of the increment (6)

Returns:

DDSDDE as a 6x6 matrix

bool abaqus_nlgeom(const double *dfgrd1)

Whether Abaqus runs the step with NLGEOM: DFGRD1 is the identity otherwise.

Parameters:: dfgrd1 – the DFGRD1 array of the UMAT interface (9 values)

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, const bool &nlgeom = true, const arma::vec &Etot_end = arma::vec())

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 (condensed when ndi < 3, select_umat_T)
dSdT – simcoon stress-temperature tangent matrix (condensed when ndi < 3)
drpldE – simcoon heat generation-strain tangent matrix (condensed when ndi < 3)
drpldT – simcoon heat generation-temperature tangent matrix
nlgeom – as smart2abaqus_M
Etot_end – as smart2abaqus_M, stored after Wm and Wt (layout in abaqus2smart_T)

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 std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &stress, 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, const int &tangent_mode = tangent_default)

Ladevèze-Le Dantec (LLD) damage model of the unidirectional ply, with coupled transverse-shear plasticity in effective stress (LLDM0).

Matrix-dominated degradation of a unidirectional fibre-reinforced ply: two scalar damage variables, \( d_{22} \) (transverse) and \( d_{12} \) (shear), and a plasticity of the matrix written in the effective (undamaged) stress. No fibre failure. The equations below are in the fibre frame, fibres along 1; props(0) puts the fibres along the global axis 1, 2 or 3, the law then runs in that frame and returns its stress, tangent and plastic strain in the global one.

Damaged elasticity. Transversely isotropic compliance \( \mathbf{S}_0 \) (moduli \( E_L, E_T, \nu_{TL}, \nu_{TT}, G_{LT} \), \( G_{TT} = E_T/(2(1+\nu_{TT})) \)), damaged as

\[ \tilde{\mathbf{S}}(d_{22}, d_{12}):\quad S_{22}, S_{33}, S_{23}, S_{66} \leftarrow \frac{S^0}{1-d_{22}}, \qquad S_{44}, S_{55} \leftarrow \frac{S^0}{1-d_{12}}, \]

the fibre terms \( S_{11}, S_{12}, S_{13} \) unchanged, and \( \boldsymbol{\sigma} = \tilde{\mathbf{L}} : (\boldsymbol{\varepsilon} - \boldsymbol{\alpha}\,\Delta\theta - \boldsymbol{\varepsilon}^p) \), \( \tilde{\mathbf{L}} = \tilde{\mathbf{S}}^{-1} \). \( \nu_{TL} \) is the MINOR Poisson ratio ( \( \nu_{TL} = \nu_{LT} E_T/E_L \), about 0.02 for carbon/epoxy), the convention of L_isotrans().

Plasticity in effective stress. With \( \tilde{\boldsymbol{\sigma}} = \mathbf{L}_0 : (\boldsymbol{\varepsilon} - \boldsymbol{\alpha}\,\Delta\theta - \boldsymbol{\varepsilon}^p) \) and \( \boldsymbol{\Theta} = \mathrm{diag}(0, A_{ts}, A_{ts}, 1, 1, 0) \) (transverse and fibre-plane shear stresses only),

\[ f = \sigma^{\mathrm{Mises}}\!\left(\boldsymbol{\Theta}\tilde{\boldsymbol{\sigma}}\right) - R(p) - \sigma_{ts,0} \leq 0, \qquad R(p) = \beta_{ts}\, p^{m_{ts}}, \]

with the associated flow \( \dot{\boldsymbol{\varepsilon}}^p = \dot p\, \partial f/\partial\tilde{\boldsymbol{\sigma}} = \dot p\, \boldsymbol{\Theta}\,\boldsymbol{\eta}(\boldsymbol{\Theta}\tilde{\boldsymbol{\sigma}}) \) ( \( \boldsymbol{\eta} \) = eta_stress()). The effective stress does not depend on the damage: the plasticity is solved first (convex cutting plane, Fischer-Burmeister), the damage then at fixed plastic strain.

Damage. The driving forces of the criterion are unilateral in the transverse direction,

\[ Y_{22} = \frac{\langle\sigma_{22}\rangle_+^2 + \langle\sigma_{33}\rangle_+^2}{2E_T(1-d_{22})^2}, \quad Y_{12} = \frac{\sigma_{12}^2 + \sigma_{13}^2}{2G_{LT}(1-d_{12})^2}, \quad Y_{ts} = \sqrt{Y_{12} + b\,Y_{22}}, \quad Y_t = \sqrt{Y_{22}}, \]

and both damage variables follow \( Y_{ts} \) (Kuhn-Tucker, \( \dot d \geq 0 \)):

\[ d_{22} + \lambda(d_{22}) = \frac{\langle Y_{ts} - Y_{22,0}\rangle_+}{Y_{22,c}}, \qquad d_{12} + \lambda(d_{12}) = \frac{\langle Y_{ts} - Y_{12,0}\rangle_+}{Y_{12,c}}, \]

\( \lambda \) the penalty of lagrange_pow_1() keeping \( d < 1 \). The two are coupled through \( Y_{ts} \) and through the stress they both relax (the damage directions \( \boldsymbol{\Lambda}_j = \partial\tilde{\mathbf{S}}/\partial d_j : \boldsymbol{\sigma} \)): the 2x2 system is full and non-symmetric. Transverse failure: when \( Y_t > Y_{22,u} \), \( d_{22} \) jumps to 0.99 and stays there; \( d_{12} \) keeps its own law.

Tangent. The tangent of the two-stage update: \( \mathrm{d}\boldsymbol{\varepsilon}^e/\mathrm{d}\boldsymbol{\varepsilon} \) from the plasticity, the flow direction held fixed over the increment (continuum: its error is first order in the increment), then the exact linearization of the active damage conditions, \( \mathbf{L}_t = \tilde{\mathbf{L}}\,\frac{\mathrm{d}\boldsymbol{\varepsilon}^e}{\mathrm{d}\boldsymbol{\varepsilon}} - \sum_j \tilde{\mathbf{L}}\boldsymbol{\Lambda}_j \otimes \frac{\mathrm{d}d_j}{\mathrm{d}\boldsymbol{\varepsilon}} \). tangent_mode is not used. L returns the damaged stiffness \( \tilde{\mathbf{L}} \).

Work split. \( W_m^{ir} \) stores the hardening energy ( \( \tfrac12 (R_n + R_{n+1})\Delta p \)); \( W_m^{d} \) the plastic dissipation \( \bar{\boldsymbol{\sigma}}:\Delta\boldsymbol{\varepsilon}^p - \bar R\,\Delta p \) and the damage dissipation \( \bar{Y}^{th}_j \Delta d_j \) with the thermodynamic forces \( Y^{th}_j = \tfrac12\,\boldsymbol{\sigma}:\boldsymbol{\Lambda}_j \); \( W_m^{r} \) the rest (overbars: increment averages). \( W_m = W_m^r + W_m^{ir} + W_m^d \) at every increment.

Material Parameters (props, 18):

Index	Symbol	Description
props[0]	axis	Fibre axis (1, 2 or 3)
props[1]	\( E_L \)	Longitudinal Young’s modulus
props[2]	\( E_T \)	Transverse Young’s modulus
props[3]	\( \nu_{TL} \)	Minor Poisson ratio \( \nu_{LT} E_T/E_L \)
props[4]	\( \nu_{TT} \)	Transverse Poisson ratio
props[5]	\( G_{LT} \)	In-plane shear modulus
props[6]	\( \alpha_L \)	Longitudinal CTE
props[7]	\( \alpha_T \)	Transverse CTE
props[8]	\( Y_{12,0} \)	Shear damage threshold
props[9]	\( Y_{12,c} \)	Shear damage characteristic force
props[10]	\( Y_{22,0} \)	Transverse damage threshold
props[11]	\( Y_{22,c} \)	Transverse damage characteristic force
props[12]	\( Y_{22,u} \)	Transverse failure force (brittle)
props[13]	\( b \)	Transverse-shear coupling of \( Y_{ts} \)
props[14]	\( A_{ts} \)	Transverse weight of the plasticity criterion
props[15]	\( \sigma_{ts,0} \)	Initial yield stress
props[16]	\( m_{ts} \)	Hardening exponent
props[17]	\( \beta_{ts} \)	Hardening coefficient

State Variables (statev, 10):

Index	Symbol	Description
statev[0]	\( T_{init} \)	Reference temperature
statev[1]	\( d_{22} \)	Transverse damage
statev[2]	\( d_{12} \)	Shear damage
statev[3]	\( p \)	Accumulated plastic strain
statev[4-9]	\( \boldsymbol{\varepsilon}^p \)	Plastic strain (Voigt, engineering shears, global frame)

// Carbon/epoxy (T300/914-like), fibres along 1
vec props = {1, 138000., 11000., 0.0223, 0.40, 5500., -0.3e-6, 30e-6,
             0.10, 0.30, 0.15, 0.40, 1.20, 2.5, 2.0, 50., 0.5, 300.};
vec statev = zeros(10);

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.

Parameters:

umat_name – Name of the constitutive model
Etot – Total strain at the beginning of the increment (6)
DEtot – Strain increment (6)
stress – Stress [in: start; out: end] (6)
Lt – Tangent modulus [output] (6x6)
L – Damaged elastic stiffness \( \tilde{\mathbf{L}} \) [output] (6x6)
DR – Rotation increment (3x3)
nprops – Number of material properties (18)
props – Material properties (table above)
nstatev – Number of state variables (10)
statev – State variables [input/output] (table above)
T – Temperature at the beginning of the increment
DT – Temperature increment
Time – Time at the beginning of the increment
DTime – Time increment
Wm – Total mechanical work [input/output]
Wm_r – Recoverable work [input/output]
Wm_ir – Stored (hardening) work [input/output]
Wm_d – Dissipated work (plasticity and damage) [input/output]
ndi – Number of direct stress components
nshr – Number of shear stress components
start – First increment flag
tnew_dt – Suggested time step ratio [output] (0.2 when a local loop does not converge)
tangent_mode – Not used: continuum tangent (see above)

void umat_damage_weibull(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &stress, 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, const int &tangent_mode = tangent_default)

void umat_external(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &stress, 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)

External mechanical UMAT subroutine wrapper.

Wrapper function for external user-defined material subroutines (UMATs) loaded as plugins at runtime. The actual constitutive model is implemented by the user through the plugin API (see umat_plugin_api.hpp).

Material parameters and state variables are user-defined and depend on the specific constitutive model implemented in the plugin.

See also

umat_plugin_ext_api for the native simcoon plugin base class

See also

umat_plugin_aba_api for Abaqus-compatible plugins

See also

umat_plugin_ans_api for ANSYS-compatible plugins

Parameters:

Etot – Total strain tensor at beginning of increment (Voigt notation: 6x1 vector)
DEtot – Strain increment tensor (Voigt notation: 6x1 vector)
stress – Stress tensor (Voigt notation: 6x1 vector) [output]
Lt – Consistent tangent modulus (6x6 matrix) [output]
L – Elastic stiffness tensor (6x6 matrix) [output]
sigma_in – Internal stress contribution for explicit solvers (6x1 vector) [output]
DR – Rotation increment matrix (3x3) 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 [output]
Wm_d – Dissipated 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 (Newton), 1=explicit (RNL), 2=dynamic implicit
tnew_dt – Suggested new time step ratio for adaptive time stepping [output]

void umat_plasticity_chaboche(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &stress, 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, const int &tangent_mode = tangent_default)

Elastic-plastic constitutive model with Chaboche kinematic hardening (J2, two Armstrong-Frederick branches, Voce isotropic term).

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} = \mathbf{X}_1 + \mathbf{X}_2 \) is the total backstress (two AF terms in this implementation)
\( 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 \( H_p = R(p) = Q \left( 1 - e^{-b p} \right) \) (the closed form of \( \dot{H}_p = b (Q - H_p) \dot{p} \), evaluated from the accumulated plastic strain \( p \) at every iterate; modular VoceHardening), where:

\( Q \) is the saturation value of the isotropic hardening stress
\( b \) 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}^{2} \left( \tfrac{2}{3} C_i - D_i\, \mathbf{X}_i : \mathbf{n} \right) + b\,(Q - H_p) \]

is the total hardening modulus combining the two kinematic contributions and the Voce isotropic term \( R'(p) = b (Q - H_p) \).

Material Parameters (props):

Fixed layout — exactly ONE Voce isotropic term and TWO Armstrong–Frederick backstresses (nprops = 10). For a variable number of terms use the modular UMAT (MODUL) with ChabocheHardening / CombinedVoceHardening.

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]	\( Q \)	Voce saturation stress	Stress	0-500 MPa
props[5]	\( b \)	Voce hardening rate	1/Strain	1-100
props[6]	\( C_1 \)	Kinematic modulus of backstress 1	Stress	10-500 GPa
props[7]	\( D_1 \)	Dynamic recovery of backstress 1	-	0-3000
props[8]	\( C_2 \)	Kinematic modulus of backstress 2	Stress	10-500 GPa
props[9]	\( D_2 \)	Dynamic recovery of backstress 2	-	0-3000

Notes on Parameter Selection:

First backstress typically has high \( C_1 \), low \( D_1 \) (long-range)
Second backstress typically has lower \( C_2 \), higher \( D_2 \) (short-range)
Set \( Q = b = 0 \) for pure kinematic hardening

State Variables (statev):

Total: \( n_{statev} = 33 \) ( \( T_{init} + p + \boldsymbol{\varepsilon}^p + \boldsymbol{a}_1 + \boldsymbol{a}_2 + \mathbf{X}_1 + \mathbf{X}_2 + H_p \)).

Index	Symbol	Description	Units
statev[0]	\( T_{init} \)	Initial/reference temperature	Temperature
statev[1]	\( p \)	Accumulated plastic strain	Strain
statev[2..7]	\( \boldsymbol{\varepsilon}^p \)	Plastic strain (Voigt, engineering shear)	Strain
statev[8..13]	\( \boldsymbol{a}_1 \)	Back-strain 1 (Voigt, engineering shear)	Strain
statev[14..19]	\( \boldsymbol{a}_2 \)	Back-strain 2 (Voigt, engineering shear)	Strain
statev[20..25]	\( \mathbf{X}_1 = \tfrac{2}{3} C_1 \boldsymbol{a}_1 \)	Backstress 1 (Voigt)	Stress
statev[26..31]	\( \mathbf{X}_2 = \tfrac{2}{3} C_2 \boldsymbol{a}_2 \)	Backstress 2 (Voigt)	Stress
statev[32]	\( H_p \)	Voce isotropic hardening stress \( R(p) \) (output)	Stress

// Example usage: 316 stainless steel (1 Voce term + 2 backstresses, fixed)
vec props(10);
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

// Isotropic hardening (Voce, Hp = Q (1 - exp(-b p)))
props(4) = 100;         // Q = 100 MPa (saturation stress)
props(5) = 10;          // b = 10 (hardening rate)

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

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

vec statev = zeros(33);  // see the state-variable table
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 stress = 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(Etot, DEtot, stress, 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 \times 6\))

See also

Idev() for deviatoric projection tensor

See also

Fischer_Burmeister() for complementarity Solver Functions

See also

eta_stress() for shifted stress 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

Integrator. The isotropic term is \( H_p = R(p) = Q(1 - e^{-bp}) \) (modular VoceHardening, statev[32] is its output) in every mode. Under tangent_mode 0-2 the return mapping is the convex cutting plane loop below (flow at each iterate); under tangent_closest_point (3, full 3D states) the kernel runs the modular PlasticityMechanism (von Mises + VoceHardening + ChabocheHardening with 2 branches) on its own props and state through closest_point_return_mapping(): backward-Euler backstresses, flow at the converged stress, and cpp_consistent_tangent() — the exact (non-symmetric, dynamic recovery) operator the cutting-plane algorithmic tangent only approximates. Non-convergence requests tnew_dt = 0.5.

Note

Stress measure. The stress returned by this model (the stress argument, written \( \boldsymbol{\sigma} \) in the relations above) is the Cauchy stress under infinitesimal strain; under finite strain the update runs in a corotational frame, so it is the rotated Kirchhoff stress \( \hat{\boldsymbol{\tau}} = \boldsymbol{Q}^{T}\boldsymbol{\tau}\,\boldsymbol{Q} \) on the frame fixed by the chosen objective rate ( \( \boldsymbol{Q} = \boldsymbol{R} \) for Green–Naghdi and \( \log_R \), the logarithmic frame for the XBM/log rate, \( \boldsymbol{F} \) for \( \log_F \)).

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 \times 1\) vector)
DEtot – Strain increment tensor (Voigt notation: \(6 \times 1\) vector)
stress – Stress tensor (Voigt notation: \(6 \times 1\) vector) [output]
Lt – Consistent tangent modulus \( \mathbf{L}_t \) ( \(6 \times 6\) matrix) [output]
L – Elastic stiffness tensor ( \(6 \times 6\) matrix) [output]
sigma_in – Internal stress contribution for explicit solvers ( \(6 \times 1\) vector) [output]
DR – Rotation increment matrix ( \(3 \times 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(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &stress, 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, const int &tangent_mode = tangent_default)

Elastic-plastic constitutive model with isotropic hardening (J2, power-law hardening).

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 \), evaluated through PowerLawHardening: exact for \( p \geq 10^{-6} \) and a \( C^1 \) quadratic blend below, because for \( m < 1 \) the exact slope \( m k p^{m-1} \) is infinite at the onset and the cutting-plane Newton cycles on the first plastic increment (committed off-surface, silently)
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 ( \(\times 2\) in Voigt)	-
statev[6]	\( \varepsilon^p_{13} \)	Plastic strain component 13 ( \(\times 2\) in Voigt)	-
statev[7]	\( \varepsilon^p_{23} \)	Plastic strain component 23 ( \(\times 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 stress = zeros(6);
mat Lt = zeros(6,6);
mat L = zeros(6,6);
mat DR = eye(3,3);
vec props = {70000, 0.3, 1e-5, 200, 500, 0.2};
vec statev = zeros(8);

umat_plasticity_iso("EPICP", Etot, DEtot, stress, Lt, L, DR, 6, props, 8, statev,
                        20.0, 0.0, 0.0, 1.0, Wm, Wm_r, Wm_ir, Wm_d,
                        3, 3, true, tnew_dt);

See also

Fischer_Burmeister_m() for the complementarity Solver Functions

See also

L_iso() for isotropic elastic stiffness construction

See also

eta_stress() for plastic flow direction computation

Note

Integrator. Under tangent_mode 0-2 the return mapping is the convex cutting plane loop below; under tangent_closest_point (3, full 3D states) the kernel runs the modular PlasticityMechanism (von Mises + PowerLawHardening) on its own props and state through closest_point_return_mapping() and returns cpp_consistent_tangent(). For this radial law the two integrators converge to the same state; the closest-point branch carries the C1 onset regularisation of the power law (PowerLawHardening) and never commits an unconverged state (non-convergence requests tnew_dt = 0.5).

Note

Stress measure. The stress returned by this model (the stress argument, written \( \boldsymbol{\sigma} \) in the relations above) is the Cauchy stress under infinitesimal strain; under finite strain the update runs in a corotational frame, so it is the rotated Kirchhoff stress \( \hat{\boldsymbol{\tau}} = \boldsymbol{Q}^{T}\boldsymbol{\tau}\,\boldsymbol{Q} \) on the frame fixed by the chosen objective rate ( \( \boldsymbol{Q} = \boldsymbol{R} \) for Green–Naghdi and \( \log_R \), the logarithmic frame for the XBM/log rate, \( \boldsymbol{F} \) for \( \log_F \)).

Note

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

Note

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

Note

The tangent modulus Lt is always computed

Parameters:

umat_name – Name of the constitutive model (EPICP)
Etot – Total strain tensor at beginning of increment (Voigt notation: \(6 \times 1\) vector)
DEtot – Strain increment tensor (Voigt notation: \(6 \times 1\) vector)
stress – Stress tensor (Voigt notation: \(6 \times 1\) vector) [output]
Lt – Consistent tangent modulus \( \mathbf{L}_t = \frac{\partial \boldsymbol{\sigma}}{\partial \boldsymbol{\varepsilon}} \) ( \(6 \times 6\) matrix) [output]
L – Elastic stiffness tensor ( \(6 \times 6\) matrix) [output]
DR – Rotation increment matrix ( \(3 \times 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)
tnew_dt – Suggested new time step size for adaptive time stepping [output]

void umat_plasticity_johnson_cook_CCP(const std::string &umat_name, 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, const int &tangent_mode = tangent_default)

Johnson-Cook elastic-viscoplastic constitutive law (EPJCK), mechanical version.

J2 (von Mises) plasticity with associated flow and the Johnson-Cook yield stress (Johnson & Cook, 1983), which multiplies a power-law strain hardening by a logarithmic strain-rate factor and a thermal-softening factor:

\[ \sigma_Y\left(p,\dot p,\theta\right) = \left(A + B\,p^{n}\right) \left(1 + C\,\ln\frac{\dot p}{\dot\varepsilon_0}\right) \left(1 - \theta^{*m}\right), \qquad \theta^{*} = \frac{\theta - \theta_{\mathrm{ref}}}{\theta_{\mathrm{melt}} - \theta_{\mathrm{ref}}} . \]

The yield function, flow rule and loading/unloading conditions are

\[ \Phi = \bar{\sigma} - \sigma_Y\left(p,\dot p,\theta\right) \leq 0, \qquad \dot{\boldsymbol{\varepsilon}}^{p} = \dot p\,\boldsymbol{\Lambda}, \quad \boldsymbol{\Lambda} = \frac{3}{2}\frac{\boldsymbol{\sigma}'}{\bar{\sigma}}, \qquad \dot p \geq 0, \quad \dot p\,\Phi = 0 , \]

with the stress given by \( \boldsymbol{\sigma} = \mathbf{L} : \left(\boldsymbol{\varepsilon} - \boldsymbol{\alpha}\,(\theta - \theta_0) - \boldsymbol{\varepsilon}^{p}\right) \).

Regularization of the Johnson-Cook law. The logarithmic rate factor is clamped at its reference value, \( \ln(\dot p/\dot\varepsilon_0) \to \max\left(\ln(\dot p/\dot\varepsilon_0), 0\right) \): below \( \dot\varepsilon_0 \) (and in particular at the onset of yielding, \( \dot p = 0 \)) the material responds with its quasi-static, reference-rate yield stress (a \( C^0 \) kink of \( \Phi \) at \( \dot p = \dot\varepsilon_0 \), the definition of the law). The homologous temperature is clamped to \( [0, 1) \) so that the law is defined below \( \theta_{\mathrm{ref}} \) and never reaches the singular melting point; for \( m < 2 \) the softening \( T^{*m} \) is continued below \( T^* = 10^{-3} \) by a \( C^2 \) cubic (johnson_cook_softening), so that its temperature derivatives, read by the entropy and the heat of the thermomechanical law, stay bounded at \( \theta_{\mathrm{ref}} \). The power-law hardening \( B p^{n} \) with \( n < 1 \) has an infinite slope at the onset, which makes the cutting-plane Newton overshoot and cycle on the first plastic increment: it is evaluated through PowerLawHardening, i.e. exact for \( p \geq 10^{-6} \) and a \( C^1 \) quadratic blend below.

Discrete update. The plastic strain rate is treated fully implicitly over the increment, \( \dot p = \Delta p / \Delta t \). The rate factor therefore depends on the unknown \( \Delta p \) and its derivative enters the scalar Jacobian of the CCP iteration:

\[ K = \frac{\partial \Phi}{\partial p} + \frac{\partial \Phi}{\partial \Delta p} = -\,n B p^{n-1}\, f_{\dot\varepsilon}\, f_\theta - \left(A + B p^{n}\right) \frac{C}{\dot p\,\Delta t}\, f_\theta , \]

so the iterates converge on the rate-dependent yield surface. With \( \Delta t = 0 \) the law degenerates to its rate-independent form. \( \Phi(\Delta p) \) is not convex across the kink \( \Delta p = \dot\varepsilon_0 \Delta t \) of the clamped rate factor, and a Newton started at \( \Delta p = 0 \) cycles between the two branches once \( \Delta t \lesssim C (A + B p^n) / (3 G \dot\varepsilon_0) \) (about \( 5\cdot10^{-5} \) s for AISI 4340). The kernel therefore decides the branch first, from \( \Phi \) at the kink: when positive, the loop starts at the kink with the right derivative of the rate factor and stays on the convex logarithmic branch, where Newton is monotone; otherwise it runs on the quasi-static branch from \( \Delta p = 0 \). Neither the unknown nor the law is altered.

Tangent operator. Assembled by compute_tangent_operator() from the converged \( \hat{B} \), \( \boldsymbol{\kappa} = \mathbf{L}:\boldsymbol{\Lambda} \) and \( \partial\Phi/\partial\boldsymbol{\sigma} \), in the three modes of parameter.hpp:

tangent_none (0): \( \mathbf{L}_t = \mathbf{L} \), for explicit integration schemes;
tangent_continuum (1): \( \mathbf{L}_t = \mathbf{L} - \hat{B}^{-1}\, \boldsymbol{\kappa}\otimes\left(\mathbf{L}:\partial\Phi/\partial\boldsymbol{\sigma}\right) \);
tangent_algorithmic (2, default): the Simo-Hughes consistent operator; the flow Hessian is \( \partial\boldsymbol{\Lambda}/\partial\boldsymbol{\sigma} \) (deta_stress) and the rate term is already part of \( \hat{B} \), so the operator is the exact Jacobian of the backward-Euler (closest-point) map at fixed \( \Delta t \); the cutting-plane state coincides with it for radial increments (J2, isotropic). It is a fixed- \( \Delta t \) derivative: at a given strain rate it tends to \( \mathbf{L} \) as \( \Delta t \to 0 \).

Energy split. The stored energy of the hardening is the reference-rate work of the thermally softened hardening stress, \( G^{ir}(\theta, p) = f_\theta(\theta)\, B p^{n+1}/(n+1) \), so the hardening force is \( A_p = -\partial G^{ir}/\partial p = -f_\theta B p^{n} \) and the rate factor is a viscous overstress acting on the dissipation:

\[ \gamma = \left(\bar{\sigma} + A_p\right)\dot p = f_\theta \left[A f_{\dot\varepsilon} + B p^{n}\left(f_{\dot\varepsilon} - 1\right)\right]\dot p \geq 0 \]

at every temperature and rate (with \( A_p = -B p^n \), EPICP’s isothermal form, the dissipation turns negative once \( f_\theta < B p^n / (A + B p^n) \)). The increments are \( \Delta W_m^{ir} = -\tfrac12 (A_p^n + A_p^{n+1})\,\Delta p \) and \( \Delta W_m^{d} = \tfrac12 (\boldsymbol{\sigma}^n + \boldsymbol{\sigma}^{n+1}):\Delta\boldsymbol{\varepsilon}^{p} + \tfrac12 (A_p^n + A_p^{n+1})\,\Delta p \), \( A_p \) being evaluated at the start and end temperatures.

Material properties (props, 11):

Index	Symbol	Description
0	\( E \)	Young’s modulus
1	\( \nu \)	Poisson’s ratio
2	\( \alpha \)	Coefficient of thermal expansion (isotropic)
3	\( A \)	Initial yield stress
4	\( B \)	Hardening coefficient
5	\( n \)	Hardening exponent
6	\( C \)	Strain-rate sensitivity
7	\( \dot\varepsilon_0 \)	Reference strain rate (1/s)
8	\( m \)	Thermal-softening exponent
9	\( \theta_{\mathrm{ref}} \)	Reference temperature (K)
10	\( \theta_{\mathrm{melt}} \)	Melting temperature (K)

Admissible ranges, not checked by the kernel: \( \dot\varepsilon_0 > 0 \), \( \theta_{\mathrm{melt}} > \theta_{\mathrm{ref}} \), \( C \geq 0 \), \( n > 0 \), \( m > 0 \) (a zero denominator or a zero reference rate propagates NaN silently). A cutting-plane loop that reaches maxiter_umat off the surface sets tnew_dt to 0.5 (step cut).

State variables (statev, 9):

Index	Symbol	Description
0	\( \theta_0 \)	Initial temperature
1	\( p \)	Accumulated plastic strain
2-7	\( \boldsymbol{\varepsilon}^{p} \)	Plastic strain (Voigt 11 22 33 12 13 23)
8	\( \dot p \)	Plastic strain rate over the last increment (output)

The same statev layout as EPICP (plastic strain at offset 2) is declared in the finite-strain convention table of umat_smart.cpp, so EPJCK is served by the log-strain / Kirchhoff box route under every objective rate.

Parameters:

umat_name – Name of the constitutive law (EPJCK)
Etot – Total strain at \( t_n \) (6x1 Voigt)
DEtot – Strain increment (6x1 Voigt)
sigma – Stress [input/output] (6x1 Voigt)
Lt – Tangent modulus [output] (6x6), see tangent_mode
L – Elastic stiffness tensor [output] (6x6)
DR – Rotation increment (3x3)
nprops – Number of material properties (11)
props – Material properties vector
nstatev – Number of state variables (9)
statev – State variables vector [input/output]
T – Temperature at \( t_n \)
DT – Temperature increment
Time – Time at \( t_n \)
DTime – Time increment (drives the plastic strain rate \( \dot p = \Delta p/\Delta t \))
Wm – Total mechanical work [input/output]
Wm_r – Recoverable mechanical work [input/output]
Wm_ir – Irrecoverable (stored) mechanical work [input/output]
Wm_d – Dissipated mechanical work [input/output]
ndi – Number of direct stress components
nshr – Number of shear stress components
start – Flag for the first increment
tnew_dt – Suggested time step ratio [output]
tangent_mode – Tangent assembly mode: tangent_none (0, elastic L for explicit integration), tangent_continuum (1) or tangent_algorithmic (2, default)

void umat_sma_mono(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &stress, 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, const int &tangent_mode = tangent_default)

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.

Elastic Symmetry Selection:

The elastic behavior is determined by the umat_name parameter:

umat_name	Symmetry	Elastic Parameters
SMAMO	Isotropic	E, nu
SMAMC	Cubic	C11, C12, C44
SMAOT	Orthotropic	E1, E2, E3, nu12, nu13, nu23, G12, G13, G23
SMATI	Transverse Isotropic	EL, ET, nuTL, nuTT, GLT

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
Supports multiple elastic symmetry types

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 \to M_n\)):**

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

**Reverse Transformation ( \(M_n \to 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 \(\to\) Orthorhombic: N = 6 variants
Cubic \(\to\) Monoclinic (NiTi): N = 12 or 24 variants
Cubic \(\to\) Tetragonal: N = 3 variants

Material Parameters (props):

The props vector layout depends on the elastic symmetry (umat_name):

SMAMO (Isotropic): 16 parameters

Index	Symbol	Description
props[0]	E	Young’s modulus
props[1]	nu	Poisson’s ratio
props[2]	alpha_iso	Coefficient of thermal expansion
props[3]	b	Slope parameter
props[4]	g	Shear strain magnitude
props[5]	Ms	Martensite start temperature
props[6]	Af	Austenite finish temperature
props[7]	nvariants	Number of martensite variants
props[8-15]	c_lambda0, p0_lambda0, …	Lagrange multiplier parameters

SMAMC (Cubic): 17 parameters

Index	Symbol	Description
props[0]	C11	Elastic constant C11
props[1]	C12	Elastic constant C12
props[2]	C44	Elastic constant C44
props[3]	alpha_iso	Coefficient of thermal expansion
props[4]	b	Slope parameter
props[5]	g	Shear strain magnitude
props[6]	Ms	Martensite start temperature
props[7]	Af	Austenite finish temperature
props[8]	nvariants	Number of martensite variants
props[9-16]	c_lambda0, p0_lambda0, …	Lagrange multiplier parameters

State Variables (statev):

Total state variables required: \( n_{statev} = 1 + N + 7 \) (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:N+6]	\( \varepsilon^{tr} \)	Macroscopic transformation strain (Voigt)	Strain
statev[N+7]	\( \xi \)	Total martensite volume fraction	-

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.

See also

L_iso() for isotropic stiffness tensor (SMAMO)

See also

L_cubic() for cubic stiffness tensor (SMAMC)

See also

L_ortho() for orthotropic stiffness tensor (SMAOT)

See also

L_isotrans() for transverse isotropic stiffness tensor (SMATI)

Note

Stress measure. The stress returned by this model (the stress argument, written \( \boldsymbol{\sigma} \) in the relations above) is the Cauchy stress under infinitesimal strain; under finite strain the update runs in a corotational frame, so it is the rotated Kirchhoff stress \( \hat{\boldsymbol{\tau}} = \boldsymbol{Q}^{T}\boldsymbol{\tau}\,\boldsymbol{Q} \) on the frame fixed by the chosen objective rate ( \( \boldsymbol{Q} = \boldsymbol{R} \) for Green–Naghdi and \( \log_R \), the logarithmic frame for the XBM/log rate, \( \boldsymbol{F} \) for \( \log_F \)).

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:

umat_name – Name of the constitutive model (SMAMO, SMAMC, SMAOT, SMATI)
Etot – Total strain tensor at beginning of increment (Voigt notation: \(6 \times 1\) vector)
DEtot – Strain increment tensor (Voigt notation: \(6 \times 1\) vector)
stress – Cauchy stress tensor (Voigt notation: \(6 \times 1\) vector) [output]
Lt – Consistent tangent modulus ( \(6 \times 6\) matrix) [output]
L – Elastic stiffness tensor ( \(6 \times 6\) matrix) [output]
DR – Rotation increment matrix ( \(3 \times 3\)) for objective integration
nprops – Number of material properties
props – Material properties vector (layout depends on umat_name)
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_sma_unified_T(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &stress, 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, const int &tangent_mode = tangent_default)

Unified phenomenological SMA model with phase transformation.

This function implements the unified phenomenological model for shape memory alloys based on Chatziathanasiou’s framework. The model describes the macroscopic thermomechanical behavior through a single scalar internal variable (martensitic volume fraction \( \xi \)) and a tensorial transformation strain \( \boldsymbol{\varepsilon}^{tr} \).

Elastic Symmetry and Criteria Selection:

The elastic behavior and transformation criteria are determined by the umat_name parameter:

umat_name	Elasticity	Criteria	Number of props
SMADI	Isotropic	Drucker	28
SMADC	Cubic	Drucker	30
SMAAI	Isotropic	Anisotropic Drucker (DFA)	35
SMAAC	Cubic	Anisotropic Drucker (DFA)	37

Effective Properties:

The effective compliance is obtained via a Reuss mixing rule:

\[\mathbf{M}_{eff} = \xi \mathbf{M}_M + (1 - \xi) \mathbf{M}_A \]

and the effective CTE is similarly interpolated.

Current Maximum Transformation Strain:

The maximum transformation strain magnitude \( H^{cur} \) evolves with stress:

\[H^{cur}(\bar{\sigma}) = H_{min} + (H_{max} - H_{min})(1 - e^{-k_1 \bar{\sigma} / \sigma_{crit}}) \]

Material Parameters (props):

SMADI (Isotropic elasticity, Drucker criteria) — 28 props:

Index	Symbol	Description	Units
0	flagT	Temperature extrapolation (0=linear, 1=smooth)	-
1	\( E_A \)	Young’s modulus of Austenite	MPa
2	\( E_M \)	Young’s modulus of Martensite	MPa
3	\( \nu_A \)	Poisson’s ratio of Austenite	-
4	\( \nu_M \)	Poisson’s ratio of Martensite	-
5	\( \alpha_A \)	CTE of Austenite	1/K
6	\( \alpha_M \)	CTE of Martensite	1/K
7	\( H_{min} \)	Minimal transformation strain magnitude	-
8	\( H_{max} \)	Maximal transformation strain magnitude	-
9	\( k_1 \)	Exponential evolution parameter for \( H^{cur} \)	-
10	\( \sigma_{crit} \)	Critical stress for \( H^{cur} \) evolution	MPa
11	\( C_A \)	Clausius-Clapeyron slope (M \( \rightarrow \) A)	MPa/K
12	\( C_M \)	Clausius-Clapeyron slope (A \( \rightarrow \) M)	MPa/K
13	\( M_{s0} \)	Martensite start temperature at zero stress	K
14	\( M_{f0} \)	Martensite finish temperature at zero stress	K
15	\( A_{s0} \)	Austenite start temperature at zero stress	K
16	\( A_{f0} \)	Austenite finish temperature at zero stress	K
17	\( n_1 \)	Martensite start smooth exponent	-
18	\( n_2 \)	Martensite finish smooth exponent	-
19	\( n_3 \)	Austenite start smooth exponent	-
20	\( n_4 \)	Austenite finish smooth exponent	-
21	\( \sigma_{caliber} \)	Calibration stress for \( C_A \) / \( C_M \)	MPa
22	\( b \)	Tension-compression asymmetry parameter (Prager)	-
23	\( n \)	Tension-compression asymmetry exponent (Prager)	-
24	\( c_\lambda \)	Penalty function exponent start point	-
25	\( p_{0\lambda} \)	Penalty function limit value	-
26	\( n_\lambda \)	Penalty function power law exponent	-
27	\( \alpha_\lambda \)	Penalty function power law parameter	-

SMADC (Cubic elasticity, Drucker criteria) — 30 props:

Index	Symbol	Description	Units
0	flagT	Temperature extrapolation (0=linear, 1=smooth)	-
1	\( E_A \)	Young’s modulus of Austenite ([100] direction)	MPa
2	\( E_M \)	Young’s modulus of Martensite ([100] direction)	MPa
3	\( \nu_A \)	Poisson’s ratio of Austenite	-
4	\( \nu_M \)	Poisson’s ratio of Martensite	-
5	\( G_A \)	Shear modulus of Austenite	MPa
6	\( G_M \)	Shear modulus of Martensite	MPa
7–29		Same as SMADI props[5–27] (common SMA parameters)

SMAAI (Isotropic elasticity, anisotropic Drucker criteria) — 35 props:

Same as SMADI (props[0–27]) followed by 7 DFA parameters:

Index	Symbol	Description
28	\( F_{dfa} \)	F parameter of DFA criteria
29	\( G_{dfa} \)	G parameter of DFA criteria
30	\( H_{dfa} \)	H parameter of DFA criteria
31	\( L_{dfa} \)	L parameter of DFA criteria
32	\( M_{dfa} \)	M parameter of DFA criteria
33	\( N_{dfa} \)	N parameter of DFA criteria
34	\( K_{dfa} \)	K parameter of DFA criteria

SMAAC (Cubic elasticity, anisotropic Drucker criteria) — 37 props:

Same as SMADC (props[0–29]) followed by 7 DFA parameters (props[30–36]).

State Variables (statev) — 17 variables:

Index	Symbol	Description	Units
0	\( T_{init} \)	Initial/reference temperature	K
1	\( \xi \)	Martensitic volume fraction	-
2	\( \varepsilon^{tr}_{11} \)	Transformation strain component 11	-
3	\( \varepsilon^{tr}_{22} \)	Transformation strain component 22	-
4	\( \varepsilon^{tr}_{33} \)	Transformation strain component 33	-
5	\( \gamma^{tr}_{12} \)	Transformation strain component 12 (engineering)	-
6	\( \gamma^{tr}_{13} \)	Transformation strain component 13 (engineering)	-
7	\( \gamma^{tr}_{23} \)	Transformation strain component 23 (engineering)	-
8	\( \xi_F \)	Forward martensitic volume fraction	-
9	\( \xi_R \)	Reverse martensitic volume fraction	-
10	\( \rho \Delta s_0 \)	Entropy difference between phases (M - A)	MPa/K
11	\( \rho \Delta E_0 \)	Internal energy difference between phases (M - A)	MPa
12	\( D \)	Stress dependence parameter for transformation limits	-
13	\( a_1 \)	Forward hardening parameter	MPa
14	\( a_2 \)	Reverse hardening parameter	MPa
15	\( a_3 \)	Equilibrium hardening parameter	MPa
16	\( Y_{0t} \)	Initial transformation critical value	MPa

References:

Chatziathanasiou, D. (2016). Ph.D. Thesis — Phenomenological constitutive modeling of shape memory alloys.
Chemisky, Y., Chatziathanasiou, D., Kumar, P., & Lagoudas, D. C. (2014). “A constitutive model for cyclic actuation of high-temperature shape memory alloys.” Mechanics of Materials, 68, 120-136.

See also

L_iso() for isotropic stiffness tensor (SMADI, SMAAI)

See also

L_cubic() for cubic stiffness tensor (SMADC, SMAAC)

See also

umat_sma_mono() for the micromechanical monocrystal SMA model

Note

Elastic convention: isotropic uses “Enu”, cubic uses “EnuG” for L_iso / L_cubic

Note

The flagT parameter controls temperature extrapolation: 0 for linear, 1 for smooth

Note

The pre-2.0 aliases SMAUT and SMANI were removed in 2.0: use SMADI and SMAAI.

Note

Convergence measure. The local Fischer-Burmeister solve reports \( \sum_j |FB_j| / Y^{crit}_j \), so each \( Y^{crit}_j \) must be a positive scale of its criterion. The signed thresholds \( Y^t_F = Y_0^t + D\,H^{cur}\,\bar{\sigma} \) and \( Y^t_R = Y_0^t + D\,\boldsymbol{\sigma}:\bar{\boldsymbol{\varepsilon}}^{T} \) are not: with a calibrated \( D < 0 \) they cross zero along a loading path (at \( \bar{\sigma} = Y_0^t / (-D\,H^{cur}) \), about 1.2 GPa for NiTi-like parameters), and a residual at the floating-point floor of \( \xi \) then reads as a non-convergence: the kernel requests step cuts no sub-increment can satisfy and the solver crawls at its minimal increment. The scales are therefore the magnitudes of the terms, \( Y^{crit}_F = |Y_0^t| + |D\,H^{cur}\,\bar{\sigma}| \) and \( Y^{crit}_R = |Y_0^t| + |D\,\boldsymbol{\sigma}:\bar{\boldsymbol{\varepsilon}}^{T}| \), identical to the thresholds wherever \( D \ge 0 \). Shared by the thermomechanical twin and by umat_sma_unified_TR.

Parameters:

umat_name – Model variant (SMADI, SMADC, SMAAI, SMAAC)
Etot – Total strain tensor at beginning of increment (Voigt notation: 6x1)
DEtot – Strain increment tensor (Voigt notation: 6x1)
stress – Cauchy stress tensor (Voigt notation: 6x1) [output]
Lt – Consistent tangent modulus (6x6) [output]
L – Elastic stiffness tensor (6x6) [output]
DR – Rotation increment matrix (3x3) for objective integration
nprops – Number of material properties
props – Material properties vector (layout depends on umat_name)
nstatev – Number of state variables (17)
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 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 ratio for adaptive time stepping [output]

void umat_sma_unified_TR(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &stress, 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, const int &tangent_mode = tangent_default)

Unified phenomenological SMA model with phase transformation and reorientation.

Extends umat_sma_unified_T with a third leading mechanism: martensite reorientation. The macroscopic state is the martensitic volume fraction \( \xi \), the transformation strain \( \boldsymbol{\varepsilon}^{tr} \), the reorientation strain accumulator \( \mathbf{E}^{Reo} \), and three back-strain channels \( \mathbf{a}_F, \mathbf{a}_R, \mathbf{a}_{reo} \) that grow with the forward / reverse / reorientation increments respectively.

Variants:

umat_name	Elasticity	Criteria	nprops
SMRDI	Isotropic	Drucker	35
SMRDC	Cubic	Drucker	37
SMRAI	Isotropic	Anisotropic Drucker (DFA)	42
SMRAC	Cubic	Anisotropic Drucker (DFA)	44

Effective back-stress and saturation-coupled reorientation surface (Form B):

\[\mathbf{X} = H^{Reo} \, (\mathbf{a} \odot \mathbf{I}_{r05}), \qquad \mathbf{a} = \mathbf{a}_F + \mathbf{a}_R + \mathbf{a}_{reo} \]

\[\Phi^{Reo} = \mathrm{Prager}\!\bigl( \boldsymbol{\sigma} - (1+\lambda_{1}^{Reo}) \, \mathbf{X} \bigr) - Y^{Reo} \]

where \( \lambda_{1}^{Reo} = \mathrm{lagrange}_{\mathrm{pow}_1}(\|\mathbf{v}^{re}\|/E_{T}^{Reo,max}, \ldots) \) is the saturation penalty on the intrinsic (per-unit-martensite) back-strain \( \mathbf{v}^{re} \).

Property layout for SMRDI (28 + 7 = 35 props):

Index	Symbol	Description
0..27	–	Same as SMADI (flagT, EA, EM, nuA, nuM, alphaA, alphaM, Hmin, Hmax, k1, sigmacrit, C_A, C_M, Ms0, Mf0, As0, Af0, n1..n4, sigmacaliber, prager_b, prager_n, c_lambda, p0_lambda, n_lambda, alpha_lambda)
28	\( Y^{Reo} \)	Stress limit for onset of reorientation
29	\( H^{Reo} \)	Reorientation kinematic hardening coefficient
30	\( E_{T}^{Reo,max} \)	Maximum reorientation back-strain magnitude
31	\( c_{\lambda Reo} \)	Reorientation penalty exponent start point
32	\( p_{0\lambda Reo} \)	Reorientation penalty limit value
33	\( n_{\lambda Reo} \)	Reorientation penalty power exponent
34	\( \alpha_{\lambda Reo} \)	Reorientation penalty power parameter

For SMRDC, SMRAI, SMRAC the same 7 reorientation parameters are appended after the SMADC / SMAAI / SMAAC property block.

State variables (statev) — 30 entries:

Per Chatziathanasiou §4.3, the internal variables of the model are \( \{\boldsymbol\sigma, T, \boldsymbol\varepsilon^F, \boldsymbol\varepsilon^R, \boldsymbol\varepsilon^{re}, \mathbf{v}^{re}, \xi^F, \xi^R\} \). There is only one back-strain in the entire model — \( \mathbf{v}^{re} \) — associated with the reorientation channel. Forward and reverse transformation use isotropic hardening only and do not carry back-strains.

Index	Description
0..16	Identical to umat_sma_unified_T (T_init, xi, ET(6), xiF, xiR, calibrated rho-Ds0, rho-DE0, D, a1, a2, a3, Y0t)
17	\( p^{TR} \) : cumulative reorientation multiplier \( \sum \Delta s^{Reo} \)
18..23	\( \mathbf{v}^{re} \) : reorientation back-strain (Voigt, strain convention with engineering shears)
24..29	\( \mathbf{E}^{Reo} \) : macroscopic reorientation strain accumulator

References:

Chatziathanasiou, D. (2016). Ph.D. Thesis.
Chatziathanasiou, Chemisky, Meraghni, Echchorfi, Patoor (2015). Smart Materials and Structures.

See also

umat_sma_unified_T for the transformation-only base model.

See also

assemble_continuum_tangent for the shared 3-mechanism tangent operator.

Note

Stress measure. The stress returned by this model (the stress argument, written \( \boldsymbol{\sigma} \) in the relations above) is the Cauchy stress under infinitesimal strain; under finite strain the update runs in a corotational frame, so it is the rotated Kirchhoff stress \( \hat{\boldsymbol{\tau}} = \boldsymbol{Q}^{T}\boldsymbol{\tau}\,\boldsymbol{Q} \) on the frame fixed by the chosen objective rate ( \( \boldsymbol{Q} = \boldsymbol{R} \) for Green–Naghdi and \( \log_R \), the logarithmic frame for the XBM/log rate, \( \boldsymbol{F} \) for \( \log_F \)).

Note

SMRDI/SMRDC/SMRAI/SMRAC route through umat_sma_unified_TR.

Note

Internal vectorial state variables are rotated as strain by DR.

Note

Active set. The local system couples the three mechanisms in one Fischer-Burmeister solve. The reorientation row has purely kinematic hardening, \( B_{22} \simeq -(1+\lambda_1^{Reo})\,H^{Reo}\,\boldsymbol{\eta}^{Reo}:\boldsymbol{\eta}^{Reo} \), and at the apex of the Drucker cone ( \( \boldsymbol{\sigma}_{eff} = 0 \): a thermal step at zero stress; or \( \xi = 0 \) with \( H^{Reo} = 0 \)) the direction \( \boldsymbol{\eta}^{Reo} \) is undefined and both the stress-flow and the back-strain terms of the row vanish, so \( B_{22} = 0 \) and the row \( -|B_{22}|\,\Delta s^{Reo} = 0 \) is singular. The transformation rows keep the isotropic term \( \partial \Phi / \partial \xi \) on their diagonal and never degenerate. A mechanism with \( \Phi_j < 0 \) and \( B_{jj} = 0 \) is therefore dropped from the iteration, \( \Delta s_j = 0 \), and the reduced system is solved: the row of an inactive mechanism with a finite diagonal decouples and gives \( \Delta s_j = 0 \) exactly, so the reduced solve returns the same increments as the full one wherever the full one is solvable.

Note

Convergence measure: same magnitude scales \( Y^{crit}_j \) as umat_sma_unified_T (see its note); \( Y^{crit}_{Reo} = Y^{Reo} \).

Note

When \( \Delta s^{Reo} = 0 \), the result is bit-identical to umat_sma_unified_T (same elastic mixing, same forward/reverse equations).

void umat_prony_Nfast(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &stress, 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, const int &tangent_mode = tangent_default)

Generalized Maxwell (Prony series) linear viscoelastic model with N branches (PRONK).

Isotropic, rate-dependent, linear viscoelastic model with thermal expansion: an equilibrium spring in parallel with \( N \) Maxwell branches, each a spring \( \mathbf{L}_i \) in series with a dashpot \( \mathbf{H}_i \).

Rheological representation:

           L_inf = L_0 - sum_i L_i   (equilibrium spring)
─────────────────────/\/\/\─────────────────────
         │                                 │
         ├──── L_1 /\/\/\──[ H_1 ]──────────┤   branch 1
         │                                 │
         │                ...              │
         │                                 │
         └──── L_N /\/\/\──[ H_N ]──────────┘   branch N

Constitutive equations (Voigt notation, engineering shear strains):

The instantaneous stiffness acts on the elastic part of the strain,

\[ \boldsymbol{\sigma} = \mathbf{L}_0 : \left( \boldsymbol{\varepsilon} - \boldsymbol{\alpha}\,(T - T_{init}) - \tilde{\boldsymbol{\varepsilon}}^{v} \right), \qquad \tilde{\boldsymbol{\varepsilon}}^{v} = \sum_{i=1}^N \mathbf{M}_0 : \mathbf{L}_i : \boldsymbol{\varepsilon}^{v}_i , \]

with \( \mathbf{L}_0 = \mathbf{L}_{iso}(E_0, \nu_0) \), \( \mathbf{M}_0 = \mathbf{L}_0^{-1} \) and \( \mathbf{L}_i = \mathbf{L}_{iso}(E_i, \nu_i) \); this is the parallel assembly \( \boldsymbol{\sigma} = \mathbf{L}_0 : (\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^{th}) - \sum_i \mathbf{L}_i : \boldsymbol{\varepsilon}^{v}_i \). The dashpot of branch \( i \) carries the branch stress \( \mathbf{A}_i = \mathbf{L}_i : (\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^{v}_i) \) and flows as

\[ \dot{\boldsymbol{\varepsilon}}^{v}_i = \mathbf{H}_i^{-1} : \mathbf{L}_i : \left( \boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^{v}_i \right), \qquad \mathbf{H}_i = 3\,\eta_{B,i}\,\mathbf{I}_{vol} + 2\,\eta_{S,i}\,\mathbf{I}_{dev} , \]

so each branch has a bulk and a shear relaxation time \( \tau_{B,i} = \eta_{B,i} / K_i \) and \( \tau_{S,i} = \eta_{S,i} / \mu_i \), with \( K_i = E_i / (3 (1 - 2\nu_i)) \) and \( \mu_i = E_i / (2 (1 + \nu_i)) \).

Instantaneous versus long-term stiffness (convention):

\( (E_0, \nu_0) \) define the instantaneous (glassy) stiffness: at \( t = 0 \) every \( \boldsymbol{\varepsilon}^{v}_i \) vanishes and \( \boldsymbol{\sigma} = \mathbf{L}_0 : \boldsymbol{\varepsilon} \). As \( t \to \infty \) each \( \boldsymbol{\varepsilon}^{v}_i \) tends to \( \boldsymbol{\varepsilon} \) and the response relaxes to

\[ \mathbf{L}_\infty = \mathbf{L}_0 - \sum_{i=1}^N \mathbf{L}_i , \]

i.e. the uniaxial relaxation modulus is the Prony series \( E(t) = E_\infty + \sum_i E_i\, e^{-t/\tau_i} \) written with \( E(0) = E_0 = E_\infty + \sum_i E_i \). The branch moduli must satisfy \( \sum_i \mathbf{L}_i < \mathbf{L}_0 \) ( \( \sum_i E_i < E_0 \) for equal Poisson ratios); nothing validates it, and a violation gives a negative long-term stiffness (the stress crosses zero during a hold). A calibration that provides the long-term modulus must pass \( E_0 = E_\infty + \sum_i E_i \).

Integration:

Backward Euler, solved in closed form (maxwell_parallel_step, linear_viscoelastic.hpp): the branches decouple, \( \boldsymbol{\varepsilon}^v_i = \mathbf{B}_i^{-1}\boldsymbol{\varepsilon}^v_{i,n} + \mathbf{C}_i\mathbf{L}_i\boldsymbol{\varepsilon}^e \) with \( \mathbf{C}_i = \Delta t\,(\mathbf{H}_i + \Delta t\,\mathbf{L}_i)^{-1} \), and the consistent tangent is

\[ \mathbf{L}_t = \mathbf{L}_0 - \sum_{i=1}^N \mathbf{L}_i \mathbf{C}_i \mathbf{L}_i \]

(exact, whatever tangent_mode but tangent_none, which returns \( \mathbf{L}_0 \)). The multipliers \( v_i \) accumulate \( \| \Delta\boldsymbol{\varepsilon}^v_i \| \). When \( \Delta t \le \) iota the branches are inactive ( \( \mathbf{C}_i = \mathbf{0} \)): the solver probes the tangent with a zero time increment at every block start and commits that answer, so a kernel returning its stationary condition there would relax one branch per block boundary.

**Material parameters (props), \( n_{props} = 4 + 4N \):**

Index	Symbol	Description	Units
props[0]	\( E_0 \)	Instantaneous Young’s modulus	Stress
props[1]	\( \nu_0 \)	Instantaneous Poisson’s ratio	-
props[2]	\( \alpha \)	Isotropic thermal expansion coefficient	1/Temperature
props[3]	\( N \)	Number of Maxwell branches	-
props[4+4(i-1)]	\( E_i \)	Young’s modulus of branch i	Stress
props[5+4(i-1)]	\( \nu_i \)	Poisson’s ratio of branch i	-
props[6+4(i-1)]	\( \eta_{B,i} \)	Bulk viscosity of branch i	Stress x Time
props[7+4(i-1)]	\( \eta_{S,i} \)	Shear viscosity of branch i	Stress x Time

**State variables (statev), \( n_{statev} = 7 + 7N \):**

Index	Symbol	Description	Units
statev[0]	\( T_{init} \)	Initial temperature	Temperature
statev[1:6]	\( \tilde{\boldsymbol{\varepsilon}}^{v} \)	Total viscous strain \( \sum_i \mathbf{M}_0 : \mathbf{L}_i : \boldsymbol{\varepsilon}^{v}_i \) (Voigt)	-
statev[7+7(i-1)]	\( v_i \)	Accumulated flow length of branch i	-
statev[8+7(i-1):13+7(i-1)]	\( \boldsymbol{\varepsilon}^{v}_i \)	Viscous strain of branch i (Voigt)	-

// One Maxwell branch: E_0 = 3000 MPa, nu_0 = 0.35, no thermal expansion,
// E_1 = 1500 MPa, nu_1 = 0.35, etaB_1 = 3000 MPa.s, etaS_1 = 1200 MPa.s
// (long-term modulus E_0 - E_1 = 1500 MPa).
vec props = {3000., 0.35, 0., 1., 1500., 0.35, 3000., 1200.};
vec statev = zeros(7 + 7*1);

See also

umat_prony_Nfast_T() thermomechanical twin (props prefixed by \( \rho, c_p \))

See also

umat_zener_Nfast() same props layout, different rheology

See also

ViscoelasticMechanism the modular (MODUL) twin of this kernel

See also

H_iso() for the bulk/shear viscosity tensor

Note

Stress measure. The stress returned by this model (the stress argument, written \( \boldsymbol{\sigma} \) in the relations above) is the Cauchy stress under infinitesimal strain; under finite strain the update runs in a corotational frame, so it is the rotated Kirchhoff stress \( \hat{\boldsymbol{\tau}} = \boldsymbol{Q}^{T}\boldsymbol{\tau}\,\boldsymbol{Q} \) on the frame fixed by the chosen objective rate ( \( \boldsymbol{Q} = \boldsymbol{R} \) for Green–Naghdi and \( \log_R \), the logarithmic frame for the XBM/log rate, \( \boldsymbol{F} \) for \( \log_F \)).

Note

Not equivalent to umat_zener_Nfast(): identical props layout, but that kernel is a generalized Kelvin chain (branches in series with the spring), not a generalized Maxwell model

Note

The thermal strain enters the stress, not the branch flow: branch i is driven by \( \boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^{v}_i \)

Note

Relaxation times should span the loading time scales; the time step should stay small against the shortest one for accuracy

Parameters:

umat_name – Name of the constitutive law (unused)
Etot – Total strain at the beginning of the increment (Voigt notation: \(6 \times 1\) vector)
DEtot – Strain increment (Voigt notation: \(6 \times 1\) vector)
stress – Stress (Voigt notation: \(6 \times 1\) vector) [input/output]
Lt – Consistent tangent modulus ( \(6 \times 6\) matrix) [output]
L – Instantaneous elastic stiffness \( \mathbf{L}_0 \) ( \(6 \times 6\) matrix) [output]
DR – Rotation increment matrix ( \(3 \times 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 the beginning of the increment
DT – Temperature increment
Time – Time at the beginning of the increment
DTime – Time increment
Wm – Total mechanical work [input/output, cumulative]
Wm_r – Recoverable (stored) work [input/output, cumulative]
Wm_ir – Irrecoverable stored work, always 0 for this model [input/output]
Wm_d – Dissipated (viscous) work [input/output, cumulative]
ndi – Number of direct stress components (typically 3)
nshr – Number of shear stress components (typically 3)
start – Flag indicating the first increment (true) or a continuation (false)
tnew_dt – Suggested new time step size for adaptive time stepping [output]
tangent_mode – Tangent operator selection (see parameter.hpp)

void umat_zener_fast(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &stress, 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, const int &tangent_mode = tangent_default)

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 \times 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

Stress measure. The stress returned by this model (the stress argument, written \( \boldsymbol{\sigma} \) in the relations above) is the Cauchy stress under infinitesimal strain; under finite strain the update runs in a corotational frame, so it is the rotated Kirchhoff stress \( \hat{\boldsymbol{\tau}} = \boldsymbol{Q}^{T}\boldsymbol{\tau}\,\boldsymbol{Q} \) on the frame fixed by the chosen objective rate ( \( \boldsymbol{Q} = \boldsymbol{R} \) for Green–Naghdi and \( \log_R \), the logarithmic frame for the XBM/log rate, \( \boldsymbol{F} \) for \( \log_F \)).

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 \times 1\) vector)
DEtot – Strain increment tensor (Voigt notation: \(6 \times 1\) vector)
stress – Stress tensor (Voigt notation: \(6 \times 1\) vector) [output]
Lt – Consistent tangent modulus ( \(6 \times 6\) matrix) [output]
DR – Rotation increment matrix ( \(3 \times 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 std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &stress, 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, const int &tangent_mode = tangent_default)

Generalized Kelvin chain (Zener) linear viscoelastic model with N units (ZENNK).

Isotropic, rate-dependent, linear viscoelastic model with thermal expansion: an instantaneous spring \( \mathbf{L}_0 \) in series with \( N \) Kelvin–Voigt units, each a spring \( \mathbf{L}_i \) in parallel with a dashpot \( \mathbf{H}_i \) (generalized Kelvin chain; \( N = 1 \) is the Poynting–Thomson / Zener solid).

Rheological representation:

                    ┌── L_1 /\/\/\──┐      ┌── L_N /\/\/\──┐
────/\/\/\──────────┤               ├──────┤               ├────
     L_0            └──[ H_1 ]──────┘ ...  └──[ H_N ]──────┘
(instantaneous)         unit 1                  unit N

Constitutive equations (Voigt notation, engineering shear strains):

The instantaneous spring carries the whole stress on the elastic part of the strain,

\[ \boldsymbol{\sigma} = \mathbf{L}_0 : \left( \boldsymbol{\varepsilon} - \boldsymbol{\alpha}\,(T - T_{init}) - \boldsymbol{\varepsilon}^{v} \right), \qquad \boldsymbol{\varepsilon}^{v} = \sum_{i=1}^N \boldsymbol{\varepsilon}^{v}_i , \]

with \( \mathbf{L}_0 = \mathbf{L}_{iso}(E_0, \nu_0) \) and \( \mathbf{L}_i = \mathbf{L}_{iso}(E_i, \nu_i) \). Each Kelvin unit sees that same stress and its dashpot carries what its spring does not,

\[ \dot{\boldsymbol{\varepsilon}}^{v}_i = \mathbf{H}_i^{-1} : \left( \boldsymbol{\sigma} - \mathbf{L}_i : \boldsymbol{\varepsilon}^{v}_i \right), \qquad \mathbf{H}_i = 3\,\eta_{B,i}\,\mathbf{I}_{vol} + 2\,\eta_{S,i}\,\mathbf{I}_{dev} , \]

so each unit has a bulk and a shear retardation time \( \tau_{B,i} = \eta_{B,i} / K_i \) and \( \tau_{S,i} = \eta_{S,i} / \mu_i \), with \( K_i = E_i / (3 (1 - 2\nu_i)) \) and \( \mu_i = E_i / (2 (1 + \nu_i)) \).

Instantaneous versus long-term response:

\( (E_0, \nu_0) \) define the instantaneous stiffness (every \( \boldsymbol{\varepsilon}^{v}_i \) vanishes at \( t = 0 \)). Under a constant stress each unit creeps to \( \boldsymbol{\varepsilon}^{v}_i = \mathbf{M}_i : \boldsymbol{\sigma} \), so the long-term response is the series compliance

\[ \mathbf{M}_\infty = \mathbf{M}_0 + \sum_{i=1}^N \mathbf{M}_i , \qquad \mathbf{M}_i = \mathbf{L}_i^{-1} , \]

i.e. the uniaxial creep compliance is \( J(t) = 1/E_0 + \sum_i (1 - e^{-t/\tau_i}) / E_i \). There is no admissibility constraint between the unit moduli and \( E_0 \) (a chain only gets more compliant), unlike the generalized Maxwell model of umat_prony_Nfast().

Integration:

Backward Euler, solved in closed form (kelvin_series_step, linear_viscoelastic.hpp): with \( \mathbf{C}_i = \Delta t\,(\mathbf{H}_i + \Delta t\,\mathbf{L}_i)^{-1} \) the implicit step is the linear system \( (\mathbf{I} + \mathbf{L}_0 \sum_i \mathbf{C}_i)\,\boldsymbol{\sigma} = \mathbf{L}_0 (\boldsymbol{\varepsilon}^e - \sum_i \mathbf{B}_i^{-1} \boldsymbol{\varepsilon}^v_{i,n}) \), and the consistent tangent is \( \mathbf{L}_t = (\mathbf{M}_0 + \sum_i \mathbf{C}_i)^{-1} \) (exact, whatever tangent_mode but tangent_none, which returns \( \mathbf{L}_0 \)). The multipliers \( v_i \) accumulate \( \| \Delta\boldsymbol{\varepsilon}^v_i \| \). When \( \Delta t \le \) iota the branches are inactive ( \( \mathbf{C}_i = \mathbf{0} \)): the solver probes the tangent with a zero time increment at every block start and commits that answer, so a kernel returning its stationary condition there would relax one branch per block boundary.

**Material parameters (props), \( n_{props} = 4 + 4N \):**

Index	Symbol	Description	Units
props[0]	\( E_0 \)	Instantaneous Young’s modulus	Stress
props[1]	\( \nu_0 \)	Instantaneous Poisson’s ratio	-
props[2]	\( \alpha \)	Isotropic thermal expansion coefficient	1/Temperature
props[3]	\( N \)	Number of Kelvin units	-
props[4+4(i-1)]	\( E_i \)	Young’s modulus of branch i	Stress
props[5+4(i-1)]	\( \nu_i \)	Poisson’s ratio of branch i	-
props[6+4(i-1)]	\( \eta_{B,i} \)	Bulk viscosity of branch i	Stress x Time
props[7+4(i-1)]	\( \eta_{S,i} \)	Shear viscosity of branch i	Stress x Time

**State variables (statev), \( n_{statev} = 7 + 7N \):**

Index	Symbol	Description	Units
statev[0]	\( T_{init} \)	Initial temperature	Temperature
statev[1:6]	\( \boldsymbol{\varepsilon}^{v} \)	Total viscous strain \( \sum_i \boldsymbol{\varepsilon}^{v}_i \) (Voigt)	-
statev[7+7(i-1)]	\( v_i \)	Accumulated flow length of branch i	-
statev[8+7(i-1):13+7(i-1)]	\( \boldsymbol{\varepsilon}^{v}_i \)	Viscous strain of branch i (Voigt)	-

// One Kelvin unit: E_0 = 3000 MPa, nu_0 = 0.35, no thermal expansion,
// E_1 = 1500 MPa, nu_1 = 0.35, etaB_1 = 3000 MPa.s, etaS_1 = 1200 MPa.s
// (long-term modulus 1/(1/E_0 + 1/E_1) = 1000 MPa).
vec props = {3000., 0.35, 0., 1., 1500., 0.35, 3000., 1200.};
vec statev = zeros(7 + 7*1);

See also

umat_zener_Nfast_T() thermomechanical twin (props prefixed by \( \rho, c_p \))

See also

umat_prony_Nfast() same props layout, different rheology

See also

umat_zener_fast() the single-unit version (ZENER: no N slot, \( n_{statev} = 8 \))

See also

H_iso() for the bulk/shear viscosity tensor

Note

Stress measure. The stress returned by this model (the stress argument, written \( \boldsymbol{\sigma} \) in the relations above) is the Cauchy stress under infinitesimal strain; under finite strain the update runs in a corotational frame, so it is the rotated Kirchhoff stress \( \hat{\boldsymbol{\tau}} = \boldsymbol{Q}^{T}\boldsymbol{\tau}\,\boldsymbol{Q} \) on the frame fixed by the chosen objective rate ( \( \boldsymbol{Q} = \boldsymbol{R} \) for Green–Naghdi and \( \log_R \), the logarithmic frame for the XBM/log rate, \( \boldsymbol{F} \) for \( \log_F \)).

Note

Not equivalent to umat_prony_Nfast(): identical props layout, but a generalized Kelvin chain (units in series) instead of a generalized Maxwell model (branches in parallel); the measured response difference for the same numbers is of order 1

Note

Retardation times should span the loading time scales; the time step should stay small against the shortest one for accuracy

Parameters:

umat_name – Name of the constitutive law (unused)
Etot – Total strain at the beginning of the increment (Voigt notation: \(6 \times 1\) vector)
DEtot – Strain increment (Voigt notation: \(6 \times 1\) vector)
stress – Stress (Voigt notation: \(6 \times 1\) vector) [input/output]
Lt – Consistent tangent modulus ( \(6 \times 6\) matrix) [output]
L – Instantaneous elastic stiffness \( \mathbf{L}_0 \) ( \(6 \times 6\) matrix) [output]
DR – Rotation increment matrix ( \(3 \times 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 the beginning of the increment
DT – Temperature increment
Time – Time at the beginning of the increment
DTime – Time increment
Wm – Total mechanical work [input/output, cumulative]
Wm_r – Recoverable (stored) work [input/output, cumulative]
Wm_ir – Irrecoverable stored work, always 0 for this model [input/output]
Wm_d – Dissipated (viscous) work [input/output, cumulative]
ndi – Number of direct stress components (typically 3)
nshr – Number of shear stress components (typically 3)
start – Flag indicating the first increment (true) or a continuation (false)
tnew_dt – Suggested new time step size for adaptive time stepping [output]
tangent_mode – Tangent operator selection (see parameter.hpp)

void get_L_elastic(phase_characteristics&)