Solid Bodies as Items for Assembly

The mechanics submodule contains classes for the generation of solid bodies. Solid body objects are supported as items of a Step and the newtonrhapson() procedure.

Solid Body

The generation of internal force vectors and stiffness matrices of solid bodies are provided as assembly-methods of a SolidBody or a SolidBodyNearlyIncompressible.

import felupe as fem

mesh = fem.Cube(n=6)
region = fem.RegionHexahedron(mesh)
field = fem.FieldContainer([fem.Field(region, dim=3)])

# a solid body
body = fem.SolidBody(umat=fem.NeoHooke(mu=1, bulk=2), field=field)

# a (nearly) incompressible solid body (to be used with quads and hexahedrons)
body = fem.SolidBodyNearlyIncompressible(umat=fem.NeoHooke(mu=1), field=field, bulk=5000)

internal_force = body.assemble.vector(field, parallel=False, jit=False)
stiffness_matrix = body.assemble.matrix(field, parallel=False, jit=False)

During assembly, several results are stored, e.g. the gradient of the strain energy density function per unit undeformed volume w.r.t. the deformation gradient (first Piola-Kirchhoff stress tensor). Other results are the deformation gradient or the fourth-order elasticity tensor associated to the first Piola-Kirchhoff stress tensor.

\[ \begin{align}\begin{aligned}\boldsymbol{F} &= \frac{\partial \boldsymbol{x}}{\partial \boldsymbol{X}}\\\boldsymbol{P} &= \frac{\partial \psi(\boldsymbol{C}(\boldsymbol{F}))}{\partial \boldsymbol{F}}\\\mathbb{A} &= \frac{\partial^2 \psi(\boldsymbol{C}(\boldsymbol{F}))}{\partial \boldsymbol{F}\ \partial \boldsymbol{F}}\end{aligned}\end{align} \]

F = body.results.kinematics
P = body.results.stress
A = body.results.elasticity

The Cauchy stress tensor, as well as the gradient and the hessian of the strain energy density function per unit undeformed volume are obtained by evaluate-methods of the solid body.

\[ \begin{align}\begin{aligned}\boldsymbol{P} &= \frac{\partial \psi(\boldsymbol{C}(\boldsymbol{F}))}{\partial \boldsymbol{F}}\\\mathbb{A} &= \frac{\partial^2 \psi(\boldsymbol{C}(\boldsymbol{F}))}{\partial \boldsymbol{F}\ \partial \boldsymbol{F}}\\\boldsymbol{\sigma} &= \frac{1}{J} \boldsymbol{P} \boldsymbol{F}^T\end{aligned}\end{align} \]

P = body.evaluate.gradient(field)
A = body.evaluate.hessian(field)
s = body.evaluate.cauchy_stress(field)

Body Force (Gravity) on a Solid Body

The generation of external force vectors or stiffness matrices of body forces acting on solid bodies are provided as assembly-methods of a SolidBodyGravity.

\[\delta W_{ext} = \int_V \delta \boldsymbol{u} \cdot \rho \boldsymbol{g} \ dV\]

body = fem.SolidBodyGravity(field=field, gravity=[9810, 0, 0], density=7.85e-9)

force_gravity = body.assemble.vector(field, parallel=False, jit=False)

Pressure Boundary on a Solid Body

The generation of force vectors or stiffness matrices of pressure boundaries on solid bodies are provided as assembly-methods of a SolidBodyPressure.

\[\delta W_{ext} = \int_{\partial V} \delta \boldsymbol{u} \cdot p \ J \boldsymbol{F}^{-T} \ d\boldsymbol{A}\]

region_pressure = fem.RegionHexahedronBoundary(
    mesh=mesh,
    only_surface=True, # select only faces on the outline
    mask=mesh.points[:, 0] == 0, # select a subset of faces on the surface
)

field_boundary = fem.FieldContainer([fem.Field(region_pressure, dim=3)])
field_boundary.link(field)

body_pressure = fem.SolidBodyPressure(field=field_boundary)

force_pressure = body_pressure.assemble.vector(
    field=field_boundary, parallel=False, jit=False
)

stiffness_matrix_pressure = body_pressure.assemble.matrix(
    field=field_boundary, parallel=False, jit=False
)

For axisymmetric problems the boundary region has to be created with the attribute ensure_3d=True.

mesh = fem.Rectangle(a=(0, 30), b=(20, 40), n=(21, 11))
region = fem.RegionQuad(mesh)

region_pressure = fem.RegionQuadBoundary(
    mesh=mesh,
    only_surface=True, # select only faces on the outline
    mask=mesh.points[:, 0] == 0, # select a subset of faces on the surface
    ensure_3d=True, # flag for axisymmetric boundary region
)

field = fem.FieldContainer([fem.FieldAxisymmetric(region)])
field_boundary = fem.FieldContainer([fem.FieldAxisymmetric(region_pressure)])
field_boundary.link(field)