[FElupe v0.1.1] - Axisymmetric Analysis with (u,p,J) Mixed-Formulation and Isotropic Hyperelasticity

[adtzlr]

This is an axisymmetric example of a mixed field formulation for the treatment of nearly-incompressibility. Matadi is used for the definition of the material behavior. The (Nastran) mesh file is located here - download and save it as conical-spring.bdf.

In the first step, the mesh is read with the help of meshio. A FElupe mesh is created with the imported points and cells.

import felupe
import meshio

conical_spring = meshio.read("conical-spring.bdf")
points = conical_spring.points
cell_type, cells = conical_spring.cells[0]

mesh = felupe.Mesh(points[:, :2], cells, cell_type)

# define cells with steel
import numpy as np

cells_steel = np.concatenate((np.arange(32), np.arange(32, 64), np.arange(64, 96)))

Next, we set up our problem.

# elements, quadrature and regions
element = felupe.Quad()
element0 = felupe.ConstantQuad()
quadrature = felupe.GaussLegendre(order=1, dim=2)
region = felupe.Region(mesh, element, quadrature)
region0 = felupe.Region(felupe.mesh.convert(mesh, order=0), element0, quadrature)

# fields
displacement = felupe.FieldAxisymmetric(region)
pressure = felupe.Field(region0)
volumeratio = felupe.Field(region0, values=1)
fields = felupe.FieldMixed((displacement, pressure, volumeratio))

# boundary conditions
fb = lambda x: x == 0.5
ft = lambda x: x == 1.5

boundaries = {}
boundaries["bottom"] = felupe.Boundary(displacement, fy=fb)
boundaries["top"] = felupe.Boundary(displacement, fy=ft, skip=(1, 0))
boundaries["move"] = felupe.Boundary(displacement, fy=ft, skip=(0, 1), value=0.2)

dof0, dof1, offsets = felupe.dof.partition(fields, boundaries)
u0ext = felupe.dof.apply(displacement, boundaries, dof0)

The original mesh is now rotated. In order to update the regions and fields, we re-define them once again.

# rotate mesh
mesh.points = felupe.mesh.rotate(
    (points, cells), angle_deg=20, axis=2, center=[-1, 0.5, 0]
)[0][:, :2]


# re-setup problem with rotated mesh
region = felupe.Region(mesh, element, quadrature)
region0 = felupe.Region(felupe.mesh.convert(mesh, order=0), element0, quadrature)
displacement = felupe.FieldAxisymmetric(region)
pressure = felupe.Field(region0)
volumeratio = felupe.Field(region0, values=1)
fields = felupe.FieldMixed((displacement, pressure, volumeratio))

# differential areas
dA = region.dV

The constitutive (hyperelastic) material behavior is defined by matadi:

from matadi import ThreeFieldVariation, MaterialHyperelastic
from matadi.math import det, trace, transpose

def W(F, C10, bulk):
    J = det(F)
    C = J ** (-2 / 3) * transpose(F) @ F
    return 2 * C10 * (trace(C) - 3) + bulk * (J - 1) ** 2 / 2

umat = ThreeFieldVariation(MaterialHyperelastic(W, C10=0.5, bulk=5000))
umat_steel = ThreeFieldVariation(MaterialHyperelastic(W, C10=40000, bulk=170000))

Finally, a newton-rhapson - procedure obtains equilibrium on all active dof's:

for iteration in range(8):

    F, p, J = fields.extract()

    P = umat.gradient([F, p, J])
    A = umat.hessian([F, p, J])

    for a in range(len(P)):
        P[a].T[cells_steel] = umat_steel.gradient([F, p, J])[a].T[cells_steel]

    for a in range(len(A)):
        A[a].T[cells_steel] = umat_steel.hessian([F, p, J])[a].T[cells_steel]

    A[1] = A[1][:, :, 0, 0]
    A[2] = A[2][:, :, 0, 0]
    A[3] = A[3][0, 0, 0, 0]
    A[4] = A[4][0, 0, 0, 0]
    A[5] = A[5][0, 0, 0, 0]

    linearform = felupe.IntegralFormMixed(P, fields, dA)
    bilinearform = felupe.IntegralFormMixed(A, fields, dA, fields)

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

    system = felupe.solve.partition(fields, K, dof1, dof0, r)
    dfields = np.split(felupe.solve.solve(*system, u0ext), offsets)

    norm = felupe.math.norm(dfields)
    print(iteration, norm)

    fields += dfields

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

Results are exported via meshio. Steel stress results are set to zero before export in order to visualize the stress state in the rubber. Note that the stress on the steel-rubber boundary is not true because it is divided by the number of cells per point.

P[0].T[cells_steel] = 0
felupe.tools.save(region, fields, r, F, P, unstack=offsets, filename="result.vtk")

Output of residuals:

0 [3.42647231e+00 2.83755098e+02 4.26773377e-02]
1 [1.57197214e-01 8.68606911e+01 8.22748048e-03]
2 [2.90977175e-02 5.34720052e+00 7.81542361e-04]
3 [1.11966963e-02 1.79743200e+00 3.15676938e-04]
4 [1.77519530e-03 3.24537239e-01 5.38398227e-05]
5 [7.40332471e-05 1.37811226e-02 1.98803111e-06]
6 [1.65640634e-07 3.08105018e-05 5.55601698e-09]

[adtzlr]

The problem in the above is example is that FElupe v0.1.1 does not understand the hessian shape of madati (u,p,J) materials.

This is the reason for the following lines of the above code in the newton-scheme:

A[1] = A[1][:, :, 0, 0]
A[2] = A[2][:, :, 0, 0]
A[3] = A[3][0, 0, 0, 0]
A[4] = A[4][0, 0, 0, 0]
A[5] = A[5][0, 0, 0, 0]

FElupe is expecting a shape of (3,3,...) for the off-diagonal parts d2WdFdp, d2WdFdJ whereas matadi produces (3,3,1,1,...). In future versions of FElupe the input should be detected and adopted accordingly.

[adtzlr]

By going through this old example it is really nice how FElupe has evolved since the last two years. I'll update this example and eventually add it to the docs.

[adtzlr]

This is a first idea for a new example with a dynamic mesh.

import felupe as fem
import numpy as np

layers = [2, 11, 2, 11, 2]
lines = [
    fem.mesh.Line(a=0, b=13, n=21).expand(n=1).translate(2, axis=1),
    fem.mesh.Line(a=0, b=13, n=21).expand(n=1).translate(2.5, axis=1),
    fem.mesh.Line(a=-0.2, b=10, n=21).expand(n=1).translate(4.5, axis=1),
    fem.mesh.Line(a=-0.2, b=10, n=21).expand(n=1).translate(5, axis=1),
    fem.mesh.Line(a=-0.4, b=7, n=21).expand(n=1).translate(6.5, axis=1),
    fem.mesh.Line(a=-0.4, b=7, n=21).expand(n=1).translate(7, axis=1),
]
faces = fem.MeshContainer(
    [
        first.fill_between(second, n=n)
        for first, second, n in zip(lines[:-1], lines[1:], layers)
    ]
)
point = lambda m: m.points[np.unique(m.cells)].mean(axis=0) - np.array([2, 0])
mask = lambda m: np.unique(m.cells)
kwargs = dict(axis=1, exponent=2.5, normalize=True)
faces.points[:] = (
    faces[1]
    .add_runouts([3], centerpoint=point(faces[1]), mask=mask(faces[1]), **kwargs)
    .points
)
faces.points[:] = (
    faces[3]
    .add_runouts([7], centerpoint=point(faces[3]), mask=mask(faces[3]), **kwargs)
    .points
)
faces.points[:21] = faces.points[21 : 2 * 21]
faces.points[-21:] = faces.points[-2 * 21 : -21]
faces.points[:] = faces[0].rotate(15, axis=2).points
faces[0].y[:21] = 2
faces[-1].y[-21:] = 8.5
mesh = fem.MeshContainer([faces.stack([1, 3]), faces.stack([0, 2, 4])])
mesh.imshow(colors=[None, "white"], show_edges=True, opacity=1)

mesh = fem.MeshContainer(
    [m.revolve(n=73, phi=270).rotate(-90, axis=1).rotate(-90, axis=2) for m in mesh]
)
mesh.imshow(colors=[None, "white"], show_edges=True, opacity=1)

[adtzlr]

Unfortunately, Read-the-docs can't handle such examples. The only way would be to share them tested with Github Actions in a separate repo? Locally, it takes less than a minute to complete.

code

import felupe as fem
import numpy as np
import pypardiso

layers = [2, 11, 2, 11, 2]
lines = [
    fem.mesh.Line(a=0, b=13, n=21).expand(n=1).translate(2, axis=1),
    fem.mesh.Line(a=0, b=13, n=21).expand(n=1).translate(2.5, axis=1),
    fem.mesh.Line(a=-0.2, b=10, n=21).expand(n=1).translate(4.5, axis=1),
    fem.mesh.Line(a=-0.2, b=10, n=21).expand(n=1).translate(5, axis=1),
    fem.mesh.Line(a=-0.4, b=7, n=21).expand(n=1).translate(6.5, axis=1),
    fem.mesh.Line(a=-0.4, b=7, n=21).expand(n=1).translate(7, axis=1),
]
faces = fem.MeshContainer(
    [
        first.fill_between(second, n=n)
        for first, second, n in zip(lines[:-1], lines[1:], layers)
    ]
)
point = lambda m: m.points[np.unique(m.cells)].mean(axis=0) - np.array([2, 0])
mask = lambda m: np.unique(m.cells)
kwargs = dict(axis=1, exponent=2.5, normalize=True)
faces.points[:] = (
    faces[1]
    .add_runouts([3], centerpoint=point(faces[1]), mask=mask(faces[1]), **kwargs)
    .points
)
faces.points[:] = (
    faces[3]
    .add_runouts([7], centerpoint=point(faces[3]), mask=mask(faces[3]), **kwargs)
    .points
)
faces.points[:21] = faces.points[21 : 2 * 21]
faces.points[-21:] = faces.points[-2 * 21 : -21]
faces.points[:] = faces[0].rotate(15, axis=2).points
faces[0].y[:21] = 2
faces[-1].y[-21:] = 8.5
container = fem.MeshContainer([faces.stack([1, 3]), faces.stack([0, 2, 4])])
container.imshow(colors=[None, "white"], show_edges=True, opacity=1)

container = fem.MeshContainer(
    [
        m.revolve(n=37, phi=180).rotate(-90, axis=1).rotate(-90, axis=2)
        for m in container
    ],
    merge=True,
    decimals=4,
)
container.imshow(colors=[None, "white"], show_edges=True, opacity=1)

mesh = container.stack()
region = fem.RegionHexahedron(mesh)
field = fem.FieldContainer([fem.Field(region, dim=3)])

regions = [fem.RegionHexahedron(m) for m in container]
fields = [fem.FieldsMixed(r, n=1) for r in regions]

boundaries, loadcase = fem.dof.uniaxial(field, clamped=True, axis=2, sym=(0, 1, 0))
solids = [
    fem.SolidBodyNearlyIncompressible(fem.NeoHooke(mu=1), fields[0], bulk=200),
    fem.SolidBody(fem.LinearElasticLargeStrain(E=2.1e5, nu=0.3), fields[1]),
]

move = fem.math.linsteps([0, -1.5], num=3)
ramp = {boundaries["move"]: move}
step = fem.Step(
    items=[
        *solids,
    ],
    ramp=ramp,
    boundaries=boundaries,
)

job = fem.Job(steps=[step])
job.evaluate(x0=field, tol=1e-1, solver=pypardiso.spsolve, parallel=True)

boundaries, loadcase = fem.dof.shear(
    field,
    axes=(0, 2),
    moves=(0, 0, -1.5),
    sym=True,
)

move = fem.math.linsteps([0, 0.5], num=5)
ramp = {boundaries["move"]: move}
step = fem.Step(
    items=[
        *solids,
    ],
    ramp=ramp,
    boundaries=boundaries,
)

job = fem.Job(steps=[step])
job.evaluate(x0=field, tol=1e-2, solver=pypardiso.spsolve, parallel=True)

plotter = fields[1].plot(show_undeformed=False, color="white", show_edges=False)
fields[0].plot(
    "Principal Values of Logarithmic Strain",
    show_undeformed=False,
    project=fem.topoints,
    plotter=plotter,
).show()