## 2 Methodology
### Research Methodology  
#### 2.2.1 Simulation Model  
This study employs the large eddy simulation (LES) technique, which is known for its effectiveness in numerically predicting the flow separation around bluff bodies, as validated in previous studies [28; 29]. The standard \(k\)-\(\epsilon\) PANS model has certain limitations in simulating flows with strong vortices and large-curvature bent wall flows. However, the Renormalization Group (RNG) \(k\)-\(\epsilon\) PANS model has shown improvements in predicting such flows [43]. To obtain more accurate results, the RNG \(k\)-\(\epsilon\) PANS turbulence model is adopted in this study, and its governing equations are expressed as follows:  
\[\frac{\partial(\rho k_{u})}{\partial t}+\frac{\partial\big{(}\rho U_{j}k_{u} \big{)}}{\partial x_{j}}=\frac{\partial}{\partial x_{j}}\Bigg{[}\alpha_{k} \bigg{(}\mu+\frac{\mu_{u}}{\sigma_{u}}\bigg{)}\frac{\partial k_{u}}{\partial x _{j}}\Bigg{]}+P_{ku}-\rho\varepsilon_{u} \tag{1}\]  
\[\frac{\partial(\rho\varepsilon_{u})}{\partial t}+\frac{\partial\big{(}\rho U_ {j}\varepsilon_{u}\big{)}}{\partial x_{j}}=\frac{\partial}{\partial x_{j}} \Bigg{[}\alpha_{k}\bigg{(}\mu+\frac{\mu_{u}}{\sigma_{u}}\bigg{)}\frac{ \partial\varepsilon_{u}}{\partial x_{j}}\Bigg{]}+C^{*}_{\epsilon 1}P_{ku}\frac{ \varepsilon_{u}}{k_{u}}-C^{*}_{\epsilon 2}\rho\frac{\varepsilon_{u}^{2}}{k_{u}} \tag{2}\]  
where \(U_{j}\) represents the resolved velocity field, \(t\) denotes time, \(\rho\) represents fluid density, \(\mu\) denotes the dynamic viscosity coefficient, \(\mu_{u}\) represents the turbulent viscosity coefficient, \(f_{k}\) represents the unresolved turbulent kinetic energy ratio, \(f_{\epsilon}\) represents the unresolved turbulent kinetic energy dissipation rate ratio, \(k_{u}\) represents the unresolved local time-averaged turbulent kinetic energy, and \(\varepsilon_{u}\) represents the unresolved local time-averaged turbulent kinetic energy dissipation rate.  
In addition,  
\[f_{k}=\frac{k_{u}}{k} \tag{3}\]  
\[f_{\epsilon}=\frac{\varepsilon_{u}}{\varepsilon} \tag{4}\]  
\[\mu_{u}=\rho C_{\mu}\frac{k_{u}^{2}}{\varepsilon_{u}} \tag{5}\]  
\[\sigma_{u}=\frac{f_{k}^{2}}{f_{\epsilon}} \tag{6}\]  
\[C^{*}_{\epsilon 2}=C^{*}_{\epsilon 1}+\frac{f_{k}}{f_{\epsilon}}(C_{\epsilon 2 }-C^{*}_{\epsilon 1}) \tag{7}\]  
Figure 1: Posture parameters and forces acting on the multi-body system during flight. Posture parameters include flight speed V, attack angle \(\phi\), angle between velocity and skis \(\alpha\), skis’ opening angle \(\lambda\), angle between body and skis \(\theta\), bending angle of upper body \(\beta\). Forces include drag, lift, yaw force, and gravity.  
\[\eta=\left(2S_{ij}\cdot S_{ij}\right)^{1/2}\frac{k_{u}}{\varepsilon_{u}}\frac{f_{ t}}{f_{k}} \tag{8}\]  
\[S_{ij}=\frac{1}{2}\left(\frac{\partial U_{i}}{\partial x_{j}}+\frac{\partial U_ {j}}{\partial x_{i}}\right) \tag{9}\]  
\[C_{c1}^{*}=C_{t1}-\frac{\eta(1-\eta/\eta_{0})}{1+\delta\eta^{3}} \tag{10}\]  
The values of the constants in the model are \(C_{\mu}=0.0845\), \(\alpha_{k}=\alpha_{\varepsilon}=1.39\), \(C_{t1}=1.42\), \(C_{t2}=1.68\), \(\eta_{0}=4.377\), and \(\delta=0.012\).  
The governing equations are discretized using the finite volume method. The coupling of pressure and velocity is solved using a consistent and coordinated approach based on the semi-implicit algorithm for pressure-linked equations (SIMPLEC). The time discretization is performed using a second-order difference scheme. The energy dissipation rate ratio and velocity terms are discretized using a second-order upwind scheme. The time step size is set to 0.0001 s, and the total duration for the integration of flow equations is about 5 s.