12
$\begingroup$

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]

enter image description here

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?

enter image description here

$\endgroup$

1 Answer 1

15
$\begingroup$

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]

enter image description here

$\endgroup$
6
  • 1
    $\begingroup$ Happy FEM-ing ;-) $\endgroup$
    – user21
    Mar 28, 2016 at 7:30
  • $\begingroup$ 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. $\endgroup$
    – henry
    Mar 28, 2016 at 7:31
  • $\begingroup$ How thick is the tip of your arrow? $\endgroup$
    – user21
    Mar 28, 2016 at 7:33
  • $\begingroup$ 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... $\endgroup$
    – henry
    Mar 28, 2016 at 7:37
  • 2
    $\begingroup$ 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. $\endgroup$
    – user21
    Mar 29, 2016 at 0:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.