Note
Go to the end to download the full example code.
Inflation of a hyperelastic balloon#
This example requires external packages.
pip install contique
With the help of contique it is possible to apply a numerical parameter continuation algorithm on any system of equilibrium equations. This advanced tutorial demonstrates the usage of FElupe in conjunction with contique. The unstable inflation of a rectangular hyperelastic balloon demonstrates this powerful approach. The deformed model and the pressure - displacement curve is plotted.
First, setup a problem in FElupe as usual (mesh, region, field, boundaries, umat, solid
and a pressure boundary). For the material definition we use the
Neo-Hooke built-in
hyperelastic material formulation, see Eq. (1).
(1)#\[\psi = \frac{\mu}{2} \left(
\text{tr}(\boldsymbol{C}) - \ln(\det(\boldsymbol{C}))
\right)\]
import contique
import numpy as np
import felupe as fem
mesh = fem.Rectangle(b=(1, 25), n=(2, 4))
region = fem.RegionQuad(mesh)
field = fem.FieldContainer([fem.FieldAxisymmetric(region, dim=2)])
boundaries = fem.dof.symmetry(field[0], axes=(0, 1))
boundaries["fix-y"] = fem.Boundary(field[0], fy=mesh.y.max(), mode="or", skip=(0, 1))
dof0, dof1 = fem.dof.partition(field, boundaries)
umat = fem.NeoHookeCompressible(mu=1)
solid = fem.SolidBodyNearlyIncompressible(umat, field, bulk=5000)
region_for_pressure = fem.RegionQuadBoundary(mesh, mask=(mesh.x == 0), ensure_3d=True)
field_for_pressure = fem.FieldContainer(
[fem.FieldAxisymmetric(region_for_pressure, dim=2)]
)
pressure = fem.SolidBodyPressure(field_for_pressure)
The next step involves the problem definition for contique. For details have a look at its README.
def fun(x, lpf, *args):
"The system vector of equilibrium equations."
field[0].values.ravel()[dof1] = x
pressure.update(lpf)
return fem.tools.fun([solid, pressure], field)[dof1]
def dfundx(x, lpf, *args):
"""The jacobian of the system vector of equilibrium equations w.r.t. the
primary unknowns."""
K = fem.tools.jac([solid, pressure], field)
return fem.solve.partition(field, K, dof1, dof0)[2]
def dfundl(x, lpf, *args):
"""The jacobian of the system vector of equilibrium equations w.r.t. the
load proportionality factor."""
pressure.update(1)
return fem.tools.fun([pressure], field)[dof1]
Next we have to init the problem and specify the initial values of unknowns (the undeformed configuration). After each completed step of the numeric continuation the results are saved.
|Step,C.| Control Component | Norm (Iter.#) | Message |
|-------|-------------------|---------------|-------------|
| 1,1 | 0+ => 0+ | 4.8e-06 ( 7#) | |
| 2,1 | 0+ => 0+ | 3.9e-06 ( 7#) | |
| 3,1 | 0+ => 0+ | 6.4e-07 ( 7#) | |
| 4,1 | 0+ => 0+ | 8.9e-05 ( 8#) | |
| 5,1 | 0+ => 0+ | 1.6e-05 ( 7#) | |
| 6,1 | 0+ => 0+ | 2.9e-05 ( 7#) | |
| 7,1 | 0+ => 0+ | 3.9e-06 ( 7#) | |
| 8,1 | 0+ => 0+ | 2.9e-06 ( 8#) | |
| 9,1 | 0+ => 0+ | 1.3e+00 ( 8#) |Failed |
| 10,1 | 0+ => 0+ | 3.4e-05 ( 8#) | |
| 11,1 | 0+ => 0+ | 5.3e-02 ( 8#) |Failed |
/home/docs/checkouts/readthedocs.org/user_builds/felupe/envs/latest/lib/python3.12/site-packages/felupe/constitution/hyperelasticity/_neo_hooke_compressible.py:162: RuntimeWarning: invalid value encountered in log
lnJ = np.log(J, out=J)
/home/docs/checkouts/readthedocs.org/user_builds/felupe/envs/latest/lib/python3.12/site-packages/felupe/constitution/hyperelasticity/_neo_hooke_compressible.py:195: RuntimeWarning: invalid value encountered in log
lnJ = np.log(J, out=J)
| 12,1 | 0+ => 3+ | 1.5e+03 ( 8#) |Failed |
| 13,1 | 0+ => 12- | 9.8e+01 ( 8#) |Failed |
| 14,1 | 0+ => 1+ | 1.1e+04 ( 8#) |Failed |
| 15,1 | 0+ => 2+ | 2.8e+02 ( 8#) |Failed |
| 16,1 | 0+ => 12+ | 3.0e+05 ( 8#) |Failed |
| 17,1 | 0+ => 1+ | 4.0e+02 ( 8#) |Failed |
| 18,1 | 0+ => 0+ | 4.4e-06 ( 6#) | |
| 19,1 | 0+ => 0+ | 2.0e-04 ( 6#) | |
| 20,1 | 0+ => 0+ | 2.4e-05 ( 7#) | |
| 21,1 | 0+ => 0+ | 3.4e-04 ( 7#) | |
| 22,1 | 0+ => 0+ | 2.3e-04 ( 7#) | |
| 23,1 | 0+ => 0+ | 5.0e-05 ( 6#) | |
| 24,1 | 0+ => 0+ | 2.3e-05 ( 6#) | |
| 25,1 | 0+ => 0+ | 2.8e-07 ( 6#) | |
| 26,1 | 0+ => 0+ | 2.6e-05 ( 6#) | |
| 27,1 | 0+ => 0+ | 7.3e-06 ( 7#) | |
| 28,1 | 0+ => 0+ | 1.1e-05 ( 7#) | |
| 29,1 | 0+ => 12+ | 1.3e+03 ( 8#) |Failed |
| 30,1 | 0+ => 0+ | 1.6e-04 ( 8#) | |
| 31,1 | 0+ => 12+ | 3.7e+01 ( 8#) |Failed |
| 32,1 | 0+ => 12+ | 7.7e-04 ( 6#) | => re-Cycle |
| 2 | 12+ => 12+ | 3.8e-08 ( 5#) | |
| 33,1 | 12+ => 12+ | 8.4e-05 ( 6#) | |
| 34,1 | 12+ => 12+ | 5.8e-04 ( 6#) | |
| 35,1 | 12+ => 12+ | 3.8e-04 ( 6#) | |
| 36,1 | 12+ => 12+ | 3.0e-03 ( 8#) |Failed |
| 37,1 | 12+ => 12+ | 1.1e-05 ( 7#) | |
| 38,1 | 12+ => 12+ | 3.3e-04 ( 7#) | |
| 39,1 | 12+ => 12+ | 7.4e-09 ( 6#) | |
| 40,1 | 12+ => 12+ | 1.7e-04 ( 6#) | |
| 41,1 | 12+ => 12+ | 6.1e-04 ( 7#) | |
| 42,1 | 12+ => 0+ | 4.3e-07 ( 6#) | => re-Cycle |
| 2 | 0+ => 0+ | 9.6e-09 ( 5#) | |
| 43,1 | 0+ => 0+ | 9.1e-05 ( 6#) | |
| 44,1 | 0+ => 0+ | 6.6e-07 ( 7#) | |
| 45,1 | 0+ => 1+ | 1.7e+03 ( 8#) |Failed |
| 46,1 | 0+ => 0+ | 1.6e-04 ( 7#) | |
| 47,1 | 0+ => 0+ | 1.2e-05 ( 7#) | |
| 48,1 | 0+ => 0+ | 6.0e-05 ( 8#) | |
| 49,1 | 0+ => 0+ | 3.2e+02 ( 8#) |Failed |
| 50,1 | 0+ => 0+ | 3.0e-04 ( 6#) | |
| 51,1 | 0+ => 0+ | 1.2e-04 ( 6#) | |
| 52,1 | 0+ => 0+ | 4.8e-05 ( 6#) | |
| 53,1 | 0+ => 0+ | 3.0e-04 ( 6#) | |
| 54,1 | 0+ => 12+ | 1.8e+04 ( 8#) |Failed |
| 55,1 | 0+ => 0+ | 7.2e+00 ( 8#) |Failed |
| 56,1 | 0+ => 0+ | 1.8e-04 ( 7#) | |
| 57,1 | 0+ => 0+ | 2.1e+01 ( 8#) |Failed |
| 58,1 | 0+ => 0+ | 7.1e-08 ( 6#) | |
| 59,1 | 0+ => 0+ | 3.5e-08 ( 6#) | |
| 60,1 | 0+ => 0+ | 3.8e-08 ( 6#) | |
| 61,1 | 0+ => 0+ | 8.2e-07 ( 6#) | |
| 62,1 | 0+ => 0+ | 6.9e-06 ( 6#) | |
| 63,1 | 0+ => 0+ | 2.7e-01 ( 8#) |Failed |
| 64,1 | 0+ => 0+ | 2.1e-07 ( 6#) | |
| 65,1 | 0+ => 0+ | 1.2e-07 ( 6#) | |
The unstable pressure-controlled equilibrium path is plotted as pressure-displacement curve.
import matplotlib.pyplot as plt
plt.plot(X[:, 0], X[:, -1], "x-", lw=3)
plt.xlabel(r"Max. Displacement $u_1(X_1=X_2=0)$ $\longrightarrow$")
plt.ylabel(r"Load-Proportionality-Factor $\lambda$ $\longrightarrow$")
The 3d-deformed configuration of the solid body is plotted.
mesh_3d = mesh.revolve(phi=90, n=6)
region_3d = fem.RegionHexahedron(mesh_3d)
values = mesh.copy(points=field[0].values).revolve(phi=90, n=6).points
u_3d = fem.Field(region_3d, values=values, dim=3)
field_3d = fem.FieldContainer([u_3d])
solid_3d = fem.SolidBodyNearlyIncompressible(umat, field_3d, bulk=5000)
solid_3d.plot("Principal Values of Cauchy Stress", project=fem.topoints).show()
Total running time of the script: (0 minutes 5.435 seconds)