# Mixed Elasticity

## Mixed Linear

The boundary-value problem (Strong form) for constitutive equation {eq}`general-mixed-linear-constitutive` may be stated as follows: Given body force $\rho \bm g$, Dirichlet boundary $\bar{\bm u}$ and applied traction $\bar{\bm t}$, find the displacement and pressure-like variables $(\bm u, p) \in \mathcal{V} \times \mathcal{Q} $ (here $\mathcal{V} = H^1(\Omega), \mathcal{Q} = L^2(\Omega) $ ), such that:

$$
  \begin{aligned}
    -\nabla \cdot \bm \sigma - \rho \bm g   &= 0, \qquad \text{in $\Omega$} \\
    -\nabla\cdot\bm u - \frac{p}{\bulk - \bulk_p} &= 0, \qquad \text{in $\Omega$} \\
    \bm{u} &= \bar{\bm u}, \quad    \text{on $\partial \Omega^{D}$} \\
    \bm \sigma \cdot \bm n &= \bar{\bm t}. \qquad \text{on $\partial \Omega^{N}$}
  \end{aligned}
$$ (mixed-linear-strong-form)

with $\bm n$ be the unit normal on the boundary and its weak formulation as:

$$
  \begin{aligned}
    \int_{\Omega}{ \nabla \bm{v} \tcolon \bm{\sigma}} \, dv
    - \int_{\partial \Omega}{\bm{v} \cdot \bar{\bm t}} \, da
    - \int_{\Omega}{\bm{v} \cdot \rho \bm{g}} \, dv
    &= 0, \quad \forall \bm v \in \mathcal V \\
    \int_{\Omega} q \left( -\nabla\cdot\bm u - \frac{p}{\bulk - \bulk_p} \right) \, dv & = 0. \quad \forall q \in \mathcal Q
  \end{aligned}
$$ (mixed-linear-weak-form)

## Mixed Hyperelasticity, initial configuration

We can state the mixed hyperelastic strong form as follow: Given body force $\rho_0 \bm g$, Dirichlet boundary $\bar{\bm u}$ and applied traction $\bar{\bm t}$, find the displacement and pressure-like variables $(\bm u, p) \in \mathcal{V} \times \mathcal{Q} $, such that:

$$
  \begin{aligned}
    -\nabla_X \cdot \bm P - \rho_0 \bm g   &= 0. \qquad \text{in $\Omega_0$} \\
    -\frac{\partial V}{\partial J} - \frac{p}{\bulk - \bulk_p} &= 0. \qquad \text{in $\Omega$} \\
    \bm{u} &= \bar{\bm u}, \quad    \text{on $\partial \Omega^{D}_0$} \\
    \bm P \cdot \bm N &= \bar{\bm t}. \qquad \text{on $\partial \Omega^{N}_0$}
  \end{aligned}
$$ (mixed-hyperelastic-strong-form)

and by multiplying the strong form {eq}`mixed-hyperelastic-strong-form` by test functions $(\bm{v}, q)$ and  integrate by parts, we can obtain the weak form for finite-strain incompressible hyperelasticity as:
find $(\bm{u}, p) \in \mathcal{V} \times \mathcal{Q} \subset H^1 \left( \Omega_0 \right) \times L^2 \left( \Omega_0 \right)$ such that

$$
  \begin{aligned}
    R_u(\bm u, p) &\coloneqq \int_{\Omega_0}{\nabla_X \bm{v} \tcolon \bm{P}} \, dV
     - \int_{\Omega_0}{\bm{v} \cdot \rho_0 \bm{g}} \, dV
     - \int_{\partial \Omega_0}{\bm{v} \cdot \bar{\bm t}} \, dA = 0, \quad \forall \bm{v} \in \mathcal{V}, \\
    R_p(\bm u, p) &\coloneqq \int_{\Omega_0} q \cdot L J \, dV = 0, \quad \forall q\in \mathcal{Q},
  \end{aligned}
$$ (mixed-hyperelastic-weak-form-initial)

where $L =  -V' - \frac{p}{\bulk - \bulk_p} $. This equation contains material/constitutive nonlinearities in defining $\bm{S}(\bm{E})$, as well as geometric nonlinearities through $\bm{P} = \bm{F}\, \bm{S}$, $\bm{E}(\bm{F})$, and the body force $\bm{g}$, which must be pulled back from the current configuration to the initial configuration.

### Newton linearization

As explained in {ref}`Newton-linearization`, we need the Jacobian form of the {eq}`mixed-hyperelastic-weak-form-initial` which is:
find $(\diff \bm{u}, \diff p) \in \mathcal{V} \times \mathcal{Q}$ such that

$$
  \begin{aligned}
    \int_{\Omega_0} \nabla_X \bm{v} \tcolon \diff \bm{P} dV &= -R_u, \quad \forall \bm{v} \in \mathcal{V}, \\
    \int_{\Omega_0} q \cdot \left( \diff L J + L \diff J \right)dV &= -R_p, \quad \forall q \in \mathcal{Q}
  \end{aligned}
$$ (mixed-jacobian-weak-form-initial)

where

$$
  \begin{aligned}
    \diff \bm{P} &= \diff \bm{F}\, \bm{S} + \bm{F} \diff \bm{S}, \\
    \diff L J + L \diff J &= J \frac{\partial L}{\partial \bm{E}} \tcolon \diff \bm{E} + J \frac{\partial L}{\partial p} \diff p + L \frac{\partial J}{\partial \bm{E}} \tcolon \diff \bm{E}
  \end{aligned}
$$ (mixed-eq-diff-P)

with $\diff \bm{F}$ and $\diff \bm E$ defined in {eq}`dE-and-dF`. The linearization of the second equation of {eq}`mixed-eq-diff-P` is

$$
  \begin{aligned}
    \diff L J + L \diff J &= \left( -V^{''} \diff J - \frac{\diff p}{\bulk - \bulk_p} \right) J + L J (\bm{C}^{-1} \tcolon \diff \bm{E}), \\
    &= \left( - J^2 \, V^{''} - J \, V' - J \frac{p}{\bulk - \bulk_p} \right) \bm{C}^{-1} \tcolon \diff \bm{E} - J \frac{\diff p}{\bulk - \bulk_p}
  \end{aligned}
$$ (mixed-eq-diff-P-2)


The linearization of the second Piola-Kirchhoff stress tensor, $\diff \bm{S}$, depends upon the material model and derived below for different hyperelastic materials.

:::{dropdown} Deriving $\diff\bm{S}$ for isochoric Neo-Hookean material

For the Neo-Hookean model {eq}`mixed-neo-hookean-stress`, we derive split

$$
\diff \bm{S} = \underbrace{\frac{\partial \bm{S}_{\text{vol}}}{\partial \bm{E}} \tcolon \diff \bm{E}}_{\diff \bm{S}_{\text{vol}}^u} + \underbrace{\frac{\partial \bm{S}_{\text{iso}}}{\partial \bm{E}} \tcolon \diff \bm{E}}_{\diff \bm{S}_{\text{iso}}} + \underbrace{\frac{\partial \bm{S}_{\text{vol}}}{\partial p} \diff p}_{\diff \bm{S}_{\text{vol}}^p},
$$ (dS-split-neo-hookean)

then,

$$
  \begin{aligned}
    \diff \bm{S}_{\text{vol}}^u &= \left( \bulk_p \diff J \, V' + \bulk_p J V^{''} \diff J - p \diff J \right) \bm{C}^{-1} + \left(\bulk_p J \, V' - p \, J \right) \diff \bm{C}^{-1}, \\
    &= \left( \bulk_p J^2 V^{''} + \bulk_p J \, V' - p J \right) (\bm{C}^{-1} \tcolon \diff \bm{E}) \bm{C}^{-1} + \left(\bulk_p J \, V' - p \, J \right) \diff \bm{C}^{-1}
  \end{aligned}
$$ (dS-vol-neo-hookean)

$$
\diff \bm{S}_{\text{vol}}^p = -J \bm{C}^{-1} \diff p,
$$ (dS-p-neo-hookean)

and

$$
  \begin{aligned}
    \diff \bm{S}_{\text{iso}} &= -\frac{2}{3}\mu J^{-2/3} (\bm{C}^{-1} \tcolon \diff \bm{E}) \left(\bm{I} - \frac{1}{3} \mathbb{{I}_1} \bm{C}^{-1} \right) \\
     &- \frac{1}{3} \mu J^{-2/3} \left( 2 \trace \diff \bm{E} \, \bm{C}^{-1} + \mathbb{{I}_1} \diff \bm{C}^{-1} \right), \\
    &= -\frac{4}{3}\mu J^{-2/3} (\bm{C}^{-1} \tcolon \diff \bm{E}) \bm{C}^{-1} \bm{E}_{\text{dev}} - \frac{1}{3} \mu J^{-2/3} \left( 2 \trace \diff \bm{E} \, \bm{C}^{-1} + \mathbb{{I}_1} \diff \bm{C}^{-1} \right)
  \end{aligned}
$$ (dS-iso-neo-hookean)

where

$$
\diff \bm{C}^{-1} = \frac{\partial \bm{C}^{-1}}{\partial \bm{E}} \tcolon \diff \bm{E} = -2 \bm{C}^{-1} \diff \bm{E} \, \bm{C}^{-1} .
$$


:::{note}
If we use single field isochoric model {eq}`isochoric-neo-hookean-stress` the linearization of the volumetric stress becomes

$$
\diff \bm{S}_{\text{vol}} = \left( \bulk J^2 V^{''} + \bulk J \, V'\right) (\bm{C}^{-1} \tcolon \diff \bm{E}) \bm{C}^{-1} + \bulk J \, V' \diff \bm{C}^{-1}.
$$ (dS-vol-isochoric-neo-hookean)

$\diff \bm{S}_{\text{vol}}^p = 0$ and we only need to solve the first equation of {eq}`mixed-hyperelastic-weak-form-initial`.

:::

:::


:::{dropdown} Deriving $\diff\bm{S}$ for isochoric Mooney-Rivlin material

For the Mooney-Rivlin model {eq}`mixed-mooney-rivlin-stress`, we derive split

$$
\diff \bm{S} = \underbrace{\frac{\partial \bm{S}_{\text{vol}}}{\partial \bm{E}} \tcolon \diff \bm{E}}_{\diff \bm{S}_{\text{vol}}^u} + \underbrace{\frac{\partial \bm{S}_{\text{iso}}}{\partial \bm{E}} \tcolon \diff \bm{E}}_{\diff \bm{S}_{\text{iso}}} + \underbrace{\frac{\partial \bm{S}_{\text{vol}}}{\partial p} \diff p}_{\diff \bm{S}_{\text{vol}}^p},
$$ (dS-split-mooney-rivlin)

where $\diff \bm{S}_{\text{vol}}^p$ and $\diff \bm{S}_{\text{vol}}^u$ are the same as {eq}`dS-p-neo-hookean`, {eq}`dS-vol-neo-hookean` (or {eq}`dS-vol-isochoric-neo-hookean` if we use single field stress {eq}`isochoric-mooney-rivlin-stress`) and the isochoric part is

$$
  \begin{aligned}
    \diff \bm{S}_{\text{iso}} &= -\frac{4}{3}(\bm{C}^{-1} \tcolon \diff \bm{E}) \left( \mu_1 J^{-2/3} + 4 \mu_2 J^{-4/3} \right) \bm{C}^{-1} \bm{E}_{\text{dev}}\\
     &- \frac{1}{3} \left( \mu_1 J^{-2/3} + 2 \mu_2 J^{-4/3} \right) \left( 2 \mathbb{{I}_1}(\diff \bm{E})\, \bm{C}^{-1} + \mathbb{{I}_1} \diff \bm{C}^{-1} \right) \\
    &- \frac{8}{3}(\bm{C}^{-1} \tcolon \diff \bm{E}) \mu_2 J^{-4/3}  \left( \mathbb{{I}_1}(\bm{E})\bm{I} - \bm{E} \right) \\
    &+2 \mu_2 J^{-4/3} \left( \mathbb{{I}_1}(\diff \bm{E})\bm{I} - \diff \bm{E} \right) \\
    & + (c_1 + c_2) \bm{C}^{-1} - \frac{4}{3} \mu_2 J^{-4/3} \left( \mathbb{{I}_1}(\bm{E}) + 2 \mathbb{{I}_2}(\bm{E})\right) \diff \bm{C}^{-1}
  \end{aligned}
$$ (dS-iso-mooney-rivlin)

where

$$
\begin{aligned}
  c_1 &= \frac{16}{9} \mu_2 J^{-4/3} (\bm{C}^{-1} \tcolon \diff \bm{E}) \left( \mathbb{{I}_1}(\bm{E}) + 2 \mathbb{{I}_2}(\bm{E}) \right). \\
  c_2 &= -\frac{4}{3} \mu_2 J^{-4/3} \left( \mathbb{{I}_1}(\diff \bm{E}) + 2 \mathbb{{I}_1}(\bm{E}) \mathbb{{I}_1}(\diff \bm{E}) - 2 \bm{E} \tcolon \diff \bm{E} \right).
\end{aligned}
$$ (coeff)

:::

:::{dropdown} Deriving $\diff\bm{S}$ for isochoric Ogden material
Similar to the Neo-Hookean model we have

$$
\diff \bm{S} = \underbrace{\frac{\partial \bm{S}_{\text{vol}}}{\partial \bm{E}} \tcolon \diff \bm{E}}_{\diff \bm{S}_{\text{vol}}^u} + \underbrace{\frac{\partial \bm{S}_{\text{iso}}}{\partial \bm{E}} \tcolon \diff \bm{E}}_{\diff \bm{S}_{\text{iso}}} + \underbrace{\frac{\partial \bm{S}_{\text{vol}}}{\partial p} \diff p}_{\diff \bm{S}_{\text{vol}}^p},
$$ (dS-split-ogden)

where $\diff \bm{S}_{\text{vol}}^p$ and $\diff \bm{S}_{\text{vol}}^u$ are similar to volumetric derivative of Neo-Hookean and Mooney-Rivlin models and the isochoric part is

$$
  \begin{aligned}
  \diff \bm{S}_{\text{iso}} &= \sum_{i=1}^3 \diff S_i^{iso} \hat{\bm{N}_i} \hat{\bm{N}_i}^T + S_i^{iso} \left( \diff \hat{\bm{N}_i} \hat{\bm{N}_i}^T + \hat{\bm{N}_i} \diff \hat{\bm{N}_i}^T\right)
  \end{aligned}
$$ (dS-iso-ogden)

For example computing $\diff S_1^{iso}$ from {eq}`mixed-ogden-stress-iso-stable1` gives

$$
  \begin{aligned}
  \diff S_1^{iso} &= -\frac{2\diff \beta_1}{\beta_1^3} \sum_{j=1}^N \frac{m_j}{3} \left[ 2\operatorname{\tt expm1}(\alpha_j \ell_1) - \operatorname{\tt expm1}(\alpha_j \ell_2) - \operatorname{\tt expm1}(\alpha_j \ell_3) \right] J^{-\alpha_j/3} \\
  &+ \frac{1}{\beta_1^2} \sum_{j=1}^N \frac{m_j \alpha_j}{3} \left[ 2 \diff \ell_1 \exp(\alpha_j \ell_1) - \diff \ell_2 \exp(\alpha_j \ell_2) - \diff \ell_3 \exp(\alpha_j \ell_3) \right] J^{-\alpha_j/3} \\
  &- \frac{1}{\beta_1^2} \sum_{j=1}^N \frac{m_j \alpha_j}{9} \left[ 2 \operatorname{\tt expm1}(\alpha_j \ell_1) - \operatorname{\tt expm1}(\alpha_j \ell_2) - \operatorname{\tt expm1}(\alpha_j \ell_3) \right] J^{-\alpha_j/3} (\bm{C}^{-1} \tcolon \diff \bm{E})
  \end{aligned}
$$ (mixed-ogden-stress-iso-linearization)

To compute $\diff \beta_i$ we differentiate $\bm{C} = \sum_{i=1}^3 \beta_i^2 \hat{\bm{N}_i} \hat{\bm{N}_i}^T$ as

$$
\diff \bm{C} = \sum_{i=1}^3 2 \beta_i \diff \beta_i \hat{\bm{N}_i} \hat{\bm{N}_i}^T + \beta_i^2 \left(\diff \hat{\bm{N}_i} \hat{\bm{N}_i}^T + \hat{\bm{N}_i} \diff \hat{\bm{N}_i}^T \right)
$$ (dC-ogden)

and used the fact that the eigenvectors are orthonormal i.e., $\langle \hat{\bm{N}_i}, \hat{\bm{N}_j} \rangle = \delta_{ij}$ so that $\langle \hat{\bm{N}_i}, \diff \hat{\bm{N}_i} \rangle = 0$, we can multiply {eq}`dC-ogden` from the left and right by $\hat{\bm{N}_i}$

$$
\langle \hat{\bm{N}_i}, \diff \bm{C} \hat{\bm{N}_i} \rangle = 2 \beta_i \diff \beta_i \Rightarrow \diff \beta_i = \frac{1}{\beta_i} \langle \hat{\bm{N}_i}, \diff \bm{E} \hat{\bm{N}_i} \rangle
$$ (dbeta)

By differentiating $\bm{E} \hat{\bm{N}_i} = \beta_i^E \hat{\bm{N}_i}$ we will arrive at

$$
  \diff \bm{E} \hat{\bm{N}_i} + \bm{E} \diff \hat{\bm{N}_i}= \diff \beta_i^E \hat{\bm{N}_i} + \beta_i^E \diff \hat{\bm{N}_i}
$$

taking the inner product of above equation with $\hat{\bm{N}_j}, \, j\neq i$ simplifies to

$$
\begin{aligned}
  &\quad \langle \hat{\bm{N}_j}, \diff \bm{E} \hat{\bm{N}_i} \rangle + \langle \diff \hat{\bm{N}_i}, \bm{E} \hat{\bm{N}_j} \rangle = \beta_i^E \langle \diff \hat{\bm{N}_i}, \hat{\bm{N}_j} \rangle \\
  &\quad \langle \hat{\bm{N}_j}, \diff \bm{E} \hat{\bm{N}_i} \rangle + \beta_j^E \langle \diff \hat{\bm{N}_i}, \hat{\bm{N}_j} \rangle = \beta_i^E \langle \diff \hat{\bm{N}_i}, \hat{\bm{N}_j} \rangle \\
  &\quad \langle \diff \hat{\bm{N}_i}, \hat{\bm{N}_j} \rangle = \frac{1}{\beta_i^E - \beta_j^E} \langle \diff \bm{E} \hat{\bm{N}_i}, \hat{\bm{N}_j} \rangle \\
  &\quad \diff \hat{\bm{N}_i} = \sum_{j \neq i} \frac{1}{\beta_i^E - \beta_j^E} \langle \diff \bm{E} \hat{\bm{N}_i}, \hat{\bm{N}_j} \rangle \hat{\bm{N}_j}
\end{aligned}
$$

and finally the linearization of $\ell_i = \log \beta_i$ is

$$
  \diff \ell_i = \frac{\diff \beta_i}{\beta_i}.
$$

which complete $\diff\bm{S}$ for Ogden model.
:::

## Mixed Hyperelasticity, current configuration

Similar to what we have shown in {eq}`initial-current-linearize`, we can write the first equation of {eq}`mixed-hyperelastic-weak-form-initial` in the current configuration (see {eq}`hyperelastic-weak-form-current` ) which yields to the Jacobian form

$$
  \begin{aligned}
    \int_{\Omega_0} \nabla_x \bm{v} \tcolon \left( \diff \bm{\tau} - \bm{\tau} \left( \nabla_x \diff \bm{u} \right)^T \right) dV &= -r_u, \quad \forall \bm{v} \in \mathcal{V}, \\
    \int_{\Omega_0} q \cdot \left( \diff L J + L \diff J \right)dV &= -r_p, \quad \forall q \in \mathcal{Q}
  \end{aligned}
$$ (mixed-jacobian-weak-form-current)

where the second equation is similar to {eq}`mixed-jacobian-weak-form-initial` and we can use {eq}`dtau`

$$
  \diff\bm\tau - \bm\tau\left( \nabla_x \diff\bm u \right)^T = \left( \nabla_x \diff\bm u \right) \bm\tau + \bm F \diff\bm S \bm F^T.
$$

to compute the linearization of the first equation in current configuration for different material models.


:::{dropdown} Representation of $\bm F \diff\bm S \bm F^T$ for isochoric Neo-Hookean material
Based on the split {eq}`dS-split-neo-hookean`, we can derive

$$
\bm{F} \diff \bm{S}_{\text{vol}}^u \bm{F}^{T} = \left( \bulk_p J^2 V^{''} + \bulk_p J \, V' - p J \right) \trace (\diff \bm \epsilon) \bm{I} - 2 \left(\bulk_p J \, V' - p \, J \right) \diff \bm \epsilon.
$$ (F-dS-vol-FTranspose-mixed)

$$
\bm{F} \diff \bm{S}_{\text{vol}}^p \bm{F}^{T} = -J \diff p \bm{I},
$$ (F-dS-p-FTranspose)

where $\diff \bm \epsilon$ is defined in {eq}`depsilon-current` and we have used {eq}`Cinv-contr-dE` and {eq}`FinvTranspose-dE-Finv`. The isochoric part can be simplified as

$$
  \begin{aligned}
  \bm{F} \diff \bm{S}_{\text{iso}} \bm{F}^{T} &= -\frac{2}{3} \mu J^{-2/3} \left( 2 \trace (\diff \bm \epsilon) \bm{e}_{\text{dev}} + \trace(\diff \bm e) \bm{I} - \mathbb{{I}_1} \diff \bm \epsilon \right),
  \end{aligned}
$$ (F-dS-iso-FTranspose-neo-hookean)

where $\mathbb{{I}_1} = \trace(\bm{b})$.
:::{note}
If we use single field isochoric model the linearization of the volumetric stress becomes

$$
\bm{F} \diff \bm{S}_{\text{vol}} \bm{F}^{T} = \left( \bulk J^2 V^{''} + \bulk J \, V'\right) \trace (\diff \bm \epsilon) \bm{I} -2  \bulk J \, V' \diff \bm \epsilon.
$$ (F-dS-vol-FTranspose)

$\diff \bm{S}_{\text{vol}}^p = 0$ and we only need to solve the first equation of {eq}`mixed-jacobian-weak-form-current`.

:::

:::

:::{dropdown} Representation of $\bm F \diff\bm S \bm F^T$ for isochoric Mooney-Rivlin material
The $\bm{F} \diff \bm{S}_{\text{vol}}^p \bm{F}^T$ and $\bm{F} \diff \bm{S}_{\text{vol}}^u \bm{F}^T$ in {eq}`dS-split-mooney-rivlin`   are the same as {eq}`F-dS-p-FTranspose`, {eq}`F-dS-vol-FTranspose-mixed` (or {eq}`F-dS-vol-FTranspose` if we use single field formulation ). The isochoric part can be derived as

$$
  \begin{aligned}
    \bm{F} \diff \bm{S}_{\text{iso}} \bm{F}^T &= -\frac{4}{3}(\trace (\diff \bm \epsilon)) \left( \mu_1 J^{-2/3} + 4 \mu_2 J^{-4/3} \right) \bm{e}_{\text{dev}}\\
     &+ \frac{2}{3} \left( \mu_1 J^{-2/3} + 2 \mu_2 J^{-4/3} \right) \left( \mathbb{{I}_1} \diff \bm \epsilon -\mathbb{{I}_1}(\diff \bm{e})\, \bm{I} \right) \\
    &- \frac{8}{3}(\trace (\diff \bm \epsilon)) \mu_2 J^{-4/3}  \left( \mathbb{{I}_1}(\bm{e})\bm{b} - \bm{b} \bm{e} \right) \\
    &+2 \mu_2 J^{-4/3} \left( \mathbb{{I}_1}(\diff \bm{e})\bm{b} - \bm{b} \diff \bm \epsilon \bm{b} \right) \\
    & + (c_1 + c_2) \bm{I} + \frac{8}{3} \mu_2 J^{-4/3} \left( \mathbb{{I}_1}(\bm{e}) + 2 \mathbb{{I}_2}(\bm{e})\right) \diff \bm \epsilon
  \end{aligned}
$$ (F-dS-iso-FTranspose-mooney-rivlin)

where

$$
\begin{aligned}
  c_1 &= \frac{16}{9} \mu_2 J^{-4/3} (\trace (\diff \bm \epsilon)) \left( \mathbb{{I}_1}(\bm{e}) + 2 \mathbb{{I}_2}(\bm{e}) \right). \\
  c_2 &= -\frac{4}{3} \mu_2 J^{-4/3} \left( \mathbb{{I}_1}(\diff \bm{e}) + 2 \mathbb{{I}_1}(\bm{e}) \mathbb{{I}_1}(\diff \bm{e}) - 2 \trace (\bm{b} \bm{e} \diff \bm \epsilon) \right).
\end{aligned}
$$ (coeff-current)

and we have used {eq}`trace-dE`, {eq}`F-dE-FTranspose`, and

$$
\begin{aligned}
  \bm{F} \bm{C}^{-1} \bm{E}_{\text{dev}} \bm{F}^T &= \bm{F}^{-T} \left( \bm{E} - \frac{1}{3} \mathbb{{I}_1}(\bm{E}) \bm{I} \right) \bm{F}^T = \bm{e}_{\text{dev}}. \\
  \bm{F} \bm{E} \bm{F}^T &= \frac{1}{2} \left( \bm{F} \bm{C} \bm{F}^T - \bm{F} \bm{F}^T \right) = \bm{b} \bm{e}. \\
  \diff \bm{E} &= \frac{1}{2} \left( \bm{F}^T \diff \bm{F} \bm{F}^{-1} \bm{F} + \bm{F}^T \bm{F}^{-T} \diff \bm{F}^T \bm{F} \right) = \bm{F}^T \diff \bm \epsilon \bm{F}. \\
  \bm{E} \tcolon \diff \bm{E} &= \trace(\bm{E} \diff \bm{E}) = \trace(\bm{E} \bm{F}^T \diff \epsilon \bm{F}) = \trace (\bm{F} \bm{E} \bm{F}^T \diff \epsilon) = \trace(\bm{b} \bm{e} \diff \bm \epsilon ).
\end{aligned}
$$

:::

:::{dropdown} Representation of $\bm F \diff\bm S \bm F^T$ for isochoric Ogden material
Based on the split {eq}`dS-split-ogden`, the $\bm{F} \diff \bm{S}_{\text{vol}}^p \bm{F}^T$ and $\bm{F} \diff \bm{S}_{\text{vol}}^u \bm{F}^T$ are the same as {eq}`F-dS-p-FTranspose`, {eq}`F-dS-vol-FTranspose-mixed` (or {eq}`F-dS-vol-FTranspose` if we use single field formulation ) and the isochoric part is

$$
  \begin{aligned}
  \bm{F} \diff \bm{S}_{\text{iso}} \bm{F}^{T} &= \sum_{i=1}^3 \beta^2_i \diff S_i^{iso} \hat{\bm{n}_i} \hat{\bm{n}_i}^T + \beta_i S_i^{iso} \left( \bm{F} \diff \hat{\bm{N}_i} \hat{\bm{n}_i}^T + \hat{\bm{n}_i} \diff \hat{\bm{N}_i}^T \bm{F}^T \right) \\
  &= \sum_{i=1}^3 \left(\diff \tau_i^{iso} - 2 \frac{\diff \beta_i}{\beta_i} \tau_i^{iso} \right) \hat{\bm{n}_i} \hat{\bm{n}_i}^T + \tau_i^{iso} \bm{A}_i
  \end{aligned}
$$ (F-dS-iso-FTranspose-ogden)

where for example $\diff \tau_1^{iso}$ can be written

$$
  \begin{aligned}
  \diff \tau_1^{iso} &= \sum_{j=1}^N \frac{m_j \alpha_j}{3} \left[ 2 \diff \ell_1 \exp(\alpha_j \ell_1) - \diff \ell_2 \exp(\alpha_j \ell_2) - \diff \ell_3 \exp(\alpha_j \ell_3) \right] J^{-\alpha_j/3} \\
  &- \sum_{j=1}^N \frac{m_j \alpha_j}{9} \left[ 2 \operatorname{\tt expm1}(\alpha_j \ell_1) - \operatorname{\tt expm1}(\alpha_j \ell_2) - \operatorname{\tt expm1}(\alpha_j \ell_3) \right] J^{-\alpha_j/3} \trace \left(\diff \bm \epsilon \right)
  \end{aligned}
$$

and

$$
\bm{A}_i = 2 \frac{\diff \beta_i}{\beta_i} \hat{\bm{n}_i} \hat{\bm{n}_i}^T + \diff \left( \hat{\bm{n}_i} \hat{\bm{n}_i}^T \right) + \left( \left( \nabla_x \diff \bm{u} \right) \hat{\bm{n}_i} \hat{\bm{n}_i}^T + \hat{\bm{n}_i} \hat{\bm{n}_i}^T \left( \nabla_x \diff \bm{u} \right)^T \right)
$$


where $\tau_i^{iso}$ is defined in {eq}`mixed-ogden-tau-iso-stable` and we have used
$\bm{F} \hat{\bm{N}_i} = \beta_i \hat{\bm{n}_i}$ and $\tau_i^{iso} = \beta_i^2 S_i^{iso}$. Note that the linearization of principal stretches, $\beta_i$, and eigenvectors, $\hat{\bm{n}_i}$, are similar to initial configuration but are written in terms of Green-Euler tensor i.e.,

$$
\diff \beta_i = \frac{1}{\beta_i} \langle \hat{\bm{n}_i}, \diff \bm{e} \hat{\bm{n}_i} \rangle.
$$

$$
\quad \diff \hat{\bm{n}_i} = \sum_{j \neq i} \frac{1}{\beta_i^e - \beta_j^e} \langle \diff \bm{e} \hat{\bm{n}_i}, \hat{\bm{n}_j} \rangle \hat{\bm{n}_j}.
$$

:::

## Perturbed Lagrange-multiplier method

In the previous sections we introduced the strong form of general mixed formulation for small and finite strain and by considering test functions $(\bm v, q)$ as derived in {eq}`mixed-linear-weak-form` and {eq}`mixed-hyperelastic-weak-form-initial`. Alternatively, we can derive the mixed $\bm u- p$ weak formulation based on minimization of two fields functional $\Pi(\bm u, p)$, which it is known as the perturbed Lagrangian approach. For the mixed linear case, we can write

$$
  \begin{aligned}
    \Pi (\bm u, p) &= \int_{\Omega} \left[ \mu \, \bm \varepsilon_{\text{dev}} \tcolon \bm \varepsilon_{\text{dev}} - p\,  \trace \bm \varepsilon + \frac{\bulk_p}{2}  \left(\trace \bm \varepsilon \right)^2 - \frac{p^2}{2 (\bulk-\bulk_p)}  \right] \, dv - \Pi_{\text{ext}} (\bm u) \\
    \Pi_{\text{ext}} (\bm u) &= \int_{\Omega}{\bm{u} \cdot \rho \bm{g}} \, dv + \int_{\partial \Omega}{\bm{u} \cdot \bar{\bm{t}}} \, da
  \end{aligned}
$$ (mixed-linear-functional)

and by invoking the stationarity of $\Pi$ with respect to $\bm u$ and $p$, we obtain

$$
  \begin{aligned}
      \int_{\Omega} \nabla \delta \bm u \tcolon \left[ 2 \mu \, \bm \varepsilon_{\text{dev}} + \left(\bulk_p \trace \bm \varepsilon - p \right) \bm{I} \right] \, dv - \bm L_{\text{ext}}(\delta \bm u) &=0, \\
      \int_{\Omega} \delta p \left( -\trace \bm \varepsilon -  \frac{p}{\bulk-\bulk_p} \right) \, dv &= 0,
  \end{aligned}
$$ (mixed-linear-weak-form-PL)

where $\bm L_{\text{ext}}(\delta \bm u) = \int_{\Omega}{\delta \bm{u} \cdot \rho \bm{g}} \, dv + \int_{\partial \Omega}{\delta \bm{u} \cdot \bar{\bm{t}}} \, da$, and $\delta \bm u$ and $\delta p$ are virtual displacement and pressure and can be seen as the test functions $\bm v, q$, i.e., $\delta \bm u = \bm v, \, \delta p = q$. It is clear that the weak form {eq}`mixed-linear-weak-form-PL` agrees with what we obtained in {eq}`mixed-linear-weak-form`. However, the hyperelastic weak form {eq}`mixed-hyperelastic-weak-form-initial`, can not be derived by minimizing any functional and its linearization is not symmetric.

To write the mixed functional $\Pi(\bm u, p)$ for hyperelastic, we need to consider the strain energy of the form

$$
  \psi \left(\bm{C} \right) = \psi_{\text{vol}}(J) + \psi_{\text{iso}}(\bar{\bm{C}}) = \frac{\bulk}{2} \left(U(J)\right)^2 + \psi_{\text{iso}}(\bar{\bm{C}})
$$ (mixed-strain-energy-PL)

we can write a two fields energy functional as

$$
  \begin{aligned}
  \Pi (\bm u, p) &= \int_{\Omega_0} \left[ \psi_{\text{iso}}(\bar{\bm{C}}) - p \, U(J) - \frac{1}{2}\frac{p^2}{\bulk - \bulk_p} + \frac{\bulk_p}{2} U^2\right] \, dV - \Pi_{\text{ext}} (\bm u) \\
  \Pi_{\text{ext}} (\bm u) &= \int_{\Omega_0}{\bm{u} \cdot \rho_0 \bm{g}} \, dV + \int_{\partial \Omega_0}{\bm{u} \cdot \bar{\bm{t}}} \, dA
  \end{aligned}
$$ (mixed-energy-functional)

Finding the stationary conditions with respect to $\bm{u}$ and $p$ by taking Gateaux derivative gives the weak form

$$
  \begin{aligned}
    \int_{\Omega_0} \bm F \underbrace{\left( \bm{S}_{\text{iso}} + (\bulk_p U - p) \, J \, U' \, \bm{C}^{-1} \right)}_{\bm S} \tcolon \nabla_X \bm v \, dV &= L_{\text{ext}} (\bm v) \\
    \int_{\Omega_0} \left(- U(J) - \frac{p}{\bulk - \bulk_p} \right) q \, dV &= 0
  \end{aligned}
$$ (mixed-hyperelastic-weak-form-initial-PL)

where, $L_{\text{ext}} (\bm v) = \int_{\Omega_0}{\bm{v} \cdot \rho_0 \bm{g}} \, dV + \int_{\partial \Omega_0}{\bm{v} \cdot \bar{\bm{t}}} \, dA $ and we have used

$$
  \begin{aligned}
    \frac{\partial \psi_{\text{iso}}}{\partial \bm u} \cdot {\delta \bm u} &= \frac{\partial \psi_{\text{iso}}}{\partial \bm E} \tcolon \frac{\partial \bm E}{ \partial \bm u}{\delta \bm u} = \bm{S}_{\text{iso}} \tcolon {\delta \bm E} = \bm{S}_{\text{iso}} \tcolon \text{sym} \left(\bm F^T \delta \bm F \right) \\
    \frac{\partial U}{\partial \bm u} \cdot {\delta \bm u} &= \frac{\partial U}{\partial J} \frac{\partial J}{\partial \bm u} {\delta \bm u} = U' {\delta J} = J \, U' \, \bm{C}^{-1} \tcolon \tcolon {\delta \bm E} = J \, U' \, \bm{C}^{-1} \tcolon \text{sym} \left(\bm F^T \delta \bm F \right)
  \end{aligned}
$$

where $\delta \bm{F} = \nabla_X \delta \bm u = \nabla_X \bm v$. The Jacobian for problem {eq}`mixed-hyperelastic-weak-form-initial-PL` can be written as

$$
  \begin{aligned}
    \int_{\Omega_0} \nabla_X \bm{v} \tcolon \left(\bm F \diff \bm{S} + \diff \bm{F} \bm{S} \right) dV &= -R_u^{PL}, \\
    \int_{\Omega_0} q \left( -\diff U - \frac{\diff p}{\bulk - \bulk_p} \right) dV &= -R_p^{PL},
  \end{aligned}
$$ (mixed-jacobian-weak-form-initial-PL)

where

$$
  \begin{aligned}
    \diff \bm{S} &= \diff \bm{S}_{\text{iso}} + \diff \bm{S}_{\text{vol}}^u  + \diff \bm{S}_{\text{vol}}^p, \\
    \diff \bm{S}_{\text{vol}}^u &= \left(\bulk_p U' \diff J \right) J U' \bm{C}^{-1} + \left(\bulk_p U - p \right) \left(\diff J \, U' \bm{C}^{-1} + J\, U^{''} \diff J \bm{C}^{-1} + J U' \diff \bm{C}^{-1} \right), \\
    &= \left[\bulk_p \left( J U' \right)^2 +  \left(\bulk_p U - p \right) \left( J U' + J^2  U^{''} \right) \right]  \left( \bm{C}^{-1} \tcolon \diff \bm E \right) \bm{C}^{-1} \\
    &+ \left(\bulk_p U - p \right) J U' \, \diff \bm{C}^{-1}\\
    \diff \bm{S}_{\text{vol}}^p &= - dp J U' \, \bm{C}^{-1} \\
    \diff U &= U' \diff J = J \, U' \, \bm{C}^{-1} \tcolon \diff \bm E,
  \end{aligned}
$$ (mixed-stationary-point-linearization)

where $\diff \bm{S}_{\text{iso}}$ is derived for Neo-Hookean, Mooney-Rivlin and Ogden in {eq}`dS-iso-neo-hookean`, {eq}`dS-iso-mooney-rivlin` and {eq}`dS-iso-ogden`. To complete the derivation we only need $U(J)$ function with condition

$$
  U(J) = 0 \quad \text{if and only if} \quad J = 1.
$$

For the volumetric strain energy function of the form given in {eq}`mixed-neo-hookean-energy`, if we choose $V(J) = \frac{1}{4} \left(J^2 - 1 - 2 \log J \right)$, from $\bulk V(J) = \frac{\bulk}{2} U^2(J)$, we will have

$$
  U(J) = \pm \sqrt{2 V} = \frac{\sign(J-1)}{\sqrt{2}} ( \underbrace{J^2 - 1 - 2 \log J }_{A})^{1/2}
$$

where the derivatives are

$$
  \begin{aligned}
    U' &= \frac{\partial U}{\partial J} = \sign(J-1) \frac{J^2 - 1}{J \sqrt{2}} A^{-1/2}, \\
    U{''} &= \frac{\partial^2 U}{\partial J^2} = \sign(J-1) \frac{1}{J^2 \sqrt{2}} \left((J^2 + 1) A^{-1/2} - (J^2 - 1)^2 A^{-3/2} \right).
  \end{aligned}
$$

## Matrix-free implementation

Similar to {ref}`matrix-free-single-field` for single field, we present mathematical formulation of matrix-free method for a general Dirichlet problem for general mixed problem $\bm{R}(\bm{u}, p) = [\bm{R}_u(\bm{u}, p) ~~ \bm{R}_p(\bm{u}, p)]^T = \bm 0$: Find $(\bm u, p) \in \mathcal{V} \times \mathcal{Q} $, such that

$$
  \begin{aligned}
    \langle\bm{v},  \bm{R}_u(\bm{u}, p) \rangle &= \int_{\Omega_0}{\bm{v} \cdot \bm{f}_0 + \nabla \bm{v} \tcolon \bm{f}_1} \, dV = 0, \quad \forall \bm{v} \in \mathcal{V}, \\
		\langle q,  \bm{R}_p(\bm{u}, p) \rangle &= \int_{\Omega_0}{q \cdot \bm{g}_0 + \nabla q \tcolon \bm{g}_1 } \, dV = 0, \quad \forall q \in \mathcal{Q},
  \end{aligned}
$$ (mixed-residual-matrix-free)

where the operators $\bm{f}_0, \bm{f}_1, \bm{g}_0, \bm{g}_1$ contain all possible sources in the problem i.e., $\bm{u}, \nabla \bm{u}, p, \nabla p$. It should be noted that the gradient in the {eq}`mixed-residual-matrix-free` depends on the configuration system and it could be with respect to initial configuration $\bm{X}$ i.e.,  $\nabla_X \bm{v}$ (Lagrangian approach) or current configuration $\bm{x}$ i.e.,  $\nabla_x \bm{v}$ (Eulerian approach).

In order to solve {eq}`mixed-residual-matrix-free` with Newton-Krylov iterative solvers we need its Jacobian form: Find $(\diff \bm{u}, \diff p) \in \mathcal{V} \times \mathcal{Q} $ such that

$$
  \begin{aligned}
    \langle\bm{v},  \bm{J}_u(\bm{u}, p) \diff \bm{U} \rangle &= \int_{\Omega_0}{\bm{v} \cdot \diff \bm{f}_0 + \nabla \bm{v} \tcolon \diff \bm{f}_1} \, dV, \\
		\langle q,  \bm{J}_p(\bm{u}, p) \diff \bm{U} \rangle &= \int_{\Omega_0}{q \cdot \diff \bm{g}_0 + \nabla q \tcolon \diff \bm{g}_1} \, dV.
  \end{aligned}
$$ (mixed-jacobian-matrix-free)

where $\diff \bm{U} = [\diff \bm{u} ~~ \diff p]^T$ and the linearization of operators $\bm{f}_0, \bm{f}_1, \bm{g}_0, \bm{g}_1$ are

$$
  \begin{aligned}
    \diff \bm{f}_i &= \frac{\partial \bm{f}_i}{\partial \bm{u}} \diff \bm{u} + \frac{\partial \bm{f}_i}{\partial \nabla \bm{u}} \nabla \diff \bm{u} + \frac{\partial \bm{f}_i}{\partial p} \diff p + \frac{\partial \bm{f}_i}{\partial \nabla p} \nabla \diff p, \quad i=0,1 \\
		\diff \bm{g}_i &= \frac{\partial \bm{g}_i}{\partial \bm{u}} \diff \bm{u} + \frac{\partial \bm{g}_i}{\partial \nabla \bm{u}} \nabla \diff \bm{u} + \frac{\partial \bm{g}_i}{\partial p} \diff p + \frac{\partial \bm{g}_i}{\partial \nabla p} \nabla \diff p, \quad i=0,1
  \end{aligned}
$$

Compare with governing equations derived in pervious sections for linear and large deformation, it is easy to verify that

* For mixed linear elasticity described in {eq}`mixed-linear-weak-form` we have

$$
  \begin{aligned}
    \bm f_0 &= -\rho \bm g, \quad \bm f_1 = \bm{\sigma}(\bm u, p), \quad \bm g_0 =-\trace \bm \varepsilon -  \frac{p}{\bulk-\bulk_p}, \quad \bm g_1 = \bm 0, \\
			\diff \bm f_1 &= \bm f_1(\diff \bm{u}, \diff p), ~ ~ \diff \bm g_0 = \bm g_0(\diff \bm{u}, \diff p)
  \end{aligned}
$$

* For mixed hyperelastic in initial configuration described in {eq}`mixed-hyperelastic-weak-form-initial` and {eq}`mixed-jacobian-weak-form-initial` we have

$$
  \begin{aligned}
		\bm f_0 &= -\rho_0 \bm g, \quad \bm f_1 = \bm F \bm S, \quad \bm g_0 = L J, \quad \bm g_1 = \bm 0, \\
		\diff \bm f_1 &= \diff \bm F \bm S + \bm F \diff \bm S, ~ ~ \diff \bm g_0 = \diff L J + L \diff J
  \end{aligned}
$$

* For mixed hyperelastic in current configuration derived in {eq}`mixed-jacobian-weak-form-current`

$$
  \begin{aligned}
		\bm f_0 &= -\rho_0 \bm g, \quad \bm f_1 = \bm \tau, \quad \bm g_0 = L J, \quad \bm g_1 = \bm 0, \\
		\diff \bm f_1 &= \diff \bm{\tau} - \bm{\tau} \left( \nabla_x \diff \bm{u} \right)^T, ~ ~ \diff \bm g_0 = \diff L J + L \diff J
  \end{aligned}
$$

* For mixed hyperelastic with variational approach in initial configuration given in {eq}`mixed-hyperelastic-weak-form-initial-PL` and {eq}`mixed-jacobian-weak-form-initial-PL`

$$
  \begin{aligned}
		\bm f_0 &= -\rho_0 \bm g, \quad \bm f_1 = \bm F \bm S, \quad \bm g_0 = -{U}(J) - \frac{p}{\bulk - \bulk_p}, \quad \bm g_1 = \bm 0, \\
		\diff \bm f_1 &= \diff \bm F \bm S + \bm F \diff \bm S, ~ ~ \diff \bm g_0 = -\diff {U}(J) - \frac{\diff p}{\bulk - \bulk_p}
  \end{aligned}
$$

We solve the above mixed equations using static and quasistatic solver in Ratel as explained in elasticity [static and quasistatic](static-quasistatic) section. For mixed elastodynamic, we only need to add the same terms (see {eq}`elastodynamics-residual` and {eq}`elastodynamics-jacobian`) to the residual and jacobian of the first equation of mixed matrix-free formulation described in {eq}`mixed-residual-matrix-free` and {eq}`mixed-jacobian-matrix-free`, respectively.