## 2 Numerical Simulation Methods
### Rigid-Flexible Coupling Equation and Numerical Solving Method
Figure 1: The flexible aircraft before and after deformation. Considering gravity as distributed force, the potential energy of the system only includes the deformation energy of the structure:  
\[U=\frac{1}{2}\eta^{T}\cdot K_{e}\cdot\eta \tag{8}\]  
In order to derive the nonlinear equations of motion the partial derivatives \(\frac{\partial L}{\partial V_{c}}\), \(\frac{\partial L}{\partial\omega}\), \(\frac{\partial L}{\partial\dot{\theta}}\), \(\frac{\partial L}{\partial\dot{\theta}}\) have to be derived first. These partial derivatives are presented here:  
\[ \begin{array}{rl}\frac{\partial L}{\partial V_{c}}&=\int_{V}\rho_{m}\big{(}V _{c}+\phi\dot{\eta}+\omega\times(p+\phi\eta)\big{)}dV\\ &=\int_{V}\rho_{m}\big{(}V_{c}+\phi\dot{\eta}-\left(\stackrel{{ \sim}}{{\rho}}+\stackrel{{\sim}}{{\phi\eta}}\right)\cdot\omega \big{)}dV\\ &=M_{\rm c}\cdot V_{c}+\underbrace{\int_{V}\rho_{m}\phi dV\cdot\dot{\eta}}_{S_ {1}}\cdot\eta-\underbrace{\int_{V}\rho_{m}\left(\stackrel{{ \sim}}{{\rho}}+\stackrel{{\sim}}{{\phi\eta}}\right)dV}_{X_{1}} \cdot\omega\\ &=M_{\rm c}\cdot V_{c}+\underbrace{S_{1}\cdot\eta+X_{1}^{T}\cdot\omega}_{X_{1} }\\ \end{array}  \tag{9}\]  
\[ \begin{array}{rl}\frac{\partial L}{\partial\omega}&=\int_{V}\rho_{m}\bigg{(} \left(\stackrel{{\sim}}{{\rho}}+\stackrel{{\sim}}{{ \phi\eta}}\right)\cdot V_{c}+\left(\stackrel{{\sim}}{{\phi\eta}} \right)\cdot\phi\dot{\eta}+\left(\stackrel{{\sim}}{{\rho}}+ \stackrel{{\sim}}{{\phi\eta}}\right)^{T}\left(\stackrel{{ \sim}}{{\rho}}+\stackrel{{\sim}}{{\phi\eta}}\right)\cdot \omega\bigg{)}dV\\ &=\int_{V}\rho_{m}\left(\stackrel{{\sim}}{{\rho}}+\stackrel{{ \sim}}{{\phi\eta}}\right)dV\cdot V_{c}+\left(\underbrace{\int_{V}\rho_{m} \widetilde{\rho}\phi dV}_{S_{2}}+\underbrace{\int_{V}\rho_{m}\left(\stackrel{{ \sim}}{{\phi\eta}}\right)\phi dV}_{X_{2}}\right)\cdot\dot{\eta}\\ &+\underbrace{\int_{V}\rho_{m}\left(\stackrel{{ \sim}}{{\rho}}+\stackrel{{\sim}}{{\phi\eta}}\right)^{T}\left( \stackrel{{\sim}}{{\rho}}+\stackrel{{\sim}}{{\phi\eta} }\right)}_{I}dV\cdot\omega\\ &=\mathrm{X}_{1}\cdot V_{c}+(S_{2}+X_{2})\cdot\dot{\eta}+J\cdot\omega\\ &\frac{\partial L}{\partial\dot{\eta}}&=\int_{V}\rho_{m}\big{(}\phi^{T}\cdot V _{c}+\phi^{T}\phi\dot{\eta}\dot{\eta}+\phi^{T}\cdot(\omega\times(\rho+\phi\eta ))\big{)}dV\\ &=\mathrm{S}_{1}^{T}\cdot V_{c}+M_{e}\cdot\dot{\eta}+\left(S_{2}^{T}+X_{2}^{T} \right)\cdot\omega\\ \end{array}  \tag{10}\]  
\[ \begin{array}{rl}\frac{\partial L}{\partial\dot{\eta}}&=\int_{V}\rho_{m}\big{(} \phi^{T}\cdot V_{c}+\phi^{T}\widetilde{\omega}^{T}\phi\cdot\dot{\eta}+\phi^{T }\widetilde{\omega}\widetilde{\rho}\cdot\omega+\phi^{T}\widetilde{\omega}^{T} \widetilde{\omega}\phi\cdot\eta\big{)}dV-K_{e}\cdot\eta\\ &=S_{1}^{T}\cdot\widetilde{\omega}^{T}\cdot V_{c}+H_{1}^{T}\cdot\dot{\eta}+H_{2 }\cdot\omega+H_{3}\cdot\eta-K_{e}\cdot\eta\end{array}  \tag{12}\]  
Calculating the time derivative of the above terms and rearranging them, the equations of motion of a deformable aircraft become the following:  
\[ \begin{array}{rl}\begin{pmatrix}M_{\rm c}&X_{1}^{T}&S_{1}\\ X_{1}&J&S_{2}+X_{2}\\ S_{1}^{T}&S_{2}^{T}+X_{2}^{T}&M_{e}\end{pmatrix}\begin{pmatrix}\dot{V}_{c}\\ \dot{\omega}\\ \dot{\eta}\end{pmatrix}\\ +\begin{pmatrix}\widetilde{\omega}M_{\rm c}&\widetilde{\omega}X_{1}^{T}&2 \widetilde{\omega S_{1}}\\ \widetilde{\omega}X_{1}&-\left(Y_{3}+Y_{3}^{T}\right)+\widetilde{V}_{c}X_{1}^{T }+\widetilde{\omega}J&\widetilde{\omega}(S_{2}+X_{2})\\ -S_{1}^{T}\widetilde{\omega}^{T}&-H_{2}&2H_{1}+C_{e}\end{pmatrix}\begin{pmatrix}V _{c}\\ \omega\\ \dot{\eta}\end{pmatrix}\\ +\begin{pmatrix}0&0&0\\ 0&0&0\\ 0&0&-H_{3}+K_{e}\end{pmatrix}\begin{pmatrix}R_{I}\\ \theta\\ \eta\end{pmatrix}=\begin{pmatrix}F\\ M\\ Q\end{pmatrix}\end{array}  \tag{13}\]The expressions for each coefficient are as follows:  
\[