Weak Form External Displacement Boundary Condition

[adtzlr]

This is an idea to apply an external displacement by a weak form item.

import felupe as fem
from felupe.math import ddot
import numpy as np

mesh = fem.Cube(n=6)
region = fem.RegionHexahedron(mesh)
field = fem.FieldContainer([fem.Field(region, dim=3)])
boundaries = fem.dof.symmetry(field[0])
umat = fem.NeoHooke(mu=1)
solid = fem.SolidBodyNearlyIncompressible(umat, field, bulk=5000)

surface = fem.FieldContainer(
    [fem.Field(fem.RegionHexahedronBoundary(mesh, mask=mesh.x == 1), dim=3)]
)


@fem.Form(v=surface)
def linearform():
    uext = np.array([-0.43, 0, 0]).reshape(3, 1, 1)

    def L(v, multiplier):
        u = surface.extract(grad=False)[0]
        return multiplier * ddot(v, u - uext)

    return [L]


@fem.Form(v=surface, u=surface)
def bilinearform():
    return [lambda v, u, multiplier: ddot(v, u) * multiplier]


move = fem.FormItem(bilinearform, linearform, kwargs=dict(multiplier=500))
step = fem.Step(items=[solid, move], boundaries=boundaries)
job = fem.Job(steps=[step])
job.evaluate()

solid.plot("Principal Values of Cauchy Stress").show()

[adtzlr]

For more details, see the Examples section in https://felupe.readthedocs.io/en/latest/felupe/assembly.html#felupe.FormItem

[adtzlr]

This approach enables a segment-to-surface contact. A simple example with a rigid sphere-shaped obstacle:

import felupe as fem
from felupe.math import dot
import numpy as np
import pypardiso

mesh = fem.Cube(n=11)
region = fem.RegionHexahedron(mesh)
field = fem.FieldContainer([fem.Field(region, dim=3)])
boundaries = fem.dof.symmetry(field[0])
umat = fem.NeoHooke(mu=1)
solid = fem.SolidBodyNearlyIncompressible(umat, field, bulk=5000)

Γ = fem.FieldContainer(
    [fem.Field(fem.RegionHexahedronBoundary(mesh, mask=mesh.x == 1), dim=3)]
)
X = fem.Field(
    fem.RegionHexahedronBoundary(mesh, mask=mesh.x == 1), values=mesh.points, dim=3
).interpolate()


def R(x, grad=False):
    "Return the (gradient of the) ramp function."
    if not grad:
        return (x + np.abs(x)) / 2
    else:
        return (np.ones_like(x) + np.sign(x)) / 2


def S(x, uext, r=np.sqrt(2), grad=False):
    if not grad:
        surface = uext + r - np.sqrt(r**2 - x[1] ** 2 - x[2] ** 2)
    else:
        surface = -(x[1] + x[2]) / np.sqrt(r**2 - x[1] ** 2 - x[2] ** 2)
    return surface


@fem.Form(v=Γ)
def linearform():

    def L(δu, uext=0, λ=1e3):
        u = Γ.extract(grad=False)[0]
        x = X + u
        w = u[0] - S(x, uext)
        return λ * dot(δu[0], R(w))

    return [L]


@fem.Form(v=Γ, u=Γ)
def bilinearform():

    def a(δu, Δu, uext=0, λ=1e2):
        u = Γ.extract(grad=False)[0]
        x = X + u
        w = u[0] - S(x, uext)
        return dot(δu[0], R(w, grad=True) * (Δu[0] - S(Δu, uext, grad=True))) * λ

    return [a]


obstacle = fem.FormItem(
    bilinearform, linearform, kwargs=dict(uext=0, λ=20), ramp_item=0
)
move = fem.math.linsteps([0, -0.4], 4)
step = fem.Step(items=[solid, obstacle], boundaries=boundaries, ramp={obstacle: move})
job = fem.Job(steps=[step])
job.evaluate(tol=1e-1, solver=pypardiso.spsolve, parallel=True)

solid.plot("Principal Values of Cauchy Stress").show()

Due to the fixed penalty multiplier, converge is not guaranted. A mixed formulation will solve this problem?

[adtzlr]

A nice addition would be an Obstacle class. This should contain the analytic surface-evaluations, its gradients as well as plot methods for viewing the model.

The above FormItem could be (re-casted in an IntegralForm representation and) provided as ObstacleContact?