Usage of mixed-field (u,p,J,T) material formulations in FElupe

[adtzlr]

Note: This is an advanced example.

The following example demonstrates the thermal shrinkage of a nearly-incompressible block of rubber.

import felupe
import matadi
from matadi.math import det, gradient

F = matadi.Variable("F", 3, 3)
p = matadi.Variable("p")
J = matadi.Variable("J")
T = matadi.Variable("T")
z = matadi.Variable("z") # add empty state variables

def nh(x, C10, bulk, alpha, T0=0):
    r"""Isotropic hyperelastic Neo-Hookean :math:`(u,p,J,T)` - material 
    formulation. Hu-Washizu total potential energy of internal forces given as
    
    ..  math::
        
        \Pi_{\text{int}}(\boldsymbol{u},p,J,T) = 
        \int_V \left( \hat{\psi}(\hat{\boldsymbol{F}}) + U(J_M) 
        + p (\det(F) - J) \right) \ dV
    
    with the multiplicative decomposed mechanical and thermal part of the 
    volume change:
    
    ..  math::
        
        J &= J_M J_T
        
        J_T &= \left( 1 + \alpha (T - T_0) \right)^3
    
    """
    
    # unpack variables
    F, p, J, T, z = x
    
    # thermal part of the volume change
    JT = (1 + alpha * (T - T0)) ** 3
    
    # isochoric and volumetric part of the strain energy function 
    W = (matadi.models.neo_hooke(det(F)**(-1 / 3) * F, C10)
      + bulk * (J / JT - 1) ** 2 / 2
      + p * (det(F) - J)
     )

    return [gradient(W, F), gradient(W, p), gradient(W, J), gradient(W, T), z]

umat = matadi.MaterialTensor(
    x=[F, p, J, T, z],
    fun=nh, 
    kwargs={"C10": 0.5, "bulk": 5000.0, "alpha": 1.5e-4},
    triu=True,
    statevars=1,
)

mesh = felupe.Cube(n=11)
region = felupe.RegionHexahedron(mesh)
field = felupe.FieldsMixed(region, n=4)
solid = felupe.SolidBody(umat, field)

boundaries = felupe.dof.symmetry(field[0])
boundaries["fixed"] = felupe.Boundary(field[0], fx=lambda x: x == 1)
boundaries["temperature"] = felupe.Boundary(
    field=field[-1], 
    fx=lambda x: True, 
    value=-150
)

dof0, dof1 = felupe.dof.partition(field, boundaries)
ext0 = felupe.dof.apply(field, boundaries, dof0)

res = felupe.newtonrhapson(items=[solid], dof0=dof0, dof1=dof1, ext0=ext0)

felupe.tools.save(
    region, 
    field, 
    filename="hyperelastic_cube_thermal_contraction.vtk",
)

Newton-Rhapson solver
=====================

| # | norm(fun) |  norm(dx) |
|---|-----------|-----------|
| 1 | 1.160e+01 | 4.743e+03 |
| 2 | 1.936e-04 | 2.837e-01 |
| 3 | 1.287e-09 | 5.771e-04 |

Solution converged in 3 iterations within 3.099 seconds.

[adtzlr]

Note: This tutorial needs FElupe >= v.2.0.0 (will be released in January 2022).

Note: This tutorial needs FElupe >= v.4.0.0 (will be released in August 2022).

Note: This tutorial needs FElupe >= v.6.0.0.