Results from 4, 8, and 16 finite element solutions to example [ex17.1]

Improved numerical solutions are obtained by considering uniform meshes of four, eight, sixteen, etc., elements. As shown in figure [fig17.10], the piecewise polynomial approximation for the axial displacement using eight uniform elements is, to the scale of the plot, very close to the exact solution. The axial force from the exact solution and from the finite element solution with eight uniform elements is shown in figure [fig17.11]. The piecewise constant axial force from the solution with eight elements is an improvement with respect to the two element model shown in figure [fig17.9].

The strain energy and the natural boundary condition at \(z=L\) are used to measure the error in the finite element solutions with respect to the exact solution. Results for uniform meshes of one to sixteen elements are listed in table [tab17.1]. The data in the table demonstrates that the strain energy converges faster than the natural boundary condition to the exact solution as the number of elements is increased.

Convergence requirements

Heubner et al. (1995, p. 85) list the following requirements for mathematical convergence of the finite element solution to the exact solution as an increasing number of smaller elements are used in the remeshng process.

• The elements must be made smaller in such a way that every point in the solution domain can always be within an element regardless of how small the element may be.

• All previous meshes must be contained in the refined meshes.

• The form of the interpolation functions must remain unchanged.

For example, the nodes in the mesh for \(\textit{M} = 2\) are also contained in the mesh for \(\textit{M} = 4\), nodes in the mesh for \(\textit{M} = 4\) are also contained in the mesh for \(\textit{M} = 8\), etc. The same linear interpolation functions ([eq17.16]) are used in each discretization.

Apparent loadings from the 8- and 16-element solutions of example [ex17.1]

As illustrated in figure [fig17.11], the axial force exhibits jumps at the nodes between neighboring elements. These jumps can be interpreted as a series of concentrated forces applied at the nodes, and these forces are called the apparent loading (Szabo and Babuska, p. 63). Let \(\underline{Q}_{k}\) denote the apparent axial force at node \(z_{k}\). From ([eq17.42]), the apparent axial force in the positive \(z\)-direction at an interior node is \[\begin{aligned} \label{eq17.44} \underline{Q}_{k}=N^{(k-1)}\left(z_{k}\right)-N^{(k)}\left(z_{k}\right) \quad k=2,3,\,\ldots,\,M-1.\end{aligned}\] The forces at the nodes computed from the jump in the internal axial force are listed in table [tab17.2] for the 8-element model and for the 16-element model. Note that the sum of these forces \(\underline{Q}_{k}\) vanishes in each model, which is a consequence of the principle work as a statement of equilibrium. At common interior nodes the force \(\underline{Q}_{k}\) is smaller in the 16-element model than in the 8-element model.

Adaptive mesh refinement beginning with the 8-element solution to example [ex17.1]

A uniform mesh may converge slowly to the exact solution with continued refinement. In practice, finite element simulations are performed on a structure where the exact solution is not known. For those structures whose exact solutions are unknown, the apparent loading from the finite element solution can be used in adaptive procedures to refine the mesh. Adaptive mesh refinement is based on assessing the relative error of the energy norm in each element between the original loading and the apparent loading (Szabo and Babuska, p. 63). The energy norm for a kinematically admissible function \(w(z)\) is denoted by \(\|w(z)\|\), and it is defined as the square root of the strain energy ([eq17.9]) (i.e., \(\|w\| \equiv \sqrt{U[w]}\)). The mesh and shape functions of the original model are not changed in the second solution of the model subject to the apparent loading. Those elements exhibiting the largest discrepancy in the energy norm between the original loading and the apparent loading are subdivided to generate a new mesh. The new mesh will not be uniform, and is called quasi-uniform. The adaptive mesh procedure is repeated with respect to the new mesh. This repeated use of mesh refinement generates a sequence of meshes. An optimum mesh is achieved when the local error is distributed uniformly through the mesh (Heubner et al., p. 514).

Consider the eight-element model (\(M= 8\)) with nine nodes. The nine nodes in the uniform mesh are \[\begin{aligned} \label{eq17.45} \left\{z_{9}\right\}=\{0,1/8,1/4,3/8,1/2,5/8,3/4,7/8,1\} L.\end{aligned}\]

Each element has the same length \(h_{k}=z_{k+1}-z_{k}=L/8\), \(k=1,2, \ldots, 8\). The dimensions of the restrained structural stiffness matrix \(\left[K_{\alpha \alpha}\right]\) is 8X8. The displacement degrees of freedom and the corresponding forces are \[\label{eq17.46} \hspace*{6pt}\left\{q_{\alpha}\right\}=\left[\begin{array}{@{}llll@{}}q_{2} & q_{3} & q_{4} & q_{5}, q_{6}, q_{7}, q_{8}, q_{9}\end{array}\right]^{T}, \textrm{ and } \left\{Q_{\alpha}\right\}=\left[Q_{2} Q_{3} Q_{4}, Q_{5}, Q_{6}, Q_{7}, Q_{8}, Q_{9}\right]^{T}.\] The nodal force vector for the original loading is \[\begin{aligned} \label{eq17.47} \left\{Q_{\alpha}\right\}=\left[(8 E A/L-c L/48) q_{1}, 0,0,0,0,0,0,0, E A \alpha \tau_{0}+F_{L}\right]^{T}=[-79,183.3,0,0,0,0,0,0,0,-14240 .]^{T}.\end{aligned}\]

The matrix equation \(\left[K_{\alpha \alpha}\right]\left\{q_{\alpha}\right\}=\left\{Q_{\alpha}\right\}\) is solved numerically to determine the displacement vector \(\left\{q_{\alpha}\right\}\). The axial displacement \(w^{(k)}\) in each of the eight elements is computed from \(\left\{q_{\alpha}\right\}\), followed by the computation of the energy norm in each element \(\left\|w^{(k)}\right\|\), \(k=1,2, \ldots, 8\).

For the apparent loading the nodal force vector from table [tab17.2] is

\[\begin{aligned} \label{eq17.48} \left\{\underline{Q}_{\alpha}\right\}=[29,906.9\;12,679.5,510.35\;2,672.57\;1,944.26\;2,750.58\;5,727.95,\;-35,260.1].\end{aligned}\] Without changing the restrained structural stiffness matrix, the equation \(\left[K_{\alpha \alpha}\right]\left\{\underline{q}_{\alpha}\right\}=\left\{\underline{Q}_{\alpha}\right\}\) is solved numerically for the displacement vector \(\left\{\underline{q}_{\alpha}\right\}\). Note that \(\left\{\underline{q}_{\alpha}\right\} \neq\left\{q_{\alpha}\right\}\). From the displacement vector \(\left\{\underline{q}_{\alpha}\right\}\) the axial displacement in the kth element \(\underline{w}^{(k)}\) is determined. The energy norm for the displacement in each element from the apparent loading is denoted by \(\left\|\underline{w}^{(k)}\right\|\), \(k=1,2, \ldots, 8\). Results for the energy norms from the original loading and apparent loading are listed in table [tab17.3]. Elements 1, 8, 2, 7, 3, 5, 4, and 6 have the largest to the smallest discrepancy in the energy norms.

Based on the data in table [tab17.3], a quasi-uniform mesh is selected to increase the number of elements in the domain where the errors in energy norm are large. For example, a mesh of fifteen nodes and fourteen elements is illustrated below: \[\begin{aligned} \label{eq17.49} \left\{z_{15}\right\}=\{0,1/24,1/12,1/8,3/16,1/4,3/8,1/2,5/8,3/4,13/16,7/8,11/12,23/24,1\} L,\end{aligned}\] where the lengths of the elements are \[\begin{aligned} \label{eq17.50} \left\{h_{14}\right\}=\left\{\frac{L}{24}, \frac{L}{24}, \frac{L}{24}, \frac{L}{16}, \frac{L}{16}, \frac{L}{8}, \frac{L}{8}, \frac{L}{8}, \frac{L}{8}, \frac{L}{16}, \frac{L}{16}, \frac{L}{24}, \frac{L}{24}, \frac{L}{24}\right\}.\end{aligned}\] The mesh with nine nodes and the mesh with fifteen nodes are depicted in figure [fig17.12]. Nodes are clustered at the beginning and end of the domain of the 14-element model where discrepancies in the energy norm were the largest. A finite element analysis of this 14-element model resulted in a strain energy of the assembly of 7,788.23 N-mm, and a natural boundary condition \(N(L)=-38{,}164.2\,\mathrm{N}\). Compared to the exact solution the percentage error in the strain energy is 0.438 percent, and the error in the natural boundary condition is 4.59 percent. Moreover, compared to the results from the 16-element model with a uniform mesh in table [tab17.1], the 14-element model has a smaller error in the strain energy and a smaller error in the natural boundary condition.

Element displacement functions and strains

The kth element is denoted by \(\Omega_{k} \equiv\left\{z \mid z_{k} \leq z \leq z_{k+1}\right\}\), where \(z_{k}<z_{k+1}\). The standard element \(\Omega_{\mathrm{st}} \equiv\{\zeta \mid-1<\zeta<1\}\) ([eq17.15]) is mapped to the kth element by eq. ([eq17.16]), and the inverse mapping is given by eq. ([eq17.17]). The lateral displacement of the kth element is denoted by \(v^{(k)}(z)\) and the rotation by \(\phi_{X}^{(k)}(z)\) . Define the generalized nodal displacements as \[\begin{aligned} \label{eq17.66} v^{(k)}\left(z_{k}\right)=q_{2 k-1} \quad \phi_{x}^{(k)}\left(z_{k}\right)=q_{2 k} \quad v^{(k)}\left(z_{k+1}\right)=q_{2 k+1} \quad \phi_{x}^{(k)}\left(z_{k+1}\right)=q_{2 k+2}.\end{aligned}\]

r103pt

See figure [fig17.14]. Admissible functions \(v(z)\) and \(\phi_{x}(z)\) must be continuous within an element and be continuous between elements. The following basis functions will be used for the beam element: \[\begin{aligned} \eta_{1}(\zeta) &=(1-\zeta)/2 \quad \eta_{2}(\zeta)=(1+\zeta)/2 \quad \eta_{3}(\zeta)=\left(\frac{1}{2} \sqrt{\frac{3}{2}}\right)\left(\zeta^{2}-1\right)\nonumber\\ \eta_{4}(\zeta)&=\left(\frac{1}{2} \sqrt{\frac{5}{2}}\right) \zeta\left(\zeta^{2}-1\right). \label{eq17.67}\end{aligned}\] The basis functions \(\eta_{1}(\zeta)\) and \(\eta_{2}(\zeta)\) are the external shape functions with the interpolation properties ([eq17.18]). The basis functions \(\eta_{3}(\zeta)\) and \(\eta_{4}(\zeta)\) are called internal shape functions. The internal shape functions vanish at the nodes \(\zeta=\mp 1\). The selection of the internal shape functions is based on Legendre polynomials. From Szabo and Babuska (1991, p. 38) the expressions for the internal shape function are \[\begin{aligned} \label{eq17.68} \eta_{j}(\zeta)=\frac{1}{\sqrt{2(2 j-1)}}\left[P_{j}(\zeta)-P_{j-1}(\zeta)\right] \quad j=3,4,\end{aligned}\] where \(P_{j}(\zeta)\) is the Legendre polynomial of degree \(j\) in \(\zeta\).² Graphs of these basis functions are shown in figure [fig17.15].

Consider the strain \(\kappa=d \phi_{x}/d z\) appearing in the bending moment ([eq17.52]), and the rotation \(\phi_{x}\) appearing in the transverse shear strain \(\psi_{y}(\zeta)\) ([eq17.53]). As discussed by Reddy (2019, p. 294), in order to avoid the numerical problem of “shear locking” in the thin beam limit where \(\psi_{y} \rightarrow 0\), a consistent interpolation procedure is employed. That is, a consistent interpolation of the shear strain \(\psi_{y}(\zeta)\) requires that the polynomial for the displacement \(v(\zeta)\) be one degree greater in \(\zeta\) than the polynomial for the rotation \(\phi_{x}(\zeta)\). The beam element developed here is capable of representing a linear distribution of the bending moment \(M_{x}(\zeta)\) and a constant shear force \(V_{y}\). Then \(d \phi_{x}/d z\) is linear in \(\zeta\), which implies the rotation is quadratic in \(\zeta\). It follows that displacement \(v(\zeta)\) is cubic in \(\zeta\). The trial functions for the kth element are \[\begin{aligned} \label{eq17.69} \begin{split}\phi_{x}^{(k)}(\zeta)&=\eta_{1}(\zeta) q_{2 k}+\eta_{2}(\zeta) q_{2 k+2}+\eta_{3}(\zeta) a_{2}^{(k)} \\ v^{(k)}(\zeta)&=\eta_{1}(\zeta) q_{2 k-1}+\eta_{2}(\zeta) q_{2 k+1}+\eta_{3}(\zeta) a_{1}^{(k)}+\eta_{4}(\zeta) a_{3}^{(k)} \end{split}.\end{aligned}\] Internal degrees of freedom \(a_{1}^{(k)}\) and \(a_{3}^{(k)}\) are displacements with dimensions of length, and \(a_{2}^{(k)}\) is a rotation in radians. The internal degrees of freedom are not associated with a particular point in \(\Omega_{\mathrm{st}}\). That is, the beam element under consideration does not have internal nodes. Equation ([eq17.69]) is written in the matrix notation as \[\begin{aligned} \label{eq17.70} \left[\begin{array}{@{}l@{}}\phi_{x}^{(k)} \\[3pt]v^{(k)}\end{array}\right]=\left[N_{q}(\zeta)\right]\left\{q^{(k)}\right\}+\left[N_{a}(\zeta)\right]\left\{a^{(k)}\right\},\end{aligned}\] where the 4X1 and the 3X1 vectors of generalized displacements are \[\begin{aligned} \label{eq17.71} \left\{q^{(k)}\right\}=\left[\begin{array}{@{}llll@{}}q_{2 k-1} & q_{2 k} & q_{2 k+1} & q_{2 k+2}\end{array}\right]^{T} \quad\left\{a^{(k)}\right\}=\left[\begin{array}{@{}lll@{}}a_{1}^{(k)} & a_{2}^{(k)} & a_{3}^{(k)}\end{array}\right]^{T}.\end{aligned}\] The shape function matrices are \[\begin{aligned} \label{eq17.72} \left[N_{q}(\zeta)\right]=\left[\begin{array}{@{}cccc@{}}0 & \eta_{1} & 0 & \eta_{2} \\\eta_{1} & 0 & \eta_{2} & 0\end{array}\right] \quad\left[N_{a}(\zeta)\right]=\left[\begin{array}{@{}ccc@{}}0 & \eta_{3} & 0 \\\eta_{3} & 0 & \eta_{4}\end{array}\right].\end{aligned}\] This beam element has seven degrees of freedom. The virtual rotation and displacement for the kth element have the same functional form as ([eq17.70]) but with different coefficients: \[\begin{aligned} \label{eq17.73} \left[\begin{array}{@{}l@{}}\bar{\phi}_{x}^{(k)} \\[3pt]\bar{v}^{(k)}\end{array}\right]=\left[N_{q}(\zeta)\right]\left\{b^{(k)}\right\}+\left[N_{a}(\zeta)\right]\left\{c^{(k)}\right\},\end{aligned}\] where \[\begin{aligned} \label{eq17.74} \left\{b^{(k)}\right\}=\left[\begin{array}{@{}lllll@{}}b_{2 k-1} & b_{2 k} & b_{2 k+1} & b_{2 k+2}\end{array}\right]^{T} \quad\left\{c^{(k)}\right\}=\left[\begin{array}{@{}lll@{}}c_{1}^{(k)} & c_{2}^{(k)} & c 3_{3}^{(k)}\end{array}\right]^{T}.\end{aligned}\] Virtual degrees of freedom \(\left[\begin{array}{@{}llll@{}}b_{2 k-1} & b_{2 k} & b_{2 k+1} & b_{2 k+2}\end{array}\right]\) are associated with the nodal displacements (external shape functions), and virtual degrees of freedom \(\left[c_{1}^{(k)} c_{2}^{(k)} c_{3}^{(k)}\right]\) are associated with the internal shape functions. Elements of the virtual displacement vectors \(\left\{b_{k}\right\}\) and \(\left\{c^{(k)}\right\}\) are independent of the trial functions ([eq17.69]).

Using the chain rule and the inverse mapping ([eq17.17]), the strains for the kth element are \[\begin{aligned} \label{eq17.75} \frac{d \phi_{x}^{(k)}}{d z}=\left(\frac{2}{h_{k}}\right) \frac{d \phi_{x}^{(k)}}{d \zeta} \quad \psi_{y}^{(k)}=\left(\frac{2}{h_{k}}\right) \frac{d v^{(k)}}{d \zeta}+\phi_{x}^{(k)},\end{aligned}\] Substitute ([eq17.70]) and ([eq17.72]) into the expressions for the strains ([eq17.75]) and write the strain-displacement relationship as \[\begin{aligned} \label{eq17.76} \left\{\varepsilon^{(k)}\right\}=\left[\begin{array}{@{}c@{}}\frac{d \phi_{x}^{(k)}}{d z} \\[3pt]\psi_{y}^{(k)}\end{array}\right]=\left[\begin{array}{@{}cccc@{}}0 & \frac{2}{h_{k}} \frac{d \eta_{1}}{d \zeta} & 0 & \frac{2}{h_{k}} \frac{d \eta_{2}}{d \zeta} \\[5pt] \frac{2}{h_{k}} \frac{d \eta_{1}}{d \zeta} & \eta_{1} & \frac{2}{h_{k}} \frac{d \eta_{2}}{d \zeta} & \eta_{2}\end{array}\right]\left[\begin{array}{@{}c@{}}q_{2 k-1} \\[3pt]q_{2 k} \\[3pt]q_{2 k+1} \\[3pt]q_{2 k+2}\end{array}\right]+\left[\begin{array}{@{}ccc@{}}0 & \frac{2}{h_{k}} \frac{d \eta_{3}}{d \zeta} & 0 \\[3pt] \frac{2}{h_{k}} \frac{d \eta_{3}}{d \zeta} & \eta_{3} & \frac{2}{h_{k}} \frac{d \eta_{4}}{d \zeta}\end{array}\right]\left[\begin{array}{@{}l@{}}a_{1}^{(k)} \\[5pt]a_{2}^{(k)} \\[5pt]a_{3}^{(k)}\end{array}\right].\end{aligned}\] In compact notation the strain-displacement relation ([eq17.76]) is \[\begin{aligned} \label{eq17.77} \left\{\varepsilon^{(k)}\right\}=\left[\varepsilon_{q}(\zeta)\right] \left\{q^{(k)}\right\}+\left[\varepsilon_{a}(\zeta)\right]\left\{a^{(k)}\right\},\end{aligned}\] where the strain-displacement matrices associated with the external shape functions and internal shape functions are \[\begin{aligned} \left[\varepsilon_{q}(\zeta)\right]= \left[\begin{array}{@{}cccc@{}} 0 & \frac{-1}{h_{k}} & 0 & \frac{1}{h_{k}} \\[5pt] \frac{-1}{h_{k}} & \frac{(1-\zeta)}{2} & \frac{1}{h_{k}} & \frac{(1+\zeta)}{2} \end{array}\right] \quad\left[\varepsilon_{a}(\zeta)\right]= \left[\begin{array}{@{}ccc@{}} 0 & \left(\frac{\sqrt{6}}{h_{k}} \zeta\right) & 0\\ \left(\frac{\sqrt{6}}{h_{k}} \zeta\right) & \left(\frac{1}{2} \sqrt{\frac{3}{2}}\left(\zeta^{2}-1\right)\right) & \left(\frac{1}{h_{k}} \sqrt{\frac{5}{2}}\left(3 \zeta^{2}-1\right)\right)\end{array}\right]. \label{eq17.78}\end{aligned}\] The virtual strains are polynomials of the same degree as strains ([eq17.77]), but with independent coefficients: \[\begin{aligned} \label{eq17.79} \left\{\bar{\varepsilon}^{(k)}\right\}=\left[\begin{array}{@{}c@{}}\frac{d \bar{\phi}_{x}^{(k)}}{d z} \\[3pt]\bar{\psi}_{y}^{(k)}\end{array}\right]=\left[\varepsilon_{q}(\zeta)\right]\left\{b^{(k)}\right\}+\left[\varepsilon_{a}(\zeta)\right]\left\{c^{(k)}\right\}.\end{aligned}\]

Principle of virtual work for a typical element

The internal virtual work ([eq17.63]) for the kth element in matrix notation is \[\begin{aligned} \label{eq17.80} B\left[\phi_{x}^{(k)}, \bar{\phi}_{x}^{(k)}, v^{(k)}, \bar{v}^{(k)}\right]=\frac{h_{k}}{2} \int_{-1}^{1}\left\{\bar{\varepsilon}^{(k)}\right\}^{T}[D]\left\{\varepsilon^{(k)}\right\} d \zeta,\end{aligned}\] where the material matrix is defined by \[\begin{aligned} \label{eq17.81} [D] \equiv\left[\begin{array}{@{}cc@{}}E I_{x x} & 0 \\0 & s_{y y}\end{array}\right].\end{aligned}\] Substitute the strain-displacement matrices ([eq17.77]) and ([eq17.79]) into ([eq17.80]) to get \[\begin{aligned} B\left[\phi_{x}^{(k)}, \bar{\phi}_{x}^{(k)}, v^{(k)}, \bar{v}^{(k)}\right] &=\frac{h_{k}}{2} \int_{-1}^{1}\left(\left\{b^{(k)}\right\}^{T}\left[\varepsilon_{q}(\zeta)\right]^{T}+\left\{c^{(k)}\right\}^{T}\left[\varepsilon_{a}(\zeta)\right]^{T}\right)\nonumber \\ &\quad\left([D]\left[\varepsilon_{q}(\zeta)\right]\left\{q^{(k)}\right\} +[D]\left[\varepsilon_{a}(\zeta)\right]\left\{a^{(k)}\right\}\right) d \zeta. \label{eq17.82}\end{aligned}\] Perform the matrix multiplications in ([eq17.82]) and arrange the result to \[\begin{aligned} \label{eq17.83} &\begin{array}{rl}B\left[\phi_{x}^{(k)}, \bar{\phi}_{x}^{(k)}, v^{(k)}, \bar{v}^{(k)}\right]\end{array}\nonumber\\ &\quad = \left\{b^{(k)}\right\}^{T}\left\{\left(\frac{h_{k}}{2} \int^{1}_{-1}\left[\varepsilon_{q}(\zeta)\right]^{T}[D]\left[\varepsilon_{q}(\zeta)\right] d \zeta\right)\left\{q^{(k)}\right\}+\left(\frac{h_{k}}{2} \int^{1}_{-1}\left[\varepsilon_{q}(\zeta)\right]^{T}[D]\left[\varepsilon_{a}(\zeta)\right]\right\}\left\{a^{(k)}\right\}\right\} \nonumber\\ &\quad\quad\!\!+ \left. \left\{c^{(k)}\right\}^{T}\left\{\left(\frac{h_{k}^{1}}{2} \int^{1}_{-1}\left[\varepsilon_{a}(\zeta)\right]^{T}[D]\left[\varepsilon_{q}(\zeta)\right] d \zeta\right)\left\{q^{(k)}\right\}+\left(\frac{h_{k}}{2} \int^{1}_{-1}\left[\varepsilon_{a}(\zeta)\right]^{T}[D]\left[\varepsilon_{a}(\zeta)\right]\right\}\left\{a^{(k)}\right\}\right\} \right\}.\end{aligned}\] The expression for the internal virtual work ([eq17.83]) is written as \[\begin{aligned} \label{eq17.84} B\left[\phi_{x}^{(k)}, \bar{\phi}_{x}^{(k)}, v^{(k)}, \bar{v}^{(k)}\right]=\left\{b_{k}\right\}^{T}\left\{\left[K_{q q}\right]\left\{q^{(k)}\right\}+\left[K_{q a}\right]\left\{a^{(k)}\right\}\right\}+\left\{c^{(k)}\right\}^{T}\left\{\left[K_{a q}\right]\left\{q^{(k)}\right\}+\left[K_{a a}\right]\left\{a^{(k)}\right\}\right\},\end{aligned}\] where the stiffness matrices are defined by \[\begin{gathered} \left[K_{q q}\right]=\frac{h_{k}}{2} \int_{-1}^{1}\left[\varepsilon_{q}(\zeta)\right]^{T}[D]\left[\varepsilon_{q}(\zeta)\right] d \zeta \quad\left[K_{q a}\right]=\frac{h_{k}}{2} \int_{-1}^{1}\left[\varepsilon_{q}(\zeta)\right]^{T}[D]\left[\varepsilon_{a}(\zeta)\right], \textrm{ and}\label{eq17.85}\\ \left[K_{a q}\right]=\frac{h_{k}}{2} \int_{-1}^{1}\left[\varepsilon_{a}(\zeta)\right]^{T}[D]\left[\varepsilon_{q}(\zeta)\right] d \zeta \quad\left[K_{a a}\right]=\frac{h_{k}}{2} \int_{-1}^{1}\left[\varepsilon_{a}(\zeta)\right]^{T}[D]\left[\varepsilon_{a}(\zeta)\right]. \label{eq17.86}\end{gathered}\] Evaluate the stiffness matrices to find \[\begin{gathered} \left[K_{q q}\right]= \left[\begin{array}{@{}cccc@{}} \frac{s_{y y}}{h_{k}} & \frac{-s_{y y}}{2} & \frac{-s_{y y}}{h_{k}} & \frac{-s_{y y}}{2} \\[5pt] \frac{-s_{y y}}{2} & \frac{E I_{x x}}{h_{k}}+\frac{h_{k} s_{y y}}{3} & \frac{s_{y y}}{2} & \frac{-E I_{xx}}{h_{k}}+\frac{h_{k} s_{y y}}{6} \\[5pt] \frac{-s_{y y}}{h_{k}} & \frac{s_{y y}}{2} & \frac{s_{y y}}{h_{k}} & \frac{s_{y y}}{2} \\[5pt] \frac{-s_{y y}}{2} & \frac{-E I_{x x}}{h_{k}}+\frac{h_{k} s_{y y}}{6} & \frac{s_{y y}}{2} & \frac{E I_{x x}}{h_{k}}+\frac{h_{k} s_{y y}}{3}\end{array}\right], \textrm{ and}\label{eq17.87}\\ \left[K_{q a}\right] =\left[K_{a q}\right]^{T} =\left[\begin{array}{@{\ }c;{2pt/2pt}c;{2pt/2pt}c@{\ }} 0 & \left(s_{y y}/\sqrt{6}\right) & 0 \\ \left(-s_{y y}/\sqrt{6}\right) & \left(-h_{k} s_{y y} /(2 \sqrt{6})\right) & 0 \\ 0 & \left(-s_{y y}/\sqrt{6}\right) & 0 \\ \left(s_{y y}/\sqrt{6}\right) & \left(-h_{k} s_{y y} /(2 \sqrt{6})\right) & 0\end{array}\right] \quad \left[K_{a a}\right]= \left[\begin{array}{@{}cc;{2pt/2pt}c@{}} \frac{2 s_{y y}}{h_{k}} & 0 & 0 \\ 0 & \frac{E I_{x x}}{h_{k}}+\frac{h_{k} s_{y y}}{5} & \frac{s_{y y}}{\sqrt{15}} \\[5pt] 0 & \frac{s_{y y}}{\sqrt{15}} & \frac{2 s_{y y}}{h_{k}}\end{array}\right]. \label{eq17.88}\end{gathered}\] Take \(\left\{b^{(k)}\right\}=\left\{q^{(k)}\right\}\) and \(\left\{c^{(k)}\right\}=\left\{a^{(k)}\right\}\) in ([eq17.84]), and multiply the result by one-half, to find the strain energy in kth element as \[\begin{aligned} \label{eq17.89} U^{(k)}&=\frac{1}{2} B\left[\phi_{x}^{(k)}, \phi_{x}^{(k)}, v^{(k)}, v^{(k)}\right]\nonumber\\ &=\frac{1}{2}\left(\left\{q^{(k)}\right\}^{T}\left[K_{q q}\right]\left\{q^{(k)}\right\}+2\left\{q^{(k)}\right\}^{T}\left[K_{q q}\right]\left\{a^{(k)}\right\}+\left\{a^{(k)}\right\}^{T}\left[K_{a a}\right]\left\{a^{(k)}\right\}\right).\end{aligned}\] From the mapping ([eq17.16]), the distributed line load intensity and the thermal gradient are evaluated in the kth element as \[\label{eq17.90} \begin{aligned} &f_{y}\left[z_{k}(\zeta)\right]=f_{y}\left[z_{k} \eta_{1}(\zeta)+z_{k+1} \eta_{2}(\zeta)\right]=f_{y}^{(k)}(\zeta) \\ &\tau_{y}\left[z_{k}(\zeta)\right]=\tau_{y}\left[z_{k} \eta_{1}(\zeta)+z_{k+1} \eta_{2}(\zeta)\right]=\tau_{y}^{(k)}(\zeta) \end{aligned}.\] The external virtual work ([eq17.65]) for the kth element is \[\begin{aligned} \label{eq17.91} F\left[\bar{\phi}_{x}^{(k)}, \bar{v}^{(k)}\right]=\left.\bar{v}^{(k)} V_{y}^{(k)}\right|_{z_{k}} ^{z_{k+1}}+\left.\bar{\phi}_{x}^{(k)} M_{x}^{(k)}\right|_{z_{k}} ^{z_{k+1}}+\int_{-1}^{1} \bar{v}^{(k)} f_{y}^{(k)}(\zeta) \frac{h_{k}}{2} d \zeta+\int_{-1}^{1}\left(E I_{x x} \frac{d \bar{\phi}_{x}^{(k)}}{d \zeta} \frac{2}{h_{k}}\right) \alpha \tau_{y}^{(k)}(\zeta)\left(\frac{h_{k}}{2}\right) d \zeta.\end{aligned}\] The boundary terms for the kth element are expressed in the form \[\begin{aligned} \label{eq17.92} \left.\bar{v}^{(k)} V_{y}^{(k)}\right|_{z_{k}} ^{z_{k+1}}+\left.\bar{\phi}_{x}^{(k)} M_{x}^{(k)}\right|_{z_{k}} ^{z_{k+1}}=b_{2 k+1} Q_{2 k+1}^{(k)}+b_{2 k-1} Q_{2 k-1}^{(k)}+b_{2 k+2} Q_{2 k+2}^{(k)}+b_{2 k} Q_{2 k}^{(k)},\end{aligned}\] where generalized forces acting on the cross sections at \(z= z_{\textrm{k}}\) and \(z=z_{\textrm{k}+1}\) are defined by \[\begin{aligned} \label{eq17.93} Q_{2k-1}^{(k)} \equiv-V_{y}\left(z_{k}\right) \quad Q_{2k}(k) \equiv-M_{x}\left(z_{k}\right) \quad Q_{2 k+1}(k) \equiv V_{y}\left(z_{k+1}\right) \quad Q_{2 k+2}^{(k)} \equiv M_{x}\left(z_{k+1}\right). \hspace*{6pt}\end{aligned}\] The integral term with respect to the line load \(f_{y}^{(k)}(\zeta)\) in ([eq17.91]) expands to \[\begin{aligned} \label{eq17.94} \int_{1-1}^{1} \bar{v}^{(k)} f_{y}^{(k)}(\zeta) \frac{h_{k}}{2} d \zeta=\frac{h_{k}}{2}\left(\int_{-1}^{1}\left(b_{2 k-1} \eta_{1}+b_{2 k+1} \eta_{2}+c_{1}^{(k)} \eta_{3}+c_{3}^{(k)} \eta_{4}\right) f_{y}^{(k)}(\zeta) d \zeta\right).\end{aligned}\] The integral term with respect to the temperature gradient \(\tau_{y}^{(k)}(\zeta)\) in ([eq17.91]) expands to \[\begin{aligned} \label{eq17.95} \int_{-1}^{1}\left(E I_{x x} \frac{d \bar{\phi}_{x}^{(k)}}{d \zeta} \frac{2}{h_{k}}\right) \alpha \tau_{y}^{(k)}(\zeta)\left(\frac{h_{k}}{2}\right) d \zeta=\int_{-1}^{1} E I_{x x}\left(b_{2 k} \frac{d \eta_{1}}{d \zeta}+b_{2 k+2} \frac{d \eta_{2}}{d \zeta}+c_{2}^{(k)} \frac{d \eta_{3}}{d \zeta}\right) \alpha \tau_{y}^{(k)}(\zeta) d \zeta.\end{aligned}\] Combining ([eq17.92]), ([eq17.94]), and ([eq17.95]), the external virtual work ([eq17.91]) is written in the matrix form \[\begin{aligned} \label{eq17.96} F\left[\bar{\phi}_{x}^{(k)}, \bar{v}^{(k)}\right]=\left\{b^{(k)}\right\}^{T}\left(\left\{Q^{(k)}\right\}+\left\{F E^{(k)}\right\}\right)+\left\{c^{(k)}\right\}^{T}\left\{F I^{(k)}\right\},\end{aligned}\] where the generalized nodal force vector is denoted by \(\left\{Q^{(k)}\right\}\), the generalized force vector associated with the external shape functions by \(\left\{F E^{(k)}\right\}\), and the generalized force vector associated with the internal shape functions by \(\left\{F I^{(k)}\right\}\). These vectors are \[\begin{aligned} \label{eq17.97} \left\{Q^{\{k\}}\right\} \equiv \left[\begin{array}{@{}c@{}} Q_{2 k-1}^{(k)} \\[4pt] Q_{2 k} \\[4pt] Q_{2 k+1}^{k} \\[4pt] Q_{2 k+2}^{(k)} \\[4pt] \end{array}\right], \left\{F E^{(k)}\right\} = \left[\begin{array}{@{}c@{}} F E_{2 k-1}^{(k)} \\[4pt] F E_{2 k}^{(k)} \\[4pt] F E_{2 k+1}^{(k)} \\[4pt] F E_{2 k+2}^{(k)} \end{array}\right] = \left[ \begin{array}{@{}c@{}} \frac{h_{k}}{2} \int^1_{-1}\eta_{1} f_{y}^{(k)} d\zeta \\[4pt] \int^1_{-1} E I_{x x} \frac{d \eta_{1}}{d \zeta} \alpha \tau_{y}^{(k)} d \zeta \\[4pt] \frac{h_{k}}{2} \int^1_{-1} \eta_{2} f_{y}^{(k)} d \zeta \\[4pt] \int^1_{-1} E I_{x x} \frac{d \eta_{2}}{d \zeta} \alpha \tau_{y}^{(k)} d \zeta \end{array} \right]\!, \left\{F I^{(k)}\right\} =\left[\begin{array}{@{}c@{}} F I_{1}^{(k)} \\[4pt] F I_{2}^{(k)} \\[4pt] F I_{3}^{(k)} \end{array}\right] = \left[\begin{array}{@{}c@{}} \frac{h_{k}}{2} \int^1_{-1} \eta_{3} f_{y}^{(k)} d \zeta \\[4pt] \int^1_{-1} E I_{x x} \frac{d \eta_{3}}{d \zeta} \alpha \tau_{y}^{(k)} d \zeta\\[4pt] \frac{h_{k}}{2} \int^1_{-1} \eta_{4} f_{y}^{(k)} d \zeta\end{array}\right].\end{aligned}\]

Equate the internal virtual work ([eq17.84]) to the external virtual work ([eq17.96]) to get \[\begin{aligned} \label{eq17.98} &\left\{b^{(k)}\right\}^{T}\left(\left[K_{q q}\right]\left\{q^{(k)}\right\}+\left[K_{q a}\right]\left\{a^{\{k\}}\right\}\right)+\left\{c^{\{k\}}\right\}^{T}\left(\left[K_{a q}\right]\left\{q^{\{k\}}\right\}+\right. {\left.\left[K_{a a}\right]\left\{a^{(k)}\right\}\right)} \nonumber\\ &\quad =\left\{b^{(k)}\right\}^{T}\left(\left\{Q^{(k)}\right\}+\left\{F E^{(k)}\right\}\right)+\left\{c^{\{k\}}\right\}{ }^{T}\left\{F I^{(k)}\right\}.\end{aligned}\] The principle of virtual work leads to the governing matrix relations for the kth element: \[\begin{gathered} \left[K_{q q}\right]\left\{q^{(k)}\right\}+\left[K_{q a}\right]\left\{a^{(k)}\right\}=\left\{Q^{(k)}\right\}+\{F E(k)\} \quad \boldsymbol{\forall}\left\{b^{(k)}\right\} \neq 0_{4 X 1}, \label{eq17.99} \\ \left[K_{a q}\right]\left\{q^{(k)}\right\}+\left[K_{a a}\right]\left\{a^{(k)}\right\}=\left\{F I^{(k)}\right\} \quad \boldsymbol{\forall}\left\{c^{(k)}\right\} \neq 0_{3 X 1}. \label{eq17.100}\end{gathered}\]

Static condensation

The internal degrees of freedom \(a_{1}^{(k)}\), \(a_{2}^{(k)}\), and \(a_{3}^{(k)}\) associated with the internal shape functions do not affect interelement continuity. To reduce the number of system equations to solve resulting from the assembly of elements, the internal degrees of freedom can be eliminated at the element level in the process called static condensation. Further, the assembly algorithm is simplified if only the nodal displacements, or external degrees of freedom are employed. The internal degrees of freedom are eliminated in terms of the nodal displacements \(\left\{q^{(k)}\right\}\) and the generalized force vector \(\left\{F I^{(k)}\right\}\) by solving ([eq17.100]). Thus, from ([eq17.100]) \[\begin{aligned} \label{eq17.101} \left\{a^{(k)}\right\}=\left[K_{a a}\right]^{-1}\left\{F I^{(k)}\right\}-\left[K_{a a}\right]^{-1}\left[K_{a q}\right]\left\{q^{(k)}\right\}=\left[K_{a a}\right]^{-1}\left\{F I^{(k)}\right\}+\left[G_{a q}\right]\left\{q^{(k)}\right\},\end{aligned}\] where \[\begin{aligned} \label{eq17.102} \left[K_{a a}\right]^{-1} =\left[\begin{array}{@{}ccc@{}} \frac{h_{k}}{2 s_{y y}} & 0 & 0 \\[4pt] 0 & \frac{6 h_{k} \mu}{E I_{x x}} & \frac{-\sqrt{\frac{3}{5}} h_{k}^{2} \mu}{E I_{x x}} \\[4pt] 0 & \frac{-\sqrt{\frac{3}{5}} h_{k}^{2} \mu}{E I_{x x}} & \frac{30 E I_{x x} h_{k}+3 h_{k}^{3} s_{y y}}{60 E I_{x x} s_{y y}+5 h_{k}^{2} s_{y y}^{2}}\end{array}\right] \left[G_{a q}\right] =\left[\begin{array}{@{}cccc@{}} 0 & \frac{h_{k}}{2 \sqrt{6}} & 0 & \frac{-h_{k}}{2 \sqrt{6}} \\[4pt] \frac{-\sqrt{6} h_{k} s_{y y} \mu}{E I_{x x}} & \frac{\sqrt{\frac{3}{2}} h_{k}^{2} s_{y y} \mu}{E I_{x x}} & \frac{\sqrt{6} h_{k} s_{y y} \mu}{E I_{x x}} & \frac{\sqrt{\frac{3}{2}} h_{k}^{2} s_{y y} \mu}{E I_{x x}} \\[4pt] \frac{h_{k}^{2} s_{y y} \mu}{\sqrt{10} E I_{x x}} & \frac{-h_{k}^{3} s_{y y} \mu}{2 \sqrt{10} E I_{x x}} & \frac{-h_{k}^{2} s_{y y} \mu}{\sqrt{10} E I_{x x}} & \frac{-h_{k}^{3} s_{y y} \mu}{2 \sqrt{10} E I_{x x}}\end{array}\right].\end{aligned}\] We define \(\left[G_{aq}\right]\) as the Guyan matrix, because the method of static condensation was first proposed by Guyan (1965), and it was also introduced by Irons (1965). In eq. ([eq17.102]) we have introduced the dimensionless factor \[\begin{aligned} \label{eq17.103} \mu=\frac{E I_{x x}}{12 E I_{x x}+h_{k}^{2} s_{y y}}.\end{aligned}\] Substitute the displacement vector \(\left\{a^{(k)}\right\}\) from ([eq17.101]) into eq. ([eq17.99]) and write the final result as \[\begin{aligned} \label{eq17.104} \left[K^{(k)}\right]\left\{q^{(k)}\right\}=\left\{Q^{(k)}\right\}+\left\{F^{(k)}\right\},\end{aligned}\] where \[\begin{aligned} \label{eq17.105} \begin{gathered} {\left[K^{(k)}\right]=\left[K_{q q}\right]+\left[K_{q a}\right]\left[G_{a q}\right]} \\ \left\{F^{(k)}\right\}=\left\{F E^{(k)}\right\}-\left[K_{q a}\right]\left[K_{a a}\right]^{-1}\left\{F I^{(k)}\right\}=\left\{F E^{(k)}\right\}+\left[G_{a q}\right]^{T}\left\{F I^{(k)}\right\}\end{gathered}.\end{aligned}\] The 4X4 matrix \(\left[K^{(k)}\right]\) is the condensed stiffness matrix for the kth element, and the 4X1 vector \(\left\{F^{(k)}\right\}\) is the condensed generalized force vector for the kth element. The generalized forces acting on the beam element separated from the nodes, and the generalized forces acting on the nodes, are depicted in figure [fig17.16].

The 4X4 condensed stiffness matrix is \[\begin{aligned} \label{eq17.106} \left[K^{(k)}\right]=\mu\left[\begin{array}{@{}cccc@{}}\frac{12 s_{y y}}{h_{k}} & -6 s_{y y} & \frac{-12 s_{y y}}{h_{k}} & -6 s_{y y} \\[4pt]-6 s_{y y} & \frac{12 E I_{x x}}{h_{k}}+4 h_{k} s_{y y} & 6 s_{y y} & \frac{-12 E I_{x x}}{h_{k}}+2 h_{k} s_{y y} \\[4pt]\frac{-12 s_{y y}}{h_{k}} & 6 s_{y y} & \frac{12 s_{y y}}{h_{k}} & 6 s_{y y} \\[4pt]-6 s_{y y} & \frac{-12 E I_{x x}}{h_{k}}+2 h_{k} s_{y y} & 6 s_{y y} & \frac{12 E I_{x x}}{h_{k}}+4 h_{k} s_{y y}\end{array}\right].\end{aligned}\]

Bending moment and shear force

Substitute vector \(\left\{a^{(k)}\right\}\) from ([eq17.101]) into the strain vector ([eq17.77]) to get \[\begin{aligned} \label{eq17.107} \left\{\varepsilon^{(k)}\right\}=\left[\varepsilon_{q}(\zeta)\right]\left\{q^{(k)}\right\}+\left[\varepsilon_{a}(\zeta)\right]\left\{a^{(k)}\right\}=\left(\left[\varepsilon_{q}(\zeta)\right]+\left[\varepsilon_{a}\right]\left[G_{a q}\right]\right)\left\{q^{(k)}\right\}+\left[\varepsilon_{a}(\zeta)\right]\left[K_{a a}\right]^{-1}\left\{F I^{(k)}\right\}.\end{aligned}\] The bending moment and shear force are determined from the matrix relation \[\begin{aligned} \label{eq17.108} \left[\begin{array}{@{}c@{}}M_{x}^{(k)} \\V_{y}^{(k)}\end{array}\right]=[D]\left(\left\{\varepsilon^{(k)}\right\}-\left[\begin{array}{@{}c@{}}\alpha \tau_{y}^{(k)} \\0\end{array}\right]\right).\end{aligned}\] Substitute the strain ([eq17.107]) into ([eq17.108]) to find the matrix relation for the moment and shear: \[\begin{aligned} \label{eq17.109} \left[\begin{array}{@{}c@{}}M_{x}^{(k)} \\V_{y}^{(k)}\end{array}\right]=[D]\left(\left[\varepsilon_{q}(\zeta)\right]+\left[\varepsilon_{a}\right]\left[G_{a q}\right]\right)\left\{q^{(k)}\right\}+[D]\left(\left[\varepsilon_{a}(\zeta)\right]\left[K_{a a}\right]^{-1}\left\{F I^{(k)}\right\}-\left[\begin{array}{@{}c@{}}\alpha \tau_{y}^{(k)} \\0\end{array}\right]\right).\end{aligned}\] The previous equation is written as \[\begin{aligned} \label{eq17.110} \left[\begin{array}{@{}c@{}}M_{x}^{(k)} \\V_{y}^{(k)}\end{array}\right]=\left[S_{q}(\zeta)\right]\left[\begin{array}{@{}c@{}}q_{2 k-1} \\q_{2 k} \\q_{2 k+1} \\q_{2 k+2}\end{array}\right]+\left[S_{F}(\zeta)\right]\left[\begin{array}{@{}l@{}}F I_{1}^{(k)} \\F I_{2}^{(k)} \\F I_{3}^{(k)}\end{array}\right]-\left[\begin{array}{@{}c@{}}E I_{x x} \alpha \tau_{y}^{(k)}(\zeta) \\0\end{array}\right],\end{aligned}\] where the stress matrices are defined as \[\begin{aligned} \left[S_{q}(\zeta)\right]&= \left[\begin{array}{@{}cccc@{}} -6 s_{y y} \mu \zeta & \left(\frac{-E I_{x x}}{h_{k}}+3 h_{k} s_{y y} \mu \zeta\right) & 6 s_{y y} \mu \zeta & \left(\frac{E I_{x x}}{h_{k}}+3 h_{k} s_{y y} \mu \zeta\right) \\ \frac{-12 \mu s_{y y}}{h_{k}} & 6 s_{y y} \mu & \frac{12 \mu s_{y y}}{h_{k}} & 6 s_{y y} \mu \end{array}\right], \textrm{ and}\label{eq17.111}\\ \left[S_{F}(\zeta)\right]&=\left[\begin{array}{@{}ccc@{}} 0 & 6 \sqrt{6} \mu \zeta & -3 \sqrt{\frac{2}{5}} h_{k} \mu \zeta \\[5pt] \sqrt{\frac{3}{2}} \zeta & \left(\frac{-\sqrt{6} h_{k} s_{y y} \mu}{E I_{x x}}\right) & \frac{3\left[20EI_{xx}(-1+3\zeta^2) + h^2_kS_{yy}(-1+5\zeta^2)\right]\mu}{2\sqrt{10}EI_{xx}} \end{array}\right].\label{eq17.112}\end{aligned}\]

Requirements of the interpolation functions.

The remeshing procedure for an increasing number of elements presented in article 1.2.2 has the old mesh embedded in the new mesh. Monotonic convergence of this sequence of finite element meshes also requires the interpolation functions to be compatible and complete (Bathe, 1982). Compatibility means the displacements must be continuous within the elements and between elements. Completeness means the displacement functions of the element must be able to represent the rigid body displacements and uniform strain states. Consequently, the beam element under consideration must be capable of representing zero strain states for rigid body displacements when the element is not subject to external loads. Note that the determinate of the stiffness matrix ([eq17.106]) is zero, since the element is not restrained against rigid body displacements. Now consider the response of the beam element under the two rigid body modes: \[\begin{aligned} \label{eq17.113} \big\{q^{(k)}\big\}_{1}=\left[\begin{array}{@{}llll@{}}1 & 0 & 1 & 0\end{array}\right]^{T} \textrm{ and } \big\{q^{(k)}\big\}_{2}=\left[\begin{array}{@{}llll@{}}h_{k}/2 & 1 & -h_{k}/2 & 1\end{array}\right]^{T}.\end{aligned}\] Generalized displacement vector \(\left\{q^{(k)}\right\}_{1}\) is a vertical displacement of the element, and the generalized displacement vector \(\left\{q^{(k)}\right\}_{2}\) is a clockwise rotation of the element about its center. In the absence of external loading the generalized force vectors \(\left\{Q^{(k)}\right\}=0_{4 X 1}\), \(\left\{F E^{(k)}\right\}=0_{4 X 1}\), \(\left\{F I^{(k)}\right\}=0_{3 X 1}\), and \(\left\{F^{(k)}\right\}=0_{4 X 1}\). The internal degrees of freedom ([eq17.101]) for the rigid body modes and no external loading also vanish; i.e., \[\begin{aligned} \label{eq17.114} \big\{a^{(k)}\big\}_{1}=\left[G_{a q}\right]\big\{q^{(k)}\big\}_{1}=0_{3 X 1} \textrm{ and } \big\{a^{(k)}\big\}_{2}=\left[G_{a q}\right]\big\{q^{(k)}\big\}_{2}=0_{3 X 1}.\end{aligned}\] Also, evaluation of eq. ([eq17.104]) for the rigid body displacements results in \(\left[K^{(k)}\right]\left\{q^{(k)}\right\}_{1}=0_{4 X 1}\) and \(\left[K^{(k)}\right]\left\{q^{(k)}\right\}_{2}=0_{4 X 1}\), which are consistent with vanishing external loads. Finally, the strains ([eq17.77]) in the element vanish for the rigid body displacements: \[\begin{aligned} \label{eq17.115} \big\{\varepsilon^{(k)}\big\}_{1}=\left[\varepsilon_{q}(\zeta)\right]\big\{q^{(k)}\big\}_{1}=0_{2 x 1} \textrm{ and } \big\{\varepsilon^{(k)}\big\}_{2}=\left[\varepsilon_{q}(\zeta)\right]\big\{q^{(k)}\big\}_{2}=0_{2 x 1}.\end{aligned}\] Therefore, the beam element satisfies part of the completeness requirement of vanishing strains under rigid body displacements. Constant strain states are demonstrated in the next two examples.

[ex17.3] A simply supported, uniform beam is shown in figure [fig17.17](a). It is subject to equal and opposite moments \(M_{a}\) at each end, and a temperature gradient through the thickness \(\tau_{y}\) that is uniform along the length of the beam. The lateral distributed load intensity \(f_{y}=0\), \(0<z<L\). The beam is modeled with one finite element, and the two nodes are \(z_{1}=0\) and \(z_{2}=L\) as shown in figure [fig17.17](b). The mapping ([eq17.16]) is \(z=\frac{1}{2}(1+\zeta) L\).

The generalized displacement vector \(\left\{q^{(1)}\right\}=\left[\begin{array}{@{}llll@{}}q_{1} & q_{2} & q_{3} & q_{4}\end{array}\right]^{T}\) is partitioned to unknown components \(\left\{q_{\alpha}\right\}\) and known components \(\left\{q_{\beta}\right\}\): \[\big\{q^{(1)}\big\}=\left[\begin{array}{@{}l@{}}\left\{q_{\alpha}\right\} \\\left\{q_{\beta}\right\}\end{array}\right], \left\{q_{\alpha}\right\}=\left[\begin{array}{@{}l@{}}q_{2} \\q_{4}\end{array}\right], \textrm{ and } \left\{q_{\beta}\right\}=\left[\begin{array}{@{}l@{}}q_{1} \\q_{3}\end{array}\right]=\left[\begin{array}{@{}l@{}}0 \\0\end{array}\right]. \label{eq17.3.a}\tag{a}\] The generalized displacement vector ([eq17.69]) is \[\left[\begin{array}{@{}l@{}}\phi_{x}^{(1)}(\zeta) \\v(1)(\zeta)\end{array}\right]=\left[\begin{array}{@{}cccc@{}}0 & \eta_{1} & 0 & \eta_{2} \\\eta_{1} & 0 & \eta_{2} & 0\end{array}\right]\left[\begin{array}{@{}c@{}}0 \\q_{2} \\0 \\q_{4}\end{array}\right]+\left[\begin{array}{@{}ccc@{}}0 & \eta_{3} & 0 \\\eta_{3} & 0 & \eta_{4}\end{array}\right]\left[\begin{array}{@{}c@{}}a_{1} \\a_{2} \\a_{3}\end{array}\right]=\left[\begin{array}{@{}c@{}}\eta_{1}(\zeta) q_{2}+\eta_{2}(\zeta) q_{4}+\eta_{3}(\zeta) a_{2} \\\eta_{3}(\zeta) a_{1}+\eta_{4}(\zeta) a_{3}\end{array}\right]. \label{eq17.3.b}\tag{b}\] Displacement boundary conditions \(v(-1)=0\) and \(v(1)=0\) are satisfied by (b). Matrix ([eq17.106]) represents the unrestrained structural stiffness matrix \(\left[K_{u}\right]\). The rows and columns of \(\left[K_{u}\right]\) are interchanged to the order 2, 4, 1, and 3 to facilitate partitioning it into submatrices: \[\left[K_{u}\right]=\mu\left[\begin{array}{@{}cccc@{}}q_{2} & q_{4} & q_{1} & q_{3} \\\frac{12 E I_{x x}}{L}+4 L s_{y y} & \frac{-12 E I_{x x}}{L}+2 L s_{y y} & -6 s_{y y} & 6 s_{y y} \\[6pt] \frac{-12 E I_{x x}}{L}+2 Ls_{y y} & \frac{12 E I_{x x}}{L}+4 L s_{y y} & -6 s_{y y} & 6 s_{y y} \\[4pt] \hdashline\\[-8pt] -6 s_{y y} & -6 s_{y y} & \frac{12 s_{y y}}{L} & \frac{-12 s_{y y}}{L} \\[6pt] 6 s_{y y} & 6 s_{y y} & \frac{-12 s_{y y}}{L} & \frac{12 s_{y y}}{L}\end{array}\right]=\left[\begin{array}{@{}l@{}}{\big[K_{\alpha \alpha}\big]\left[K_{\alpha \beta}\right]} \\{\left[K_{\beta \alpha}\right]\left[K_{\beta \beta}\right]}\end{array}\right]. \label{eq17.3.c}\tag{c}\] The restrained structural stiffness matrix is \[\left[K_{\alpha \alpha}\right]=\frac{E I_{x x}}{12 E I_{x x}+L^{2} s_{y y}}\left[\begin{array}{@{}ll@{}}\frac{12 E I_{x x}}{L}+4 L s_{y y} & \frac{-12 E I_{x x}}{L}+2 L s_{y y} \\[6pt] \frac{-12 E I_{x x}}{L}+2 L s_{y y} & \frac{12 E I_{x x}}{L}+4 L s_{y y}\end{array}\right], \label{eq17.3.d}\tag{d}\] where \(\mu\) is defined in eq. ([eq17.103]). From ([eq17.97]) actions \(\left.F I_{1}^{(1)}=F I_{2}^{(1)}=F I_{2}^{(1)}\right)=0\), since \(f_{y}=0\) and \(\tau_{y}\) is spatially uniform in \(\zeta\). Hence, from ([eq17.105]) the generalized nodal force vector is \(\left\{F^{(1)}\right\}=\left\{F E^{(1)}\right\}\), which is evaluated from ([eq17.97]). The generalized nodal force vector \(\left\{Q^{(1)}\right\}=\left[Q_{1}^{(1)} Q_{2}^{(1)} Q_{3}^{(1)} Q_{4}^{(1)}\right]^{T}\) is partitioned into known components \(\left\{Q_{0,}\right\}\) and unknown components \(\left\{Q_{\beta}\right\}\): \[\left\{Q_{\alpha}\right\}=\left[\begin{array}{@{}l@{}} Q_{2} \\Q_{4}\end{array}\right]=\left[\begin{array}{@{}c@{}}-M_{a} \\M_{a}\end{array}\right], \textrm{ and } \left\{Q_{\beta}\right\}=\left[\begin{array}{@{}l@{}}Q_{1} \\Q_{3}\end{array}\right]. \label{eq17.3.e}\tag{e}\] The generalized force vector \(\left\{F^{(1)}\right\}=\left[\begin{array}{@{}llll@{}}F_{1}^{(1)} & F_{2}^{(1)} & F_{3}^{(1)} & F_{4}^{(1)}\end{array}\right]^{T}\) ([eq17.97]) is prescribed and is partitioned as \[\left\{F_{\alpha}\right\}=\left[\begin{array}{@{}c@{}}F_{2}^{(1)}\\[4pt] F_{4}^{(1)}\end{array}\right]=\left[\begin{array}{@{}l@{}} \int\limits^{1}_{-1}\frac{\left(-E I_{x x} \alpha \tau_{y}\right)}{2} d\zeta \\ \int\limits^{1}_{-1}\frac{\left(E I_{x x} \alpha \tau_{y}\right)}{2} d\zeta \end{array}\right]=\left[\begin{array}{@{}c@{}}-E I_{x x} \alpha \tau_{y} \\E I_{x x} \alpha \tau_{y}\end{array}\right], \textrm{ and } \left\{F_{\beta}\right\}=\left[\begin{array}{@{}l@{}}F_{1}^{(1)} \\[4pt]F_{3}^{(1)}\end{array}\right]=\left[\begin{array}{@{}l@{}}0 \\0\end{array}\right]. \label{eq17.3.f}\tag{f}\] The matrix equation to determine the unknown displacements is \(\left[K_{\alpha \alpha}\right]\left\{q_{\alpha}\right\}=\left\{Q_{\alpha}\right\}+\left\{F_{\alpha}\right\}\), or \[\frac{E I_{x x}}{12 E I_{x x}+L^{2} s_{y y}}\left[\begin{array}{@{}ll@{}} \frac{12 E I_{x x}}{L}+4 L s_{y y} & \frac{-12 E I_{x x}}{L}+2 L s_{y y} \\[4pt] \frac{-12 E I_{x x}}{L}+2 L s_{y y} & \frac{12 E I_{x x}}{L}+4 L s_{y y} \end{array}\right]\left[\begin{array}{@{}l@{}}q_{2} \\[4pt] q_{4}\end{array}\right]=\left[\begin{array}{@{}c@{}}-M_{a}-E I_{x x} \alpha \tau_{y} \\[4pt] M_{a}+E I_{x x} \alpha \tau_{y}\end{array}\right]. \label{eq17.3.g}\tag{g}\] The solution for \(\left\{q_{\alpha}\right\}\) from eq. ([eq17.3.g]) is \[\left\{q_{\alpha}\right\}=\left[\begin{array}{@{}l@{}}q_{2} \\q_{4}\end{array}\right]=\frac{L\left(M_{a}+E I_{x x} \alpha \tau_{y}\right)}{2 E I_{x x}}\left[\begin{array}{@{}c@{}}-1 \\1\end{array}\right]. \label{eq17.3.h}\tag{h}\] The rotations are equal magnitude and of the opposite sense. The reactive force vector is determined from \(\left[K_{\beta \alpha}\right]\left\{q_{\alpha}\right\}+\left[K_{\beta \beta}\right]\left\{q_{\beta}\right\}=\left\{Q_{\beta}\right\}+\left\{F_{\beta}\right\}\), which is \[\frac{E I_{x x}}{12 E I_{x x}+L^{2} s_{y y}}\left[\begin{array}{@{}cc@{}}-6 s_{y y} & -6 s_{y y} \\6 s_{y y} & 6 s_{y y}\end{array}\right] \frac{L\left(M_{a}+E I_{x x} \alpha \tau_{y}\right)}{2 E I_{x x}}\left[\begin{array}{@{}c@{}}-1 \\1\end{array}\right]=\frac{L\left(M_{a}+E I_{x x} \alpha \tau_{y}\right)}{2\left(12 E I_{x x}+L^{2} s_{y y}\right)}\left[\begin{array}{@{}c@{}}6 s_{y y}-6 s_{y y} \\-6 s_{y y}+6 s_{y y}\end{array}\right]=\left[\begin{array}{@{}l@{}}0 \\0\end{array}\right]. \label{eq17.3.i}\tag{i}\] Hence, the reactive force vector \(\left\{Q_{\beta}\right\}\) vanishes. The vector \(\left\{a^{(1)}\right\}\) associated with internal shape function is determined from ([eq17.101]) and ([eq17.102]). These computations are \[\left\{a^{(1)}\right\}=\left[K_{a a}\right]^{-1}\left\{F I^{(1)}\right\}+\left[G_{a q}\right]\left\{q_{k}\right\}=0_{3 x 1}+\left[\begin{array}{@{}cccc@{}} 0 & \frac{L}{2 \sqrt{6}} & 0 & \frac{-L}{2 \sqrt{6}} \\[7pt] \frac{-\sqrt{6} L s_{y y}\mu}{E I_{x x}} & \frac{\sqrt{\frac{3}{2}} L^{2} \mu}{E I_{x x}} & \frac{\sqrt{6} L s_{y y} \mu}{E I_{x x}} & \frac{\sqrt{\frac{3}{2}} L^{2} \mu}{E I_{x x}} \\[7pt] \frac{L^{2} s_{y y} \mu}{\sqrt{10} E I_{x x}} & \frac{-L^{3} s_{y y} \mu}{2 \sqrt{10} E I_{x x}} & \frac{-L^{2} s_{y y} \mu}{\sqrt{10} E I_{x x}} & \frac{-L^{3} s_{y y} \mu}{2 \sqrt{10}EI_{xx}} \end{array}\right]\left[\begin{array}{@{}c@{}}0 \\ q_{2} \\ 0 \\ q_{4}\end{array}\right]. \label{eq17.3.j}\tag{j}\] Since \(q_{2}=-q_{4}\), we get \[\left\{a^{(1)}\right\}=\left[\begin{array}{@{}c@{}}\frac{-L}{\sqrt{6}} \\0 \\0\end{array}\right] q_{4}=\left[\begin{array}{@{}c@{}}\frac{-L^{2}\left(M_{a}+E I_{x x} \alpha \tau_{y}\right)}{2 \sqrt{6} E I_{x x}} \\[4pt]0 \\0\end{array}\right]. \label{eq17.3.k}\tag{k}\] The generalized displacement vector in eq. ([eq17.3.b]) is \[\left[\begin{array}{@{}l@{}}\phi_{x}^{(1)} \\v^{(1)}\end{array}\right]=\left[\begin{array}{@{}c@{}}\frac{L\left(M_{a}+E I_{x x} \alpha \tau_{y}\right) \zeta}{2 E I_{x x}} \\[4pt]\frac{L^{2}\left(M_{a}+E I_{x x} \alpha \tau_{y}\right)\left(1-\zeta^{2}\right)}{8 E I_{x x}}\end{array}\right]. \label{eq17.3.l}\tag{l}\] From ([eq17.110]) the bending moment and shear force are \[\left[\begin{array}{@{}c@{}}M_{x} \\V_{y}\end{array}\right]=\left[S_{q}\right]\left[\begin{array}{@{}c@{}}0 \\q_{2} \\0 \\q_{4}\end{array}\right]+\left[S_{F}\right]\left[\begin{array}{@{}l@{}}0 \\0 \\0\end{array}\right]-\left[\begin{array}{@{}c@{}}E I_{x x} \alpha \tau_{y} \\0\end{array}\right]=\left[\begin{array}{@{}c@{}}S_{q}\end{array}\right]\left[\begin{array}{@{}c@{}}0 \\-1 \\0 \\1\end{array}\right] q_{4}-\left[\begin{array}{@{}c@{}}E I_{x x} \alpha \tau_{y} \\0\end{array}\right]. \label{eq17.3.m}\tag{m}\] The matrix multiplication of \(\left[S_{q}\right]\) times the 4X1 vector in eq. ([eq17.3.m]) is \[\begin{aligned} &\left[\begin{array}{@{}c@{}cc@{}c@{}} -6 s_{y y}\mu\zeta & \left(\frac{-E I_{x x}}{L}+3 L s_{y y} \mu \zeta\right) & 6 s_{y y}\mu\zeta &\left(\frac{E I_{x x}}{L}+3 L s_{y y} \mu \zeta\right) \\ \frac{-12\mu s_{y y}}{L} & 6 s_{y y}\mu & \frac{12\mu s_{y y}}{L} & 6 s_{y y}\mu \end{array}\right] \left[\begin{array}{@{}c@{}}0 \\-1 \\0 \\1\end{array}\right] \nonumber\\ &\quad=\left[\begin{array}{@{}c@{}} -\left(-\frac{EI_{x x}}{L}+3\mu\zeta L s_{yy}\right)+\left(\frac{E I_{x x}}{L}+3 \mu \zeta L s_{y y}\right) \\ 0 \end{array}\right]. \label{eq17.3.n}\tag{n}\end{aligned}\] The final result for the bending moment and shear force is \[\left[\begin{array}{@{}c@{}} M_{x} \\ V_{y} \end{array}\right]=\left[\begin{array}{@{}c@{}} \frac{2 E I_{x x}}{L} \\ 0 \end{array}\right] \frac{L\left(M_{a}+E I_{x x} \alpha \tau_{y}\right)}{2 E I_{x x}}-\left[\begin{array}{@{}c@{}} E I_{x x} \alpha \tau_{y} \\ 0 \end{array}\right]=\left[\begin{array}{@{}c@{}} M_{a} \\0 \end{array}\right]. \label{eq17.3.o}\tag{o}\]

The bending moment is constant in the element and it is equal to the prescribed external moment \(M_{a}\). The shear force vanishes within the element. If \(M_{a}=0\) and \(\tau_{y} \neq 0\), then bending moment \(M_{x}=0\). However, the rotation and displacement are not zero; i.e., \[\left[\begin{array}{@{}l@{}}\phi_{x} \\v\end{array}\right]=\left.\left[\begin{array}{@{}c@{}}\left(L \alpha \tau_{y}\right) \zeta/2 \\\frac{L^{2}}{8} \alpha \tau_{y}\left(1-\zeta^{2}\right)\end{array}\right]\right|_{M_{a}=0}. \label{eq17.3.p}\tag{p}\]

This one-element solution is the same as the exact solution.

[ex17.4] A cantilever, uniform beam of length L is subject to a vertical force \(F_{a}\) at its free end as shown in figure [fig17.18]. The distributed load intensity \(f_{y}=0\), and the through-the-thickness temperature gradient \(\tau_{y}(z)=0\), \(0<z<L\). One element models the entire beam as shown in figure [fig17.17](b). Since \(f_{y}=0\) and \(\tau_{y}(z)=0\), prescribed actions ([eq17.97])

l137pt

\(\left\{F E^{(1)}\right\}=0_{4 X 1}\) and \(\left\{F I^{(1)}\right\}=0_{3 X 1}\). It follows from ([eq17.105]) that \(\left\{F^{(1)}\right\}=0_{4 X 1}\). The 4X1 generalized nodal displacement vector \(\left\{q^{(1)}\right\}\) and the 4X1 generalized force vector \(\left\{Q^{(1)}\right\}\) are partitioned into known and unknown components as follows: \[\begin{aligned} \left\{q_{\alpha}\right\} &=\left[\begin{array}{@{}ll@{}}q_{3} & q_{4}\end{array}\right]^{T}, \left\{Q_{\alpha}\right\}=\left[\begin{array}{@{}ll@{}}Q_{3} & Q_{4}\end{array}\right]^{T}=\left[\begin{array}{@{}ll@{}}F_{a} & 0\end{array}\right]^{T},\nonumber\\ \left\{q_{\beta}\right\} &=\left[\begin{array}{@{}ll@{}}q_{1} & q_{2}\end{array}\right]^{T}=0_{2 X 1}, \textrm{ and } \left\{Q_{\beta}\right\}=\left[\begin{array}{@{}ll@{}}Q_{1} & Q_{2}\end{array}\right]^{T}. \label{eq17.4.a}\tag{a}\end{aligned}\]

Matrix ([eq17.106]) represents the unrestrained structural stiffness matrix \(\left[K_{u}\right]\). The rows and columns of \(\left[K_{u}\right]\) are interchanged to the order 3, 4, 1, and 2 to facilitate partitioning it into submatrices: \[\left[K_{u}\right]=\mu \begin{array}{@{}c@{}} \begin{array}{@{}cccc@{}} q_{3} & q_{4} & q_{1} & q_{2} \\[5pt] \phantom{6 s_{y y}} & \phantom{\frac{12 E I_{x x}}{L}+4 L s_{y y}} & \phantom{-6 s_{y y}} & \phantom{\frac{-12 E I_{x x}}{L}+2 L s_{y y}} \\ \end{array} \\[-12pt] \left[\begin{array}{@{\ }cc;{2pt/2pt}cc@{\ }} \frac{12 s_{y y}}{L} & 6 s_{y y} & \frac{12 s_{y y}}{L} & 6 s_{y y} \\ 6 s_{y y} & \frac{12 E I_{x x}}{L}+4 L s_{y y} & -6 s_{y y} & \frac{-12 E I_{x x}}{L}+2 L s_{y y} \\[3pt]\hdashline & & & \\[-8pt] \frac{12^{\vphantom{M}} s_{y y}}{L} & -6 s_{y y} & \frac{12 s_{y y}}{L} & -6 s_{y y} \\6 s_{y y} & \frac{-12 E I_{x x}}{L}+2 L s_{y y} & -6 s_{y y} & \frac{12 E I_{x x}}{L}+4 L s_{y y}\end{array}\right] \end{array} =\left[\begin{array}{@{}l@{}}{\left[K_{\alpha \alpha}\right]\left[K_{\alpha \beta}\right]} \\{\left[K_{\beta \alpha}\right]\left[K_{\beta \beta}\right]}\end{array}\right]. \label{eq17.4.b}\tag{b}\] The matrix equation to determine the unknown displacements is \(\left[K_{\alpha \alpha}\right]\left\{q_{\alpha}\right\}+\left[K_{\alpha\beta}\right]\left\{q_{\beta}\right\}=\left\{Q_{\alpha}\right\}\), or \[\frac{E I_{x x}}{12 E I_{x x}+L^{2} s_{y y}}\left[\begin{array}{@{}cc@{}}\frac{12 s_{y y}}{L} & 6 s_{y y} \\[4pt] 6 s_{y y} & \frac{12 E I_{x x}}{L}+4 L s_{y y}\end{array}\right]\left[\begin{array}{@{}l@{}}q_{3} \\q_{4}\end{array}\right]=\left[\begin{array}{@{}c@{}}F_{a} \\0\end{array}\right]. \label{eq17.4.c}\tag{c}\] The solution of eq. ([eq17.4.c]) for the displacements is \[\left[\begin{array}{@{}l@{}}q_{3} \\q_{4}\end{array}\right]=\left[\begin{array}{@{}c@{}}\frac{F_{a} L^{3}}{3 E I_{x x}}+\frac{F_{a} L}{s_{y y}} \\[4pt] \frac{-F_{a} L^{2}}{2 E I_{x x}}\end{array}\right]. \label{eq17.4.d}\tag{d}\] The matrix formulation to determine the generalized reactive force vector is \(\left[K_{\beta \alpha}\right]\left\{q_{\alpha}\right\}+\left[K_{\beta \beta}\right]\left\{q_{\beta}\right\}=\left\{Q_{\beta}\right\}\), or \[\left[\begin{array}{@{}l@{}}Q_{1} \\Q_{2}\end{array}\right]=\mu\left[\begin{array}{@{}cc@{}}\frac{-12 s_{y y}}{L} & -6 s_{y y} \\6 s_{y y} & \frac{-12 E I_{x x}}{L}+2 L s_{y y}\end{array}\right]\left[\begin{array}{@{}c@{}}\frac{L^{2}}{3 E I_{x x}}+\frac{1}{s_{y y}} \\[4pt] \frac{-L}{2 E I_{x x}}\end{array}\right] F_{a} L=\mu L F_{a}\left[\begin{array}{@{}c@{}}\frac{-\left(12 E I_{x x}+L^{2} s_{y y}\right)}{E I_{x x} L} \\[4pt] \frac{12 E I_{x x}+L^{2} s_{y y}}{E I_{x x}}\end{array}\right]=\left[\begin{array}{@{}l@{}}-F_{a} \\L F_{a}\end{array}\right]. \label{eq17.4.e}\tag{e}\] The generalized displacement vector \(\left\{a^{(1)}\right\}\) associated with internal shape function is determined from ([eq17.101]) and ([eq17.102]). The result is \[\left\{a^{(1)}\right\}=\left[K_{a a}\right]^{-1}\left\{F I^{(1)}\right\}+\left[G_{a q}\right]\left\{q_{k}\right\}=0_{3 X 1}+\left[G_{a q}\right] \left[\begin{array}{@{}l@{}}0 \\0 \\q_{3} \\q_{4}\end{array}\right]= \left[\begin{array}{@{}c@{}}\frac{F_{a} L^{3}}{4 \sqrt{6} E I_{x x}} \\[6pt]\frac{F_{a} L^{2}}{2 \sqrt{6} E I_{x x}} \\[6pt]\frac{-F_{a} L^{3}}{12 \sqrt{10} E I_{x x}}\end{array}\right]. \label{eq17.4.f}\tag{f}\] The generalized displacement vector ([eq17.70]) is \[\left[\begin{array}{@{}l@{}}\phi_{x}^{(1)} \\[4pt]v^{(1)}\end{array}\right]=\left[\begin{array}{@{}cccc@{}}0 & \eta_{1} & 0 & \eta_{2} \\\eta_{1} & 0 & \eta_{2} & 0\end{array}\right]\left[\begin{array}{@{}c@{}}0 \\0 \\q_{3} \\ q_{4} \end{array}\right]+\left[\begin{array}{@{}ccc@{}} 0 & \eta_{3} & 0 \\ \eta_{3} & 0 & \eta_{4} \end{array}\right]\left[\begin{array}{@{}l@{}} a_{1}^{(1)} \\[4pt] a_{2}^{(1)} \\[4pt] a_{3}^{(1)} \end{array}\right]=\left[\begin{array}{@{}c@{}} \frac{F_{a} L^{2}}{8 E I_{x x}}\left(-3-2 \zeta+\zeta^{2}\right) \\[6pt] \frac{F_{a} L(1+\zeta)\left[24 E I_{x x}+L^{2} s_{y y}\left(5+4 \zeta-\zeta^{2}\right)\right]}{48 E I_{x x} s_{y y}} \end{array}\right]. \label{eq17.4.g}\tag{g}\] The vector of the bending moment and shear force ([eq17.110]) is \[\left[\begin{array}{@{}c@{}}M_{x} \\V_{y}\end{array}\right]=\left[S_{q}\right]\left[\begin{array}{@{}c@{}}0 \\0 \\q_{3} \\q_{4}\end{array}\right]=\left[\begin{array}{@{}c@{}}\left(6 \mu s_{y y} \zeta\right)\left(\frac{F_{a} L^{3}}{3 E I_{x x}}+\frac{F_{a} L}{s_{y y}}\right)+\left(\frac{E I_{x x}}{L}+3 \mu L_{k} s_{y y} \zeta\right)\left(\frac{-F_{a} L^{2}}{2 E I_{x x}}\right) \\[6pt] \left(12 \mu s_{y y}\right)\left(\frac{F_{a} L^{3}}{3 E I_{x x}}+\frac{F_{a} L}{s_{y y}}\right)+\left(6 \mu s_{y y}\right)\left(\frac{-F_{a} L^{2}}{2 E I_{x x}}\right)\end{array}\right]=\left[\begin{array}{@{}c@{}}\frac{F_{a} L(-1+\zeta)}{2} \\[4pt]F_{a}\end{array}\right]. \label{eq17.4.h}\tag{h}\] The shear force is constant in the beam and is equal to the prescribed external force \(F_{\textrm{a}}\). For \(F_{a}>0\) the bending moment decreases linearly from zero at the free end to a minimum (\(-F_{a} L\)) at the clamped end. The finite element solution for transverse bending is the same as the exact solution of the governing boundary value problem.

[ex17.5]The wing spar of a light airplane described in example [ex6.6] on page  is modeled as a cantilever beam as shown in figure [fig17.19](a). In a symmetric maneuver of the airplane the total lift \(L=12{,}000\) lb., and the lift acting on each wing is L/2. Assume the airload acts along the locus of shear centers so that spar bends without twist in torsion. The airload is distributed elliptically over the wing, so that the airload intensity \(f_{y}(z)\) per unit span is given as \[f_{y}(z)=\frac{2 L}{\pi z_{m a x}} \sqrt{1-\left(\frac{z}{z_{m a x}}\right)^{2}} \quad 0 \leq z \leq z_{\max}, \label{eq17.5.a}\tag{a}\] where \(z\) is the spanwise coordinate, \(z=0\) at the root, and \(z=z_{\max}=120\,\mathrm{in}\) at the tip of the wing spar. The transverse temperature gradient \(\tau_{y}(z)=0\), \(0<z<z_{\max}\). Data taken from example [ex6.6] are: the modulus of elasticity \(E=10.5 \times 10^{6}\,\mathrm{lb}./\text { in. }^{2}\), the second area moment of the cross section about the \(x\)-axis \(I_{x x}=101.619\,\mathrm{in}^{4}\), and the transverse shear coefficient \(s_{y y}=2.4278 \times 10^{6}\,\mathrm{lb}\)

The finite element model shown in figure [fig17.19](b) has three nodes \(\left\{z_{3}\right\}=\{0,60,120\} \text { in. }\), and two equal length elements \(h_{1}=h_{2}=60\,\mathrm{in}\). The stiffness matrices ([eq17.106]) for each element are \[\begin{gathered} \left[K^{(1)}\right]=\left[\begin{array}{@{}cccc@{}}q_{1} & q_{2} & q_{3} & q_{4} \\24,048 & -721,441 & -24,048 & -721,441 \\-721,441 & 3.94266 \times 10^{7} & 721,441 & 3.85991 \times 10^{6} \\-24,048 & 721,441 & 24,048 & 721,441 \\-721,441 & 3.85991 \times 10^{6} & 721,441 & 3.94266 \times 10^{7}\end{array}\right] \textrm{ and}\label{eq17.5.b}\tag{b}\\ \left[K^{(2)}\right]=\left[\begin{array}{@{}cccc@{}}q_{3} & q_{4} & q_{4} & q_{5} \\24,048 & -721,441 & -24,048 & -721,441 \\-721,441 & 3.94266 \times 10^{7} & 721,441 & 3.85991 \times 10^{6} \\-24,048 & 721,441 & 24,048 & 721,441 \\-721,441 & 3.85991 \times 10^{6} & 721,441 & 3.94266 \times 10^{7}\end{array}\right]. \label{eq17.5.c}\tag{c}\end{gathered}\] The assembly process effects the matrix addition \(\left[K_{u}\right]=\left[K^{(1)}\right]+\left[K^{(2)}\right]\) to get the unrestrained structural stiffness matrix as \[\left[K_{u u}\right]=\left[\begin{array}{@{}cccccc@{}}q_{1} & q_{2} & q_{3} & q_{4} & q_{5} & q_{6} \\24,048 & -721,441 & -24,048 & -721441 . & 0 & 0 \\-721,441 & 3.94266 \times 10^{7} & 721441 & 3.85991 \times 10^{6} & 0 & 0 \\[3pt] \hdashline\\[-8pt] -24,048 & 721,441 & 48,096.1 & 0 & -24,048 & -721441 . \\[3pt]\hdashline\\[-8pt] 721,441 & 3.85991 \times 10^{6} & 0 & 7.88531 \times 10^{7} & 721,441 & 3.85991 \times 10^{6} \\0 & 0 & -24,048 & 721,441 & 24,048 & 721,441 \\0 & 0 & -721,441 & 3.85991 \times 10^{6} & 721,441 & 3.94266 \times 10^{7}\end{array}\right]. \label{eq17.5.d}\tag{d}\] The generalized nodal displacement vector \(\left[\begin{array}{@{}llllll@{}}q_{1} & q_{2} & q_{3} & q_{4} & q_{5} & q_{6}\end{array}\right]^{T}\) is partitioned into unknown components \(\left\{q_{\alpha}\right\}\) and known components \(\left\{q_{\beta}\right\}\): \[\left\{q_{\alpha}\right\}=\left[\begin{array}{@{}llll@{}}q_{3} & q_{4} & q_{5} & q_{6}\end{array}\right]^{T} \quad\left\{q_{\beta}\right\}=\left[\begin{array}{@{}ll@{}}q_{1} & q_{2}\end{array}\right]^{T}=0_{2 X 1}. \label{eq17.5.e}\tag{e}\] We partition the unrestrained structural stiffness in eq. ([eq17.5.d]) in the form \[\begin{aligned} \left[K_{u}\right]=\left[\begin{array}{@{}ll@{}}{\left[K_{\beta \beta}\right]} & {\left[K_{\beta \alpha}\right]} \\[4pt]{\left[K_{\alpha \beta}\right]} & {\left[K_{\alpha \alpha}\right]}\end{array}\right], \textrm{ where}\label{eq17.5.f}\tag{f}\end{aligned}\] \[\begin{aligned} \left[K_{\alpha \alpha}\right]&= \left[\begin{array}{@{}cccc@{}} q_{3} & q_{4} & q_{5} & q_{6} \\[4pt]48,096.1 & 0 & -24,048 & -721,441 \\[4pt]0 & 7.88531 \times 10^{7} & 721,441 & 3.85991 \times 10^{6} \\[4pt]-24,048 & 721,441 & 24,048 & 721,441 \\[4pt]-721,441 & 3.85991 \times 10^{5} & 721,441 & 3.94266 \times 10^{7}\end{array}\right],\nonumber\\ \left[K_{\alpha \beta}\right] &= \begin{array}{@{}c@{}} \begin{array}{@{}c@{\qquad}c@{}} q_{1} & q_{2} \end{array}\\[4pt] \left[\begin{array}{@{}cc@{}} -24,048 & 721,441 \\[4pt] -721,441 & 3.85991 \times 10^{6} \\[4pt] 0 & 0 \\[4pt] 0 & 0 \end{array}\right]\\ \end{array} =\left[K_{\beta \alpha}\right]^{T}\!, \label{eq17.5.g}\tag{g}\end{aligned}\] and \[\left[K_{\beta \beta}\right]=\left[\begin{array}{@{}cc@{}}q_{1} & q_{2} \\ 24,048 & -721,441 \\-721,441 & 3.94266 \times 10^{7}\end{array}\right]. \label{eq17.5.h}\tag{h}\] The generalized forces acting at the nodes of elements 1 and 2 are shown in figure [fig17.20]. The external generalized forces acting on elements \(\Omega_1\) and \(\Omega_2\) from the distributed airload ([eq17.5.a]) are computed by numerical integration because of the complexity of the integrands ([eq17.97]). (Refer to article 1.3.6 for details on numerical integration.) The external generalized forces acting on each element are determined from eq. ([eq17.105]). The transpose of the Guyan matrix used to compute the generalized force vectors \(\left\{F^{(1)}\right\}\) and \(\left\{F^{(2)}\right\}\) in eq. ([eq17.105]) is the same for each element in this example. It is given by \[\left[G_{a q}^{(1)}\right]^{T}=\left[G_{a q}^{(2)}\right]^{T}=\left[\begin{array}{@{}ccc@{}} 0 & -0.016526 & 0.128289 \\12.2474 & 0.496859 & -3.84866 \\ 0 & 0.016562 & -0.128289 \\ -12.2474 & 0.496859 & -3.84866\end{array}\right]. \label{eq17.5.i}\tag{i}\] The results for the generalized forces are \[\begin{gathered} \left\{F E^{(1)}\right\} \equiv\left[\begin{array}{@{}c@{}}1,869\,\mathrm{lb}. \\0 \\1,784.98\,\mathrm{lb}. \\0\end{array}\right] \quad\left\{F I^{(1)}\right\} \equiv\left[\begin{array}{@{}c@{}}-1,499.03\,\mathrm{lb}. \\0 \\266.3992\,\mathrm{lb}.\end{array}\right] \quad\left\{F^{(1)}\right\}=\left[\begin{array}{@{}c@{}}F_{1}^{(1)} \\F_{2}^{(1)} \\F_{3}^{(1)} \\F_{4}^{(1)}\end{array}\right]=\left[\begin{array}{@{}c@{}}1,872.39\,\mathrm{lb}. \\-18,460.9\,\mathrm{lb}.\mathrm{-in}. \\1,781.6\,\mathrm{lb}. \\18,257.7\,\mathrm{lb}.\mathrm{-in}.\end{array}\right], \label{eq17.5.j}\tag{j} \\ \left\{F E^{(2)}\right\} \equiv\left[\begin{array}{@{}c@{}}1,384.05\,\mathrm{lb}. \\0 \\961.96\,\mathrm{lb}. \\0\end{array}\right] \quad\left\{F I^{(2)}\right\} \equiv\left[\begin{array}{@{}c@{}}-991.799\,\mathrm{lb}. \\0 \\125.595\,\mathrm{lb}.\end{array}\right], \textrm{ and } \left\{F^{(2)}\right\}=\left[\begin{array}{@{}c@{}}F_{3}^{(2)} \\F_{4}^{(2)} \\F_{5}^{(2)} \\F_{6}^{(2)}\end{array}\right]=\left[\begin{array}{@{}c@{}}1,400.17\,\mathrm{lb}. \\-12,630.4\,\mathrm{lb}\mathrm{-in}. \\945.848\,\mathrm{lb}. \\11,663.6\,\mathrm{lb}.\mathrm{-in}.\end{array}\right]. \label{eq17.5.k}\tag{k}\end{gathered}\] Assembly of the two elements is shown in figure [fig17.19](b). The external generalized force vectors for the assembly are \[\left[\begin{array}{@{}c@{}}F_{1} \\ F_{2} \\ F_{3} \\ F_{4} \\ F_{5} \\ F_{6}\end{array}\right]=\left[\begin{array}{@{}c@{}}F_{1}^{1)} \\ F_{2}^{(1)} \\ F_{3}^{(1)}+F_{3}^{(2)} \\ F_{4}^{(1)}+F_{4}^{(2)} \\ F_{5}^{(2)} \\ F_{6}^{(2)}\end{array}\right]=\left[\begin{array}{@{}c@{}}1,872.39\,\mathrm{lb}. \\ -18,460.9\,\mathrm{lb}.\,\mathrm{in}. \\ 3,181.77\,\mathrm{lb}. \\ 5,627.30\,\mathrm{lb}.\,\mathrm{in}. \\ 945.848\,\mathrm{lb}. \\ 11,663.6\,\mathrm{lb}.\,\mathrm{in}.\end{array}\right]\left[\begin{array}{@{}c@{}}Q_{1} \\ Q_{2} \\ Q_{3} \\ Q_{4} \\ Q_{5} \\ Q_{6}\end{array}\right]=\left[\begin{array}{@{}c@{}}Q_{1}^{(1)} \\[4pt] Q_{2}^{(2)} \\[4pt] Q_{3}^{(1)}+Q_{3}^{(2)} \\[4pt] Q_{4}^{(1)}+Q_{4}^{(2)} \\[4pt] Q_{5}^{(2)} \\[4pt] Q_{6}^{(2)}\end{array}\right]. \label{eq17.5.l}\tag{l}\] The 6X1 generalized nodal force vector \(\left[\begin{array}{@{}llllll@{}}Q_{1} & Q_{2} & Q_{3} & Q_{4} & Q_{5} & Q_{6}\end{array}\right]^{T}\) is partitioned into known components \(\left\{Q_{\alpha}\right\}\) and unknown components \(\left\{Q_{\beta}\right\}\) as follows: \[\left\{Q_{\alpha}\right\}=\left[\begin{array}{@{}llll@{}} Q_{3} & Q_{4} & Q_{5} & Q_{6}\end{array}\right]^{T}=0_{4 X 1}, \textrm{ and } \left\{Q_{\beta}\right\}=\left[Q_{1}\ \ Q_{2}\right]^{T}. \label{eq17.5.m}\tag{m}\] There are no generalized point forces acting at nodes 2 and 3. Generalized point forces acting at node 1 are reactive. The generalized forces from the airload are also partitioned as \[\begin{gathered} \left\{F_{\alpha}\right\}=\left[\begin{array}{@{}llll@{}}F_{3} & F_{4} & F_{5} & F_{6}\end{array}\right]^{T}=\left[\begin{array}{@{}llll@{}}3,181.76 & 5,627.37 & 945.848 & 11,663.6\end{array}\right]^{T}, \textrm{ and}\label{eq17.5.n}\tag{n}\\ \left\{F_{\beta}\right\}=\left[\begin{array}{@{}ll@{}}F_{1} & F_{2}\end{array}\right]^{T}=[1,872.39-18,460.9]^{T}. \label{eq17.5.o}\tag{o}\end{gathered}\]

The matrix equation to determine the unknown generalized displacements is \(\left[K_{\alpha \alpha}\right]\left\{q_{\alpha}\right\}+\left[K_{\alpha \beta}\right]\left\{q_{\beta}\right\}=\left\{Q_{\alpha}\right\}+\left\{F_{\alpha}\right\}\).The solution for the unknown displacement vector is \[\left\{q_{\alpha}\right\}=\left[\begin{array}{@{}llll@{}}0.447103 & -0.0091821 & 1.06555 & -0.0101218\end{array}\right]^{T}. \label{eq17.5.p}\tag{p}\] The matrix equation to determine the unknown forces is \(\left[K_{\beta \alpha}\right]\left\{q_{\alpha}\right\}+\left[K_{\beta \beta}\right]\left\{q_{\beta}\right\}=\left\{Q_{\beta}\right\}+\left\{F_{\beta}\right\}\), or \[\left[\begin{array}{@{}cccc@{}}-24,048 & -721,441 & 0 & 0 \\721,441 & 3.85991 \times 10^{6} & 0 & 0\end{array}\right]\left[\begin{array}{@{}c@{}}0.447103 \\-0.0091821 \\1.06555 \\-0.0101218\end{array}\right]+\left[\begin{array}{@{}cc@{}}24,048 & -721,441 \\-721,441 & 3.94266 \times 10^{7}\end{array}\right]\left[\begin{array}{@{}l@{}}0 \\0\end{array}\right]=\left[\begin{array}{@{}l@{}}Q_{1} \\Q_{2}\end{array}\right]+\left[\begin{array}{@{}c@{}}1,872.39 \\-18,460.9\end{array}\right]. \label{eq17.5.q}\tag{q}\] Perform the matrix algebra in eq. ([eq17.5.q]) to find that the generalized reactive forces are \[Q_{1}=-6{,}000.0\,\mathrm{lb}. \quad Q_{2}=305,577\,\mathrm{lb} \cdot-\mathrm{in}. \label{eq17.5.r}\tag{r}\] The external generalized forces acting on the two-element model of wing spar are shown in the free body diagram in figure [fig17.21].

Gaussian integration

Consider the integral to compute the external force component \(F E_{2k+1}^{(k)}\) in eq. ([eq17.97]) for element two in example [ex17.5]. The expression for this force component is \[\begin{aligned} \label{eq17.116} F E_{5}^{(2)}=\frac{h_{2}}{2} \int_{-1}^{1} \frac{1}{2}(1+\zeta) \frac{2 L}{\pi z_{max }} \sqrt{1-\left(\frac{z^{(2)}(\zeta)}{z_{max }}\right)^{2}} d \zeta,\end{aligned}\] where \[\begin{aligned} \label{eq17.117} z^{(2)}=\eta_{1}(\zeta)(60 .)+\eta_{2}(\zeta)(120)=90+30 \zeta.\end{aligned}\] Numerical evaluation of eq. ([eq17.116]) leads to \[\begin{aligned} \label{eq17.118} F E_{5}^{(2)}=\int_{-1}^{1}\left[954.930(1+\zeta) \sqrt{0.4357-0.375 \zeta-0.0625 \zeta^{2}}\right] d \zeta.\end{aligned}\] We carry out Gaussian integration of eq. ([eq17.118]) after a discussion of the method.

The method of Gaussian quadrature is to approximate the integral \(I=\int_{-1}^{1} f(\zeta) d \zeta\) by \[\begin{aligned} \label{eq17.119} I \approx I_{\text {appr }}=\sum_{i=1}^{n s} w_{i} f\left(\zeta_{j i}\right),\end{aligned}\] where \(\zeta_{i}\) are the abscissas of the Legendre polynomial \(P_{i}(\zeta)\), and \(w_{i}\) are the weight factors. The abscissas, or roots, of the Legendre polynomial are symmetrically located in the interval \(-1 \leq \zeta \leq 1\). The weight factors are determined such that a polynomial of degree \(p\) is exactly equal to the sum in ([eq17.119]). Consider a polynomial of second degree \(f(\zeta)=c_{0}+c_{1} \zeta+c_{2} \zeta^{2}\). The exact integral is \[\begin{aligned} \label{eq17.120} I=\int_{-1}^{1}\left(c_{0}+c_{1} \zeta+c_{2} \zeta^{2}\right) d \zeta=2 c_{0}+\frac{2}{3} c_{2}.\end{aligned}\] The exact integral is determined by two coefficients, \(c_0\) and \(c_2\). The Legendre polynomial for \(n=2\) is \(P_{2}(\zeta)=\left(-1+3 \zeta^{2}\right)/2\). The roots of \(P_{2}(\zeta)=0\) are \[\begin{aligned} \label{eq17.121} -1/\sqrt{3} \textrm{ and } 1 /(\sqrt{3}).\end{aligned}\] Equation ([eq17.119]) for \(n=2\) is \[\begin{aligned} \label{eq17.122} I=I_{\text {appr. }}=w_{1} f\left(\frac{-1}{\sqrt{3}}\right)+w_{2} f\left(\frac{1}{\sqrt{3}}\right).\end{aligned}\] Evaluating eq. ([eq17.122]) we get \[\begin{aligned} \label{eq17.123} I=I_{\text {appr. }}=2\left(w_{1}+w_{2}\right) c_{0}+\frac{\left(w_{2}-w_{1}\right)}{\sqrt{3}} c_{1}+\frac{1}{3}\left(w_{1}+w_{2}\right).\end{aligned}\] Equate like coefficients in ([eq17.120]) and ([eq17.123]) to get three linear equations for the weight factors \(w\)1, and \(w\)2: \[\begin{aligned} \label{eq17.124} 2=w_{1}+w_{2} \quad 0=\left(w_{2}-w_{1}\right) /(\sqrt{3}) \quad \frac{2}{3}=\frac{1}{3}\left(w_{1}+w_{2}\right).\end{aligned}\] The solution of eq. ([eq17.124]) is \(w_{1}=w_{2}=1\). Note that the weight factors of the positive root and the corresponding negative root are the same because the term with the odd power of \(\zeta\) does not contribute to I ([eq17.120]). Thus, a polynomial of degree 2 can be integrated exactly by evaluating it at \(\mp 1/\sqrt{3}\) and multiplying by the appropriate equal weight factors.

Consider the cubic polynomial \(f(\zeta)=c_{0}+c_{1} \zeta+c_{2} \zeta^{2}+c_{3} \zeta^{3}\) whose exact integral is \(I=2 c_{0}+(2/3) c_{2}\). The Legendre polynomial for \(n\) = 3 is \(P_{3}(\zeta)=\left(-3 \zeta+5 \zeta^{3}\right)/2\). The roots of \(P_{3}(\zeta)=0\) are \(0\), and \(\mp \sqrt{3/5}\). Equation ([eq17.119]) for \(n=3\) is \[\begin{aligned} \label{eq17.125} I=I_{\text {appr. }}=w_{1} f(0)+w_{2} f(-\sqrt{3/5})+w_{2} f(\sqrt{3/5})=c_{0}\left(w_{1}+w_{2}\right)+c_{2}\left(6 w_{2}/ 5\right).\end{aligned}\] Equating like coefficients between the last equation and the exact integral we find the weight factors \(w_{1}=8/9\) and \(w_{2}=5/9\). In general, a polynomial of degree \(p\) is integrated exactly if \(n \geq(p+1)/2\).

Now consider the integrand in the definite integral ([eq17.118]). The integrand is \[\begin{aligned} \label{eq17.126} f(\zeta)=954.930(1+\zeta) \sqrt{0.4357-0.375 \zeta-0.0625 \zeta^{2}}.\end{aligned}\] Function \(f(\zeta)\) is continuous in the interval \(|\zeta| \leq 1\) as shown in figure [fig17.22]. Its Taylor’s formula is an infinite series \(f(\zeta)=\sum_{m=0}^{\infty} a_{m} \zeta^{m}\), where the \(a_{m}\) are real numbers. Consequently if \(f(\zeta)\) is replaced by its Taylor series, then the integrand is a polynomial of infinite degree. The integral of a polynomial of infinite degree implies an infinite number of abscissas and weights for Gaussian quadrature, which, of course, is not practical. Only a finite number of abscissas and weights in Gaussian quadrature are considered in the sequence of numerical integrations to follow.

Gaussian quadrature of the integrand ([eq17.119]) for \(n=2\) is \[\begin{aligned} \label{eq17.127} I_{\text {appr. }}=w_{1} f(-1/\sqrt{3})+w_{1} f(1/\sqrt{3})=(1)(321.154)+(1)(673.889)=995.043\,\mathrm{lb}\end{aligned}\] Gaussian quadrature of the integrand ([eq17.119]) for \(n=3\) is \[\begin{aligned} \label{eq17.128} I_{\text {appr. }}=w_{1} f(0)+w_{2} f\left(-\sqrt{\frac{3}{5}}\right)+w_{2} f\left(\sqrt{\frac{3}{5}}\right)=\left(\frac{8}{9}\right)(631.627)+\left(\frac{5}{9}\right)(178.857)+\left(\frac{5}{9}\right)(560.829)=972.382\,\mathrm{lb}\end{aligned}\] For \(n=2\), 3, 4, 5, 6, 7, 8, 9, and 10 the abscissa and weight factors are listed in table [tab17.4]. Approximate values of the integrals of ([eq17.126]) are listed in table [tab17.5] based on the data in table [tab17.4]. The results in table [tab17.5] show the integrals are slowly decreasing as the number of terms in the Gaussian integration increase. Apparently, the approximate values of the integrals are asymptotically approaching the value of 961.96 lb., which is the value of the integral of ([eq17.126]) computed from the function NIntegrate[\(f(\zeta), \{\zeta, -1, 1\}\)] in Mathematica.

Element displacement functions and strains

The kth element is denoted by \(\Omega_{k} \equiv\left\{z \mid z_{k} \leq z \leq z_{k+1}\right\}\), where \(Z_{k}<Z_{k+1}\). The standard element \(\Omega_{\text {st }} \equiv\{\zeta \mid-1<{\zeta<1}\}\) ([eq17.15]) is mapped to the kth element by ([eq17.16]), and the inverse mapping is given by ([eq17.17]). The lateral displacement of the kth element is denoted by \(v^{(k)}(z)\) and the clockwise rotation by \(-\left(\frac{d v^{(k)}}{d z}\right)\). Define the generalized nodal displacements as \[\begin{aligned} \label{eq17.139} v^{(k)}\left(z_{k}\right)=q_{2 k-1}\quad -\!\left.\left(\frac{d v^{(k)}}{d z}\right)\right|_{z=z_{k}}=q_{2 k} \quad v^{(k)}\left(z_{k+1}\right)=q_{2 k+1} \quad-\!\left.\left(\frac{d v^{(k)}}{d z}\right)\right|_{z=z_{k+1}}=q_{2 k+2}.\end{aligned}\]

r106pt

See figure [fig17.23]. Admissible functions \(v(z)\) and its derivative must be continuous within an element and be continuous between elements. Interpolation functions that satisfy these continuity requirements are Hermite cubic interpolation functions, which are denoted by \(\phi_{i}(\zeta)\),. \(i=1,2,3,4\). These interpolation functions are \[\begin{aligned} \label{eq17.140} \begin{array}{@{}ll@{}}\phi_{1}(\zeta)=\left(2-3 \zeta+\zeta^{3}\right)/4 & \phi_{2}(\zeta)=\left(h_{k}/ 8\right)\left(-1+\zeta+\zeta^{2}-\zeta^{3}\right) \\[4pt]\phi_{3}(\zeta)=\left(2+3 \zeta-\zeta^{3}\right)/4 & \phi_{4}(\zeta)=\left(h_{k}/8\right)\left(1+\zeta-\zeta^{2}-\zeta^{3}\right)\end{array}.\end{aligned}\] Let the derivative \(\frac{d \phi_{i}}{d \zeta}\) be denoted by \(\phi_{i}{ }^{\prime}\). The interpolation properties of the Hermite cubic functions are \[\begin{aligned} \label{eq17.141} \begin{array}{@{}cccc@{}}\phi_{1}(-1)=1 & \phi_{1}^{\prime}(-1)=0 & \phi_{1}(1)=0 & \phi_{1}^{\prime}(1)=0 \\[4pt]\phi_{2}(-1)=0 & \phi_{2}^{\prime}(-1)=-h_{k}/2 & \phi_{2}(1)=0 & \phi_{2}^{\prime}(1)=0 \\[4pt]\phi_{3}(-1)=0 & \phi_{3}^{\prime}(-1)=0 & \phi_{3}(1)=1 & \phi_{3}^{\prime}(1)=0 \\[4pt]\phi_{4}(-1)=0 & \phi_{4}^{\prime}(-1)=0 & \phi_{4}(1)=0 & \phi_{4}^{\prime}(1)=-h_{k}/2\end{array}.\end{aligned}\]

Graphs of Hermite cubic functions are shown in figure [fig17.24]. The displacement for element \(\Omega_{\textrm{k}}\) is \(v^{(k)}(\zeta)=[N(\zeta)]\left\{q^{(k)}\right\}\), where \([N(\zeta)]\) is the 1X4 matrix of interpolation functions, and \(\left\{q^{(k)}\right\}\) is the 4X1 vector of the generalized displacements at the nodes. The interpolation matrix and nodal displacement vector are \[\begin{aligned} \label{eq17.142} [N(\zeta)]=\left[\phi_{1}(\zeta) \phi_{2}(\zeta) \phi_{3}(\zeta) \phi_{4}(\zeta)\right], \textrm{ and } \left\{q^{(k)}\right\}=\left[\begin{array}{@{}llll@{}} q_{2 k-1} & q_{2 k} & q_{2 k+1} & q_{2 k+2}\end{array}\right]^{T}.\end{aligned}\] The virtual displacement is \(\bar{v}(\zeta)=[N(\zeta)]\{b\}\), where \(\{b\}=\left[\begin{array}{@{}llll@{}}b_{1} & b_{2} & b_{3} & b_{4}\end{array}\right]^{T}\). The second derivative of the displacement, or curvature \(\kappa\) of the centroidal axis in bending, is expressed as \[\begin{gathered} \kappa=-\left(\frac{d^{2} v}{d z^{2}}\right)=-\left(\frac{d^{2} v}{d \zeta^{2}}\right)\left(\frac{2}{h_{k}}\right)^{2}=-\left(\frac{2}{h_{k}}\right)^{2} \frac{d^{2}}{d \zeta^{2}}[N(\zeta)]\left\{q^{(k)}\right\}=\left[N_{\kappa}(\zeta)\right]\big\{q^{(k)}\big\}, \textrm{ where}\label{eq17.143}\\ \left[N_{\kappa}(\zeta)\right]=\frac{1}{h_{k}^{2}}\left[\begin{array}{@{}llll@{}} -6 \zeta & h_{k}(-1+3 \zeta) & 6 \zeta & h_{k}(1+3 \zeta)\end{array}\right]. %\left[-6 \zeta h_{k}(-1+3 \zeta) 6 \zeta h_{k}(1+3 \zeta)\right]. \label{eq17.144}\end{gathered}\] The the virtual curvature is \(\bar{\kappa}=-\left(\frac{d^{2} \bar{v}^{(k)}}{d z^{2}}\right)=\left[N_{\kappa}(\zeta)\right]\{b\}=\{b\}^{T}\left[N_{\kappa}(\zeta)\right]^{T}\). The bilinear form evaluates as \[\begin{aligned} \label{eq17.145} B[v, \bar{v}]=\int_{0}^{L}\left(\bar{\kappa} E I_{x x} \kappa\right) d z=\{b\}^{T}\left(\int_{-1}^{1}\left[N_{\kappa}(\zeta)\right]^{T} E I_{x x}\left[N_{\kappa}(\zeta)\right] \frac{h_{k}}{2} d \zeta\right)\big\{q^{(k)}\big\}=\{b\}^{T}[K]\big\{q^{(k)}\big\}.\end{aligned}\] The element stiffness matrix is \[\begin{aligned} \label{eq17.146} [K]=\frac{E I_{x x}}{h_{k}^{3}}\left[\begin{array}{@{}cccc@{}}12 & -6 h_{k} & -12 & -6 h_{k} \\-6 h_{k} & 4 h_{k}^{2} & 6 h_{k} & 2 h_{k}^{2} \\-12 & 6 h_{k} & 12 & 6 h_{k} \\-6 h_{k} & 2 h_{k}^{2} & 6 h_{k} & 4 h_{k}^{2}\end{array}\right].\end{aligned}\]

Referring to Fig. [fig17.16] on page , we express the boundary terms in the external virtual work ([eq17.138]) as \[\begin{aligned} \left.\bar{v} V_{y}\right|_{z_{k+1}}-\left.\bar{v} V_{y}\right|_{z_{k}}+\left.\left(-\frac{d \bar{v}}{d z}\right) M_{x}\right|_{z_{k+1}}-\left.\left(-\frac{d \bar{v}}{d z}\right) M_{x}\right|_{z_{k}}&=b_{3} Q_{2 k+1}-b_{1}\left(-Q_{2 k-1}\right)+b_{4} Q_{2 k+2}-b_{2}\left(-Q_{2 k}\right) \nonumber \\ &=\{b\}^{T}\big\{Q^{(k)}\big\}. \label{eq17.147}\end{aligned}\] The prescribed distributed load terms in the external virtual work are \[\begin{aligned} \label{eq17.148} \int_{0}^{L}\hspace*{-2.5pt} f_{y} \bar{v} d z+\int_{0}^{L}\left[E I_{x x}\left(-\frac{d^{2} \bar{v}}{d z^{2}}\right)\right] \alpha \tau_{y} d z=\{b\}^{T}\left(\int_{-1}^{1} f_{y}^{(k)}[N(\zeta)]^{T} \frac{h_{k}}{2} d \zeta+\int_{-1}^{1} E I_{x x}\left[N_{k}(\zeta)\right]^{T} \alpha \tau_{y}^{(k)}(\zeta) \frac{h_{k}}{2} d \zeta\right),\end{aligned}\] where the distributed loading in the element is \[\begin{aligned} \label{eq17.149} f_{y}^{(k)}=f_{y}\left[\eta_{1}(\zeta) z_{k}+\eta_{2}(\zeta) z_{k+1}\right] \quad \tau_{y}^{(k)}=\tau_{y}\left[\eta_{1}(\zeta) z_{k}+\eta_{2}(\zeta) z_{k+1}\right].\end{aligned}\] The final result for the external virtual work for element \(\Omega_{\textrm{k}}\) is \[\begin{gathered} F[\bar{v}]=\{b\}^{T}\left(\left\{Q^{k}\right\}+\big\{F^{(k)}\big\}\right), \textrm{ where}\label{eq17.150}\\ \big\{Q^{(k)}\big\}=\left[\begin{array}{@{}llll@{}}Q_{2 k-1} & Q_{2 k} & Q_{2 k+1} & Q_{2 k+2}\end{array}\right]^{T} \textrm{ and } \left\{F^{(k)}\right\}=\left[\begin{array}{@{}llll@{}}F_{2 k-1} &F_{2 k} &F_{2 k+1} &F_{2 k+2}\end{array}\right]^{T}. \label{eq17.151}\end{gathered}\] Explicit expressions for the components of the generalized force vector from the distributed load are \[\begin{gathered} F_{2 k-1}=\int_{-1}^{1}\left[f_{y}^{(k)}(\zeta)\right] \phi_{1}(\zeta) \frac{h_{k}}{2} d \zeta+\int_{-1}^{1}\left[E I_{x x}\left(\frac{-6 \zeta}{h_{k}^{2}}\right)\right]\left[\alpha \tau_{y}^{(k)}(\zeta)\right] \frac{h_{k}}{2} d \zeta, \label{eq17.152} \\ F_{2 k}=\int_{-1}^{1}\left[f_{y}^{(k)}(\zeta)\right] \phi_{2}(\zeta) \frac{h_{k}}{2} d \zeta+\int_{-1}^{1}\left[E I_{x x} \frac{(-1+3 \zeta)}{h_{k}}\right]\left[\alpha \tau_{y}^{(k)}(\zeta)\right] \frac{h_{k}}{2} d \zeta, \label{eq17.153} \\ F_{2 k+1}=\int_{-1}^{1}\left[f_{y}^{(k)}(\zeta)\right] \phi_{3}(\zeta) \frac{h_{k}}{2} d \zeta+\int_{-1}^{1}\left[E I_{x x}\left(\frac{6 \zeta}{h_{k}^{2}}\right)\right]\left[\alpha \tau_{y}^{(k)}(\zeta)\right] \frac{h_{k}}{2} d \zeta, \textrm{ and}\label{eq17.154} \\ F_{2 k+2}=\int_{-1}^{1}\left[f_{y}^{(k)}(\zeta)\right] \phi_{4}(\zeta) \frac{h_{k}}{2} d \zeta+\int_{-1}^{1}\left[E I_{x x} \frac{(1+3 \zeta)}{h_{k}}\right]\left[\alpha \tau_{y}^{(k)}(\zeta)\right] \frac{h_{k}}{2} d \zeta. \label{eq17.155}\end{gathered}\]

[sec17.5] Bathe, K-J. Finite Element Procedures in Engineering Analysis. Englewood Cliffs, NJ: Prentice-Hall, Inc., 1982, p. 167.

Guyan, R.J. “Reduction of Stiffness and Mass Matrices.” AIAA Journal 3, no.2, (1965): 380.

Huebner, K. H., E. A.Thornton, and T. G. Byrom. The Finite Element Method for Engineers. 3d ed.New York: John Wiley & Sons, Inc., 1995.[Huebner]

Irons, B,M. “Structural Eigenvalue Problems: Elimination of Unwanted Variables.” AIAA Journal 3, no.5 (1965): 961–962

Reddy, J. N. An Introduction to the Finite Element Method, 4th ed. New York: McGraw Hill-Education, 2019.[Reddy]

Szabo, B., and I. Babuska. Finite Element Analysis. New York: John Wiley & Sons, Inc., 1991.[SzaboBabuska]