Solid Mechanics

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. For a general three dimensional solid object these equations are:

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

where the number in brackets indicates the number of equations. We see there are 15 equations for 15 unknown variables (6 strains, 6 stresses and 3 displacements). In two dimensions this simplifies to 8 equations with 2 displacements, 3 stresses, and 3 strains.

Strain-displacement equations

Under the action of applied forces, a point in the solid originally at $(x,y,z)$ moves to position $(X,Y,Z)$. This movement can be described completely by the displacement vector \begin{equation} \boldsymbol{u(\boldsymbol{x})} = \{u_x(x,y,z) \quad u_y(x,y,z) \quad u_z(x,y,z)\}, \end{equation} the brackets $\{\}$ indicating that this is a column vector.

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_{xy} = -\frac{2P}{\pi} \frac{y^3}{\left(x^2+y^2\right)^2}, \end{equation} \begin{equation} \sigma_{yy} = -\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,v)$ can be determined. The displacements are given by \begin{align} u &= -\frac{P}{4\pi\mu}\left((\kappa-1)\theta - \frac{2xy}{r^2}\right), \label{eq:dispx}\\ v &= -\frac{P}{4\pi\mu}\left((\kappa+1)\log r - \frac{2y^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 as is defined as \[ \kappa = \begin{cases} 3 - 4\nu & \qquad \text{plane strain}, \\ \cfrac{3 - \nu}{1+\nu} & \qquad \text{plane stress}, \end{cases} \] and $\mu$ is the well known 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}).

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} \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".

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}