Nitsche Method : Fully tied 3D Contact

Introduction

In this article, we take advantage of this capability to solve an optimization problem that fully constrains three-dimensional displacements using the mortar method.

This approach offers a key advantage — it enables the treatment of various combinations of bodies without altering their mesh structures.

Algorithm

Nitche Method

Instead of enforcing the contact constraint through a Lagrange multiplier (which introduces extra unknowns), Nitsche’s method weakly imposes the constraint by modifying the weak form with consistent terms.

The key idea is to:

• – Add penalty-like terms to stabilize the interface,
• – Add consistent adjoint terms to maintain symmetry and consistency.

This results in an augmented bilinear form on the contact surface.

In the general form of force equilibrium equation, we will use below weak form

\int_{\Omega} \boldsymbol{\varepsilon}(\boldsymbol{v}) : \boldsymbol{\sigma}(\boldsymbol{u}) , d\Omega = \int_{\Omega} \boldsymbol{v} \cdot \boldsymbol{f} , d\Omega
+
\int_{\Gamma_t} \boldsymbol{v} \cdot \bar{\boldsymbol{t}} , d\Gamma

And we would obtain linear equation from the weak form

\begin{align}
\mathbf{K}u =F_{N}+F_{body}\\
\begin{bmatrix}
\mathbf{K}_1 & \mathbf{0} \\
\mathbf{0} & \mathbf{K}_2
\end{bmatrix}
\begin{bmatrix}
\mathbf{u}_1 \\[4pt]
\mathbf{u}_2
\end{bmatrix}
=
\begin{bmatrix}
\mathbf{f}_1 \\[4pt]
\mathbf{f}_2
\end{bmatrix}
\end{align}

In the Nitsche method, a penalty-like matrix is added to the stiffness matrix obtained from the above weak formulation by using bilinear form on the contact surface as it is mentioned above,

\begin{aligned}
a_\Gamma(u_1, u_2; v_1, v_2)
&=\\
\int_{\Gamma}
\Big[
\frac{\alpha}{h} (u_1 - u_2) \cdot (v_1 - v_2)\\
- \frac{1}{2} \big( (\boldsymbol{t}_1(u_1) + \boldsymbol{t}_2(u_2)) \cdot (v_1 - v_2) \big)\\
- \frac{1}{2} \big( (\boldsymbol{t}_1(v_1) + \boldsymbol{t}_2(v_2)) \cdot (u_1 - u_2) \big)
  \Big] ds
  \end{aligned}
 

\boldsymbol{t}_i(u_i) = \boldsymbol{C} : \boldsymbol{\varepsilon}(u_i)  \boldsymbol{n},
\qquad
\boldsymbol{\varepsilon}(u_i) = \tfrac{1}{2}(\nabla u_i + \nabla u_i^\mathsf{T})\\

where

\boldsymbol{C} is the elasticity tensor,
\boldsymbol{n} is the unit normal vector on the interface ( \Gamma ),
h is the characteristic element size,
and \alpha is the penalty parameter.

From this form, we would get K_{Nitsche}.

(K+K_{Nitsche})u=F_{N}+F_{body}

The meaning of 1st term is that it prevents jump of displacement on the interface. Additionally, the meaning of 2nd and 3rd terms are kind of energy (like compliance).

Implementation with scikit-fem

Code

import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import minres, spsolve

import skfem
from skfem import FacetBasis
from skfem.helpers import dot, mul, sym_grad
from skfem.supermeshing import elementwise_quadrature, intersect

from common import (
    assemble_top_pressure_load,
    make_default_meshes,
    make_elasticity_form,
    make_vector_p1_bases,
)

# ----------------------------
# 0) Meshes / bases / material
# ----------------------------
mesh_top, mesh_bot = make_default_meshes()
elem, basis_top, basis_bot = make_vector_p1_bases(mesh_top, mesh_bot)
a_elastic, E, nu, mu, lam = make_elasticity_form(E=210e9, nu=0.30)

K1 = a_elastic.assemble(basis_top).tocsr()
K2 = a_elastic.assemble(basis_bot).tocsr()

# ----------------------------
# 1) Supermesh on contact facets (non-matching)
# ----------------------------
m1t, orig1 = mesh_top.trace("contact", mtype=skfem.MeshTri, project=lambda p: p[[0, 1]])
m2t, orig2 = mesh_bot.trace("contact", mtype=skfem.MeshTri, project=lambda p: p[[0, 1]])

m12, t1, t2 = intersect(m1t, m2t)
quad1 = elementwise_quadrature(m1t, m12, t1)
quad2 = elementwise_quadrature(m2t, m12, t2)

fb_top = FacetBasis(mesh_top, elem, facets=orig1[t1], quadrature=quad1)
fb_bot = FacetBasis(mesh_bot, elem, facets=orig2[t2], quadrature=quad2)

# ----------------------------
# 2) Symmetric Nitsche for tied contact: u_top = u_bot
# ----------------------------
gamma = 20.0
alpha0 = lam + 2.0 * mu


@skfem.BilinearForm
def a_nitsche_tied(u1, u2, v1, v2, w):
    n = w.n

    ju = u1 - u2
    jv = v1 - v2

    eps_u1 = sym_grad(u1)
    eps_u2 = sym_grad(u2)
    eps_v1 = sym_grad(v1)
    eps_v2 = sym_grad(v2)

    t1 = 2.0 * mu * mul(eps_u1, n) + lam * np.trace(eps_u1) * n
    t2 = 2.0 * mu * mul(eps_u2, n) + lam * np.trace(eps_u2) * n
    s1 = 2.0 * mu * mul(eps_v1, n) + lam * np.trace(eps_v1) * n
    s2 = 2.0 * mu * mul(eps_v2, n) + lam * np.trace(eps_v2) * n

    penalty = gamma * alpha0 / w.h

    return (
        penalty * dot(ju, jv)
        - 0.5 * dot(t1 + t2, jv)
        - 0.5 * dot(s1 + s2, ju)
    )


B = skfem.asm(a_nitsche_tied, fb_top * fb_bot).tocsr()
K = sp.bmat([[K1, None], [None, K2]], format="csr") + B

# ----------------------------
# 3) External load + Dirichlet
# ----------------------------
F_top, area, pressure = assemble_top_pressure_load(mesh_top, elem, total_force=100.0)

F = np.zeros(K.shape[0])
F[:basis_top.N] = F_top

D_bot = basis_bot.get_dofs("dirichlet").all()
D = D_bot + basis_top.N

Kc, Fc, x0, I = skfem.condense(K, F, D=D)

try:
    u_red = spsolve(Kc, Fc)
except Exception:
    u_red, info = minres(Kc, Fc, rtol=1e-10, maxiter=20000)
    if info != 0:
        raise RuntimeError(f"MINRES did not converge, info={info}")

u = np.zeros(K.shape[0])
u[I] = u_red

res = K @ u - F
print("K:", K.shape, "nnz:", K.nnz)
print("||res||2 =", np.linalg.norm(res))
print("relative ||res||2 / ||F||2 =", np.linalg.norm(res) / (np.linalg.norm(F) + 1e-30))


# ----------------------------
# 4) Export combined VTU
# ----------------------------
def export_combined_vtu(mesh_top, mesh_bot, u, basis_top, file_path="combined.vtu"):
    import meshio

    n_top = mesh_top.p.shape[1]
    n_bot = mesh_bot.p.shape[1]
    points = np.vstack([mesh_top.p.T, mesh_bot.p.T])
    cells = [("tetra", np.vstack([mesh_top.t.T, mesh_bot.t.T + n_top]))]

    u_top = u[:basis_top.N].reshape(n_top, 3)
    u_bot = u[basis_top.N:].reshape(n_bot, 3)
    u_all = np.vstack([u_top, u_bot])

    node_tag = np.zeros(n_top + n_bot, dtype=np.int32)
    top_nodes_contact = np.unique(mesh_top.facets[:, mesh_top.boundaries["contact"]])
    bot_nodes_contact = np.unique(mesh_bot.facets[:, mesh_bot.boundaries["contact"]]) + n_top
    node_tag[np.hstack([top_nodes_contact, bot_nodes_contact])] = 1

    bot_nodes_dir = np.unique(mesh_bot.facets[:, mesh_bot.boundaries["dirichlet"]]) + n_top
    node_tag[bot_nodes_dir] = 2

    top_nodes_force = np.unique(mesh_top.facets[:, mesh_top.boundaries["force"]])
    node_tag[top_nodes_force] = 3

    meshio.Mesh(
        points=points,
        cells=cells,
        point_data={"u": u_all, "node_tag": node_tag},
    ).write(file_path)

    print(f"wrote {file_path}")


export_combined_vtu(mesh_top, mesh_bot, u, basis_top, "tied_contact-nitsche.vtu")