Engineering at Alberta Courses » FEA in One Dimension

Finite Element Analysis: FEA in One Dimension

As shown in the previous chapter, the approximate methods involve assuming a particular form for the unknown variables such that the approximate solution has a finite number of unknown parameters. The approximate solution, or often termed trial function, that were used in the previous section were often polynomials which are continuous and differentiable along the whole domain (Figure 1a). The finite element analysis, however, involves using piecewise linear (piecewise affine) or piecewise nonlinear functions for the approximate solution (Figure 1b). Afterwards, a weak formulation of the problem is solved using the Virtual Work method (which as shown previously is equivalent to the Galerkin method).

Figure 1. The final displacement function as (a) the sum of continuously differentiable group
of functions versus (b) the sum of a group of piecewise linear functions

Interpolation Order

In traditional finite element analysis, the values of the displacements at specific points (nodes) across the domain are the main unknowns to be calculated using the method. The displacements at intermediate locations between the nodes are interpolated according to a chosen interpolation function. The interpolation functions are classified according to their differentiability. A $C_{n}$ interpolation function is an interpolation function that can be differentiated $n$ times while an interpolation function that is discontinuous across nodes is termed $C_{-1}$ . Figure 2 shows three nodes (node 1, node 2 and node 3) and three different interpolation functions for the displacement at intermediate locations between the nodes. Figure 2 illustrates discontinuous $C_{-1}$ , continuous $C_{0}$ , and once differentiable $C_{1}$ interpolation functions.

Figure 2. Examples of interpolation order across nodes

One Dimensional Linear Elements

To illustrate the finite element method, we will start by solving the same example that was solved before using the Galerkin method but employing a finite element approximation. Using a four-piecewise linear trial function, find the approximate displacement function of the shown bar. Assume that the bar is linear elastic with Young’s modulus $E$ and cross-sectional area $A$ and that the small strain tensor is the appropriate measure of strain. Ignore the effect of Poisson’s ratio.
PL1
The different aspects of the finite element analysis method will be presented through solving this example as follows.

Shape Functions

In this example, a piecewise $C_0$ linear trial function for the displacement that satisfies the essential boundary condition $@X_1=0:u=0$ is imposed. The displacement is interpolated between the four nodes and thus has the following form:

$u=u_1N_1+u_2N_2+u_3N_3+u_4N_4$

where $\forall i\leq 4:u_i$ are multipliers that happen to be the displacements at the nodes 1, 2, 3, and 4 respectively, while $N_i$ are the interpolation functions that have a value of 1 at node $i$ and decrease linearly to zero at the neighbouring nodes (Figure 3). These functions are traditionally termed the “hat” functions, for obvious reasons. They are also referred to as the shape functions since they define the shape of the deformation of the model. Except when they vanish, the shape functions have the form:

$\begin{split} N_1 & =\bigg\{\begin{matrix} \frac{4X_1}{L} & 0\leq X_1\leq \frac{L}{4}\\ \frac{4}{L}\left(\frac{L}{2}-X_1\right) & \frac{L}{4}\leq X_1\leq \frac{L}{2} \end{matrix}\\ N_2 & =\bigg\{\begin{matrix} \frac{4}{L}\left(X_1-\frac{L}{4}\right) & \frac{L}{4}\leq X_1\leq \frac{L}{2}\\ \frac{4}{L}\left(\frac{3L}{4}-X_1\right) & \frac{L}{2}\leq X_1\leq \frac{3L}{4} \end{matrix}\\ N_3 & =\bigg\{\begin{matrix} \frac{4}{L}\left(X_1-\frac{L}{2}\right) & \frac{L}{2}\leq X_1\leq \frac{3L}{4}\\ \frac{4}{L}\left(L-X_1\right) & \frac{3L}{4}\leq X_1\leq L \end{matrix}\\ N_4 & =\frac{4}{L}\left(X_1-\frac{3L}{4}\right) \end{split}$

The derivatives of the shape functions are very simple and have the form:

$\begin{split} \frac{\mathrm{d}N_1}{\mathrm{d}X_1} & =\bigg\{\begin{matrix} \frac{4}{L} & 0\leq X_1\leq \frac{L}{4}\\ \frac{-4}{L} & \frac{L}{4}\leq X_1\leq \frac{L}{2} \end{matrix}\\ \frac{\mathrm{d}N_2}{\mathrm{d}X_1} & =\bigg\{\begin{matrix} \frac{4}{L}& \frac{L}{4}\leq X_1\leq \frac{L}{2}\\ \frac{-4}{L}& \frac{L}{2}\leq X_1\leq \frac{3L}{4} \end{matrix}\\ \frac{\mathrm{d}N_3}{\mathrm{d}X_1}& =\bigg\{\begin{matrix} \frac{4}{L} & \frac{L}{2}\leq X_1\leq \frac{3L}{4}\\ \frac{-4}{L} & \frac{3L}{4}\leq X_1\leq L \end{matrix}\\ \frac{\mathrm{d}N_4}{\mathrm{d}X_1} & =\frac{4}{L} \end{split}$

Figure 3. One dimensional piecewise linear interpolation functions or “hat” functions

Solution using the Galerkin Method

The Galerkin method, involves equating the integral of the weighted residuals to zero. The weight functions are chosen as the functions $\phi_i(X_1)$ which happen to be equal to the shape functions $N_i$ . After integration by parts as shown previously, the Galerkin method produces four equations of the form:

$\int_0^L \! EA\frac{\mathrm{d}N_i}{\mathrm{d}X_1}\left(\frac{\mathrm{d}u}{\mathrm{d}X_1}\right) \,\mathrm{d}X_1=\int_0^L \! N_i cX_1 \,\mathrm{d}X_1$

There are four trial or shape functions $N_i$ . Each appears only within a sub-interval and is equal to zero at the remaining part of the domain. Thus, the integration for each equation is only performed on the part where $N_i$ does not vanish.

Equation 1: $N_1$ vanishes every where except for $0\leq X_1\leq \frac{L}{2}$ . In this part of the domain, $N_3$ and $N_4$ vanish. Therefore, for the first equation we have:

$\begin{split} u & =u_1N_1+u_2N_2\\ \frac{\mathrm{d}u}{\mathrm{d}X_1} & =u_1\frac{\mathrm{d}N_1}{\mathrm{d}X_1}+u_2\frac{\mathrm{d}N_2}{\mathrm{d}X_1} \end{split}$

Therefore, the first equation has the form:

$u_1\int_0^{\frac{L}{2}} \! EA\frac{\mathrm{d}N_1}{\mathrm{d}X_1}\frac{\mathrm{d}N_1}{\mathrm{d}X_1} \,\mathrm{d}X_1 +u_2\int_0^{\frac{L}{2}} \! EA\frac{\mathrm{d}N_1}{\mathrm{d}X_1}\frac{\mathrm{d}N_2}{\mathrm{d}X_1} \,\mathrm{d}X_1 =\int_0^{\frac{L}{2}} \! N_1cX_1 \,\mathrm{d}X_1$

The domain of integration can be updated to exclude the sub-domains when the integrand is equal to zero:

$u_1\int_0^{\frac{L}{2}} \! EA\frac{\mathrm{d}N_1}{\mathrm{d}X_1}\frac{\mathrm{d}N_1}{\mathrm{d}X_1} \,\mathrm{d}X_1 +u_2\int_{\frac{L}{4}}^{\frac{L}{2}} \! EA\frac{\mathrm{d}N_1}{\mathrm{d}X_1}\frac{\mathrm{d}N_2}{\mathrm{d}X_1} \,\mathrm{d}X_1 =\int_0^{\frac{L}{2}} \! N_1cX_1 \,\mathrm{d}X_1$

Equation 2: $N_2$ vanishes every where except for $\frac{L}{4}\leq X_1\leq \frac{3L}{4}$ . In this part of the domain, $N_4$ vanishes. Therefore, for the second equation we have:

$\begin{split} u & =u_1N_1+u_2N_2+u_3N_3\\ \frac{\mathrm{d}u}{\mathrm{d}X_1} & =u_1\frac{\mathrm{d}N_1}{\mathrm{d}X_1}+u_2\frac{\mathrm{d}N_2}{\mathrm{d}X_1}+u_3\frac{\mathrm{d}N_3}{\mathrm{d}X_1} \end{split}$

Therefore, the second equation has the form:

$u_1\int_{\frac{L}{4}}^{\frac{3L}{4}} \! EA\frac{\mathrm{d}N_2}{\mathrm{d}X_1}\frac{\mathrm{d}N_1}{\mathrm{d}X_1} \,\mathrm{d}X_1 +u_2\int_{\frac{L}{4}}^{\frac{3L}{4}} \!EA \frac{\mathrm{d}N_2}{\mathrm{d}X_1}\frac{\mathrm{d}N_2}{\mathrm{d}X_1} \,\mathrm{d}X_1+u_3\int_{\frac{L}{4}}^{\frac{3L}{4}} \! EA\frac{\mathrm{d}N_2}{\mathrm{d}X_1}\frac{\mathrm{d}N_3}{\mathrm{d}X_1} \,\mathrm{d}X_1 =\int_{\frac{L}{4}}^{\frac{3L}{4}} \! N_2cX_1 \,\mathrm{d}X_1$

The domain of integration can be updated to exclude the sub-domains when the integrand is equal to zero:

$u_1\int_{\frac{L}{4}}^{\frac{L}{2}} \! EA\frac{\mathrm{d}N_2}{\mathrm{d}X_1}\frac{\mathrm{d}N_1}{\mathrm{d}X_1} \,\mathrm{d}X_1 +u_2\int_{\frac{L}{4}}^{\frac{3L}{4}} \! EA\frac{\mathrm{d}N_2}{\mathrm{d}X_1}\frac{\mathrm{d}N_2}{\mathrm{d}X_1} \,\mathrm{d}X_1+u_3\int_{\frac{L}{2}}^{\frac{3L}{4}} \! EA\frac{\mathrm{d}N_2}{\mathrm{d}X_1}\frac{\mathrm{d}N_3}{\mathrm{d}X_1} \,\mathrm{d}X_1 =\int_{\frac{L}{4}}^{\frac{3L}{4}} \! N_2cX_1 \,\mathrm{d}X_1$

Equation 3: $N_3$ vanishes every where except for $\frac{L}{2}\leq X_1\leq L$ . In this part of the domain, $N_1$ vanishes. Therefore, for the third equation we have:

$\begin{split} u & =u_2N_2+u_3N_3+u_4N_4\\ \frac{\mathrm{d}u}{\mathrm{d}X_1} & =u_2\frac{\mathrm{d}N_2}{\mathrm{d}X_1}+u_3\frac{\mathrm{d}N_3}{\mathrm{d}X_1}+u_4\frac{\mathrm{d}N_4}{\mathrm{d}X_1} \end{split}$

Therefore, the third equation after excluding the sub-domains when the intgrand is equal to zero has the form:

$u_2\int_{\frac{L}{2}}^{\frac{3L}{4}} \! EA\frac{\mathrm{d}N_3}{\mathrm{d}X_1}\frac{\mathrm{d}N_2}{\mathrm{d}X_1} \,\mathrm{d}X_1 +u_3\int_{\frac{L}{2}}^{L} \! EA\frac{\mathrm{d}N_3}{\mathrm{d}X_1}\frac{\mathrm{d}N_3}{\mathrm{d}X_1} \,\mathrm{d}X_1+u_4\int_{\frac{3L}{4}}^{L} \! EA\frac{\mathrm{d}N_3}{\mathrm{d}X_1}\frac{\mathrm{d}N_4}{\mathrm{d}X_1} \,\mathrm{d}X_1 =\int_{\frac{L}{2}}^{L} \! N_3cX_1 \,\mathrm{d}X_1$

Equation 4: $N_4$ vanishes every where except for $\frac{3L}{4}\leq X_1\leq L$ . In this part of the domain, $N_1$ and $N_2$ vanish. Therefore, for the fourth equation we have:

$\begin{split} u & =u_3N_3+u_4N_4\\ \frac{\mathrm{d}u}{\mathrm{d}X_1} & =u_3\frac{\mathrm{d}N_3}{\mathrm{d}X_1}+u_4\frac{\mathrm{d}N_4}{\mathrm{d}X_1} \end{split}$

Therefore, the fourth equation after excluding the sub-domains when the intgrand is equal to zero has the form:

$u_3\int_{\frac{3L}{4}}^{L} \! EA\frac{\mathrm{d}N_4}{\mathrm{d}X_1}\frac{\mathrm{d}N_3}{\mathrm{d}X_1} \,\mathrm{d}X_1 +u_4\int_{\frac{3L}{4}}^{L} \! EA\frac{\mathrm{d}N_4}{\mathrm{d}X_1}\frac{\mathrm{d}N_4}{\mathrm{d}X_1} \,\mathrm{d}X_1=\int_{\frac{3L}{4}}^{L} \! N_4cX_1 \,\mathrm{d}X_1$

The above four equations can be written in matrix form as follows:

$\begin{split} \left(\begin{matrix} \int_0^{\frac{L}{2}} \! EA\frac{\mathrm{d}N_1}{\mathrm{d}X_1}\frac{\mathrm{d}N_1}{\mathrm{d}X_1} \,\mathrm{d}X_1 & \int_{\frac{L}{4}}^{\frac{L}{2}} \! EA\frac{\mathrm{d}N_1}{\mathrm{d}X_1}\frac{\mathrm{d}N_2}{\mathrm{d}X_1} \,\mathrm{d}X_1 & 0 & 0\\ \int_{\frac{L}{4}}^{\frac{3L}{4}} \!EA \frac{\mathrm{d}N_2}{\mathrm{d}X_1}\frac{\mathrm{d}N_1}{\mathrm{d}X_1} \,\mathrm{d}X_1 & \int_{\frac{L}{4}}^{\frac{3L}{4}} \! EA\frac{\mathrm{d}N_2}{\mathrm{d}X_1}\frac{\mathrm{d}N_2}{\mathrm{d}X_1} \,\mathrm{d}X_1 & \int_{\frac{L}{4}}^{\frac{3L}{4}} \! EA\frac{\mathrm{d}N_2}{\mathrm{d}X_1}\frac{\mathrm{d}N_3}{\mathrm{d}X_1} \,\mathrm{d}X_1 & 0\\ 0& \int_{\frac{L}{2}}^{\frac{3L}{4}} \! EA\frac{\mathrm{d}N_3}{\mathrm{d}X_1}\frac{\mathrm{d}N_2}{\mathrm{d}X_1} \,\mathrm{d}X_1 & \int_{\frac{L}{2}}^{L} \! EA\frac{\mathrm{d}N_3}{\mathrm{d}X_1}\frac{\mathrm{d}N_3}{\mathrm{d}X_1} \,\mathrm{d}X_1 & \int_{\frac{3L}{4}}^{L} \! EA\frac{\mathrm{d}N_3}{\mathrm{d}X_1}\frac{\mathrm{d}N_4}{\mathrm{d}X_1} \,\mathrm{d}X_1 \\ 0&0&\int_{\frac{3L}{4}}^{L} \! EA\frac{\mathrm{d}N_4}{\mathrm{d}X_1}\frac{\mathrm{d}N_3}{\mathrm{d}X_1} \,\mathrm{d}X_1 &\int_{\frac{3L}{4}}^{L} \! EA\frac{\mathrm{d}N_4}{\mathrm{d}X_1}\frac{\mathrm{d}N_4}{\mathrm{d}X_1} \,\mathrm{d}X_1 \end{matrix} \right) \left(\begin{array}{c}u_1\\u_2\\u_3\\u_4 \end{array}\right)\\ =\left(\begin{array}{c}f_1\\f_2\\f_3\\f_4 \end{array}\right) \end{split}$

where $\forall i\leq 4:f_i=\int_{0}^{L} \! N_icX_1 \,\mathrm{d}X_1$ . The integration on both sides is relatively simple to calculate. The final equations have the following form:

$\left(\begin{matrix} \frac{8}{L} & \frac{-4}{L}& 0 & 0\\ \frac{-4}{L} & \frac{8}{L} & \frac{-4}{L} & 0\\ 0& \frac{-4}{L} & \frac{8}{L} & \frac{-4}{L} \\ 0&0&\frac{-4}{L} & \frac{4}{L} \end{matrix} \right) \left(\begin{array}{c}u_1\\u_2\\u_3\\u_4 \end{array}\right)=\frac{cL^2}{EA}\left(\begin{array}{c}\frac{1}{16}\\\frac{1}{8}\\\frac{3}{16}\\\frac{11}{96}\end{array}\right)$

which when solved gives the solution:

$\left(\begin{array}{c}u_1\\u_2\\u_3\\u_4 \end{array}\right)=\frac{cL^3}{EA}\left(\begin{array}{c}\frac{47}{384}\\\frac{11}{48}\\\frac{39}{128}\\\frac{1}{3}\end{array}\right)$

Comparison with the Exact Solution

The exact solution as shown previously is known and is equal to:

$u_{\mbox{exact}}=\frac{cL^2}{2EA}X_1-\frac{c}{6EA}X_1^3$

The exact displacement at the points 1, 2, 3, and 4 are given by:

$\begin{split} {u_{\mbox{exact}}}_1&=\frac{47cL^3}{384EA}\\ {u_{\mbox{exact}}}_2&=\frac{11cL^3}{48EA}\\ {u_{\mbox{exact}}}_3&=\frac{39cL^3}{128EA}\\ {u_{\mbox{exact}}}_4&=\frac{cL^3}{3EA} \end{split}$

The finite element analysis approximation gave exact solution for the displacements of the points 1, 2, 3, and 4. However, in between those points, the finite element analysis approximation following linear interpolation. As shown in the figure below when setting $c=1$ , $E=1$ , $A=1$ , and $L=1$ units, the displacement values are almost accurate. However, the linear interpolation of the displacements means that the strain is constant in the area between points. Therefore, the stresses are constant between points and in fact, the stress distribution is discontinuous!
FEA4

View Mathematica Code

uexact = c*L^2/2/EA*x - c/6/EA*x^3;
N1 = Piecewise[{{4 x/L, 0 <= x <= L/4}, {4/L (L/2 - x), L/4 <= x <= L/2}}, 0];
N2 = Piecewise[{{4/L (x - L/4), L/4 <= x <= L/2}, {4/L (3 L/4 - x), L/2 <= x <= 3 L/4}}, 0];
N3 = Piecewise[{{4/L (x - L/2), L/2 <= x <= 3 L/4}, {4/L (L - x), 3 L/4 <= x <= L}}, 0];
N4 = Piecewise[{{4/L (x - 3 L/4), 3 L/4 <= x <= L}}, 0];
u1 = c*L^3/EA*(47/384);
u2 = c*L^3/EA*(11/48);
u3 = c*L^3/EA*(39/128);
u4 = c*L^3/EA*(1/3);
u = u1*N1 + u2*N2 + u3*N3 + u4*N4;
uexactp = uexact /. {c -> 1, EA -> 1, L -> 1}
sigmaexact = D[uexactp, x]
up = u /. {c -> 1, EA -> 1, L -> 1}
(*Assuming Young's modulus = 1 unit *)
sigmap = D[up, x]
Plot[{uexactp, up}, {x, 0, 1}, AxesLabel -> {"x", "displacement"}, PlotLegends -> {uexact, u_FEA}]
Plot[{sigmaexact, sigmap}, {x, 0, 1}, AxesLabel -> {"x", "sigma_11"},  PlotLegends -> {"sigma_exact", "sigma_FEA"}]

Solution Using the Virtual Work Method

As stated previously, the Galerkin method and the virtual work method are equivalent. The principle of virtual work states that, from a state of equilibrium, when a virtual admissible (compatible) displacement is applied, the increment in the work done by the external and body forces during the application of the virtual displacement is equal to the increment in the deformation energy stored. We will adopt the traditional matrix multiplication convention with differentiation between row and column vectors. The displacement function has the form:

$u=u_1N_1+u_2N_2+u_3N_3+u_4N_4=\left<\begin{array}{cccc}N_1 & N_2 & N_3 & N_4\end{array}\right>\left(\begin{array}{c}u_1\\u_2\\u_3\\u_4\end{array}\right)$

The stress component $\sigma_{11}$ due to the assumed trial displacement function has the following form:

$\begin{split} \sigma_{11}&=E\varepsilon_{11}=E\frac{\mathrm{d}u}{\mathrm{d}X_1}=E\sum_{i=1}^4u_i\frac{\mathrm{d}N_i}{\mathrm{d}X_1}\\ &=E\left<\begin{array}{cccc}\frac{\mathrm{d}N_1}{\mathrm{d}X1} & \frac{\mathrm{d}N_2}{\mathrm{d}X1} & \frac{\mathrm{d}N_3}{\mathrm{d}X1} & \frac{\mathrm{d}N_4}{\mathrm{d}X1}\end{array}\right>\left(\begin{array}{c}u_1\\u_2\\u_3\\u_4\end{array}\right) \end{split}$

The virtual displacement can be set by choosing virtual parameters $u_i^*$ :

$u^*=u_1^*N_1+u_2^*N_2+u_3^*N_3+u_4^*N_4$

Setting:

$\begin{split} N&=\left<\begin{array}{cccc}N_1 & N_2 & N_3 & N_4\end{array}\right>\\ B&=\left<\begin{array}{cccc}\frac{\mathrm{d}N_1}{\mathrm{d}X1} & \frac{\mathrm{d}N_2}{\mathrm{d}X1} & \frac{\mathrm{d}N_3}{\mathrm{d}X1} & \frac{\mathrm{d}N_4}{\mathrm{d}X1}\end{array}\right>\\ u_e&=\left(\begin{array}{c}u_1\\u_2\\u_3\\u_4\end{array}\right)\\ u_e^*&=\left<\begin{array}{cccc}u_1^* & u_2^*& u_3^* & u_4^*\end{array}\right>\\ \end{split}$

then the equations can be rewritten as follows:

$\begin{split} u&=Nu_e\\ u^*&=u_e^*N^T\\ \varepsilon_{11}&=Bu_e\\ \varepsilon_{11}^*&=u_e^*B^T\\ \sigma_{11}&=EBu_e \end{split}$

The internal virtual work has the form:

$IVW=\int_0^{L} \! A\varepsilon_{11}^*\sigma_{11} \,\mathrm{d}X_1=u_e^* \int_0^{L} \! EAB^TB \,\mathrm{d}X_1 u_e$

The external virtual work has the form:

$EVW= u_e^*\int_0^{L} \! cX_1 N^T \,\mathrm{d}X_1$

Equating the internal and external virtual work, the following is obtained:

$\begin{split} &{\small\left<\begin{array}{cccc}u_1^* & u_2^*& u_3^* & u_4^*\end{array}\right> \left(\begin{matrix} \int_0^{\frac{L}{2}} \! \frac{\mathrm{d}N_1}{\mathrm{d}X_1}\frac{\mathrm{d}N_1}{\mathrm{d}X_1} \,\mathrm{d}X_1 & \int_{\frac{L}{4}}^{\frac{L}{2}} \! \frac{\mathrm{d}N_1}{\mathrm{d}X_1}\frac{\mathrm{d}N_2}{\mathrm{d}X_1} \,\mathrm{d}X_1 & 0 & 0\\ \int_{\frac{L}{4}}^{\frac{3L}{4}} \! \frac{\mathrm{d}N_2}{\mathrm{d}X_1}\frac{\mathrm{d}N_1}{\mathrm{d}X_1} \,\mathrm{d}X_1 & \int_{\frac{L}{4}}^{\frac{3L}{4}} \! \frac{\mathrm{d}N_2}{\mathrm{d}X_1}\frac{\mathrm{d}N_2}{\mathrm{d}X_1} \,\mathrm{d}X_1 & \int_{\frac{L}{4}}^{\frac{3L}{4}} \! \frac{\mathrm{d}N_2}{\mathrm{d}X_1}\frac{\mathrm{d}N_3}{\mathrm{d}X_1} \,\mathrm{d}X_1 & 0\\ 0& \int_{\frac{L}{2}}^{\frac{3L}{4}} \! \frac{\mathrm{d}N_3}{\mathrm{d}X_1}\frac{\mathrm{d}N_2}{\mathrm{d}X_1} \,\mathrm{d}X_1 & \int_{\frac{L}{2}}^{L} \! \frac{\mathrm{d}N_3}{\mathrm{d}X_1}\frac{\mathrm{d}N_3}{\mathrm{d}X_1} \,\mathrm{d}X_1 & \int_{\frac{3L}{4}}^{L} \! \frac{\mathrm{d}N_3}{\mathrm{d}X_1}\frac{\mathrm{d}N_4}{\mathrm{d}X_1} \,\mathrm{d}X_1 \\ 0&0&\int_{\frac{3L}{4}}^{L} \! \frac{\mathrm{d}N_4}{\mathrm{d}X_1}\frac{\mathrm{d}N_3}{\mathrm{d}X_1} \,\mathrm{d}X_1 &\int_{\frac{3L}{4}}^{L} \! \frac{\mathrm{d}N_4}{\mathrm{d}X_1}\frac{\mathrm{d}N_4}{\mathrm{d}X_1} \,\mathrm{d}X_1 \end{matrix} \right) \left(\begin{array}{c}u_1\\u_2\\u_3\\u_4 \end{array}\right)}\\ &=\frac{1}{EA}\left<\begin{array}{cccc}u_1^* & u_2^*& u_3^* & u_4^*\end{array}\right>\left(\begin{array}{c}f_1\\f_2\\f_3\\f_4 \end{array}\right) \end{split}$

where $\forall i\leq 4:f_i=\int_{0}^{L} \! N_icX_1 \,\mathrm{d}X_1$ . Note that the factor $EA$ is constant in this problem and was pulled out as a common factor for simplicity.
Since the coefficients $u_i^*$ are arbitrary, the coefficients of each $u_i^*$ on both sides of the equality are equal. Therefore, the same equations of the Galerkin method are retrieved.

Elements, Nodes, and Degrees of Freedom

In the finite element method, the domain is discretized into points which are called nodes and the regions between the nodes are called elements. The approximate solution involves assuming a trial function for the displacements interpolated between the values of the displacements the nodes. The displacement function is formed by multiplying the displacements of the nodes by shape functions. The displacements of the nodes are termed the degrees of freedom of the system. The right-hand side of the matrix equation formed is equivalent to a nodal force vector and is formed by lumping the distributed load and forming equivalent concentrated loads at the nodes.

For this particular example the solution is exact at the selected nodes. However, the solution within the elements will be linearly interpolated between the two values on each side and thus is not exact. The value of the gradient of the displacement at each node does not exist, as the displacement function is not differentiable at the nodes. In the finite element analysis method, the stresses that are related to the strains (the gradients of the displacements) are usually averaged at the nodes, and the degree of variation of the stresses between the different elements could be an indication of the accuracy of the solution. The matrix multiplied by the degrees of freedom is, in conservative systems, symmetric and is termed the stiffness matrix.

When the virtual work method was applied, the integration was done on the whole domain. However, it is possible to perform the integration locally on the elements. The equations of the virtual work can be written as follows:

$\begin{split} &\int_0^{\frac{L}{4}} \! A\varepsilon_{11}^*\sigma_{11} \,\mathrm{d}X_1+\int_{\frac{L}{4}}^{\frac{L}{2}} \! A\varepsilon_{11}^*\sigma_{11} \,\mathrm{d}X_1+\int_{\frac{L}{2}}^{\frac{3L}{4}} \! A\varepsilon_{11}^*\sigma_{11} \,\mathrm{d}X_1+\int_{\frac{3L}{4}}^{L} \! A\varepsilon_{11}^*\sigma_{11} \,\mathrm{d}X_1\\ &= \int_0^{\frac{L}{4}} \! u^*cX_1 \,\mathrm{d}X_1+\int_{\frac{L}{4}}^{\frac{L}{2}} \! u^*cX_1 \,\mathrm{d}X_1+\int_{\frac{L}{2}}^{\frac{3L}{4}} \! u^*cX_1 \,\mathrm{d}X_1+\int_{\frac{3L}{4}}^{L} \! u^*cX_1 \,\mathrm{d}X_1 \end{split}$

Because the integration can be done locally on the elements, it is possible to isolate the elements and perform the integration on each element separately (Figure 4). On each element, the shape functions or interpolation functions are associated with the end nodes. For example, element e2 connects nodes 1 and 2 with two shape functions,
$N_1$ and $N_2$ . What is remarkable about the method is the similarity between all the elements shown in (Figure 4). Thus, a local stiffness matrix for each element can be developed, and then, the global stiffness matrix can be easily assembled by combining all the local stiffness matrices.

Element 1:
On element 1, the displacement is given by: $u=u_1N_1$ and the virtual displacement is given by $u^*=u_1^*N_1$ . By considering the first element, the first term (denoted $IVW_1$ ) in the internal virtual work equation has the following form:

$IVW_1=u_1^* \int_0^{\frac{L}{4}} \! EA\frac{\mathrm{d}N_1}{\mathrm{d}X_1}\frac{\mathrm{d}N_1}{\mathrm{d}X_1} \,\mathrm{d}X_1 u_1=u_1^*K^{e1}_{11} u_1$

$K^{e1}_{11}$ is the local stiffness entry for element 1 corresponding to the degree of freedom $u_1$ .

The first term in the external virtual work has the form:

$EVW_1=u_1^*\int_0^{\frac{L}{4}} \! N_1cX_1 \,\mathrm{d}X_1=u_1^* f^{e1}_1$

where $f^{e1}_1$ is the nodal force to be applied at node 1 corresponding to lumping the distributed load acting on element 1.

Element 2:
On element 2, the displacement is given by: $u=u_1N_1+u_2N_2$ and the virtual displacement is given by $u^*=u_1^*N_1+u_2^*N_2$ . By considering the second element, the second term in the internal virtual work equation has the following form:

$\begin{split} IVW_2&= \left<\begin{array}{cc}u_1^* & u_2^*\end{array}\right> \left(\begin{matrix} \int_{\frac{L}{4}}^{\frac{L}{2}} \! EA\frac{\mathrm{d}N_1}{\mathrm{d}X_1}\frac{\mathrm{d}N_1}{\mathrm{d}X_1} \,\mathrm{d}X_1 & \int_{\frac{L}{4}}^{\frac{L}{2}} \! EA\frac{\mathrm{d}N_1}{\mathrm{d}X_1}\frac{\mathrm{d}N_2}{\mathrm{d}X_1} \,\mathrm{d}X_1\\ \int_{\frac{L}{4}}^{\frac{L}{2}} \! EA\frac{\mathrm{d}N_2}{\mathrm{d}X_1}\frac{\mathrm{d}N_1}{\mathrm{d}X_1} \,\mathrm{d}X_1 & \int_{\frac{L}{4}}^{\frac{L}{2}} \! EA\frac{\mathrm{d}N_2}{\mathrm{d}X_1}\frac{\mathrm{d}N_2}{\mathrm{d}X_1} \,\mathrm{d}X_1 \end{matrix} \right) \left(\begin{array}{c}u_1\\u_2\end{array}\right)\\ &=\left<\begin{array}{cc}u_1^* & u_2^*\end{array}\right> \left(\begin{matrix} K^{e2}_{11} & K^{e2}_{12}\\K^{e2}_{21} & K^{e2}_{22} \end{matrix}\right)\left(\begin{array}{c}u_1\\u_2\end{array}\right) \end{split}$

where $K^{e2}_{ij}$ are the local stiffness matrix entries for element 2 corresponding to the degrees of freedom $u_1$ and $u_2$ .
The second term in the external virtual work has the form:

$\begin{split} EVW_2&=\left<\begin{array}{cc}u_1^* & u_2^*\end{array}\right>\left(\begin{array}{c}\int_\frac{L}{4}^{\frac{L}{2}} \! N_1cX_1 \,\mathrm{d}X_1\\\int_\frac{L}{4}^{\frac{L}{2}} \! N_2cX_1 \,\mathrm{d}X_1\end{array}\right)\\ &=\left<\begin{array}{cc}u_1^* & u_2^*\end{array}\right>\left(\begin{array}{c}f^{e2}_1\\f^{e2}_2\end{array}\right) \end{split}$

where $f^{e2}_1$ and $f^{e2}_2$ are the nodal forces to be applied at nodes 1 and 2 corresponding to lumping the distributed load acting on element 2.

Element 3:
On element 3, the displacement is given by: $u=u_2N_2+u_3N_3$ and the virtual displacement is given by $u^*=u_2^*N_2+u_3^*N_3$ . By considering the third element, the third term in the internal virtual work equation has the following form:

$\begin{split} IVW_3&=\left<\begin{array}{cc}u_2^* & u_3^*\end{array}\right> \left(\begin{matrix} \int_{\frac{L}{2}}^{\frac{3L}{4}} \! EA\frac{\mathrm{d}N_2}{\mathrm{d}X_1}\frac{\mathrm{d}N_2}{\mathrm{d}X_1} \,\mathrm{d}X_1 & \int_{\frac{L}{2}}^{\frac{3L}{4}} \! EA\frac{\mathrm{d}N_2}{\mathrm{d}X_1}\frac{\mathrm{d}N_3}{\mathrm{d}X_1} \,\mathrm{d}X_1\\ \int_{\frac{L}{2}}^{\frac{3L}{4}} \! EA\frac{\mathrm{d}N_3}{\mathrm{d}X_1}\frac{\mathrm{d}N_2}{\mathrm{d}X_1} \,\mathrm{d}X_1 & \int_{\frac{L}{2}}^{\frac{3L}{4}} \! EA\frac{\mathrm{d}N_3}{\mathrm{d}X_1}\frac{\mathrm{d}N_3}{\mathrm{d}X_1} \,\mathrm{d}X_1 \end{matrix}\right) \left(\begin{array}{c}u_2\\u_3\end{array}\right)\\ &=\left<\begin{array}{cc}u_2^* & u_3^*\end{array}\right> \left(\begin{matrix} K^{e3}_{11} & K^{e3}_{12}\\K^{e3}_{21} & K^{e3}_{22} \end{matrix}\right)\left(\begin{array}{c}u_2\\u_3\end{array}\right) \end{split}$

where $K^{e3}_{ij}$ are the local stiffness matrix entries for element 3 corresponding to the degrees of freedom $u_2$ and $u_3$ .
The third term in the external virtual work has the form:

$\begin{split} EVW_3&=\left<\begin{array}{cc}u_2^* & u_3^*\end{array}\right>\left(\begin{array}{c}\int_\frac{L}{2}^{\frac{3L}{4}} \! N_2cX_1 \,\mathrm{d}X_1\\\int_\frac{L}{2}^{\frac{3L}{4}} \! N_3cX_1 \,\mathrm{d}X_1\end{array}\right)\\ &=\left<\begin{array}{cc}u_2^* & u_3^*\end{array}\right>\left(\begin{array}{c}f^{e3}_1\\f^{e3}_2\end{array}\right) \end{split}$

where $f^{e3}_1$ and $f^{e3}_2$ are the nodal forces to be applied at nodes 2 and 3 corresponding to lumping the distributed load acting on element 3.

Element 4:
On element 4, the displacement is given by: $u=u_3N_3+u_4N_4$ and the virtual displacement is given by $u^*=u_3^*N_3+u_4^*N_4$ . By considering the fourth element, the fourth term in the internal virtual work equation has the following form:

$\begin{split} IVW_4&=\left<\begin{array}{cc}u_3^* & u_4^*\end{array}\right> \left(\begin{matrix} \int_{\frac{3L}{4}}^{L} \! EA\frac{\mathrm{d}N_3}{\mathrm{d}X_1}\frac{\mathrm{d}N_3}{\mathrm{d}X_1} \,\mathrm{d}X_1 & \int_{\frac{3L}{4}}^{L} \! EA\frac{\mathrm{d}N_3}{\mathrm{d}X_1}\frac{\mathrm{d}N_4}{\mathrm{d}X_1} \,\mathrm{d}X_1\\ \int_{\frac{3L}{4}}^{L} \! EA\frac{\mathrm{d}N_4}{\mathrm{d}X_1}\frac{\mathrm{d}N_3}{\mathrm{d}X_1} \,\mathrm{d}X_1 & \int_{\frac{3L}{4}}^{L} \! EA\frac{\mathrm{d}N_4}{\mathrm{d}X_1}\frac{\mathrm{d}N_4}{\mathrm{d}X_1} \,\mathrm{d}X_1 \end{matrix}\right)\left(\begin{array}{c}u_3\\u_4\end{array}\right)\\ &=\left<\begin{array}{cc}u_3^* & u_4^*\end{array}\right> \left(\begin{matrix} K^{e4}_{11} & K^{e4}_{12}\\K^{e4}_{21} & K^{e4}_{22} \end{matrix}\right)\left(\begin{array}{c}u_3\\u_4\end{array}\right) \end{split}$

where $K^{e4}_{ij}$ are the local stiffness matrix entries for element 4 corresponding to the degrees of freedom $u_3$ and $u_4$ .
The fourth term in the external virtual work has the form:

$\begin{split} EVW_4&=\left<\begin{array}{cc}u_3^* & u_4^*\end{array}\right>\left(\begin{array}{c}\int_\frac{3L}{4}^{L} \! N_3cX_1 \,\mathrm{d}X_1\\\int_\frac{3L}{4}^{L} \! N_4cX_1 \,\mathrm{d}X_1 \end{array}\right)\\ &=\left<\begin{array}{cc}u_3^* & u_4^*\end{array}\right>\left(\begin{array}{c}f^{e4}_1\\f^{e4}_2\end{array}\right) \end{split}$

where $f^{e4}_1$ and $f^{e4}_2$ are the nodal forces to be applied at nodes 3 and 4 corresponding to lumping the distributed load acting on element 4.

Global Stiffness Matrix
The global stiffness matrix can be assembled by the equation:

$IVW_1+IVW_2+IVW_3+IVW_4=EVW_1+EVW_2+EVW_3+EVW_4$

The resulting equations (global stiffness matrix and global nodal forces vector) considering the arbitrariness of the multipliers $u_i^*$ :

$\left(\begin{matrix} K^{e1}_{11}+K^{e2}_{11} & K^{e2}_{12} & 0 & 0\\ K^{e2}_{21} & K^{e2}_{22}+K^{e3}_{11} & K^{e3}_{12} & 0\\ 0& K^{e3}_{21}& K^{e3}_{22}+K^{e4}_{11} & K^{e4}_{12}\\ 0& 0& K^{e4}_{21} & K^{e4}_{22} \end{matrix}\right)\left(\begin{array}{c}u_1\\u_2\\u_3\\u_4\end{array}\right)=\left(\begin{array}{c}f^{e1}_1+f^{e2}_1\\f^{e2}_2+f^{e3}_1\\f^{e3}_2+f^{e4}_1\\f^{e4}_2\end{array}\right)$

It is important to note that because the assumed displacement function satisfies the essential boundary conditions, the structure is stable, and the equations can be solved directly. In general, the structure under consideration would have five degrees of freedom and would have a $5\times 5$ global stiffness matrix and a $5\times 1$ global nodal forces vector. The nodal forces vector corresponding to the displacement of the first node where the displacement is prescribed would have an unknown reaction. The stiffness matrix can then be reduced into a $4\times 4$ matrix which can be solved.

Figure 4. Discretizing the domain into elements and nodes

Properties of the Stiffness Matrix

In the example above, the final equations had the following form:

$Ku=F$

where $K\in\mathbb{M}^n$ is the $n\times n$ global stiffness matrix, $u\in\mathbb{R}^n$ is the vector of degrees of freedom while $F\in\mathbb{R}^n$ is the nodal forces vector. In general, the global stiffness matrix $K$ of an elastic structure formed using the finite element analysis method whether the problems has one, two, or three dimensions has the following properties:

Sparsity: The stiffness matrix formed in the traditional finite element analysis is sparse, i.e., contains many zero entries. This property is due to the choice of piecewise-linear interpolation functions. The shape functions (or the interpolation functions) are chosen to be “Hat functions” that take values of 1 at each node and decrease linearly or in a parabolic fashion across the neighboring elements connected to that node. Nodes that are not connected together with any elements have a corresponding zero entry in the global stiffness matrix.
Invertibility: (See the section on Invertible linear maps) The stiffness matrix of a structure is a linear map between $u\in\mathbb{R}^n$ and $F \in\mathbb{R}^n$ . Such map is invertible if and only if the map is injective (one to one), which means that the only element in the kernel of the map is the zero vector. Otherwise, if the map is not injective, then there are two nonzero distinct displacement vectors $v$ and $w\in\mathbb{R}^n$ that can be produced by the same force vector $F \in\mathbb{R}^n$ . I.e., $K(v-w)=F-F=0$ , then the nonzero vector displacement $v-w$ can be applied to the structure without any external forces! If no boundary conditions are applied to a structure, then naturally rigid body motions and rigid body rotations can be applied without any external forces. Therefore, a stiffness matrix of a structure that is not restrained is not invertible and its determinant is equal to zero. Once enough boundary conditions are applied to the structure, then the only element of the kernel of the linear map is the zero element, the stiffness matrix is invertible, and its determinant is not equal to zero. Notice, however, that when a structure is not well supported and does not have the correct boundary conditions, the stiffness matrix is termed “ill conditioned”. Some finite element analysis software will still be able to find a possible solution, depending on the solution method. A user should be cautious in such cases since the solution could involve very large numbers for the displacement that would render the solution inaccurate.
Symmetry: The stiffness matrix of an elastic structure is symmetric. The symmetry is a direct consequence of elasticity. Elasticity implies that there is an energy function such that $F_i=\frac{\partial W}{\partial u_i}$ . Therefore, $K_{ij}=\frac{\partial F_i}{\partial u_j}=\frac{\partial^2 W}{\partial u_i\partial u_j}$ implying that $K$ is symmetric.
Positive and Semi-Positive Definiteness: (See the corresponding section in linear vector maps) A stable structure will deform in the direction of increasing external forces. Mathematically speaking, the dot product between the displacement and the force vectors has to be positive. If the structure is well restrained, i.e., $K$ is invertible, then $K$ is positive definite and any displacement field in the structure will produce positive energy. In mathematical terms, this imples: $\forall u \in\mathbb{R}^n:F.u=Ku.u > 0$ . However, if the structure is not well restrained, then $K$ is not invertible and $\exists u\in\mathbb{R}^n$ such that $Ku = O$ . If the material of the structure is elastic, then the global stiffness matrix of the structure is semi-positive definite. Most finite element analysis software will check the invertibility of the stiffness matrix by finding the smallest eigenvalue of the stiffness matrix. If one of the eigenvalues of the stiffness matrix is zero, then the stiffness matrix is not invertible. If one of the eigenvalues is negative, then the stiffness matrix is not positive definite, since positive definiteness ensures positive eigenvalues (See the corresponding section in linear vector maps). A direct consequence of the positive definiteness is that the diagonal entries of the stiffness matrix have to be positive; since a diagonal element represents the force applied at the corresponding degree of freedom to produce unit displacement at this particular degree of freedom, then the requirement that the energy is positive ensures that the force and the displacement are in the same direction.
Any vertical column of the stiffness matrix is equal to the vector of nodal forces required to move the corresponding degree of freedom by one unit displacement while keeping the remaining degrees of freedom equal to zero.

One Dimensional Quadratic Elements

In the previous section, the domain of a one-dimensional problem was discretized using piecewise linear functions. Another approximation is the C0 piecewise quadratic interpolation functions, where on each element, the displacement is quadratic. If one wishes to obtain better accuracy, either a finer mesh of linear elements is to be used, or quadratic elements are to be used instead of the linear elements (Figure 5).

Figure 5. Linear vs. quadratic one dimensional elements

To define a quadratic one dimensional element, an additional inner node is introduced. Similar to the linear elements, the displacement function is continuous across nodes but not differentiable at the boundaries between elements. In general, the quadratic form of the displacement is:

$u=a_0+a_1X_1+a_2X_1^2$

Figure 6 shows the displacement function of a quadratic element with three nodes. An alternative form of the displacement is:

$u=u_1N_1+u_2N_2+u_3N_3$

where $u_1$ , $u_2$ , and $u_3$ are the displacements of the nodes 1, 2, and 3, while $N_1$ , and $N_2$ , and $N_3$ are the shape functions associated with these nodes (Figure 6). $a_0$ , $a_1$ , and $a_2$ are the generalized degrees of freedom of the element while $u_1$ , $u_2$ , and $u_3$ are the nodal degrees of freedom. To find the shape functions $N_1$ , $N_2$ , and $N_3$ , the following three equations are solved:

$\begin{split} @X_1=0&:u=u_1\Rightarrow a_0=u_1\\ @X_1=\frac{L}{2}&:u=u_2\Rightarrow a_0+a_1\frac{L}{2}+a_2\frac{L^2}{4}=u_2\\ @X_1=L&:u=u_3\Rightarrow a_0+a_1L+a_2L^2=u_3 \end{split}$

Therefore,

$\begin{split} a_0&=u_1\\ a_1&=-\frac{3u_1-4u_2+u_3}{L}\\ a_2&=\frac{2(u_1-2u_2+u_3)}{L^2} \end{split}$

Therefore, the displacement function in terms of the nodal degrees of freedom has the form:

$u=u_1-\frac{3u_1-4u_2+u_3}{L}X_1+\frac{2(u_1-2u_2+u_3)}{L^2}X_1^2$

Rearranging:

$u=u_1\left(\frac{(L-2X_1)(L-X_1)}{L^2}\right)+u_2\left(\frac{4(L-X_1)X_1}{L^2}\right)+u_3\left(\frac{-(L-2X_1)X_1}{L^2}\right)$

Therefore, the shape functions are:

$\begin{split} N_1&=\frac{(L-2X_1)(L-X_1)}{L^2}\\ N_2&=\frac{4(L-X_1)X_1}{L^2}\\ N_3&=\frac{-(L-2X_1)X_1}{L^2} \end{split}$

Figure 6. One dimensional quadratic element displacement function and shape functions

It can be checked that the sum of the shape functions is always equal to 1, which allows for the rigid-body motion or constant displacement of an element to be modelled (Why?). The local stiffness matrix and the local nodal forces vectors can be established similar to the previous section. First, the following matrices are defined:

$\begin{split} N&=\left<\begin{array}{cccc}N_1 & N_2 & N_3\end{array}\right>\\ B&=\left<\begin{array}{cccc}\frac{\mathrm{d}N_1}{\mathrm{d}X1} & \frac{\mathrm{d}N_2}{\mathrm{d}X1} & \frac{\mathrm{d}N_3}{\mathrm{d}X1} \end{array}\right>\\ u_e&=\left(\begin{array}{c}u_1\\u_2\\u_3\end{array}\right)\\ \end{split}$

then the displacement $u$ , $\varepsilon_{11}$ , and $\sigma_{11}$ can be written in terms of the nodal degrees of freedom vector $u_e$ as follows:

$\begin{split} u&=Nu_e\\ \varepsilon_{11}&=Bu_e\\ \sigma_{11}&=EBu_e \end{split}$

If $E$ is Young’s modulus and $A$ is the cross sectional area of the element, then the local stiffness matrix of the element is a $3\times 3$ matrix and has the form:

$K^e=\int_0^{L} \! EA B^TB \, \mathrm{d}X_1$

If $E$ and $A$ are constant, then $K^e$ has the form:

$K^e=\frac{EA}{3L}\left(\begin{matrix}7 & -8 & 1\\-8 &16 & -8\\1 & -8 & 7\end{matrix}\right)$

The nodal forces vector is a 3 dimensional vector and has the form:

$f^e=\int_0^{L} \! p N^T \, \mathrm{d}X_1$

where $p$ is the distributed load per unit length acting on the element. If $p$ is constant, then the nodal forces vector have the form:

$f^e=pL\left(\begin{array}{c}\frac{1}{6}\\\frac{2}{3}\\\frac{1}{6}\end{array}\right)$