Constitutive modeling

To handle large deformations it is possible to use the Kirchhoff stress and logarithmic strains as a basis, and use small strain constitutive models with an additive strain split. What this means is that out of the box any strain->stress mapping that works with small strain will work in a large deformation context. Plasticity may be adapted using an exponential mapping.

To model stress rate equations, we can't use this method. To solve this an objective stress rate must be used.

We can update our logarithmic strain by first taking calculating the velocity gradient \(L\)

\[ L = \nabla v \]

The stretch rate can be found from either from this, or directly from the grid:

\[ D = \frac{1}{2}(L + L^T) \]\[ D = \frac{1}{2}(\nabla v + (\nabla v)^T) \]

We can then find our deformation gradient update from this update.

\[ \Delta F = I + L \]\[ F_{n+1} = \Delta F F_n \]

By taking our engineering logarithmic strain, we can calculate the left Cauchy-Green strain tensor \(b_n\):

\[ \varepsilon_{n+1} = exp(2b_{n+1}) \]

We can then update the left Cauchy-Green strain by the formula:

\[ b_{n+1} = \Delta F b_n (\Delta F)^T \]

From our updated strain \(b_{n+1}\), we can re-derive our engineering logarithmic strain measures.

\[ \varepsilon_{n+1} = \frac{1}{2}\ln(b_{n+1}) \]

This is fine for constitutive modelling of hyperelastic materials, and plastic materials - but if you want to write constitutive models in terms of a stress rate this won't work.

If you want to define a viscoelastic material - where the stress at every step is not defined entirely by the strain we need an integral form of stress update.

If we have a constitutive function \(C\) of the form:

\[ \dot\sigma = C(\sigma,\dot\varepsilon) \]

We will find that as the material moves and rotates through the world we will get stresses pointing in the wrong direction. To solve this we use a pair of so-called "objective" measures of stress and strain \(\mathring{\sigma},\mathring{\varepsilon}\).

There are many objective stress rates, with varying degrees of complexity. A very good choice of objective stress rate keeping with our logarithmic strains is the corrotational logarithmic stress rate, as it turns out the strain-rate measure \(\mathring{\varepsilon}\) is the natural Eulieran stretch measure \(D\)!.

\[ \mathring{\sigma} = \dot{\sigma} + \Omega^{log} \sigma - \sigma\Omega^{log}\ \]

Where we use the tensor factor \(\Omega^{log}\) which somehow captures the rotation of the stress \(\sigma\) over each timestep. For posterity the form of this tensor is defined by the eigenprojectections and the stretch-rate:

\[ P^i = v_i \otimes v_i\\ \]\[ \Omega^{log}_{ij} = w_{ij} + \sum^{n}_{A,B=1,A\neq B}{\left(\frac{\lambda_A + \lambda_B}{\lambda_A - \lambda_B} + \frac{2}{\ln{\lambda_A} - \ln{\lambda_B}}\right)} P^A_{ik} D_{kl} P^B_{lm} \]

By utilising this adjustments factor we may still keep our rate-based small-strain constitutive formulations in a large deformation context but understanding that we are using logarithmic strain under the hood.

Models

From our logarithmic strain, we may implement small strain theory constitutive models (additive strain split) in a large deformation system with ease. Each constitutive model is implemented as a different class, found in the src/particles file.

Linear elasticity

The simplest model we may implement is isotropic linear elasticity with material non-linearity. We state there is a linear relationship between the work conjugate pair of logarithmic strain and Kirchhoff stress.

\[ \tau = [D^e] \varepsilon \]

Material points with elastic behaviour are derived from the particle-elastic class. The class contains \(E\) and \(\nu\), and \([D^e]\) variables. When either of \(E\) or \(\nu\) are updated, the elastic matrix is updated alongside it.

This is made easy by CLOS's ability to hook into setting parameters:

(defun update-elastic-matrix (particle)
  (with-accessors ((de mp-elastic-matrix)
                    (E  mp-E)
                    (nu mp-nu))
      particle
    (setf de (cl-mpm/constitutive::linear-elastic-matrix E nu))))

(defmethod (setf mp-E) :after (value (p particle-elastic))
  (update-elastic-matrix p))

We define a function (defun) for updating the elastic matrix, and a generic method (defmethod) that runs after a material points mp-E slot has been set. On the caller side, this happens seamlessly.

(setf (mp-e mp) 10d0) -> [D^e] updated
(setf (mp-damage mp) 1d0) -> No special behaviour

Linear viscoelasticity

We may use the corotational rate of Kirchhoff stress to write a simple compressible viscoelastic model.

\[ \dot{\tau} = G \mathring{\varepsilon}^D - \frac{G}{\eta} \tau^D \]

Where we have extended the Maxwell viscoelastic model to only apply on the deviatoric stress. This is then implemented with a non-objective integration algorithm as for an explicit-dynamic code the rotation per timestep is small enough to be accurate.

\[ \tau_{n+1} = \tau_n + (G \Delta t \mathring{\varepsilon}^D - \Delta t \frac{G}{\eta} \tau^D + \Delta t \tau \Omega^{log} - \Delta t \Omega^{log}\tau) \]

Material points with viscoelastic behaviour are derived from the particle-viscoelastic class - itself inheriting from the particle-elastic class.

Continuum damage

Modeling fracture numerically is a very complex and finicky topic, and often comes with many issues that are hard to overcome.

One of the most powerful and easy to understand frameworks is continuum damage mechanics, where we attempt to model fracture in some bulk averaged way.

Instead of trying to represent every micro-crack in a continuum discreetly, we lump them together into a state parameter damage \(d\). The materials state of decay can be represented with a value of 0 (undamaged) to 1 (completely damaged). One way of thinking of this parameter is the ratio of area of the material that has now become a crack/void.

\[ d = \frac{A_D}{A} \]

This parameter then effects the average properties of a material, such as through the constitutive model of a damaging elastic material.

Undamaged:

\[ \tau = [D^e] \varepsilon \]

Damaged:

\[ \overline{\tau} = (1-d)[D^e] \varepsilon \]

Where we represent the damaged stress with an over-bar.

To evolve damage, we can prescribe a relationship between an unbounded history dependent variable \(\kappa\) and the bounded damage \(d\).

\[ d(\kappa) = 1 - e^{\frac{\kappa - \varepsilon_0}{\varepsilon_f - \varepsilon_0}} \]

Where \(\varepsilon_0\) is the initiation strain - directly limited to the failure stress, and \(\varepsilon_f\) is the failure strain and related to the softening.

Damage evolution

For a quasi-static simulation we may evolve \(\kappa\) as the maximum of our damage criteria \(\chi\) over time:

\[ \kappa_{n+1} = max[\chi,\kappa_{n}] \]

For dynamic fracture simulations with lots of transient stresses we can take a delayed damage approach, and use a similar approach to viscoplastic regularization.

We introduce a characteristic time \(\tau\) that slows damage updating, leading to a PDE that must be time-integrated.

\[ \dot{\kappa} = \frac{\langle\chi - \kappa\rangle}{\tau} \]

Damage criteria

The damage criteria determines how the damage evolves.

Maximum principal stress

The maximum principal stress is the simplest damage criteria:

\[ \chi = \sigma_I \]