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. 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.

Parameters:

stress – Abaqus stress array
ddsdde – Abaqus mechanical tangent array
ddsddt – Abaqus stress-temperature tangent array
drplde – Abaqus heat generation-strain tangent array
drpldt – Abaqus heat generation-temperature tangent
stran – Abaqus total strain array
dstran – Abaqus strain increment array
time – Abaqus time array
dtime – Time increment
temperature – Temperature at start of increment
Dtemperature – Temperature increment
nprops – Number of material properties
props – Abaqus material properties array
nstatev – Number of state variables
statev – Abaqus state variables array
ndi – Number of direct stress components
nshr – Number of shear stress components
drot – Abaqus rotation increment matrix
sigma – [out] simcoon stress vector
dSdE – [out] simcoon mechanical tangent matrix
dSdT – [out] simcoon stress-temperature tangent matrix
drpldE – [out] simcoon heat generation-strain tangent matrix
drpldT – [out] simcoon heat generation-temperature tangent matrix
Etot – [out] simcoon total strain vector
DEtot – [out] simcoon strain increment vector
T – [out] simcoon temperature
DT – [out] simcoon temperature increment
Time – [out] simcoon total time
DTime – [out] simcoon time increment
props_smart – [out] simcoon properties vector
Wm – [out] simcoon mechanical work quantities vector
Wt – [out] simcoon thermal work quantities vector
statev_smart – [out] simcoon state variables vector
DR – [out] simcoon rotation matrix
start – [out] true if this is the first increment

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

Transfer simcoon output to Abaqus format (mechanical, simple).

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

Parameters:

stress – [out] Abaqus stress array
ddsdde – [out] Abaqus tangent operator array
statev – [out] Abaqus state variables array
ndi – Number of direct stress components
nshr – Number of shear stress components
sigma – simcoon stress vector
statev_smart – simcoon state variables vector
Wm – simcoon work quantities vector
Lt – simcoon tangent matrix

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

Full transfer from simcoon to Abaqus format (mechanical).

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

Parameters:

stress – [out] Abaqus stress array
ddsdde – [out] Abaqus tangent operator array
stran – [out] Abaqus total strain array
dstran – [out] Abaqus strain increment array
time – [out] Abaqus time array
dtime – [out] Time increment
temperature – [out] Temperature
Dtemperature – [out] Temperature increment
nprops – Number of material properties
props – [out] Abaqus material properties array
nstatev – Number of state variables
statev – [out] Abaqus state variables array
ndi – Number of direct stress components
nshr – Number of shear stress components
drot – [out] Abaqus rotation increment matrix
sigma – simcoon stress vector
Lt – simcoon tangent matrix
Etot – simcoon total strain vector
DEtot – simcoon strain increment vector
T – simcoon temperature
DT – simcoon temperature increment
Time – simcoon total time
DTime – simcoon time increment
props_smart – simcoon properties vector
Wm – simcoon work quantities vector
statev_smart – simcoon state variables vector
DR – simcoon rotation matrix
start – simcoon start flag

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

Transfer simcoon output to Abaqus format (thermomechanical).

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

Parameters:

stress – [out] Abaqus stress array
ddsdde – [out] Abaqus mechanical tangent array
ddsddt – [out] Abaqus stress-temperature tangent array
drplde – [out] Abaqus heat generation-strain tangent array
drpldt – [out] Abaqus heat generation-temperature tangent
rpl – [out] Abaqus heat generation rate
statev – [out] Abaqus state variables array
ndi – Number of direct stress components
nshr – Number of shear stress components
sigma – simcoon stress vector
statev_smart – simcoon state variables vector
r – simcoon heat generation rate
Wm – simcoon mechanical work quantities vector
Wt – simcoon thermal work quantities vector
dSdE – simcoon mechanical tangent matrix
dSdT – simcoon stress-temperature tangent matrix
drpldE – simcoon heat generation-strain tangent matrix
drpldT – simcoon heat generation-temperature tangent matrix

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

Transfer Ansys USERMAT input arrays to simcoon format (mechanical).

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

Parameters:

stress – Ansys stress array (ncomp)
dstran – Ansys strain increment array (ncomp)
sedEl – Specific elastic strain energy
sedPl – Specific plastic strain energy
epseq – Equivalent plastic strain
statev – Ansys state variables array
props – Ansys material properties array
Time – Total time
DTime – Time increment
temperature – Temperature
Dtemperature – Temperature increment
ncomp – Number of stress/strain components
nprops – Number of material properties
nstatev – Number of state variables
sigma – [out] simcoon stress vector
DEtot – [out] simcoon strain increment vector
Wm – [out] simcoon work quantities vector
statev_smart – [out] simcoon state variables vector
props_smart – [out] simcoon properties vector
T – [out] simcoon temperature
DT – [out] simcoon temperature increment

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

Transfer simcoon output to Ansys USERMAT format (mechanical).

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

Parameters:

stress – [out] Ansys stress array
ddsdde – [out] Ansys tangent operator array
sedEl – [out] Specific elastic strain energy
sedPl – [out] Specific plastic strain energy
statev – [out] Ansys state variables array
ncomp – Number of stress/strain components
sigma – simcoon stress vector
Lt – simcoon tangent matrix
Wm – simcoon work quantities vector
statev_smart – simcoon state variables vector

void umat_damage_LLD_0(const 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) anisotropic damage model with coupled plasticity for composite materials.

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

Transversely isotropic elastic behavior
Anisotropic damage in matrix-dominated modes (transverse tension, in-plane shear)
Coupled transverse-shear plasticity
Thermodynamically consistent damage evolution based on energy release rate
Separate damage variables for different failure modes
No fiber failure (matrix-dominated damage only)

Material Symmetry:

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

Direction 1: Fiber direction (longitudinal, no damage)
Directions 2-3: Transverse plane (isotropic, subject to damage)

Damage Variables:

Two independent scalar damage variables characterize the material state:

\( d_{22} \in [0,1] \): Transverse damage (matrix cracking perpendicular to fibers)
\( d_{12} \in [0,1] \): Shear damage (matrix/interface damage in fiber-matrix plane)

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

Effective Stress Concept:

The effective stress acting on the undamaged material configuration is:

\[\tilde{\boldsymbol{\sigma}} = \mathbf{M}^{-1} : \boldsymbol{\sigma} \]

where \( \mathbf{M} \) is the damage effect tensor:

\[\mathbf{M} = \text{diag}(1, 1-d_{22}, 1-d_{22}, 1-d_{12}, 1-d_{12}, 1) \]

Damaged Elastic Stiffness:

The elastic stiffness degrades with damage:

\[\mathbf{L}(d_{22}, d_{12}) = \mathbf{L}_0 : \mathbf{M} \]

where \( \mathbf{L}_0 \) is the undamaged transversely isotropic stiffness.

Specific moduli degradation:

\( E_2 = E_{2,0} (1 - d_{22}) \) (transverse Young’s modulus)
\( E_3 = E_{3,0} (1 - d_{22}) \) (out-of-plane Young’s modulus)
\( G_{12} = G_{12,0} (1 - d_{12}) \) (in-plane shear modulus)
\( G_{13} = G_{13,0} (1 - d_{12}) \) (out-of-plane shear modulus)
\( E_1 \) (fiber direction) remains constant (no fiber damage)

Thermodynamic Forces:

The energy release rates driving damage evolution are:

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

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

Damage Evolution Laws:

**Transverse Damage ( \(d_{22}\)):**

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

where:

\( Y_{22,0} \) is the damage initiation threshold
\( Y_{22,c} \) is the characteristic energy release rate
\( Y_{22,u} \) is the ultimate energy release rate (failure)
\( b \) is the damage evolution exponent

**Shear Damage ( \(d_{12}\)):**

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

where:

\( Y_{12,0} \) is the shear damage initiation threshold
\( Y_{12,c} \) is the characteristic shear energy release rate

Coupled Transverse-Shear Plasticity:

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

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

where:

\( A_{ts} \) is the transverse-shear coupling parameter
\( \sigma_{ts,0} \) is the initial yield stress
\( \alpha_{ts}, \beta_{ts} \) are hardening parameters
\( p_{ts} \) is the accumulated plastic strain

Plastic Flow Rule:

Associative flow in the transverse-shear plane:

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

Material Parameters (props):

Index	Symbol	Description	Units	Typical Range (CFRP)
props[0]	axis	Fiber orientation axis (1, 2, or 3)	-	1
props[1]	\( E_L \)	Longitudinal Young’s modulus	Stress	100-200 GPa
props[2]	\( E_T \)	Transverse Young’s modulus	Stress	5-15 GPa
props[3]	\( \nu_{TL} \)	Major Poisson’s ratio	-	0.25-0.35
props[4]	\( \nu_{TT} \)	Transverse Poisson’s ratio	-	0.35-0.45
props[5]	\( G_{LT} \)	In-plane shear modulus	Stress	3-8 GPa
props[6]	\( \alpha_L \)	Longitudinal CTE	1/Temperature	-0.5 to 0 e-6 /K
props[7]	\( \alpha_T \)	Transverse CTE	1/Temperature	20-40 e-6 /K
props[8]	\( Y_{12,0} \)	Shear damage initiation threshold	Energy/Volume	0.05-0.2 MPa
props[9]	\( Y_{12,c} \)	Characteristic shear energy release	Energy/Volume	0.1-0.5 MPa
props[10]	\( Y_{22,0} \)	Transverse damage initiation	Energy/Volume	0.1-0.3 MPa
props[11]	\( Y_{22,c} \)	Characteristic transverse energy	Energy/Volume	0.2-0.6 MPa
props[12]	\( Y_{22,u} \)	Ultimate transverse energy	Energy/Volume	0.5-2.0 MPa
props[13]	\( b \)	Transverse damage exponent	-	1-5
props[14]	\( A_{ts} \)	Transverse-shear coupling	-	1-3
props[15]	\( \sigma_{ts,0} \)	Initial yield stress	Stress	30-80 MPa
props[16]	\( \alpha_{ts} \)	Linear hardening parameter	Stress	0-500 MPa
props[17]	\( \beta_{ts} \)	Quadratic hardening parameter	Stress	0-5000 MPa

State Variables (statev):

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

Index	Symbol	Description	Units
statev[0]	\( T_{init} \)	Initial/reference temperature	Temperature
statev[1]	\( d_{22} \)	Transverse damage variable	-
statev[2]	\( d_{12} \)	Shear damage variable	-
statev[3]	\( p_{ts} \)	Accumulated plastic strain	Strain
statev[4]	\( \varepsilon^p_{11} \)	Plastic strain component 11	Strain
statev[5]	\( \varepsilon^p_{22} \)	Plastic strain component 22	Strain
statev[6]	\( \varepsilon^p_{33} \)	Plastic strain component 33	Strain
statev[7]	\( \varepsilon^p_{12} \)	Plastic strain component 12 (engineering)	Strain
statev[8]	\( \varepsilon^p_{13} \)	Plastic strain component 13 (engineering)	Strain
statev[9]	\( \varepsilon^p_{23} \)	Plastic strain component 23 (engineering)	Strain

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

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

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

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

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

See also

L_isotrans() for transversely isotropic stiffness construction

See also

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

See also

rotate_strain() for strain tensor rotation

See also

rotate_stress() for stress tensor rotation

References:

Ladevèze, P., & Le Dantec, E. (1992). “Damage modelling of the elementary ply for laminated composites.” Composites Science and Technology, 43(3), 257-267.
Allix, O., & Ladevèze, P. (1992). “Interlaminar interface modelling for the prediction of delamination.” Composite Structures, 22(4), 235-242.
Ladevèze, P. (1992). “A damage computational method for composite structures.” Computers & Structures, 44(1-2), 79-87.
Matzenmiller, A., Lubliner, J., & Taylor, R. L. (1995). “A constitutive model for anisotropic damage in fiber-composites.” Mechanics of Materials, 20(2), 125-152.
Pinho, S. T., et al. (2012). “Material and structural response of polymer-matrix fibre-reinforced composites.” Journal of Composite Materials, 46(19-20), 2313-2341.

Note

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 LLD model is specifically designed for unidirectional fiber composites

Note

Damage is irreversible and monotonically increasing

Note

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

Note

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

Note

Parameter identification requires multiple test configurations:

Note

- Transverse tension for \(d_{22}\) parameters

Note

- In-plane shear for \(d_{12}\) parameters

Note

- Off-axis tests for plasticity coupling

Note

Convergence requires small load steps once damage initiates

Note

Material axes must be properly oriented relative to global coordinates

Parameters:

Etot – Total strain tensor at beginning of increment (Voigt notation: \(6 \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]
L – Damaged 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 damage [output]
Wm_d – Dissipated work (damage + plasticity) [output]
ndi – Number of direct stress components (typically 3)
nshr – Number of shear stress components (typically 3)
start – Flag indicating first increment (true) or continuation (false)
solver_type – Solver type: 0=implicit, 1=explicit, 2=dynamic implicit
tnew_dt – Suggested new time step size for adaptive time stepping [output]

void umat_damage_weibull(const 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_CCP(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 solved by the Convex Cutting Plane (CCP) algorithm.

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

J2 (von Mises) plasticity with associative flow rule
Multiple Armstrong-Frederick backstresses for nonlinear kinematic hardening
Optional isotropic hardening with Voce law
Dynamic recovery (recall terms) in backstress evolution
Convex Cutting Plane algorithm for return mapping
Thermal expansion effects
Consistent tangent modulus for implicit FE analysis

Constitutive Equations:

The yield function is defined as:

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

where:

\( \sigma_{eq}(\boldsymbol{\eta}) = \sqrt{\frac{3}{2} \boldsymbol{\eta}_{dev} : \boldsymbol{\eta}_{dev}} \) is the von Mises equivalent stress
\( \boldsymbol{\eta} = \boldsymbol{\sigma} - \mathbf{X} \) is the shifted (effective) stress tensor
\( \mathbf{X} = \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, integrated in rate form \( \dot{H}_p = b \left( Q - H_p \right) \dot{p} \) (a single term, so \( R(p) = Q \left( 1 - e^{-b p} \right) \) exactly under monotonic flow): 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 (integrated in rate form \( \dot{H}_p = b (Q - H_p) \dot{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 \)	Integrated Voce isotropic hardening stress	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, rate form dHp = b (Q - Hp) dp)
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_CCP(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

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

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_CCP(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 solved by the Convex Cutting Plane (CCP) algorithm.

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

J2 (von Mises) plasticity with associative flow rule
Isotropic hardening following a power law: \( H_p = k \cdot p^m \)
Convex Cutting Plane algorithm for return mapping
Thermal expansion effects
Consistent tangent modulus for implicit FE analysis

Constitutive Equations:

The yield function is defined as:

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

where:

\( \sigma_{eq} = \sqrt{\frac{3}{2} \mathbf{s} : \mathbf{s}} \) is the von Mises equivalent stress
\( \mathbf{s} \) is the deviatoric stress tensor
\( H_p(p) = k \cdot p^m \) is the isotropic hardening function
\( p = \int \sqrt{\frac{2}{3} \dot{\boldsymbol{\varepsilon}}^p : \dot{\boldsymbol{\varepsilon}}^p} \, dt \) is the accumulated plastic strain
\( \sigma_Y \) is the initial yield stress

The plastic flow rule (associative plasticity):

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

Evolution of accumulated plastic strain:

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

Convex Cutting Plane Algorithm:

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

Find \( \Delta p \geq 0 \) such that \( \Phi(\boldsymbol{\sigma}, p) \leq 0 \) and \( \Delta p \cdot \Phi = 0 \)
This is solved using the Fischer-Burmeister function for robust convergence
The method provides a consistent tangent modulus for implicit finite element analysis

Material Parameters (props):

The material properties vector must contain 6 constants:

Index	Symbol	Description	Units
props[0]	\( E \)	Young’s modulus	Stress
props[1]	\( \nu \)	Poisson’s ratio	-
props[2]	\( \alpha \)	Isotropic thermal expansion coefficient	1/Temperature
props[3]	\( \sigma_Y \)	Initial yield stress	Stress
props[4]	\( k \)	Hardening parameter	Stress
props[5]	\( m \)	Hardening exponent	-

State Variables (statev):

The state variables vector contains 8 internal variables:

Index	Symbol	Description	Units
statev[0]	\( T_{init} \)	Initial temperature	Temperature
statev[1]	\( p \)	Accumulated plastic strain	-
statev[2]	\( \varepsilon^p_{11} \)	Plastic strain component 11	-
statev[3]	\( \varepsilon^p_{22} \)	Plastic strain component 22	-
statev[4]	\( \varepsilon^p_{33} \)	Plastic strain component 33	-
statev[5]	\( \varepsilon^p_{12} \)	Plastic strain component 12 ( \(\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_CCP("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

See also

L_iso() for isotropic elastic stiffness construction

See also

eta_stress() for plastic flow direction computation

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_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

Legacy aliases SMAUT and SMANI map to SMADI and SMAAI respectively

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

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)

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

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

Constitutive Equations:

The relaxation modulus is expressed as a Prony series:

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

where:

\( E_\infty \) is the long-term (equilibrium) modulus
\( E_i \) is the modulus of the i-th Maxwell element
\( \tau_i \) is the relaxation time of the i-th element
\( N \) is the number of Maxwell elements

Incremental Form:

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

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

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

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

Bulk and Shear Decomposition:

The model handles volumetric and deviatoric responses independently:

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

Consistent Tangent Modulus:

For implicit finite element analysis:

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

Material Parameters (props):

For N Maxwell elements, the props vector contains:

Index	Symbol	Description	Units
props[0]	\( N \)	Number of Maxwell elements	-
props[1]	\( E_\infty \)	Equilibrium Young’s modulus	Stress
props[2]	\( \nu_\infty \)	Equilibrium Poisson’s ratio	-
props[3]	\( \alpha \)	Thermal expansion coefficient	1/Temperature
props[4]	\( E_1 \)	Modulus of 1st Maxwell element	Stress
props[5]	\( \tau_1 \)	Relaxation time of 1st element	Time
…	…	…	…
props[4+2(i-1)]	\( E_i \)	Modulus of i-th element	Stress
props[5+2(i-1)]	\( \tau_i \)	Relaxation time of i-th element	Time

State Variables (statev):

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

Index	Symbol	Description	Units
statev[0]	\( T_{init} \)	Initial temperature	Temperature
statev[1:6]	\( \mathbf{q}_1 \)	Internal stresses of 1st element (Voigt)	Stress
statev[7:12]	\( \mathbf{q}_2 \)	Internal stresses of 2nd element (Voigt)	Stress
…	…	…	…
statev[1+6(i-1):6i]	\( \mathbf{q}_i \)	Internal stresses of i-th element (Voigt)	Stress

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

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

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

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

See also

L_iso() for isotropic elastic stiffness construction

See also

Ivol() for volumetric projection tensor

See also

Idev() for deviatoric projection tensor

References:

Ferry, J. D. (1980). Viscoelastic Properties of Polymers (3rd ed.). Wiley.
Simo, J. C., & Hughes, T. J. R. (1998). Computational Inelasticity. Springer.
Park, S. W., & Schapery, R. A. (1999). “Methods of interconversion between linear viscoelastic material functions. Part I—A numerical method based on Prony series.” International Journal of Solids and Structures, 36(11), 1653-1675.

Note

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 Prony series provides an efficient representation for viscoelastic relaxation

Note

Relaxation times should span the expected loading time scales

Note

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

Note

The model assumes small strains and linear viscoelasticity

Note

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

Parameters:

Etot – Total strain tensor at beginning of increment (Voigt notation: \(6 \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_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 Zener (Standard Linear Solid) viscoelastic model with N parallel branches.

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

Rheological Representation:

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

Constitutive Equations:

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

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

where:

\( E_\infty \) is the long-term (equilibrium) modulus
\( E_i \) is the modulus of the i-th Maxwell element
\( \tau_i = \eta_i / E_i \) is the relaxation time of the i-th element
\( N \) is the number of Maxwell elements

Instantaneous Modulus:

At \( t = 0 \):

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

Stress Decomposition:

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

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

where \( \mathbf{q}_i \) are the internal stress-like variables for each Maxwell element.

Internal Variable Evolution:

Each internal variable evolves according to:

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

Incremental Form:

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

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

Consistent Tangent Modulus:

For implicit finite element analysis:

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

Material Parameters (props):

For N Maxwell elements, the props vector contains:

Index	Symbol	Description	Units
props[0]	\( N \)	Number of Maxwell elements	-
props[1]	\( E_\infty \)	Equilibrium Young’s modulus	Stress
props[2]	\( \nu \)	Poisson’s ratio (assumed constant)	-
props[3]	\( \alpha \)	Thermal expansion coefficient	1/Temperature
props[4]	\( E_1 \)	Modulus of 1st Maxwell element	Stress
props[5]	\( \tau_1 \)	Relaxation time of 1st element	Time
…	…	…	…
props[4+2(i-1)]	\( E_i \)	Modulus of i-th element	Stress
props[5+2(i-1)]	\( \tau_i \)	Relaxation time of i-th element	Time

State Variables (statev):

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

Index	Symbol	Description	Units
statev[0]	\( T_{init} \)	Initial temperature	Temperature
statev[1:6]	\( \mathbf{q}_1 \)	Internal stresses of 1st element (Voigt)	Stress
statev[7:12]	\( \mathbf{q}_2 \)	Internal stresses of 2nd element (Voigt)	Stress
…	…	…	…
statev[1+6(i-1):6i]	\( \mathbf{q}_i \)	Internal stresses of i-th element (Voigt)	Stress

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

References:

Zener, C. (1948). Elasticity and Anelasticity of Metals. University of Chicago Press.
Ferry, J. D. (1980). Viscoelastic Properties of Polymers (3rd ed.). Wiley.
Simo, J. C., & Hughes, T. J. R. (1998). Computational Inelasticity. Springer.
Park, S. W., & Schapery, R. A. (1999). “Methods of interconversion between linear viscoelastic material functions.” Int. J. Solids Struct., 36(11), 1653-1675.

See also

umat_zener_fast() for single Maxwell element version

See also

umat_prony_Nfast() for equivalent Prony series formulation

Note

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

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

Note

Typically 5-10 Maxwell elements are sufficient for most polymers

Note

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

Note

The model assumes small strains and linear viscoelasticity

Note

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

Parameters:

Etot – Total strain tensor at beginning of increment (Voigt notation: \(6 \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 get_L_elastic(phase_characteristics&)