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]
)
in second line;Γ
defined but not used; why are operatorsDiv
andGrad
madeInactive
?; a rectangular mesh might work better than a triangular one. $\endgroup$Inactive
prevents the operators in question from executing. Also, why isNDSolveValue
executed twice with essentially the same arguments? I recommend that you clean up the code in order to attract responses. Best wishes. $\endgroup$