## 3 Stress Analysis by the HDQM  
The system of differential equations is solved numerically by using the HDQM. Hence, the formulations in terms of the harmonic differential quadrature are given in this section. The numbers of grid points in \(x\) and \(z\) direction are denotes by \(N\) and \(M\). Since the method of modification of weighting coefficient-3 (MMWC-3) [16; 32] is used to implement \(w_{,x}(\pm L/2,\pm c)=0\) for simplicity, grid V is used in \(x\)-direction for reliability consideration. It was observed that the spurious complex eigenvalues of the modified weighting coefficient matrix did not exist if grid V was used [19]. The non-dimensional form of grid V is defined as  
\[\xi_{1}=-1,\ \xi_{i}\ =\ (N-2)\ \ \text{Gauss points}\ (i=2,3,...,N-1),\ \xi_{N}=1 \tag{9}\]  
The actual grid coordinates in \(x\)-direction are \(x_{i}=L\xi_{i}/2\), \((i=1,2,...,N)\). For convenience and accuracy considerations, grid III is used in \(z\)-direction. Its non-dimensional form is defined as  
\[\eta_{j}=-\cos\left[(j-1)\pi/(M-1)\right](j=1,2,...,M) \tag{10}\]  
The actual grid coordinates in \(z\)-direction are \(z_{j}=c\eta_{j}\), \((j=1,2,...,M)\). In terms of harmonic differential quadrature, Equation (2) at all inner grid points are given by  
\[\left\{ \begin{array}{l}Q_{11}\sum\limits_{k=2}^{N-1}B_{lk}^{x}u_{ik}+Q_{66} \sum\limits_{j=1}^{M}B_{ij}^{z}u_{jl}+(Q_{13}+Q_{66})\sum\limits_{j=1}^{M} \sum\limits_{k=2}^{N-1}A_{ij}^{z}A_{lk}^{x}w_{jk}=0\\ Q_{66}\sum\limits_{k=2}^{N-1}B_{lk}^{x}w_{ik}+Q_{33}\sum\limits_{j=1}^{M}B_{ij }^{z}w_{jl}+(Q_{13}+Q_{66})\sum\limits_{j=1}^{M}\sum\limits_{k=2}^{N-1}A_{ij}^ {z}A_{lk}^{x}u_{jk}=0\\ \\ (i=2,3,...,M-1;l=2,3,...,N-1)\end{array} \right. \tag{11}\]  
where \(A_{lk}^{x}\) and \(A_{ij}^{z}\) are weighting coefficients of the first-order derivative with respect to the variable of \(x\) or \(z\), \(B_{lk}^{x}\) and \(B_{ij}^{z}\) are weighting coefficients of the second-order derivative with respect to the variable of \(x\) or \(z\), respectively. The fixed boundary conditions, i.e., \(u(\pm L/2,z)=w(\pm L/2,z)=0\), have been applied. Explicit formulas to compute the weighting coefficients are available [16; 23; 24]. \(A_{lk}^{x}\) and \(A_{ij}^{z}\) can be computed by  
\[A_{lk}^{x}=\left\{ \begin{array}{ll}\dfrac{2p}{L}\prod\limits_{i=1}^{N}\sin p (\xi_{l}-\xi_{i})\left/\prod\limits_{i\
eq k}^{N}\sin p(\xi_{k}-\xi_{i})&(l \
eq k)\\ \\ \dfrac{2p}{L}\sum\limits_{i\
eq l}^{N}\dfrac{\cos p(\xi_{l}-\xi_{i})}{\sin p( \xi_{l}-\xi_{i})}&(l=k)\end{array} \right. \tag{12}\]  
and  
\[A_{ij}^{z}=\left\{ \begin{array}{ll}\dfrac{p}{c}\prod\limits_{k=l}^{M}\sin p (\eta_{i}-\eta_{k})\left/\prod\limits_{k=l}^{M}\sin p(\eta_{j}-\eta_{k})&(i \
eq j)\\ \\ \dfrac{p}{c}\sum\limits_{k=l}^{M}\dfrac{\cos p(\eta_{i}-\eta_{k})}{\sin p( \eta_{i}-\eta_{k})}&(i=j)\end{array} \right. \tag{13}\]where \(\underline{k},\eta\in[-1,1]\) and \(p\) can be 0.5 [23], \(\pi/8\)[16], \(\pi/4\)[24; 26; 31] or any other values, usually less than 1.0. In present study, \(p\) = 0.25. It is not difficult to show that the weighting coefficients are approximately the same as the ones of the DQM when \(N\) and \(M\) are large enough or \(p\) takes a very small number. \(B_{lk}^{x}\) and \(B_{ij}^{z}\) can be conveniently computed by  
\[B_{lk}^{x}=\sum\limits_{j=1}^{N}A_{lj}^{x}A_{jk}^{x},\ B_{ij}^{z}=\sum\limits_ {k=1}^{M}A_{ik}^{z}A_{kj}^{z} \tag{14}\]  
In terms of harmonic differential quadrature, Equations (3) and (4) are expressed by  
\[ \begin{split} f_{b}E_{1}^{b}\left(\sum\limits_{k=2}^{N-1}B_{lk}^ {x}u_{1k}+0.5f_{b}\sum\limits_{k=2}^{N-1}C_{lk}^{x}w_{1k}\right)+Q_{66}\left( \sum\limits_{k=2}^{N-1}A_{lk}^{x}w_{1k}+\sum\limits_{j=1}^{M}A_{ij}^{z}u_{jl }\right)=0\\ (l=2,3,...,N-1)\end{split}  \tag{15}\]  
and  
\[