[FElupe v0.1.1] A Cube (Arbitrary Order Lagrange Hexahedron) with Isotropic Hyperelasticity under Uniaxial Compression

[adtzlr]

This is an example with an arbitrary order lagrange hexahedron element in combination with isotropic hyperelasticity under uniaxial compression, glued to a rigid end plate.

As always, import FElupe (and matadi).

import felupe
import matadi

First, we define a plot function (will be used later).

import numpy as np
import matplotlib.pyplot as plt

def plot(mesh, field=None):
    "Plot deformed mesh."

    # init figure
    fig, ax = plt.subplots(subplot_kw={"projection": "3d"})

    # undeformed coordinates
    x = mesh.points

    # deformed coordinates
    if field is not None:
        x += field.values

    # reshape coordinates for surface representation
    xx = x[:, 0].reshape(*([order + 1] * 3))
    yy = x[:, 1].reshape(*([order + 1] * 3))
    zz = x[:, 2].reshape(*([order + 1] * 3))

    # surface properties
    surface_props = {"color": "1.0", "edgecolor": "k", "linewidth": 1.5}

    # plot surfaces
    ax.plot_surface(xx[0, :, :], yy[0, :, :], zz[0, :, :], **surface_props)
    ax.plot_surface(xx[-1, :, :], yy[-1, :, :], zz[-1, :, :], **surface_props)
    #
    ax.plot_surface(xx[:, 0, :], yy[:, 0, :], zz[:, 0, :], **surface_props)
    ax.plot_surface(xx[:, -1, :], yy[:, -1, :], zz[:, -1, :], **surface_props)
    #
    ax.plot_surface(xx[:, :, 0], yy[:, :, 0], zz[:, :, 0], **surface_props)
    ax.plot_surface(xx[:, :, -1], yy[:, :, -1], zz[:, :, -1], **surface_props)

    # set limits
    b = max([max(mesh.points[:, a]) for a in range(3)])
    ax.set_xlim(0, b)
    ax.set_ylim(0, b)
    ax.set_zlim(0, b)

    # turn off grid and axis
    ax.grid(False)
    ax.axis("off")

    return fig, ax

The parameters of this example are the order of the lagrange polynomials and the length, width and height of the mesh.

order = 3

L = 4
B = 4
H = 1

Next, we create our mesh (in this case the points of one element):

# init mesh
xp = np.linspace(0, H / 2, order + 1)
yp = np.linspace(0, B / 2, order + 1)
zp = np.linspace(0, L / 2, order + 1)
xpv = [xp, yp, zp]
xv = np.meshgrid(*xpv[::-1], indexing="ij")

points = np.vstack([x.ravel() for x in xv[::-1]]).T
cells = np.arange(len(points)).reshape(1, -1)

mesh = felupe.Mesh(points, cells)

# plot undeformed element
fig, ax = plot(mesh)

This will be followed by setting up the problem

# setup problem
element = felupe.ArbitraryOrderLagrangeElement(order=order, ndim=3)
quadrature = felupe.GaussLegendre(order=order, dim=3)
region = felupe.Region(mesh, element, quadrature)
field = felupe.Field(region, dim=3)

and the definition of the material behavior:

# material definition
umat = matadi.MaterialHyperelastic(
    matadi.models.third_order_deformation,
    C10=0.3,
    C01=0.2,
    C11=0.02,
    C20=-0.05,
    C30=0.015,
    bulk=50,
)

A Newton-Rhapson scheme is used to obtain the state of equilibrium. Due to the nearly incompressible material parameters small increments of external prescribed displacements have to be applied.

# external displacement W/H
WH = np.linspace(0, -0.2, 11)

for increment, move in enumerate(WH[1:] * H / 2):

    # loadcase
    boundaries, dof0, dof1, u0ext = felupe.dof.uniaxial(field, right=H / 2, move=move)

    for iteration in range(8):

        F = field.extract()
        P = umat.gradient([F])[0]
        A = umat.hessian([F])[0]

        linearform = felupe.IntegralForm(P, field, region.dV, grad_v=True)
        bilinearform = felupe.IntegralForm(
            A, field, region.dV, field, grad_v=True, grad_u=True
        )

        r = linearform.assemble(parallel=True).toarray()[:, 0]
        K = bilinearform.assemble(parallel=True)

        system = felupe.solve.partition(field, K, dof1, dof0, r)
        dfield = felupe.solve.solve(*system, u0ext)

        norm = felupe.math.norm(dfield)
        print(
            "Inc. %d," % (1 + increment),
            "u_ext %1.2f," % (2 * move),
            "iter. %d," % (1 + iteration),
            "norm %1.2e" % norm,
        )

        field += dfield

        if np.sum(norm) < 1e-10:
            break

    print(
        "\nReaction Force X = %1.4f\n\n" % (4 * r[boundaries["move"].dof].sum() / L / B)
    )

The residuals per increment and iteration are printed out.

Inc. 1, u_ext -0.02, iter. 1, norm 2.35e-01
Inc. 1, u_ext -0.02, iter. 2, norm 2.58e-03
Inc. 1, u_ext -0.02, iter. 3, norm 8.32e-06
Inc. 1, u_ext -0.02, iter. 4, norm 7.71e-11

Reaction Force X = -0.1855


Inc. 2, u_ext -0.04, iter. 1, norm 2.39e-01
Inc. 2, u_ext -0.04, iter. 2, norm 2.67e-03
Inc. 2, u_ext -0.04, iter. 3, norm 9.49e-06
Inc. 2, u_ext -0.04, iter. 4, norm 1.08e-10
Inc. 2, u_ext -0.04, iter. 5, norm 5.64e-16

Reaction Force X = -0.3864


Inc. 3, u_ext -0.06, iter. 1, norm 2.44e-01
Inc. 3, u_ext -0.06, iter. 2, norm 2.76e-03
Inc. 3, u_ext -0.06, iter. 3, norm 1.10e-05
Inc. 3, u_ext -0.06, iter. 4, norm 1.54e-10
Inc. 3, u_ext -0.06, iter. 5, norm 5.91e-16

Reaction Force X = -0.6033


Inc. 4, u_ext -0.08, iter. 1, norm 2.48e-01
Inc. 4, u_ext -0.08, iter. 2, norm 2.84e-03
Inc. 4, u_ext -0.08, iter. 3, norm 1.34e-05
Inc. 4, u_ext -0.08, iter. 4, norm 2.46e-10
Inc. 4, u_ext -0.08, iter. 5, norm 6.20e-16

Reaction Force X = -0.8375


Inc. 5, u_ext -0.10, iter. 1, norm 2.52e-01
Inc. 5, u_ext -0.10, iter. 2, norm 2.91e-03
Inc. 5, u_ext -0.10, iter. 3, norm 1.79e-05
Inc. 5, u_ext -0.10, iter. 4, norm 4.88e-10
Inc. 5, u_ext -0.10, iter. 5, norm 1.01e-15

Reaction Force X = -1.0903


Inc. 6, u_ext -0.12, iter. 1, norm 2.57e-01
Inc. 6, u_ext -0.12, iter. 2, norm 2.99e-03
Inc. 6, u_ext -0.12, iter. 3, norm 2.74e-05
Inc. 6, u_ext -0.12, iter. 4, norm 1.27e-09
Inc. 6, u_ext -0.12, iter. 5, norm 7.89e-16

Reaction Force X = -1.3639


Inc. 7, u_ext -0.14, iter. 1, norm 2.62e-01
Inc. 7, u_ext -0.14, iter. 2, norm 3.15e-03
Inc. 7, u_ext -0.14, iter. 3, norm 4.52e-05
Inc. 7, u_ext -0.14, iter. 4, norm 3.93e-09
Inc. 7, u_ext -0.14, iter. 5, norm 1.06e-15

Reaction Force X = -1.6613


Inc. 8, u_ext -0.16, iter. 1, norm 2.66e-01
Inc. 8, u_ext -0.16, iter. 2, norm 3.40e-03
Inc. 8, u_ext -0.16, iter. 3, norm 7.38e-05
Inc. 8, u_ext -0.16, iter. 4, norm 1.18e-08
Inc. 8, u_ext -0.16, iter. 5, norm 2.11e-15

Reaction Force X = -1.9867


Inc. 9, u_ext -0.18, iter. 1, norm 2.70e-01
Inc. 9, u_ext -0.18, iter. 2, norm 3.75e-03
Inc. 9, u_ext -0.18, iter. 3, norm 1.14e-04
Inc. 9, u_ext -0.18, iter. 4, norm 3.11e-08
Inc. 9, u_ext -0.18, iter. 5, norm 1.16e-14

Reaction Force X = -2.3446


Inc. 10, u_ext -0.20, iter. 1, norm 2.74e-01
Inc. 10, u_ext -0.20, iter. 2, norm 4.25e-03
Inc. 10, u_ext -0.20, iter. 3, norm 1.65e-04
Inc. 10, u_ext -0.20, iter. 4, norm 7.14e-08
Inc. 10, u_ext -0.20, iter. 5, norm 5.67e-14

Reaction Force X = -2.7408

Finally, the deformed body is plotted.

[adtzlr]

Hint: VTK export is not supported for arbitrary order lagrange elements due to different cell-point connectivity between FElupe and VTK.