Implementing viscoelasticity with state variables

[bhaveshshrimali]

Hi @adtzlr ,

I really like the matadi approach. It is performant, as well as offloads the tangent calculations to AD. I am starting to think about viscoelaticity and how it can be implemented within felupe. Consider the strain energy:

def ψ(x, mu1, mu2, alph1, alph2, a1, a2, m1, m2, K1, K2, bta1, bta2, eta0, etaInf, λ, dt):
    F, p, ζ = x[0], x[1], x[2]

    J = det(F)
    Js = (λ + p + sqrt((λ + p) ** 2 + 4 * λ * (mu1 + mu2 + m1 + m2))) / (2 * λ)

    C = F.T @ F
    Cvinv = inv(ζ)
    C_Cvinv = C @ Cvinv

    Iv1 = trace(ζ)

    I1 = trace(C)
    Ie1 = trace(C_Cvinv)

    W_I1 = 3 ** (1 - alph1) / 2 / alph1 * mu1 * (I1**alph1 - 3**alph1) + 3 ** (
        1 - alph2
    ) / 2 / alph2 * mu2 * (I1**alph2 - 3**alph2)
    W_Ie1 = 3 ** (1 - a1) / 2 / a1* m1 * (Ie1**a1- 3**a1) + 3 ** (
        1 - a2
    ) / 2 / a2* m2* (Ie1**a2 - 3**a2)
    ψs = p * (J - Js) + W_I1 - (mu1 + mu2 + m1 + m2) * log(Js) + λ / 2 * (Js - 1) ** 2

    ζ_new = ζ + rk_update() # this update of the state variable

    return gradient(ψs, F), gradient(ψs, p), ζ_new

where the update function for the internal state variables ζ (if I understand these correctly) can be an explicit (like Runge-Kutta) or implicit (backward euler). The latter would need a Newton solver embedded within the constitutive relation and hence I don't want to get into that yet (or maybe its easily implementable). For now, consider that it is fully explicit, i.e. knowing at time-step n the strain (Cn) the internal vairable (ζn), we can fully determine the update ζ given the increment dt.

def rk_update(C, Cv, Cn, Cvn):
    pass

Concretely, I am planning to implement the following constitutive behavior

with

and the ODE G(t, Cv) can be discretized using a time-discretization of your choice (in this case, say RK4 or RK5)

Any thoughts on how this could be achieved within felupe with/without using matadi.

[bhaveshshrimali]

And the constitutive model is here: https://pamies.cee.illinois.edu/Publications_files/CRM_2016.pdf

Specifically, equations (36)-(37) along with (40)-(41)

[adtzlr]

Hi @bhaveshshrimali,

you have to split / concatenate the state variables Cn and Cv (hopefully vertsplit is correct here, I'm not 100% sure right now):

Edit: I think you need only one state variable in this model, Cn isn't necessary. Thus, you may skip the vertsplit/vertcat stuff.

Code:

import matadi as mat
from matadi.math import (
    det,
    trace,
    inv,
    sqrt,
    gradient,
    log,
    astensor,
    asvoigt,
    vertsplit,
    vertcat,
)

F = mat.Variable("F", 3, 3)
p = mat.Variable("p", 1, 1)
ζ = mat.Variable("ζ", 12, 1)

def ψ(
    x, mu1, mu2, alph1, alph2, a1, a2, m1, m2, K1, K2, bta1, bta2, eta0, etaInf, λ, dt
):
    F, p, ζ = x[0], x[1], x[2]

    Cn, Cvn = vertsplit(ζ, [0, 6, 12])
    Cn = astensor(Cn)
    Cvn = astensor(Cvn)

    J = det(F)
    Js = (λ + p + sqrt((λ + p) ** 2 + 4 * λ * (mu1 + mu2 + m1 + m2))) / (2 * λ)

    C = F.T @ F

    Cvinv = inv(Cv)
    C_Cvinv = C @ Cvinv

    Iv1 = trace(Cv)

    I1 = trace(C)
    Ie1 = trace(C_Cvinv)

    W_I1 = 3 ** (1 - alph1) / 2 / alph1 * mu1 * (I1**alph1 - 3**alph1) + 3 ** (
        1 - alph2
    ) / 2 / alph2 * mu2 * (I1**alph2 - 3**alph2)
    W_Ie1 = 3 ** (1 - a1) / 2 / a1 * m1 * (Ie1**a1 - 3**a1) + 3 ** (
        1 - a2
    ) / 2 / a2 * m2 * (Ie1**a2 - 3**a2)
    ψs = p * (J - Js) + W_I1 - (mu1 + mu2 + m1 + m2) * log(Js) + λ / 2 * (Js - 1) ** 2

    Cv = Cvn + rk_update(C, Cn, Cvn)  # this update of the state variable

    ζ_new = vertcat(asvoigt(C), asvoigt(Cv))

    return gradient(ψs, F), gradient(ψs, p), ζ_new

(this is still an incomplete function, I just added the state variable handling)

CasADi (the backend of matadi) also comes with an ode function. I can't help you here as I've never used it but hopefully their docs give you some idea. Otherwise proceed it as you described it - use your own update function of state variables for the evolution equation. For my understanding, you have to do it at the beginning (like in my code block).

Edit: I see, due to the explicit Runge-Kutta update the update relies on the old state variables.

[adtzlr]

I tried to implement it with pure FElupe using felupe.Hyperelastic - but I failed somewhere.

If you're interested, here is the code. Syntax is a bit different compared to matadi.

import numpy as np
import felupe as fem
import tensortrax.math as tm


def viscoelastic_two_potential(C, ζn, dt, μ, α, m, a, η0, ηinf, β, K):
    # extract old state variables
    Cvn = tm.special.from_triu_1d(ζn, like=C)

    # third invariant of the right Cauchy-Green deformation tensor
    # first invariant of distortional part of the right Cauchy-Green deformation tensor
    I3 = tm.linalg.det(C)
    I1 = tm.trace(C) * I3 ** (-1 / 3)

    def invariants(Cv):
        "Invariants of the elastic and viscous parts of the deformation."
        Ce = C @ tm.linalg.inv(Cv)

        I1e = I3 ** (-1 / 3) * tm.trace(Ce)
        I1v = tm.trace(Cv)
        I2e = I3 ** (-2 / 3) * (I1e**2 - tm.trace(Ce @ Ce)) / 2

        return I1e, I1v, I2e

    def G(Cv):
        "Evolution equation."
        I1e, I1v, I2e = invariants(Cv)

        J2Neq = (I1e**2 / 3 - I2e) * tm.sum(
            [3 ** (1 - a[r]) * m[r] * I1e ** (a[r] - 1) for r in [0, 1]]
        ) ** 2
        u = [3 ** (1 - a[r]) * m[r] * I1e ** (a[r] - 1) for r in [0, 1]]
        ηK = ηinf + (η0 - ηinf + K[0] * (I1v ** β[0] - 3 ** β[0])) / (
            1 + (K[1] * J2Neq) ** β[1]
        )

        return tm.sum(u) / ηK * (C - I1e / 3 * tm.linalg.inv(Cv))

    # update the state variables
    k1 = G(Cvn)
    k2 = G(Cvn + k1 * dt / 2)
    k3 = G(Cvn + (3 * k1 + k2) * dt / 16)
    k4 = G(Cvn + k3 * dt / 2)
    k5 = G(Cvn + 3 * (-k2 + 2 * k3 + 3 * k4) * dt / 16)
    k6 = G(Cvn + (k1 + 4 * k2 + 6 * k3 - 12 * k4 + 8 * k5) * dt / 7)
    Cv = Cvn + dt / 90 * (7 * k1 + 32 * k3 + 12 * k4 + 32 * k5 + 7 * k6)

    I1e, I1v, I2e = invariants(Cv)

    # strain energy functions of the equilibrium and non-equilibrium parts
    p = [0, 1]
    ψEq = [3 ** (1 - α[r]) / (2 * α[r]) * μ[r] * (I1 ** α[r] - 3 ** α[r]) for r in p]
    ψNEq = [3 ** (1 - a[r]) / (2 * a[r]) * m[r] * (I1e ** a[r] - 3 ** a[r]) for r in p]

    return tm.sum(ψEq) + tm.sum(ψNEq), tm.special.triu_1d(Cv)


λ = fem.math.linsteps([1, 2, 1.5], num=[10, 5])
dλ = np.abs(np.diff(λ)).mean()
dλdt = 0.00023

umat = fem.Hyperelastic(
    viscoelastic_two_potential,
    nstatevars=6,
    dt=dλ / dλdt,
    μ=[1.08, 0.0170],
    α=[0.26, 7.68],
    m=[1.57, 0.59],
    a=[-10, 7.53],
    η0=2.11,
    ηinf=0.1,
    β=[3.0, 1.929],
    K=[442, 1289.49],
)

mesh = fem.Cube(n=2)
region = fem.RegionHexahedron(mesh)
field = fem.FieldContainer([fem.Field(region, dim=3)])
boundaries, loadcase = fem.dof.uniaxial(field, clamped=False)

solid = fem.SolidBodyNearlyIncompressible(umat, field, bulk=5000)
solid.results.statevars[[0, 3, 5]] += 1.0

step = fem.Step(items=[solid], ramp={boundaries["move"]: λ - 1}, boundaries=boundaries)
job = fem.CharacteristicCurve(steps=[step], boundary=boundaries["move"])
job.evaluate()
fig, ax = job.plot(
    xlabel="Displacement $u$ in mm $\longrightarrow$",
    ylabel="Normal Force $F$ in N $\longrightarrow$",
)
img = solid.imshow("Principal Values of Cauchy Stress")

[bhaveshshrimali]

Thanks a lot, Andreas!! I tried the simplest patch test, i.e. Cube(n=1) with the material parameters for VHB. It does run but the reaction force output is incorrect. I will try to get to debug this over the weekend once I have more time.

My approach to debugging would be to spell out the timesteps explicitly (which is currently masked behind job.evaluate()), and within each time-step the Jacobian matrix and residual vector explicitly, such that I can step it through the debugger to see if the state-variable update and/or the material tangent is incorrect. If you have some notes on how that could look like, that would be really helpful. Otherwise, I'll get to it during the weekend

Once this works out, then it makes FElupe in my opinion edge past all the currently available packages in terms of performance and ease of implementation.

import numpy as np
import felupe as fem
import tensortrax.math as tm
from matplotlib import pyplot as plt


def viscoelastic_two_potential(C, ζn, dt, μ, α, m, a, η0, ηinf, β, K):
    # extract old state variables
    Cvn = tm.special.from_triu_1d(ζn, like=C)

    # third invariant of the right Cauchy-Green deformation tensor
    # first invariant of distortional part of the right Cauchy-Green deformation tensor
    I3 = tm.sqrt(tm.linalg.det(C))
    I1 = tm.trace(C)  # * I3 ** (-1 / 3)

    def invariants(Cv):
        "Invariants of the total, elastic and viscous parts of the deformation."
        Ce = C @ tm.linalg.inv(Cv)
        CvinvC = tm.linalg.inv(Cv) @ C
        I1e = tm.trace(Ce)
        I1v = tm.trace(Cv)
        I2e = (I1e**2 - tm.trace(CvinvC @ Ce)) / 2

        return I1e, I1v, I2e

    def G(Cv):
        "Evolution equation."
        I1e, I1v, I2e = invariants(Cv)

        J2Neq = (I1e**2 / 3 - I2e) * tm.sum(
            [3 ** (1 - a[r]) * m[r] * I1e ** (a[r] - 1) for r in [0, 1]]
        ) ** 2
        u = [3 ** (1 - a[r]) * m[r] * I1e ** (a[r] - 1) for r in [0, 1]]
        ηK = ηinf + (η0 - ηinf + K[0] * (I1v ** β[0] - 3 ** β[0])) / (
            1 + (K[1] * J2Neq) ** β[1]
        )

        return tm.sum(u) / ηK * (C - I1e / 3 * Cv)

    # update the state variables
    k1 = G(Cvn)
    k2 = G(Cvn + k1 * dt / 2)
    k3 = G(Cvn + (3 * k1 + k2) * dt / 16)
    k4 = G(Cvn + k3 * dt / 2)
    k5 = G(Cvn + 3 * (-k2 + 2 * k3 + 3 * k4) * dt / 16)
    k6 = G(Cvn + (k1 + 4 * k2 + 6 * k3 - 12 * k4 + 8 * k5) * dt / 7)
    Cv = Cvn + dt / 90 * (7 * k1 + 32 * k3 + 12 * k4 + 32 * k5 + 7 * k6)

    I1e, I1v, I2e = invariants(Cvn)

    # strain energy functions of the equilibrium and non-equilibrium parts
    p = [0, 1]
    ψEq = [3 ** (1 - α[r]) / (2 * α[r]) * μ[r] * (I1 ** α[r] - 3 ** α[r]) for r in p]
    ψNEq = [3 ** (1 - a[r]) / (2 * a[r]) * m[r] * (I1e ** a[r] - 3 ** a[r]) for r in p]

    return tm.sum(ψEq) + tm.sum(ψNEq), tm.special.triu_1d(Cv)


λ = fem.math.linsteps([1, 2, 1.5], num=[10, 5])
dλ = np.abs(np.diff(λ)).mean()
dλdt = 0.03

umat = fem.Hyperelastic(
    viscoelastic_two_potential,
    nstatevars=6,
    dt=dλ / dλdt,
    μ=[13.54, 1.08],
    α=[1.0, -2.474],
    m=[5.42, 20.78],
    a=[-10, 1.948],
    η0=7014.0,
    ηinf=0.1,
    β=[1.852, 0.26],
    K=[3507.0, 1.0],
)

mesh = fem.Cube(n=1)
region = fem.RegionHexahedron(mesh)
field = fem.FieldContainer([fem.Field(region, dim=3)])
boundaries, loadcase = fem.dof.uniaxial(field, clamped=False)

solid = fem.SolidBodyNearlyIncompressible(umat, field, bulk=5000)
solid.results.statevars[[0, 3, 5]] += 1.0

step = fem.Step(items=[solid], ramp={boundaries["move"]: λ - 1}, boundaries=boundaries)
job = fem.CharacteristicCurve(steps=[step], boundary=boundaries["move"])
job.evaluate()
fig, ax = job.plot(
    xlabel="Displacement $u$ in mm $\longrightarrow$",
    ylabel="Normal Force $F$ in N $\longrightarrow$",
)
plt.show()

[bhaveshshrimali]

Some initial thoughts. I would still need to incorporate Cn in the update function, since the rhs G(t, Cv) when evaluated at tn, tn+Δt/2, tn + Δt/4, tn+3/4*Δt would depend on Cn, Cvn and C. I'll think on how to do that.

Thanks again for this. Super-helpful!

[adtzlr]

Matadi could also be a good choice here. I just experimented with FElupe to have an example or a new hyperelastic model formulation to add. Pure FElupe uses Tensortrax - it's self-written AD. For simple strain energy functions it works really well. However, I have the feeling that it is super slow for things like the Runge-Kutta update of state variables. I have a few more things in mind to show you.

This discussion is great because it gives me some ideas what has to be documented to get this to work!

[adtzlr]

With matadi, there is still something wrong with the implementation. Please note that for a SolidBodyNearlyIncompressible it is assumed that a strain energy function on the distortional part of the deformation $\hat{\psi}(\boldsymbol{F})$ has to be provided. If you don't enforce this, stresses are not zero at the undeformed state. You'd have to add something like ...- (mu1 + mu2 + m1 + m2) * log(det(F)... as you did. But then the isochoric-volumetric split is not done the right way.

$$\psi(\boldsymbol{F}, p, \bar{J}) = \hat{\psi}(\boldsymbol{F}) + U(\bar{J}) + p (J - \bar{J})$$

where $U(\bar{J})=K (\bar{J}-1)^2 / 2$ is added by SolidBodyNearlyIncompressible. This means that the added bulk-modulus is not the final bulk modulus as $\psi(\boldsymbol{F})$ contains some hydrostatic part. Beside that, this could also lead to volumetric locking (haven't tried that yet). The solution for compressible materials is to use fem.ThreeFieldVariation and fem.SolidBody.

$$\psi(\bar{\boldsymbol{F}}, p, \bar{J}) = {\psi}(\bar{\boldsymbol{F}}) + p (J - \bar{J})$$

with

$$ \bar{\boldsymbol{F}} = \left ( \frac{\bar{J}}{\det(\boldsymbol{F})} \right)^{-1/3} \boldsymbol{F} $$

Or, sometimes the default displacement-based formulation is also fine.

Code

import numpy as np
import felupe as fem
from pypardiso import spsolve

import matadi as mat
from matadi.math import (
    det,
    trace,
    inv,
    gradient,
)

F = mat.Variable("F", 3, 3)
ζ = mat.Variable("ζ", 3, 3)


def ψ(x, dt, μ, α, m, a, η0, ηinf, β, K):
    F, ζn = x[0], x[1]
    
    # extract old state variables
    Cvn = ζn
    
    J = det(F)
    C = (F.T @ F) * J ** (-2 / 3)
    I1 = trace(C)

    def invariants(Cv):
        "Invariants of the elastic and viscous parts of the deformation."
        Ce = C @ inv(Cv)
        I1e = trace(Ce)
        I1v = trace(Cv)
        I2e = (I1e**2 - trace(Ce @ Ce)) / 2
        return I1e, I1v, I2e

    I1e, I1v, I2e = invariants(Cvn)
    
    def G(Cv):
        "Evolution equation."
        I1e, I1v, I2e = invariants(Cv)
        J2Neq = (I1e**2 / 3 - I2e) * sum(
            [3 ** (1 - a[r]) * m[r] * I1e ** (a[r] - 1) for r in [0, 1]]
        ) ** 2
        u = [3 ** (1 - a[r]) * m[r] * I1e ** (a[r] - 1) for r in [0, 1]]
        ηK = ηinf + (η0 - ηinf + K[0] * (I1v ** β[0] - 3 ** β[0])) / (
            1 + (K[1] * J2Neq) ** β[1]
        )

        return sum(u) / ηK * (C - I1e / 3 * Cv)
    
    # a) Euler-explicit
    # I1e, I1v, I2e = invariants(Cvn)
    # J2Neq = (I1e**2 / 3 - I2e) * sum(
    #     [3 ** (1 - a[r]) * m[r] * I1e ** (a[r] - 1) for r in [0, 1]]
    # ) ** 2
    # u = [3 ** (1 - a[r]) * m[r] * I1e ** (a[r] - 1) for r in [0, 1]]
    # ηK = ηinf + (η0 - ηinf + K[0] * (I1v ** β[0] - 3 ** β[0])) / (
    #     1 + (K[1] * J2Neq) ** β[1]
    # )
    # Cv = Cvn + dt * (sum(u) / ηK * (C - I1e / 3 * Cvn))

    # b) Runge-Kutta
    # k1 = G(Cvn)
    # k2 = G(Cvn + k1 * dt / 2)
    # k3 = G(Cvn + (3 * k1 + k2) * dt / 16)
    # k4 = G(Cvn + k3 * dt / 2)
    # k5 = G(Cvn + 3 * (-k2 + 2 * k3 + 3 * k4) * dt / 16)
    # k6 = G(Cvn + (k1 + 4 * k2 + 6 * k3 - 12 * k4 + 8 * k5) * dt / 7)
    # Cv = Cvn + dt / 90 * (7 * k1 + 32 * k3 + 12 * k4 + 32 * k5 + 7 * k6)
    
    # c) disabled
    Cv = Cvn
    
    I1e, I1v, I2e = invariants(Cv)

    # strain energy functions of the equilibrium and non-equilibrium parts
    p = [0, 1]
    ψEq = [3 ** (1 - α[r]) / (2 * α[r]) * μ[r] * (I1 ** α[r] - 3 ** α[r]) for r in p]
    ψNEq = [3 ** (1 - a[r]) / (2 * a[r]) * m[r] * (I1e ** a[r] - 3 ** a[r]) for r in p]

    return gradient(sum(ψEq) + sum(ψNEq), F), Cv


umat = mat.MaterialTensor(
    x=[F, ζ],
    fun=ψ,
    statevars=1,
    kwargs=dict(
        dt=500,
        μ=[1.08, 0.0170],
        α=[0.26, 7.68],
        m=[1.57, 0.59],
        a=[-10, 7.53],
        η0=2.11,
        ηinf=0.1,
        β=[3.0, 1.929],
        K=[442, 1289.49],
    ),
    compress=True,
)
view = fem.ViewMaterialIncompressible(
    umat,
    ux=fem.math.linsteps([1, 1.5], num=[100]),
    bx=None,
    ps=None,
    statevars=np.eye(3).reshape(3, 3, 1, 1)
)
view.plot(
    show_kwargs=False,
)

mesh = fem.Cube(n=11)
region = fem.RegionHexahedron(mesh)
field = fem.FieldContainer([fem.Field(region, dim=3)])
boundaries, loadcase = fem.dof.uniaxial(field, clamped=False)

solid = fem.SolidBodyNearlyIncompressible(umat, field, bulk=5000)
solid.results.statevars.T[:] = np.eye(3)

move = fem.math.linsteps([0, 1], num=10)
step = fem.Step(items=[solid], ramp={boundaries["move"]: move}, boundaries=boundaries)
job = fem.CharacteristicCurve(steps=[step], boundary=boundaries["move"])
job.evaluate(solver=spsolve, tol=1e-2)
fig, ax = job.plot(
    xlabel="Displacement $u$ in mm $\longrightarrow$",
    ylabel="Normal Force $F$ in N $\longrightarrow$",
)
solid.imshow("Principal Values of Cauchy Stress")

The equilibrium part works as expected.

The analytic uniaxial incompressible deformation

via fem.ViewMaterialIncompressible

and the solution of FElupe

What is missing

• The viscoelastic evolution equation (update of state variables), see Implementing viscoelasticity with state variables #666 (comment)

[adtzlr]

A very basic evolution equation (inspired by [1])

m = 1.0
Cv = Cvn + dt * m * C
Cv *= det(Cv) ** (-1 / 3)

completes with no errors. So the error is definitelty in the G(t, Cv)-function.

[adtzlr]

Some initial thoughts. I would still need to incorporate Cn in the update function, since the rhs G(t, Cv) when evaluated at tn, tn+Δt/2, tn + Δt/4, tn+3/4*Δt would depend on Cn, Cvn and C. I'll think on how to do that.

Hmm. I don't find Δt on the rhs of G(t, Cv) if I use the Runge-Kutta method. Δt is only used in the update of Cv(t=tn). The only quantities needed are the old viscoelastic Cv(t=tn) and the time increment Δt. Could you point me to C(t=tn)? I must have overlooked that.

[bhaveshshrimali]

Correct, but the rhs is also a funciton of C, and hence the intermediate evaluations should also, in principle, depend on it. E.g. Backward Euler would look like:

Cv_n+1 = Cv_n + Δt⋅G(C, Cv_n+1)

And an RK5 update could look like

let,
C₀₅ = C(tn + dt/2),
C₀₂₅ = C(tn + dt/4),
C₀₇₅ = C(tn + 3*dt/4),
C = (tn + dt) # not known -- updated during Newton iterations

where I would linearly interpolate for the intermediate time-steps (between the current newton iterate and Cn)

and correspondingly,

k1 = G(Cn, Cvn)
k2 = G(C₀₅, Cvn + k1 * dt / 2)
k3 = G(C₀₂₅, Cvn + (3 * k1 + k2) * dt / 16)
k4 = G(C₀₅, Cvn + k3 * dt / 2)
k5 = G(C₀₇₅, Cvn + 3 * (-k2 + 2 * k3 + 3 * k4) * dt / 16)
k6 = G(C, Cvn + (k1 + 4 * k2 + 6 * k3 - 12 * k4 + 8 * k5) * dt / 7)
Cv = Cvn + dt / 90 * (7 * k1 + 32 * k3 + 12 * k4 + 32 * k5 + 7 * k6)

Hence, although it is fully explicit in Cvn, it does depend on C, as described here: here

[adtzlr]

If you need more than one state variable, split and concatenate the vector of state variables. From/to Voigt-Notation is optional, but reduces the number of needed state variables for symmetric tensors.

from matadi.math import gradient, astensor, asvoigt, vertcat, vertsplit

F = mat.Variable("F", 3, 3)
ζ = mat.Variable("ζ", 12)

def ψ(x, **kwargs):
    F, ζn = x[0], x[1]

    # extract old state variables
    ζ1n, ζ2n = vertsplit(ζn, [0, 6, 12])
    Cvn = astensor(ζ1n)
    Cn = astensor(ζ2n)

    # ...

    return gradient(W, F), vertcat(asvoigt(Cv), asvoigt(C))

Without From/To Voigt-Notation:

ζ = mat.Variable("ζ", 6, 3)

def ψ(x, **kwargs):
    F, ζn = x[0], x[1]
    Cvn, Cn = vertsplit(ζn, [0, 3, 6])

    return gradient(W, F), vertcat(Cv, C)

[bhaveshshrimali]

Correct, I got that part from your first comment. I guess my question would then be, is Cn=C(t=tn), because that is what I need. And this assignment would have to be done after every converged solution is obtained.

I would have to step in through a debugger to understand the intrinsics of how the stiffness calculation works inside FElupe.

[adtzlr]

I guess my question would then be, is Cn=C(t=tn), because that is what I need. And this assignment would have to be done after every converged solution is obtained.

When using a SolidBody (or the nearly-incompressible variant), normally you don't have to take care of the state-variables update mechanism.

I still haven't figured out yet how to treat local Newton iterations...

Below follows a brief overview on technical details.

State Variables in FElupe

The state variables are stored in a private (temporary) attribute SolidBody.results._statevars for each stress evaluation, where umat.gradient([F, ζ], **kwargs) is called inside SolidBody.evaluate.gradient(field). The code is here:

felupe/src/felupe/mechanics/_solidbody.py

Lines 293 to 300 in a79926f

	gradient = self.umat.gradient(
	[*self.results.kinematics, self.results.statevars], *args, **kwargs
	)
	self.results.gradient = gradient[0]
	self.results.stress, self.results._statevars = gradient[:-1], gradient[-1]
	return self.results.stress

After each converted Newton iteration, the state variables are written to SolidBody.results.statevars.

felupe/src/felupe/mechanics/_helpers.py

Lines 72 to 74 in a79926f

	def update_statevars(self):
	if self._statevars is not None:
	self.statevars = self._statevars

This is done inside the default check()-function of newtonrhapson(), which determines whether the solution has converged or not.

felupe/src/felupe/tools/_newton.py

Lines 188 to 190 in a79926f

	if success and items is not None:
	for item in items:
	[item.results.update_statevars() for item in items]

I know, this is a bit tricky to understand - but it works great! 🚀

Automatic Differentiation using matadi (with CasADi as backend)

With matadi, things are a bit more difficult due to AD. You don't operate on NumPy arrays directly inside the strain-energy function but on CasADi symbolic variables instead. The numeric values are mapped and evaluated on the symbolic variables.

[adtzlr]

The $\beta$-parameter seems to make the model unstable. Or I made a typo somewhere. With $\beta_r=1$ it works. If I add a small number like 1e-6 in $\eta_K$, the code is much more stable, even with the original parameters.

def G(Cv):
    # ....
    ηK = ηinf + (η0 - ηinf + K[0] * (I1v ** β[0] - 3 ** β[0])) / (
                1 + (K[1] * J2Neq + 1e-6) ** β[1]
            )
    # ...

Due to the Runge-Kutta method, small increments are necessary. Newton-iterations would definitely improve this. Shouldn't be too hard to implement.

Edit: See #741 for details on implementing Newton iterations.

Code (with matadi)

import numpy as np
import felupe as fem
from pypardiso import spsolve

import matadi as mat
from matadi.math import (
    det,
    trace,
    inv,
    gradient,
)

F = mat.Variable("F", 3, 3)
ζ = mat.Variable("ζ", 3, 3)


def ψ(x, dt, μ, α, m, a, η0, ηinf, β, K):
    F, ζn = x[0], x[1]
    
    # extract old state variables
    Cvn = ζn
    
    J = det(F)
    C = (F.T @ F) * J ** (-2 / 3)
    I1 = trace(C)

    def invariants(Cv):
        "Invariants of the elastic and viscous parts of the deformation."
        Ce = C @ inv(Cv)
        I1e = trace(Ce)
        I1v = trace(Cv)
        I2e = (I1e**2 - trace(Ce @ Ce)) / 2
        return I1e, I1v, I2e

    I1e, I1v, I2e = invariants(Cvn)
    
    def G(Cv):
        "Evolution equation."
        I1e, I1v, I2e = invariants(Cv)
        J2Neq = (I1e**2 / 3 - I2e) * sum(
            [3 ** (1 - a[r]) * m[r] * I1e ** (a[r] - 1) for r in [0, 1]]
        ) ** 2
        u = [3 ** (1 - a[r]) * m[r] * I1e ** (a[r] - 1) for r in [0, 1]]
        ηK = ηinf + (η0 - ηinf + K[0] * (I1v ** β[0] - 3 ** β[0])) / (
            1 + (K[1] * J2Neq + 1e-6) ** β[1]
        )

        return sum(u) / ηK * (C - I1e / 3 * Cv)

    k1 = G(Cvn)
    k2 = G(Cvn + k1 * dt / 2)
    k3 = G(Cvn + (3 * k1 + k2) * dt / 16)
    k4 = G(Cvn + k3 * dt / 2)
    k5 = G(Cvn + 3 * (-k2 + 2 * k3 + 3 * k4) * dt / 16)
    k6 = G(Cvn + (k1 + 4 * k2 + 6 * k3 - 12 * k4 + 8 * k5) * dt / 7)
    Cv = Cvn + dt / 90 * (7 * k1 + 32 * k3 + 12 * k4 + 32 * k5 + 7 * k6)

    I1e, I1v, I2e = invariants(Cv)

    # strain energy functions of the equilibrium and non-equilibrium parts
    p = [0, 1]
    ψEq = [3 ** (1 - α[r]) / (2 * α[r]) * μ[r] * (I1 ** α[r] - 3 ** α[r]) for r in p]
    ψNEq = [3 ** (1 - a[r]) / (2 * a[r]) * m[r] * (I1e ** a[r] - 3 ** a[r]) for r in p]

    return gradient(sum(ψEq) + sum(ψNEq), F), Cv


umat = mat.MaterialTensor(
    x=[F, ζ],
    fun=ψ,
    statevars=1,
    kwargs=dict(
        dt=0.1,
        μ=[1.08, 0.0170],
        α=[0.26, 7.68],
        m=[1.57, 0.59],
        a=[-10, 7.53],
        η0=2.11,
        ηinf=0.1,
        β=[3.0, 1.929],
        K=[442, 1289.49],
    ),
    compress=True,
)
view = fem.ViewMaterialIncompressible(
    umat,
    ux=fem.math.linsteps([1, 1.5, 1], num=[200]),
    bx=None,
    ps=None,
    statevars=np.eye(3).reshape(3, 3, 1, 1)
)
view.plot(
    show_kwargs=False,
)

mesh = fem.Cube(n=2)
region = fem.RegionHexahedron(mesh)
field = fem.FieldContainer([fem.Field(region, dim=3)])
boundaries, loadcase = fem.dof.uniaxial(field, clamped=False)

solid = fem.SolidBodyNearlyIncompressible(umat, field, bulk=5000)
solid.results.statevars.T[:] = np.eye(3)

move = fem.math.linsteps([0, 0.5, 0, 0.5], num=200)
step = fem.Step(items=[solid], ramp={boundaries["move"]: move}, boundaries=boundaries)
job = fem.CharacteristicCurve(steps=[step], boundary=boundaries["move"])
job.evaluate(solver=spsolve, tol=1e-2)
fig, ax = job.plot(
    xlabel="Displacement $u$ in mm $\longrightarrow$",
    ylabel="Normal Force $F$ in N $\longrightarrow$",
)
solid.imshow("Principal Values of Cauchy Stress")

[bhaveshshrimali]

Thanks, Andreas! That beta-instability is a great catch.!! This is awesome. My only thought about the constitutive update is that G currently just evaluates everything at t (or C) along with the intermediates for Cv (Cvn, Cvn + k1*dt/2., ....). For the material tangent that should suffice, but the update itself would be slightly inaccurate. I'll think about this more later in the weekend.

Here is how it looks for a comparison against the material behavior for a completely incomrpessible solid (VHB4910 in the linked paper above).

Code

import numpy as np
import felupe as fem
from pypardiso import spsolve
from pandas import read_excel
from matplotlib import pyplot as plt

plt.rc("text", usetex=True)
plt.rc("text.latex", preamble=r"\usepackage{amsmath,amsfonts,amssymb,bm}")

import matadi as mat
from matadi.math import (
    det,
    trace,
    inv,
    gradient,
)

F = mat.Variable("F", 3, 3)
ζ = mat.Variable("ζ", 3, 3)


def ψ(x, dt, μ, α, m, a, η0, ηinf, β, K):
    F, ζn = x[0], x[1]

    # extract old state variables
    Cvn = ζn

    J = det(F)
    C = (F.T @ F) * J ** (-2 / 3)
    I1 = trace(C)

    def invariants(Cv):
        "Invariants of the elastic and viscous parts of the deformation."
        Ce = C @ inv(Cv)
        I1e = trace(Ce)
        I1v = trace(Cv)
        I2e = (I1e**2 - trace(Ce @ Ce)) / 2
        return I1e, I1v, I2e

    # I1e, I1v, I2e = invariants(Cvn)

    def G(Cv):
        "Evolution equation."
        I1e, I1v, I2e = invariants(Cv)
        J2Neq = (I1e**2 / 3 - I2e) * sum(
            [3 ** (1 - a[r]) * m[r] * I1e ** (a[r] - 1) for r in [0, 1]]
        ) ** 2
        u = [3 ** (1 - a[r]) * m[r] * I1e ** (a[r] - 1) for r in [0, 1]]
        ηK = ηinf + (η0 - ηinf + K[0] * (I1v ** β[0] - 3 ** β[0])) / (
            1 + (K[1] * J2Neq + 1e-10) ** β[1]
        )

        return sum(u) / ηK * (C - I1e / 3 * Cv)

    k1 = G(Cvn)
    k2 = G(Cvn + k1 * dt / 2)
    k3 = G(Cvn + (3 * k1 + k2) * dt / 16)
    k4 = G(Cvn + k3 * dt / 2)
    k5 = G(Cvn + 3 * (-k2 + 2 * k3 + 3 * k4) * dt / 16)
    k6 = G(Cvn + (k1 + 4 * k2 + 6 * k3 - 12 * k4 + 8 * k5) * dt / 7)
    Cv = Cvn + dt / 90 * (7 * k1 + 32 * k3 + 12 * k4 + 32 * k5 + 7 * k6)

    I1e, I1v, I2e = invariants(Cv)

    # strain energy functions of the equilibrium and non-equilibrium parts
    p = [0, 1]
    ψEq = [3 ** (1 - α[r]) / (2 * α[r]) * μ[r] * (I1 ** α[r] - 3 ** α[r]) for r in p]
    ψNEq = [3 ** (1 - a[r]) / (2 * a[r]) * m[r] * (I1e ** a[r] - 3 ** a[r]) for r in p]

    return gradient(sum(ψEq) + sum(ψNEq), F), Cv


dt_vals = 0.5

umat = mat.MaterialTensor(
    x=[F, ζ],
    fun=ψ,
    statevars=1,
    kwargs=dict(
        dt=dt_vals,
        μ=[13.54, 1.08],
        α=[1.0, -2.474],
        m=[5.42, 20.78],
        a=[-10.0, 1.948],
        η0=7014.0,
        ηinf=0.1,
        β=[1.852, 0.26],
        K=[3507.0, 1.0],
    ),
    compress=True,
)
# view = fem.ViewMaterialIncompressible(
#     umat,
#     ux=fem.math.linsteps([1, 1.5, 1], num=[200]),
#     bx=None,
#     ps=None,
#     statevars=np.eye(3).reshape(3, 3, 1, 1)
# )
# view.plot(
#     show_kwargs=False,
# )

mesh = fem.Cube(n=3)
region = fem.RegionHexahedron(mesh, quadrature=fem.GaussLegendre(order=3, dim=3))
field = fem.FieldContainer([fem.Field(region, dim=3)])
boundaries, loadcase = fem.dof.uniaxial(field, clamped=False)

solid = fem.SolidBodyNearlyIncompressible(umat, field, bulk=5000)
solid.results.statevars.T[:] = np.eye(3)

# move = fem.math.linsteps([0, 2.0, 0], num=200)

ldot = 0.05
lfinal = 3.0
tfinal = (lfinal - 1.0) / ldot
num_steps = tfinal / dt_vals
move = fem.math.linsteps([0.0, lfinal - 1.0, 0.358], num=int(num_steps) + 1)

step = fem.Step(items=[solid], ramp={boundaries["move"]: move}, boundaries=boundaries)
job = fem.CharacteristicCurve(steps=[step], boundary=boundaries["move"])
job.evaluate(
    solver=spsolve,
    tol=1e-3,
    filename="visco.xdmf",
    newton_verbose=True,
    verbose=True,
    parallel=True,
)
fig, ax = job.plot(
    xlabel="Displacement $u$ in mm $\longrightarrow$",
    ylabel="Normal Force $F$ in N $\longrightarrow$",
)

plt.close()
line = ax.lines[0]

xdata, ydata = line.get_xdata(), line.get_ydata()

data_analytical = read_excel(
    "S11_vs_l11_matrix_RK5_dl.xlsx", index_col=None, header=None
).values

fig, ax3 = plt.subplots(1, 1, figsize=(4, 4))
ax3.plot(1.0 + xdata, ydata, color="blue", label=r"FElupe: $\dot{\lambda} = 0.05$/s")
ax3.plot(
    data_analytical[:, 0],
    data_analytical[:, 1] / 1.0e3,
    color="red",
    linestyle="--",
    label="Material Behavior: Uniaxial",
)
ax3.set_xlabel("$\lambda$", fontsize=14)
ax3.set_ylabel("$S_{un}$", fontsize=14)
ax3.legend(fontsize=14, edgecolor="k", fancybox=False)
ax3.grid(linestyle="--")
fig.tight_layout()
fig.savefig("visco.png")
plt.close()

Comparison with material behavior (solving the ODE for Cv)

[adtzlr]

Great @bhaveshshrimali! 🚀

If I look at Fig. 5b of the paper, the peak stress is near 70 MPa. If I use more increments in FElupe, the peak stress increases to approx. 69 MPa.

[bhaveshshrimali]

This is already looking fantastic. Thanks a lot!

[adtzlr]

Hi @bhaveshshrimali,

I have something for you 🎁 🎁 . I know we've already discussed this question in detail and we found a solution, but I managed to implement the Newton-Rhapson iterations - finally! Unfortunately not with matadi, but with tensortrax. It's possible to use quite large increments with this constitutive model formulation.

The trick is to un-couple the evolution() equation from the AD tensors. Hence, all arguments must be simple arrays. That's why C.x is passed instead of C. Internally, Cvn is a Tensor to evaluate the Jacobian of the evolution equation, but totally isolated from the strain energy density. FElupe's newtonrhapson() function is used to update Cv. The solve-argument is a new math-helper math.solve_2d(A, b) to flatten the second-order tensor Cv and move the batch-dimensions to the front for numpy.linalg.solve.

(Unfortunately I still have no idea how to handle that with matadi.)

Code (this requires FElupe >=8.4.0)

from tensortrax import function, jacobian
from tensortrax.math import trace
from tensortrax.math.linalg import inv
from tensortrax.math.special import dev, from_triu_1d, triu_1d

import felupe as fem


@fem.constitution.isochoric_volumetric_split
def viscoelastic_two_potential(C, ζn, dt, μ, α, m, a, η0, ηinf, β, K):
    # extract old state variables
    Cvn = from_triu_1d(ζn, like=C)

    # first invariant of the right Cauchy-Green deformation tensor
    I1 = trace(C)

    def invariants(Cv, C):
        "Invariants of the elastic and viscous parts of the deformation."
        Ce = C @ inv(Cv)

        I1e = trace(Ce)
        I1v = trace(Cv)
        I2e = (I1e**2 - trace(Ce @ Ce)) / 2

        return I1e, I1v, I2e

    def evolution(Cv, Cvn, C):
        "Evolution equation."
        I1e, I1v, I2e = invariants(Cv, C)
        J2Neq = (I1e**2 / 3 - I2e) * sum(
            [3 ** (1 - a[r]) * m[r] * I1e ** (a[r] - 1) for r in [0, 1]]
        ) ** 2
        u = [3 ** (1 - a[r]) * m[r] * I1e ** (a[r] - 1) for r in [0, 1]]
        ηK = ηinf + (η0 - ηinf + K[0] * (I1v ** β[0] - 3 ** β[0])) / (
            1 + (K[1] * J2Neq + 1e-10) ** β[1]
        )
        G = sum(u) / ηK * dev(C @ inv(Cv)) @ Cv
        return G - (Cv - Cvn) / dt

    # update the state variables
    Cv = fem.newtonrhapson(
        x0=Cvn,
        fun=function(evolution, ntrax=C.ntrax),
        jac=jacobian(evolution, ntrax=C.ntrax),
        solve=fem.math.solve_2d,
        args=(Cvn, C.x),
        verbose=0,
        tol=1e-1,
    ).x

    I1e, I1v, I2e = invariants(Cv, C)

    # strain energy functions of the equilibrium and non-equilibrium parts
    p = [0, 1]
    ψEq = [3 ** (1 - α[r]) / (2 * α[r]) * μ[r] * (I1 ** α[r] - 3 ** α[r]) for r in p]
    ψNEq = [3 ** (1 - a[r]) / (2 * a[r]) * m[r] * (I1e ** a[r] - 3 ** a[r]) for r in p]

    return sum(ψEq) + sum(ψNEq), triu_1d(Cv)


dt_vals = 4.0

umat = fem.Hyperelastic(
    viscoelastic_two_potential,
    nstatevars=6,
    dt=dt_vals,
    μ=[13.54, 1.08],
    α=[1.0, -2.474],
    m=[5.42, 20.78],
    a=[-10, 1.948],
    η0=7014.0,
    ηinf=0.1,
    β=[1.852, 0.26],
    K=[3507.0, 1.0],
    # parallel=True,
)

mesh = fem.Cube(n=2)
region = fem.RegionHexahedron(mesh)
field = fem.FieldContainer([fem.Field(region, dim=3)])
boundaries, loadcase = fem.dof.uniaxial(field, clamped=False)

solid = fem.SolidBodyNearlyIncompressible(umat, field, bulk=5000)
solid.results.statevars[[0, 3, 5]] += 1.0

ldot = 0.05
lfinal = 3.0
tfinal = (lfinal - 1.0) / ldot
num_steps = tfinal / dt_vals
move = fem.math.linsteps([0.0, lfinal - 1.00, 0.358], num=int(num_steps) + 1)

step = fem.Step(items=[solid], ramp={boundaries["move"]: move}, boundaries=boundaries)
job = fem.CharacteristicCurve(steps=[step], boundary=boundaries["move"])
job.evaluate()
fig, ax = job.plot(
    xlabel="Displacement $u$ in mm $\longrightarrow$",
    ylabel="Normal Force $F$ in N $\longrightarrow$",
    marker="x",
)
img = solid.imshow("Principal Values of Cauchy Stress")

If you're interested, we could add this to the examples section in the docs!? It's running fast in that basic setting, so that should be fine for the read-the-docs CI/CD-pipeline. I'd be great if you could open a PR with a new example as examples/ex13_nonlinear-viscoelasticity-newton.py. No need to polish/finish the code, just copy the above code-block and I'll do the rest. I'd just like you to get listed as contributor because you helped a lot on implementing this model!

[bhaveshshrimali]

That's great @adtzlr ! This looks fantastic. Let me take a look in a bit and get started on the example. I appreciate you helping out. FEluPE is in fact a pretty robut and efficient code mechanics and constitutive modeling. Sorry I was traveling the past couple weeks, hence the delay in my response.