Specifying NeumannValue for Axial load instead of transverse load (Plane Stress situation)

Question

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).

This axial loading situation is depicted by this lovely figure:

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["DifferentialEquations`InterpolatingFunctionAnatomy`"]
Needs["NDSolve`FEM`"]

xload = 1000;
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))), 
         0}, {0, -(Y/(1 - \[Nu]^2))}}.Inactive[Grad][
        v[x, y], {x, y}]), {x, y}]};


uif, vif} = 
  NDSolveValue[{planeStressOperator[Ey, n] == {0, 
      NeumannValue[-xload, x == span]}, 
    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]

Answer 1

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["DifferentialEquations`InterpolatingFunctionAnatomy`"]
xload = 1000;
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] == 
    {0, NeumannValue[-xload, x == span]}, 
     DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, x == 0]}, 
     {u, v}, {x, 0, span}, {y, 0, width}];
{meshcanti} = InterpolatingFunctionCoordinates[usolcanti[[2]]];
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[[2]][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.