Revision as of 13:00, 8 November 2016

Basic equations of elasticity

To determine the distribution of static stresses and displacements in a solid body we must obtain a solution (either analytically or numerically) to the basic equations of the theory of elasticity, satisfying the boundary conditions on forces and/or displacements. The equations thus form a boundary value problem. For a general three dimensional solid object the equations governing its behavior are:

• equations of equilibrium (3)
• strain-displacement equations (6)
• stress-strain equations (6)

where the number in brackets indicates the number of equations. The equations of equilibrium are three tensor partial differential equations for the balance of linear momentum, the strain-displacement equations are relations that stem from infinitesimal strain theory and the stress-strain equations are a set of linear algebraic constitutive relations (3D Hooke's law). In two dimensions these equations simplify to 8 equations (2 equilibrium, 3 strain-displacement, and 3 stress-strain).

A large amount of confusion regarding linear elasticity originates from the many different notations used in this field. For this reason we first provide an overview of the different notations that are used.

Direct tensor form

In direct tensor form (independent of coordinate system) the governing equations are:

• Equation of motion (an expression of Newton's second law):

\[ \boldsymbol{\nabla} \cdot \boldsymbol{\sigma} + \boldsymbol{F} = \rho \ddot{\boldsymbol{u}} \]

• Strain-displacement equations:

\[ \boldsymbol{\varepsilon} = \frac{1}{2}\left[ \boldsymbol{\nabla} \boldsymbol{u} + (\boldsymbol{\nabla} \boldsymbol{u})^T \right] \]

• Stress-displacement equations (constitutive equations). For a linear elastic material this is Hooke's law:

\[ \boldsymbol{\sigma} = \boldsymbol{C} : \boldsymbol{\varepsilon} \]

where $\boldsymbol{\sigma}$ is the Cauchy stress tensor, $\boldsymbol{\varepsilon}$ is the infinitesimal strain tensor, $\boldsymbol{u}$ is the displacement vector, $\boldsymbol{C}$ is the fourth order stiffness tensor, $\boldsymbol{F}$ is the body force per unit volume (a vector quantity), $\rho$ is the mass density, $\boldsymbol{\nabla}(\bullet)$ is the gradient operator, $\boldsymbol{\nabla}\cdot(\bullet)$ is the divergence operator, $(\bullet)^T$ represents a transpose, $\ddot{(\bullet)}$ is the second derivative with respect to time, and $\boldsymbol{A}:\boldsymbol{B}$ is the inner product of two second order tensors (a tensor contraction).

Cartesian coordinate form

Using the Einstein summation convention (implied summations over repeated indexes) the equations are:

• Equation of motion (an expression of Newton's second law):

\[ \sigma_{ji,j} + F_i = \rho \partial_{tt} u_i \]

• Strain-displacement equations:

\[ \varepsilon_{ij} = \frac{1}{2}\left(u_{j,i}+u_{i,j}\right) \]

• Stress-displacement equations (constitutive equations). For a linear elastic material this is Hooke's law:

\[ \sigma_{i,j} = C_{ijkl} \varepsilon_{kl}, \] where $i,j = 1,2,3$ represent, respectively, $x$, $y$ and $z$, the $(\bullet),j$ subscript is a shorthand for partial derivative $\partial(\bullet)/\partial x_j$ and $\partial_{tt}$ is shorthand notation for $\partial^2/\partial t^2$, $\sigma_{ij}=\sigma_{ji}$ is the Cauchy stress tensor (with 6 independent components), $F_i$ are the body forces, $\rho$ is the mass density, $u_i$ is the displacement, $\varepsilon_{ij} = \varepsilon_{ji}$ is the strain tensor (also with 6 independent components), and, finally $C_{ijkl}$ is the fourth-order stiffness tensor that due to symmetry requirements $C_{ijkl} = C_{klij} = C_{jikl} =C_{ijkl}$ can be reduced to 21 different elements.

Matrix-vector (FEM) notation

• Equation of motion:

\[ \boldsymbol{L}^T\boldsymbol{\sigma} + \boldsymbol{F} = \rho\ddot{\boldsymbol{u}} \]

• Strain-displacement equations:

\[ \boldsymbol{\varepsilon} = \boldsymbol{L}\boldsymbol{u} \]

• Stress-displacement equations (constitutive equations). For a linear elastic material this is Hooke's law:

\[ \boldsymbol{\sigma} = \boldsymbol{C}\boldsymbol{\varepsilon} \] where $\boldsymbol{\sigma}$ is the Cauchy stress tensor represented in vector form (6 components), $\boldsymbol{L}$ is a differential operator matrix (size $3 \times 6$), $\bullet^T$ is transpose of a matrix $\bullet$, $\boldsymbol{F}$ is the body force vector, $\rho$ is the mass density, $\boldsymbol{u}$ is the displacement vector, $\boldsymbol{\varepsilon}$ is the strain tensor represented in vector form (6 components), and $\boldsymbol{C}$ is the symmetric stress-strain matrix (size $6 \times 6$ with 21 material constants $C_{ij} = C_{ji}$). Certain literature prefers using the symbol $\boldsymbol{D}$ instead of $\boldsymbol{C}$ for the stress-strain matrix. The symbol $C$ is then used as the "compliance tensor" that relates the strains to the stresses, e.g. $\boldsymbol{\varepsilon} = \boldsymbol{C}\boldsymbol{\sigma}$.

To solve the basic equations of elasticity two approaches exist according to the boundary conditions of the boundary value problem. In the displacement formulation the displacements are prescribed everywhere at the boundaries and the stresses and strains are eliminated from the formulation. The other possible option is that the surface tractions are prescribed everywhere on the surface boundary. The equations of elasticity are then manipulated to leave the stresses as the unknown to be solved for. This approach is known as the stress formulation and will not be considered here further.

Displacement formulation

The goal of the displacement formulation is to eliminate the strains and stresses from the basic equations of elasticity and leave only the displacements as the unknown to be solved for. The first step is to substitute the strain-displacement equations into Hooke's law, eliminating the strains as unknowns: \begin{equation}\label{eq:3d_hooke} \sigma_{ij} = \lambda\delta_{ij}\varepsilon_{kk} + 2\mu\varepsilon_{ij} \quad \Rightarrow \quad \begin{cases} \sigma_{xx} = \lambda (u_{x,x}+u_{y,y}+u_{z,z}) + 2\mu u_{x,x}\\ \sigma_{yy} = \lambda (u_{x,x}+u_{y,y}+u_{z,z}) + 2\mu u_{y,y}\\ \sigma_{zz} = \lambda (u_{x,x}+u_{y,y}+u_{z,z}) + 2\mu u_{z,z}\\ \sigma_{xy} = \mu (u_{x,y} + u_{y,x})\\ \sigma_{yz} = \mu (u_{y,z} + u_{z,y})\\ \sigma_{zx} = \mu (u_{z,x} + u_{x,z}) \end{cases} \end{equation} where $\lambda$ and $\mu$ are Lamé parameters. The formula above can also be written more concisely as \[ \sigma_{ij} = \lambda \delta_{ij} u_{k,k} + \mu(u_{i,j}+u_{j,i}). \] The next step is to substitute these equations into the equilibrium equation \[ \sigma_{ij,j} + F_i = \rho \partial_{tt} u_i \quad \Rightarrow \quad \begin{cases} \frac{\partial \sigma_{xx}}{\partial x} + \frac{\partial \sigma_{xy}}{\partial y} + \frac{\partial \sigma_{xz}}{\partial z} + F_x = \rho \frac{\partial^2 u_x}{\partial t^2} \\ \frac{\partial \sigma_{yx}}{\partial x} + \frac{\partial \sigma_{yy}}{\partial y} + \frac{\partial \sigma_{yz}}{\partial z} + F_y = \rho \frac{\partial^2 u_y}{\partial t^2} \\ \frac{\partial \sigma_{zx}}{\partial x} + \frac{\partial \sigma_{zy}}{\partial y} + \frac{\partial \sigma_{zz}}{\partial z} + F_z = \rho \frac{\partial^2 u_z}{\partial t^2} \end{cases} \] resulting in \begin{align*} \frac{\partial}{\partial x}\left( \lambda \left( \frac{\partial u_x}{\partial x} + \frac{\partial u_y}{\partial y} + \frac{\partial u_z}{\partial z} \right) + 2\mu \frac{\partial u_x}{\partial x} \right) + \frac{\partial}{\partial y}\left(\mu\left(\frac{\partial u_x}{\partial y}+\frac{\partial u_y}{\partial x}\right)\right) + \frac{\partial}{\partial z}\left(\mu\left(\frac{\partial u_x}{\partial z}+\frac{\partial u_z}{\partial x}\right)\right) + F_x &= \rho \frac{\partial^2 u_x}{\partial t^2} \\ \frac{\partial}{\partial x}\left(\mu\left(\frac{\partial u_y}{\partial x}+\frac{\partial u_x}{\partial y}\right)\right) + \frac{\partial}{\partial y}\left(\lambda \left( \frac{\partial u_x}{\partial x} + \frac{\partial u_y}{\partial y} + \frac{\partial u_z}{\partial z} \right) + 2\mu \frac{\partial u_y}{\partial y}\right) + \frac{\partial}{\partial z}\left(\mu\left(\frac{\partial u_y}{\partial z}+\frac{\partial u_z}{\partial y}\right)\right) + F_y &= \rho \frac{\partial^2 u_y}{\partial t^2}\\ \frac{\partial}{\partial x}\left(\mu\left(\frac{\partial u_z}{\partial x}+\frac{\partial u_x}{\partial z}\right)\right) + \frac{\partial}{\partial y}\left(\mu\left(\frac{\partial u_z}{\partial y}+\frac{\partial u_y}{\partial z}\right)\right) + \frac{\partial}{\partial z}\left(\lambda \left( \frac{\partial u_x}{\partial x} + \frac{\partial u_y}{\partial y} + \frac{\partial u_z}{\partial z} \right) + 2\mu \frac{\partial u_z}{\partial z}\right) + F_z &= \rho \frac{\partial^2 u_z}{\partial t^2} \end{align*} Using the assumption that Lamé parameters $\lambda$ and $\mu$ are constant we can rearrange to produce: \begin{align*} (\lambda + \mu) \frac{\partial}{\partial x}\left(\frac{\partial u_x}{\partial x}+\frac{\partial u_y}{\partial y}+\frac{\partial u_z}{\partial z}\right) + \mu\left(\frac{\partial^2 u_x}{\partial x^2}+\frac{\partial^2 u_x}{\partial y^2}+\frac{\partial^2 u_x}{\partial z^2}\right) + F_x &= \rho \frac{\partial^2 u_x}{\partial t^2} \\ (\lambda + \mu) \frac{\partial}{\partial y}\left(\frac{\partial u_x}{\partial x}+\frac{\partial u_y}{\partial y}+\frac{\partial u_z}{\partial z}\right) + \mu\left(\frac{\partial^2 u_y}{\partial x^2}+\frac{\partial^2 u_y}{\partial y^2}+\frac{\partial^2 u_y}{\partial z^2}\right) + F_y &= \rho \frac{\partial^2 u_y}{\partial t^2} \\ (\lambda + \mu) \frac{\partial}{\partial z}\left(\frac{\partial u_x}{\partial x}+\frac{\partial u_y}{\partial y}+\frac{\partial u_z}{\partial z}\right) + \mu\left(\frac{\partial^2 u_z}{\partial x^2}+\frac{\partial^2 u_z}{\partial y^2}+\frac{\partial^2 u_z}{\partial z^2}\right) + F_z &= \rho \frac{\partial^2 u_z}{\partial t^2} \\ \end{align*}

We can easily see that only the displacements are left in these equations. The equations obtaines in this manner are known as the Navier-Cauchy equations.

Navier-Cauchy equations

The Navier or Navier-Cauchy equations describe the dynamics of a solid through the displacement vector field $\b{u}$. The equations can also be expressed concisely in vector form as follows \begin{equation}\label{eq:navier_mu_lambda} \rho \frac{\partial^2 \b{u}}{\partial t^2} = (\lambda + \mu) \nabla(\nabla \cdot \b{u}) + \mu \nabla^2 \b{u} + \b{F} \end{equation} where $\mu$ and $\lambda$ are Lamé constants, $\rho$ is the object density and $\b{F}$ are the external forces. In certain cases we may prefer to use the Young modulus $E$ and the Poisson ratio $\nu$ instead of the Lamé constants. In this case the Navier equations are \begin{equation}\label{eq:navier_young_poisson} \rho \frac{\partial^2 \b{u}}{\partial t^2} = \frac{E}{2(1+\nu)}\left(\nabla^2 \b{u} + \frac{1}{1-2\nu}\nabla\left(\nabla \cdot \b{u}\right)\right) + \b{F} \end{equation}

Two-dimensional stress distributions

Many problems in elasticity can be simplified as two-dimensional problems described by plane theory of elasticity. In general there are two types of problems we may encounter in plane analysis: plane stress and plane strain. The first problem arises in analysis of thin plates loaded in the plane of the plate, while the second is used for elongated bodies of constant cross section subject to uniform loading.

Plane stress

Plane stress distributions build on the assumption that the normal stress and shear stresses directed perpendicular to the $x$-$y$ plane are assumed zero: \begin{equation}\label{eq:pstress_assump} \sigma_{zz} = \sigma_{zx} = \sigma_{zy} = 0. \end{equation} It is also assumed that the stress components do not vary through the thickness of the plate (the assumptions do violate some compatibility conditions, but are still sufficiently accurate for practical applications if the plate is thin).

Using (\ref{eq:pstress_assump}) the three-dimensional Hooke's law can be reduced to: \begin{equation}\label{eq:planestress} \boldsymbol{\sigma} = \boldsymbol{C}\boldsymbol{\varepsilon} \end{equation} in matrix form where for isotropic materials, we have \begin{equation}\label{eq:planestressmatrix} \boldsymbol{C} = \frac{E}{1-\nu^2} \begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1-\nu}{2} \end{bmatrix} \qquad \text{(Plane stress)} \end{equation} and \begin{equation}\label{eq:2d_stress} \b{\sigma} = \begin{bmatrix} \sigma_{xx} \\ \sigma_{yy} \\ \sigma_{xy}\end{bmatrix}, \end{equation} \begin{equation}\label{eq:2d_strain} \b{\varepsilon} = \begin{bmatrix} \varepsilon_{xx} \\ \varepsilon_{yy} \\ \gamma_{xy} \end{bmatrix} = \begin{bmatrix} \varepsilon_{xx} \\ \varepsilon_{yy} \\ 2\varepsilon_{xy} \end{bmatrix} = \begin{bmatrix} \frac{\partial u_x}{\partial x} \\ \frac{\partial u_y}{\partial y} \\ \frac{\partial u_y}{\partial x} + \frac{\partial u_x}{\partial y} \end{bmatrix}. \end{equation} Note that the strain tensor uses the so-called engineering shear strain $\gamma_{ij}$. Under the plane stress assumption the Navier equations are given as: \begin{equation}\label{eq:navier_plane_stress} \rho \frac{\partial^2 \b{u}}{\partial t^2} = \frac{E}{2(1+\nu)} \left( \b{\nabla}^2 \b{u} + \frac{1+\nu}{1-\nu}\b{\nabla}\left( \b{\nabla}\cdot \b{u}\right) \right) + \b{F} \end{equation} or in component notation \begin{align} \rho \frac{\partial^2 u_x}{\partial t^2} &= \frac{E}{2(1+\nu)} \left( \frac{\partial^2 u_x}{\partial x^2} + \frac{\partial^2 u_x}{\partial y^2} \right) + \frac{E}{2(1-\nu)} \frac{\partial}{\partial x}\left( \frac{\partial u_x}{\partial x} + \frac{\partial u_y}{\partial y} \right) \\ \rho \frac{\partial^2 u_y}{\partial t^2} &= \frac{E}{2(1+\nu)} \left( \frac{\partial^2 u_y}{\partial x^2} + \frac{\partial^2 u_y}{\partial y^2} \right) + \frac{E}{2(1-\nu)} \frac{\partial}{\partial y}\left( \frac{\partial u_x}{\partial x} + \frac{\partial u_y}{\partial y} \right) \end{align}

Plane strain

The plane strain problem arises in analysis of walls, dams, tunnels where one dimension of the structure is very large in comparison to the other two dimensions ($x$- and $y$- coordinates). It is also appropriate for small-scale problems such as bars and rollers compressed by forces normal to their cross section. In all such problems the body may be imagined as a prismatic cylinder with one dimension much larger that the other two. The applied forces act in the $x$-$y$ plane and do not vary in the $z$ direction, leading to the assumption \begin{equation} \frac{\partial}{\partial z} = u_z = 0. \end{equation}

With the above assumption it follows immediately that \begin{equation} \varepsilon_{zz} = \varepsilon_{zx} = \varepsilon_{zy} = 0. \end{equation} The three dimensional Hooke's law can now be reduced to \begin{equation}\label{eq:planestrain} \boldsymbol{\sigma} = \boldsymbol{C}\boldsymbol{\varepsilon} \end{equation} where the matrix $C$ is given by \begin{equation}\label{eq:matrixplanestrain} \boldsymbol{C} = \frac{E}{(1+\nu)(1-2\nu)} \begin{bmatrix} 1-\nu & \nu & 0 \\ \nu & 1-\nu & 0 \\ 0 & 0 & \frac{1-2\nu}{2} \end{bmatrix} \qquad \text{(Plane strain)} \end{equation} The vectors $\boldsymbol{\sigma}$ and $\boldsymbol{\varepsilon}$ are the same as above for plane stress and are given in (\ref{eq:2d_stress}) and (\ref{eq:2d_strain}), respectively. In the case of plane strain we have additonal non-zero components of the stress tensor: \begin{equation}\label{eq:sigmazz} \sigma_{zz} = \nu(\sigma_{xx}+\sigma_{yy}) \end{equation} \begin{equation} \sigma_{yz} = \sigma_{zx} = 0 \end{equation} The reason $\sigma_{zz}$ is not included in the matrix stress-strain equation (\ref{eq:planestrain}) is because it is linearly dependent on the normal stresses $\sigma_{xx}$ and $\sigma_{zz}$.

The Navier-Cauchy equation that follows from the plane strain assumption can be found in (\ref{eq:navier_young_poisson}). For purposes of completeness we provide here the equations in component notation: \begin{align} \rho \frac{\partial^2 u_x}{\partial t^2} &= \frac{E}{2(1+\nu)} \left( \frac{\partial^2 u_x}{\partial x^2} + \frac{\partial^2 u_x}{\partial y^2} \right) + \frac{E}{2(1+\nu)(1-2\nu)} \frac{\partial}{\partial x}\left( \frac{\partial u_x}{\partial x} + \frac{\partial u_y}{\partial y} \right) \\ \rho \frac{\partial^2 u_y}{\partial t^2} &= \frac{E}{2(1+\nu)} \left( \frac{\partial^2 u_y}{\partial x^2} + \frac{\partial^2 u_y}{\partial y^2} \right) + \frac{E}{2(1+\nu)(1-2\nu)} \frac{\partial}{\partial y}\left( \frac{\partial u_x}{\partial x} + \frac{\partial u_y}{\partial y} \right) \end{align}

Plane strain with Lamé constants

Alternatively the stiffness tensor $\b{C}$ in the stress-strain equation $\b{\sigma} = \b{C}\b{\varepsilon}$ may be expressed in terms of $\mu$ and $\lambda$: \begin{equation} \boldsymbol{C} = \begin{bmatrix} 2\mu + \lambda & \lambda & 0 \\ \lambda & 2\mu + \lambda & 0 \\ 0 & 0 & 2\mu \end{bmatrix} \qquad \text{(Plane strain)} \end{equation} A word of caution applies here. When using the stiffness tensor in this form the strain tensor must be given as: \[ \b{\varepsilon} = \{\varepsilon_{xx} \quad \varepsilon_{yy} \quad \varepsilon_{xy}\}^T \] WITHOUT the factor 2 in front of $\varepsilon_{xy}$ (see equation (\ref{eq:2d_strain})).

Connection between plane stress and plane strain

For isotropic materials with elastic modulus $E$ and Poisson's ratio $\nu$ it is possible to go from plane stress to plane strain, or vice-versa, by replacing $E$ and $\nu$ in the stress-strain matrix with a fictitious modulus $E^*$ and fictitious Poisson ratio $\nu^*$. This allows us to "reuse" a plane stress program to solve plane strain or again vice-versa (as long as the material is isotropic). A few exercises on this topic are given at this link (page 13).

To go from plane stress ($s$) to plane strain ($n$) insert the fictitious quantities \begin{equation}\label{eq:ston1} E_n^* = \frac{E_s}{1-\nu_s^2}, \end{equation} \begin{equation}\label{eq:ston2} \nu_n^* = \frac{\nu_s}{1-\nu_s}. \end{equation}

In the opposite case we want to go from plane strain ($n$) to plane stress ($s$) we can use: \begin{equation}\label{eq:ntos1} E_s^* = \frac{E_n(1+2\nu_n)}{(1+\nu_n)^2}, \end{equation} \begin{equation}\label{eq:ntos2} \nu_s^* = \frac{\nu_n}{1+\nu_n}. \end{equation}

Point contact on a 2D half-plane

A starting point to solve problems in contact mechanics is to understand the effect of a point-load applied to a homogeneous, linear elastic, isotropic half-plane. This problem may be defined either as plane stress or plain strain (for solution with FreeFem++ we have choosen the latter). The traction boundary conditions for this problem are: \begin{equation}\label{eq:bc} \sigma_{xy}(x,0) = 0, \quad \sigma_{yy}(x,y) = -P\delta(x,y) \end{equation} where $\delta(x,y)$ is the Dirac delta function. Together these boundary conditions state that there is a singular normal force $P$ applied at $(x,y) = (0,0)$ and there are no shear stresses on the surface of the elastic half-plane.

The analytical relations for the stresses can be found from the Flamant solution (stress distributions in a linear elastic wedge loaded by point forces a the tip. When the "wedge" is flat we get a half-plane. The derivation uses polar coordinates.) and are given as: \begin{equation} \sigma_{xx} = -\frac{2P}{\pi} \frac{x^2y}{\left(x^2+y^2\right)^2}, \end{equation} \begin{equation} \sigma_{yy} = -\frac{2P}{\pi} \frac{y^3}{\left(x^2+y^2\right)^2}, \end{equation} \begin{equation} \sigma_{xy} = -\frac{2P}{\pi} \frac{xy^2}{\left(x^2+y^2\right)^2}, \end{equation} for some point $(x,y)$ in the half-plane. From this stress field the strain components and thus the displacements $(u_x,u_y)$ can be determined. The displacements are given by \begin{align} u_x &= -\frac{P}{4\pi\mu}\left((\kappa-1)\theta - \frac{2xy}{r^2}\right), \label{eq:dispx}\\ u_y &= -\frac{P}{4\pi\mu}\left((\kappa+1)\log r + \frac{2x^2}{r^2}\right), \label{eq:dispy} \end{align} where $$r = \sqrt{x^2+y^2}$$ and $$\tan \theta = \frac{x}{y}.$$ The symbol $\kappa$ is known as Dundars constant and is defined as \[ \kappa = \begin{cases} 3 - 4\nu & \quad \text{(Plane strain)}, \\ \cfrac{3 - \nu}{1+\nu} & \quad \text{(Plane stress)}. \end{cases} \] The last remaining symbol is $\mu$ which represents the shear modulus (sometimes also denoted with $G$).

Numerical solution with FreeFem++

Due to the known analytical solution the point-contact problem can be used for benchmarking numerical PDE solvers in terms of accuracy (as well as computational efficiency). The purpose of this section is to compare the numerical solution obtained by FreeFem++ with the analytical solution, as well as provide a reference numerical solution for the C++ library developed in our laboratory.

For purposes of simplicity we limit ourselves to the domain $(x,y) \in \Omega = [-1,1] \times[-1,-0.1]$ and prescribe Dirichlet displacement on the boundaries $\Gamma_D$ from the known analytical solution (\ref{eq:dispx}, \ref{eq:dispy}). This way we avoid having to deal with the Dirac delta traction boundary condition (\ref{eq:bc}). The problem can be described as find $\boldsymbol{u(\boldsymbol{x})}$ that satisfies \begin{equation} \boldsymbol{\nabla}\cdot\boldsymbol{\sigma}= 0 \qquad \text{on }\Omega \end{equation} and \begin{equation} \boldsymbol{u} = \boldsymbol{u}_{\text{analytical}} \qquad \text{on }\Gamma_D \end{equation} where $\boldsymbol{u}_\text{analytical}$ is given in equations (\ref{eq:dispx}) and (\ref{eq:dispy}).

To solve the point-contact problem in FreeFem++ we must first provide the weak form of the balance equation: \begin{equation*} \boldsymbol{\nabla}\cdot\boldsymbol{\sigma} + \boldsymbol{b} = 0. \end{equation*} The corresponding weak formulation is \begin{equation}\label{eq:weak} \int_\Omega \boldsymbol{\sigma} : \boldsymbol{\varepsilon}(\boldsymbol{v}) \, d\Omega - \int_\Omega \boldsymbol{b}\cdot\boldsymbol{v}\,d\Omega = 0 \end{equation} where $:$ denotes the tensor scalar product (tensor contraction), i.e. $\boldsymbol{A}:\boldsymbol{B} =\sum_{i,j} A_{ij}B_{ij}$. The vector $\boldsymbol{v}$ is the test function or so-called "virtual displacement".

Equation (\ref{eq:weak}) can be handed to FreeFem++ with the help of Voigt or Mandel notation, that reduces the symmetric tensors $\boldsymbol{\sigma}$ and $\boldsymbol{\varepsilon}$ to vectors. The benefit of Mandel notation is that it allows the tensor scalar product to be performed as a scalar product of two vectors. For this reason we create the following macros:

macro u [ux,uy] // displacements
macro v [vx,vy] // test function
macro b [bx,by] // body forces
macro e(u) [dx(u[0]),dy(u[1]),(dx(u[1])+dy(u[0]))/2] // strain (for post-processing)
macro em(u) [dx(u[0]),dy(u[1]),sqrt(2)*(dx(u[1])+dy(u[0]))/2] // strain in Mandel notation
macro A [[2*mu+lambda,mu,0],[mu,2*mu+lambda,0],[0,0,2*mu]] // stress-strain matrix

The weak form (\ref{eq:weak}) can then be expressed naturally in FreeFem++ syntax as

int2d(Th)((A*em(u))'*em(v)) - int2d(Th)(b'*v)

Stress and displacement fields

Convergence studies

For the convergence between analytical and numerical solution we vary the number of nodes by increasing the grid size in both $x$- and $y$- directions simultaneously by powers of two from $2^2$ (16 nodes all together) to $2^7$ (16384 nodes all together). The $L^2$ error norm is used to measure the "difference" between solutions. Since the displacements are the variable that obtain from FreeFem++, we use the displacement magnitude $|\boldsymbol{u}| = \sqrt{u_x^2+u_y^2}$ to define our $L^2$-error norm. The exact equation we have used is \begin{equation} L^2\text{-norm} = \sqrt{\frac{\int_\Omega (|\boldsymbol{u_{\text{numerical}}}|-|\boldsymbol{u_{\text{analytical}}}|)^2d\Omega}{\int_\Omega|\boldsymbol{u_{\text{analytical}}}|^2d\Omega}}. \end{equation} Results are shown in Figure 1.

Contact between parallel cylinders

References

• Theory of matrix structural analysis