Question 1 of 3: Tapered bar modelled with constant-area elements

Problem 1: Tapered bar modelled with constant-area elements (equal value)

Question text not reproduced: the examination questions are © Engineers and Geoscientists BC. Open the official past paper (linked at the top of this page) to read the question, then follow the worked solution below.

Given. A straight bar of length $L$ whose cross-sectional area grows linearly from $A_0$ at the free tip to $A_1$ at the built-in end, loaded by a single axial force $P$ at the tip.

Find. The axial displacement and stress fields predicted by (i) a single constant-area bar element, (ii) two constant-area bar elements, each element taking its area at its own mid-length, and a quantitative comparison of both against the exact solution of the stated boundary-value problem.

[Figure not reproduced: Figure 1 (redrawn). The tapered bar: free tip at $x = 0$ with area $A_0 = 2$ in$^2$, built-in end at $x = L = 20$ in with area $A_1 = 3$ in$^2$, tip load $P = 1000$ lb. See the official exam paper.]

Approach. Replace the taper by one, then two, prismatic two-node bar elements whose areas are sampled at the element mid-lengths; assemble and solve the resulting chain, which is a series of axial springs of stiffness $A_eE/L_e$; then integrate the governing equation in closed form to obtain the exact field and compare.

1.1 One constant-area element

1. Sample the area at the element mid-length. The instruction "evaluate the area at the center of each element length" fixes the element area, for the single element spanning the whole bar, at $x = L/2 = 10\ \text{in}$: $$A^{(1)}=A_0+\left(A_1-A_0\right)\frac{L/2}{L}=2+\frac{10}{20}=2.500\ \text{in}^{2}.$$ This is simply the average of the two end areas, because the taper is linear.
2. Form the element stiffness. A two-node bar element of area $A$, modulus $E$ and length $L_e$ has $$[k]=\frac{AE}{L_e}\begin{bmatrix}1 & -1\\ -1 & 1\end{bmatrix},$$ so the scalar axial stiffness of this element is $$k^{(1)}=\frac{A^{(1)}E}{L}=\frac{(2.500)(10\times10^{6})}{20}=1.250\times10^{6}\ \text{lb/in}.$$ Node 1 sits at the loaded tip and node 2 at the wall.
3. Apply the boundary conditions and solve. Node 2 is built in, $u_2 = 0$, which leaves the single equation $k^{(1)}u_1 = F_1$ with $F_1 = -P$ (the load pulls the tip in the $-x$ sense). Hence $$u_1=\frac{-P}{k^{(1)}}=\frac{-1000}{1.250\times10^{6}}=\boxed{-8.000\times10^{-4}\ \text{in}}$$ which is an elongation of $8.000\times10^{-4}$ in, the tip moving away from the wall.
4. Recover the element stress. A constant-strain bar element has a single strain, and therefore a single stress, across its whole length: $$\sigma^{(1)}=E\,\varepsilon=E\,\frac{u_2-u_1}{L}=\left(10\times10^{6}\right)\frac{0-(-8.000\times10^{-4})}{20}=\boxed{400.0\ \text{psi}}.$$ The same number follows directly from equilibrium, $\sigma = P/A^{(1)} = 1000/2.500 = 400.0$ psi, because a bar element carrying only end loads transmits the full applied force. This identity is the quickest check available on any stepped-bar model.

1.2 Two constant-area elements

1. Sample the two element areas. Each element is $L_e = 10$ in long, so their mid-lengths are at $x = 5$ in and $x = 15$ in: $$A^{(1)}=2+\frac{5}{20}=2.250\ \text{in}^{2},\qquad A^{(2)}=2+\frac{15}{20}=2.750\ \text{in}^{2}.$$ Their average is again 2.500 in$^2$, the one-element value — the refinement redistributes material without changing the total.
2. Form the two element stiffnesses. With $L_e = 10$ in, $$k^{(1)}=\frac{(2.250)(10\times10^{6})}{10}=2.250\times10^{6}\ \text{lb/in},\qquad k^{(2)}=\frac{(2.750)(10\times10^{6})}{10}=2.750\times10^{6}\ \text{lb/in}.$$ Halving the element length roughly doubles each stiffness, which is why the two-element chain is not simply the one-element result repeated.
3. Assemble and reduce. Numbering the tip node 1, the mid node 2 and the wall node 3, the global system is $$\begin{bmatrix}k^{(1)} & -k^{(1)} & 0\\ -k^{(1)} & k^{(1)}+k^{(2)} & -k^{(2)}\\ 0 & -k^{(2)} & k^{(2)}\end{bmatrix}\begin{Bmatrix}u_1\\u_2\\u_3\end{Bmatrix}=\begin{Bmatrix}-P\\0\\R\end{Bmatrix}.$$ Striking out the third row and column for $u_3 = 0$ leaves a two-equation system whose second equation, $-k^{(1)}u_1+(k^{(1)}+k^{(2)})u_2 = 0$, expresses the fact that node 2 carries no applied load.
4. Solve the reduced system. Because the chain is in series with a single tip load, every element carries the same internal force $P$, so it is quickest to work back from the wall: $$u_2=\frac{-P}{k^{(2)}}=\frac{-1000}{2.750\times10^{6}}=\boxed{-3.6364\times10^{-4}\ \text{in}},$$ $$u_1=u_2-\frac{P}{k^{(1)}}=-3.6364\times10^{-4}-\frac{1000}{2.250\times10^{6}}=\boxed{-8.0808\times10^{-4}\ \text{in}}.$$ Substituting these back into the first global equation returns $-P$ exactly, confirming the solve.
5. Recover the two element stresses. Each element again has a single constant stress, $$\sigma^{(1)}=\frac{P}{A^{(1)}}=\frac{1000}{2.250}=\boxed{444.4\ \text{psi}},\qquad \sigma^{(2)}=\frac{P}{A^{(2)}}=\frac{1000}{2.750}=\boxed{363.6\ \text{psi}}.$$ The stress field is therefore a two-step staircase with a jump of $444.4 - 363.6 = 80.8$ psi at the shared node, where the exact field is perfectly smooth. Element forces $A^{(1)}\sigma^{(1)} = A^{(2)}\sigma^{(2)} = 1000$ lb confirm equilibrium in each element separately.

1.3 Comparison with the exact solution

The exact field follows from integrating the governing equation directly, and it is worth doing before commenting on the numbers, because it supplies the yardstick.

1. Integrate the governing equation. With $A(x) = A_0 + (A_1-A_0)x/L$ linear in $x$, $$u(x)=\int\frac{P}{EA(x)}\,dx=\frac{PL}{E\left(A_1-A_0\right)}\ln A(x)+C,$$ and imposing $u(L) = 0$ to fix the constant gives the closed form $$\boxed{u(x)=\frac{PL}{E\left(A_1-A_0\right)}\ln\frac{A(x)}{A_1}}.$$ The logarithm is the signature of a linear taper; a stepped model can only ever approximate it by straight segments.
2. Evaluate the exact displacements at the two nodes of interest. With $A_1 - A_0 = 1\ \text{in}^2$, the leading coefficient is $PL/(E\Delta A) = (1000)(20)/(10\times10^{6})=2.000\times10^{-3}$ in, so $$u(0)=2.000\times10^{-3}\ln\frac{2}{3}=\boxed{-8.1093\times10^{-4}\ \text{in}},\qquad u(10)=2.000\times10^{-3}\ln\frac{2.5}{3}=-3.6464\times10^{-4}\ \text{in}.$$
3. Write the exact stress field. Equilibrium of any cut gives the axial force as $P$ everywhere, so $$\sigma(x)=\frac{P}{A(x)}=\frac{1000}{2+x/20}\ \text{psi},$$ which falls smoothly from $\sigma(0) = 500.0$ psi at the tip through $\sigma(10) = 400.0$ psi to $\sigma(L) = 333.3$ psi at the wall. The strain field is $\varepsilon(x) = \sigma(x)/E$, ranging from $50.0$ to $33.3$ microstrain.
4. Tabulate the displacement errors. Comparing tip values, $$e_1=\frac{-8.0000-(-8.1093)}{-8.1093}=-1.35\ \%,\qquad e_2=\frac{-8.0808-(-8.1093)}{-8.1093}=-0.35\ \%,$$ and at the interior node the two-element model is in error by $-0.28$ per cent. Halving the element size has reduced the tip error by a factor of $1.348/0.351 = 3.84$, essentially the factor of four expected of a second-order method.

The two comparisons are best seen graphically. The displacement plot shows the piecewise-linear finite-element fields chording the exact logarithmic curve, always on the stiff side; the stress plot shows the element-constant staircase straddling the exact hyperbola.

1.4 Comments on the results

Four observations carry the marks here, and all four are visible in the two plots above.

The displacements converge quickly and monotonically from the stiff side. A two-node bar element enforces a linear displacement field, so the finite-element solution can only ever be a chord of the true curve between nodes. Because the true compliance is $\int_0^L dx/A(x)$ and $1/A$ is a convex function of $x$ when $A$ is linear, sampling $A$ at the element mid-point systematically under-estimates that integral. The model is therefore stiffer than the bar it represents, and both tip displacements come out short of the exact value — by 1.35 per cent with one element and 0.35 per cent with two. The near-fourfold error reduction on halving the mesh is the expected $O(h^{2})$ behaviour for this mid-point sampling.

The stresses are far less accurate than the displacements, and they are discontinuous. Displacements are the primary unknowns and are computed from an energy statement; stresses are recovered by differentiating an approximate field, which loses one order of accuracy. With one element the model reports 400.0 psi everywhere, against a true range of 500.0 to 333.3 psi — a 20 per cent under-prediction at the critical tip section. Two elements narrow that worst-case error to 11 per cent (444.4 against 500.0) but introduce an unphysical 80.8 psi jump at the shared node, since nothing in the formulation enforces stress continuity between elements.

Where the element stress is right, it is right for a reason. Each element stress coincides exactly with the true stress at that element's mid-point — 444.4 psi at $x = 5$ in, 363.6 psi at $x = 15$ in and 400.0 psi at $x = 10$ in for the single element. That is not a coincidence: sampling the area at mid-length makes the mid-point the superconvergent stress point of the element, and it is the reason the question insists on that particular sampling rule. In practice this is why post-processors report stresses at interior sampling points rather than at nodes.

The engineering conclusion. Neither mesh is fit for a strength check as it stands: the governing section is the smallest one, at the tip, and both models under-predict its stress. If this bar were being sized, one would either refine the mesh strongly toward the tip, extrapolate the mid-point stresses back to $x = 0$, or — far better — use a tapered (variable-area) bar element whose formulation integrates $A(x)$ exactly and returns the logarithmic displacement field at once. The displacement result, by contrast, is already adequate for a stiffness or deflection-limit check with only two elements.