## Theoretical Modeling
### Peters' Finite State Inflow Theory for Aerodynamic Modeling  
The unsteady aerodynamic loads used in the current framework are based on the Peters' 2D finite state inflow theory [45]. The model calculates aerodynamic loads on a thin airfoil section undergoing large motions in an incompressible flow. The aerodynamic loads per unit span calculated about the aerodynamic center (A.C.) are given as:  
\[l_{ac}=\pi\rho b^{2}(\dot{x}+\dot{y}\dot{\alpha}-d\ddot{\alpha})+c_{la}\rho b \dot{y}^{2}\left[\frac{\dot{z}}{\dot{y}}+(\frac{1}{2}b-d)\frac{\dot{\alpha}}{ \dot{y}}-\frac{\lambda_{0}}{\dot{y}}\right]+\rho b\dot{y}^{2}(c_{l0}+c_{l\delta}\delta) \tag{38}\]  
\[m_{ac}=-\pi\rho b^{3}\left[\frac{1}{2}\ddot{z}+\dot{y}\dot{\alpha}+(\frac{1}{ 8}b-\frac{1}{2}d)\ddot{\alpha}\right]+2\rho b^{2}\dot{y}^{2}(c_{m0}+c_{m\delta}\delta) \tag{39}\]  
\[d_{ac}=\rho b\dot{y}^{2}c_{d0} \tag{40}\]  
where \(\rho\) is the air density, \(b\) is the semichord, \(d\) is the distance of the mid-chord in front of the reference axis (here, it is the shear center (C)), and \(\delta\) is the trailing-edge flap deflection angle. \(c_{la}\) is the lift curve slope and \(c_{l\delta}\) and \(c_{m\delta}\) are the lift and moment slopes due to flap deflection, respectively. \(c_{l0}\), \(c_{m0}\) and \(c_{d0}\) are the lift, moment and drag coefficients for zero angle of attack (AOA), respectively. Furthermore, \(\lambda_{0}\) is the inflow parameter, which accounts for induced flow due to free vorticity.  
\[\lambda_{0}=\frac{1}{2}\mathbf{b}^{T}\lambda \tag{41}\]  
where \(\lambda\) is the vector of inflow states [45], which is given by:  
\[\mathbf{A}\dot{\lambda}+\frac{\dot{y}}{b}\lambda=\mathbf{c}[\ddot{z}+\dot{y} \dot{x}+b(\frac{1}{2}-d)\dot{\alpha}] \tag{42}\]  
where \(\mathbf{A}\), \(\mathbf{b}\) and \(\mathbf{c}\) are constants given as:  
\[\mathbf{A}=\mathbf{D}+\mathbf{d}\mathbf{b}^{T}+\mathbf{c}\mathbf{d}^{T}+ \frac{1}{2}\mathbf{c}\mathbf{b}^{T}\quad b_{i}=(-1)^{i-1}\frac{(N+i)!}{(N-i)!} \frac{1}{(i!)^{2}}\quad c_{i}=\frac{2}{i}\]  
\[d_{1}=0\quad d_{i}=\frac{1}{2}\ (i\
eq 1)\quad D_{ij}=\frac{1}{2i}\ (i=j+1)\quad D _{ij}=-\frac{1}{2i}\ (i=j-1);\quad D_{ij}=0\ (i\
eq j\pm 1)\]  
where \(N\) is the number of finite states.  
To transfer the loads from A.C. to the wing reference axis (here, it is shear center (C)), one may use  
\[l_{ra}=l_{ac};\quad m_{ra}=m_{ac}+(\frac{1}{2}b+d)l_{ac};\quad d_{ra}=d_{ac} \tag{43}\]  
Furthermore, the aerodynamic loads are rotated to the body coordinate system, and then the global coordinate system, which yields:  
\[\mathbf{F}_{a}=\mathbf{U}_{e}\begin{bmatrix}\cos(\alpha)&0&-\sin(\alpha)\\ 0&1&0\\ \sin(\alpha)&0&\cos(\alpha)\end{bmatrix}\begin{bmatrix}\cos(\beta)&\sin(\beta )&0\\ -\sin(\beta)&\cos(\beta)&0\\ 0&0&1\end{bmatrix}\begin{bmatrix}-d_{ra}\\ 0\\ -l_{ra}\end{bmatrix}l \tag{44}\]\[\mathbf{M}_{a}=\mathbf{U}_{e}\begin{bmatrix}\cos(\alpha)&0&-\sin(\alpha)\\ 0&1&0\\ \sin(\alpha)&0&\cos(\alpha)\end{bmatrix}\begin{bmatrix}\cos(\beta)&\sin(\beta) &0\\ -\sin(\beta)&\cos(\beta)&0\\ 0&0&1\end{bmatrix}\begin{bmatrix}0\\ m_{ra}\\ 0\end{bmatrix}l \tag{45}\]  
where \(l\) is the local span, which is here the length of the corresponding element.  
Finally, the elemental force and moment are allocated equally to the two nodes of each element.