Hyperelasticity

Find approximate solution to following non-linear system of PDEs

\[ \begin{align}\begin{aligned}\mathbf{u}_t &= \mathbf{v} &&\quad\text{ in }\Omega\times(0, T),\\\mathbf{v}_t &= \operatorname{div} (J \mathbb{T} \mathbb{F}^{-\top}) &&\quad\text{ in }\Omega\times(0, T),\\J^2 - 1 &= \begin{cases} 0 & \text{incompressible case} \newline -p/\lambda & \text{compressible case} \end{cases} &&\quad\text{ in }\Omega\times(0, T),\\\mathbf{u} = \mathbf{v} &= 0 &&\quad\text{ on }\Gamma_\mathrm{D}\times(0, T),\\J \mathbb{T} \mathbb{F}^{-\top} \mathbf{n} &= \mathbf{g} &&\quad\text{ on }\Gamma_\mathrm{N}\times(0, T),\\J \mathbb{T} \mathbb{F}^{-\top} \mathbf{n} &= 0 &&\quad\text{ on }\partial\Omega\backslash(\Gamma_\mathrm{D}\cup\Gamma_\mathrm{N})\times(0, T),\\\mathbf{u} = \mathbf{v} &= 0 &&\quad\text{ on }\Omega\times\{0\}\end{aligned}\end{align} \]

where

\[ \begin{align}\begin{aligned}\mathbb{F} &= \mathbb{I} + \nabla\mathbf{u},\\J &= \det{\mathbb{F}},\\\mathbb{B} &= \mathbb{F}\,\mathbb{F}^\top,\\\mathbb{T} &= -p\mathbb{I} + \mu (\mathbb{B-I})\end{aligned}\end{align} \]

using \(\theta\)-scheme discretization in time and arbitrary discretization in space with data

\[ \begin{align}\begin{aligned}\Omega &= \begin{cases} (0, 20) \times (0, 1) & \text{in 2D} \newline \text{lego brick } 10 \times 2 \times 1H & \text{in 3D} \end{cases}\\\Gamma_\mathrm{D} &= \begin{cases} \left\{ x=0 \right\} & \text{in 2D} \newline \left\{ x = \inf_{\mathbf{x}\in\Omega}{x} \right\} & \text{in 3D} \end{cases}\\\Gamma_\mathrm{N} &= \begin{cases} \left\{ x=20 \right\} & \text{in 2D} \newline \left\{ x = \sup_{\mathbf{x}\in\Omega}{x} \right\} & \text{in 3D} \end{cases}\\T &= 5,\\\mathbf{g} &= \begin{cases} J \mathbb{F}^{-\top} \Bigl[\negthinspace\begin{smallmatrix}0\newline100t\end{smallmatrix}\Bigr] & \text{in 2D} \newline J \mathbb{F}^{-\top} \Bigl[\negthinspace\begin{smallmatrix}0\newline0\newline100t\end{smallmatrix}\Bigr] & \text{in 3D} \end{cases}\\\mu &= \frac{E}{2(1+\nu)},\\\lambda &= \begin{cases} \infty & \text{incompressible case} \newline \frac{E\nu}{(1+\nu)(1-2\nu)} & \text{compressible case} \end{cases}\\E &= 10^5,\\\nu &= \begin{cases} 1/2 & \text{incompressible case} \newline 0.3 & \text{compressible case} \end{cases}\end{aligned}\end{align} \]

Mesh file of lego brick lego_beam.xml. Within shell download by

wget 

Task 1

Discretize the equation in time using the Crank-Nicolson scheme and derive a variational formulation of the problem. Consider discretization using P1/P1/P1 mixed element.

Task 2

Build 2D mesh:

mesh = RectangleMesh(Point(x0, y0), Point(x1, y1), 100, 5, 'crossed')

Prepare facet function marking \(\Gamma_\mathrm{N}\) and \(\Gamma_\mathrm{D}\) and plot it to check its correctness.

Hint

You can get coordinates of \(\Gamma_\mathrm{D}\) by something like x0 = mesh.coordinates()[:, 0].min() for lego mesh. Analogically for \(\Gamma_\mathrm{N}\).

Task 3

Define Cauchy stress and variational formulation of the problem.

Hint

Get geometric dimension by gdim = mesh.geometry().dim() to be able to write the code independently of the dimension.

Task 4

Prepare a solver and write simple time-stepping loop. Use time step \(\Delta t=\tfrac14\).

Prepare a solver by:

problem = NonlinearVariationalProblem(F, w, bcs=bcs, J=J)
solver = NonlinearVariationalSolver(problem)
solver.parameters['newton_solver']['relative_tolerance'] = 1e-6
solver.parameters['newton_solver']['linear_solver'] = 'mumps'

to increase the tolerance reasonably and employ powerful sparse direct solver MUMPS.

Prepare nice plotting of displacement by:

plot(u, mode="displacement")

Manipulate the plot how shown in the Matplotlib note.

Task 4

Solve the compressible 2D problem.

Solve the incompressible 2D problem.

Task 5

Solve the 3D compressible problem. Use time step \(\Delta t=\tfrac12\).

Load mesh by:

mesh = Mesh('lego_beam.xml')

Use the following optimization:

# Limit quadrature degree
dx = dx(degree=4)
ds = ds(degree=4)

You can also try to run the 3D problem in parallel:

# Disable plotting
export MPLBACKEND=template
export DOLFIN_NOPLOT=1

# Run the code on <np> processors
mpirun -n <np> python <script>.py

Task 6

Plot computed displacement \(u\) in Paraview using Warp by vector filter.

Reference solution

Show/Hide Code

Download Code

from dolfin import *
import matplotlib.pyplot as plt
import os


def solve_elasticity(facet_function, E, nu, dt, T_end, output_dir):
    """Solves elasticity problem with Young modulus E, Poisson ration nu,
    timestep dt, until T_end and with output data going to output_dir.
    Geometry is defined by facet_function which also defines rest boundary
    by marker 1 and traction boundary by marker 2."""

    # Get mesh and prepare boundary measure
    mesh = facet_function.mesh()
    gdim = mesh.geometry().dim()
    dx = Measure("dx")
    ds = Measure("ds", subdomain_data=facet_function, subdomain_id=2)

    # Limit quadrature degree
    dx = dx(degree=4)
    ds = ds(degree=4)

    # Build function space
    element_v = VectorElement("P", mesh.ufl_cell(), 1)
    element_s = FiniteElement("P", mesh.ufl_cell(), 1)
    mixed_element = MixedElement([element_v, element_v, element_s])
    W = FunctionSpace(mesh, mixed_element)
    info("Num DOFs {}".format(W.dim()))

    # Prepare BCs
    bc0 = DirichletBC(W.sub(0), gdim*(0,), facet_function, 1)
    bc1 = DirichletBC(W.sub(1), gdim*(0,), facet_function, 1)
    bcs = [bc0, bc1]

    # Define constitutive law
    def stress(u, p):
        """Returns 1st Piola-Kirchhoff stress and (local) mass balance
        for given u, p."""
        mu = Constant(E/(2.0*(1.0 + nu)))
        F = I + grad(u)
        J = det(F)
        B = F * F.T
        T = -p*I + mu*(B-I) # Cauchy stress
        S = J*T*inv(F).T # 1st Piola-Kirchhoff stress
        if nu == 0.5:
            # Incompressible
            pp = J-1.0
        else:
            # Compressible
            lmbd = Constant(E*nu/((1.0 + nu)*(1.0 - 2.0*nu)))
            pp = 1.0/lmbd*p + (J*J-1.0)
        return S, pp

    # Timestepping theta-method parameters
    q = Constant(0.5)
    dt = Constant(dt)

    # Unknowns, values at previous step and test functions
    w = Function(W)
    u, v, p = split(w)
    w0 = Function(W)
    u0, v0, p0 = split(w0)
    _u, _v, _p = TestFunctions(W)

    I = Identity(W.mesh().geometry().dim())

    # Balance of momentum
    S, pp = stress(u, p)
    S0, pp0 = stress(u0, p0)
    F1 = (1.0/dt)*inner(u-u0, _u)*dx \
       - ( q*inner(v, _u)*dx + (1.0-q)*inner(v0, _u)*dx )
    F2a = inner(S, grad(_v))*dx + pp*_p*dx
    F2b = inner(S0, grad(_v))*dx + pp0*_p*dx
    F2 = (1.0/dt)*inner(v-v0, _v)*dx + q*F2a + (1.0-q)*F2b

    # Traction at boundary
    F = I + grad(u)
    bF_magnitude = Constant(0.0)
    bF_direction = {2: Constant((0.0, 1.0)), 3: Constant((0.0, 0.0, 1.0))}[gdim]
    bF = det(F)*dot(inv(F).T, bF_magnitude*bF_direction)
    FF = inner(bF, _v)*ds

    # Whole system and its Jacobian
    F = F1 + F2 + FF
    J = derivative(F, w)

    # Initialize solver
    problem = NonlinearVariationalProblem(F, w, bcs=bcs, J=J)
    solver = NonlinearVariationalSolver(problem)
    solver.parameters['newton_solver']['relative_tolerance'] = 1e-6
    solver.parameters['newton_solver']['linear_solver'] = 'mumps'

    # Extract solution components
    u, v, p = w.split()
    u.rename("u", "displacement")
    v.rename("v", "velocity")
    p.rename("p", "pressure")

    # Create files for storing solution
    vfile = XDMFFile(os.path.join(output_dir, "velo.xdmf"))
    ufile = XDMFFile(os.path.join(output_dir, "disp.xdmf"))
    pfile = XDMFFile(os.path.join(output_dir, "pres.xdmf"))

    # Prepare plot window
    fig = plt.figure()
    fig.show()

    # Time-stepping loop
    t = 0
    while t <= T_end:
        t += float(dt)
        info("Time: {}".format(t))

        # Increase traction
        bF_magnitude.assign(100.0*t)

        # Prepare to solve and solve
        w0.assign(w)
        solver.solve()

        # Store solution to files and plot
        ufile.write(u, t)
        vfile.write(v, t)
        pfile.write(p, t)
        fig.clear()
        plot(u, mode="displacement")
        fig.canvas.draw()

    # Close files
    vfile.close()
    ufile.close()
    pfile.close()


def geometry_2d(length):
    """Prepares 2D geometry. Returns facet function with 1, 2 on parts of
    the boundary."""
    n = 5
    x0 = 0.0
    x1 = x0 + length
    y0 = 0.0
    y1 = 1.0
    mesh = RectangleMesh(Point(x0, y0), Point(x1, y1), int((x1-x0)*n), int((y1-y0)*n), 'crossed')
    boundary_parts = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
    left  = AutoSubDomain(lambda x: near(x[0], x0))
    right = AutoSubDomain(lambda x: near(x[0], x1))
    left .mark(boundary_parts, 1)
    right.mark(boundary_parts, 2)
    return boundary_parts


def geometry_3d():
    """Prepares 3D geometry. Returns facet function with 1, 2 on parts of
    the boundary."""
    mesh = Mesh('lego_beam.xml')
    gdim = mesh.geometry().dim()
    x0 = mesh.coordinates()[:, 0].min()
    x1 = mesh.coordinates()[:, 0].max()
    boundary_parts = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
    left  = AutoSubDomain(lambda x: near(x[0], x0))
    right = AutoSubDomain(lambda x: near(x[0], x1))
    left .mark(boundary_parts, 1)
    right.mark(boundary_parts, 2)
    return boundary_parts


if __name__ == '__main__':
    parameters['std_out_all_processes'] = False

    solve_elasticity(geometry_2d(20.0), 1e5, 0.3, 0.25, 5.0, 'results_2d_comp')
    solve_elasticity(geometry_2d(20.0), 1e5, 0.5, 0.25, 5.0, 'results_2d_incomp')
    solve_elasticity(geometry_2d(80.0), 1e5, 0.3, 0.25, 5.0, 'results_2d_long_comp')
    solve_elasticity(geometry_3d(),     1e5, 0.3, 0.50, 5.0, 'results_3d_comp')