Viscoelasticity with local Newton-iterations

[adtzlr]

As a follow-up on #666, this is finite-strain viscoelasticity (probably the most simple model formulation) with local Newton-iterations to solve the evolution equation

• using tensortrax.
• using matadi.

This could enable e.g. an implementation of the Bergstrom-Boyce model?

It's basically straigh-forward

1. define the evolution equation,
2. create the function + jacobian,
3. use a solver which places the batch-axis according to numpy.linalg.solve and
4. use the Newton-Rhapson solver inside the strain-energy function.

import numpy as np
import felupe as fem

import tensortrax as tr
from tensortrax.math import trace
from tensortrax.math.special import from_triu_1d, triu_1d, dev
from tensortrax.math.linalg import inv

mesh = fem.Rectangle(n=6)
region = fem.RegionQuad(mesh)
field = fem.FieldContainer([fem.FieldPlaneStrain(region, dim=2)])

boundaries, loadcase = fem.dof.uniaxial(field, clamped=True)


def local_solve(A, b):
    B = np.einsum("...i->i...", b.reshape(9, np.prod(b.shape[2:])))
    a = np.einsum("ij...->...ij", A.reshape(9, 9, np.prod(A.shape[4:])))
    return np.einsum("...i->i...", np.linalg.solve(a, B)).reshape(*b.shape)


@fem.constitution.isochoric_volumetric_split
def finite_strain_viscoelastic(C, Cin, mu, eta, dtime):
    # update of state variables by evolution equation
    Cin = from_triu_1d(Cin, like=C)

    def evolution(Ci, Cin, C, mu, eta, dtime):
        return mu / eta * dev(C @ inv(Ci)) @ Ci - (Ci - Cin) / dtime

    fun = tr.function(evolution, ntrax=C.ntrax)
    jac = tr.jacobian(evolution, ntrax=C.ntrax)

    Ci = fem.newtonrhapson(
        Cin, fun, jac, local_solve, args=(Cin, C.x, mu, eta, dtime), verbose=0
    ).x

    # first invariant of elastic part of right Cauchy-Green deformation tensor
    I1 = trace(C @ inv(Ci))

    # strain energy function and state variable
    return mu / 2 * (I1 - 3), triu_1d(Ci)


umat = fem.Hyperelastic(
    finite_strain_viscoelastic, mu=1.0, eta=1.0, dtime=0.1, nstatevars=6
)
solid = fem.SolidBodyNearlyIncompressible(umat, field, bulk=5000)
solid.results.statevars[[0, 3, 5]] += 1

move = fem.math.linsteps([0, 0.3, 0], num=11)
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 $d_1$ in mm $\longrightarrow$",
    ylabel="Normal Force $F_1$ in N $\longrightarrow$",
)
solid.plot("Principal Values of Cauchy Stress").show()

[adtzlr]

see also #742 and the newly added test:

felupe/tests/test_constitution_newton.py

Lines 42 to 66 in bd2d6fc

	@fem.constitution.isochoric_volumetric_split
	def finite_strain_viscoelastic_newton(C, Cin, mu, eta, dtime):
	"Finite Strain Viscoelastic material formulation."
	def evolution(Ci, Cin, C, mu, eta, dtime):
	"Viscoelastic evolution equation."
	return mu / eta * dev(C @ inv(Ci)) @ Ci - (Ci - Cin) / dtime
	# update of state variables by evolution equation
	Cin = from_triu_1d(Cin, like=C)
	Ci = fem.newtonrhapson(
	x0=Cin,
	fun=tr.function(evolution, ntrax=C.ntrax),
	jac=tr.jacobian(evolution, ntrax=C.ntrax),
	solve=fem.math.solve_2d,
	args=(Cin, C.x, mu, eta, dtime),
	verbose=0,
	tol=1e-2,
	).x
	# first invariant of elastic part of right Cauchy-Green deformation tensor
	I1 = trace(C @ inv(Ci))
	# strain energy function and state variable
	return mu / 2 * (I1 - 3), triu_1d(Ci)

[adtzlr]

Now that the local Newton-iterations are working as expected, let's try to implement the viscoelastic microsphere model. This is really nice because the evolution equation is derived by the constraint equation. Newton's method is applied independently on the affine micro-stretches. Hence, the nd-solver is called with n=0 and the micro-stretches are treated as additional elementwise-operating trailing axis.

import felupe as fem
import tensortrax.math as tm
import tensortrax as tr
from functools import partial


def visco(C, statevars, **kwargs):

    def W(λ, ζn, μ, τ, δ, Δt):
        εn = tm.array(ζn[:21], like=λ)

        def evolution(ε, εn, λ):
            ψ = μ / 2 * (tm.log(λ) - ε) ** 2
            ϕ = δ / (τ * μ * (1 + δ)) * (τ * μ * tm.abs(ε - εn) / Δt) ** ((1 + δ) / δ)
            return ψ + ϕ

        ε = fem.newtonrhapson(
            x0=εn.x,
            fun=tr.jacobian(evolution, ntrax=1 + λ.ntrax),
            jac=tr.hessian(evolution, ntrax=1 + λ.ntrax),
            solve=partial(fem.math.solve_nd, n=0),
            args=(εn.x, λ.x),
            verbose=0,
        ).x
        ψ = μ / 2 * (tm.log(λ) - ε) ** 2
        return ψ, ε

    ψf, statevars1 = (
        fem.constitution.hyperelasticity.models.microsphere.affine_stretch_statevars(
            C, statevars[:21], f=W, kwargs=kwargs
        )
    )

    ψc, statevars2 = (
        fem.constitution.hyperelasticity.models.microsphere.affine_tube_statevars(
            C, statevars[21:], f=W, kwargs=kwargs
        )
    )

    return ψf + ψc, tm.stack([*statevars1, *statevars2])


umat = fem.MaterialAD(
    visco,
    **{"μ": 6, "τ": 0.1, "δ": 1, "Δt": 10.0},
    nstatevars=42,
    parallel=True,
)
ux = fem.math.linsteps(
    [1, 1.5, 1],
    num=20,
)
ax = umat.plot(
    incompressible=True,
    ux=ux,
    ps=ux,
    bx=(ux - 1) / 2 + 1,
)

[adtzlr]

Sadly, with the parameters from the paper, Newton's method fails.

import felupe as fem
import tensortrax.math as tm
import tensortrax as tr
import pypardiso
from functools import partial


def visco(C, statevars, **kwargs):

    def W(λ, ζn, μ, τ, δ, Δt):
        εn = tm.array(ζn, like=λ)

        def evolution(ε, εn, λ):
            ψ = μ / 2 * (tm.log(λ) - ε) ** 2
            ϕ = δ / (τ * μ * (1 + δ)) * (τ * μ * tm.abs(ε - εn) / Δt) ** ((1 + δ) / δ)
            return ψ + ϕ * Δt

        ε = fem.newtonrhapson(
            x0=εn.x,
            fun=tr.jacobian(evolution, ntrax=1 + λ.ntrax),
            jac=tr.hessian(evolution, ntrax=1 + λ.ntrax),
            solve=partial(fem.math.solve_nd, n=0),
            args=(εn.x, λ.x),
            verbose=0,
            tol=1e-2,
        ).x
        ψ = μ / 2 * (tm.log(λ) - ε) ** 2
        return ψ, ε

    w1 = fem.constitution.hyperelasticity.models.microsphere.affine_stretch_statevars
    w2 = fem.constitution.hyperelasticity.models.microsphere.affine_tube_statevars
    dt = 0.1
    ψf1, sv1 = w1(
        C, statevars[:21], f=W, kwargs={"μ": 6.16, "τ": 0.1, "δ": 2.97, "Δt": dt}
    )
    ψf2, sv2 = w1(
        C, statevars[21:42], f=W, kwargs={"μ": 1.35, "τ": 1e1, "δ": 2.46, "Δt": dt}
    )
    ψf3, sv3 = w1(
        C, statevars[42:63], f=W, kwargs={"μ": 1.76, "τ": 1e3, "δ": 24.46, "Δt": dt}
    )
    ψc1, sv4 = w2(
        C, statevars[63:84], f=W, kwargs={"μ": 7.19, "τ": 0.1, "δ": 1.02, "Δt": dt}
    )
    ψc2, sv5 = w2(
        C, statevars[84:105], f=W, kwargs={"μ": 1.13, "τ": 1e1, "δ": 2.05, "Δt": dt}
    )
    ψc3, sv6 = w2(
        C, statevars[105:126], f=W, kwargs={"μ": 1.63, "τ": 1e3, "δ": 16.47, "Δt": dt}
    )

    return ψf1 + ψf2 + ψf3 + ψc1 + ψc2 + ψc3, tm.stack(
        [*sv1, *sv2, *sv3, *sv4, *sv5, *sv6]
    )


umat = fem.MaterialAD(
    visco,
    nstatevars=126,
    parallel=True,
)
ux = fem.math.linsteps(
    [1.01, 1.5, 1],
    num=100,
)
ax = umat.plot(
    incompressible=True,
    ux=ux,
    ps=ux,
    bx=(ux - 1) / 2 + 1,
)