Modeling and Motion Simulation for A Flying-Wing Underwater Glider with A Symmetrical Airfoil

The flying-wing underwater glider (UG), shaped as a blended wing body, is a new type of underwater vehicle and still requires further research. The shape layout and the configuration of the internal actuators of the flying-wing UG are different from those of “legacy gliders” which have revolving bodies, and these two factors strongly affect the dynamic performance of the vehicle. Considering these differences, we propose a new configuration of the internal actuators for the flying-wing UG and treat the flying-wing UG as a multi-body system when establishing its dynamic model. In this paper, a detailed dynamic model is presented using the Newton-Euler method for the flying-wing UG. Based on the full dynamic model, the effect of the internal actuators on the steady gliding motion of vehicle is studied theoretically, and the relationship between the state parameters of the steady gliding motion and the controlled variables is obtained by solving a set of equilibrium equations. Finally, the behaviors of two classical motion modes of the glider are analyzed based on the simulation. The simulation results demonstrate that the motion performance of the proposed flying-wing UG is satisfactory.


Introduction
Underwater gliders are a class of buoyancy-driven autonomous underwater vehicles (AUVs). They glide through the ocean by changing their net buoyancy with the buoyancy control system and controlling their attitude with internal sliding masses or rudder. The buoyancy control system and internal sliding masses are known as internal actuators. Underwater gliders have been broadly applied in environment monitoring, resource exploration and military field because of long endurance, long distance, low noise and low cost.
In 1989, Henry Stommel (1989) published a revolutionary article about buoyancy driven floats for oceanography. Since then, some outstanding gliders, such as the Slocum glider (Webb et al., 2001), the Spray glider (Sherman et al., 2001), the Seaglider (Eriksen et al., 2001) and PETREL (Wang et al., 2011), have been developed and widely applied in various fields. It is noteworthy that the shape layouts of all these "legacy gliders" are a revolving body with two wings, and the most common revolving body is a cylinder with semi-elliptical head and tail, such as the Slocum glider. Based on this unique shape, these legacy gliders are similar in the configuration of internal actuators.
The motion performance of the underwater glider relies on the external forces exerted on the glider, and the shape of the glider is an important factor in determining the hydrodynamic forces and moments. Because the cylindrical hull hardly generates any lift force in water, the lift-to-drag ratio of the glider with a cylindrical body is low and this is one of the main reasons why the motion performance of these legacy gliders is limited . To achieve higher hydrodynamic efficiency, a flying-wing UG shaped as a blended wing body was proposed. Jenkins et al. (2003) fully studied the feasibility of the flying-wing UG, and found that the maximum lift-to-drag ratio of their model reached 17, which was much higher than that of the legacy gliders. After that, a flying-wing UG named XRay was developed, and the lift-to-drag ratio of the XRay exceeded 10 (D' Spain et al., 2007). Then, a new flying-wing UG named ZRay was designed, and the lift-to-drag ratio of ZRay exceeded 35 (D'Spain, 2009). Using computational fluid dynamics (CFD) and model-scale testing, Haase et al. (2017) did a research on predicting the lift-to-drag ratio and assessing the static stability of a blended wing-body glider in pitch and yaw. Based on their study, a new flying-wing UG named Sun Ray Glider was designed and a series of underwater tests were conducted. These studies have demonstrated that the flying-wing UG has better hydrodynamic efficiency compared with the legacy glider and has received great attention. However, the shape of the flying-wing UG is more complicated than that of the legacy glider and facing great challenges in the design and manufacturing process. The published research on the flying-wing UG is not abundant yet and thus requires further research.
In recent years, the Northwestern Polytechnical University has conducted some researches on the shape design of the flying-wing UG and made some achievements. We defined the longitudinal section of the flying-wing UG with a symmetrical airfoil which was chosen from the National Advisory Committee for Aeronautics (NACA) airfoils. Then this two-dimensional (2D) shape was used consistently throughout the full span of vehicle, gradually tapering to either wingtip. Based on the geometric model of the shape described above, the shape of a flying-wing UG was optimized using the lift-to-drag ratio as the optimization target (Sun et al., 2015). However, the result of this method loses large volume. Then the maximum gliding range was used as the target to optimize the shape of flying-wing UG (Sun et al., 2017). Besides, we actively carried out a series of studies for the flying-wing UG including modeling, analysis, and experiments.
Dynamic modeling and motion simulation always play an important role in the study of underwater vehicles. Based on Euler-Langrange system, a nonlinear six degrees of freedom (DOF) dynamic model of the marine vehicles, which has been widely used for modeling and control of the traditional AUVs, was described in Fossen (1994). Leonard and Graver (2001) modelled a laboratory-scale underwater glider as a multi-body system, composed of a rigid body and several point masses. Then Graver (2005) explained the dynamic model for the underwater glider in details. This model has been widely used for the legacy gliders. Utilizing the Gibbs and Appell equations, the mathematical model for a thermal UG was proposed to take the internal actuators into consideration (Wang and Wang, 2009). Using the empirical approach, Hussain et al. (2011) developed the mathematical model of an underwater glider and presented the identification on net buoyancy, depth and pitching angle of the glider. Wang et al. (2011) treated a hybrid-driven underwater glider as a multi-rigid-body system and developed the dynamic model based on momentum theorem. Zhang et al. (2013) did some research on the steady spiraling motion of the legacy glider, and a recursive algorithm was proposed to solve the complex dynamic equations of the steady spiraling motion. Taking into consideration of the effect of a nonuniform flow on the dynamics of the glider, Fan and Woolsey (2014) presented a nonlinear multi-body dynamic model for an underwater glider operating in an unsteady, nonuniform flow. Considering the hull deformation and seawater density variation, Yang et al. (2017) established a dynamic mod-el for the deep-sea underwater glider.
From these studies on the legacy gliders and AUVs, we can see that dynamic models of different underwater vehicles are not the same. The differences between these dynamic models include (but are not limited to) varied vehicle configuration, hydrodynamic profile and working environment. In addition to the hydrodynamic profile, the flyingwing UG and the legacy gliders are also different in the configuration of their internal actuators. First, the roll control system of the legacy glider is usually a rotatable eccentric mass subject to the constraints of its cylinder body, while the flying-wing UG has enough space to set other forms of the roll control system. For example, in Sun et al. (2017), a new kind of internal actuators layout of the flying-wing UG was proposed, where two tanks of buoyancy-regulating are symmetrically placed to control the buoyancy and roll of the glider. Secondly, the buoyancy control system of the flyingwing UG is much heavier than that of the legacy gliders and would change the position of its center of gravity (CG), so the buoyancy control system cannot be treated as a mass with a fixed location. Therefore, it is essential to establish a detailed dynamic model for the flying-wing UG.
In this paper, considering that the sliding mass is much better than the buoyancy control system for roll control, we propose a new configuration of the internal actuators for the flying-wing UG in Section 2. A detailed dynamic model for the flying-wing UG is built using the Newton-Euler method in Section 3. When establishing the dynamic model, the flying-wing UG is treated as a multi-body system to consider the effect of the internal actuators. Based on this detailed dynamic model, the effect of the internal actuators on the steady gliding motion is studied in Section 4. Further simulation and analysis on two classical motion modes is conducted in Section 5. Note that the shape of the glider we used is based on the shape optimized by Sun et al. (2017). A brief summary is given in Section 6.

Configuration of a flying-wing underwater glider
As shown in Fig. 1a, the shape of flying-wing UG is totally different from that of the legacy gliders with a torpedo shape. This unique shape layout allows for the carriage of larger payloads (Jenkins and D'Spain, 2016). Sun et al. (2017) proposed to use two tanks which were symmetrically placed to control the buoyancy and roll of the flyingwing UG. Considering that the buoyancy can be easily affected by the seawater density variation and the density of seawater varies with temperature, salinity and pressure, which vary with different seas (Yang et al., 2017). The sliding mass is much better than the buoyancy control system for roll control and a new configuration of the internal actuators for the flying-wing UG as below is proposed.
The flying-wing UG we study in this paper measures 8.5 m in wing span and has a volume of 3.53 m 3 . As shown in Fig. 1b, the vehicle is driven by a buoyancy control system m rb and is steered by two sliding masses. The buoyancy control system could change the net buoyancy of the glider by pumping water in or out of the reservoir, and the sliding masses are constrained to move along the suspension system. The redistribution of these internal actuators will affect the attitude of the glider. In light of the configuration of underwater glider, many researchers model the underwater glider as a multi-body system. Based on the same assumption, we treat the flying-wing UG as a multi-body system that consists of a rigid body (mass ) and three moving point masses. Here, we denote the buoyancy control system as a point mass . Note that, the point mass has a variable position and we will discuss it in detail later. The moving mass which slides along the vehicle centerline for pitch control is denoted as the pitch control mass . Different from legacy gliders, the roll control system of the flying-wing UG here is composed of a movable mass which moves along the wing of the glider and a fixed mass which is symmetrically placed to keep balance ( ). For the convenience of analyzing, the roll control system is treated as a moving particle and denoted as the roll control mass in this work. Then the total vehicle mass is Considering that the vehicle is operated in shallow water, the change of seawater density can be ignored for buoyancy control. The glider displaces a fixed volume of fluid . Thus, the net mass of the vehicle is 3 Dynamic modeling for a flying-wing underwater glider

Frames of reference and kinematics
{i, j, k} To describe the motion of the multi-body system, three coordinate frames, i.e., the inertial frame, the body frame and the flow frame, are defined as shown in Fig. 2 and Fig. 3. Definition of the three frames is the same as that of the common underwater vehicles (Fossen, 1994;Graver, 2005;Zhang et al., 2013). The inertial frame, represented by an orthonormal triad , is a reference frame which is fixed m rb in the inertial space and its origin is fixed on the sea surface. The body frame is fixed with the rigid body (mass ). The origins of the body frame and flow frame sit at the glider's center of buoyancy (CB), which is located at the geometric center of the glider. According to the above assumptions, the CB of the glider is fixed with the vehicle.
In light of the definition of the reference frames, the rotation matrix between the inertial frame and body frame can be expressed as: where is a rotation matrix and describes the body frame relative to inertial frame. The matrix has the property like , where is the rotation mat-    rix from inertial frame to body frame, and are Euler angles. The notation is and is . Similarly, the rotation matrix from the flow frame to the body frame is where the rotational matrix also has the property like , the angle is defined as attack angle and the angle is defined as sideslip angle.
Based on the three defined frames, the kinematic equation of the glider can be presented as (Fossen, 1994): where denotes vector of position and attitude of the glider, and , are the position and orientation vectors, respectively. represents velocity vector of the glider, where is the translational velocity of the origin of the body frame and is the angular velocity of the body frame.
is the transformation matrix to transform vectors from the body frame to the inertial frame, where and can be expressed as: Expand Eq. (5) in detail: t· tan(·) where the notation is .
For bottom-heavy design, which is similar to other underwater vehicles, there is an offset from the CB to the CG of the flying-wing UG. The offset vector from the origin of the body frame to the CG of the rigid body (mass ) with respect to the body frame can be denoted as . The CG of the rigid body is moving with a translational velocity : Similarly, the position of three moving point masses in body frame can be expressed as , , . The kinematic equations of these internal actuators are is the standard Euclidean basis for .
3.2 External forces and moments acting on the glider The external forces and moments acting on the glider include gravity, buoyancy, hydrodynamic forces and addedmass forces. The dynamic equations of the multi-body system would be expressed in the body frame. Hence, the coordinates of the forces expressed in the body frame would be convenient, and this is what we will discuss in this section.

Gravity and buoyancy
The gravitational and buoyant forces of the glider expressed in the inertial frame can be described as: (10) Based on the relationship between the body frame and inertial frame, we can obtain The gravity of pitch control mass, roll control mass, the buoyancy control system and the rigid body (mass ) can be written as with respect to the body frame, where the subscript * denotes p, r, b, and rb, respectively.
As the CB of the glider coincides with the origin of the body frame, the buoyancy will not generate any moment on the flying-wing UG with respect to the body frame, but the gravity will. Therefore, the restoring force and moment in body frame can be expressed by the following vectors where is the offset vector from the origin of the body frame to the CG of the glider with respect to the body frame and .

Hydrodynamic forces and moments v V
According to the definition of the flow frame, the relationship between the attack angle, the sideslip angle, the velocity and the magnitude of the velocity is V= The hydrodynamic forces and moments we discuss in this section include hydrodynamic position forces and damping forces. The hydrodynamic forces in the flow frame are calculated with a linear approximation and can be expressed as: where is the density of seawater, represents the characteristic area of the glider, is the length of the glider, and are linear hydrodynamic coefficients. In addition, , and .

B W R
The rotation matrix maps the hydrodynamic force and moment from the flow frame to the body frame as:

Added-mass forces
Generally, the effect of the fluid surrounding the glider on the glider cannot be neglected and the effect is defined as added mass forces. Referring to Fossen (1994), the energy approach based on the Kirchhoff's equations is commonly used to derive the expression of the added mass forces. The kinetic energy of the surrounding fluid is written in the quadratic form: where is the system inertia matrix of added mass terms, and we have . is the added mass matrix, is the added inertia matrix, and is the added cross term.
Then, using the Kirchhoff's equations for Eq. (19), the added-mass forces acting on the glider with respect to the body frame can be calculated as: ω) The added Coriolis-centripetal matrix is calculated as: where character denote the skew-symmetric matrix satisfying for any .

Dynamics
(mass m rb )

Rigid body
The flying-wing UG is treated as a multi-body system. The dynamic equations of the rigid body (mass ) are similar to those of the conventional AUV. As in Fossen (1994) and Fossen (2002), using the Kirchhoff's equations or based on the Newton-Euler method, the dynamic equations of the rigid body (mass ) can be expressed in the body frame as: is a vector of the total forces acting on the rigid body (mass ) without the added mass forces and is a vector of the total moments acting on the rigid body (mass ) without moments caused by the added mass forces. The added mass forces and moments are put on the left of the equations. The matrices are: where is the identity matrix, is the moment of inertia of the mass with respect to the body frame and . is constant in time and can be measured by the 3D modeling software. The rigid body (mass ) is a submerged vehicle, so the external forces and moments discussed in Section 3.2 all act on the rigid body (mass ). Besides these external forces and moments, there are forces and moments acting on the rigid body by the three point masses, as calculated by where , and are respectively the forces acting upon the rigid body mass by the three point masses with respect to the body frame.

Pitch control mass and the roll control mass
The pitch control mass which can slide along theaxis of the body frame is treated as a moving particle. The rigid body (mass ) exerts a force on the pitch control mass . The sum of all forces acting on the particle is defined as , and can be expressed as: x v pr =ρ px e The body frame is a moving frame, and the pitch control mass can slide along the -axis with the velocity relative to the body frame. The convected velocity of the pitch control mass is . According to the Newtonian mechanics, the relative acceleration , convected acceleration and Coriolis acceleration can be calculated as: 27) According to the theorem of composition of the acceleration of a particle, the absolute acceleration of the pitch control mass is +ρ px e 1 + 2ω ×ρ px e 1 .
(28) Using the Newton's second laws, we finally yields the dynamic equations of the pitch control mass in matrix form: However, to obtain the complete dynamic equations of the multi-body system in the later work, we need to supplement a set of equations by premultiplying to each term on the both sides of Eq. (29) by matrix . Combining the result with Eq. (29) yields The matrices used in the above equation are: (32) Similar to the pitch control mass, the dynamic equations of the roll control mass can be expressed as: The matrices are:

Buoyancy control system
x For simplicity, we consider the case of the piston-cylinder arrangements to model the dynamic equations of the buoyancy control system. Let denote the position of the piston face and denote the length of the piston cylinder, see Fig. 4. The other parameters of the buoyancy control system are shown in Table 1. When , the vehicle is neutrally buoyant. When the buoyancy control system pumps the water into the piston cylinder and , the vehicle tends to dive down. When the buoyancy control system pumps water from the piston cylinder to outside and , the glider is buoyant and tends to climb up. Therefore, we treated the buoyancy control system as a particle with variable mass which can move along the -axis.
According to general equation of variable-mass motion, we can obtain the dynamic equation of the point mass : where is the force caused by the water pumping out or pumping in. However, this force is very small compared with the gravity, so can be neglected in our model.
ll m b where and are too small and are zero during most of the time when the glider is working, so they can be neglected while the variety of the position of the mass can not. By substituting Eq. (37) into Eq. (36), we obtain the dynamic equations of the buoyancy control system in the body  ZHAO Liang et al. China Ocean Eng., 2019, Vol. 33, No. 3, P. 322-332 327 frame: Similar to the pitch control mass, we supplement a set of equations and combine the result with Eq. (38), and we obtain: The matrices used in the above equation are:

Complete dynamic equations of the multi-body system
By combining Eq. (22), Eq. (30), Eq. (33) and Eq. (39), we can finally obtain the complete dynamic equations of the multi-body system: , .

Modeling and analysis of the steady gliding motion
The gliding motion is the principal working pattern of the glider. By changing its net buoyancy with the buoyancy control system and controlling the pitching angle using the pitch control mass, the vehicle could follow a sawtooth motion pattern in the vertical plane (as shown in Fig. 5). In this section, we will discuss the steady gliding motion and investigate the relationship between the controlled variables of glider and its motion performance. Our aim is to identify general trends to provide guidelines for the designs of the vehicle and controller.

Equation of motion in the vertical plane
When the vehicle is gliding in the vertical plane, the roll angle , the yaw angle , the translational velocity and the angular velocity are all zero. Therefore, the kinematic equations (Eq. (7)) can be simplified as: The buoyancy control system and the pitch control mass move with a low and uniform velocity relative to the glider, and they are static in most of the time during the gliding motion. Hence, the higher order terms of state variables in the dynamic model can be omitted, and the dynamic equations (Eq. (41)) can be simplified in the vertical plane as: . When the glider is gliding in a stable state in the vertical plane, the acceleration and angular velocities are all zero (Fig. 6). Then, the equilibrium equations in the vertical plane can be expressed in the flow frame as: where , and . ξ Noting from the Fig. 6, the glide angle can be defined as: (47) tξ = sξ/cξ s 2 ξ + c 2 ξ = 1 Using the fact that and , we can simplify Eqs. (44) and (45) as: Considering the fact that , Eq. (46) can be rewritten as: In our case, is 0.07 m and the pitch control mass is 200 kg. The mass of the buoyancy control system is initially set at 76.5 kg, which is denoted as . Thus, the relationship between the point mass and the net mass can be described by . From Fig. 4, the relationship between and the net mass is where .

Motion analysis of the steady gliding motion
The steady gliding motion of the flying-wing UG in the longitudinal plane can be completely described by Eqs. (47)-(51). Examining these five equations closely and counting the number of the unknown, dependent variables in these equations, we can find that there are seven unknown variables and in these equations. Note that, and are controlled variables and the others are state variables. The controlled variables are set by us, and these parameters can be calculated by MATLAB through the equilibrium equations. The relationship between the controlled variables and the other state variables can be estimated by the generalized least square method and the results are shown in Fig. 7. Referring to Fig. 5, the glide ratio is defined as The results indicate that in the steady gliding motion of the flying-wing UG, the attack angle first decreases fast and then slowly, and the glide angle shows an approximately linear increasing trend when the pitch control mass moves forward. The variation of glide ratio is similar to that of the attack angle. As the buoyancy control system putting water into the system and the pitch control mass moving forward, the velocity increases.
Referring to Fig. 8a, when the net mass is 25.5 kg, the horizontal velocity of the glider in steady gliding motion shows an approximately linear increasing trend as the pitch control mass moves forward, while the glide ratio decreases fast. When the pitch control mass is fixed ( =100 mm), we find that the horizontal velocity of glider first increases rapidly and then approximates linear growth as the buoyancy control system puts the water into the system, while the glide ratio remains nearly steady; as shown in Fig. 8b. From these results we can find that if the effect of the buoyancy control system on the glide ratio can be neglected, the buoyancy control system can be treated as a Fig. 7. Relationship between state variables and controlled variables in the steady gliding motion. ZHAO Liang et al. China Ocean Eng., 2019, Vol. 33, No. 3, P. 322-332 329 particle with a fixed position.

Numerical simulation
The actual motion of the flying-wing UG can be divided into gliding motion and spiraling motion. Based on the dynamic equations we derived, numerical simulation is conducted to evaluate the motion performance of the flyingwing UG in these two classical motion modes. The main purpose of this section is to analyze the basic motion performance of the glider, and the effect of ocean currents on the glider is ignored.

Numerical simulation of the gliding motion
The initial conditions of the flying-wing UG are set as =[0 m, 0 m, 10 m, 0 rad, 0 rad, 0 rad, 0 m/s, 0 m/s, 0 m/s, 0 rad/s, 0 rad/s, 0 rad/s] T . Setting the controlled variables and as 25.5 ∆m ρ px kg and 98.8 mm respectively, the glider dives down in the water. By changing the control variables and to -25.5 kg and -66.3 mm, the glider climbs up. The simulation time is assumed to be 8000 s.
As illustrated in Fig. 9, the glider follows a sawtooth motion pattern with a maximum velocity of 0.583 m/s. The pitch and attack angles are kept at -3.526° and 1.735° when the vehicle moves down, and the glide path ratio of the glider exceeds 10.

Numerical simulation of the spiraling motion
Other than the gliding motion, another classical motion mode of the glider is the spiraling motion. When the vehicle enters a spiraling motion, the accelerations are all zero and the roll angle, the pitch angle, the translational velocity and the angular velocity are all constant. Besides, the angular velocity expressed in the inertial frame is parallel to the   The initial conditions in this simulation are the same as the gliding motion simulation. By setting the controlled variables , and as 25.5 kg, 98.8 mm and 100 mm respectively, the flying-wing UG is made to enter the spiral-m r ing motion. In our case, the roll control mass is 200 kg and the simulation time is assumed to be 8000 s. Fig. 10 shows that the flying-wing UG we study performs a spiraling motion with a turning diameter of 692.7 m. The glider moves at a velocity of 0.6457 m/s, and the attack angle is 1.5137°.

Conclusions
In this paper, a flying-wing UG with a symmetrical airfoil is treated as a multi-body system and a nonlinear dynamic model of the glider is developed using Newton-Euler method. Considering the differences between the legacy gliders and the flying-wing UG, the buoyancy control system of the flying-wing UG is treated as a moving particle with variable mass, and the roll control mass is modeled as a moving particle which is similar to the pitch control mass. Then, a series of theoretical analyses of the steady gliding motion is presented to find the relationship between the state parameters of the steady gliding motion and the controlled variables. The results can provide guidelines for the designs of the vehicle and controller. Based on the dynamic model we derive, the dynamic simulation of the glider in gliding motion and spiraling motion is created. The simulation results show that the glide path ratio of the glider we discuss exceeds 10 when the glider glides with a maximum velocity of 0.583 m/s. These results demonstrate that the flying-wing UG we designed has an excellent motion performance.
The analysis of the steady spiraling motion is more complex than the steady gliding motion. Therefore, the relationship between the state parameters of the steady spiraling motion and the controlled variables is not studied in this work. Besides, one main limitation of this work is that the effect of ocean currents on the vehicle dynamic is ignored. In future work, the stability and control algorithms of the flying-wing UG will be discussed with consideration of the impact of ocean currents.