# Beam deformation due to localized force distribution

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["NDSolveFEM"]
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?

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["NDSolveFEM"]
mesh = uif["ElementMesh"];
Show[{mesh["Wireframe"["MeshElement" -> "BoundaryElements"]],
ElementMeshDeformation[mesh, {uif, vif}][
"Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]},
ImageSize -> 300]


• Happy FEM-ing ;-) – user21 Mar 28 '16 at 7:30
• Haha ! Thanks. Is there also an option to assign multiple individual seperated point forces as I drew in my question ? Thank you for your help. – henry Mar 28 '16 at 7:31
• How thick is the tip of your arrow? – user21 Mar 28 '16 at 7:33
• Actually, there is no thickness. Those are just point forces. Of Course, one could assign a very thin thickness, but I thought of maybe using a Kronecker Delta within the Weak Form of the FEM to add those Point Forces... – henry Mar 28 '16 at 7:37
• Physically there is no point force (at least on the scale of beams) So you'd need to have some length (area) over which this force applies. What you could to is to translate the point force into a displacement and use a DirichletCondition to enforce that displacement. – user21 Mar 29 '16 at 0:59