2D Elasticity problem

Keywords: PDE:Elasticity Method:FEM FE:Lagrange Solver:Direct

The elasticity problem illustrates how to use vector unknown in XLiFE++:

\[\begin{split}\DeclareMathOperator{\sigmaop}{\sigma} \DeclareMathOperator{\divop}{div} \begin{cases} -\divop (\sigmaop \mathbf{u}) -\omega^2 \mathbf{u} = \mathbf{f} & \mathrm{in\;}\Omega \\ \sigmaop \mathbf{u} * \mathbf{n}=0 & \mathrm{\;on\;}\partial \Omega \end{cases}\end{split}\]

For homogeneous isotropic material:

\[\DeclareMathOperator{\epsilonop}{\varepsilon} \DeclareMathOperator{\sigmaop}{\sigma} \DeclareMathOperator{\divop}{div} \sigmaop \mathbf{u} = \lambda (\divop \mathbf{u}) \mathbb{I} + 2 \mu \epsilonop \mathbf{u}, \quad \epsilonop_{ij} \mathbf{u} = \partial_i u_j.\]

The variational formulation in \(V=(H^1(\Omega))^3\) is to find:

\[\DeclareMathOperator{\epsilonop}{\varepsilon} \DeclareMathOperator{\divop}{div} \mathbf{u} \in V \mid \forall \mathbf{v} \in V, \lambda \int_\Omega \divop \mathbf{u} \divop \mathbf{\bar{v}} + 2 \mu \int_\Omega \epsilonop \mathbf{u} : \epsilonop \mathbf{\bar{v}} - \omega^2 \int_\Omega \mathbf{u} \cdot \mathbf{\bar{v}} = \int_{\Omega} \mathbf{f} \cdot \mathbf{\bar{v}}\]

This is implemented as follows:

#include "xlife++.h"
using namespace xlifepp;

// data function
Vector<Real> f(const Point& P, Parameters& pa = defaultParameters)
{  Vector<Real> F(2, 0.); F(2)=-0.005; return F;}

int main(int argc, char** argv)
{
  init(argc, argv, _lang=en); // mandatory initialization of xlifepp

  // mesh rectangle
  Rectangle rect(_center=Point(0., 0.), _xlength=20, _ylength=2, _nnodes=Numbers(50, 5), _domain_name="Omega", _side_names=Strings("", "", "", "Gamma"));
  Mesh mesh2d(rect, _shape=triangle, _generator=gmsh);
  Domain omega=mesh2d.domain("Omega"), Gamma=mesh2d.domain("Gamma");
  // create P1 Lagrange interpolation
  Space V(_domain=omega, _interpolation=P1, _name="V");
  Unknown u(V, _name="u", _dim=2);  TestFunction v(u, _name="v");
  // create bilinear form and linear form
  Real lambda=167.06, mu=56,67, omg2=0, rho=7.86;
  BilinearForm auv = lambda*intg(omega, div(u)*div(v)) + 2*mu*intg(omega, epsilon(u) % epsilon(v)) - omg2*intg(omega, u|v);
  LinearForm fv=intg(omega, f|v);
  EssentialConditions ecs= (u|Gamma=0);
  TermMatrix A(auv, ecs, _name="A");
  TermVector B(fv, _name="B");
  //solve linear system AX=B using direct method
  TermVector U=directSolve(A, B);
  thePrintStream<<U;
  saveToFile("U", U, _format=vtu);

  // create the deformation of the mesh
  for (number_t i=0;i<mesh2d.nbOfNodes();i++)
    mesh2d.nodes[i] += U.evaluate(mesh2d.nodes[i]).value<std::vector<Real> >();
  saveToFile("Ud", mesh2d, _format=msh);

  return 0;
}

Fig. 31 Displacement and modulus of the solution of the elasticity 2D problem