LargeDeformations: Rigid Body Rotation

This page is based on a Jupyter notebook.

In this benchmark we test a basic kinematic feature of the finite strain / large deformations implementation in OpenGeoSys. An element is subjected to a rigid body rotation. The expected result is a stress- and strain-free motion. This is confirmed for the finite-strain implementation while the small-strain implementation shows phantom strains and stresses.

Basic definitions

The small deformation code uses the linearized strain tensor:

\begin{align} \boldsymbol{\epsilon} = \frac{1}{2} \left( \text{grad}, \mathbf{u} + \text{grad},^\text{T} \mathbf{u} \right) \end{align}

while the large deformation code is set up in a Total Lagrangian formulation and rests on Green-Lagrange strains:

\begin{align} \mathbf{E} = \frac{1}{2} \left( \text{Grad}, \mathbf{U} + \text{Grad}^\text{T} \mathbf{U} + \text{Grad}^\text{T} \mathbf{U},\text{Grad}, \mathbf{U}\right) \end{align}

A rigid body rotation in 2D can be described by

\begin{align} \mathbf{u} = [X_1 (\cos \vartheta - 1) - X_2 \sin \vartheta] \mathbf{E}_1 + [X_1 \sin \vartheta + X_2 (\cos \vartheta - 1)]\mathbf{E}_2 \end{align}

While this yields $\mathbf{E} = \mathbf{0}$ in finite deformation kinematics, we obtain a linearized strain tensor coordinate matrix with

\begin{align} [\boldsymbol{\epsilon}]_{ij} = \left( \begin{array}{cc} \cos \vartheta - 1 & 0 \ 0 & \cos \vartheta - 1 \end{array} \right) \end{align}

For both cases the OGS’s linear elastic model is used to compute stresses. In the TL formulation this amounts to a Saint-Venant-Kirchhoff model. It thus suffices to illustrate the behaviour of strain values.

import os… (click to toggle)

import os
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import ogstools as ot
import pyvista as pv

# Some plot settings
plt.style.use("seaborn-v0_8-deep")
plt.rcParams["lines.linewidth"] = 2.0
plt.rcParams["lines.color"] = "black"
plt.rcParams["legend.frameon"] = True
plt.rcParams["font.family"] = "serif"
plt.rcParams["legend.fontsize"] = 14
plt.rcParams["font.size"] = 14
plt.rcParams["axes.spines.right"] = False
plt.rcParams["axes.spines.top"] = False
plt.rcParams["axes.spines.left"] = True
plt.rcParams["axes.spines.bottom"] = True
plt.rcParams["axes.axisbelow"] = True
plt.rcParams["figure.figsize"] = (8, 6)

out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))… (click to toggle)

out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
if not out_dir.exists():
    out_dir.mkdir(parents=True)

model_s = ot.Project(… (click to toggle)

model_s = ot.Project(
    input_file="square_1e0.prj", output_file=f"{out_dir}/square_1e0_small.prj"
)
model = ot.Project(input_file="square_1e0.prj", output_file="square_1e0.prj")

model_s.replace_text("SMALL_DEFORMATION", xpath="./processes/process/type")… (click to toggle)

model_s.replace_text("SMALL_DEFORMATION", xpath="./processes/process/type")
model_s.replace_text(
    "StandardElasticityBrick",
    xpath="./processes/process/constitutive_relation/behaviour",
)
model_s.replace_text("square_1e0_small", xpath="./time_loop/output/prefix")
model_s.remove_element(
    xpath="./processes/process/secondary_variables/secondary_variable[@internal_name='deformation_gradient']"
)
model_s.remove_element(
    xpath="./processes/process/secondary_variables/secondary_variable[@internal_name='volume_ratio']"
)
model_s.remove_element(xpath=".//vtkdiff[field='deformation_gradient']")
model_s.remove_element(xpath=".//vtkdiff[field='volume_ratio']")
model_s.write_input()

model.run_model(logfile=f"{out_dir}/out.txt", args=f"-o {out_dir}")… (click to toggle)

model.run_model(logfile=f"{out_dir}/out.txt", args=f"-o {out_dir}")
model_s.run_model(logfile=f"{out_dir}/out.txt", args=f"-o {out_dir} -m .")

OGS finished with project file square_1e0.prj.
Execution took 0.21545720100402832 s
Project file written to output.

OGS finished with project file /var/lib/gitlab-runner/builds/vZ6vnZiU/0/ogs/build/release-all/Tests/Data/LargeDeformation/RigidBody/RigidBody/square_1e0_small.prj.
Execution took 0.20690560340881348 s
Project file written to output.

pv.set_plot_theme("document")… (click to toggle)

pv.set_plot_theme("document")
pv.set_jupyter_backend("static")  # comment out for interactive graphics

reader = pv.get_reader(f"{out_dir}/square_1e0.pvd")… (click to toggle)

reader = pv.get_reader(f"{out_dir}/square_1e0.pvd")
reader_s = pv.get_reader(f"{out_dir}/square_1e0_small.pvd")

reader.set_active_time_value(0.0)… (click to toggle)

reader.set_active_time_value(0.0)
reader_s.set_active_time_value(0.0)
mesh = reader.read()[0]  # nulltes Gitter lesen
mesh_s = reader_s.read()[0]  # nulltes Gitter lesen

points = mesh.point_data["epsilon"].shape[0]

xs = mesh.points[:, 0]… (click to toggle)

xs = mesh.points[:, 0]
ys = mesh.points[:, 1]

def ploteps(time, angle, quantity="epsilon"):… (click to toggle)

def ploteps(time, angle, quantity="epsilon"):
    reader.set_active_time_value(time)
    reader_s.set_active_time_value(time)
    mesh = reader.read()[0]  # nulltes Gitter lesen
    mesh_s = reader_s.read()[0]  # nulltes Gitter lesen
    eps_vec = mesh.point_data[quantity][:, 0]
    eps_vec_s = mesh_s.point_data[quantity][:, 0]

    print(
        "Expected: %.2f (small strain) and 0 (large strain)"
        % (np.cos(np.deg2rad(angle)) - 1)
    )

    sargs = {
        "title": "small deformation, " + str(angle) + "°",
        "title_font_size": 15,
        "label_font_size": 15,
        "n_labels": 2,
        "position_x": 0.2,
        "position_y": 0.85,
        "fmt": "%.1f",
        "width": 0.6,
    }

    p = pv.Plotter(shape=(1, 2), border=False)
    p.subplot(0, 0)
    p.add_mesh(
        mesh,
        scalars=eps_vec_s,
        show_edges=False,
        show_scalar_bar=True,
        colormap="RdBu_r",
        scalar_bar_args=sargs,
    )
    # p.show_bounds(ticks="outside", xlabel="", ylabel="")
    # p.add_axes()
    p.view_xy()
    p.camera.zoom(1.2)

    sargs1 = {
        "title": "large deformation, " + str(angle) + "°",
        "title_font_size": 15,
        "label_font_size": 15,
        "n_labels": 2,
        "position_x": 0.2,
        "position_y": 0.85,
        "fmt": "%.1f",
        "width": 0.6,
    }

    p.subplot(0, 1)
    p.add_mesh(
        mesh,
        scalars=eps_vec,
        show_edges=False,
        show_scalar_bar=True,
        colormap="RdBu_r",
        scalar_bar_args=sargs1,
    )
    # p.show_bounds(ticks="outside", xlabel="", ylabel="")
    # p.add_axes()
    p.view_xy()
    p.camera.zoom(1.2)

    p.window_size = [800, 400]
    p.show()

We plot the normal strain in the $x$-direction for both kinematic formulations on the undeformed configuration as it undergoes a 360° rotation and find our expectations confirmed.

ploteps(0, 0)

Expected: 0.00 (small strain) and 0 (large strain)

/var/lib/gitlab-runner/builds/vZ6vnZiU/0/ogs/build/release-all/.venv/lib/python3.12/site-packages/pyvista/plotting/mapper.py:669: RuntimeWarning: All-NaN axis encountered
  clim = [np.nanmin(scalars), np.nanmax(scalars)]

png

ploteps(0.5, 45)

Expected: -0.29 (small strain) and 0 (large strain)

png

ploteps(1, 90)

Expected: -1.00 (small strain) and 0 (large strain)

png

ploteps(2, 180)

Expected: -2.00 (small strain) and 0 (large strain)

png

ploteps(3, 270)

Expected: -1.00 (small strain) and 0 (large strain)

png

ploteps(4, 360)

Expected: 0.00 (small strain) and 0 (large strain)

png

This article was written by Thomas Nagel. If you are missing something or you find an error please let us know.
Generated with Hugo 0.122.0 in CI job 504124 | Last revision: March 1, 2023