Difference between revisions of "Linear elasticity"

Latest revision as of 18:03, 4 December 2022

Go back to Examples.

On this page we showcase some basic examples from linear elasticity. To read more about the governing equations, refer to the Solid Mechanics page and examples therein, which are considered in more detail.

All examples here will be the solutions of the Cauchy-Navier equation $$ (\lambda + \mu) \nabla(\nabla \cdot \vec{u}) + \mu \nabla^2 \vec{u} = 0. $$

Cantilever beam

Consider a beam of dimensions $L \times D$ having a narrow rectangular cross section. The beam occupies a region of $[0, L] \times [-D/2, D/2]$. The beam is bent by a force $P$ applied at the end $x = 0$ and the other end of the beam is fixed at $x = L$, as illustrated below.

The stresses in such a beam are given as: \begin{equation} \sigma_{xx} = -\frac{Pxy}{I}, \sigma_{yy} = 0, \sigma_{xy} = -\frac{P}{2I}\left(\frac{D^2}{4} - y^2 \right), \label{eq:sxy} \end{equation} where $I = D^3/12$ is the moment of inertia.

The exact solution in terms of the displacements in $x$ and $y$ direction is \begin{align}\label{eq:beam_a1} u_x(x,y) = u(x, y) &= -\frac{Py}{6EI}\left(3(x^2-L^2) -(2+\nu)y^2 + 6 (1+\nu) \frac{D^2}{4}\right) \\ \label{eq:beam_a2} u_y(x,y) = v(x, y) &= \frac{P}{6EI}\left(3\nu x y^2 + x^3 - 3L^2 x + 2L^3\right) \end{align} where $E$ is Young's modulus and $\nu$ is the Poisson ratio. More details can be found in File:Cantilever beam.nb.

The solution of the cantilever beam problem is illustrated in the cantilever_beam.cpp file. The problem is defined by the following pyhsical parameters:

    const double E = 72.1e9;
    const double nu = 0.33;
    const double P = 1000;
    const double D = 5;
    const double L = 30;

and the problem itself is specified as follows:

    for (int i : domain.interior()) {
        (lam+mu)*op.graddiv(i) + mu*op.lap(i) = 0.0;
    }
    for (int i : domain.types() == RIGHT) {
        double y = domain.pos(i, 1);
        op.value(i) = {(P*y*(3*D*D*(1+nu) - 4*(2+nu)*y*y)) / (24.*E*I), -(L*nu*P*y*y) / (2.*E*I)};
    }
    for (int i : domain.types() == LEFT) {
        double y = domain.pos(i, 1);
        op.traction(i, lam, mu, {-1, 0}) = {0, -P*(D*D - 4*y*y) / (8.*I)};
    }
    for (int i : domain.types() == TOP) {
        op.traction(i, lam, mu, {0, 1}) = 0.0;
    }
    for (int i : domain.types() == BOTTOM) {
        op.traction(i, lam, mu, {0, -1}) = 0.0;
    }

The plot of the numerical solution produced by the accompanying cantilever_beam.m is shown below.

Point contact 2D

Next we consider the point contact in 2D with the known analytical solution $$\vec{u}(x, y) = \left[-\frac{P}{4 \pi \mu } \left(\frac{2 \mu}{\lambda +\mu } \text{atan2}(y, x)+\frac{2 x y}{x^2+y^2}\right),-\frac{P}{4 \pi \mu } \left(\frac{y^2-x^2}{x^2+y^2}-\frac{(\lambda +2 \mu ) \log \left(x^2+y^2\right)}{\lambda +\mu}\right)\right]$$

More details can be found in File:Point contact2d.nb and on the Point contact page.

The problem is solved numerically on the domain $\Omega = [-1, -1] \times [1, -0.1]$ with Dirichlet boundary conditions on a scattered node set. See the point_contact2d.cpp file for reference.

    for (int i : domain.interior()) {
        (lam+mu)*op.graddiv(i) + mu*op.lap(i) = 0.0;
    }
    for (int i : domain.boundary()) {
        op.value(i) = analytical(domain.pos(i));
    }

The obtained solution plotted with point_contact2d.m is shown below.

Point contact 3D

We also present a 3-D example of linear elasticity point contact solved on a scattered node set with variable density. Due to singularity at the origin, we choose the domain $\Omega = [-a, -b]^3$ with $a=1$ and $b=0.1$. The problem details are available in File:Point contact3d.nb.

See the point_contact2d.cpp file for reference.

The assembly of the matrix looks exactly the same as is 2-D.

    for (int i : domain.interior()) {
        (lam+mu)*op.graddiv(i) + mu*op.lap(i) = 0.0;
    }
    for (int i : domain.boundary()) {
        op.value(i) = analytical(domain.pos(i));
    }

The computation of stresses also follows mathematical definition $\sigma = \lambda \operatorname{tr}(\varepsilon) I + 2\mu \varepsilon, \varepsilon = \frac12 (\nabla \vec{u} + (\nabla \vec{u}^\mathrm{T}))$.

    VectorField<double, 6> stress(N);
    for (int i = 0; i < N; ++i) {
        Eigen::Matrix3d grad = eop.grad(u, i);
        Eigen::Matrix3d eps = 0.5*(grad + grad.transpose());
        Eigen::Matrix3d s = lam * eps.trace() * Eigen::Matrix3d::Identity(3, 3) + 2*mu*eps;
        stress[i][0] = s(0, 0);
        stress[i][1] = s(1, 1);
        stress[i][2] = s(2, 2);
        stress[i][3] = s(0, 1);
        stress[i][4] = s(0, 2);
        stress[i][5] = s(1, 2);
    }

The obtained solution plotted with point_contact3d.m is shown below.

Go back to Examples.