## 3 Stress Analysis by the HDQM
\begin{split}\frac{E_{1}^{b}f_{b}^{3}}{12}\sum\limits_{k=2}^{N-1} D_{lk}^{x}w_{1k}-Q_{13}\sum\limits_{k=2}^{N-1}A_{lk}^{x}u_{1k}-Q_{33}\sum \limits_{j=1}^{M}A_{1j}^{z}w_{jl}\\ -\frac{f_{b}Q_{66}}{2}\left(\sum\limits_{k=2}^{N-1}B_{lk}^{x}w_{ 1k}+\sum\limits_{j=1}^{M}\sum\limits_{k=2}^{N-1}A_{1j}^{z}A_{lk}^{x}u_{jk} \right)=0\\ (l=2,3,...,N-1)\end{split}  \tag{16}\]  
where  
\[C_{lk}^{x}=\sum\limits_{j=1}^{N}A_{lj}^{x}\widetilde{B}_{jk}^{x},\ D_{lk}^{x} =\sum\limits_{j=1}^{N}B_{lj}^{x}\widetilde{B}_{jk}^{x} \tag{17}\]  
Note that boundary conditions \(u_{1j}=w_{1j}=0\) (\(j=1\), \(N\)) have been applied. Besides, Equations (15) and (16) are not applied at the two corner points (\(l=1\), \(N\)). In formulating \(C_{lk}^{x}\) and \(D_{lk}^{x}\), \(w_{x}=0\) at nodes 1 and \(N\) has been built-in by using MMWC-3 [16; 32]. The MMWC-3 is very simple and only \(\widetilde{B}_{jk}^{x}\) is computed slightly different from Equation (14), namely,  
\[\left\{ \begin{array}{l}\widetilde{B}_{jk}^{x}=\sum\limits_{l=2}^{N-1}A_{jl }^{x}A_{lk}^{x}\ (i=1,N;j=1,2,...,N)\\ \widetilde{B}_{jk}^{x}=B_{jk}^{x}\ (i=2,3,...,N-1;j=1,2,...,N)\end{array} \right. \tag{18}\]  
Similarly, in terms of harmonic differential quadrature, Equations (5) and (6) are expressed as  
\[ \begin{split} f_{t}E_{1}^{t}\left(\sum\limits_{k=2}^{N-1}B_{lk}^ {x}u_{Mk}-0.5f_{t}\sum\limits_{k=2}^{N-1}C_{lk}^{x}w_{Mk}\right)-Q_{66}\left( \sum\limits_{k=2}^{N-1}A_{lk}^{x}w_{Mk}+\sum\limits_{j=1}^{M}A_{Mj}^{z}u_{jl }\right)=0\\ (l=2,3,...,N-1)\end{split}  \tag{19}\]  
and  
\[ \begin{split}\frac{E_{1}^{t}f_{t}^{3}}{12}\sum\limits_{k=2}^{N-1 }D_{lk}^{x}w_{Mk}+Q_{13}\sum\limits_{k=2}^{N-1}A_{lk}^{x}u_{Mk}+Q_{33}\sum \limits_{j=1}^{M}A_{Mj}^{z}w_{jl}\\ -\frac{f_{t}Q_{66}}{2}\left(\sum\limits_{k=2}^{N-1}B_{lk}^{x}w_{Mk}+ \sum\limits_{j=1}^{M}\sum\limits_{k=2}^{N-1}A_{Mj}^{z}A_{lk}^{x}u_{jk}\right)= q_{Ml}\end{split}  \tag{20}\]  
\[ \begin{split}(l=2,3,...,N-1)\end{split} \]where  
\[q_{MI}=P\sum\limits_{k=2}^{N-1}l_{l}(L_{p}\xi_{k}/2)H_{k}/H_{l}L \tag{21}\]  
Note that boundary conditions \(u_{Mj}=w_{Mj}=0\) (\(j=1\), \(N\)) have been applied. Besides, Equations (19) and (20) are not applied at the two corner points (\(l=1\), \(N\)) and Equation (21) is obtained by using the method of work equivalent load. The method is similar to the one proposed by the senior author [16] (pp. 98-99, 125-127) but is more rigorous. Essential to the method is multiplying both sides of Equation (6) by \(w^{\varepsilon}(x)/2\) and then integrating the resultant equation from \(-L/2\) to \(L/2\) numerically by Gauss quadrature. Differential equations at all integration points are then obtained. Detailed derivations for obtaining the work equivalent load \(q_{MI}\) are included in the Appendix A.  
Equations (11), (15), (16), (19) and (20) contain \(2(N-2)\times M\) algebraic equations, which can be written in following matrix form,  
\[\left[K\right]\left\{U\right\}=\left\{F\right\} \tag{22}\]  
Or in a partitioned form as  
\[\left[ \begin{array}{cc}\left[k_{uu}\right]&\left[k_{uu}\right]\\ \left[k_{uu}\right]&\left[k_{ww}\right]\end{array} \right]\left\{ \begin{array}[ ]{c}\left\{u\right\}\\ \left\{w\right\}\end{array} \right\}=\left\{ \begin{array}{c}\left\{f_{u} \right\}\\ \left\{f_{w}\right\}\end{array} \right\} \tag{23}\]  
where \(\left[K\right]\) is a non-symmetric full matrix. Since the fixed boundary conditions have been applied already, displacement vector \(\left\{U\right\}\) can be obtained as  
\[\left\{U\right\}=\left[K\right]^{-1}\left\{F\right\} \tag{24}\]  
Once the displacement vector \(\left\{U\right\}\) is obtained, the stress components at any grid point of the core can be calculated by  
\[\left\{