Third Medium Contact with quadratic quads #923

adtzlr · 2025-01-02T22:20:36Z

adtzlr
Jan 2, 2025
Maintainer

This example shows how to use the third medium contact method with quadratic quads and automatic differentiation for the strain energy function. This may be used as a starting point for topolgy optimization including contact.

code

r"""
Third Medium Contact
--------------------

.. topic:: Frictionless contact method simulated by a third medium [1]_.

   * create a mesh container with multiple meshes
   
   * define multiple solid bodies and create a top-level field
   
   * create an animation of the deformed solid

This contact method uses a third medium for two solid contact bodies with a a Hessian-
based regularization [2]_. First, let's create sub meshes with quadratic quad cells for
the solid body. All sub meshes are merged by stacking the meshes of the
:class:`mesh container <felupe.MeshContainer>` into a :class:`mesh <felupe.Mesh>`.
"""

import numpy as np
import pypardiso
import felupe as fem
import felupe.constitution.tensortrax as mat
import tensortrax.math as tm

t = 0.1
L = 1
H = 0.5

nt, nL, nH = (5, 36, 13)

a = fem.Rectangle(a=(0, 0), b=(t, t), n=nt)
b = a.translate(move=H - t, axis=1)
c = fem.Rectangle(a=(t, 0), b=(L, t), n=(nL, nt))
d = c.translate(move=H - t, axis=1)
e = fem.Rectangle(a=(0, t), b=(t, H - t), n=(nt, nH))
body = fem.MeshContainer([a, b, c, d, e], merge=True).stack().add_midpoints_edges()

# %%
# This is repeated for the medium, where two meshes are added to a mesh container. Extra
# cells are added on the right edge to improve convergence.
f = fem.Rectangle(a=(t, t), b=(L, H - t), n=(nL, nH))
g = fem.Rectangle(a=(L, 0), b=(L + t / (nt - 1), H), n=(2, 2 * (nt - 1) + nH))
medium = fem.MeshContainer([f, g], merge=True).stack().add_midpoints_edges()

# %%
# Both stacked meshes are added to a top-level mesh container. For each mesh of the
# container, quad regions and plane-strain vector-fields for the displacements are
# created. The Neo-Hookean isotropic hyperelastic material formulation from the solid
# body is also used for the isotropic part of the background material but with scaled
# down material parameters.
G = 5 / 14
K = 5 / 3
kr = 1e-6
gamma0 = 1e-6


def neo_hooke(C, mu, bulk):
    I3 = tm.linalg.det(C)
    I1 = tm.trace(C)
    return bulk * tm.log(I3) ** 2 / 8 + mu * (I3 ** (-1 / 3) * I1 - 3) / 2


container = fem.MeshContainer([body, medium], merge=True, decimals=6)
regions = [fem.RegionQuadraticQuad(m, hess=True) for m in container.meshes]
fields = [fem.FieldPlaneStrain(r, dim=2).as_container() for r in regions]
umats = [
    mat.Hyperelastic(neo_hooke, mu=G, bulk=K),
    mat.Hyperelastic(neo_hooke, mu=G * gamma0, bulk=K * gamma0),
]
solids = [fem.SolidBody(umat, f) for umat, f in zip(umats, fields)]

# %%
# A top-level plane-strain field is created on a stacked mesh. This field includes all
# unknowns and is required for the selection of prescribed degrees of freedom as well
# as for Newton's method. The left end edge is fixed and the vertical displacement is
# prescribed on the outermost top-right point.
region = fem.RegionQuadraticQuad(container.stack())
field = fem.FieldPlaneStrain(region, dim=2).as_container()

bounds = dict(
    fixed=fem.Boundary(field[0], fx=0),
    move=fem.Boundary(field[0], fx=L, fy=H, mode="and", skip=(1, 0)),
)

# %%
# The so-called HuHu-regularization is created by two weak-:func:`forms <felupe.Form>`,
# which are derived from the regularization potential, see Eq.
# :eq:`huhu-regularization` [3]_.
#
# .. math::
#    :label: huhu-regularization
#
#    \Pi = \frac{1}{2} \int_V
#        \mathbb{H}(\boldsymbol{u})~\vdots~\mathbb{H}(\boldsymbol{u})~dV
#
from felupe.math import dddot, hess


@fem.Form(v=fields[1], u=fields[1], kwargs={"kr": kr * K * L**2})
def bilinearform():
    return [lambda v, u, kr: kr * dddot(hess(u), hess(v))]


@fem.Form(v=fields[1], kwargs={"kr": kr * K * L**2})
def linearform():
    u = fields[1][0]
    return [lambda v, kr: kr * dddot(hess(v), hess(u)[:2, :2, :2])]


regularization = fem.FormItem(
    bilinearform=bilinearform,
    linearform=linearform,
    kwargs={"kr": kr * K * L**2},
)

# %%
# The prescribed displacement is ramped up to the maximum value and released until zero.
move = fem.math.linsteps([0, 1, 0], num=50)
step = fem.Step(
    items=[*solids, regularization],
    ramp={bounds["move"]: move * -0.4 * L},
    boundaries=bounds,
)

# %%
# The deformation is captured in a GIF-file. The top-level field has to be passed as the
# ``x0``-argument for :func:`Newton's method <felupe.newtonrhapson>`, which is called on
# evaluation. After all frames are recorded, it is important to ``close()`` the plotter.
plotter = container.plot(
    colors=["grey", "white"],
    line_width=3,
    off_screen=True,
    show_edges=False,
)
plotter.open_gif("third-medium-contact.gif", fps=25)


def record(stepnumber, substepnumber, substep, plotter):
    "Update the mesh of the plotter and write a frame."
    if substepnumber in np.arange(len(move), step=1):
        for m in plotter.meshes:
            m.points[:, :2] = region.mesh.points + field[0].values
        plotter.write_frame()


job = fem.Job([step], callback=record, plotter=plotter)
job.evaluate(x0=field, solver=pypardiso.spsolve, maxiter=25)

plotter.close()

# %%
# References
# ~~~~~~~~~~
#
# .. [1] https://en.wikipedia.org/wiki/Third_medium_contact_method
#
# .. [2] G. L. Bluhm, O. Sigmund, and K. Poulios, "Internal contact modeling for finite
#        strain topology optimization", Computational Mechanics, vol. 67, no. 4.
#        Springer Science and Business Media LLC, pp. 1099–1114, Mar. 04, 2021. |DOI-2|.
#
# .. [3] A. H. Frederiksen, O. Sigmund, and K. Poulios, "Topology optimization of self-
#        contacting structures", Computational Mechanics, vol. 73, no. 4. Springer
#        Science and Business Media LLC, pp. 967–981, Oct. 07, 2023. |DOI-3|.
#
# .. |DOI-2| image:: https://zenodo.org/badge/DOI/10.1007/s00466-021-01974-x.svg
#    :target: https://www.doi.org/10.1007/s00466-021-01974-x
#
# .. |DOI-3| image:: https://zenodo.org/badge/DOI/10.1007/s00466-023-02396-7.svg
#    :target: https://www.doi.org/10.1007/s00466-023-02396-7
#

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Third Medium Contact with quadratic quads #923

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{editor}}'s edit

{{editor}}'s edit

Uh oh!

Replies: 0 comments

Select a reply

Uh oh!

Third Medium Contact with quadratic quads #923

Uh oh!

Uh oh!

adtzlr Jan 2, 2025 Maintainer

Replies: 0 comments

adtzlr
Jan 2, 2025
Maintainer