19 Constitutive Models
19.1 Abstract Interface for Materials
Material models need to implement a number of different functionalities to interface with different parts of the codebase.
19.1.1 Scalar Stiffness Estimation
The general form of material stiffness is a rank-four tensor that may depend nonlinearly on the unknown solution. This would add significant complexity, computational expense, and nonlinearity in situations where a term of a formulation only needs to broadly scale with material stiffness up to some constant. Examples of such situations are penalty regularizations of Dirichlet boundary conditions, material interfaces, and mechanical contact, or stabilization terms to relieve locking or otherwise improve solution quality. To be compatible with such penalty and stabilized methods, materials must implement an interface that returns a constant (i.e., solution-independent) estimate of material stiffness. We currently assume that this estimate takes the form of an estimated shear modulus and estimated bulk modulus . These estimates can be combined into an estimate of Young’s modulus using a standard conversion formula: For linear isotropic elasticity, these estimates have the obvious choices of the nominal shear and bulk moduli, i.e., and . For more general isotropic materials, it is usually sufficient to linearize about the undeformed reference configuration, which will result in an isotropic stiffness with unique shear and bulk moduli. Some general formulas for estimating effective bulk and shear moduli of isotropic hyperelastic models are given in [80], where they are used to define pressure stabilization parameters.
Evaluation of stiffness estimates in the reference configuration reduces nonlinearity of stabilization and/or penalty terms, but corresponds to an assumption that stiffness does not dramatically change in deformed configurations. Ignoring significant softening, e.g., loss of shear stiffness during perfectly-plastic yielding, is largely benign when estimating penalty parameters, since softening tends to be localized such that the (locally) overestimated stiffness would remain appropriate elsewhere and not be a limiting influence on critical time step size. However, the ideal design of these estimates for materials that are highly anisotropic or whose stiffness increases significantly with deformation (e.g., Fung-type materials in soft tissue mechanics) remains an open problem.
19.2 Isotropic neo-Hookean material model
The isotropic neo-Hookean material is a two parameter hyperelastic material model that is identical to the linear elastic model in the small deformation limit.
19.3 Elastoplasticity
The following elastoplastic material model is an isotropic, nonlinear, rate-independent model, based largely on the model in Chapter 9 of Simo & Hughes’ book [101] (see also [100], [99], [31]). It relies on the concept of an intermediate stress-free configuration whose kinematics are described through a multiplicative decomposition of the deformation gradient into plastic and elastic parts:
As is standard for metals, it is further assumed that the plastic deformation is isochoric, i.e., .
With each incremental increase in load, a trial state is defined. The trial state assumes no additional plastic flow since the previous increment.
- This assumption is tested using the yield criterion, which defines a yield surface in stress-space:
If the trial state lies inside the yield surface, the assumption of no plastic flow holds.
Otherwise, the trial state must be updated in a way that returns it to the closest point on the yield surface; this is called the Return Mapping or Radial Return algorithm [101]. This is performed by approximating the plastic flow rules as finite difference time derivatives and algorithmically integrating them over the pseudo-time step. The result provides updates for the elastic and plastic strain measures and the stress at the current step.
An overview of the equations and algorithms involved are included in the following subsections.
19.3.1 Elastoplastic theory
Elastic-plastic strain decomposition:
The deformation gradient is decomposed into plastic and elastic parts. The elastic left Cauchy-Green tensor is defined as expected: as is the plastic right Cauchy-Green tensor: For the energy density function used here, it is convenient to define volume preserving terms, denoted by a bar:
Hyperelastic constitutive law:
A neo-Hookean energy density function is used in which the energy associated with the volumetric and volume-preserving deformation are separated into two parts: where This leads to the following equation for the Kirchhoff stress:
von Mises yield condition
While the von Mises yield condition is often written in terms of the Kirchhoff stress, , for plasticity models, it is defined here as depending on the Cauchy stress, to better reflect yield stresses and hardening curves obtained by experiment: where is the yield stress and is the hardening curve as a function of the equvialent plastic strain, .
Flow rule
Evolution of the hardening parameter
Kuhn-Tucker loading/unloading conditions
Consistency condition
19.3.2 Radial return algorithm
The radial return algorithm specifies the update for stress and strain variables at each load increment in the case of plastic loading. It is developed by expressing the flow rule with a backward difference approximation of the time derivative: where has been rewritten as . In general, such a time discretization will lead to a nonlinear system of equations that must be solved iteratively, where the terms at step are known and those at are unknown. However, this particular flow rule and strain energy function allow a closed form update to be derived, expressed in terms of “trial” variables that depend only on the previous time step. The resulting updates are shown in the following algorithm.
◾ Radial return algorithm
1: //Update the current configuration
2:
3:
4:
5:
6: //Compute the damage/elastic predictor
7:
8:
9:
10:
11: //Check for plastic loading
12:
13: if then
14: Set
15:
16: return , ,
17: else
18: Proceed to return mapping
19: end if
20: //Return-mapping
21:
22: Solve:
23:
24:
25:
26: //Update the Kirchhoff stress
27:
28: //Update the intermediate configuration
29: , where is solved for iteratively to satisfy the equation
30: return , ,
Remark: The update (see, e.g., [101]) does not maintain the relation that is implied by isochoric plastic deformation. Instead, the following update is used: where is solved for iteratively to satisfy the equation . Since does not depend on the trace of , this update does not affect the elastoplastic moduli.
19.3.2.1 Precision improvement
To improve the precision of the plasticity model, we store the elastic left Cauchy-Green deformation tensor that does not include the identity tensor, denoted herein as . Then, we update the trial deformation that also does not include , based on from the previous step.
The deformation gradient is where is the displacement gradient, and is the identity matrix. Also define
The inverse of can be expressed as where
Split the elastic left Cauchy-Green deformation tensor into the identity and the part without the identity:
At time step , the trial deformation is obtained as where which is equivalent to where
With these, we get where
19.3.3 Consistent spatial elastoplastic moduli
The spatial moduli tensor can be expressed through the following relationship (using indicial notation): The consistent algorithmic elastoplastic moduli are derived by substituting the algorithmically updated stress for in the preceding equation and differentiating: The resulting moduli are included below.
Spatial elasticity tensor for hyperelastic model
Scaling factors
Consistent algorithmic moduli
The final term of the moduli, , comes from basing the yield condition on the Cauchy stress instead of the Kirchhoff stress.
Furthermore, both the terms and in the consistent tangent do not have major symmetry. This is because, for the strain energy and yield functions used here, the flow rule is not associative. If required to be symmetric due to, for example, an iterative solver, these terms may be symmetrized, but the nonlinear convergence rate will be negatively affected.
Finally, note that does not reduce to when is zero. This has implications for the first iteration of each nonlinear solve, when is initialized as . On that iteration, will be less than or (nearly) equal to zero, which will lead to the elastic moduli being used for the nonlinear update. However, if the material will indeed undergo plastic flow at that timestep, using the elastic moduli instead of the elastoplastic moduli on the first iteration can have a significantly negative effect on convergence. To mitigate this problem, a check is made on the first nonlinear iteration such that if is nearly zero (or if plastic flow occurred on the previous step), then it is temporarily assumed that plastic flow will continue on the current step and the elastoplastic moduli are used. For subsequent iterations of each nonlinear solve, the extra check is not needed and the elastoplastic moduli are only used if .
19.3.4 Material failure
A material failure happens when is met, for which a material failure indicator is introduced as follows: is a stress triaxility-dependent failure plastic strain. In other words, the failure occurs when reaches , given a constant stress triaxiality .
Metals do not typically experience failure in the compression-dominent regime but do present a complex failure locus in the shear-dominant and tension-dominant regimes. A simplified failure locus can take the following form. The cutoff triaxiality is a meterial specific parameter and typically around -1/3.
When the material failure is initiated, the stress is set to zero. Numerically, the stress at an integration point with is gradually degraded over a range , where is a small value.
19.4 Thermal expansion
The most general kinematic framework for elasticity with thermal expansion is a multiplicative decomposition of the deformation gradient, where is the deformation gradient due to an elastic deformation superposed on a thermal expansion with deformation gradient where is the thermal stretch ratio. The definition of the thermal stretch ratio depends on what strain measure is used to interpret thermal strains. If thermal strain is understood in the sense of engineering strain, i.e., the ratio of length change to initial length, then the thermal stretch ratio is given by where is the (linear) thermal expansion coefficient, and is the temperature change from some reference temperature about which is calibrated. If depends on temperature, it is important to note that this definition assumes a secant coefficient of thermal expansion, whose definition depends on the choice of reference temperature. If thermal strain is instead understood in the sense of logarithmic (also known as Hencky) strain, then the thermal stretch ratio is given by This is the default behavior, because a logarithmic-strain interpretation of thermal strain arises more naturally when computing the secant coefficient of thermal expansion via integration of a given tangent coefficient. Using the logarithmic-strain interpretation, we can compute from a given tangent coefficient using the formula where is the reference temperature and is the tangent coefficient of thermal expansion. In particular, a constant tangent coefficient implies a constant secant coefficient under the logarithmic-strain interpretation, but this is no longer true under the engineering-strain interpretation. In the small-strain setting of linearized kinematics, there is no distinction between the engineering-strain and logarithmic-strain definitions of the thermal stretch ratio, because for .
To account for thermal expansion’s effects on the equilibrium of an elastic solid, the elastic potential energy in an existing constitutive model can be modified by replacing the full deformation gradient with the elastic one in the definition of the energy functional. Stresses and material moduli then follow via differentiation in the usual way. In elastoplastic analyses, the multiplicative decomposition becomes , i.e., the elastoplastic deformation is applied on top of thermal expansion, taking the place of in the preceding discussion of thermoelasticity.
19.4.1 Thermal stresses in linear elasticity
If we assume that both thermal expansion and stretches are small, then we can approximate the multiplicative decomposition by an additive decomposition of the Green–Lagrange strain, with defined from in the usual way, and Further assuming small rotations, we can make the usual substitution of the small engineering strain for Green–Lagrange strain to get Plugging this into a linear-elastic isotropic constitutive model, we get the following expression for stress: where refers to the stress understood as a linear function of engineering strain and we formally convert the thermal strain to a "thermal stress" given by in which the "thermal pressure" is given by Thus, in the weak form of the problem, thermal strain can be added to an existing implementation of linear elasticity by adding a term to the weak form of the problem, in a manner closely analogous to a body force.
19.4.2 Approximate thermal pressure in finite strain theory
The formal conversion of a thermal expansion strain into a thermal stress load relies on the linearity of the mapping , which is no longer applicable in the cases of finite deformation and/or a nonlinear constitutive model. However, in the case where volume changes are small, and the material’s pressure response is a function of the deformation gradient Jacobian , i.e., where is the contribution of volume change to the elastic energy density, we can still approximate the effect of thermal expansion through a thermal pressure of the form where refers to the effective tangent bulk modulus evaluated in the reference configuration. The effective bulk modulus can be defined as in hyperelastic models where the contribution of volume change to energy density is entirely determined by . Subtracting from the diagonal of the Cauchy stress provides a good approximation of thermal expansion at finite strains when both and are close to one. It does not impose any limits on rotations or shear deformation. This approximation simplifies the implementation of unsteady displacement/pressure formulations where the pressure obeys a rate equation. A direct implementation of thermal pre-strain in this setting would require access to time derivative information about the thermal deformation along with solution of an auxiliary problem to initialize the pressure. Displacement/pressure formulations are primarily used in cases where volume changes are small, thus satisfying the kinematic assumption behind the approximation. One may also choose to use this approximation in static displacement/pressure formulations, for consistent behavior of thermal expansion across multi-procedure analyses involving both static and dynamic procedures. Because this formulation is derived from an assumption of small thermal expansion strains, it does not make any distinction between the engineering-strain and logarithmic-strain interpretations of thermal stretch.
19.5 Viscoelasticity
Viscoelasticity is conceptualized as an overstress, i.e., a stress contribution added on top of the stress due to a given equilibrium constitutive model that governs the steady-state stress-strain behavior. In practice, one often uses a sum of multiple such overstresses with different parameters, e.g., a Prony series. In this section, we discuss the formulation and discretization from [91] of a single viscoelastic overstress term parameterized by a shear modulus and timescale, with the understanding that several overstresses can be superposed.
19.5.1 Viscoelastic theory
Viscoelasticity relies on a multiplicative decomposition of the deformation gradient, where is an inelastic deformation that lags behind the full deformation, catching up over some timescale if the deformation is held fixed. When the system is not in an equilibrium with , we have a non-identity elastic part of the deformation given by , which leads to a viscous stress that is computed by plugging its corresponding elastic right Cauchy-Green tensor into a hyperelastic potential. In particular we use an isotropic neo-Hookean elastic model to generate the viscous stress, which can be expressed as the following reference-configuration second Piola-Kirchhoff tensor: where is the viscous shear modulus, is the inelastic right Cauchy-Green tensor, and the operators and are defined as and for an arbitrary rank-two tensor . Note that the use of a neo-Hookean potential to generate does not prevent this viscous stress from being applied as an overstress to an equilibrium stress response coming from an unrelated (i.e., non-neo-Hookean) constitutive model. Meanwhile, is a local history variable that evolves according to where is the viscous timescale.
19.5.2 Time discretization
If the inelastic right Cauchy-Green tensor’s evolution equation is integrated with a modified backward Euler method that enforces a volume-preserving inelastic deformation as a correction at each step, we obtain the following explict formula for at time level , in terms of its previous value and the -level full right Cauchy-Green tensor: This is extended to time integrators other than backward Euler by letting and be whatever time levels the displacement is evaluated at in consecutive time steps, e.g., in a midpoint rule or and in a generalized- method. In other words, the backward Euler method is applied between displacement evaluations, similar to a standard implementation of the radial return mapping for plasticity.