## 2 Theoretical Model
### Shell Formulation
\begin{array}{cccccccc}\frac{F_{\tau}}{H_{1}}&0&0 &0&0&0&0&0&0\\ 0&\frac{F_{\tau}}{H_{2}}&0&0&0&0&0&0&0\\ 0&0&\frac{F_{\tau}}{H_{1}}&\frac{F_{\tau}}{H_{2}}&0&0&0&0&0\\ 0&0&0&0&\frac{F_{\tau}}{H_{1}}&0&\frac{\partial F_{\tau}}{\partial\zeta}&0&0\\ 0&0&0&0&0&\frac{F_{\tau}}{H_{2}}&0&\frac{\partial F_{\tau}}{\partial\zeta}&0\\ 0&0&0&0&0&0&0&\frac{\partial F_{\tau}}{\partial\zeta}\end{array} \right] \tag{17}\]  
The stress results for each order \(\tau\) of kinematic expansion can be calculated over the generalized strains (11) as  
\[\mathbf{S}^{(\tau)}=\sum_{\eta=0}^{N+1}\mathbf{A}^{(\tau\eta)}\varepsilon^{(\eta)} \tag{18}\]  
where the algebraic vector \(\mathbf{S}^{(\tau)}=\mathbf{S}^{(\tau)}(\alpha_{1},\alpha_{2},t)\) is taken as  
\[\mathbf{S}^{(\tau)}=\left[ \begin{array}{cccccccc}N_{1}^{(\tau)}&N_{2}^{( \tau)}&N_{12}^{(\tau)}&N_{21}^{(\tau)}&T_{1}^{(\tau)}&T_{2}^{(\tau)}&P_{1}^{( \tau)}&P_{2}^{(\tau)}&S_{3}^{(\tau)}\end{array} \right]^{T} \tag{19}\]  
The constitutive operator \(\mathbf{A}^{(\tau\eta)}\) is defined as follows  
\[\mathbf{A}^{(\tau\eta)}=\left[ \begin{array}{cccccccc}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&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&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\bar{\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)}^{(\bar{\tau}\eta)}&A_{44(00)}^{( \tau\bar{\eta})}&0\\ 0&0&0&0&A_{45(10)}^{(\bar{\tau}\eta)}&A_{55(01)}^{(\bar{\tau}\eta)}&A_{45(00)}^{( \bar{\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)}^{(\bar{\tau} \bar{\eta})}\\ \end{array} \right] \tag{20}\]where the elastic constant in (20) can be evaluated according to the following expressions  
\[ \begin{array}{l}A^{(\tau\eta)}_{nm\ (pq)}=\sum\limits_{k=1}^{l}\int\limits_{ \xi_{k}}^{\xi_{k+1}}\overline{B}^{(k)}_{nm}F_{\eta}F_{\tau}\frac{H_{1}H_{2}}{H_ {1}^{2}H_{2}^{2}}d\zeta\\ A^{(\bar{\tau}\eta)}_{nm\ (pq)}=\sum\limits_{k=1}^{l}\int\limits_{\xi_{k}}^{\xi_{k+1}} \overline{B}^{(k)}_{nm}F_{\eta}\frac{\partial F_{\tau}}{\partial\xi}\frac{H_{1 }H_{2}}{H_{1}^{2}H_{2}^{2}}d\zeta\\ A^{(\tau\bar{\eta})}_{nm\ (pq)}=\sum\limits_{k=1}^{l}\int\limits_{\xi_{k}}^{\xi_{k+1}} \overline{B}^{(k)}_{nm}\frac{\partial F_{\tau}}{\partial\xi}F_{\tau}\frac{H_{ 1}H_{2}}{H_{1}^{2}H_{2}^{2}}d\zeta\\ A^{(\bar{\tau}\bar{\eta})}_{nm\ (pq)}=\sum\limits_{k=1}^{l}\int\limits_{\xi_{k}}^{\xi_{k+1}} \overline{B}^{(k)}_{nm}\frac{\partial F_{\eta}}{\partial\xi}\frac{F_{\tau}}{ \partial\xi}\frac{H_{1}H_{2}}{H_{1}^{2}H_{2}^{2}}d\zeta\end{array}  \tag{21}\]  
for \(\tau,\eta=0,1,2,\ldots,N,N+1\), \(n,m=1,2,3,4,5,6\), and \(p,q=0,1,2\). The superscripts \(\eta,\tau\) on the left-hand side of (21) denote the use of the corresponding thickness functions \(F_{\eta},F_{\tau}\), where the tilde denotes that the derivative of the corresponding thickness function with respect to \(\zeta\) has been used. Due to the dependence of all the quantities \(F_{\eta},F_{\tau},H_{1},H_{2}\) in (21) on the thickness coordinate \(\zeta\), the numerical method must be introduced to solve the problem. The notation \(\overline{B}^{(k)}_{nm}\) is related to the elastic constants of the material \(\overline{C}^{(k)}_{nm}\) as follows  
\[