# Application of Finite Layer Method in Pavement Structural Analysis
## 2 Numerical Solution of EasyFEM
### FE Solution of Displacements
\begin{array}{l}[B_{1}^{mn}\ B_{2}^{mn}\cdots B_{k+1}^{mn}]=\left[ \begin{array}{cc}E^{mn}N_{1}(y)&E^{mn}N_{2}(y)&\cdots E^{mn}N_{k+1}(y)\end{array} \right]\end{array}  \tag{6}\]  
is called the strain-displacement matrix. The relation between stresses and strains is used to express the total potential energy with nodal displacements \(\{d\}\):  
\[\Pi(\{d\})=\int_{V}\frac{1}{2}([B]\{d\})^{T}[D]([B]\{d\})dV-\int_{V}([N]\{d\})^ {T}\{b\}dV-\int_{S}([N]\{d\})^{T}\{f\}dS. \tag{7}\]  
The nodal displacements \(\{d\}\) which correspond to the minimum of the functional \(\Pi\) are determined by the condition that the first variation of \(\Pi\) should be zero, i.e.,  
\[\delta\Pi(\{d\})=0, \tag{8}\] \[\delta\Pi(\{d\})=\int_{V}\{\delta d\}^{T}[B]^{T}[D][B]\{d\}dV-\int_ {V}\{\delta d\}^{T}[N]^{T}\{b\}dV-\int_{S}\{\delta d\}^{T}[N]^{T}\{f\}dS\.\]  
Then, the arbitrariness of \(\{\delta d\}\) leads the following equilibrium equation for FE:  
\[\int_{V}[B]^{T}[D][B]dV\{d\}-\int_{V}[N]^{T}\{b\}dV-\int_{S}[N]^{T}\{f\}dS\ =0, \tag{9}\]  
which is commonly given in the following form:  
\[[k]\{d\}=\{f\}, \tag{10}\] \[[k]=\int_{V}[B]^{T}[D][B]dV\ \text{and}\ \{f\}=\int_{V}[N]^{T}\{b\}dV+\int_{S}[N]^{T}\{f\}dS,\,begin\]  
where \([k]\) is the element stiffness matrix; and \(\{f\}\) is the load vector. From Equations (4) and (10), the stiffness matrix of one element includes:  
\[ \begin{array}{l}I_{1}=\int_{0}^{a}\int_{0}^{c}\cos\frac{m\pi x}{a}\sin\frac{ n\pi z}{c}\sin\frac{p\pi x}{a}\sin\frac{q\pi z}{c}\cdot dxdz,\\ I_{2}=\int_{0}^{a}\int_{0}^{c}\cos\frac{m\pi x}{a}\sin\frac{n\pi z}{c}\cos\frac{ p\pi x}{a}\sin\frac{q\pi z}{c}\cdot dxdz,\\ I_{3}=\int_{0}^{a}\int_{0}^{c}\cos\frac{m\pi x}{a}\cos\frac{n\pi z}{c}\cos\frac{ p\pi x}{a}\cos\frac{q\pi z}{c}\cdot dxdz.\end{array}  \tag{11}\]  
The integrals exhibit orthogonal properties, which ensure that:  
\[I_{1}=I_{2}=I_{3}=\left\{ \begin{array}{l}\frac{ac}{4},\ \ for\ m=p\ and\ n=q\\ 0,\ \ for\ m\
eq p\ or\ n\
eq q\end{array} \right\}. \tag{12}\]This means that the matrix \((k^{mnpq})^{\rm e}\) is diagonal. In other words, the non-zero values are only located on the diagonal, where \(m=p\) and \(n=q\). Thus, the stiffness matrix can be reduced to:  
\[\left(k^{mumm}_{gk}\right)^{e}=\frac{ac}{4}\int_{length}\left(B^{mnT}_{g}D^{bm }_{k}\right)dy\qquad m=1,\ 2\ldots M;\ n=1,2\ldots N, \tag{13}\]  
where \(g\) and \(k\) represent the nodes of the element, respectively. The length corresponds to the length of the element. A typical term for the force vector becomes:  
\[(f^{mn})^{e}=\int\int\int_{vol}(N^{mn})^{T}\{p\}^{mm}dxdydz. \tag{14}\]  
By assembling the stiffness matrix of each element to the global domain, the global linear system is achieved as follows:  
\[\left[ \begin{array}{cccc}K^{1111}&&&\\ &K^{1212}&&\\ &&\ddots&\\ &&&K^{MNMN}\end{array} \right]\left\{ \begin{array}{c}U^{11}\\ U^{12}\\ \vdots\\ U^{MN}\end{array} \right\}+\left\{ \begin{array}{c}F^{11}\\ F^{12}\\ \vdots\\ F^{MN}\end{array} \right\}=0. \tag{15}\]  
Equation (15) shows that the large system of equations splits up into \(M\times N\) separate problems:  
\[K^{mumm}U^{mn}+F^{mn}=0, \tag{16}\]  
where \(K\) is the global stiffness matrix; \(U\) is the global displacement vector; and \(F\) is the global loading vector.