Minimal example for plane elastostatics¶
This is a minimal example for plane elasticity. It will be used as point of departure for the remaining examples in this part.
We consider a plane problem of linear elasticity. Considering a square with given dimension \(h\), i.e., a domain \(\Omega = (0,h) \times (0,h)\), the governing field equations in \(\Omega\) are
with stress tensor \(\sigma\), strain tensor \(e\), displacement vector \(\mathbf{u}\), elasticity tensor \(C\), and volume force vector \(\mathbf{f}\).
In view of problems to be addressed later on, the square consists of three horizontal stripes. At the bottom \(\partial \Omega_\mathrm{B}\), the square is completely fixed, whereas at the left boundary (\(\partial \Omega_\mathrm{L}\)), only the horizontal displacement is blocked. In addition, a vertical displacement \(\bar{u}_2\) is prescribed at the top (\(\partial \Omega_\mathrm{T}\)). Furthermore, a volume force can be applied.
The corresponding variational problem reads
with test function \(\mathbf{v} \in V\) and \(\mathbf{v}=\mathbf{0}\) at \(\partial \Omega_\mathrm{B}\), \(v_1 = 0\) at \(\partial \Omega_\mathrm{L}\) as well as \(v_2 = 0\) at \(\partial \Omega_\mathrm{B}\), but otherwise arbitrary. Since, \(\mathbf{u}\) and \(\mathbf{v}\) are vector functions, \(V\) is a vector function space.
Note
You can download the complete python script
homPlaneElast.py
In the following, the main parts of the script are discussed in more detail. First, the necessary modules are imported. Different ways to import module or just specific functions of modules exist. Here, explicit imports are preferred. The DOLFINx version is printed, because there can be changes regarding the implementation depending on the version. At the moment DOLFINx version: 0.8.0 is used.
# test homogeneous plane elasticity
# ubuntu Release 18.04.6 LTS (Bionic Beaver) 64-bit
# docker version 24.0.2, build cb74dfc
import numpy as np
import dolfinx.fem.petsc
import ufl
from dolfinx import fem, io, mesh, default_scalar_type
from mpi4py import MPI
from petsc4py import PETSc
print("===============================================")
print(f"DOLFINx version: {dolfinx.__version__}")
print(f"based on GIT commit: {dolfinx.git_commit_hash}")
print(f"of https://github.com/FEniCS/dolfinx/")
print("===============================================")
Next, some general variables are defined, among them the names for the files used for storing the results.
#---------------------------------------------------------------------
# DATA SPECIFICATION
#---------------------------------------------------------------------
nr_elemsX1 = 9
nr_elemsX2 = 9
strip_height = 1.0
plate_height = 3.0*strip_height
YoungMod = 1.e5
PoissonNr = 0.3
VolumeForce = 0.0
res_filename = "HomPlaneElastResults.xdmf"
A mesh of quadrilateral elements defined on a rectangular domain is generated using a built-in function. Afterwards, a compatible vector function space \(V\) is defined. The vector character is explicit due to (gdim, ), since gdim=2.
#---------------------------------------------------------------------
# MESH
#---------------------------------------------------------------------
my_mesh = mesh.create_rectangle( comm=MPI.COMM_WORLD,
points=((0.0, 0.0),
(plate_height, plate_height)),
n=(nr_elemsX1, nr_elemsX2),
cell_type= mesh.CellType.quadrilateral)
gdim = my_mesh.geometry.dim
V = fem.functionspace(my_mesh, ("Lagrange", 1, (gdim, )))
Implementing homogeneous and inhomogeneous Dirichlet boundary conditions involves a number of steps:
Define the values for the components of the displacement vector to be prescribed.
Define functions to identify the corresponding boundary nodes (vertices) by means of a criterion.
Pass through the boundary facets (here edges) checking if the vertices of a facet fulfill the criterion.
Identify the corresponding degrees of freedom (dof’s). If only specific displacement components should be prescribed, the corresponding subspaces must be specified.
Assign the prescribed values to the dof’s, specifying again the corresponding subspaces if necessary.
Join all Dirichlet boundary conditions in a list.
#---------------------------------------------------------------------
# SUBDOMAINS DIRICHLET BC'S
#---------------------------------------------------------------------
u_Bottom = np.array([0, 0], dtype=default_scalar_type)
u_Top = -0.5
u_Left = 0.0
def Bottom(x):
return np.isclose(x[1], 0)
def Left(x):
return np.isclose(x[0], 0)
def Top(x):
return np.isclose(x[1], plate_height)
fdim = my_mesh.topology.dim - 1
Bottom_facets = mesh.locate_entities_boundary(my_mesh, fdim, Bottom)
Left_facets = mesh.locate_entities_boundary(my_mesh, fdim, Left)
Top_facets = mesh.locate_entities_boundary(my_mesh, fdim, Top)
Bottom_dofs = fem.locate_dofs_topological(V, fdim, Bottom_facets)
Left_dofs = fem.locate_dofs_topological(V.sub(0), fdim, Left_facets)
Top_dofs = fem.locate_dofs_topological(V.sub(1), fdim, Top_facets)
bcBottom = fem.dirichletbc( u_Bottom, Bottom_dofs, V)
bcLeft = fem.dirichletbc(u_Left, Left_dofs, V.sub(0))
bcTop = fem.dirichletbc(u_Top , Top_dofs , V.sub(1))
bcs = [bcBottom, bcLeft, bcTop]
Constituve relations are defined for plane stress and plane strain, working eventually with Lame’s constants based on Young’s modulus and Poisson’s number.
#---------------------------------------------------------------------
# CONSTITUTIVE LAW
#---------------------------------------------------------------------
E = fem.Constant(my_mesh, YoungMod)
nu = fem.Constant(my_mesh, PoissonNr)
model = "plane_stress"
mu_c = E/2/(1+nu)
lmbda_c = E*nu/(1+nu)/(1-2*nu)
if model == "plane_stress":
lmbda_c = 2*mu_c*lmbda_c/(lmbda_c+2*mu_c)
def epsilon(u):
return ufl.sym(ufl.grad(u))
# Stress function
def sigma(u):
return lmbda_c*ufl.nabla_div(u)*ufl.Identity(2) + 2*mu_c*epsilon(u)
The definition of the variational problem is rather straight forward, since the syntax is almost identical with the FEniCs syntax.
#---------------------------------------------------------------------
# VARIATIONAL PROBLEM
#---------------------------------------------------------------------
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(my_mesh, default_scalar_type((0, VolumeForce)))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx
petsc_opts={"ksp_type": "preonly", "pc_type": "lu"}
problem = fem.petsc.LinearProblem(a, L, bcs, petsc_options=petsc_opts)
uh = problem.solve()
After solving the problem, the results can be written into one or more files to be processed, for instantce with ParaView. Displacement results are given at the nodes of the mesh. Therefore, they can be written directly into the result file(s). Stresses on the other hand are computed at integration points and the results must therefore be processed before writing them into the output file(s).
#---------------------------------------------------------------------
# RESULTS
#---------------------------------------------------------------------
V0 = fem.functionspace(my_mesh, ("Lagrange", 1, (gdim, gdim)))
sig_exp = fem.Expression(sigma(uh), V0.element.interpolation_points())
sig = fem.Function(V0, name="Stress")
sig.interpolate(sig_exp)
out_ASCII = False
if out_ASCII:
encoding = dolfinx.io.XDMFFile.Encoding.ASCII
else:
encoding = dolfinx.io.XDMFFile.Encoding.HDF5
xdmffile = dolfinx.io.XDMFFile(MPI.COMM_WORLD, res_filename, "w", encoding)
xdmffile.write_mesh(my_mesh)
uh.name = "Displacement"
xdmffile.write_function(uh)
#sig.name = "Stress"
#xdmffile.write_function(sig)
xdmffile.close()