In this tutorial the implementation of linear elasticity problems in 2D is shown. We consider a thin plate with a hole at its centre. The plate is fixed at its left face and and there is some applied force on its right face. This is shown in the following figure. First we will give a short overview of the equations of linear elasticity in three dimensions and then reduce the problem to two dimensions.
The equation to be solved is
for the unknown displacement where .
The applied body forces are described by the mapping , called the density of the applied body force per unit volume. The tensor is the stress tensor, which will be described later.
On a subset the body is fixed and we get the homogenous Dirichlet boundary condition
On the subset we have applied surface forces, which are given by a mapping , called the density of the applied surface force per unit area. From that we get the boundary condition
where is the unit outer normal vector along .
The linear stress tensor is given by (also known as Hooke's Law)
with the linear strain tensor
Reduction to 2-dimensional elasticity
For the reduction to a 2-dimensional problem there are two kinds of reduction techniques: the plane stress state and the plane strain state. The plane stress state is a good approximation for very thin bodies. In this case one can assume that all components of connected to the -direction vanish, i.e.
By throwing away the component of in the -direction and expressing by and , we can write the material law from the first section in the form
For the plane strain state one can assume that all components of connected to the -direction vanish. Again by throwing away the component of in the -direction one can then write the material law in the form
Here we will use the plane stress state.
Variational Formulation
We choose
where now . Then we have the bilinearform
with . Here we use for matrices .
Finally, the variational formulation is: Find such that
in the case of the plane stress state and
in the case of the plane strain state. We will use this splitting of the bilinear form for computing the system matrix later.
are included. In this example we use a cig-file as input. This file is create with a Python script and contains geometry data, the densities of the applied forces, material constants and values for the h- and p-refinement. For using that cig-file as input, we write a method that first sets some default values and then reads the file.
MeshEPS< Real > drawMeshEPS(concepts::Mesh &msh, std::string filename, const Real scale=100, const Real greyscale=1.0, const unsigned int nPoints=2)
We load the attribute of the edges with the Dirichlet boundary condition from the input file and set the boundary condition for all edges with this attribute.
void setupDivDiv(vectorial::BilinearForm< F, typename concepts::Realtype< F >::type > &bf, const concepts::ElementFormulaContainer< F, typename concepts::Realtype< F >::type > frm=concepts::ElementFormulaContainer< F, typename concepts::Realtype< F >::type >())
The right hand side vectors are computed from the applied volume and surface forces. We compute these vectors for each component of the forces and concatenate them to one vector for the vectorial system.
Now we write the solution (the displacement) and the resulting strains and stresses to a matlab binary file. Before doing that we recompute the shape functions on equidistant points.
The output of the solutions of the program is saved in a file called elasticity2D_stress.mat. For regarding the output (for example the stress ) in MATLAB we can type the following in the MATLAB command window
First we will give a short overview of the equations of linear elasticity in three dimensions and then reduce the problem to two dimensions.![\[
- {\rm div} \:\sigma(u) = f \quad \text{in } \Omega,
\]](form_73_dark.png)
![\[
u = 0 \quad \text{on } \Gamma_0.
\]](form_79_dark.png)
![\[
\sigma(u) n = g \quad \text{on } \Gamma_1,
\]](form_82_dark.png)
![\[
\sigma(u) = \lambda\; {\rm trace} (\varepsilon(u)) + 2\mu \varepsilon(u).
\]](form_85_dark.png)
![\[
\varepsilon_{ij}(u) = \frac{1}{2} (\partial_j u_i + \partial_i u_j ).
\]](form_86_dark.png)
![\[
\sigma_{i3}(u) = \sigma_{3i}(u) = 0 \quad i=1,2,3.
\]](form_88_dark.png)
![\[
\sigma^{\rm 2D} (u) = 2\mu \frac{\lambda}{\lambda+2\mu}\; {\rm trace}\left(\varepsilon^{\rm 2D}(u)\right) + 2\mu \varepsilon^{\rm 2D}(u).
\]](form_92_dark.png)
![\[
\sigma^{\rm 2D} (u) = \lambda\; {\rm trace}\left(\varepsilon^{\rm 2D}(u)\right) + 2\mu \varepsilon^{\rm 2D}(u).
\]](form_94_dark.png)
![\[
V := \{ v \in H^1(\Omega)^2\,|\,\text{with } v = 0 \text{ on } \Gamma_0 \}
\]](form_95_dark.png)
![\[
B(u,v) = \int_\Omega \sigma^{\rm 2D}(u):\varepsilon^{\rm 2D}(v) \;{\rm d}x
\]](form_97_dark.png)
![\[
2\mu \frac{\lambda}{\lambda+2\mu} \int_\Omega {\rm div}(u){\rm div}(v) \;{\rm d}x + 2\mu \int_\Omega \varepsilon^{\rm 2D}(u): \varepsilon^{\rm 2D}(v) \;{\rm d}x = \int_\Omega f \cdot v \;{\rm d}x + \int_{\Gamma_1} g \cdot v \;{\rm d}s \quad \forall v \in V
\]](form_101_dark.png)
![\[
\lambda \int_\Omega {\rm div}(u){\rm div}(v) \;{\rm d}x + 2\mu \int_\Omega \varepsilon^{\rm 2D}(u): \varepsilon^{\rm 2D}(v) \;{\rm d}x = \int_\Omega f \cdot v \;{\rm d}x + \int_{\Gamma_1} g \cdot v \;{\rm d}s \quad \forall v \in V
\]](form_102_dark.png)
1.9.8