# Shock Smoothing

## Artificial bulk viscosity

The history of the artificial shock viscosity can be traced back to the seminal work by {cite}`Neumann-Richtmyer1950` who proposed that a viscous term $q$ be added to an otherwise inviscid fluid's momentum equation for shocks propagating in one dimension, i.e.,

$$
\begin{aligned}
	q := c_0^2 \rho (\Delta x)^2\Big(\frac{\partial\dot{x}}{\partial x}\Big)^2\, ,
\end{aligned}
$$ (Shock-NeumannRichtmeyerQ)

where $q = 0$ for expanding motions ($\partial\dot{x}/\partial x \geq 0$), $x$ is the coordinate in the direction of motion, $\Delta x$ is the grid spacing, $\rho$ is the material density and $c_0$ is a constant $\approx 2$. {cite}`Landshoff1955` introduced a $q$ term that was linear in the velocity gradient, i.e.,

$$
\begin{aligned}
	q := c_L\rho\Delta x c\Big|\frac{\partial\dot{x}}{\partial x}\Big|\, ,
\end{aligned}
$$ (Shocks-LandshoffQ)

where $c_L$ is a constant $\approx 1$ and $c$ is the local sound speed, i.e., the square root of the local P-wave modulus divided by the local density. For further reading on the subject, we refer the interested reader to reviews by {cite}`Benson2007` and {cite}`Margolin2022`, with details on thermodynamical validity discussed in {cite}`Mattson2014`. In this work, we use the artificial viscosity given by {cite}`Wilkins1980` (which is also used in LS-DYNA {cite}`DYNA-Theory`) which combines linear and quadratic terms into one form, i.e.,

$$
\begin{aligned}
	q = \begin{cases} \rho l\Big(C_0 l (\nabla_x\cdot\dot{\bm{u}})^2 - C_1 c (\nabla_x\cdot\dot{\bm{u}})\Big) & \text{if }(\nabla_x\cdot\dot{\bm{u}}) < 0\\[0.5cm] 0 & \text{if } (\nabla_x\cdot\dot{\bm{u}})\geq 0\, ,\end{cases}
\end{aligned}
$$ (Shocks-WilkinsQ)

where $l$ is a characteristic length scale (in 3-D, this reduces to the cubed root of the element volume where the shock is occurring), $C_0$ and $C_1$ are constants typically taken to be 1.5 and 0.06, respectively.

The extent to which the shock viscosity effects the solution depends on the magnitude of the constants $C_0$ and $C_1$. The latter is responsible for damping out oscillations behind the front (i.e., the weak [elastic] waves produced by compression of the material, irrespective of the strength of the shock front) and the former is responsible for damping out oscillations near or ahead of the front (corresponding to stronger shocks). Large values of either constant tend to overly smooth out the shock over multiple elements; examples of shock-like loadings with varying ranges of the constants that demonstrate this phenomenon are shown in {cite}`Irwin2023c`.

The shock viscosity is essentially an opposing pressure term added on to the solid stress, such that the augmented solid Cauchy stress is written as

$$
\begin{aligned}
	\tilde{\bm{\sigma}} = \bm{\sigma}^{\rm{dev}} - (P + q)\bm{1}\, ,
\end{aligned}
$$

which, using the common definitions for deviatoric stress $\bm{\sigma}^{\rm{dev}}$ and hydrostatic pressure $P$, we may also write as,

$$
\begin{aligned}
	\tilde{\sigma}_{ij} = \begin{cases} \sigma_{ij} - \underbrace{\rho h[C_0 h(\nabla_x\cdot\dot{\bm{u}})^2 - C_1 c (\nabla_x\cdot\dot{\bm{u}})]}_{:= q}\delta_{ij} & \text{if }(\nabla_x\cdot\dot{\bm{u}}) < 0 \text{ (in compression)}\\[0.5cm] \sigma_{ij} & \text{if } (\nabla_x\cdot\dot{\bm{u}})\geq 0 \text{ (in tension)}\, .\end{cases}
\end{aligned}
$$

The weak formulation for a single-phase hyperelastic material is defined as

$$
    \begin{aligned}
        r(\bm u) \coloneqq \int_{\Omega_0}{\nabla_x \bm{v} \tcolon \bm{\tau}} \, dV
      - \int_{\Omega_0}{\bm{v} \cdot \rho_0 \bm{g}} \, dV
      - \int_{\Omega_0}{\bm{v} \cdot \rho_0\ddot{\bm{u}}} \, dV\\
      - \int_{\partial \Omega_0}{\bm{v} \cdot \bar{\bm t}} \, dA
         = 0, \quad \forall \bm{v} \in \mathcal{V}.
    \end{aligned}
$$

where we have included the dynamic term since presumably for static and quasistatic problems, the wave speeds induced by loading are too low to warrant inclusion of a shock viscosity to stabilize the numerical solution. Note that when shock viscosity is enabled (during compression), the following term is added to the residual above:

$$
    \int_{\Omega_0}{\nabla_x \bm{v} \tcolon (-Jq)\bm{I}} \, dV\, .
$$

It is convenient to write the quantity $Jq$ as

$$
    q_\tau := Jq = J^{1/3}\rho_0 h_0(J^{1/3}C_0h_0[\nabla_x\cdot\dot{\bm{u}}]^2 - C_1c[\nabla_x\cdot\dot{\bm{u}}])\, ,
$$

where we have used

$$
    h = (dv)^{1/3} = (JdV)^{1/3} = J^{1/3}h_0\, .
$$

The following term must also be added to the jacobian when the material is under compression:

$$
    \int_{\Omega_0}\nabla_x \bm{v} \tcolon [-\diff q_\tau\bm{I} + q_\tau\bm{I}(\nabla_x\diff\bm{u})^T] \, dV\, ,
$$

:::{dropdown} Derivation of $\diff q_\tau$

Define

$$
    q_\tau = J^{1/3} \rho_0 h_0 (A - B) \quad \text{where} \quad A = J^{1/3} C_0 h_0 (\nabla_x\cdot\dot{\bm{u}})^2\, , \, B = C_1 c (\nabla_x\cdot\dot{\bm{u}})\, .
$$

Then,

$$
    \diff q_\tau = J^{1/3} \rho_0 h_0 (\diff A - \diff B) + \frac{q_\tau}{J^{1/3}} \diff(J^{1/3})\, ,
$$

where

$$
    \begin{aligned}
        &\diff A &&= \diff(J^{1/3}) \, C_0 h_0 (\nabla_x\cdot\dot{\bm{u}})^2 + J^{1/3} C_0 h_0 \diff([\nabla_x\cdot\dot{\bm{u}}]^2)\\
        &\phantom{\diff A} &&= J^{1/3} C_0 h_0 (\nabla_x\cdot\dot{\bm{u}}) \left(\frac{(\nabla_x\cdot\dot{\bm{u}}) \trace\diff\bm{\epsilon}}{3} + 2 \mathrm{shift}_v (\nabla_x\cdot\diff\bm{u} - \nabla_x\bm{u} \tcolon (\nabla_x\diff\bm{u})^T)\right)\, ,\\
        &\diff B &&= C_1(\diff c \, (\nabla_x\cdot\dot{\bm{u}}) + c \diff(\nabla_x\cdot\dot{\bm{u}}))\\
        &\phantom{\diff B} &&= C_1 c \left(\frac{\trace\diff\bm{\epsilon}}{2} (\nabla_x\cdot\dot{\bm{u}}) + (\nabla_x\cdot\diff\bm{u} - \nabla_x\bm{u} \tcolon (\nabla_x\diff\bm{u})^T)\right)\, ,
    \end{aligned}
$$

and where we have used

$$
    \begin{aligned}
        &\diff(\nabla_x\cdot\dot{\bm{u}}) &&= \mathrm{shift_v}\diff(\nabla_x\cdot\bm{u}) = \mathrm{shift_v}(\nabla_x\cdot\diff\bm{u} - \nabla_x\bm{u}\tcolon(\nabla_x\diff\bm{u})^T)\, ,\\
        &\diff([\nabla_x\cdot{\bm{u}}]^2) &&= 2(\nabla_x\cdot\dot{\bm{u}})\diff(\nabla_x\cdot\dot{\bm{u}})\, ,\\
        &\diff J^{1/3} &&= \frac{1}{3J^{2/3}}\diff J = \frac{J}{3J^{2/3}}\trace\diff\bm{\epsilon} = \frac{J^{1/3}}{3}\trace\diff\bm{\epsilon}\, ,\\
        &\diff c &&= \diff\left(\frac{M}{\rho}\right) = -\frac{c}{2\rho}\diff\rho = -\frac{c}{2\rho}\diff\left(\frac{\rho_0}{J}\right) = \frac{c}{2\rho}\frac{\rho_0}{J^2}\diff J= \frac{c}{2}\trace\diff\bm{\epsilon}\, .
    \end{aligned}
$$

Thus,

$$
    \begin{aligned}
        \diff q_\tau &= J^{1/3} \rho_0 h_0 (\diff A - \diff B) + \frac{q_\tau \trace\diff\bm{\epsilon}}{3}\, .
    \end{aligned}
$$
:::

### Extension to poromechanics
For poromechanics, we assume the shock viscosity acts only open the solid phase, such that the solid skeleton effective stress is

$$
\begin{aligned}
	\tilde{\sigma}_{ij(E)}^\rs = \begin{cases} \sigma_{ij(E)}^\rs - \underbrace{\rho^\rs h[C_0 h (\nabla_x\cdot\dot{\bm{u}})^2 - C_1 c_\rs (\nabla_x\cdot\dot{\bm{u}})]}_{:= q}\delta_{ij} & \text{if }\nabla_x\cdot\dot{\bm{u}} < 0 \text{ (in compression)}\\[0.5cm] \sigma_{ij(E)}^\rs & \text{if } \nabla_x\cdot\dot{\bm{u}}\geq 0 \text{ (in tension)}\, ,\end{cases}
\end{aligned}
$$

where $c_\rs$ denotes the local sound speed through the solid constituent,

$$
  c_\rs := \sqrt{\frac{M_d}{\rho^\rs}}\, ,
$$

with $M_d$ the P-wave modulus of the drained mixture (solid skeleton).

Presently the poromechanics implementation is a work in progress.