Mixed-field hyperelasticity with quadratic triangles

A 90° section of a plane-strain circle is subjected to frictionless uniaxial compression by a vertically moved rigid top plate. A mixed-field formulation is used with quadratic triangles.

from functools import partial

import numpy as np

import felupe as fem

A 90° section of a circle with quadratic triangles is created. The midpoints are shifted to the outer radius. An additional point, used as center- (control-) point, is added to the mesh.

mesh = fem.Circle(n=6, sections=[0]).triangulate().add_midpoints_edges()
mask = np.isclose(mesh.x**2 + mesh.y**2, 1, atol=0.05)
mesh.points[mask] /= np.linalg.norm(mesh.points[mask], axis=1).reshape(-1, 1)
mesh.add_points([0, 1.1])
mesh.clear_points_without_cells()

Let’s create a region for quadratic triangles and a mixed-field container with two dual fields, one for the pressure and another one for the volume ratio. The dual fields are disconnected.

region = fem.RegionQuadraticTriangle(mesh)
field = fem.FieldsMixed(
    region, n=3, values=(0.0, 0.0, 1.0), planestrain=True, disconnect=True
)

# create a nearly-incompressible hyperelastic solid body and the rigid top plate
umat = fem.NearlyIncompressible(material=fem.NeoHooke(mu=1), bulk=5000)
solid = fem.SolidBody(umat=umat, field=field)
top = fem.MultiPointContact(
    field=field,
    points=np.arange(mesh.npoints)[np.isclose(mesh.x**2 + mesh.y**2, 1)],
    centerpoint=-1,
    skip=(1, 0),
)
mesh.plot(nonlinear_subdivision=4, plotter=top.plot(line_width=5, opacity=1)).show()

Static Scene

Interactive Scene

/home/docs/checkouts/readthedocs.org/user_builds/felupe/envs/stable/lib/python3.12/site-packages/felupe/view/_scene.py:251: PyVistaFutureWarning: The default value of `algorithm` for the filter
`UnstructuredGrid.extract_surface` will change in the future. It currently defaults to
`'dataset_surface'`, but will change to `None`. Explicitly set the `algorithm` keyword to
silence this warning.
  surface = surface.extract_surface(
/home/docs/checkouts/readthedocs.org/user_builds/felupe/envs/stable/lib/python3.12/site-packages/felupe/view/_scene.py:290: PyVistaFutureWarning: The default value of `algorithm` for the filter
`UnstructuredGrid.extract_surface` will change in the future. It currently defaults to
`'dataset_surface'`, but will change to `None`. Explicitly set the `algorithm` keyword to
silence this warning.
  .extract_surface(nonlinear_subdivision=nonlinear_subdivision)

A step is used containts the solid body and the rigid top plate as items. The rigid vertical movement of the top plate is applied in a ramped manner.

boundaries = fem.dof.symmetry(field[0])
boundaries["move"] = fem.Boundary(field[0], fy=1.1, skip=(1, 0))
move = fem.math.linsteps([0, -0.4], num=4)
ramp = {boundaries["move"]: move}
step = fem.Step(items=[solid, top], ramp=ramp, boundaries=boundaries)
job = fem.Job(steps=[step]).evaluate()

The maximum principal values of the Cauchy stress tensor are plotted. The cell-based means are projected to the mesh-points.

solid.plot(
    "Principal Values of Cauchy Stress",
    nonlinear_subdivision=4,
    plotter=top.plot(line_width=5, opacity=1),
    project=partial(fem.project, mean=True),
).show()

Static Scene

Interactive Scene

/home/docs/checkouts/readthedocs.org/user_builds/felupe/envs/stable/lib/python3.12/site-packages/felupe/view/_scene.py:251: PyVistaFutureWarning: The default value of `algorithm` for the filter
`UnstructuredGrid.extract_surface` will change in the future. It currently defaults to
`'dataset_surface'`, but will change to `None`. Explicitly set the `algorithm` keyword to
silence this warning.
  surface = surface.extract_surface(
/home/docs/checkouts/readthedocs.org/user_builds/felupe/envs/stable/lib/python3.12/site-packages/felupe/view/_scene.py:290: PyVistaFutureWarning: The default value of `algorithm` for the filter
`UnstructuredGrid.extract_surface` will change in the future. It currently defaults to
`'dataset_surface'`, but will change to `None`. Explicitly set the `algorithm` keyword to
silence this warning.
  .extract_surface(nonlinear_subdivision=nonlinear_subdivision)