I am trying to simulate a bar/link element under axial tensile load. The link is of unit thickness and I may use a plane stress model. A plane stress model is an approximation that may be used for thin plates (length and width are one or more orders larger than the thickness). The link/bar element has an axial load, $R$ (top half of figure). This element deforms only by elongation or contraction but no bending. (source: Source: https://en.wikiversity.org/wiki/Nonlinear_finite_elements/Axially_loaded_bar#/media/File:AxialBar.png)

The dimensions of the link are as follows: span ($L$)=20, width=10, thickness=1. The material properties are arbitrarily assigned to Young's modulus (Ey)=$200 \times 10^9$, Poisson's ration ($\nu$)=0.3.

Analytical solution:

There exists a simple analytical solution to find the axial elongation of this link using $\delta = P L / (Area \times Ey)$. If a load of $P=1000$ units acts, axially we can find $\delta=10 \times 10^{-9}$.

If interested, the governing equation to solve this is quite simple:

$Area \times Ey \frac{d^2u}{dx^2} = 0$ with boundary conditions $u(0)=0$ and $u'(x=L)=F/(Area \times Ey)$. The first boundary condition is clearly a Dirichlet while the second boundary condition shows promise as Neumann (to my meager understanding of what Mathematica says).

Mathematica assisted FEM solution:

That was the analytical solution that can be found in a trice. Now, I head to the Mathematica solution with the use of Finite Element Methods.

The Plane Stress Model for structural analysis (FEMDocumentation/tutorial/SolvingPDEwithFEM) has a plane stress operator that may be used (potentially) on any geometry (standard shapes like rectangles to arbitrary shapes).

I define my domain as:

\[CapitalOmega] = ImplicitRegion[True, {{x, 0, span}, {y, 0, width}}];


The Dirichlet boundary conditions that say that there is no deflection where the bar is fixed are defined as:

DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, x == 0]}


How do I define the NeumannValue for the axial load? The example only shows a NeumannValue for the bending of a beam (transverse load) as:

NeumannValue[-1, x == 5]

I have tried a few permutations or combinations with the NeumannValue but to no avail. One example is NeumannValue[-1, x==5 && y] and this turned out to be an error.

My full Mathematica code

Needs["DifferentialEquationsInterpolatingFunctionAnatomy"]
Needs["NDSolveFEM"]

span = 20;
width = 10;
Ey = 200*10^9;
n = 33/100;
iorder = 5;
cellmeasure = 10;
\[CapitalOmega] =
ImplicitRegion[True, {{x, 0, span}, {y, 0, width}}];
planeStressOperator[
Y_, \[Nu]_] := {Inactive[
Div][({{0, -((Y \[Nu])/(1 - \[Nu]^2))}, {-((Y (1 - \[Nu]))/(2 \
(1 - \[Nu]^2))), 0}}.Inactive[Grad][v[x, y], {x, y}]), {x, y}] +
Inactive[
Div][({{-(Y/(1 - \[Nu]^2)),
0}, {0, -((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2)))}}.Inactive[
Grad][u[x, y], {x, y}]), {x, y}],
Inactive[
Div][({{0, -((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2)))}, {-((Y \
\[Nu])/(1 - \[Nu]^2)), 0}}.Inactive[Grad][u[x, y], {x, y}]), {x, y}] +
Inactive[
Div][({{-((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2))),
v[x, y], {x, y}]), {x, y}]};

uif, vif} =
NDSolveValue[{planeStressOperator[Ey, n] == {0,
DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, x == 0]}, {u,
v}, {x, 0, span}, {y, 0, width},
Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> cellmeasure},
"IntegrationOrder" -> iorder}}];


The deflection should be found by uif[span,width] vif[span, width]

• This appears to be a 1D problem. Why is the transverse dimension included? Jul 9, 2015 at 2:29
• @bbgodfrey The transverse direction, although included should not affect the results. The object here is that the plane stress operator could be used for this example. If you look at modern FEM software, this axial problem includes transverse direction but only an axial load. The question here is how to model the axial load in terms of Neumann or other values mma understands. Jul 9, 2015 at 11:29
• A few questions/comments: Extra ) in second line; Γ defined but not used; why are operators Div and Grad made Inactive?; a rectangular mesh might work better than a triangular one. Jul 9, 2015 at 12:39
• @bbgodfrey I have not used the Dirichlet condition gamma although I could have. I'm not sure about the "inactive", I forget what it does. Yes a rectangular mesh works better but I want to compare rect to tri :) Jul 9, 2015 at 12:44
• Inactive prevents the operators in question from executing. Also, why is NDSolveValue executed twice with essentially the same arguments? I recommend that you clean up the code in order to attract responses. Best wishes. Jul 9, 2015 at 12:49

I believe that NeumannValue[-xload, x == span], as given in the question, is correct, as can be seen from the following simple test, for which the solution is obvious.

test[Y_, ν_] := Inactive[Div][({{Y, 0}, {0, Y}}.Inactive[Grad][u[x, y], {x, y}]), {x, y}];
sol = NDSolveValue[{test[Ey, n] == NeumannValue[-xload, x == span],
u[0, y] == 0}, u, {x, 0, span}, {y, 0, width}];
Plot3D[sol[x, y], {x, 0, span}, {y, 0, width}, PlotRange -> All] For completeness, note that the code in the question can be reduced in size to

Needs["DifferentialEquationsInterpolatingFunctionAnatomy"]
span = 20;
width = 10;
Ey = 200*10^9;
n = 33/100;
planeStressOperator[Y_, ν_] :=
{Inactive[Div][({{0, -((Y ν)/(1 - ν^2))}, {-((Y (1 - ν))/(2 (1 - ν^2))), 0}}.
Inactive[Grad][v[x, y], {x, y}]), {x, y}] +
Inactive[Div][({{-(Y/(1 - ν^2)), 0}, {0, -((Y (1 - ν))/(2 (1 - ν^2)))}}.
Inactive[Grad][u[x, y], {x, y}]), {x, y}],
Inactive[Div][({{0, -((Y (1 - ν))/(2 (1 - ν^2)))}, {-((Y ν)/(1 - ν^2)), 0}}.
Inactive[Grad][u[x, y], {x, y}]), {x, y}] +
Inactive[Div][({{-((Y (1 - ν))/(2 (1 - ν^2))), 0}, {0, -(Y/(1 - ν^2))}}.
Inactive[Grad][v[x, y], {x, y}]), {x, y}]};
usolcanti = NDSolveValue[{planeStressOperator[Ey, n] ==
DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, x == 0]},
{u, v}, {x, 0, span}, {y, 0, width}];
{meshcanti} = InterpolatingFunctionCoordinates[usolcanti[]];
meshcanti["Wireframe"] Needs["NDSolveFEM"] and (in this particular case) Method are unnecessary in Ver 10, and eliminating ImplicitRegion leads to a rectangular mesh rather than a triangular one.

Plot3D[usolcanti[][x, y], {x, 0, span}, {y, 0, width}] The results are unchanged from those produced by the code in the question but differ from that of the test case, because the differential equations are different.

• I suppose I am unable to understand what the x and y directions are. Is x along the thickness of this bar element and y along its span? If so, should the load be an +xload? (of course the difference between - and + loads would be the sign of the deflection alone) Jul 10, 2015 at 14:19
• @drN Thank you for accepting the answer. x and y are as you defined them in your question, with x the length along the bar, and y its thickness. Jul 10, 2015 at 18:16
• Yes, I just realized that. Now I am on the hunt for a generalized Hookes law that describes a cantilever with transverse load (original mma example). Would you mind putting the final statement in your answer "The results are unchanged from those produced by the code in the question but differ from that of the test case, because the differential equations are different." in bold? This is vital. I was under the impression that the plane stress operator in full would still solve the linear deflection situation while it solves only the transverse load condition. Jul 10, 2015 at 18:18