# Finite elements/Axial bar finite element solution

## Axially loaded bar: The Finite Element Solution

The finite element method is a type of Galerkin method that has the following advantages:

1. The functions ${\displaystyle \varphi _{i}}$ are found in a systematic manner.
2. The functions ${\displaystyle \varphi _{i}}$ are chosen such that they can be used for arbitrary domains.
3. The functions ${\displaystyle \varphi _{i}}$ are piecewise polynomials.
4. The functions ${\displaystyle \varphi _{i}}$ are non-zero only on a small part of the domain.

As a result, computations can be done in a modular manner that is suitable for computer implementation.

### Discretization

The first step in the finite element approach is to divide the domain into elements and nodes, i.e., to create the finite element mesh.

Let us consider a simple situation and divide the rod into 3 elements and 4 nodes as shown in Figure 6.

 Figure 6. Finite element mesh and basis functions for the bar.

### Shape functions

The functions ${\displaystyle \varphi _{i}}$ have special characteristics in finite element methods and are generally written as ${\displaystyle N_{i}}$ and are called basis functions, shape functions, or interpolation functions.

Therefore, we may write

{\displaystyle {\begin{aligned}K_{ij}&=\int _{0}^{L}{\cfrac {dN_{j}}{dx}}AE{\cfrac {dN_{i}}{dx}}~dx\\f_{j}&=\int _{0}^{L}N_{j}\mathbf {q} ~dx+\left.N_{j}{\boldsymbol {R}}\right|_{x=L}~.\end{aligned}}}

The finite element basis functions are chosen such that they have the following properties:

1. The functions ${\displaystyle N_{i}}$ are bounded and continuous.
2. If there are ${\displaystyle n}$ nodes, then there are ${\displaystyle n}$ basis functions - one for each node. There are four basis functions for the mesh shown in Figure 6.
3. Each function ${\displaystyle N_{i}}$ is nonzero only on elements connected to node ${\displaystyle i}$.
4. ${\displaystyle N_{i}}$ is 1 at node ${\displaystyle i}$ and zero at all other nodes.

### Stiffness matrix

Let us compute the values of ${\displaystyle K_{ij}}$ for the three element mesh. We have

{\displaystyle {\begin{aligned}K_{ij}&=\int _{0}^{L}{\cfrac {dN_{j}}{dx}}AE{\cfrac {dN_{i}}{dx}}~dx\\&=\int _{0}^{x_{2}}{\cfrac {dN_{j}}{dx}}AE{\cfrac {dN_{i}}{dx}}~dx+\int _{x_{2}}^{x_{3}}{\cfrac {dN_{j}}{dx}}AE{\cfrac {dN_{i}}{dx}}~dx+\int _{x_{3}}^{L}{\cfrac {dN_{j}}{dx}}AE{\cfrac {dN_{i}}{dx}}~dx\\&\equiv \int _{\Omega _{1}}{\cfrac {dN_{j}}{dx}}AE{\cfrac {dN_{i}}{dx}}~dx+\int _{\Omega _{2}}{\cfrac {dN_{j}}{dx}}AE{\cfrac {dN_{i}}{dx}}~dx+\int _{\Omega _{3}}{\cfrac {dN_{j}}{dx}}AE{\cfrac {dN_{i}}{dx}}~dx\end{aligned}}}

The components of ${\displaystyle \mathbf {K} }$ are

{\displaystyle {\begin{aligned}K_{11}&=\int _{\Omega _{1}}{\cfrac {dN_{1}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx+\int _{\Omega _{2}}{\cfrac {dN_{1}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx+\int _{\Omega _{3}}{\cfrac {dN_{1}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx\\K_{12}&=\int _{\Omega _{1}}{\cfrac {dN_{2}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx+\int _{\Omega _{2}}{\cfrac {dN_{2}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx+\int _{\Omega _{3}}{\cfrac {dN_{2}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx\\K_{13}&=\int _{\Omega _{1}}{\cfrac {dN_{3}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx+\int _{\Omega _{2}}{\cfrac {dN_{3}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx+\int _{\Omega _{3}}{\cfrac {dN_{3}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx\\K_{14}&=\int _{\Omega _{1}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx+\int _{\Omega _{2}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx+\int _{\Omega _{3}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx\\K_{22}&=\int _{\Omega _{1}}{\cfrac {dN_{2}}{dx}}AE{\cfrac {dN_{2}}{dx}}~dx+\int _{\Omega _{2}}{\cfrac {dN_{2}}{dx}}AE{\cfrac {dN_{2}}{dx}}~dx+\int _{\Omega _{3}}{\cfrac {dN_{2}}{dx}}AE{\cfrac {dN_{2}}{dx}}~dx\\K_{23}&=\int _{\Omega _{1}}{\cfrac {dN_{3}}{dx}}AE{\cfrac {dN_{2}}{dx}}~dx+\int _{\Omega _{2}}{\cfrac {dN_{3}}{dx}}AE{\cfrac {dN_{2}}{dx}}~dx+\int _{\Omega _{3}}{\cfrac {dN_{3}}{dx}}AE{\cfrac {dN_{2}}{dx}}~dx\\K_{24}&=\int _{\Omega _{1}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{2}}{dx}}~dx+\int _{\Omega _{2}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{2}}{dx}}~dx+\int _{\Omega _{3}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{2}}{dx}}~dx\\K_{33}&=\int _{\Omega _{1}}{\cfrac {dN_{3}}{dx}}AE{\cfrac {dN_{3}}{dx}}~dx+\int _{\Omega _{2}}{\cfrac {dN_{3}}{dx}}AE{\cfrac {dN_{3}}{dx}}~dx+\int _{\Omega _{3}}{\cfrac {dN_{3}}{dx}}AE{\cfrac {dN_{3}}{dx}}~dx\\K_{34}&=\int _{\Omega _{1}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{3}}{dx}}~dx+\int _{\Omega _{2}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{3}}{dx}}~dx+\int _{\Omega _{3}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{3}}{dx}}~dx\\K_{44}&=\int _{\Omega _{1}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{4}}{dx}}~dx+\int _{\Omega _{2}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{4}}{dx}}~dx+\int _{\Omega _{3}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{4}}{dx}}~dx~.\end{aligned}}}

The matrix ${\displaystyle \mathbf {K} }$ is symmetric, so we don't need to explicitly compute the other terms.

From Figure 6, we see that ${\displaystyle N_{1}}$ is zero in elements 3 and 4, ${\displaystyle N_{2}}$ is zero in element 4, ${\displaystyle N_{3}}$ is zero in element 1, and ${\displaystyle N_{4}}$ is zero in elements 1 and 2. The same holds for ${\displaystyle dN_{i}/dx}$.

Therefore, the coefficients of the ${\displaystyle \mathbf {K} }$ matrix become

{\displaystyle {\begin{aligned}K_{11}&=\int _{\Omega _{1}}{\cfrac {dN_{1}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx~;~K_{12}=\int _{\Omega _{1}}{\cfrac {dN_{2}}{dx}}AE{\cfrac {dN_{1}}{dx}}~dx~;~K_{13}=0~;~K_{14}=0\\K_{22}&=\int _{\Omega _{1}}{\cfrac {dN_{2}}{dx}}AE{\cfrac {dN_{2}}{dx}}~dx+\int _{\Omega _{2}}{\cfrac {dN_{2}}{dx}}AE{\cfrac {dN_{2}}{dx}}~dx~;~K_{23}=\int _{\Omega _{2}}{\cfrac {dN_{3}}{dx}}AE{\cfrac {dN_{2}}{dx}}~dx~;~K_{24}=0\\K_{33}&=\int _{\Omega _{2}}{\cfrac {dN_{3}}{dx}}AE{\cfrac {dN_{3}}{dx}}~dx+\int _{\Omega _{3}}{\cfrac {dN_{3}}{dx}}AE{\cfrac {dN_{3}}{dx}}~dx~;~K_{34}=\int _{\Omega _{3}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{3}}{dx}}~dx\\K_{44}&=\int _{\Omega _{3}}{\cfrac {dN_{4}}{dx}}AE{\cfrac {dN_{4}}{dx}}~dx~.\end{aligned}}}

We can simplify our calculation further by letting ${\displaystyle N_{k}^{e}}$ be the shape functions over an element ${\displaystyle e}$. For example, the shape functions over element ${\displaystyle 2}$ are ${\displaystyle N_{1}^{2}}$ and ${\displaystyle N_{2}^{2}}$ where the local nodes ${\displaystyle 1}$ and ${\displaystyle 2}$ correspond to global nodes ${\displaystyle 2}$ and ${\displaystyle 3}$, respectively. Then we can write,

{\displaystyle {\begin{aligned}K_{11}&=\int _{\Omega _{1}}{\cfrac {dN_{1}^{1}}{dx}}AE{\cfrac {dN_{1}^{1}}{dx}}~dx~;~K_{12}=\int _{\Omega _{1}}{\cfrac {dN_{2}^{1}}{dx}}AE{\cfrac {dN_{1}^{1}}{dx}}~dx~;~K_{13}=0~;~K_{14}=0\\K_{22}&=\int _{\Omega _{1}}{\cfrac {dN_{2}^{1}}{dx}}AE{\cfrac {dN_{2}^{1}}{dx}}~dx+\int _{\Omega _{2}}{\cfrac {dN_{1}^{2}}{dx}}AE{\cfrac {dN_{1}^{2}}{dx}}~dx~;~K_{23}=\int _{\Omega _{2}}{\cfrac {dN_{2}^{2}}{dx}}AE{\cfrac {dN_{1}^{2}}{dx}}~dx~;~K_{24}=0\\K_{33}&=\int _{\Omega _{2}}{\cfrac {dN_{2}^{2}}{dx}}AE{\cfrac {dN_{2}^{2}}{dx}}~dx+\int _{\Omega _{3}}{\cfrac {dN_{1}^{3}}{dx}}AE{\cfrac {dN_{1}^{3}}{dx}}~dx~;~K_{34}=\int _{\Omega _{3}}{\cfrac {dN_{2}^{3}}{dx}}AE{\cfrac {dN_{1}^{3}}{dx}}~dx\\K_{44}&=\int _{\Omega _{3}}{\cfrac {dN_{2}^{3}}{dx}}AE{\cfrac {dN_{2}^{3}}{dx}}~dx~.\end{aligned}}}

Let ${\displaystyle K_{kl}^{e}}$ be the part of the value of ${\displaystyle K_{ij}}$ that is contributed by element ${\displaystyle e}$. The indices ${\displaystyle kl}$ are local and the indices ${\displaystyle ij}$ are global. Then,

{\displaystyle {\begin{aligned}K_{11}&=K_{11}^{1}~;~K_{12}=K_{12}^{1}~;~K_{13}=0~;~K_{14}=0\\K_{22}&=K_{22}^{1}+K_{11}^{2}~;~~K_{23}=K_{12}^{2}~;~~K_{24}=0\\K_{33}&=K_{22}^{2}+K_{11}^{3}~;~K_{34}=K_{12}^{3}\\K_{44}&=K_{22}^{3}\\\end{aligned}}}

We can therefore see that if we compute the stiffness matrices over each element and assemble them in an appropriate manner, we can get the global stiffness matrix ${\displaystyle \mathbf {K} }$.

#### Stiffness matrix for two-noded elements

For our problem, if we consider an element ${\displaystyle e}$ with two nodes, the local hat shape functions have the form

${\displaystyle N_{1}^{e}(\mathbf {x} )={\cfrac {\mathbf {x} _{2}-\mathbf {x} }{h_{e}}}~;~~N_{2}^{e}(\mathbf {x} )={\cfrac {\mathbf {x} -\mathbf {x} _{1}}{h_{e}}}}$

where ${\displaystyle h_{e}}$ is the length of the element.

Then, the components of the element stiffness matrix are

{\displaystyle {\begin{aligned}K_{11}^{e}&=\int _{x_{1}}^{x_{2}}{\cfrac {dN_{1}^{e}}{dx}}AE{\cfrac {dN_{1}^{e}}{dx}}~dx=\int _{x_{1}}^{x_{2}}\left(-{\cfrac {1}{h_{e}}}\right)AE\left(-{\cfrac {1}{h_{e}}}\right)~dx={\cfrac {AE}{h_{e}}}\\K_{12}^{e}&=\int _{x_{1}}^{x_{2}}{\cfrac {dN_{2}^{e}}{dx}}AE{\cfrac {dN_{1}^{e}}{dx}}~dx=\int _{x_{1}}^{x_{2}}\left({\cfrac {1}{h_{e}}}\right)AE\left(-{\cfrac {1}{h_{e}}}\right)~dx=-{\cfrac {AE}{h_{e}}}\\K_{22}^{e}&=\int _{x_{1}}^{x_{2}}{\cfrac {dN_{2}^{e}}{dx}}AE{\cfrac {dN_{2}^{e}}{dx}}~dx=\int _{x_{1}}^{x_{2}}\left({\cfrac {1}{h_{e}}}\right)AE\left({\cfrac {1}{h_{e}}}\right)~dx={\cfrac {AE}{h_{e}}}\end{aligned}}}

In matrix form,

${\displaystyle {\mathbf {K} ^{e}={\cfrac {AE}{h_{e}}}{\begin{bmatrix}1&-1\\-1&1\end{bmatrix}}}}$

The components of the global stiffness matrix are

{\displaystyle {\begin{aligned}K_{11}&={\cfrac {AE}{h_{1}}}~;~~K_{12}=-{\cfrac {AE}{h_{1}}}~;~~K_{13}=0~;~~K_{14}=0\\K_{22}&={\cfrac {AE}{h_{1}}}+{\cfrac {AE}{h_{2}}}~;~~K_{23}=-{\cfrac {AE}{h_{2}}}~;~~K_{24}=0\\K_{33}&={\cfrac {AE}{h_{2}}}+{\cfrac {AE}{h_{3}}}~;~~K_{34}=-{\cfrac {AE}{h_{3}}}\\K_{44}&={\cfrac {AE}{h_{3}}}\\\end{aligned}}}

In matrix form,

${\displaystyle {\mathbf {K} =AE{\begin{bmatrix}{\cfrac {1}{h_{1}}}&-{\cfrac {1}{h_{1}}}&0&0\\&{\cfrac {1}{h_{1}}}+{\cfrac {1}{h_{2}}}&-{\cfrac {1}{h_{2}}}&0\\&&{\cfrac {1}{h_{2}}}+{\cfrac {1}{h_{3}}}&-{\cfrac {1}{h_{3}}}\\{\text{Symm.}}&&&{\cfrac {1}{h_{3}}}\\\end{bmatrix}}}}$

Similarly, for the load vector ${\displaystyle \mathbf {f} }$, we have

{\displaystyle {\begin{aligned}f_{j}&=\int _{0}^{L}N_{j}\mathbf {q} ~dx+\left.N_{j}{\boldsymbol {R}}\right|_{x=L}\\&=\int _{\Omega _{1}}N_{j}\mathbf {q} ~dx+\int _{\Omega _{2}}N_{j}\mathbf {q} ~dx+\int _{\Omega _{3}}N_{j}\mathbf {q} ~dx+\left.N_{j}{\boldsymbol {R}}\right|_{x=L}\end{aligned}}}

The components of the load vector are

{\displaystyle {\begin{aligned}f_{1}&=\int _{\Omega _{1}}N_{1}\mathbf {q} ~dx+\int _{\Omega _{2}}N_{1}\mathbf {q} ~dx+\int _{\Omega _{3}}N_{1}\mathbf {q} ~dx+\left.N_{1}{\boldsymbol {R}}\right|_{x=L}\\f_{2}&=\int _{\Omega _{1}}N_{2}\mathbf {q} ~dx+\int _{\Omega _{2}}N_{2}\mathbf {q} ~dx+\int _{\Omega _{3}}N_{2}\mathbf {q} ~dx+\left.N_{2}{\boldsymbol {R}}\right|_{x=L}\\f_{3}&=\int _{\Omega _{1}}N_{3}\mathbf {q} ~dx+\int _{\Omega _{2}}N_{3}\mathbf {q} ~dx+\int _{\Omega _{3}}N_{3}\mathbf {q} ~dx+\left.N_{3}{\boldsymbol {R}}\right|_{x=L}\\f_{4}&=\int _{\Omega _{1}}N_{4}\mathbf {q} ~dx+\int _{\Omega _{2}}N_{4}\mathbf {q} ~dx+\int _{\Omega _{3}}N_{4}\mathbf {q} ~dx+\left.N_{4}{\boldsymbol {R}}\right|_{x=L}\end{aligned}}}

Once again, since ${\displaystyle N_{1}}$ is zero in elements 2 and 3, ${\displaystyle N_{2}}$ is zero in element 3, ${\displaystyle N_{3}}$ is zero in element 1, and ${\displaystyle N_{4}}$ is zero in elements 1 and 2, we have

{\displaystyle {\begin{aligned}f_{1}&=\int _{\Omega _{1}}N_{1}\mathbf {q} ~dx+\left.N_{1}{\boldsymbol {R}}\right|_{x=L}\\f_{2}&=\int _{\Omega _{1}}N_{2}\mathbf {q} ~dx+\int _{\Omega _{2}}N_{2}\mathbf {q} ~dx+\left.N_{2}{\boldsymbol {R}}\right|_{x=L}\\f_{3}&=\int _{\Omega _{2}}N_{3}\mathbf {q} ~dx+\int _{\Omega _{3}}N_{3}\mathbf {q} ~dx+\left.N_{3}{\boldsymbol {R}}\right|_{x=L}\\f_{4}&=\int _{\Omega _{3}}N_{4}\mathbf {q} ~dx+\left.N_{4}{\boldsymbol {R}}\right|_{x=L}\end{aligned}}}

Now, the boundary ${\displaystyle \mathbf {x} =L}$ is at node 4 which is attached to element 3. The only non-zero shape function at this node is ${\displaystyle N_{4}}$. Therefore, we have

{\displaystyle {\begin{aligned}f_{1}&=\int _{\Omega _{1}}N_{1}\mathbf {q} ~dx\\f_{2}&=\int _{\Omega _{1}}N_{2}\mathbf {q} ~dx+\int _{\Omega _{2}}N_{2}\mathbf {q} ~dx\\f_{3}&=\int _{\Omega _{2}}N_{3}\mathbf {q} ~dx+\int _{\Omega _{3}}N_{3}\mathbf {q} ~dx\\f_{4}&=\int _{\Omega _{3}}N_{4}\mathbf {q} ~dx+\left.N_{4}{\boldsymbol {R}}\right|_{x=L}\end{aligned}}}

In terms of element shape functions, the above equations can be written as

{\displaystyle {\begin{aligned}f_{1}&=\int _{\Omega _{1}}N_{1}^{1}\mathbf {q} ~dx=f_{1}^{1}\\f_{2}&=\int _{\Omega _{1}}N_{2}^{1}\mathbf {q} ~dx+\int _{\Omega _{2}}N_{1}^{2}\mathbf {q} ~dx=f_{2}^{1}+f_{1}^{2}\\f_{3}&=\int _{\Omega _{2}}N_{2}^{2}\mathbf {q} ~dx+\int _{\Omega _{3}}N_{1}^{3}\mathbf {q} ~dx=f_{2}^{2}+f_{1}^{3}\\f_{4}&=\int _{\Omega _{3}}N_{2}^{3}\mathbf {q} ~dx+\left.N_{2}^{3}{\boldsymbol {R}}\right|_{x=L}=f_{2}^{3}+{\boldsymbol {R}}\end{aligned}}}

The above shows that the global load vector can also be assembled from the element load vectors if we use finite element shape functions.

#### Load vector for two-noded elements

Using the linear shape functions discussed earlier and replacing ${\displaystyle \mathbf {q} }$ with ${\displaystyle a\mathbf {x} }$, the components of the element load vector ${\displaystyle \mathbf {f} ^{e}}$ are

{\displaystyle {\begin{aligned}f_{1}^{e}&=\int _{x_{1}}^{x_{2}}N_{1}^{e}a\mathbf {x} ~dx=\int _{x_{1}}^{x_{2}}\left({\cfrac {x_{2}-x}{h_{e}}}\right)a\mathbf {x} ~dx={\cfrac {a}{h_{e}}}\left({\cfrac {x_{2}(x_{2}^{2}-x_{1}^{2})}{2}}-{\cfrac {x_{2}^{3}-x_{1}^{3}}{3}}\right)\\f_{2}^{e}&=\int _{x_{1}}^{x_{2}}N_{2}^{e}a\mathbf {x} ~dx=\int _{x_{1}}^{x_{2}}\left({\cfrac {x-x_{1}}{h_{e}}}\right)a\mathbf {x} ~dx=-{\cfrac {a}{h_{e}}}\left({\cfrac {x_{1}(x_{2}^{2}-x_{1}^{2})}{2}}-{\cfrac {x_{2}^{3}-x_{1}^{3}}{3}}\right)\end{aligned}}}

In matrix form, the element load vector is written

${\displaystyle {\mathbf {f} ^{e}={\cfrac {a}{h_{e}}}{\begin{bmatrix}{\cfrac {x_{2}(x_{2}^{2}-x_{1}^{2})}{2}}-{\cfrac {x_{2}^{3}-x_{1}^{3}}{3}}\\{\cfrac {x_{2}^{3}-x_{1}^{3}}{3}}-{\cfrac {x_{1}(x_{2}^{2}-x_{1}^{2})}{2}}\end{bmatrix}}}}$

Therefore, the components of the global load vector are

{\displaystyle {\begin{aligned}f_{1}&={\cfrac {a}{h_{1}}}\left({\cfrac {x_{2}(x_{2}^{2}-x_{1}^{2})}{2}}-{\cfrac {x_{2}^{3}-x_{1}^{3}}{3}}\right)\\f_{2}&={\cfrac {a}{h_{1}}}\left({\cfrac {x_{2}^{3}-x_{1}^{3}}{3}}-{\cfrac {x_{1}(x_{2}^{2}-x_{1}^{2})}{2}}\right)+{\cfrac {a}{h_{2}}}\left({\cfrac {x_{3}(x_{3}^{2}-x_{2}^{2})}{2}}-{\cfrac {x_{3}^{3}-x_{2}^{3}}{3}}\right)\\f_{3}&={\cfrac {a}{h_{2}}}\left({\cfrac {x_{3}^{3}-x_{2}^{3}}{3}}-{\cfrac {x_{2}(x_{3}^{2}-x_{2}^{2})}{2}}\right)+{\cfrac {a}{h_{3}}}\left({\cfrac {x_{4}(x_{4}^{2}-x_{3}^{2})}{2}}-{\cfrac {x_{4}^{3}-x_{3}^{3}}{3}}\right)\\f_{4}&={\cfrac {a}{h_{3}}}\left({\cfrac {x_{4}^{3}-x_{3}^{3}}{3}}-{\cfrac {x_{3}(x_{4}^{2}-x_{3}^{2})}{2}}\right)+{\boldsymbol {R}}\end{aligned}}}

### Displacement trial function

Recall that we assumed that the displacement can be written as

${\displaystyle \mathbf {u} _{h}(\mathbf {x} )=a_{1}\varphi _{1}(\mathbf {x} )+a_{2}\varphi _{2}(\mathbf {x} )+\dots +a_{n}\varphi _{n}(\mathbf {x} )=\sum _{i=1}^{n}a_{i}\varphi _{i}(\mathbf {x} )~.}$

If we use finite element shape functions, we can write the above as

${\displaystyle \mathbf {u} _{h}(\mathbf {x} )=a_{1}~N_{1}(\mathbf {x} )+a_{2}~N_{2}(\mathbf {x} )+\dots +a_{n}~N_{n}(\mathbf {x} )=\sum _{i=1}^{n}a_{i}~N_{i}(\mathbf {x} )}$

where ${\displaystyle n}$ is the total number of nodes in the domain. Also, recall that the value of ${\displaystyle N_{i}}$ is 1 at node ${\displaystyle i}$ and zero elsewhere. Therefore, we have

{\displaystyle {\begin{aligned}u_{1}&:=\mathbf {u} _{h}(\mathbf {x} _{1})=a_{1}~N_{1}(\mathbf {x} _{1})+a_{2}~N_{2}(\mathbf {x} _{1})+\dots +a_{n}~N_{n}(\mathbf {x} _{1})=a_{1}\\u_{2}&:=\mathbf {u} _{h}(\mathbf {x} _{2})=a_{1}~N_{1}(\mathbf {x} _{2})+a_{2}~N_{2}(\mathbf {x} _{2})+\dots +a_{n}~N_{n}(\mathbf {x} _{2})=a_{2}\\\vdots \\u_{n}&:=\mathbf {u} _{h}(\mathbf {x} _{n})=a_{1}~N_{1}(\mathbf {x} _{n})+a_{2}~N_{2}(\mathbf {x} _{n})+\dots +a_{n}~N_{n}(\mathbf {x} _{n})=a_{n}\\\end{aligned}}}

Therefore, the trial function can be written as

${\displaystyle \mathbf {u} _{h}(\mathbf {x} )=u_{1}~N_{1}(\mathbf {x} )+u_{2}~N_{2}(\mathbf {x} )+\dots +u_{n}~N_{n}(\mathbf {x} )=\sum _{i=1}^{n}u_{i}~N_{i}(\mathbf {x} )}$

where ${\displaystyle u_{i}}$ are the nodal displacements.

### Finite element system of equations

If all the elements are assumed to be of the same length ${\displaystyle h}$, the finite element system of equations (${\displaystyle \mathbf {K} \mathbf {a} =\mathbf {f} }$) can then be written as

${\displaystyle {\cfrac {AE}{h}}{\begin{bmatrix}1&-1&0&0\\-1&2&-1&0\\0&-1&2&-1\\0&0&-1&1\end{bmatrix}}{\begin{bmatrix}u_{1}\\u_{2}\\u_{3}\\u_{4}\end{bmatrix}}={\cfrac {a}{h}}{\begin{bmatrix}\left({\cfrac {x_{2}(x_{2}^{2}-x_{1}^{2})}{2}}-{\cfrac {x_{2}^{3}-x_{1}^{3}}{3}}\right)\\\left({\cfrac {x_{2}^{3}-x_{1}^{3}}{3}}-{\cfrac {x_{1}(x_{2}^{2}-x_{1}^{2})}{2}}\right)+\left({\cfrac {x_{3}(x_{3}^{2}-x_{2}^{2})}{2}}-{\cfrac {x_{3}^{3}-x_{2}^{3}}{3}}\right)\\\left({\cfrac {x_{3}^{3}-x_{2}^{3}}{3}}-{\cfrac {x_{2}(x_{3}^{2}-x_{2}^{2})}{2}}\right)+\left({\cfrac {x_{4}(x_{4}^{2}-x_{3}^{2})}{2}}-{\cfrac {x_{4}^{3}-x_{3}^{3}}{3}}\right)\\\left({\cfrac {x_{4}^{3}-x_{3}^{3}}{3}}-{\cfrac {x_{3}(x_{4}^{2}-x_{3}^{2})}{2}}\right)+{\boldsymbol {R}}{\cfrac {h}{a}}\end{bmatrix}}}$

### Essential boundary conditions

To solve this system of equations we have to apply the essential boundary condition ${\displaystyle \mathbf {u} =0}$ at ${\displaystyle \mathbf {x} =0}$. This is equivalent to setting ${\displaystyle u_{1}=0}$. The reduced system of equations is

${\displaystyle {\cfrac {AE}{h}}{\begin{bmatrix}2&-1&0\\-1&2&-1\\0&-1&1\end{bmatrix}}{\begin{bmatrix}u_{2}\\u_{3}\\u_{4}\end{bmatrix}}={\cfrac {a}{h}}{\begin{bmatrix}\left({\cfrac {x_{2}^{3}-x_{1}^{3}}{3}}-{\cfrac {x_{1}(x_{2}^{2}-x_{1}^{2})}{2}}\right)+\left({\cfrac {x_{3}(x_{3}^{2}-x_{2}^{2})}{2}}-{\cfrac {x_{3}^{3}-x_{2}^{3}}{3}}\right)\\\left({\cfrac {x_{3}^{3}-x_{2}^{3}}{3}}-{\cfrac {x_{2}(x_{3}^{2}-x_{2}^{2})}{2}}\right)+\left({\cfrac {x_{4}(x_{4}^{2}-x_{3}^{2})}{2}}-{\cfrac {x_{4}^{3}-x_{3}^{3}}{3}}\right)\\\left({\cfrac {x_{4}^{3}-x_{3}^{3}}{3}}-{\cfrac {x_{3}(x_{4}^{2}-x_{3}^{2})}{2}}\right)+{\boldsymbol {R}}{\cfrac {h}{a}}\end{bmatrix}}}$

This system of equations can be solved for ${\displaystyle u_{2}}$, ${\displaystyle u_{3}}$, and ${\displaystyle u_{4}}$. Let us do that.

Assume that ${\displaystyle A}$, ${\displaystyle E}$, ${\displaystyle L}$, ${\displaystyle a}$, and ${\displaystyle R}$ are all equal to 1. Then ${\displaystyle x_{1}=0}$, ${\displaystyle x_{2}=1/3}$, ${\displaystyle x_{3}=2/3}$, ${\displaystyle x_{4}=1}$, and ${\displaystyle h=1/3}$. The system of equations becomes

${\displaystyle {\begin{bmatrix}2&-1&0\\-1&2&-1\\0&-1&1\end{bmatrix}}{\begin{bmatrix}u_{2}\\u_{3}\\u_{4}\end{bmatrix}}={\begin{bmatrix}0.037\\0.074\\0.383\end{bmatrix}}\qquad \implies \qquad {\begin{bmatrix}u_{1}\\u_{2}\\u_{3}\\u_{4}\end{bmatrix}}={\begin{bmatrix}0.0\\0.494\\0.951\\1.333\end{bmatrix}}}$

### Computing element strains and stresses

From the above, it is clear that the displacement field within an element ${\displaystyle e}$ is given by

${\displaystyle \mathbf {u} ^{e}=u_{1}^{e}N_{1}^{e}(\mathbf {x} )+u_{2}^{e}N_{2}^{e}(\mathbf {x} )~.}$

Therefore, the strain within an element is

${\displaystyle {\boldsymbol {\varepsilon }}^{e}={\frac {\partial \mathbf {u} ^{e}}{\partial x}}=u_{1}^{e}{\frac {\partial N_{1}^{e}}{\partial x}}+u_{2}^{e}{\frac {\partial N_{2}^{e}}{\partial x}}~.}$

In matrix notation,

${\displaystyle {\boldsymbol {\varepsilon }}^{e}=\mathbf {B} ^{e}\mathbf {u} ^{e}=\left[{\frac {\partial N_{1}^{e}}{\partial x}}~~{\frac {\partial N_{2}^{e}}{\partial x}}\right]{\begin{bmatrix}u_{1}^{e}\\u_{2}^{e}\end{bmatrix}}=\left[-{\cfrac {1}{h}}~~{\cfrac {1}{h}}\right]{\begin{bmatrix}u_{1}^{e}\\u_{2}^{e}\end{bmatrix}}}$

The stress in the element is given by

${\displaystyle {\boldsymbol {\sigma }}^{e}=E{\boldsymbol {\varepsilon }}^{e}~.}$

For our discretization, the element stresses are

{\displaystyle {\begin{aligned}{\boldsymbol {\sigma }}^{1}&=1.48\\{\boldsymbol {\sigma }}^{2}&=1.37\\{\boldsymbol {\sigma }}^{3}&=1.15\end{aligned}}}

A plot of this solution is shown in Figure 7.

 Figure 7(a). FEM vs exact solutions for displacements of an axially loaded bar.
 Figure 7(b). FEM vs exact solutions for stresses in an axially loaded bar.

#### Matlab code

The finite element code (Matlab) used to compute this solution is given below.

 function AxialBarFEM A = 1.0; L = 1.0; E = 1.0; a = 1.0; R = 1.0; e = 3; h = L/e; n = e+1; for i=1:n node(i) = (i-1)*h; end for i=1:e elem(i,:) = [i i+1]; end K = zeros(n); f = zeros(n,1); for i=1:e node1 = elem(i,1); node2 = elem(i,2); Ke = elementStiffness(A, E, h); fe = elementLoad(node(node1),node(node2), a, h); K(node1:node2,node1:node2) = K(node1:node2,node1:node2) + Ke; f(node1:node2) = f(node1:node2) + fe; end f(n) = f(n) + R; Kred = K(2:n,2:n); fred = f(2:n); d = inv(Kred)*fred; dsol = [0 d']; fsol = K*dsol'; sum(fsol) figure; p0 = plotDisp(E, A, L, R, a); p1 = plot(node, dsol, 'ro--', 'LineWidth', 3); hold on; legend([p0 p1],'Exact','FEM'); for i=1:e node1 = elem(i,1); node2 = elem(i,2); u1 = dsol(node1); u2 = dsol(node2); [eps(i), sig(i)] = elementStrainStress(u1, u2, E, h); end figure; p0 = plotStress(E, A, L, R, a); for i=1:e node1 = node(elem(i,1)); node2 = node(elem(i,2)); p1 = plot([node1 node2], [sig(i) sig(i)], 'r-','LineWidth',3); hold on; end legend([p0 p1],'Exact','FEM'); function [p] = plotDisp(E, A, L, R, a) dx = 0.01; nseg = L/dx; for i=1:nseg+1 x(i) = (i-1)*dx; u(i) = (1/(6*A*E))*(-a*x(i)^3 + (6*R + 3*a*L^2)*x(i)); end p = plot(x, u, 'LineWidth', 3); hold on; xlabel('x', 'FontName', 'palatino', 'FontSize', 18); ylabel('u(x)', 'FontName', 'palatino', 'FontSize', 18); set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18); function [p] = plotStress(E, A, L, R, a) dx = 0.01; nseg = L/dx; for i=1:nseg+1 x(i) = (i-1)*dx; sig(i) = (1/(2*A))*(-a*x(i)^2 + (2*R + a*L^2)); end p = plot(x, sig, 'LineWidth', 3); hold on; xlabel('x', 'FontName', 'palatino', 'FontSize', 18); ylabel('\sigma(x)', 'FontName', 'palatino', 'FontSize', 18); set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18); function [Ke] = elementStiffness(A, E, h) Ke = (A*E/h)*[[1 -1];[-1 1]]; function [fe] = elementLoad(node1, node2, a, h) x1 = node1; x2 = node2; fe1 = a*x2/(2*h)*(x2^2-x1^2) - a/(3*h)*(x2^3-x1^3); fe2 = -a*x1/(2*h)*(x2^2-x1^2) + a/(3*h)*(x2^3-x1^3); fe = [fe1;fe2]; function [eps, sig] = elementStrainStress(u1, u2, E, h) B = [-1/h 1/h]; u = [u1; u2]; eps = B*u sig = E*eps;