[WIP] Hypervisco

[cmhamel]

This pull request implements a simple isotropic hyperviscoelastic model under the assumption of incompressible viscous flow with an arbitrary number of prony terms. The model form looks like the following

$\psi\left(\mathbf{F}, \mathbf{F}^e_1, \mathbf{F}^e_2, ..., \mathbf{F}^e_n, \theta\right) = \psi^{eq}\left(\mathbf{F}, \theta\right) + \sum_{i=1}^{n}\psi^{neq}_i\left(\mathbf{F}^e_i, \theta\right),$
where $\psi^{eq}$ is taken to be the strain energy of the hyperfoam model, and the $\psi^{neq}_i$'s are $n$ independent strain energies for different branches of a standard viscoelastic Prony series

Currently the equilibrium branch is a NeoHookean model, but we can generalize this later. The non-equilibrium branch currently has the following model form
$\psi^{neq}_i = G_i|\dev\mathbf{E}_i^e|^2$

Kinematics
$\mathbf{F} = \mathbf{F}^e_i\mathbf{F}^v_i,$

Flow rule update
$\dot{\mathbf{F}}^v_i = \mathbf{D}^v_i\mathbf{F}^v_i$

where the above is
$\mathbf{D}^v_i = \dot{\gamma}^v_i\left(\frac{\dev\mathbf{M}^e_i}{2\bar{M}_i}\right)$

For now a simple viscous shearing rate suitable for polymers in the rubbery state is used and is shown below

$\dot{\gamma}^v_i = \frac{\bar{M}_i}{G_i\tau_i}$

Constitutive model update [WIP]

Start with standard exponential rule
$\mathbf{F}^v_{n+1} = \exp\left(\Delta t\mathbf{D}^v_{n + 1}\right)\mathbf{F}^v_n,$

Do the usual stuff to show

$\mathbf{E}^e_{n + 1} = \mathbf{E}^e_{trial} -\Delta t\mathbf{D}^v_{n + 1}.$

$2G\mathbf{E}^e_{n + 1} = 2G\mathbf{E}^e_{trial} -2G\Delta t\mathbf{D}^v_{n + 1}.$

For this particular for of the viscous shearing rate we get lucky and have an analytic solution

$\mathbf{M}^e_{n + 1} = \mathbf{M}^e_{trial} - 2G\Delta t\left(\frac{1}{2G\tau}\mathbf{M}^e_{n + 1}\right).$

$\left(1 + \frac{\Delta t}{\tau}\right)\mathbf{M}^e_{n + 1} = \mathbf{M}^e_{trial},$

$\mathbf{M}^e_{n + 1} = \frac{\mathbf{M}^e_{trial}}{\left(1 + \frac{\Delta t}{\tau}\right)}.$

What I'm unsure about is what the variational update looks like...

Currently, I'm simply return this
$\psi_{eq} + \psi_{neq} + \frac{1}{2}\Delta t \eta \dot{\gamma}_v^2$

where $\eta = G\tau$

Also in the pull request a patched up some interface issues with mechanics functions related to a time step.

[cmhamel]

Leaving it as a draft so we can iterate on what model form/features will suite all of our needs

[btalamini]

Late comment. Can I ask, why did dt get moved into the factory function?

[cmhamel]

If I remember correctly, things weren't properly connected so dt wasn't making its way to the material models. Maybe I misunderstood the interface as it was though?

[btalamini]

The J2 model had this function so that we could put the finite strain and the infinitesimal strain versions in one implementation. For this model, which only makes sense as a finite deformation model, we might as well do all of the update in the _compute_state_new() function which will make it easier to read. Maybe this will speed things up a little, too.

[btalamini]

FWIW, If we end up implementing multiple Prony terms, I'd suggest trying with a jax.lax.fori loop, or jax.lax.scan. I think it might be more straightforward than vmap.

[ralberd]

Sorry, I kind of rushed and merged this in. Should we create a separate pull request to address Brandon's comments?