Wave equation

2022-12-04T19:08:48Z

E6WikiAdmin:

{{Box-round|title= Related papers |

[https://e6.ijs.si/ParallelAndDistributedSystems/publications/62417923.pdf J. Močnik Berljavac, P. Mishra, J. Slak, G. Kosec; RBF-FD analysis of 2D time-domain acoustic wave propagation in heterogeneous media, Computers & Geosciences, 2021 [DOI: 10.1016/j.cageo.2021.104796]]
}}

__TOC__

Go back to [[Medusa#Examples|Examples]].

== 2D wave equation with Dirichlet boundary conditions ==
We are solving a problem analogous to a clamped circular membrane being forced to oscillate by fixing its centre to a rigid oscillating rod. The domain of calculation is for this reason an annulus bounded by the point of contact with the circular source in the centre and the clamped edge of the circular membrane.

[[File:wave_2d_domain_geo.jpg|300px]]

Consider the time dependent solution to 2D wave equation on annulus shaped domain

<math>
\begin{align*}
\frac{ \partial^2 u}{\partial t^2} &= c^2 \nabla^2 u &&\text{in } \Omega, \\
u(r,t=0)&=0 &&\text{in } \Omega, \\
u &= 0 &&\text{on } \partial \Omega_O,\\
u &= f(t) &&\text{on } \partial \Omega_I,
\end{align*}
</math>

where $\partial\Omega_I$ denotes the inner and $\partial\Omega_O$ the outer boundary of the domain. Through the boundary condition on the inner boundary the source is introduced to the problem as a function of time

<math>
\begin{align*}
f(t)= u_o \sin \omega_o t.
\end{align*}
</math>

First the domain is constructed by subtracting a smaller circle domain from a larger one. Boundaries of the domain are populated in the same step.

<syntaxhighlight lang="c++" line>
// // identifier to be added to nodes on the inner boundary
int CENTRE = -10;

// build circle domain
BallShape<Vec2d> domain({0, 0}, outer_radius);
auto discretization = domain.discretizeBoundaryWithStep(dx);

// build source domain
BallShape<Vec2d> empty({0, 0}, inner_radius);
auto discretization_empty = empty.discretizeBoundaryWithStep(dx, CENTRE);

// substract the source domain
discretization -= discretization_empty;
</syntaxhighlight>

Next the domain is populated with nodes in accordance with the desired density function. Once the domain is constructed and discretized, the support (neighborhood) of $n$ nodes is found for each node.

<syntaxhighlight lang="c++" line>

GeneralFill<Vec2d> fill;
fill.seed(fill_seed);
discretization.fill(fill, fill_density);

// find support
FindClosest find_support(n);
discretization.findSupport(find_support);
</syntaxhighlight>

The density function "fill_density" takes as argument the coordinates inside the domain and returns the desired distance between neighboring nodes for that location.
Now the domain preparation is complete, as domain is constructed and populated with nodes with support relations established. We continue by constructing the approximation engine.
<syntaxhighlight lang="c++" line>
WLS<Monomials<Vec2d>, GaussianWeight<Vec2d>, ScaleToClosest> approx(m - 1 , sigma);
</syntaxhighlight>

As we are solving a time propagation problem implicitly, we will have to solve a system of equations for every time step. This is equivalent to solving a matrix equation. For this reason a space matrix is constructed. Next the wave equation is turned into code. Medusa encodes this information inside the just constructed matrix $M$. $E1$ and $E0$ represent the known solutions on previous two time steps.

<syntaxhighlight lang="c++" line>
SparseMatrix<double> M(domain_size, domain_size);
// Set equation on interior
for (int i : interior) {
-(v*v*dt*dt) * op.lap(i) + op.value(i)= 2*E1(i)-E0(i);
}
</syntaxhighlight>
Same is done also for the boundary nodes. Since dirichlet boundary conditions are employed the values of the Right Hand Side ($rhs$) of matrix equation are just set to prescribed values ($0$ on the outer boundary and $f(t)$ on inner boundary).
<syntaxhighlight lang="c++" line>
// Set boundary conditions
for (int i : boundary) {
if (discretization.types()[i] == CENTRE) {
op.value(i) = A * std::sin(omega*dt*2); // step
} else {
op.value(i) = 0.0;
}
}
</syntaxhighlight>
At each time step the matrix equation is solved using an Eigen linear equation solver, in this case, the BiCGStab algorithm. $E0$ and $E1$ are updated and the source is set to new value.

<syntaxhighlight lang="c++" line>
// Solve matrix system
VectorXd E2 = solver.solve(rhs);

// Update previous states
E0 = E1;
E1 = E2;

// Update rhs on interior
for (int i : interior) {
rhs(i) = 2*E1(i)-E0(i);
}

// Update rhs on boundary
for (int i : (discretization.types() == CENTRE)) {
rhs(i) = A * std::sin(omega*dt*tt); // step
}
</syntaxhighlight>
The animation of solution is presented below. The complete exampe is avalible at [https://gitlab.com/e62Lab/medusa/tree/dev/examples/wave_equation wave_equation_2D.cpp].

[[File:EM_wave_2d_dt_e-6_dens40_omega128_v8 2X speed.mp4|400px]]

== 2D time-domain acoustic wave propagation in heterogeneous media ==

Go back to [[Medusa#Examples|Examples]].

Wave equation

2022-12-04T19:06:57Z

E6WikiAdmin:

{{Box-round|title= Related papers |

[https://e6.ijs.si/ParallelAndDistributedSystems/publications/62417923.pdf J. Močnik Berljavac, P. Mishra, J. Slak, G. Kosec; RBF-FD analysis of 2D time-domain acoustic wave propagation in heterogeneous media, Computers & Geosciences, 2021 [DOI: 10.1016/j.cageo.2021.104796]]
}}

__TOC__

Go back to [[Medusa#Examples|Examples]].

== 2D wave equation with Dirichlet boundary conditions ==
We are solving a problem analogous to a clamped circular membrane being forced to oscillate by fixing its centre to a rigid oscillating rod. The domain of calculation is for this reason an annulus bounded by the point of contact with the circular source in the centre and the clamped edge of the circular membrane.

[[File:wave_2d_domain_geo.jpg|300px]]

Consider the time dependent solution to 2D wave equation on annulus shaped domain

<math>
\begin{align*}
\frac{ \partial^2 u}{\partial t^2} &= c^2 \nabla^2 u &&\text{in } \Omega, \\
u(r,t=0)&=0 &&\text{in } \Omega, \\
u &= 0 &&\text{on } \partial \Omega_O,\\
u &= f(t) &&\text{on } \partial \Omega_I,
\end{align*}
</math>

where $\partial\Omega_I$ denotes the inner and $\partial\Omega_O$ the outer boundary of the domain. Through the boundary condition on the inner boundary the source is introduced to the problem as a function of time

<math>
\begin{align*}
f(t)= u_o \sin \omega_o t.
\end{align*}
</math>

First the domain is constructed by subtracting a smaller circle domain from a larger one. Boundaries of the domain are populated in the same step.

<syntaxhighlight lang="c++" line>
// // identifier to be added to nodes on the inner boundary
int CENTRE = -10;

// build circle domain
BallShape<Vec2d> domain({0, 0}, outer_radius);
auto discretization = domain.discretizeBoundaryWithStep(dx);

// build source domain
BallShape<Vec2d> empty({0, 0}, inner_radius);
auto discretization_empty = empty.discretizeBoundaryWithStep(dx, CENTRE);

// substract the source domain
discretization -= discretization_empty;
</syntaxhighlight>

Next the domain is populated with nodes in accordance with the desired density function. Once the domain is constructed and discretized, the support (neighborhood) of $n$ nodes is found for each node.

<syntaxhighlight lang="c++" line>

GeneralFill<Vec2d> fill;
fill.seed(fill_seed);
discretization.fill(fill, fill_density);

// find support
FindClosest find_support(n);
discretization.findSupport(find_support);
</syntaxhighlight>

The density function "fill_density" takes as argument the coordinates inside the domain and returns the desired distance between neighboring nodes for that location.
Now the domain preparation is complete, as domain is constructed and populated with nodes with support relations established. We continue by constructing the approximation engine.
<syntaxhighlight lang="c++" line>
WLS<Monomials<Vec2d>, GaussianWeight<Vec2d>, ScaleToClosest> approx(m - 1 , sigma);
</syntaxhighlight>

As we are solving a time propagation problem implicitly, we will have to solve a system of equations for every time step. This is equivalent to solving a matrix equation. For this reason a space matrix is constructed. Next the wave equation is turned into code. Medusa encodes this information inside the just constructed matrix $M$. $E1$ and $E0$ represent the known solutions on previous two time steps.

<syntaxhighlight lang="c++" line>
SparseMatrix<double> M(domain_size, domain_size);
// Set equation on interior
for (int i : interior) {
-(v*v*dt*dt) * op.lap(i) + op.value(i)= 2*E1(i)-E0(i);
}
</syntaxhighlight>
Same is done also for the boundary nodes. Since dirichlet boundary conditions are employed the values of the Right Hand Side ($rhs$) of matrix equation are just set to prescribed values ($0$ on the outer boundary and $f(t)$ on inner boundary).
<syntaxhighlight lang="c++" line>
// Set boundary conditions
for (int i : boundary) {
if (discretization.types()[i] == CENTRE) {
op.value(i) = A * std::sin(omega*dt*2); // step
} else {
op.value(i) = 0.0;
}
}
</syntaxhighlight>
At each time step the matrix equation is solved using an Eigen linear equation solver, in this case, the BiCGStab algorithm. $E0$ and $E1$ are updated and the source is set to new value.

<syntaxhighlight lang="c++" line>
// Solve matrix system
VectorXd E2 = solver.solve(rhs);

// Update previous states
E0 = E1;
E1 = E2;

// Update rhs on interior
for (int i : interior) {
rhs(i) = 2*E1(i)-E0(i);
}

// Update rhs on boundary
for (int i : (discretization.types() == CENTRE)) {
rhs(i) = A * std::sin(omega*dt*tt); // step
}
</syntaxhighlight>
The animation of solution is presented below. The complete exampe is avalible at [https://gitlab.com/e62Lab/medusa/tree/dev/examples/wave_equation wave_equation_2D.cpp].

[[File:EM_wave_2d_dt_e-6_dens40_omega128_v8 2X speed.mp4|400px]]

Go back to [[Medusa#Examples|Examples]].

Wave equation

2022-12-04T19:04:48Z

E6WikiAdmin:

{{Box-round|title= Related papers |

[https://e6.ijs.si/ParallelAndDistributedSystems/publications/62417923.pdf J. Močnik Berljavac, P. Mishra, J. Slak, G. Kosec; RBF-FD analysis of 2D time-domain acoustic wave propagation in heterogeneous media, Computers & Geosciences, 2021 [DOI: 10.1016/j.cageo.2021.104796]]
}}

__TOC__

Go back to [[Medusa#Examples|Examples]].

== 2D wave equation with Dirichlet boundary conditions ==
We are solving a problem analogous to a clamped circular membrane being forced to oscillate by fixing its centre to a rigid oscillating rod. The domain of calculation is for this reason an annulus bounded by the point of contact with the circular source in the centre and the clamped edge of the circular membrane.

[[File:wave_2d_domain_geo.jpg|300px]]

Consider the time dependent solution to 2D wave equation on annulus shaped domain

<math>
\begin{align*}
\frac{ \partial^2 u}{\partial t^2} &= c^2 \nabla^2 u &&\text{in } \Omega, \\
u(r,t=0)&=0 &&\text{in } \Omega, \\
u &= 0 &&\text{on } \partial \Omega_O,\\
u &= f(t) &&\text{on } \partial \Omega_I,
\end{align*}
</math>

where $\partial\Omega_I$ denotes the inner and $\partial\Omega_O$ the outer boundary of the domain. Through the boundary condition on the inner boundary the source is introduced to the problem as a function of time

<math>
\begin{align*}
f(t)= u_o \sin \omega_o t.
\end{align*}
</math>

First the domain is constructed by subtracting a smaller circle domain from a larger one. Boundaries of the domain are populated in the same step.

<syntaxhighlight lang="c++" line>
// // identifier to be added to nodes on the inner boundary
int CENTRE = -10;

// build circle domain
BallShape<Vec2d> domain({0, 0}, outer_radius);
auto discretization = domain.discretizeBoundaryWithStep(dx);

// build source domain
BallShape<Vec2d> empty({0, 0}, inner_radius);
auto discretization_empty = empty.discretizeBoundaryWithStep(dx, CENTRE);

// substract the source domain
discretization -= discretization_empty;
</syntaxhighlight>

Next the domain is populated with nodes in accordance with the desired density function. Once the domain is constructed and discretized, the support (neighborhood) of $n$ nodes is found for each node.

<syntaxhighlight lang="c++" line>

GeneralFill<Vec2d> fill;
fill.seed(fill_seed);
discretization.fill(fill, fill_density);

// find support
FindClosest find_support(n);
discretization.findSupport(find_support);
</syntaxhighlight>

The density function "fill_density" takes as argument the coordinates inside the domain and returns the desired distance between neighboring nodes for that location.
Now the domain preparation is complete, as domain is constructed and populated with nodes with support relations established. We continue by constructing the approximation engine.
<syntaxhighlight lang="c++" line>
WLS<Monomials<Vec2d>, GaussianWeight<Vec2d>, ScaleToClosest> approx(m - 1 , sigma);
</syntaxhighlight>

As we are solving a time propagation problem implicitly, we will have to solve a system of equations for every time step. This is equivalent to solving a matrix equation. For this reason a space matrix is constructed. Next the wave equation is turned into code. Medusa encodes this information inside the just constructed matrix $M$. $E1$ and $E0$ represent the known solutions on previous two time steps.

<syntaxhighlight lang="c++" line>
SparseMatrix<double> M(domain_size, domain_size);
// Set equation on interior
for (int i : interior) {
-(v*v*dt*dt) * op.lap(i) + op.value(i)= 2*E1(i)-E0(i);
}
</syntaxhighlight>
Same is done also for the boundary nodes. Since dirichlet boundary conditions are employed the values of the Right Hand Side ($rhs$) of matrix equation are just set to prescribed values ($0$ on the outer boundary and $f(t)$ on inner boundary).
<syntaxhighlight lang="c++" line>
// Set boundary conditions
for (int i : boundary) {
if (discretization.types()[i] == CENTRE) {
op.value(i) = A * std::sin(omega*dt*2); // step
} else {
op.value(i) = 0.0;
}
}
</syntaxhighlight>
At each time step the matrix equation is solved using an Eigen linear equation solver, in this case, the BiCGStab algorithm. $E0$ and $E1$ are updated and the source is set to new value.

<syntaxhighlight lang="c++" line>
// Solve matrix system
VectorXd E2 = solver.solve(rhs);

// Update previous states
E0 = E1;
E1 = E2;

// Update rhs on interior
for (int i : interior) {
rhs(i) = 2*E1(i)-E0(i);
}

// Update rhs on boundary
for (int i : (discretization.types() == CENTRE)) {
rhs(i) = A * std::sin(omega*dt*tt); // step
}
</syntaxhighlight>
The animation of solution is presented below. The complete exampe is avalible at [https://gitlab.com/e62Lab/medusa/tree/dev/examples/wave_equation wave_equation_2D.cpp]. This example is covered in more detail in linked document [[:File:Meshless solution to wave equation.pdf]].

[[File:EM_wave_2d_dt_e-6_dens40_omega128_v8 2X speed.mp4|400px]]

Go back to [[Medusa#Examples|Examples]].

Solid Mechanics

2022-12-04T17:33:16Z

E6WikiAdmin: /* References */

= Numerical examples =

The following pages describe several case studies in the field of elastostatics and elastodynamics. The subpages include the problem description (equations and boundary conditions), analytical solutions if obtainable and numerical results for given parameter values.

*[[Point contact]]
*[[Cantilever beam]]
*[[Hertzian contact]]
*[[FWO case | Fretting fatigue contact]]

=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 [https://en.wikipedia.org/wiki/Linear_elasticity#Stress_formulation '''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 [https://en.wikipedia.org/wiki/Lam%C3%A9_parameters 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 stress with Lamé constants ====

The plane stress stiffness tensor (\ref{eq:planestressmatrix}) may also be expressed in terms of Lamé constants $\lambda$ and $\mu$ by substituting $E$ and $\nu$ with
\[E = \frac{\mu(3\lambda+2\mu)}{\lambda+\mu}\]
\[\nu = \frac{\lambda}{2(\lambda+\mu)}\]

With these substitutions the Navier equations can be written in their component form as
\begin{align}
\rho \frac{\partial^2 u_x}{\partial t^2} &= \mu \left( \frac{\partial^2 u_x}{\partial x^2} + \frac{\partial^2 u_x}{\partial y^2} \right) + \left(\frac{2 \lambda \mu}{\lambda + 2 \mu}+\mu\right) \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} &= \mu \left( \frac{\partial^2 u_y}{\partial x^2} + \frac{\partial^2 u_y}{\partial y^2} \right) + \left(\frac{2 \lambda \mu}{\lambda + 2 \mu}+\mu\right) \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 & \mu^*
\end{bmatrix} \qquad \text{(Plane strain)}
\end{equation}
(*Note that this matrix requires the strain tensor to be defined using the engineering shear strain $\gamma_{xy} = 2\varepsilon_{xy}$. If we define the strain tensor as $\b{\varepsilon} = \{\varepsilon_{xx}\;\varepsilon_{yy}\;\varepsilon_{xy}\}^T$ then matrix component $C_{33}$ must be $2\mu$.)

The Navier equations in vector form can be found in (\ref{eq:navier_mu_lambda}). In component notation the equations are:
\begin{align}
\rho \frac{\partial^2 u_x}{\partial t^2} &= \mu \left( \frac{\partial^2 u_x}{\partial x^2} + \frac{\partial^2 u_x}{\partial y^2} \right) + (\lambda+\mu) \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} &= \mu \left( \frac{\partial^2 u_y}{\partial x^2} + \frac{\partial^2 u_y}{\partial y^2} \right) + (\lambda+\mu) \frac{\partial}{\partial y}\left( \frac{\partial u_x}{\partial x} + \frac{\partial u_y}{\partial y} \right)
\end{align}

=== 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 [http://www.colorado.edu/engineering/cas/courses.d/IFEM.d/IFEM.Ch14.d/IFEM.Ch14.pdf at this link (page 13)].

To go from plane stres'''s''' ($s$) to plane strai'''n''' ($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}

<div class="toccolours mw-collapsible mw-collapsed">
'''Substitution from plane stress to plane strain'''
<div class="mw-collapsible-content">Note: this derivation is shown just to confirm the above formulas. In a numerical program, we keep the stress-strain matrix and just use the above formulas (\ref{eq:ston1}) and (\ref{eq:ston2}) to update the values of $E$ and $\nu$. Also note we have omitted the indexes ($s$) and ($n$) in the following derivations.

We start with the plane stress matrix with the inserted fictitious values $E^*$ and $\nu^*$
\[
\frac{E^*}{1-{\nu^*}^2}
\begin{bmatrix}
1 & \nu^* & 0 \\
\nu^* & 1 & 0 \\
0 & 0 & \frac{1}{2}(1-\nu^*)
\end{bmatrix}
\]
and make the above substitutions (\ref{eq:ston1}) and (\ref{eq:ston2}) leading to:
\[
\frac{E}{\left(1-\nu^2\right)\left(1-\left(\frac{\nu}{1-\nu}\right)^2\right)}
\begin{bmatrix}
1 & \frac{\nu}{1-\nu} & 0 \\
\frac{\nu}{1-\nu} & 1 & 0 \\
0 & 0 & \frac{1}{2}(1-\frac{\nu}{1-\nu})
\end{bmatrix}.
\]
We can then use the rule to convert sums of squares into products as well as bring the factor $1/(1-v)$ out of the matrix:
\[
\frac{E}{\left(1-\nu\right)\left(1+\nu\right)\left(1-\frac{\nu}{1-\nu}\right)\left(1+\frac{\nu}{1-\nu}\right)\left(1-\nu\right)}
\begin{bmatrix}
1-\nu & \nu & 0 \\
\nu & 1-\nu & 0 \\
0 & 0 & \frac{1}{2}(1-2\nu)
\end{bmatrix}.
\]
By joining some of the factors to a common denominator and rearranging leads to
\[
\frac{E\left(1-\nu\right)\left(1-\nu\right)}{\left(1-\nu\right)\left(1+\nu\right)\left(1-2\nu\right)1\left(1-\nu\right)}
\begin{bmatrix}
1-\nu & \nu & 0 \\
\nu & 1-\nu & 0 \\
0 & 0 & \frac{1}{2}(1-2\nu)
\end{bmatrix}.
\]
The final step is canceling the factors that occur in both the numerator and denominator
\[
\frac{E}{\left(1+\nu\right)\left(1-2\nu\right)}
\begin{bmatrix}
1-\nu & \nu & 0 \\
\nu & 1-\nu & 0 \\
0 & 0 & \frac{1}{2}(1-2\nu)
\end{bmatrix}
\]
which leads exactly to the relationship for plane strain given in (\ref{eq:planestrain}).
</div>
</div>

In the opposite case we want to go from plane strai'''n''' ($n$) to plane stres'''s''' ($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}

<div class="toccolours mw-collapsible mw-collapsed">
'''Substitution from plane strain to plane stress'''
<div class="mw-collapsible-content">
Note that we have again omitted the indexes ($s$) and ($n$) since our wish is expressed by the bold title.

We start with the plane strain matrix with the inserted fictitious values $E^*$ and $\nu^*$
\[
\frac{E^*}{\left(1+\nu^*\right)\left(1-2\nu^*\right)}
\begin{bmatrix}
1-\nu^* & \nu^* & 0 \\
\nu^* & 1-\nu^* & 0 \\
0 & 0 & \frac{1}{2}(1-2\nu^*)
\end{bmatrix}
\]
and make the above substitutions (\ref{eq:ntos1}) and (\ref{eq:ntos2}) leading to:
\[
\frac{E(1+2\nu)}{\left(1+\nu\right)^2\left(1+\frac{\nu}{1+\nu}\right)\left(1-2\frac{\nu}{1+\nu}\right)}
\begin{bmatrix}
1-\frac{\nu}{1+\nu} & \frac{\nu}{1+\nu} & 0 \\
\frac{\nu}{1+\nu} & 1-\frac{\nu}{1+\nu} & 0 \\
0 & 0 & \frac{1}{2}(1-2\frac{\nu}{1+\nu})
\end{bmatrix}.
\]
Writing some of the sums with common denominators and rearranging leads to
\[
\frac{E(1+2\nu)(1+\nu)(1+\nu)}{\left(1+\nu\right)^2\left(1+2\nu\right)\left(1-\nu\right)}
\begin{bmatrix}
\frac{1}{1+\nu} & \frac{\nu}{1+\nu} & 0 \\
\frac{\nu}{1+\nu} & \frac{1}{1+\nu} & 0 \\
0 & 0 & \frac{1}{2}\left(\frac{1-\nu}{1+\nu}\right)
\end{bmatrix}.
\]
We can also bring the factor $1/(1+\nu)$ out from the matrix components. After canceling all factors that occur in both denominator and numerator we are left with
\[
\frac{E}{\left(1+\nu\right)\left(1-\nu\right)}
\begin{bmatrix}
1 & \nu & 0 \\
\nu & 1 & 0 \\
0 & 0 & \frac{1}{2}(1-\nu)
\end{bmatrix}.
\]
Rewriting the product of sums as a sum of squares gives us the matrix for plane stress in (\ref{eq:planestress}).
</div>
</div>

=== Principal stresses ===
This are the stresses when we rotate our coordinate system so that there are no shear components.
The principal stresses for plane stress conditions are given as eigenvalues of $\sigma$ with
\[
\sigma_{1,2} = \frac12 (s_{xx} + s_{zz}) \pm \sqrt{\frac{1}{4}(s_{xx}-s_{zz})^2 + s_{xz}^2}.
\]
The principal shear stress or principal stress difference is the stress where only shear stresses are present. In case of plane stress it is calculated as
\[
\tau_1 = \sqrt{\frac{1}{4}(s_{xx}-s_{zz})^2 + s_{xz}^2} = \frac{|\sigma_1 - \sigma_2|}{2}.
\]

When the stress tensor is z 3x3 matrix, the formulas for eigenvalues are more complicated.

=== Von Mises stress ===

A more detailed explanation of the von Mises stress can be found on [https://en.wikipedia.org/wiki/Von_Mises_yield_criterion Wikipedia].

In materials science and engineering the von Mises yield criterion can be also formulated in terms of the von Mises stress or equivalent tensile stress, $\sigma_\mathrm{v}$. This is a scalar value of stress that can be computed from the Cauchy stress tensor $\boldsymbol{\sigma}$. In this case, a material is said to start yielding when the von Mises stress reaches a value known as yield strength, $\sigma_\mathrm{y}$.

Mathematically the von Mises yield criterion can be expressed as:
\[
\sigma _\mathrm{v}=\sigma_\mathrm{y}={\sqrt {3J_{2}}},
\]
where $J_2$ is the [https://en.wikipedia.org/wiki/Cauchy_stress_tensor#Stress_deviator_tensor second deviatoric stress invariant].

In the case of general '''plane stress''' ($\sigma_{33}=0$,$\sigma_{23}=0$, $\sigma_{31}=0$) the von Mises stress can be calculated as:
\[
\sigma_\mathrm{v} = \sqrt{\sigma_{11}^2 - \sigma_{11}\sigma_{22}+\sigma_{22}^2 + 3\sigma_{12}^2}.
\]

= computational solid mechanics =

The goal in computational solid mechanics is to solve the equation of motion for a solid body
\[\rho \frac{\partial^2 \b{u}}{\partial t^2} = \nabla\cdot\b{\sigma} + \b{F}\]

on the given domain $\Gamma = \Gamma_N \cup \Gamma_D$ with Dirichlet (essential) boundary conditions
\[ \b{u} = \bar{\b{u}} \quad \text{ on } \Gamma_D,\]
and traction (natural) boundary conditions
\[ \b{\sigma}\cdot\b{n} = \bar{\b{t}},\quad \text{ on } \Gamma_N,\]
where the quantities with the bar sign $\bar{}$ indicate prescribed displacements and surface tractions, respectively. Another type of boundary condition are mixed (Robyn) type boundary conditions. In this case the displacement might be given in one direction, and the traction will be given for the other direction.
In case the right-hand side value is zero, the boundary conditions are called homogeneous. Homogeneous Dirichlet boundary conditions mean the the body is fully fixed (constrained), while zero traction boundary conditions allow the surface to move freely.

Using Hooke's law and the strain-displacement relation, the equation of motion can be transformed into the Navier-Cauchy equations
\[\rho \frac{\partial^2 \b{u}}{\partial t^2} = \mu\nabla^2\b{u} + (\lambda+\mu)\nabla(\nabla\cdot\b{u}) + \b{F}.\]
The solution of this equation with appropriate BC's will give us the desired displacement values across the domain.

To obtain the steady state solution we have two possible strategies:
# Add a (linear) dampening term to the Navier equation, and simulate the dynamic behaviour until the motion ceases: \[\rho \frac{\partial^2 \b{u}}{\partial t^2} + \eta_C \frac{\partial \b{u}}{\partial t} = \mu\nabla^2\b{u} + (\lambda+\mu)\nabla(\nabla\cdot\b{u}) + \b{F}.\]
# Directly solve for steady state: \[\mu\nabla^2\b{u} + (\lambda+\mu)\nabla(\nabla\cdot\b{u}) + \b{F} = 0.\]

Displacement boundary conditions are very simple to implement with both approaches. The desired values are simply prescribed for dynamic simulation with explicit methods, In case of implicit methods where we need to solve a system of equations the prescribed displacements are placed into the right-hand side of the system, and the value 1 is placed in the diagonal of the matrix (for a given node).

For the traction boundary condition

\[\begin{bmatrix}
\sigma_{xx} & \sigma_{xy} \\
\sigma_{yx} & \sigma_{yy}
\end{bmatrix} \cdot
\begin{bmatrix}
n_x \\
n_y
\end{bmatrix} =
\begin{bmatrix}
\bar{t}_x \\
\bar{t}_y
\end{bmatrix}\]

we first use Hooke's law and the strain-displacement relations to express the traction boundary in terms of displacements:

\[\begin{bmatrix}
(2\mu+\lambda)\partial_x u_x + \lambda \partial_y u_y & \mu(\partial_x u_y + \partial_y u_x) \\
\mu(\partial_x u_y + \partial_y u_x) & \lambda \partial_x u_x + (2\mu + \lambda) \partial_y u_y
\end{bmatrix} \cdot
\begin{bmatrix}
n_x \\
n_y
\end{bmatrix} =
\begin{bmatrix}
\bar{t}_x \\
\bar{t}_y
\end{bmatrix}\]

=References=
* Slak J., Kosec G.; Refined meshless local strong form solution of Cauchy-Navier equation on an irregular domain. Engineering analysis with boundary elements. 2018;11 ; [http://comms.ijs.si/~gkosec/data/papers/31107623.pdf manuscript]
* Theory of matrix structural analysis

'''This list might be incomplete, check also reference on the Main Page.'''