Beam deformation due to localized force distribution

Question

I would like to simulate the deformation of a beam which is fixed on one end and has some localized forces applied along the top.

I have been inspired by this Mathematica example: http://www.wolfram.com/mathematica/new-in-10/pdes-and-finite-elements/compute-a-plane-strain-deformation.html

First let's define the planeStrainOperator:

The planeStrainOperator represents the system of differential equations that needs to be solved to find the displacements at quilibrium. The standard continuum mechanics are derived from the Navier-Stokes equation.

The displacements

$u{}_{1}=u_{1}(x_{1},x_{2})$

$u{}_{2}=u_{2}(x_{1},x_{2})$

The deformations

$\epsilon{}_{11}=\partial u_{1}(x_{1},x_{2})/\partial x_{1}$

$\epsilon{}_{22}=\partial u_{2}(x_{1},x_{2})/\partial x_{2}$

$\epsilon{}_{12}=\epsilon{}_{21}=1/2(\partial u_{1}(x_{1},x_{2})/\partial x_{2}+\partial u_{2}(x_{1},x_{2})/\partial x_{1})$

The stresses on each plane of a 2D rectangle ($\nu$ is the Poisson's Ratio, E the Young's modul)

$\sigma{}_{11}=\frac{E}{(1+\nu)*(1-2*\nu)}*(\epsilon_{11}(1-\nu)+\nu*\epsilon_{22})$

$\sigma{}_{22}=\frac{E}{(1+\nu)*(1-2*\nu)}*(\epsilon_{22}(1-\nu)+\nu*\epsilon_{11})$

$\sigma{}_{12}=\sigma{}_{21}=\frac{E}{(1+\nu)}*\epsilon_{12}$

The quilibrium equations

$(1)\frac{\partial\sigma_{11}}{\partial x_{1}}+\frac{\partial\sigma_{12}}{\partial x_{2}}=0$

$(2)\frac{\partial\sigma_{21}}{\partial x_{1}}+\frac{\partial\sigma_{22}}{\partial x_{2}}=0$

In a vector form:

$(1)\nabla.\left(\begin{array}{c} \sigma_{11}\\ \sigma_{12} \end{array}\right)=0$

$(2)\nabla.\left(\begin{array}{c} \sigma_{21}\\ \sigma_{22} \end{array}\right)=0$

We now only need to insert the definitions given above to receive the planeStrainOperator.

$(1)\nabla.\left(\begin{array}{cc} \frac{(1-\nu)*E}{(1-2*\nu)*(1+\nu)} & 0\\ 0 & \frac{-E}{2*(1+\nu)} \end{array}\right).\nabla\left(u_{1}\right)+\nabla.\left(\begin{array}{cc} 0 & \frac{-\nu*E}{(1-2\nu)*(1+\nu)}\\ \frac{-E}{2*(1+\nu)} & 0 \end{array}\right).\nabla\left(u_{2}\right)=0$

$(2)\nabla.\left(\begin{array}{cc} 0 & \frac{-E}{2*(1+\nu)}\\ \frac{-\nu E}{(1-2*\nu)*(1+\nu)} & 0 \end{array}\right).\nabla\left(u_{1}\right)+\nabla.\left(\begin{array}{cc} \frac{-E}{2*(1+\nu)} & 0\\ 0 & \frac{-(1-\nu)*E}{(1-2*\nu)*(1+\nu)} \end{array}\right).\nabla\left(u_{2}\right)=0$

In the Mathematica example, Young's modulus is represented by a Y.

planeStrainOperator[Y_, ν_] := 
  {Inactive[Div][({{0, -((Y ν)/((1 - 2 ν) (1 + ν)))}, {-(Y/(2 (1 + ν))), 0}} . 
     Inactive[Grad][v[x, y], {x, y}]), {x, y}] + 
   Inactive[Div][({{-((Y (1 - ν))/((1 - 2 ν) (1 + ν))), 0}, {0, -(Y/(2 (1 + ν)))}} . 
     Inactive[Grad][u[x, y], {x, y}]), {x, y}], 
   Inactive[Div][({{0, -(Y/(2 (1 + ν)))}, {-((Y ν)/((1 - 2 ν) (1 + ν))), 0}} . 
     Inactive[Grad][u[x, y], {x, y}]), {x, y}] + 
   Inactive[Div][({{-(Y/(2 (1 + ν))), 0}, {0, -((Y (1 - ν))/((1 - 2 ν) (1 + ν)))}} . 
     Inactive[Grad][v[x, y], {x, y}]), {x, y}]};

This system of differential equation is then solved over a rectangular boundary region.

The boundary conditions are imposed.

1. One fixed end ($x = 0$): Dirchlet boundary condition -> Displacement = 0

2. One lose end ($x = 5$): Neumann boundary condition -> Imposed force of 1 Unit (e.g.. Newton)

{uif, vif} = 
  NDSolveValue[
    {planeStrainOperator[10^3, 33/100] == {0, NeumannValue[-1., x == 5]}, 
    DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, x == 0]}, 
    {u, v}, {x, 0, 5}, {y, 0, 1}];

Plot

Needs["NDSolve`FEM`"]
mesh = uif["ElementMesh"];
Show[{
  mesh["Wireframe"["MeshElement" -> "BoundaryElements"]],
  ElementMeshDeformation[mesh, {uif, vif}][
    "Wireframe"[
      "ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]},
  ImageSize -> 300]

What I would like to do, is to calculate the Deformation, due to some localized forces applied along the top... but I don't know how?

Answer 1

Here is an example with a force that pushes down between the stretch defined by 2<=x<=3 && y==1 and a fixation at the point y==1 && x==5. Applying a force over a stretch is done by a NeumannValue and setting a fixed displacement on a point is done by a DirichletCondition.

{uif, vif} = 
  NDSolveValue[{planeStrainOperator[10^3, 33/100] == {0, 
      NeumannValue[-5., y == 1 && 2 <= x <= 3]}, 
    DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, x == 0], 
    DirichletCondition[v[x, y] == 0.1, y == 1 && x == 5]}, {u, v}, {x,
     0, 5}, {y, 0, 1}];


Needs["NDSolve`FEM`"]
mesh = uif["ElementMesh"];
Show[{mesh["Wireframe"["MeshElement" -> "BoundaryElements"]], 
  ElementMeshDeformation[mesh, {uif, vif}][
   "Wireframe"[
    "ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]},
  ImageSize -> 300]