Research on the water-entry attitude of a submersible aircraft

Background The water entry of a submersible aircraft, which is transient, highly coupled, and nonlinear, is complicated. After analyzing the mechanics of this process, the change rate of every variable is considered. A dynamic model is build and employed to study vehicle attitude and overturn phenomenon during water entry. Experiments are carried out and a method to organize experiment data is proposed. The accuracy of the method is confirmed by comparing the results of simulation of dynamic model and experiment under the same condition. Results Based on the analysis of the experiment and simulation, the initial attack angle and angular velocity largely influence the water entry of vehicle. Simulations of water entry with different initial and angular velocities are completed, followed by an analysis, and the motion law of vehicle is obtained. To solve the problem of vehicle stability and control during water entry, an approach is proposed by which the vehicle sails with a zero attack angle after entering water by controlling the initial angular velocity. With the dynamic model and optimization research algorithm, calculation is performed, and the optimal initial angular velocity of water-entry is obtained. Conclusions The outcome of simulations confirms that the effectiveness of the propose approach by which the initial water-entry angular velocity is controlled.


2015)
, changing of cavity shape (He et al. 2012a, b;Li et al. 2010), and free surface (Han et al. 2014;Lu et al. 2008;Zhang et al. 2012). The typical methods adopted were numerical simulation and experiments. In He et al. (2012), the volume-of-fluid method coupled with the dynamic mesh method was used to simulate the movement of the vertical water-entry body. Wang and Shi (2012) built the oblique water-entry numerical model of airborne missile to simulate the initial water-entry hydroballistics. Park and Jung (2003) utilized the numerical analysis method to study the impact force and ricochet behavior of high-speed water-entry bodies. In Hu and Liu (2014a, b), numerical simulation and experiment were respectively used to study the characteristics of flat-bottomed body entering water, and the results indicated the impact forces mainly depend on water-entry velocities. In Liang et al. (2014), the wing root pivot joint's radial load of a submersible airplane, which imitates the locomotion of gannet's Morus plunge diving, was studied by implementing a test device named Mimic-Gannet. The outcomes revealed the wing load characteristic during the plunge of gannet. May (1975) collected valuable information and experiment data on water entry and systematically analyzed variables that may affect water-entry trajectory. Zhang et al. (2011), Gu et al. (2012) concluded that the head style is crucial to cavities and water-entry stabilization by studying the water entry of objects with different head styles.
Numerical simulation and experiment were the main approaches used in studies on water entry. Most studies focused on fixed work conditions and special phenomenon, neglecting the change of attitude and control during water entry. Nguyen et al. Nguyen (2014) developed an implicit algorithm based on a dual-time pseudo-compressibility method to compute water impact forces on bodies. And in Nguyen et al. (2016), a 6 DOF rigid body motion model and a moving Chimera grid scheme were developed for multiphase flow around moving bodies with application to the water entry problem.
With theoretical analysis methods, a submersible vehicle is taken as the research object in this study. Under the condition of slow velocity, omitting change of free surface, water jet, and cavity, the force mechanism of the object is analyzed to establish a dynamic model of water entry. Finally, the motion law, attitude change, and control of the object during water entry are studied using the dynamic model and optimum searching algorithm.

Shape model
In this study, the peaked arch shape model is used; the model is 2 m long, with a vertex angle of 30° and a uniform mass distribution. It is designed as a linear cutting tail. Figure 1 illustrates the physical model of the vehicle. The computational formula of the radius is shown as the Eq. (1).

Force analysis
The body axis coordinate system, which is selected as reference coordinate system, is built with the origin at the center of the mass (Fig. 2). The water-entry angle is θ. Considering the minimal effect of air force, the vehicle only interacts with the weight W, buoyancy B, and fluid force F during water entry. The weight W remains unchanged, whereas the buoyancy B and fluid force F change along with x a , which is the water-entry distance of the vehicle.

Weight W
During the entire water-entry process, the weight W remains unchanged.
In the Eq. (2), R(x) is the radius of the vehicle, ρ is the density and V is the volume. The center-of-gravity position is: (2)

Buoyancy B
The magnitude and location of the buoyancy B changes along with the water-entry distance x a . The buoyancy B is assumed to be unaffected by the variability of the free surface. The buoyancy and location can be calculated as follows: In Eqs. (4, 5), ρ 0 is the density of water, L in is the length of the vehicle in the water during water entry, and the data range is [L − x a , L]. The distance between the center of mass and the center of buoyancy is:

Fluid force F
During water entry, the fluid force F is difficult to analyze. To analyze the fluid force F expediently in this study, it is divided into ideal and viscous fluid for calculation.

➀ Ideal fluid F i
In terms of the theory of momentum and moment of momentum, the vector and moment of ideal force on the vehicle can be obtained as follows: In above equation, F i is the ideal fluid force on the vehicle, M i is the ideal fluid torque on vehicle, Q f and K f are the ideal fluid force momentum and the moment of momentum on the vehicle, respectively, which can be obtained by the added mass (Yan 2005), ω is the rotation angular velocity, and v is the velocity vector.
According to the theory of potential flow, Q f and K f can be expressed as follows: In above equation, 11 , 12 , . . . , 66 is the added mass, which depends only on the shape of the wetted area.
During water entry, the wetted volume of the vehicle and the added mass change along with time (Liao et al. 2012). Therefore, the variation of added mass must be fully considered. In this study, the effect of the changing added mass change and the rate of change are accounted for. Based on Eqs. (7) and (8) and the definition of added mass (Yan 2005), the ideal fluid force equation can be obtained as follows: By only considering the vertical plane, Eq. (10) can be simplified as follows: In which, = Given the symmetrical characteristic of the vehicle, λ 12 = λ 16 = λ 21 = λ 61 = 0 and λ 26 = λ 62 can be obtained. Then, Eq. (11) can be simplified as: where F ix and F iy are the ideal fluid force on the vehicle at X and Y directions, respectively, M iz is the ideal fluid torque of the vehicle at Z direction; v x and v y are the vehicle's velocity component at X and Y directions, respectively, ω z is the angular velocity of the vehicle at Z direction, and λ is a term of the added mass. When the vehicle enters water, the diving volume and added mass increase; therefore, the profile analysis method is adopted to calculate the added mass (Лoгвинooвч 1969).
The added mass λ 11 of slender body is extremely small (Liao et al. 2012), therefore λ 11 is set to zero. At the same time, ˙ 11 equals zero. By differentiating both sides of Eq. (13), the change rate of added mass can be obtained as follows: In Eq. (14), dx a is the change rate of the vehicle's front part that entered the free surface of water, and dx a = v x dt. R is the radius of the vehicle, as shown in Eq. (1).
➁ Viscous fluid force F μ Experiment and theoretical analysis results indicate that the motion drag coefficient of underwater vehicle is closely related to the velocity, Reynold number, and attack angle (Michael and Jerry 1976). The effect of the varying Reynold number can be ignored here because the vehicle's speed is assumed low and does not vary in a large scale.
The computational fluid dynamics (CFD) method can be used to calculate the drag F μx , lift F μy , moment of force M μz in a given speed, and the attack angle when the vehicle is full immersion. In this study, ANSYS CFX 14.0 is used to calculate the drag F μx , lift F μy , and moment of force M μz .
As seen in Fig. 3, the pressure nephogram of the vehicle when the fluid speed is 20 m/s and the attack angle is zero. The calculation step is set to 0.01 ms, and the end time is 0.01 s. That is, the result is derived by calculating 1000 steps. Then, the drag F μx , lift F μy , and moment of force M μz can be obtained.
The viscous fluid dynamic coefficient C d , C l , and m z can be calculated with Eq. (15). (13)

Fig. 3 The pressure nephogram of the vehicle
The fluid dynamic coefficient under different speeds and attack angles can be calculated through the same method. However, a huge workload is required to calculate the coefficient under every speed and attack angle, which is time consuming. The table lookup interpolation method is adopted to reduce the workload. Considering that the fluid dynamic coefficient is minimally influenced by speed, the speeds 1, 5, 10, 20 m/s and the attack angles 0°, 10°, 20°, …, 90°, are chosen. Figure 4 displays the result.
From Fig. 4, it can be seen that, the drag coefficient C d increased with the increase of the attack angle, while the lift coefficient C l and the moment coefficient m z increase at first when α < 50°, then decrease when α > 50°. The changing law of C l and m z resulted from the shape of the vehicle and the turbulent flow theory.
According to the calculation by CFD, a one-to-one correspondence exists between the viscous fluid dynamic coefficient and the changing relationship of the attack angle and speed. Therefore, three two-dimensional interpolation tables of the viscous fluid dynamic coefficient can be established, in which the line represents the speed and the row represents the attack angle. Based on the current speed and angle, the viscous fluid dynamic coefficient is calculated by interpolation. Then, the viscous fluid force during water entry can be calculated. The viscous fluid dynamic coefficient can be calculated by the equation below.
The fluid dynamic coefficient calculated is under the condition that the vehicle is in full immersion. In practice, the vehicle is partly immersed when crossing the water surface. Given that the fluid dynamic coefficient increases as the wetted area increases, the relationship between them can be assumed linear. The viscous fluid force during water entry can be obtained by: In the equation, S 0 is the superficial area of the vehicle, and S is the wetted area of the vehicle. After the vehicle completely immerses into water, the pitch moment caused by viscous fluid influence is:

Dynamic model
According to the force analysis above, a dynamic model of water entry is established below (Xu et al. 2015).
In the equation, J is the moment of inertia of the vehicle, which can be obtained by the moment of inertia theorem of the disc shape and parallel axis formula.
The differential equation set, which can describe the water-entry process of vehicle consists of Eqs. (1)-(20). Based on the given initial value condition, the equation set can be employed with the classic four-order Runge-Kutta iteration algorithm. The initial value condition includes the initial speed (v x0 , v y0 ), initial water-entry angle θ 0 , initial rotational angular speed ω z0 , and initial distance x a0 , which is between the head peak of the vehicle and water surface along the axial direction.

Experiment model
To verify the accuracy and suitability of the dynamic model, a water-entry experiment was performed on the model. Figure 5 depicts the size and physical diagram of experiment model. The model was composed of head, stern, and body. The head and stern were made of aluminum. The body was made of standard steel tube. The head and stern were assembled to the body with screw thread. The waterproofing gasket was used to guarantee the air tightness. The cylindrical quirk in the head and stern reduced the weight of the experiment model. An acceleration sensor was placed in the cylindrical quirk in the back head. The model weighed 65.1 kg, the volume was 0.052 m 2 , the mean density was 1.252 kg/m 3 , and the rotary inertia was 33.56 kg m 2 . In addition, the barycenter is nearer to the head and 1.08 m to the stern.

Experiment process
On-shore and underwater high-speed cameras were used during the experiment. The on-shore high-speed camera filmed the trajectory from the launch to the water entry. The underwater high speed camera filmed the trajectory of the under-water attitude. An acceleration sensor was fixed on the model head to record the acceleration of the entire process. The experiment model was launched from the second floor to the pool by a high-pressure launch set. Different water-entry experiments can be realized by changing the experiment conditions, such as angle of launcher support and pressure.

Experiment data processing
1. Data processing of on-shore high speed camera The on-shore high-speed camera was set to 1000 fps to film the air trajectory of the experiment model. The video was then sampled every 50 fps (Δt = 0.05 s). The head coordinate (X head , Y head ) and stern coordinate (X tail , Y tail ) of the launcher point were obtained. The linear logistic algorithm was adopted to estimate the dynamic state of the experiment model. The detail procedures were as follows: ➀ The center-of-gravity position was 1.08 m from the stern throughout the experiment. The coordinates of the center-of-gravity position and slope angle were derived from Eqs. (21) and (22). Figure 6 depicts the air trajectory of the model. The red dots denote the center-of-gravity position. Fig. 6 Air trajectory of vehicle model ➁ Linear regression algorithm was adopted to calculate the initial state of the air trajectory (X 0 , v x0 , Y 0 , v y0 , θ 0 , ω 0 ). The model considered gravity regardless of air resistance. The center-of-gravity position and slope angle were: ➂ At the water-entry moment T = t in , horizontal velocity can be denoted as V x_in = v x0 , vertical velocity V y_in = v y0 + gt in , and span rate ω in = ω 0 . The coordinate of the center-of-gravity position and slope angle were directly read from the pictures. Then, the dynamic state (X in , v x_in , Y in , v y_in , θ in , ω in ) was obtained.

Data processing of under-water high speed camera
The underwater high-speed camera filmed the water-entry process. The vision field of underwater high-speed camera was narrow, and the time of motion in the vision field was short. The video was sampled every 5 fps, that is, Δt = 0.005 s, to accurately show the change of dynamic state when the vehicle model entered the water. When the launch angle was 30° and the initial pressure was 0.7 Mpa, the underwater high-speed camera shot the picture at 0.04 s after the vehicle model entered the water. The light area showed the bubbles caused by the water entry. The dark area showed the water. The vehicle model, bubbles, and the inverted reflection were easily distinguished in water, as depicted in Fig. 7. Figure 8 depicts the change of the model's attitude over time.
Illustrating the vehicle model in the completely same scene was difficult because the vision of underwater high-speed camera was narrow. To obtain the attitude and centerof-gravity position during water entry, the pictures were processed as follows:  Fig. 9. ➁ The slope angle was acquired through the easily recognizable front profile, as shown in Fig. 10. The first three pictures were neglected because of the limited water-entry time and vague slope angle. ➂ The position and attitude were obtained according to the outline size and slope angle. Then, the center-of-gravity positions were obtained according to the propositional size. Figure 11 depicts the change of attitude. The yellow line represents the outline, and the red line indicates the water surface.

Data processing of acceleration sensor
The acceleration sensor, as depicted in Fig. 12, recorded the change of impact loading from the launch to the complete water-entry process. This study investigated the water

Model verification
To verify the accuracy of the proposed dynamic equation, experiment and simulation were carried out under the same condition, and the results were compared. The water-entry state was not under control because of the restriction of experiment condition. With the pressure of 0.7 Mpa and launch angle of 30°, the initial water-entry state parameters of different launcher angles were obtained using the proposed experiment data processing. The horizontal velocity, vertical velocity, span rate, attack angle, and slope angle at the water-entry moment were derived as 21.7, 19.48 m/s, 25.43°/s, 18.02°, and 59.93°, respectively. Table 1 shows the change of attitude and center-of-gravity position over time in the experiment.
The initial simulation condition was set to the same value as the initial water-entry state parameter in the experiment. The simulation time was set to 0.15 s. Figure 13 displays the results of experiment and simulation. The blue solid box shows the vision field  Fig. 15 shows the comparison of the change of axial direction (X axis) and radial direction (Y axis) loading. In Fig. 15, when the vehicle entered the water, the radial direction loading of the vehicle plummeted to −25 g at around 0.07 s. It indicates that when the vehicle entered the water, the radial direction would withstand a big overload at around 0.07 s. Figures 13,  14 and 15 reveal that the simulation was consistent with the experiment, thus the proposed dynamic equations were reasonable. Nevertheless, some errors were encountered. The main causes were as follows: 1. The simulation neglected the effects of perturbation of the water surface, water-entry cavity, and splashing. 2. The model was treated as rigid body in simulation, which was inconsistent with reality. 3. The viscous hydrodynamic coefficients were calculated by interpolation, which caused the errors. 4. Errors also existed during data acquisition and processing.  In addition, there are some limitations of the dynamic model: (1) the model can be applied for low-velocity problem only; (2) the model can be applied to axisymmetric body or symmetric flow only; (3) usually, the added mass method is only stable when the body mass is larger than the added liquid mass due to the moving body.
From the results of experiment and simulation, overturning can easily happen. The overturn phenomenon is strongly connected with the initial water-entry angle of attack and rotational angular velocity.

Model simulation
Using the proposed dynamical model, other simulations were run to reveal the relationship among the overturn, and angle of attack, and rotational angular velocity.

Influence of initial angle of attack
The water-entry velocity was reset to 20 m/s, water-entry to 60°, and rotational angular to 0. The initial water-entry angle of attack was set to 0°, ±5°, ±10°, and ±20°, considering the structure tolerance of the water-entry impact force (Yang et al. 2001). The simulation results are shown in Figs. 16, 17, 18, and 19.
The results showed that when α < 0°, the increase of the absolute angle of attack value decreased water-entry depth, horizontal displacement, and slope angle. It also accelerated the increase of negative rotational angular velocity, negative angle of attack, and slope angle. The vehicle rotated downwards and even overturned. When α ≥ 0°, the decrease of initial angle of attack value decreased water-entry depth and horizontal displacement and accelerated the decrease of slope angle. It accelerated the increase of positive rotational angular velocity and positive angle of attack. The vehicle rotated downwards and even overturned. The decrease of initial water-entry angle of attack When the vehicle reached the water surface, the fluid torque influenced the motion of vehicle and increased with high initial angle of attack. As the vehicle entered water, the buoyancy and buoyancy center changed. The misalignment between center-of-gravity position and buoyant center resulted to another torque. Both torques were affected simultaneously. When α < 0°, the initial angle of attack and fluid torque were both negative. The vehicle began to rotate downwards. The increase of slope angle caused the negative increase of the angle of attack and rotational angular velocity. Both torques accelerated the change of angle of attack and slope angle when the vehicle entered the water completely. The overturn occurred downwards when the slope angle exceeded the threshold. When α ≥ 0°, the initial angle of attack and fluid torque were both positive. The vehicle began to rotate upwards. The decrease of slope angle caused the positive

Influence of initial rotational angular velocity
The water-entry velocity was reset to 20 m/s, water-entry to 60°, and initial water-entry angle of attack to 0°. The initial water-entry rotational angular velocity was set to 0°/s, ±10°/s, ±30°/s, and ±45°/s. The simulation results are shown in Figs. 20, 21, 22, and 23. The results showed that initial rotational angular velocity affected the water-entry process. The effects of the two torques were the same as those mentioned above. When ω ≥ −10°/s, the increase of initial rotational angular velocity accelerated the decrease of water-entry depth, horizontal displacement, and slope angle. It also accelerated the increase of positive rotational angular velocity and angle of attack. The vehicle rotated upwards and even overturned. When ω ≤ −30°/s, the decrease of initial rotational angular velocity accelerated the decrease of water-entry depth, horizontal displacement, rotational angular velocity, and angle of attack. It also accelerated the increase of slope angle. The vehicle rotated downwards and even overturned. Choosing the proper water-entry initial rotational angular velocity is crucial to prevent the overturn phenomenon.

Water-entry control
Water entry was difficult to control because the process was transient and complex. The initial angle of attack and initial rotational angular velocity influenced the water-entry trajectory. The vehicle was more stable with a smaller angle of attack after entering the water completely. Vehicle attitude should be controlled before water entry to guarantee stability after entering the water completely. Attitude control mechanism could adjust the initial rotational angular velocity to achieve the navigation with zero-attack after entering the water completely. The dynamic model proposed in "Establishment of dynamical model" section is a multivariate differential equation set. Finding an analytical solution is difficult, and, thus, this problem can be solved by converting it to constraint optimization problem.
where function U(ω) is the proposed dynamic model in "Establishment of dynamical model" section, ω z-begin is the controllable initial rotational angular velocity, and α is the angle of attack. When an initial rotational angular velocity ω z-begin was given, we can use the dynamic model proposed in "Establishment of dynamical model" section to calculate the attack angle α. In the constraint condition, T is the total time to complete the waterentry process. The condition can be avoided by deriving a substantial calculated ω z-begin , which may cause a larger π angle variation of the inclination angle before the vehicle completely enters the water.
In Eq. (26), ω z-begin is the input variable of U(ω), which allows |α| to be minimum. It is also a unimodal function according to the real vehicle motion. The golden section algorithm (Chen 2005) was adopted to search for the optimal value.
Step 5 k = k + 1, go to step 2. Table 2 shows the optimal rotational angular velocity when water-entry velocity was 20 m/s, slope angle was 60°, and angle of attack were 0°, ±5°, and ±10°. The results are shown in Figs. 24,25,26,and 27. Simulation results showed that, with the optimal initial rotational angular velocity, the slope angle changed slightly. After entering the water, the angle of attack and rotational angular velocity approached zero quickly. The trajectory was stable and straight; thus, the overturn phenomenon was avoided, which was beneficial for the subsequent underwater control.