## 3 Higher-Order Equivalent Layer-Wise Approach
It should be noticed that \(\mathbf{A}^{(\tau\eta)}\) represents the stiffness matrix of the considered shell structure. In particular, if a laminated composite structure made of \(l\) orthotropic layers is considered, the constitutive operator takes the following aspect  
\[\mathbf{A}^{(\tau\eta)}=\left[ \begin{array}{cccccc}A_{11(20)}^{(\tau\eta)}&A _{12(11)}^{(\tau\eta)}&A_{16(20)}^{(\tau\eta)}&A_{16(11)}^{(\tau\eta)}&0&0&0&A_{ 13(10)}^{(\tau\bar{\eta})}\\ A_{12(11)}^{(\tau\eta)}&A_{22(02)}^{(\tau\eta)}&A_{26(11)}^{(\tau\eta)}&A_{26(02 )}^{(\tau\eta)}&0&0&0&A_{23(01)}^{(\tau\bar{\eta})}\\ A_{16(20)}^{(\tau\eta)}&A_{26(11)}^{(\tau\eta)}&A_{66(20)}^{(\tau\eta)}&A_{66(11 )}^{(\tau\eta)}&0&0&0&0&A_{36(10)}^{(\tau\bar{\eta})}\\ A_{16(11)}^{(\tau\eta)}&A_{26(02)}^{(\tau\eta)}&A_{66(11)}^{(\tau\eta)}&A_{66(02 )}^{(\tau\eta)}&0&0&0&0&A_{36(01)}^{(\tau\bar{\eta})}\\ 0&0&0&0&A_{44(20)}^{(\tau\eta)}&A_{45(11)}^{(\tau\eta)}&A_{44(10)}^{(\tau\bar{ \eta})}&A_{45(10)}^{(\tau\bar{\eta})}&0\\ 0&0&0&0&A_{45(11)}^{(\tau\eta)}&A_{55(02)}^{(\tau\eta)}&A_{45(01)}^{(\tau\bar{ \eta})}&A_{55(01)}^{(\tau\bar{\eta})}&0\\ 0&0&0&0&A_{44(10)}^{(\bar{\tau}\eta)}&A_{45(01)}^{(\tau\eta)}&A_{44(00)}^{(\tau \bar{\eta})}&A_{45(00)}^{(\tau\bar{\eta})}&0\\ 0&0&0&0&A_{45(10)}^{(\tau\eta)}&A_{55(01)}^{(\tau\eta)}&A_{45(00)}^{(\tau\bar{ \eta})}&0\\ A_{13(10)}^{(\bar{\tau}\eta)}&A_{23(01)}^{(\bar{\tau}\eta)}&A_{36(10)}^{(\bar {\tau}\eta)}&A_{36(01)}^{(\bar{\tau}\eta)}&0&0&0&0&A_{33(00)}^{(\tau\bar{\eta}) }\end{array} \right] \tag{29}\]  
Each term collected in the stiffness matrix can be evaluated by means of the following definitions  
\[ \begin{array}{c}A_{nm\,(pq)}^{(\tau\eta)}=\sum\limits_{k=1}^{l}\int\limits_ {\xi_{k}}^{\xi_{k+1}}\overline{B}_{nm}^{(k)}F_{\eta}F_{\tau}\frac{H_{1}H_{2}} {H_{1}^{p}H_{2}^{d}}d\zeta\\ A_{nm\,(pq)}^{(\bar{\tau}\eta)}=\sum\limits_{k=1}^{l}\int\limits_{\xi_{k}}^{\xi_{ k+1}}\overline{B}_{nm}^{(k)}F_{\eta}\frac{\partial F_{\tau}}{\partial\zeta} \frac{H_{1}H_{2}}{H_{1}^{p}H_{2}^{d}}d\zeta\\ A_{nm\,(pq)}^{(\tau\bar{\eta})}=\sum\limits_{k=1}^{l}\int\limits_{\xi_{k}}^{\xi_{ k+1}}\overline{B}_{nm}^{(k)}\frac{\partial F_{\eta}}{\partial\zeta}F_{\tau}\frac{H_{1}H_{ 2}}{H_{1}^{p}H_{2}^{d}}d\zeta\\ A_{nm\,(pq)}^{(\bar{\tau}\bar{\eta})}=\sum\limits_{k=1}^{l}\int\limits_{\xi_{k}}^{ \xi_{k+1}}\overline{B}_{nm}^{(k)}\frac{\partial F_{\eta}}{\partial\zeta}\frac{ \partial F_{\tau}}{\partial\zeta}\frac{H_{1}H_{2}}{H_{1}^{p}H_{2}^{d}}d\zeta \end{array}  \tag{30}\]for \(\tau,\eta=0,1,2,\ldots,N,N+1\), \(n,m=1,2,3,4,5,6\), and \(p,q=0,1,2\). It should be noticed that the symbols \(\tau,\eta\) specify the thickness function order, whereas \(\widetilde{\tau},\widetilde{\eta}\) denote that the derivatives of the corresponding thickness functions \(F_{\tau}\), \(F_{\eta}\) with respect to \(\zeta\) must be evaluated. It is pointed out that this evaluation can be easily performed since the definition of the thickness functions is shown in extended form in (12) and (18), respectively. Quantities \(\overline{B}_{nm}^{(k)}\) are defined as follows  
\[\overline{B}_{nm}^{(k)}=\overline{E}_{nm}^{(k)} \tag{31}\]  
for \(n,m=1,2,3,6\) and  
\[\overline{B}_{nm}^{(k)}=\kappa\overline{E}_{nm}^{(k)} \tag{32}\]  
for \(n,m=4,5\). This distinction is performed to introduce the shear correction factor \(\kappa=1/\chi\) that has to be included in the model if the chosen order of kinematic expansion does not provide a parabolic profile of the through-the-thickness shear stresses. Typically, a constant value is used to define the shear correction factor \(\kappa\). In the present paper, the value \(\chi=1.2\) is used for this purpose. Nevertheless, it should be noticed that a parabolic function can be chosen as shear correction function as illustrated in the work by [PERSON] and [PERSON] [92]. On the other hand, quantities \(\overline{E}_{nm}^{(k)}\), for \(n,m=1,2,3,4,5,6\), represent the elastic constant of the material of the \(k\)-th layer. They can assume different meanings depending on the case.