Plasticity-based material models

Drucker-Prager model - DruckerPrager

The Drucker-Prager plasticity model (Contributed by Simon Rolshoven, LSC, FENAC, EPFL.) is an isotropic elasto-plastic model based on a yield function

\[f \left(\bsig, \tauY\right) = F\left(\bsig\right) - \tauY\]

with the pressure-dependent equivalent stress

\[F\left(\bsig\right) = \alpha I_1 + \sqrt{J_2}\]

As usual, \(\bsig\) is the stress tensor, \(\tauY\) is the yield stress under pure shear, and \(I_1\) and \(J_2\) are the first invariant and second deviatoric invariant of the stress tensor. The friction coefficient \(\alpha\) is a positive parameter that controls the influence of the pressure on the yield limit, important for cohesive-frictional materials such as concrete, soils or other geomaterials. Regarding Mohr-Coulomb plasticity, relation to cohesion, \(c\), and the angle of friction, \(\theta\), exists for the Drucker-Prager model

\[\begin{split}\alpha = \frac{2\sin\theta}{(3-\sin\theta)\sqrt{3}}\\ \tauY = \frac{6c\cos\theta}{(3-\sin\theta)\sqrt{3}}\end{split}\]

The flow rule is derived from the plastic potential

\[g\left(\bsig\right) = \alphaPsi I_1 + \sqrt{J_2}\]

where \(\alphaPsi\) is the dilatancy coefficient. An associated model with \(\alphaPsi=\alpha\) would overestimate the dilatancy of concrete, so the dilatancy coefficient is usually chosen smaller than the friction coefficient.

The model is described by the equations

\[\begin{split}\bsig & = &\mbf{D} : \left(\eps - \epsp \right) \\ \tauY & = & h(\kappa) \\ \epspd = \dot{\lambda} \frac{\partial g}{\partial \bsig} & = & \dot{\lambda} \left( \alphaPsi \mbf{\delta} + \frac{\mbf{s}}{2\sqrt{J_2}} \right) \\ \dot{\kappa} & = & \sqrt{\frac{2}{3}} \; \| \epspd \|\end{split}\]

and

\[\dot{\lambda} \ge 0, \;\;\; f \left(\bsig, \tauY\right) \le 0, \;\;\; \dot{\lambda}\, f \left(\bsig, \tauY\right) = 0\]

which represent the linear elastic law, hardening law, evolution laws for plastic strain and hardening variable, and the loading-unloading conditions.

In the above, \(\mbf{D}\) is the elastic stiffness tensor, \(\eps\) is the strain tensor, \(\epsp\) is the plastic strain tensor, \(\lambda\) is the plastic multiplier, \(\mbf{\delta}\) is the unit second-order tensor, \(\mbf{s}\) is the deviatoric stress tensor, \(\kappa\) is the hardening variable, and a superior dot marks the derivative with respect to time.

The flow rule has the form given in Eq.~(eq:flowRule) at all points of the conical yield surface with the exception of its vertex, located on the hydrostatic axis.

For the present model, the evolution of the hardening variable can be explicitly linked to the plastic multiplier. Substituting the flow rule (eq:flowRule) into Eq.~(eq:hardening) and computing the norm leads to

\[\dot\kappa = k \dot \lambda\]

with a constant parameter \(k = \sqrt{1/3 + 2 \alphaPsi^2}\), so the hardening variable is proportional to the plastic multiplier. For \(\alpha=\alphaPsi=0\), the associated \(J_2\)-plasticity model is recovered as a special case.

In the simplest case of linear hardening, the hardening function is a linear function of \(\kappa\), given by

(7)\[h(\kappa) = \tau_0 + H E \kappa\]

where \(\tau_0\) is the initial yield stress, and \(H\) is the hardening modulus normalized with the elastic modulus. Alternatively, an exponential hardening function

\[h(\kappa) = \tau_{\mathrm{limit}} + \left( \tau_0 - \tau_{\mathrm{limit}} \right) \mathrm{e}^{-\kappa/\kappac}\]

can be used for a more realistic description of hardening.

The stress-return algorithm is based on the Newton-iteration. In plasticity, this is commonly called Closest-Point-Projection (CPP), and it generally leads to quadratic convergence. The implemented algorithm is convergent in any stress case, but in the vicinity of the vertex region, quadratic convergence might be lost because of insufficient regularity of the yield function.

The algorithmic tangent stiffness matrix is implemented for both the regular case and the vertex region. Generally, the error decreases quadratically (of course only asymptotically). Again, in the vicinity of the vertex region, quadratic convergence might be lost due to insufficient regularity. Furthermore, the tangent stiffness matrix does not always exist for the vertex case. In these cases, the elastic stiffness is used instead. It is generally safer (but slower) to use the elastic stiffness if you encounter any convergence problems, especially if your problem is tension-dominated.

Table 11 DP material - summary

Description	DP material
Record Format	DruckerPrager num (in) d (rn) tAlpha (rn) E (rn) n (rn) alpha (rn) alphaPsi (rn) ht (in) iys (rn) lys (rn) hm (rn) kc (rn) [yieldtol (rn)]
Parameters	num material model number d material density tAlpha thermal dilatation coefficient E Young modulus n Poisson ratio alpha friction coefficient alphaPsi dilatancy coefficient ht hardening type, 1: linear hardening, 2: exponential hardening iys initial yield stress in shear, \(\tau_0\) lys limit yield stress for exponential hardening, \(\tau_{\mathrm{limit}}\) hm hardening modulus normalized with E-modulus (!) kc \(\kappa_c\) for the exponential softening law yieldtol tolerance of the error in the yield criterion, default value 1.e-14 newtonIter maximum number of iterations in \(\lambda\) search, default value 30
Supported modes	3dMat, PlaneStrain, 3dRotContinuum
Tests/Examples	sm/DruckerPrager_01.in

Drucker-Prager model with tension cut-off and isotropic damage - DruckerPragerCut

The Drucker-Prager plasticity model with tension cut-off is a multisurface model, appropriate for cohesive-frictional materials such as concrete loaded both in compression and tension. The plasticity model is formulated for isotropic hardening and enhanced by isotropic damage, which is driven by the cumulative plastic strain. The model can be used only in the small-strain context, with additive split of the strain tensor into the elastic and plastic parts.

The basic equations include the additive decomposition of strain into elastic and plastic parts,

\[\veps = \veps_{\rm{e}} + \veps_{\rm{p}},\]

the stress strain law

\[\vsig = (1-\omega)\bar{\vsig}=(1-\omega)\mbf{D} :(\veps-\veps_{\rm{p}}),\]

the definition of the yield function in terms of the effective stress,

\[f \left(\bar{\vsig},\kappa \right) = \alpha I_1 + \sqrt{J_2} - \tau_Y,\]

the flow rule

\[g \left( \bar{\vsig} \right) = \alphaPsi I_1 + \sqrt{J_2},\]

the linear hardening law

\[\tau_Y(\kappa) = \tau_0 + H\kappa,\]

where \(\tau_0\) represents the initial yield stress under pure shear, the damage law

(8)\[\omega(\kappa) = \omega_c(1-\mbox{e}^{-a\kappa}),\]

where \(\omega_c\) is critical damage and \(a\) is a positive dimensionless parameter. More detailed descriptioin of some parameters is in Section Drucker-Prager model - DruckerPrager.

The dilatancy coefficient \(\alphaPsi\) controls flow associativeness; if \(\alphaPsi=\alpha\), an associate model is recovered, which overestimates the dilatancy of concrete. Hence, the dilatancy coefficient is usually chosen smaller, \(\alphaPsi \leq \alpha\), and the non-associated model is formulated.

Table 12 Drucker Prager material with tension cut-off - summary

Description	Drucker Prager material with tension cut-off
Record Format	DruckerPragerCut num (in) d (rn) tAlpha (rn) E (rn) n (rn) tau0 (rn) alpha (rn) [alphaPsi (rn)] [H (rn)] [omega_crit (rn)] [a (rn)] [yieldtol (rn)] [NewtonIter (in)]
Parameters	num material model number d material density tAlpha thermal dilatation coefficient E Young modulus n Poisson ratio tau0 initial yield stress in shear \(\tau_0\) alpha friction coefficient alphaPsi dilatancy coefficient, equals to alpha by default H hardening modulus (can be negative in the case of plastic softening), 0 by default omega_crit critical damage in damage law ((8)), 0 by default a exponent in damage law ((8)), 0 by default yieldtol tolerance of the error in the yield criterion, default value 1.e-14 newtonIter maximum number of iterations in \(\lambda\) search, default value 30
Supported modes	1dMat, 3dMat, PlaneStrain

Mises plasticity model with isotropic damage - MisesMat

This model is appropriate for the description of plastic yielding in ductile material such as metals, and it can also cover the effects of void growth. The model uses the Mises yield condition (in terms of the second deviatoric invariant, \(J_2\)), associated flow rule, linear isotropic hardening driven by the cumulative plastic strain, and isotropic damage, also driven by the cumulative plastic strain. The model can be used in the small-strain context, with additive split of the strain tensor into the elastic and plastic parts, or in the large-strain context, with multiplicative split of the deformation gradient and with yield condition formulated in terms of Kirchhoff stress (which is the true Cauchy stress multiplied by the Jacobian).

Small-strain formulation:

The small-strain version of hardening Mises plasticity can be combined with isotropic damage. The basic equations include the additive decomposition of strain into elastic and plastic parts,

(9)\[\veps = \veps_{\rm{e}} + \veps_{\rm{p}},\]

the stress strain law

(10)\[\vsig = (1-\omega)\bar{\vsig}=(1-\omega)\mbf{D} :(\veps-\veps_{\rm{p}}),\]

the definition of the yield function in terms of the effective stress,

(11)\[f(\bar{\vs},\kappa) = \sqrt{\frac{3}{2}\bar{\vs}:\bar{\vs}}-\sigma_Y(\kappa) = \sqrt{3 J_2(\bar{\vsig})}-\sigma_Y(\kappa),\]

the incremental definition of cumulative plastic strain

(12)\[\dot{\kappa} = \| \epspd \|,\]

the linear hardening law (for htype = 0)

(13)\[\sigma_Y(\kappa) = \sigma_0 + H\kappa,\]

the evolution law for the plastic strain

\[\epspd = \dot{\lambda}\frac{\partial f}{\partial \bar{\vs}},\]

the loading-unloading conditions

(14)\[\dot{\lambda} > 0 \qquad f(\bar{\vs},\kappa)\leq 0 \qquad \dot{\lambda} f(\bar{\vs},\kappa) = 0,\]

and the damage law

(15)\[\omega(\kappa) = \omega_c(1-\mbox{e}^{-a\kappa}),\]

In the equations above, \(\veps\) is the strain tensor, \(\veps_{\rm{e}}\) is the elastic strain tensor, \(\veps_{p}\) is the plastic strain tensor, \(\mbf{D}\) is the elastic stiffness tensor, \(\vsig\) is the nominal stress tensor, \(\bar{\vsig}\) is the effective stress tensor, \(\bar{š}\) is the effective deviatoric stress tensor, \(\sigma_Y\) is the magnitude of stress at yielding under uniaxial tension (or compression), \(\kappa\) is the cumulated plastic strain, \(H\) is the hardening modulus, \(\lambda\) is the plastic multiplier, \(\omega\) is the damage variable, \(\omega_c\) is critical damage and \(a\) is a positive dimensionless parameter.

Large-strain formulation

is based on the introduction of an intermediate local configuration, with respect to which the elastic response is characterized. This concept leads to a multiplicative decomposition of deformation gradient into elastic and plastic parts:

\[\mbf{F} = \mbf{F}^e\mbf{F}^p.\]

The stress-evaluation algorithm can be based on the classical radial return mapping; see [SimoHughes] for more details. Damage is not yet implemented in the large-strain version of the model.

The model description and parameters are summarized in Table 13.

Table 13 Mises plasticity – summary.

Description	Mises plasticity model with isotropic hardening
Record Format	MisesMat in d (rn) E (rn) n (rn) sig0 (rn) H (rn) [htype (in)] [h_eps (ra)] [h(eps) (ra)] omega_crit (rn)a (rn)
Parameters	material number d material density E Young’s modulus n Poisson’s ratio sig0 initial yield stress in uniaxial tension (compression) (Required if htype = 0, which is default) H hardening modulus (can be negative in the case of plastic softening) (Required if htype = 0, which is default) htype hardening type (Optional parameter. Default = 0) h_eps array of plastic strains (Required if htype = 1) h(eps) array of yield stresses (Required if htype = 1) omega_crit critical damage in damage law ((15)) a exponent in damage law ((15))
Supported modes	1dMat, PlaneStrain, 3dMat, 3dMatF
Tests/Examples	sm/Mises01.in, sm/Mises02.in, sm/bond_link_2.in, sm/libeam3dboundary.in

VTKxml output can report Mises stress, which equals to \(\sqrt{3J_2}\). When no hardening/softening exists, Mises stress reaches values up to given uniaxial yield stress sig0.

Mises plasticity model with isotropic damage, nonlocal - MisesMatNl, MisesMatGrad

The small-strain version of the model is regularized by the over-nonlocal formulation with damage driven by a combination of local and nonlocal cumulated plastic strain,

(16)\[\hat{\kappa} = (1-m)\kappa + m\bar{\kappa},\]

where \(m\) is a dimensionless material parameter (typically \(m>1\)) and \(\bar{\kappa}\) is the nonlocal cumulated plastic strain, which is evaluated either using the integral approach, or using the implicit gradient approach.

Integral nonlocal formulation

One possible regularization technique is based on the integral definition of nonlocal cumulated plastic strain

(17)\[\bar{\kappa}(x) = \int\limits_V \alpha(x,s)\kappa(s)\,{\rm d}s\]

The nonlocal weight function is usually defined as

\[\alpha(x,s) = \frac{\alpha_0(\Vert x-s\Vert )}{\int\limits_V\alpha_0(\Vert x-t\Vert )\,{\rm d}t}\]

where

(18)\[\begin{split}\alpha_0(r) = \begin{cases} \left(1-\frac{r^2}{R^2}\right)^2 &\text{if $r<R$}\\ \\ 0 & \text{if $r \ge R$} \end{cases}\end{split}\]

is a nonnegative function, for \(r<R\) monotonically decreasing with increasing distance \(r=\Vert x-s\Vert\), and \(V\) denotes the domain occupied by the investigated material body. The key idea is that the damage evolution at a certain point depends not only on the cumulated plastic strain at that point, but also on points at distances smaller than the interaction radius \(R\), considered as a new material parameter.

Implicit gradient formulation

The gradient formulation can be conceived as the differential counterpart to the integral formulation. The nonlocal cumulated plastic strain is computed from a Helmholtz-type differential equation

(19)\[\bar{\kappa} - l^2\nabla^2\bar{\kappa} = \kappa\]

with homogeneous Neumann boundary condition

(20)\[\frac{\partial\bar{\kappa}}{\partial n} = 0.\]

In ((19)), \(l\) is the length scale parameter and \(\nabla\) is the Laplace operator.

The model description and parameters are summarized in Table 14 and Table 15. Note that the internal length parameter r has the meaning of the radius of interaction \(R\) for the integral version (and thus has the dimension of length) but for the gradient version it has the meaning of the coefficient \(l^2\) multiplying the Laplacean in (19), and thus has the dimension of length squared.

Table 14 Nonlocal integral Mises plasticity – summary.

Description	Nonlocal Mises plasticity with isotropic hardening
Record Format	MisesMatNl in d (rn) E (rn) n (rn) sig0 (rn) H (rn) omega_crit (rn)a (rn)r (rn)m (rn)[wft (in)][scalingType (in)]
Parameters	material number d material density E Young’s modulus n Poisson’s ratio sig0 initial yield stress in uniaxial tension (compression) H hardening modulus omega_crit critical damage a exponent in damage law r nonlocal interaction radius \(R\) from eq. ((18)) m over-nonlocal parameter wft selects the type of nonlocal weight function (see Section Nonlocal isotropic damage model for tensile failure - Idmnl1): ; 1 - default, quartic spline (bell-shaped function with bounded support) ; 2 - Gaussian function ; 3 - exponential function (Green function in 1D) ; 4 - uniform averaging up to distance \(R\) ; 5 - uniform averaging over one finite element ; 6 - special function obtained by reducing the 2D exponential function to 1D (by numerical integration) scalingType selects the type of scaling of the weight function (e.g. near a boundary; see Section Nonlocal isotropic damage model for tensile failure - Idmnl1): ; 1 - default, standard scaling with integral of weight function in the denominator ; 2 - no scaling (the weight function normalized in an infinite body is used even near a boundary) ; 3 - Borino scaling (local complement)
Supported modes	1dMat, PlaneStrain, 3dMat

Table 15 Gradient-enhanced Mises plasticity – summary.

Description	Gradient-enhanced Mises plasticity with isotropic damage
Record Format	MisesMatGrad in d (rn) E (rn) n (rn) sig0 (rn) H (rn) omega_crit (rn)a (rn)r (rn)m (rn)
Parameters	material number d material density E Young’s modulus n Poisson’s ratio sig0 initial yield stress in uniaxial tension (compression) H hardening modulus omega_crit critical damage a exponent in damage law r internal length scale parameter \(l^2\) from eq. ((19)) m over-nonlocal parameter
Supported modes	1dMat, PlaneStrain, 3dMat

Rankine plasticity model with isotropic damage and its nonlocal formulations - RankMat, RankMatNl, RankMatGrad

This model has a very similar structure to the model described in Section Mises plasticity model with isotropic damage - MisesMat, but is based on the Rankine yield condition. It is available in the small-strain version only, and so far exclusively for plane stress analysis. The basic equations (9)–(10) and (12)–(15) remain valid, and the yield function (11) is redefined as

(21)\[f(\bar{\vsig},\kappa) = \max_{I}\bar{\sigma_I}-\sigma_Y(\kappa)\]

where \(\bar{\sigma_I}\) are the principal values of the effective stress tensor \(\bar{\vsig}\). The hardening law can have either the linear form (13), or the exponential form

(22)\[\sigma_Y(\kappa) = \sigma_0 + \Delta\sigma_Y\left(1-\exp(-H\kappa/\Delta\sigma_Y)\right),\]

where \(H\) is now the initial plastic modulus and \(\Delta\sigma_Y\) is the value of yield stress increment asymptotically approached as \(\kappa\rightarrow\infty\). In damage law (15), parameter \(\omega_c\) is always set to 1. If the plastic hardening is linear, the user can specify either the exponent \(a\) from (15), or the dissipation per unit volume, \(g_f\), which represents the area under the stress-strain diagram (and parameter \(a\) is then determined automatically such that the area under the diagram has the prescribed value). For exponential plastic hardening, the evaluation of \(a\) from \(g_f\) is not properly implemented and it is better to specify \(a\) directly.

The model description and parameters are summarized in Table 16–Table 18. Note that the default value of parameter \(m\) is equal to 1 for the integral model but for the gradient model it is equal to 2. Also note that the internal length parameter :param{r} has the meaning of the radius of interaction \(R\) for the integral version (and thus has the dimension of length) but for the gradient version it has the meaning of the coefficient \(l^2\) multiplying the Laplacean in (19), and thus has the dimension of length squared.

For the gradient model, it is possible to specify parameter negligible_damage, which sets the minimum value of damage that is considered as nonzero.

The approximate solution of Helmholtz equation \((`\)implicitGradient``)` can lead to very small but nonzero nonlocal kappa at some points that are actually elastic. If such small values are positive, they lead to a very small but nonzero damage. If this is interpreted as “loading”, the tangent terms are activated, but damage will not actually grow at such points and the convergence rate is slowed down. It is better to consider such points as elastic. By default, negligible_damage is set to 0, but it is recommended to set it e.g.to 1.e-6.

Table 16 Rankine plasticity – summary.

Description	Rankine plasticity with isotropic hardening and damage
Record Format	RankMat in d (rn) E (rn) n (rn) plasthardtype (in) sig0 (rn) H (rn) delSigY (rn) yieldtol (rn) (gf (rn) \(\|\) a (rn))
Parameters	material number d material density E Young’s modulus n Poisson’s ratio plasthardtype type of plastic hardening (0=linear=default, 1=exponential) sig0 initial yield stress in uniaxial tension (compression) H initial hardening modulus (default value 0.) delSigY final increment of yield stress (default value 0., needed only if plasthardtype=1) yieldtol relative tolerance in the yield condition gf dissipation per unit volume a exponent in damage law ((15))
Supported modes	PlaneStress

Table 17 Nonlocal integral Rankine plasticity – summary.

Description	Nonlocal Rankine plasticity with isotropic hardening and damage
Record Format	RankMatNl in d (rn) E (rn) n (rn) plasthardtype (in) sig0 (rn) H (rn) delSigY (rn) yieldtol (rn) (gf (rn) \(\|\) a (rn)):elemparam:r{rn}m (rn)[wft (in)][scalingType (in)]
Parameters	material number d material density E Young’s modulus n Poisson’s ratio plasthardtype type of plastic hardening (0=linear=default, 1=exponential) sig0 initial yield stress in uniaxial tension (compression) H initial hardening modulus (default value 0.) delSigY final increment of yield stress (default value 0.) yieldtol relative tolerance in the yield condition gf dissipation per unit volume a exponent in damage law ((15)) r internal length scale parameter \(l^2\) from eq. ((19)) m over-nonlocal parameter (default value 1.) wft selects the type of nonlocal weight function (see Section Nonlocal isotropic damage model for tensile failure - Idmnl1): ; 1 - default, quartic spline (bell-shaped function with bounded support) ; 2 - Gaussian function ; 3 - exponential function (Green function in 1D) ; 4 - uniform averaging up to distance \(R\) ; 5 - uniform averaging over one finite element ; 6 - special function obtained by reducing the 2D exponential function to 1D (by numerical integration) scalingType selects the type of scaling of the weight function (e.g. near a boundary; see Section Nonlocal isotropic damage model for tensile failure - Idmnl1): ; 1 - default, standard scaling with integral of weight function in the denominator ; 2 - no scaling (the weight function normalized in an infinite body is used even near a boundary) ; 3 - Borino scaling (local complement)
Supported modes	PlaneStress

Table 18 Gradient-enhanced Rankine plasticity – summary.

Description	Gradient-enhanced Rankine plasticity with isotropic hardening and damage
Record Format	RankMatGrad in d (rn) E (rn) n (rn) plasthardtype (in) sig0 (rn) H (rn) delSigY (rn) yieldtol (rn) (gf (rn) \(\|\) a (rn)):elemparam:r{rn}m (rn)negligible_damage (rn)
Parameters	material number d material density E Young’s modulus n Poisson’s ratio plasthardtype type of plastic hardening (0=linear=default, 1=exponential) sig0 initial yield stress in uniaxial tension (compression) H hardening modulus (default value 0.) delSigY final increment of yield stress (default value 0.) yieldtol relative tolerance in the yield condition gf dissipation per unit volume a exponent in damage law ((15)) r internal length scale parameter \(l\) from eq. ((19)) m over-nonlocal parameter (default value 2.) negligible_damage optional parameter (default value 0.)
Supported modes	PlaneStress

Perfectly plastic material with Mises yield condition - Steel1

This is an older model, kept here for compatibility with previous versions. It uses Mises plasticity condition with no hardening and under small strain only. The model description and parameters are summarized in Tab.~:samp:ref{Steel1_table}. All its features are included in the model described in Section~:samp:ref{sec:misplast}.

Table 19 Perfectly plastic material with Mises condition – summary

Description	Perfectly plastic material with Mises condition
Record Format	Steel1 num (in) d (rn) E (rn) n (rn) tAlpha (rn) Ry (rn)
Parameters	num material model number d material density E Young modulus n Poisson ratio tAlpha thermal dilatation coefficient Ry uniaxial yield stress
Supported modes	3dMat, PlaneStress, PlaneStrain, 1dMat, 2dPlateLayer, 2dBeamLayer, 3dShellLayer 3dBeam, PlaneStressRot
Tests/Examples	benchmark/sm/steel1.in

Composite plasticity model for masonry - Masonry02

Masonry is a composite material made of bricks and mortar. Nonlinear behavior of both components should be considered to obtain a realistic model able to describe cracking, slip, and crushing of the material. The model is based on paper by Lourenco and Rots [Rots]. It is formulated on the basis of softening plasticity for tension, shear, and compression (see fig.(Fig. 1)). Numerical implementation is based on modern algorithmic concepts such as implicit integration of the rate equations and consistent tangent stiffness matrices.

Fig. 1 Composite yield surface model for masonry

The approach used in this work is based on the idea of concentrating all the damage in the relatively weak joints and, if necessary, in potential tension cracks in the bricks. The joint interface constitutive model should include all important damage mechanisms. Here, the concept of interface elements is used. An interface element allows to incorporate discontinuities in the displacement field and its behavior is described in terms of a relation between the tractions and relative displacement across the interface. In the present work, these quantities will be denoted as \(\sig\), generalized stress, and \(\e\), generalized strain. For 2D configuration, \(\sig=\{\sigma, \tau\}^T\) and \(\e=\{u_n, u_s\}^T\), where \(\sigma\) and \(\tau\) are the normal and shear components of the traction interface vector and \(n\) and \(s\) subscripts distinguish between normal and shear components of displacement vector.

Fig. 2 Modeling strategy for masonry

The elastic response is characterized in terms of the elastic constitutive matrix \(\mbf{D}\) as

\[\sig = \mbf{D} \e\]

For a 2D configuration, \(\mbf{D} = \text{diag}\{k_n, k_s\}\). The terms of the elastic stiffness matrix can be obtained from the properties of both masonry and joints as

\[k_n = \frac{E_b E_m}{t_m (E_b - E_m)};\; k_s = \frac{G_b G_m}{t_m (G_b - G_m)}\]

where \(E_b\) and \(E_m\) are Young’s moduli, \(G_b\) and \(G_m\) are shear moduli for brick and mortar, and \(t_m\) is the thickness of the joint. One should note that there is no contact algorithm assumed between bricks, this means that the overlap of neighboring units will be visible. On the other hand, the interface model includes a compressive cap, where the compressive inelastic behavior of masonry is lumped.

Tension mode

In the tension mode, the exponential softening law is assumed (see Tensile behavior of proposed model (f_t=0.2 MPa, G_f^I=0.018 N/mm)). The yield function has the following form

\[f_1(\sig, \kappa_1) = \sigma - f_t(\kappa_1)\]

where the yield value \(f_t\) is defined as

(23)\[f_t = f_{t0} \exp\left(-\frac{f_{t0}}{G^I_f} \kappa_1\right)\]

Fig. 3 Tensile behavior of proposed model (\(f_t=0.2\) MPa, \(G_f^I=0.018\) N/mm)

The \(f_{t0}\) represents the tensile strength of the joint or interface; and \(G^I_f\) is the mode-I fracture energy. For the tension mode, the associated flow hypothesis is assumed.

Shear mode

For the shear mode a Coulomb friction envelope is used. The yield function has the form

\[f_2(\sig,\kappa_2) = \vert\tau\vert+\sigma\tan\phi(\kappa_2)-c(\kappa_2)\]

According to [Rots] the variations of friction angle \(\phi\) and cohesion \(c\) are assumed as

(24)\[\begin{split}c &= c_0\exp\left(-\frac{c_0}{G^{II}_f}\kappa_2\right)\\ \tan\phi &= \tan\phi_0+(\tan\phi_r-\tan\phi_0)\left(\frac{c_0-c}{c_0}\right)\end{split}\]

where \(c_0\) is initial cohesion of joint, \(\phi_0\) initial friction angle, \(\phi_r\) residual friction angle, and \(G^{II}_f\) fracture energy in mode II failure. A non-associated plastic potential \(g_2\) is considered as

\[g_2=\vert\tau\vert+\sigma\tan\Phi-c\]

Fig. 4 Shear behavior of proposed model for different confinement levels in MPa (\(c_0=0.8\ \rm{MPa},\ \tan\phi_0=1.0,\ \tan\phi_r=0.75,{\rm and}\ G_f^{II}=0.05\ {N/mm}\))

Coupling of tension/shear modes

The tension and Coulomb friction modes are coupled with isotropic softening. This means that the percentage of softening in the cohesion is assumed to be the same as on the tensile strength

\[\dot\kappa_1=\lambda_1+\frac{G^I_f}{G^{II}_f}\frac{c_0}{f_{t0}}\lambda_2;\ \ \dot\kappa_2=\frac{G^{II}_f}{G^I_f}\frac{f_{t0}}{c_0}\lambda_1+\lambda_2\]

This follows from ((23)) and ((24)). However, in the corner region, when both yield surfaces are activated, such approach will lead to a non-acceptable penalty. For this reason a quadratic combination is assumed

\[\dot\kappa_1=\sqrt{(\lambda_1)^2+\left(\frac{G^I_f}{G^{II}_f}\frac{c_0}{f_{t0}}\lambda_2\right)^2};\ \ \dot\kappa_2=\sqrt{\left(\frac{G^{II}_f}{G^I_f}\frac{f_{t0}}{c_0}\lambda_1\right)^2+(\lambda_2)^2}\]

Cap mode

For the cap mode, an ellipsoid interface model is used. The yield condition is assumed as

\[f_3(\sig, \kappa_3) = C_{nn}\sigma^2+C_{ss}\tau^2 + C_n\sigma-\bar{\sigma}^2(\kappa_3)\]

where \(C_{nn},\ C_{ss}\) and \(C_n\) are material model parameters and \(\bar{\sigma}\) is yield value, originally assumed in the following form of hardening/softening law [Rots]

(25)\[\begin{split}\bar{\sigma}_1(\kappa_3) &= \bar{\sigma}_i+(\bar{\sigma}_p-\bar{\sigma}_i)\sqrt{\frac{2\kappa_3}{\kappa_p}-\frac{\kappa_3^2}{\kappa_p^2}};\;\;\kappa_3\in(0,\kappa_p)\\ \bar{\sigma}_2(\kappa_3) &= \bar{\sigma}_p+(\bar{\sigma}_m-\bar{\sigma}_p)\left(\frac{\kappa_3-\kappa_p}{\kappa_m-\kappa_p}\right)^2;\;\;\kappa_3\in(\kappa_p, \kappa_m)\\ \bar{\sigma}_3(\kappa_3) &= \bar{\sigma}_r+(\bar{\sigma}_m-\bar{\sigma}_r)\exp\left(m\frac{\kappa_3-\kappa_m}{\bar{\sigma}_m-\bar{\sigma}_r}\right);\;\;\kappa_3\in(\kappa_m, \infty)\end{split}\]

with \(m=2(\bar{\sigma}_m-\bar{\sigma}_p)/(\kappa_m-\kappa_p)\). The hardening/softening law ((25)) is shown in fig.(Fig. 5). Note that the curved diagram is a \(C^1\) continuous \(\sigma-\kappa_3\) relation. The energy under the load-displacement diagram can be related to a “compressive fracture energy”.

The original hardening law ((25).1) exhibits indefinite slope for \(\kappa_3=0\), which can cause the problems with numerical implementation. This has been overcomed by replacing this hardening law with parabolic equation given by

\[\bar{\sigma}_1(\kappa_3) = \bar{\sigma}_i-2*(\bar{\sigma}_i-\bar{\sigma}_p)*\frac{\kappa_3}{\kappa_p}+(\bar{\sigma}_i-\bar{\sigma}_p)\frac{\kappa_3}{\kappa_p}\]

An associated flow and strain hardening hypothesis are being considered. This yields

\[\dot\kappa_3=\lambda_3\sqrt{(2C_{nn}\sigma+C_n)*(2C_{nn}\sigma+C_n) + (2C_{ss}\tau)*(2C_{ss}\tau)}\]

Fig. 5 Hardening/softening law for cap mode

The model parameters are summarized in Table 20. There is one algorithmic issue, that follows from the model formulation. Since the cap mode hardening/softening is not coupled to hardening/softening of shear and tension modes the it may happen that when the cap and shear modes are activated, the return directions become parallel for both surfaces. This should be avoided by adjusting the input parameters accordingly (one can modify dilatancy angle, for example).

Table 20 Composite model for masonry - summary.

Description	Composite plasticity model for masonry
Record Format	Masonry02 num (in) d (rn) E (rn) n (rn) ft0 (rn) gfi (rn) gfii (rn) kn (rn) ks (rn) c0 (rn) tanfi0 (rn) tanfir (rn) tanpsi (rn) si (rn) sp (rn) sm (rn) sr (rn) kp (rn) km (rn) kr (rn) cnn (rn) css (rn) cn (rn)
Parameters	num material model number d material density E Young modulus n Poisson ratio ft0 tensile strength gfi fracture energy for mode I gfii fracture energy for mode II kn joint elastic property ks joint elastic property c0 initial cohesion tanfi0 initial friction angle tanfir residual friction angle tanpsi dilatancy (si, sp, sm, sr) cap parameters \(\{\bar{\sigma}_i, \bar{\sigma}_p, \bar{\sigma}_m, \bar{\sigma}_r\}\) (kp, km,kr) cap parameters \(\{\kappa_p, \kappa_m, \kappa_r\}\) cnn,:param:css,:param:cn cap mode parametrs
Supported modes	_2dInterface

Nonlinear elasto-plastic material model for concrete plates and shells - Concrete2

Nonlinear elasto-plastic material model with hardening. Takes into account uniaxial stress + transverse shear in concrete layers with transverse stirrups. Can be used only for 2d plates and shells with layered cross section and together with explicit integration method (stiffness matrix is not provided). The model description and parameters are summarized in Table 21.

Table 21 Nonlinear elasto-plastic material model for concrete - summary.

Description	Nonlinear elasto-plastic material model for concrete plates and shells
Record Format	Concrete2 num (in) d (rn) E (rn) n (rn) SCCC (rn) SCCT (rn) EPP (rn) EPU (rn) EOPU (rn) EOPP (rn) SHEARTOL (rn) IS_PLASTIC_FLOW (in) IFAD (in) STIRR_E (rn) STIRR_Ft (rn) STIRR_A (rn) STIRR_TOL (rn) STIRR_EREF (rn) STIRR_LAMBDA (rn)
Parameters	num material model number d material density E Young modulus n Poisson ratio SCCC pressure strength SCCT tension strength EPP threshold effective plastic strain for softening in compression EPU ultimate eff. plastic strain EOPP threshold volumetric plastic strain for softening in tension EOPU ultimate volumetric plastic strain SHEARTOL threshold value of the relative shear deformation (psi**2/eef) at which shear is considered in layers. For lower relative shear deformations the transverse shear remains elastic decoupled from bending. default value SHEARTOL = 0.01 IS_PLASTIC_FLOW indicates that plastic flow (not deformation theory) is used in pressure IFAD State variables will not be updated, otherwise update state variables STIRR_E Young modulus of stirrups STIRR_R stirrups uniaxial strength = elastic limit STIRR_A stirrups area/unit length (beam) or /unit area (shell) STIRR_TOL stirrups tolerance of equilibrium in the z direction (=0 no iteration) STIRR_EREF stirrups reference strain rate for Peryzna’s material STIRR_LAMBDA coefficient for that material (stirrups) SHTIRR_H isotropic hardening factor for stirrups
Supported modes	3dShellLayer, 2dPlateLayer