Distributed force on hyperelastic object

[KaiWolven]

Hi,
I've been trying for some time to apply deformations to hyperelastic objects, which is quite easy using Felupe. However, know that I plan on more complex stuff, I wanted to apply forces directly to the top surface of an object, specifically at Z.max(). I noticed the existence of the fem.SolidBodyForce function, but I'm not sure if it can be used in this case or if there are any other possible options.
Thank you.

[adtzlr]

Hi,

SolidBodyForce is for things like gravity (body force). Do you want to apply a concentrated force on a surface? You could have a look on Ex. 8 which uses a MultiPointConstraint. Consider the following examples:

Note

It's important to remove the MPC - center point from mesh.points_without_cells. In this array, all hanging points are collected and ignored to avoid errors. But we would like to include the degrees-of-freedom of this center point in this example.

import felupe as fem
import numpy as np

mesh = fem.Cube(n=6)
mesh.update(points=np.vstack([mesh.points, [2, 0.5, 0.5]]))
mesh.points_without_cells = np.array([], dtype=int)

region = fem.RegionHexahedron(mesh)
field = fem.FieldContainer([fem.Field(region, dim=3)])
solid = fem.SolidBody(umat=fem.NeoHooke(mu=1, bulk=2), field=field)

fixed = dict(fixed=fem.Boundary(field[0], fx=0))
points = np.arange(mesh.npoints)[mesh.x == 1.0]
mpc = fem.MultiPointConstraint(field, points=points, centerpoint=-1)
force = fem.PointLoad(field, points=[-1])

table = fem.math.linsteps([0, 1], num=5, axis=0, axes=3)
ramp = {force: 1 * table}
step = fem.Step(items=[solid, mpc, force], ramp=ramp, boundaries=fixed)

job = fem.Job(steps=[step]).evaluate()
plotter = solid.plot("Principal Values of Cauchy Stress", plotter=mpc.plot())
plotter.show()

It is also possible to apply the MPC only in the direction of the force, i.e.

import felupe as fem
import numpy as np

mesh = fem.Cube(n=6)
mesh.update(points=np.vstack([mesh.points, [2, 0.5, 0.5]]))
mesh.points_without_cells = np.array([], dtype=int)

region = fem.RegionHexahedron(mesh)
field = fem.FieldContainer([fem.Field(region, dim=3)])
solid = fem.SolidBody(umat=fem.NeoHooke(mu=1, bulk=2), field=field)

boundaries = dict(
    left=fem.Boundary(field[0], fx=0),
    right=fem.Boundary(field[0], fx=2, skip=(1, 0, 0)),
)
points = np.arange(mesh.npoints)[mesh.x == 1.0]
mpc = fem.MultiPointConstraint(field, points=points, centerpoint=-1, skip=(0, 1, 1))
force = fem.PointLoad(field, points=[-1])

table = fem.math.linsteps([0, 1], num=5, axis=0, axes=3)
ramp = {force: table}
step = fem.Step(items=[solid, mpc, force], ramp=ramp, boundaries=boundaries)

job = fem.Job(steps=[step]).evaluate()
plotter = solid.plot("Principal Values of Cauchy Stress", plotter=mpc.plot())
plotter.show()

If you know your force distribution, you can of course also apply the force on the mesh-points directly.

import felupe as fem
import numpy as np

mesh = fem.Cube(n=6)
region = fem.RegionHexahedron(mesh)
field = fem.FieldContainer([fem.Field(region, dim=3)])
solid = fem.SolidBody(umat=fem.NeoHooke(mu=1, bulk=2), field=field)

boundaries = dict(left=fem.Boundary(field[0], fx=0))
points = np.arange(mesh.npoints)[mesh.x == mesh.x.max()]
force = fem.PointLoad(field, points=points)

table = fem.math.linsteps([0, 1], num=5, axis=0, axes=3)
ramp = {force: 1 / len(force.points) * table}
step = fem.Step(items=[solid, force], ramp=ramp, boundaries=boundaries)

job = fem.Job(steps=[step]).evaluate()
solid.plot("Principal Values of Cauchy Stress").show()

Hope that helps!

[KaiWolven]

The main problem with tetra meshes is that I'm not being able to read externally generated ones (I'm using Gmsh to transform already existing stl files), modify the number of elements and operate with it inside Felupe.

[adtzlr]

Okay, I understand. There is

import felupe

container = fem.mesh.read("my_mesh.msh")
mesh = container.meshes[0]

Does that function work for you?

[KaiWolven]

If i start with a closed tetra-mesh like this:

After applying that function, it seems to read just fine, both in .gms and .vtk:

However, it won't continue working after the fem.RegionTetraMINI(mesh) bit. Printing the mesh class after reading gives "class 'felupe.mesh._mesh.Mesh'", which seems fine. I've tried many approaches when generating the tetra-mesh, but still cannot get it to work properly. I've also tried other Region types in case TetraMINI wasn't compatible. Moreover, if i reduce the opacity when plotting, there seems to be no inner meshing; I donnot know if its just the visualization of there and actual problem when reading the file:

Hope this explains the problem a little bit better.

[adtzlr]

It seems that PyVista doesn't show points inside the volume if you set the opacity to a value lower than 1. I think you'd have to modify the mesh before to do so.

Here is an example. Let's take a gmsh-file, e.g. cube6.geo. First, create the msh-file within gmsh and then follow this code to use it in FElupe.

code (requires the pypardiso-package)

import felupe as fem
import pypardiso

mesh = fem.mesh.read("cube6.msh")[-1]
region = fem.RegionTetraMINI(mesh.add_midpoints_volumes(cell_type="tetra"))
field = fem.FieldsMixed(region, n=3)

boundaries, loadcase = fem.dof.uniaxial(field, clamped=True)

umat = fem.NearlyIncompressible(fem.NeoHooke(mu=1), bulk=5000)
solid = fem.SolidBody(umat, field)

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

job = fem.CharacteristicCurve(steps=[step], boundary=boundaries["move"])
job.evaluate(solver=pypardiso.spsolve)

fig, ax = job.plot(
    xlabel="Displacement $d_1$ in mm $\longrightarrow$",
    ylabel="Normal Force $F_1$ in N $\longrightarrow$",
)
solid.plot("Principal Values of Cauchy Stress").show()

[KaiWolven]

Thanks for the info! It seems to work just fine with a tetra mesh. The problem was the definition I had for the region part.

I just have one last question: in that example, how would you apply other hyperelastic models appart from NeoHooke, such as Mooney-Rivlin or Ogden? I've tried with fem.Hyperelastic and other previous felupe definitions, but it seems it is not compatible with this workflow.

[adtzlr]

That's great @KaiWolven! 🎉 🎉

To use other hyperelastic materials from FElupe, use

umat = fem.NearlyIncompressible(
    fem.Hyperelastic(fem.mooney_rivlin, C10=0.4, C01=0.1), bulk=5000
)
solid = fem.SolidBody(umat, field)

or with a custom strain energy density function:

import tensortrax.math as tm

def my_mooney_rivlin(C, C10, C01):
    I1 = tm.trace(C)
    I2 = (I1**2 - tm.trace(C @ C)) / 2
    I3 = tm.linalg.det(C)
    return C10 * (I3**(-1/3) * I1 - 3) + C01 * (I3**(-2/3) * I2 - 3)
    
umat = fem.NearlyIncompressible(
    fem.Hyperelastic(my_mooney_rivlin, C10=0.4, C01=0.1), bulk=5000
)
solid = fem.SolidBody(umat, field)

You have to wrap your model in felupe.NearlyIncompressible(..., bulk=5000) to add a volumetric strain energy density function and to create a mixed-field formulation (displacement, pressure, volume-ratio). SolidBodyNearlyIncompressible(..., bulk=5000) does the same, but in a more efficient way for quads and hexahedrons.