Research on Wave-Added Resistance and Longitudinal Stability Characteristics of Amphibious Aircraft in Rule Wave

Abstract

Assessing the safety of amphibious aircraft hinges significantly on two key factors: wave-added resistance and motion stability during takeoff and landing on water surfaces. To tackle this, we employed the Reynolds-averaged Navier–Stokes (RANS) equations solved via the finite volume method. We utilized the volume fraction method to accurately capture the free surface and employed the overset grid technique to manage the relative motion between the aircraft and the liquid surface. Our approach involves establishing a numerical simulation scheme to investigate the water-planing motion of amphibious aircraft across varying wave heights, wavelengths, speeds, and center-of-gravity positions. The computational findings demonstrate a close match between calculated forces and aircraft motion compared to experimental values. Notably, we observed pronounced nonlinearity in wave-added resistance. Under high sea conditions, operating in a short-wavelength environment or with a rearward center-of-gravity position proves advantageous for reducing wave-added resistance. Conversely, poor longitudinal stability is evident during planing in long waves.

1. Introduction

Amphibious aircraft [1] is a special type of aircraft that can take off in both land and water environments, offering advantages that conventional aircraft cannot match. The water-landing process of an amphibious aircraft consists of two stages: water contact and planing. Due to the use of a hull with a curved or variable incidence angle section on the bottom of the aircraft fuselage and the presence of floats below the wings, wave-added resistance and the nonlinear motion response induced by waves at high speeds are highly sensitive to water surface environmental parameters. Hence, wave-added resistance and the nonlinear motion response induced by waves are highly sensitive to surface environmental parameters at high speeds, thus significantly impacting stability and safety during takeoff and landing operations. Many scholars have conducted research on amphibious aircraft using numerical simulation methods. Qu Qiulin et al. [2,3] simulated the forced landing process of the ARJ21 aircraft and the NACA TN 2929 model on calm water using the FLUENT software with the FVM + VOF method. Guo Baodong et al. [4,5] studied the water dynamic performance of various aircraft rigid models with different aerodynamic layouts based on the Fluent software. Sun Peicheng [6] used the general coupling algorithm based on MSC. Dytran to simulate the water impact during water contact of amphibious aircraft and compared the results with model test data, demonstrating the feasibility of the algorithm. Cai Yufeng [7] studied the water-planing motion of amphibious aircraft on a calm water surface using the STAR-CCM+ software with the overset grid method and the volume fraction method, verifying the effectiveness of the approach. Duan Xupeng [8] studied the unsteady water dynamics and flow phenomena during the takeoff process of amphibious aircraft using OpenFOAM. The research showed that the yaw and pitch motions of the aircraft are mainly determined by hydrodynamics, while aerodynamics play a supporting role; Huang [9,10,11] conducted experiments on the hydrodynamic models of amphibious aircraft and studied their motion response characteristics, further advancing the research on the hydrodynamic performance of amphibious aircraft. Li Xinying [12,13,14] addressed the challenges of large motion response and strong nonlinear flow field in the high-speed water planing of amphibious aircraft. The study proposed a numerical simulation method based on the traditional dynamic mesh technique called “state estimation–accurate computation” and verified the feasibility of the numerical simulation method. Chen Mo [15] developed a numerical simulation scheme using the RANS method to study the motion-response characteristics and typical water-surface skipping behavior of waterborne aircraft during water planing in waves. Hu Kaiye [16] studied the motion response of amphibious aircraft models in regular waves using numerical simulation techniques. The study analyzed the planing resistance, bottom pressure, and free surface wave profiles of the amphibious aircraft under different sea conditions. Chinvorarat Sinchai [17] conducted performance analysis of a light amphibious aircraft using MATLAB. Shi Linfei [18] analyzed the effects of wave and speed variations on the motion response and stability of amphibious aircraft during planing and takeoff using CFD methods and overlapping grid techniques. Zhou Hui [19] studied the effects of different wave factors and speeds on the motion response of amphibious aircraft using the CGFDM method.

Research on amphibious aircraft primarily concentrates on validating the feasibility and accuracy of numerical simulations, alongside studying motion-response characteristics during takeoff and landing. However, scant attention has been given to investigating wave-added resistance and longitudinal stability during high-speed planing in waves. Unlike planing on calm water, taking off in waves can lead to increased resistance and motion response, potentially causing deceleration and instability, thus significantly impacting the well-being of onboard personnel. This poses a critical limitation on operational capabilities and can result in discomfort or, in severe cases, failed takeoff. Understanding the patterns of increased resistance and motion response in waves is essential to determine safe takeoff conditions. To address this, our study employed a scheme that combines overlapping grids with dynamic meshes, integrating the finite volume method (FVM) with the volume of fluid method (VOF). This setup allowed us to simulate the water planing of amphibious aircraft in head waves, enabling us to investigate wave-added resistance and longitudinal stability characteristics under various conditions, including planing speeds, center-of-gravity positions, wave heights, and wavelengths. Ultimately, our research aims to provide valuable insights for enhancing the safety assessment of takeoff and landing operations of amphibious aircraft.

2. Numerical Methods

The STAR-CCM+ version is 15.04. During computation, we applied the finite volume method to discretize the momentum equation while employing a segregated solver to handle the equations in the time domain. To accurately capture the free surface, we utilized the volume of fluid (VOF) method managed through the HRIC scheme. Our computational environment incorporates the effects of gravity, with standard atmospheric pressure serving as the reference pressure for initialization. For turbulence modeling, we adopted the SST k-ω turbulence model, supplemented with the All y+ Wall Treatment wall function for boundary layer treatment. Discretization of convective terms was carried out using a second-order upwind scheme.

Although actual sea conditions typically involve irregular waves, we initially simplified our study by assuming regular wave conditions to establish reference conclusions. To simulate the high-speed planing motion of amphibious aircraft in waves, we employed the boundary wave-making method to simulate first-order regular waves. At the inlet boundary, flow velocity and wave parameters were specified, with only pitch and heave motions of the aircraft being considered during calculation. Prior to computation, the model’s attitude was adjusted to correspond with the initial inclination and draught reflecting the static floating state. The free surface position was set according to a fixed coordinate system to ensure the aircraft remained in static floating equilibrium under the initial conditions. Throughout the calculation process, the flow solver and six-degree-of-freedom motion solver provided data on the aircraft’s force, moment, attitude, and displacement.

2.1. Control Equations

In this study, the computational flow field is considered to be unsteady and incompressible. The numerical model was established based on the fundamental equations of fluid mechanics, specifically the Navier–Stokes equations. The control equations for the three-dimensional, incompressible, unsteady Reynolds-averaged Navier–Stokes (RANS) method are as follows [20]:

$$\frac{\partial {\overline{u}}_i}{\partial x_i}=0$$

(1)

$$\rho \frac{\partial {\overline{u}}_i}{\partial t}+\rho {\overline{u}}_j\frac{\partial {\overline{u}}_i}{\partial x_j}=\rho {\overline{F}}_i-\frac{\partial \overline{p}}{\partial x_i}+\frac{\partial }{\partial x_j}(\mu \frac{\partial {\overline{u}}_i}{\partial x_j}-\rho \overline{{u_i}^{\prime }{u_j}^{\prime }})$$

(2)

In the equation, $\overline{u_i}$ represents the time-averaged velocity, ${u_i}^{\prime }$ represents the turbulent velocity fluctuations, and $-\rho \overline{{u_i}^{\prime }{u_j}^{\prime }}$ represents the Reynolds stress.

2.2. Free Surface Treatment Methods

Accurately capturing and tracking the free surface is crucial for simulating the motion of amphibious aircraft in waves and studying wave-added drag. The volume of fluid (VOF) method was thus employed to capture the free surface between water and air. Its fundamental principle involves constructing and tracking the interface between the two immiscible fluids by using the volume fraction function within each control volume at any given time. This method is particularly suitable for strong nonlinear problems like amphibious aircraft, where splashing phenomena occur.

In the VOF model [21], the distribution and position of the gas–liquid interface are represented by the phase volume fraction α_i. It is defined as follows:

$${\alpha }_i=\frac{V_i}{V}$$

(3)

In the equation, V_i represents the volume occupied by phase i within a grid cell, and V represents the volume of the grid cell itself. In the VOF model, the sum of the volume fractions of all phases within a grid cell is constrained to be equal to 1.

The presence of different phases within a grid cell can be determined based on the values of the phase volume fractions in the grid cell. Taking the example of the current study, where the fluid domain consists of a water phase and a gas phase, if α_i = 0, it indicates that the grid cell contains only the gas phase; if α_i = 1, it indicates that the grid cell contains only the water phase; if 0 < α_i < 1, it indicates that the grid cell represents a region of the free surface where both the water and gas phases are present.

2.3. Turbulence Models

The SST (Menter) k-ω turbulence model [22] considers the transport characteristics of turbulent shear stress and is capable of accurately predicting flow separation points and separation regions caused by adverse pressure gradients. Therefore, it has a significant advantage in predicting complex turbulent flows with separation. The equations for the SST k-ω turbulence model are as follows:

$$\left\{\begin{array}{l}\frac{d}{dt}{\displaystyle \underset{V}{\int }\rho kdV+{\displaystyle \underset{A}{\int }\rho k\left(v-v_g\right)da={\displaystyle \underset{A}{\int }\left(\mu +{\sigma }_k{\mu }_t\right)\nabla k\cdot da+{\displaystyle \underset{V}{\int }\left[{\gamma }_{eff}{G}_k-{\gamma }^{\prime }\rho {\beta }^{*}{f}_{{\beta }^{*}}\left(\omega k-{\omega }_0{k}_0\right)+S_k\right]}dV}}}\\ \frac{d}{dt}{\displaystyle \underset{V}{\int }\rho \omega dV+{\displaystyle \underset{A}{\int }\rho \omega \left(v-v_g\right)da={\displaystyle \underset{A}{\int }\left(\mu +{\sigma }_{\omega }{\mu }_t\right)\nabla \omega \cdot da+{\displaystyle \underset{V}{\int }\left[G_k-\rho \beta f_{\beta }\left({\omega }^2-{\omega }_0^2\right)+S_{\omega }\right]dV}}}}\end{array}\right.$$

(4)

In the equation, V represents the volume, μ is the dynamic viscosity, and σ_k and σ_ω are model coefficients. f_β* and f_β are the free shear correction factor and the turbulent elongation correction factor, respectively. S_k and S_ω are the source terms. k represents turbulent kinetic energy, and ω represents specific turbulent dissipation rate.

3. Numerical Simulation Schemes

3.1. Geometric Models

The present study focused on a specific amphibious aircraft with a hull design featuring a deep V-shaped “long aft body + twisted forebody” structure. The midsection of the hull incorporates a step structure, while the bow of the hull includes wave-dissipating grooves designed to enhance the aircraft’s sea keeping and protect the engines from fluid splashing. The model has a high length-to-width ratio and comprises a single wing structure, T-shaped tail structure, dual-side float structures, and four sets of propeller engines. The model, as shown in Figure 1, includes components such as floats, flaps, ailerons, horizontal stabilizers, and elevators. To optimize computational efficiency, simplifications were made for components such as propellers and float supports. The key parameters of the aircraft model are presented in

Table 1. Key parameters of the scaled-down model aircraft.

Fuselage Length/m	Wingspan/m	Wing Area/m2	Weight/kg	Planing Step Height/m
4.278	4.56	1.75	87.12	0.23
Deadrise angle/°	Beam/m	Draft/m	Flap deflection/°	Rudder deflection/°
22	0.38	0.185	45	0

, with a scale ratio of 1:8.5.

3.2. Computational Domain and Grid Partition

The computational domain and model positioning are shown in Figure 2. Due to the symmetry of the model, only pitch and heave motions were considered, so a half model was used for the simulation. The grid was divided using a structured grid to ensure both computational accuracy and efficiency. For the simulation of the model’s pitch and heave motions, an overlapping grid approach was employed. This involved setting up overlapping grid regions around the aircraft model and the surrounding area, allowing the aircraft to pitch and heave within that region. The basic principle of the overlapping grid technique is to divide the computational domain into multiple subdomains and use domain-sharing methods for information exchange between different subdomains. The advantage of this approach is that information transfer is carried out through domain sharing rather than boundary sharing, which reduces the difficulty of generating subdomain grids. Compared to dynamically deforming grids, overlapping grids provide better grid quality and do not require grid deformation due to aircraft motion. Regardless of the extent of aircraft motion, the grid quality around the aircraft hull is effectively maintained without significant grid deformation.

In this scheme, the velocity inlet boundary of the computational domain was positioned at a distance of 2L from the aircraft’s nose, while the pressure outlet boundary was located at a distance of 3L from the aircraft’s tail. The side velocity inlet boundary was located 2L away from the symmetry plane, the top velocity inlet boundary was situated 1L away from the aircraft keel, and the bottom velocity inlet boundary was positioned 2L away from the keel. The background grid region had a length of 2L, a width of 0.75L, and a height of 1.2L. The overlapping domain had a length of 1.2L, a width of 0.5L, and a height of 0.5L. The calculation scheme adopted two sets of coordinate systems. Under the fixed coordinate system O-XYZ, O is the leading edge position of the nose, the OX axis points to the initial sliding direction, the OY axis points to the right side, and the OZ axis is perpendicular to the right hand rule. In the random coordinate system o-xyz, o is the center of gravity of the aircraft, where the ox axis is pointing from the tail to the nose, the oy axis is pointing from the wing root to the wing tip, and the oz axis is vertically upward according to the right-hand rule.

4. Calculation Results and Discussion

4.1. Verification of Grid and Time Step Independence

The high-speed water planing of amphibious aircraft on the water surface belongs to a strongly nonlinear gas–water two-phase flow problem. Especially during high-speed planing in waves, phenomena such as wave generation, splashing, significant variations in hull pressure gradients, flow separation at the step, wake formation, and Kelvin wake field occur. To accurately simulate the physical phenomena around the aircraft’s hull during water planing and achieve accurate solutions for aerodynamics and hydrodynamics, conducting spatial and temporal sensitivity analysis is the most direct approach to reducing numerical simulation discretization errors. It ensures that the true values of the simulation results are independent of spatial and temporal resolutions.

In terms of grid resolution, it is necessary to refine the grids around the region of free surface surrounding the aircraft as well as the hull’s sliding surface, wings, horizontal stabilizers, and other geometric surfaces experiencing forces. Based on the characteristics of the existing grid scheme, three sets of grid refinement schemes, namely coarse, medium, and fine, were employed. Isotropic refinement was utilized for volumetric grids. It can be observed that the grid quantity increased exponentially with different grid refinement schemes, and the grid sizes are presented in

Table 2. Different grid refinement scheme sizes.

Grid Refinement Schemes	Grid Refinement Schemes for the Region of Free Surface	Grid Refinement Schemes for the Force Components	Grid Size
Coarse grid	15‰L	7.5‰L	6.53 million
Medium grid	7.5‰L	3.75‰L	8.25 million
Fine grid	3.75‰L	1.875‰L	14.36 million

. Regarding the time resolution, six different schemes were used for comparison: 0.01 s, 0.007 s, 0.004 s, 0.003 s, 0.002 s, and 0.001 s.

The grid refinement schematic for the amphibious aircraft is shown in Figure 3 and Figure 4.

For the high-speed water-planing condition of the amphibious aircraft model, with the longitudinal position of the center of gravity at 26% MAC (mean aerodynamic chord) and a velocity of 10 m/s, a comparative calculation was conducted. The model’s static floating state corresponds to a longitudinal pitch angle of 2.94°. Prior to the calculation, the model’s attitude was adjusted to match the pitch angle and draft corresponding to the static floating state.

The results of the calculations for different grid refinement schemes are presented in

Table 3. The comparison of the calculation results for different grid refinement schemes.

	Test Values	Coarse Grid	Deviation (%)	Medium Grid	Deviation (%)	Fine Grid	Deviation (%)
Resistance (N)	195.910	171.534	12.440	185.822	5.150	189.200	3.430
Pitch angle (°)	8.700	7.911	9.070	8.273	4.910	8.430	3.100
Heave (m)	140.000	145.689	−4.060	142.304	−1.650	141.110	−0.790

. During the calm water-planing process of the amphibious aircraft, the accuracy of drag, pitch angle, and lift/drag calculations improved with the increase in grid density in critical areas such as the free surface and force components. In comparison with experimental results, the fine grid scheme yielded a drag deviation of 3.430% compared to the experimental value, a pitch angle deviation of 3.100%, and a lift/drag deviation of 0.790%. This indicates that the grid refinement in critical areas such as the free surface and force components effectively enhances the simulation accuracy of physical quantities like velocity and pressure gradients around the aircraft, enabling accurate solution of the six-degree-of-freedom motion equilibrium equations for the amphibious aircraft under dynamic balance conditions. However, in terms of computational efficiency, the fine-grid scheme takes three times longer than the medium-grid scheme and coarse-grid scheme due to the significant increase in grid size. The difference in computational time between the medium-grid scheme and coarse-grid scheme is relatively small. Moreover, the deviations in drag, pitch angle, and lift/drag calculations for the medium-grid scheme are all within 10%. Therefore, the medium-grid scheme satisfied the research requirements in terms of both accuracy and efficiency.

The calculation results under different time steps are shown in

Table 4. The comparison of results obtained with different time step sizes.

	Test Values	t = 0.001 s	Deviation (%)	t = 0.002 s	Deviation (%)	t = 0.003 s	Deviation (%)
Resistance (N)	195.91	190.822	2.60	189.078	3.49	188.19	3.94
Pitch angle (°)	8.7	8.573	1.46	8.564	1.56	8.396	3.49
Heave (mm)	140	140.704	−0.50	141.358	−0.97	141.634	1.17
	Test values	t = 0.004 s	Deviation (%)	t = 0.007 s	Deviation (%)	t = 0.01 s	Deviation( %)
Resistance (N)	195.91	185.822	5.15	177.334	9.48	162.9	16.85
Pitch angle (°)	8.7	8.273	4.91	7.95	8.62	7.446	14.41
Heave (mm)	140	142.304	1.65	143.91	2.79	144.85	3.46

. The time step also has a significant effect on the calculation accuracy of the force and attitude of the aircraft. Reducing the time step is beneficial to improve the calculation accuracy of drag, pitch, and heave. When t ≤ 0.004 s, although further reducing the time step can improve the numerical simulation accuracy to a certain extent, the calculation efficiency will be greatly reduced. Therefore, t = 0.004 s was used in subsequent calculations.

4.2. Wave-Added Resistance Characteristics during Water Planing

The typical wave heights (h) were chosen as 0.035 m, 0.106 m, and 0.224 m, respectively. Additionally, the influence of different center-of-gravity (CG) positions on the results was investigated. The longitudinal positions of the three CGs were set at 23.5MAC (G1), 26% MAC(G2), and 28.5% MAC(G3).

The formula for calculating wave-added resistance is as follows:

$$R_{aw}=\overline{R_t}-R_{cw}$$

(5)

In the equation, R_αw is wave-added resistance; $\overline{R_t}$ presents average value of total resistance during a stable period in wave navigation conditions; R_cw is stable value of total resistance in calm water conditions at the same speed.

To facilitate comparison, we nondimensionalize R_αw as C_αw:

$$C_{aw}=\frac{R_{aw}}{\rho \mathrm{g}{A}^2(B^2/L)}$$

(6)

In the equation, ρ is the water density, g is the acceleration of gravity, A is the amplitude, B is the width of the aircraft fuselage, and L is the length of the aircraft fuselage.

$$F{r}_{\nabla }=\frac{v}{\sqrt{g{\nabla }^{1/3}}}$$

(7)

In the equation, $v$ is the speed of aircraft, and $\nabla$ is displacement volume.

4.2.1. The Influence of Different Wave Heights on Wave-Added Resistance Characteristics

Through numerical simulations, the time–history curves of resistance, heave, and pitch angle under different wave heights for a center of gravity at G2 and a planing speed of 8 m/s are shown in Figure 5. After data processing, the variation curves of wave-added resistance coefficient, heave amplitude, and pitch angle amplitude with wavelength were obtained and are depicted in Figure 6. As the wavelength increases, the amplitudes of heave and pitch angle response in high-speed planing conditions of the aircraft also increase. However, after λ/L > 1.87, the increase in heave response slows down, and the pitch angle slightly decreases. Nevertheless, under different wave heights, higher wave heights lead to more intense motion responses, indicating a higher risk for the aircraft. The wave-added resistance characteristics of the amphibious aircraft during high-speed planing vary significantly with wavelength, showing distinct differences from conventional displacement-type ships. When λ/L ≤ 1.4, the wave-added resistance coefficient decreases with increasing wavelength. The wave-added resistance coefficient is significantly higher for the 0.035 m compared to 0.106 m and 0.224 m, while the difference between the 0.106 m and 0.224 m is relatively small. When λ/L > 1.4, both 0.035 m (wave height) and 0.106 m (wave height) exhibit further reduction in the wave-added resistance coefficient with increasing wavelength, showing a “resistance reduction” characteristic. The variation of the wave-added resistance coefficient for 0.035 m (wave height) does not follow a clear pattern, while the change in the wave-added resistance coefficient for the 0.106 m becomes less pronounced with increasing wavelength. The 0.224 m (wave height) still exhibits an “increased resistance” characteristic, with a minimum value observed at λ/L = 1.64. This might be due to the fact that amphibious aircraft, compared to conventional ships, have higher speeds and lower draft as well as the influence of lift generated by components such as wings and tail fins. During high-speed planing in waves, parts of the aircraft’s body may lift off the water surface (Figure 7 and Figure 8). At this time, the aircraft is only subject to aerodynamic drag and splashing resistance. Additionally, the complex configuration of the aircraft, including components such as floaters, further contributes to the lack of clear regularities in the variation of wave-added resistance.

4.2.2. The Influence of Center-of-Gravity Position and Planing Speed on Wave-Added Resistance Characteristics

In Figure 9, Figure 10 and Figure 11, the wave-added resistance, heave, and pitch angle amplitudes of the amphibious aircraft under different speeds and wavelengths are shown in 0.224 m (wave height), with the center-of-gravity positions at G1, G2, and G3. In high wave height, different speeds exhibit different motion-response characteristics and wave-added resistance characteristics. When v = 4 m/s (Fr_∇ < 3), the amphibious aircraft is in a transitional sailing state. When λ/L ≤ 1.4, both the heave and pitch angle amplitudes increase rapidly with increasing wavelength, and the differences among different center-of-gravity positions are not significant. When λ/L > 1.4, the increase in heave amplitude becomes less pronounced with increasing wavelength, and shifting the center of gravity aft reduces the magnitude of heave increase. The pitch angle amplitude initially increases and then decreases with increasing wavelength, and a more aft center of gravity results in a smaller pitch angle for the same wavelength. Under different center-of-gravity positions, the wave-added resistance exhibits different trends with increasing wavelength. When λ/L ≤ 0.94, the wave-added resistance decreases rapidly with increasing wavelength for the G1 and G2 configurations. The G3 configuration, with a more aft center of gravity, shows a significant reduction in wave-added resistance. When λ/L > 0.94, the wave-added resistance tends to decrease with increasing wavelength.

When v = 8 m/s (Fr_∇ > 3), the amphibious aircraft is in the planing state, and the magnitudes of heave and pitch angle variations increase significantly compared to v = 4 m/s, resulting in more pronounced motion response. The heave amplitude of the aircraft monotonically increases with increasing wavelength, and for the same wavelength, a more aft center of gravity results in a larger heave amplitude. The pitch angle variations for the G1 and G2 configurations follow similar patterns with wavelength, and their amplitudes are also comparable. When λ/L ≤ 1.87, the pitch angle increases with increasing wavelength; when λ/L > 1.87, the pitch angle decreases with increasing wavelength. The aft center of gravity has a significant influence on the pitch angle and exhibits different variation patterns. When λ/L ≤ 1.17, the pitch angle decreases with increasing wavelength; when λ/L > 1.17, the pitch angle monotonically increases with increasing wavelength. Under different center-of-gravity positions, the wave-added resistance exhibits maximum values. When λ/L ≤ 1.87, the wave-added resistance coefficient increases with increasing wavelength for the G1 configuration, while it decreases first and then increases for the G2 configuration. When λ/L > 1.87, the wave-added resistance coefficient decreases with increasing wavelength for the G1 configuration, while it decreases first and then increases for the G2 configuration. The G3 configuration has the maximum wave-added resistance coefficient at λ/L = 2.34. When λ/L ≤ 2.34, the wave-added resistance coefficient decreases first, then increases monotonically with increasing wavelength. When λ/L > 2.34, the wave-added resistance coefficient decreases with increasing wavelength.

4.3. Study of Longitudinal Stability during Water Surface Planing in Waves

When the amphibious aircraft is planing at high speed on the calm water surface, the dynamic load of the fuselage is too concentrated on the tail planing surface, resulting in longitudinal instability, i.e., “dolphin movement”. This longitudinal instability is related to the velocity (v), the displacement ($\Delta$), and the longitudinal position of the center of gravity (X_g). Compared with the calm water surface, the wave water surface disturbance will cause the aircraft inclination angle and draft to change dramatically, which may cause the dynamic load to concentrate to the tail and aggravate the occurrence of longitudinal instability. Longitudinal instability is a state that loses balance in the longitudinal direction. It includes the state of gliding in calm water (dolphin movement). In waves, the aircraft will jump under the influence of waves. The aircraft may still be able to maintain balance after jumping. It may also capsize. This is not only related to geometric or dynamic characteristics but also to wave conditions.

$$C_B=\frac{\Delta }{\frac{1}{2}\rho v^2{B}^2}$$

(8)

In the equation, C_B is the dynamic load factor, $\Delta$ is the displacement, and B is the beam width.

Numerical simulations were conducted for water surface planing of amphibious aircraft under two different wave conditions, with a significant wave height (h) of 0.035 m. The simulations considered varying center-of-gravity positions and planing speeds. The wavelength-to-ship length ratio (λ/L) was set at 0.701 and 1.403 for the respective wave conditions. The simulations yielded average resistance values, heave amplitudes, and pitch angle amplitudes for different planing speeds and center-of-gravity positions. Based on the results of the pitch angle amplitudes, a longitudinal motion stability boundary curve was plotted.

The time–history curves of the pitch angle for different center-of-gravity positions under longitudinal stable and longitudinal unstable conditions, corresponding to a significant wave height (h) of 0.035 m and wavelength-to-ship length ratios (λ/L) of 0.701 and 1.403, are shown in Figure 12 and Figure 13, respectively. At a planing speed of 15 m/s for longitudinal stability and 16 m/s for longitudinal instability, the sliding states corresponding to the moments of maximum and minimum pitch angles are depicted in Figure 14 and Figure 15, respectively. In the unstable state, the aircraft exhibits a severe motion response with large pitch angle and heave amplitudes, potentially causing the aircraft to lift off the water surface.

The longitudinal motion stability boundary curves for water surface planing under two wave conditions, with a significant wave height (h) of 0.035 m and wavelength-to-ship length ratios (λ/L) of 0.701 and 1.403, were plotted based on the longitudinal stability critical curve equation. The calculated results for the dimensionless longitudinal stable speed limits at different center-of-gravity positions were processed and are presented in

Table 5. Dimensionless calculation results of the upper and lower limits of water surface planing speed combinations.

	λ/L = 0.701			λ/L = 1.403
xg/m	1.908	2.015	2.122	1.908	2.015	2.122
v1/(m·s−1)	9	15	15	9	8	14
v2/(m·s−1)	10	16	16	10	9	15
CB1/(xg/B)	0.003494	0.001191	0.001131	0.003494	0.004187	0.001298
CB2/(xg/B)	0.002830	0.001046	0.000994	0.002830	0.003309	0.001131

.

The longitudinal stability boundary line for water surface planing of amphibious aircraft is generally in the form of an exponential shape and can be represented by the equation shown in Equation (9) [23].

$${\left[\frac{C_B}{x_g/B}\right]}_{\mathrm{critical}}=\mathrm{m}F{r_{\nabla }}^n$$

(9)

The exponential fitting function was transformed into a linear polynomial fitting function for solving, and the longitudinal stability boundary curve equation for water surface planing with a wavelength-to-ship length ratio (λ/L) of 0.701 was obtained as given:

$${\left[\frac{C_B}{x_g/B}\right]}_{\mathrm{critical}}=\frac{0.086}{F{r_{\nabla }}^{2.1591}}$$

(10)

The equation for the longitudinal stability boundary curve of water surface planing with a wavelength-to-ship length ratio (λ/L) of 1.403 is:

$${\left[\frac{C_B}{x_g/B}\right]}_{\mathrm{critical}}=\frac{0.1536}{F{r_{\nabla }}^{2.5937}}$$

(11)

Based on the longitudinal stability critical curve equation, the longitudinal stability boundary curves for wavelength-to-ship length ratios (λ/L) of 0.701 and 1.403 were plotted as shown in Figure 16. The region above the boundary curves represents the unstable area, while the region below the boundary curves represents the stable area. It can be observed that the longitudinal stability boundary curve corresponding to a wavelength-to-ship length ratio of 0.701 is generally located above the curve, corresponding to a ratio of 1.403, indicating a larger stable region. This suggests that the longitudinal stability of the amphibious aircraft during high-speed water surface planing in waves with a wavelength-to-ship length ratio of 0.701 is better than that during high-speed water surface planing in waves with a ratio of 1.403.

5. Conclusions

This paper employs a numerical simulation method to replicate the wave and water surface navigation motion of amphibious aircraft equipped with a full appendage structure. We delved into analyzing and summarizing the impact of wave height, center-of-gravity position, and planing speed on the wave-added resistance of amphibious aircraft. Additionally, we herein propose a longitudinal stability critical curve for high-speed planing in waves under typical working conditions. This critical curve serves as a crucial reference for gauging power demand and conducting safety assessments for amphibious aircraft during planing and takeoff operations. Furthermore, this paper draws the following conclusions:

(1)

The finite volume method and the fluid volume fraction method, combined with overlapping grid techniques, were employed to effectively simulate the forces and motions experienced by the amphibious aircraft during high-speed water surface planing. The calculations considered the hydrodynamic and aerodynamic forces as well as their coupling effects, yielding reasonably accurate results;

(2)

During high-speed water surface planing of amphibious aircraft, the phenomenon of water skipping may occur, and wave resistance exhibits strong nonlinear characteristics. For the studied amphibious aircraft in this paper, as the wave height increases (λ/L ≤ 1.4), the wave resistance decreases. Shifting the center of gravity rearward is beneficial for reducing wave resistance, and higher speeds lead to a more pronounced reduction in wave resistance;

(3)

When the amphibious aircraft is planing at high speed in waves, the aircraft will jump under the influence of waves. The aircraft may still be able to maintain balance after jumping or may lose longitudinal stability and capsize. The influence of wavelength is greater than that of the longitudinal position of the center of gravity. The high-speed motion response in long waves is more intense than that in short waves, which is not conducive to maintaining motion stability.