Locomotor transition: how squid jet from water to air

The amazing multi-modal locomotion of flying squid helps to achieve fast-speed migration and predator-escape behavior. Observation of flying squid has been rarely reported in recent years, since it is challenging to clearly record the flying squid’s aquatic-aerial locomotion in a marine environment. The existing reports of squid-flying events are rare and merely record the in-air motion. Therefore, the water-air locomotor transition of flying squid is still unknown. This paper proposes the idea of using CFD to simulate the process of the flying squid (Sthenoteuthis oualaniensis (S. oualaniensis)) launching from water into air. The results for the first time reveal the flow field information of squid in launching phase and show the kinematic parameters of flying squid in quantification. Both a trailing jet and pinch-off vortex rings are formed to generate launching thrust, and the formation number Lω/Dω is 5.22, demonstrating that the jet strategy is to produce greater time-averaged thrust rather than higher propulsion efficiency. The results also indicate that the maximum flying speed negatively correlates with the launch angle, indicating that a lower launch angle could result in a larger flying speed for the flying squid to escape. These findings explore the multi-modal locomotion of flying squid from a new perspective, helping to explain the trade-off strategy of water-to-air transition, and further enhance the performance of aquatic-aerial vehicles.

The mantle length l t The total body lengtḣ m (t) The instantaneous mass transfer at time t M b The body mass M e (t) The ejected water mass at time t M t (t) The total squid mass (contains water) at time t Q The value of Q-criterion r m The greatest mantle radius r m-max The maximum r m r m-r The resting mantle radius (r m at rest) s o The cross-sectional area of orifice t The time T The jet cycle T j The thrust generated by jet momentum flux T p The thrust generated by over-pressure impulse T a The thrust generated by added mass variation v i , v j , v l (i, j , l = 1, 2, 3) The fluid velocity vector components v jet (t) The jet velocity at time t v s The squid velocity V The fluid velocity V max The maximum water volume in the mantle cavity α o The orifice attack angle δ The mantle thickness Locomotor transition: how squid jet from water to air

Introduction
The multi-modal locomotion of marine animals and seabirds, including plunge-diving [1,2] and water-air launching [3], helps to obtain prey and avoid predators. Squid are special compared with the water-air transition animals, such as flying fish [4,5], gannet [6][7][8][9], cormorants [10], etc. The squid's strong mantle cavity [11], the soft and flexible arms, and the muscular lateral fins work in coordination to provide efficient propulsion [12]. Therefore, squids can reach the fastest speed (around 8 m s −1 ) known among aquatic invertebrates [13] and the speed of flying squid can reach up to 43.9-49.7 BL s −1 in the air as measured by body length (BL) [14], while the speed limitation of common fish is 25 BL s −1 [15]. Therefore, squid has always been recognized as excellent swimmer and escaper in the ocean. The study of the morphology and motion mechanism of squid has been a hot topic [16][17][18][19][20]. Different from jet propulsion underwater, during which the mantle alternately fills with and ejects water, squid flying activity is technically rocket propulsion because they cannot refill the mantle cavity in the air [21]. Kinetic energy with high power density can be generated by rocket propulsion, thus assisting the squid's aerial flight locomotion. The researchers are always interested in the squid flying behavior [22][23][24][25]. In 2004, Macia et al documented fourteen events of squid flight from 1892 to 2003, identifying that at least 8 squid species could fly and use fins and arms to control flight attitude [25]. The above studies are all qualitative reports and studies via photos or visual observations. Some other researchers have also tried to study the flight behavior of squid quantitatively by sequential photography or video (stacks.iop. org/BB/15/036014/mmedia). In 2013, Muramatsu et al took sequential photographs in the Northwest Pacific (35°34.0′N, 146°19.3′E) and defined flight process based on photographic observation of oceanic squid (Ommastrephidae) [14]. According to the analysis of these photographs, the average mantle length of the squid was around 60.1% of body length, and the speed after launching ranged from 8.8 to 11.2 m s −1 (43.4-49.7 BL s −1 ), the total flying distance was around 26.4-33.5 m. In 2010, Dr kubodera and Dr Ohizumi recorded the whole process of squids' flight behavior in the air, which was believed to be the only video recording squid flight [26]. By analyzing the video, Odor et al calculated the speed, acceleration, flight distance, and other kinematic parameters of the squid in the aerial flight phase in 2015 [26]. However, higher quality research resources are still scarce because it is difficult to capture the flight behavior of squid by taking videos or sequential photos in a marine environment.
According to the biological observations, S. oualaniensis, commonly known as the purpleback flying squid, which mainly distributes in Indian Ocean and Pacific Ocean, is considered to be a stronger squid with flight capabilities [27]. In general, the squid multimodal locomotion can be divided into six phases, i.e. roaming, accelerating, launching, jetting, gliding, and diving as shown in figure 1. However, current researches mainly focused on speed, acceleration, thrust, nozzle angle, mantle volume [28,29], or vortex propulsion [30][31][32][33][34] etc, through biological observation on underwater jetting or aerial flight separately [14,26]. To the best of our knowledge, research on squid launching phase is almost blank due to the difficulties to capture the water-air transition movement. Benefiting from the development of CFD technology, most of issues in biomimetic area that are impossible to solve by traditional methods now can be investigated through numerical simulation [35][36][37][38][39]. The CFD technology can supply efficient numerical experiments to simulate the motion of research objects under the condition of the theoretical models developed properly [37]. Therefore, CFD methodology is a feasible way to simulate the multi-modal locomotion of flying squid in launching phase.
In this paper, in order to quantitatively investigate the kinematics and propulsion during the squid launching phase as shown in figure 1(c), several numerical experiments were developed to simulate the squid motion. A 3D squid model was designed based on the biological morphology of S. oualaniensis. Initial parameters used in the simulation, including jet velocity, squid mass, water cavity volume, etc, were properly determined according to references. This paper investigates the jet strategy of the squid in launching phase. The complete kinematic process of the flying squid's water-air launching is reproduced, analyzed and revealed quantitatively for the first time. The influence of different launch angles is studied. By examining and reconstructing the morphology of squid, as well as investigating squid's launching behavior through CFD simulation, this paper offers a new perspective on multi-modal locomotion research of flying squid, thus promoting research methodology of such marine creatures and further enhancing the performance of aquatic-aerial vehicles [40].
The remainder of this paper is organized as follows. Section 2 presents the biological prototype, initial parameters and three launch angles are introduced. In section 3, the numerical methodology is illustrated. In section 4, the experiment results are discussed. At last, the conclusions are drawn in section 5.

Biological prototype
There are few direct observations of squid flying events. Ommastrephidae including Ommatrephes bartrami, Todarodes sagittatus and S. oualaniensis are commonly reported in these existing records. S. oualaniensis, known as the purpleback flying squid, has excellent flying ability [26]. In a biological observation record, the speed of S. oualaniensis ranged from 8.8 to 11.2 m s −1 (43.4-49.7 BL s −1 ) during flight and the flying distance was estimated to be 26.4-33.5 m [14]. In addition, S. oualaniensis, the average speed of which can reach 6.76 m s −1 in jetting phase (figure 1(d)), flies faster and farther compared to S. pteropus and Dosidicus gigas [26]. Therefore, to make the research more representative, the S. oualaniensis was chosen as experimental research object.
For S. oualaniensis living in the eastern tropical Pacific Ocean, the mantle length and mass are related to its age [41]. Therefore, the mantle length of S. oualaniensis can be described as, where l m is mantle length of the squid, Age denotes squid's life time, the unit of which is day. In addition, the body mass of S. oualaniensis can be presented as, where M b is body mass of the squid. From equations (1) and (2), we set the Age of our biological prototype as 84 days, then the mantle length and body mass of the prototype can be calculated as 148 mm and 131.5 g, respectively. Besides, total length l t of squids is around 1.5 times longer than mantle length, so the total length of our prototype is 222 mm [28]. The schematics of biological prototype's key geometry parameters are shown in figure 2. Studies have shown that the squid's orifice diameter D o , funnel orifice angle of attack α o and orifice shape vary during launching phase [28,29], which are necessary for posture adjustments and balance keeping of the body. Furthermore, to produce enough thrust, water in the mantle cavity of the squid is constantly ejected, so the squid's volume decreases until the water jetting finishes. However, it can be known that soon after the startup of jetting, D o and α o of squid tend to be a stable value, and the orifice shape turns into a circle [29]. Since the squid orifice is not coincident with body axis but on the bottom edge, if the direction of the jet is completely parallel to the body axis (α o = 0°), the jet thrust will generate a large pitch moment. It Then a small amount of water is sucked into the mantle cavity and ejected backwards. At the same time, the arms stroke backwards, helping to provide a stronger propulsive force to accelerate. (c) Launching: the squid rolls down their fins, folds their arms, inhales sufficient water, and breaks the water surface by jet propulsion. (d) Jetting: the squid spreads their arms and fins to adjust flight attitude and to provide lift, and continue jetting until the water in the mantle is exhausted. (e) Gliding: the squid stops jetting water, keeps its fins and arms spread, and glides utilizing the kinetic energy obtained from jetting. (f) Diving: the squid rolls down its fins and folds the arms again when approaching the sea surface, streamlines the shape to minimize water resistance, and then enters the water. is reported that during squid high speed underwater motion, the orifice attack angle is approaching −3° [28] which can help reducing the pitch moment. Therefore, we set orifice diameter and angle of attack as fixed values, where D o is set as 4 mm, and α o as −3° with respect to body axis. In addition, for a squid of which l m is 0.135 m, the greatest mantle radius ranges from 9.4 mm to 16.0 mm during jetting [28]. The body deformation of squid is mainly reflected on the mantle radius variation with time. Because of the large drag coefficient of water, the influence of the mantle radius changes on the squid launching phase mainly occurs underwater. The squid mass varies from 0.25 kg to 0.23 kg during the underwater motion, i.e. the mantle radius varies from 16 mm to 14 mm. Therefore, quite little amount of water is ejected during the underwater motion period. The influence of the body geometry variation on launching phase is tiny. To simplify calculation and ensure accurate convergence of results, we assume that the mantle radius of the squid remains unchanged during the CFD simulation, thus the greatest radius of our proto type r m was set as 15 mm. Table 1 shows the specific parameters of the biological prototype. Figure 3 shows the 3D geometry model of squid constructed in SolidWorks 2016 (Dassault Systemes, Vélizy-Villacoublay, France). The shape of the entire model refers to the streamlined body of the S. oualaniensis, and the morphology and posture are designed exactly following a live flying squid during launching phase. In nature, the mantle lip opens under muscular control to refill water into the mantle cavity, during which the orifice is closed [42]. After refilling, the squid closes its mantle lip and folds the fins and arms around the mantle, then gets ready to launch. To minimize the resistance, the shape of the squid turns into a spindle in launching phase as shown in figure 3. Meanwhile, the fins on the tail are tightly wrapped around the mantle. As the squid's body orientation is tail forward during launching, this behavior greatly reduces the resistance. In addition, the ten arms of the squid close up along the head direction. When scared or threatened by predators, flying squid will break the sea surface by jetting water backwards in a posture shown in figure 1(c) [14]. The 3D geometry model was established according to the body parameters shown in table 1.

3D geometry model and coordinate system definition
The coordinate system is presented in figure 3. O-xyz is the ground coordinate system, where O is the end point of the squid arm at the beginning. The xOz plane is parallel to the water surface, and y -axis is defined according to the right-hand rule from z-axis to x-axis. C-uvw is the body coordinate, and C is the squid's center of mass. The uCv plane, the symmetrical plane of the geometry model, is coincident with the xOy plane. The u-axis is along the body axis toward the tail, and w-axis points the same direction as z-axis. So v can be defined following the right-hand rule from w-axis to u-axis. Therefore, C-uvw can be obtained by transforming O-xyz by rotating along the direction z with angle of γ (launch angle, refer to table 2).  Bioinspir. Biomim. 15 (2020) 036014

Jet flow parameters
For fluid simulation, several input parameters are necessary, including jet-flow speed, squid's body weight with water, etc. The speed of squid multi-modal locomotion varies with different states, such as roaming and escape jetting. When roaming underwater, the swim propulsion is generated by fins and water jets, and the squid's mantle cavity inhales the appropriate amount of water and sprays backwards slowly [12]. However, for the escape jetting state, the mantle cavity expands extremely to store a large amount of water and jet in a high speed, which can provide a considerable thrust [29]. For a complete launch process of the squid, jet-flow speed, and body weight (contains water) will change over time.
Water volume ( V max , cm 3 ) in the mantle cavity for the squid is estimated by treating the squid mantle as a single cylinder with a variable diameter and length less than total mantle length. With corrections on mantle muscle and viscera, the variation of the cylinder radius can be used to simulate the uneven deformation of the mantle cavity [28], where V max is maximum water volume in the mantle cavity, r m-max represents maximum r m , r m-r denotes resting mantle radius (r m of the squid measured from dead specimen), δ is mantle thickness, and l m means mantle length. M b represents squid mass, and d s is squid density. In our study, M b = 0.131 kg; l m = 0.148 m; r m-r = 0.0138 m; r m-max = 0.03 m; d s = 1055 kg m −3 . Therefore, for the squid modal adopted in our experiment, V max is 0.127 kg derived by equation (3).
The jet velocity can be determined by the mantle cavity geometry change with time. However, capturing such change during squid's launching phase is quite difficult. Therefore, the flying squid jet velocity used in simulation was obtained by referring to the squid underwater escape jet velocity law curve, because the water-to-air launch of squid can also be regarded as an escape behavior [29]. According to [29], the squid jet velocity is calculated by time varying flow and jet area. The result of the research shows that a large initial peak exists in the velocity curve. After the peak, speed decreases sharply to a small value, then increases stably to a large one, and then gradually decreases to zero, which is like a sine curve. The peak is mainly caused by the pressure generated in the mantle before jetting. That is, at the moment of orifice opening, the overpressure in the mantle cavity generates an instantaneous impulse of jetting flow. Then, the change of the velocity curve afterwards is controlled by the contraction of the squid mantle muscle. The relationship between jet velocity v jet and maximum mantle water volume V max can be expressed as follows, where ρ w is the water density, s o represents the crosssectional area of orifice, and T denotes jet cycle which is estimated as 1s according to the research of O'dor et al by biological observation [26]. Then, the jet velocity v jet can be calculated by the following equation, where V is the water volume in the mantle cavity. Therefore, the instantaneous jet velocity can be obtained as shown in figure 4(a). Then, the jet flow was added on the orifice of the squid model by UDF (user defined function) in Fluent. When the squid moves throughout the computational space, the jet region will move with the orifice to simulate the squid's waterjet process. In addition, the body weight of the squid is also changing as the mantle cavity jets water continuously, which will affect the inertia and gravity of the squid. The ejected water mass M e can be obtained by integrating the jet velocity, which can be presented as, Furthermore, the change in squid total mass is determined by the amount of ejected water, so the squid total body mass can be expressed as, where M t is the squid total mass (contains water) at time t. Thus M t and M e changing over time during launching phase can be described as figures 4(b) and (c), respectively. Then the squid total weight M t attribute is defined in UDF according to figure 4(b). Therefore, the Fluent 6-DOF solver will calculate effect of M t loss on the squid motion.

Launch angles and assumptions
Flying squid has excellent flying skills through years of biological evolution and has established a perfect flying strategy trade-off mechanism. Flying squid using their fins, arms and probably their fin flaps to adjust flying posture like flying fish during the jetting and gliding phases [21]. For the launching phase, because of instantaneity and arms and fins being folded, the launch angle is the main factor that dominates the movement of the flying squid [43]. Therefore, the influence of different launch angles Mantle cavity V max (ml) 30 15 a Velocity profile f (t) refers to figure 4(a).
on the launching phase was investigated. As shown in figure 5, the simulation cases were divided into three conditions by the different launch angle γ, i.e. the angle between the body axis of the squid and the water surface which were set to 30°, 45° and 60°, respectively. The COMs (center of mass) of the squid prototype in the three cases were arranged on the arc with O as the center. Meanwhile, the squid body axis in the three cases converged at the point O. In addition, all the three cases share the same biological prototype, jet propulsion parameters, and environmental setup as shown in table 2.
By analyzing existing records of flying squid events, the squid launches out of water in a short time. To reduce pitch moment induced by lateral loads (especially when moving underwater), the body longitudinal axis and velocity vector of a launching vehicle are consistent during each launching phase [44,45]. The same principle applies to the squid motion as well. If the squid body axis has a small angle with the advancing direction underwater, the squid, the shape of which is similar to a spindle in the launching phase, will generate a great pitch moment causing the body to rise or bow, resulting in uncontrollable pitching. During the squid launching process, the squid longitudinal axis and the motion direction are also highly consistent by self-adjustment, which can be observed by videos. Therefore, under the condition of following the real motion of flying squid, the squid pitch motion was constrained during simulation, under which condition the change of the center of gravity can be ignored.
In addition, the underwater buoyancy was calculated by Fluent and was considered constant. As mentioned in section 2.1, the squid mass varies from 0.25 kg to 0.23 kg during the underwater motion, i.e. quite little amount of water is ejected during the water motion. The squid body buoyancy variation is within 0.1 N (squid model tonnage is 0.24 kg), which is about 3% of the average value of T j . Therefore, the influence of the body buoyancy variation on launching motion is tiny. To simplify calculation and ensure rapid convergence of results, the mantle radius of the squid was  The launch angles are set as 30°, 45° and 60°, respectively. For these three conditions, the profile of jet velocity, squid mass, and ejected water mass are identical and all subject to parameters mentioned in figure 4. treated as unchanged, i.e. the buoyancy was derived by hydrodynamic calculation of Fluent and remain unchanged.

Thrust generation
The value of the thrust has a direct impact on the movement of the squid. Jet momentum flux(predicted by slug model [46]), orifice over-pressure [32] and added mass variation [47][48][49][50] contribute to produce the squid thrust.
Assuming that the orifice ejects fluid in a quiescent environment, the thrust generated by jet momentum flux can be predicted by slug model [46]. The jet velocity in the orifice area can be described by the mass transfer flow as [51], where ṁ is the instantaneous mass transfer in the orifice, A denotes the orifice area. And the thrust generated by the jet momentum flux T j (t) is presented as, where D o can be referred from table 1, ρ w = 10 3 kg m −3 , and mass transfer can be expressed as During the simulation, the manually added jet flow cannot directly generate the force on the squid body, the thrust T j was applied to squid body through UDF after external calculation according to equation (9). However, researches have shown that the slug model ignores the pressure contribution to vortex ring formation because it makes no distinction between the trailing jet and the vortex ring [32]. The total impulse is mainly determined by the jet momentum flux term and over-pressure at the orifice exit plane. The total impulse can be presented as, where I j is impulse due to jet momentum flux, and this part of thrust can be derived by equation (9). I p denotes the impulse caused by the over-pressure at the orifice exit plane, which is defined by, p ∞ represents the ambient pressure at (x, r) → (∞, 0). The thrust generated by over-pressure impulse, T p , can be calculated as follows, Equations (10)- (12) are already embedded in Fluent solver program. Therefore, T p was computed real-time during simulation in Fluent.
In addition, previous studies have shown that added mass variation play an important role in thrust generation. If the object passes through a fluid of which density cannot be ignored, the motion equation can be described as follows according to [47], where m a is the fluid added mass, and f d is resistive force. It is obvious that the size and the shape variations of the body are important in determining the propulsive performance. It can also be known that recovery of added-mass energy enables a shrinking rocket in a dense inviscid flow to achieve greater escape speed than an identical rocket in a vacuum or liquid with less density. Substitute the body geometry and jet flow parameters mentioned in sections 2.1 and 2.3 into equation (13), the maximum value of the thrust generated by added mass variation T a is about 0.031 N during underwater part of the squid launching phase, and when t > 0.05 s, T a is less than 0.01 N. However, the maximum value T j is about 9.2 N as calculated by equation (9) and the maximum ratio of T a /T j = 0.52% at t = 0.02 s. In addition, as discussed in section 2.1, most of the water inside squid mantle cavity is ejected during aerial motion. It is also known that if the object is operating in fluid of a lower density, the external forces may be neglected [47]. Therefore, to simplify the simulation, the thrust generated by added mass variation is neglected.

Solver
Choosing a suitable and efficient CFD solver can expand the research scope of the experiment and improve the simulation efficiency and accuracy. According to the different discrete principles of the governing equation, the numerical solution of CFD produces several methods, finite difference method (FDM), finite element method (FEM), finite volume method (FVM), etc. Some discrete methods (such as FDM) only satisfy the integral conservation when the mesh is extremely fine. Therefore, the discrete equations derived by FVM are the most commonly used numerical methods due to the better accuracy in the case of coarse meshes and better integral conservation property. For fluid flow, the conservative governing equation of the physical quantity φ in the control volume can be expressed as follows, where ρ is the fluid density, V denotes the velocity, Γ represents the diffusion coefficient, and S φ is the source term. The FVM is to integrate the governing equation shown in equation (14) within the control volume. That is, to solve the governing equation in integral form as shown below, where Ω V is the control volume. The above equation can be simply described as follows, where φ t , φ A , φ г and φ s are net increase in physical quantity caused by transient, convection, diffusion and source terms, respectively. By this method, the nonlinear partial differential equation is transformed into an algebraic equation group on the control volume, and then the algebraic equations are solved to obtain the solution of the flow field. Currently, commercial CFD software with good versatility such as Fluent, CFX, STAR-CD, etc all adopt the finite volume method. Since the Fluent has a built-in six-degree-of-freedom (6-DOF) solver, the movement of the squid can be calculated based on hydrodynamic/aerodynamic forces, gravity or other moments acting on the squid in the flow field. In addition, with the Fluent dynamic mesh technology and the multiphase flow model, the multi-modal locomotion of the flying squid can be well simulated by Fluent. Moreover, since the movement of flying squid is a 3D process, the Fluent 16 3D (single-precision) using the finite volume method is selected as the solver.

3D meshing and boundary conditions
Since the body-fitted grid makes the computational domain coincide with the mesh boundary, it can help to improve the convergence speed of the numerical solution and the resolution. Therefore, the bodyfitted coordinate (PCBFC) is used for meshing. The unstructured tetrahedral mesh is used, considering its effectiveness when adapting to complex solution domain with irregular shapes and curved boundaries, especially for 3D problems. Generating the mesh, specifying the boundary conditions, etc are all done by the preprocessor GAMBIT.
The rectangular cube computational domain grid is shown in figure 6(a). The purple slice perpendicular to the water surface coincides with xOy plane in the ground coordinate system. The green slice passes through the squid center of mass, and is parallel to the xOz plane. The blue area at the center of the entire computational domain is the initial location of the squid body. The 2D magnified view of the slice when launch angle γ = 45° is shown in figure 6(b). The slice is located in the xOy plane passing through the body axis of the squid. Since the gradient of physical quanti ty in the normal direction and the tangential direction of the region close to the squid surface changes drastically, the mesh near the wall of the squid is densified. The total number of grid is around 5.21 million.
During the simulation, grid independence has been verified. The grid independence was tested by comparing the drags of the squid moving in the air under v s = 5 m s −1 and in the water under v s = 1.02 m s −1 , respectively. Three different grid densities of 2612 346 (coarse), 5209 678 (medium) and 8024 826 (fine) nodes were used for calculation. In underwater and air motion, the drag differences between medium grid and the fine grid are 1.94% and 2.4% (within 5% [35]), respectively. To save the computational consumption, the medium grid density was selected to implement the simulation.
The boundary condition together with the initial condition is called the definite condition. Only when the boundary condition is set, the solution of the flow field exists and is unique. The squid flies into air from water, and there is a multiphase flow. Therefore, the pressure outlet boundary condition is set at the waterair interface. For subsonic flow, at the water-air boundary, set the gauge pressure to 0 (relative pressure). To make the simulation easier to converge, set the backflow direction to normal to boundary. The backflow turbulence parameters on the boundary are defined by the intensity I and viscosity ratio µ t /µ L . Turbulence intensity I is defined as follows, where u′ is the root-mean-square of the turbulent velocity fluctuations and U is the mean velocity. Set the backflow turbulent to 1% at the boundary. The turbulent viscosity ratio µ t /µ L represents the ratio of the turbulent viscosity µ t to the laminar viscosity µ L , and the turbulent viscosity ratio is set to 10 in this simulation. For the near-wall flow, a wall function is used to simulate, and no slip condition is applied to the wall.

Dynamic mesh
For the simulation of moving and deforming flow, the moving reference frame (MRF), sliding mesh, and dynamic mesh methods are often used in Fluent. The dynamic mesh method updates and reconstructs the grid of the computational domain before each time step of iteration based on the motion of the object to calculate the unsteady FSI (fluid structure interaction) problem. Therefore, the dynamic mesh method provides an effective solution for the flying squid multi-locomotion simulation. The core of the dynamic mesh is the mesh generation. Firstly, the initial mesh is defined. Then, the Spring-Based Smoothing model is used to calculate the dynamic process of the mesh to ensure that the mesh topology does not change and no interpolation is needed, thus ensuring the accuracy of the calculation. In addition, the Spring-Based Smoothing model has good compatibility with tetrahedral meshes. In the dynamic mesh calculation, the time derivative term of the conservation governing equation shown in equation (15) can be expressed as a first-order backward difference format, where n and n + 1 are different time layers. Ω n+1 V on the n + 1 layer can be expressed as, where dΩ V /dt is the time derivative of the control volume, which can be calculated by the following formula according to the conservation law of the grid, where − → u g is the grid velocity, n f denotes the face number of the control volume, − → A j represents the area vector of the face j , and the dot product − → u gj · − → A j is calculated by, δΩ Vj is the volume swept by the face j in the time interval ∆t. In addition, the movement of the squid is given by the 6-DOF solver provided with the dynamic mesh in Fluent. The squid weight attribute ( figure 4(b)) and T j were added via UDF (F UDF ). When the simulation object has weight properties, the 6-DOF solver will calculate gravitational effect on the movement of the squid. The buoyancy, T p and other hydrodynamic forces were derived by hydrodynamic calculation of Fluent (F Fluent ), which will also be considered in the 6-DOF solver. Therefore, the 6-DOF solver calculates the position and attitude of the squid based on forces/ moments acting on the squid given by UDF (F UDF ) and hydrodynamic forces calculated by Fluent (F Fluent ). The governing equation for solving the translational motion of the center of gravity in the ground coordinate system is, where ˙ V presents the translational acceleration vector of the center of gravity, m is the mass, and F denotes the force vector caused by the external force.
The governing equation of the squid's rotational motion in the body coordinate system is: where L is the inertia tensor, − → M B is the moment vector of the squid, and − → ω B is the angular velocity vector. Converting − → M G in the inertial coordinate system form to the body coordinate system is expressed as: where R represents the transformation matrix, (25) where C λ = cos (λ), S λ = sin (λ), and angles ϕ, θ, and ψ are Euler angles rotating around the x, y , and z axes, respectively. The translational velocity and the rotational angular velocity can be calculated iteratively by equations (22) and (23), respectively, using which to update the grids position and attitude.
The fluid force applied on the squid is obtained by integrating the pressure and shear stress on the surface of the squid. Dynamic mesh with the 6-DOF model greatly simplifies and speeds up the calculation process.

Turbulence model
Although the Navier-Stokes equation can describe the details of turbulent motion, it consumes a large amount of computational resources due to the multiscale vortex motion and mathematically strong nonlinearity. For example, the number of grids required by direct numerical simulation (DNS) is about Re 9 4 (Re is Reynolds number) and integration steps are above 10 5 , which have strict requirements for computing speed and memory space. Since the simulation process only needs to be concerned with the change of the average flow field caused by turbulence, the Reynolds averaged Navier-Stokes (RANS) method is used to express the fluid motion as a sum of the statistical average value and the fluctuating value. For Newtonian fluids, after averaging, the governing equation becomes where u i , u j , and u k denote the Reynolds mean velocity components with the average sign omitted, ρ is the density, p is the pressure, µ is the dynamic viscosity, u i and u j are the components of turbulent velocity fluctuations, and δ ij is the Kroneck tensor component.
Here the equation introduces the correlation term ρu i u j (Reynolds stress term) of the turbulent fluctuating value. In order to solve the equations, it is necessary to make assumptions about the Reynolds stress by turbulence viscous model. The two-equation model standard k-epsilon (k − ε) model has a wide range of application, and it has low requirements on grid accuracy, and has self-contained wall correction, which has good robustness. Therefore, standard k − ε model is adopted as the turbulence model [52]. The k − ε turbulence model can be described as follows.
In addition, the value of Y + , which is the dimensionless distance from the centroid of the first layer grid to the wall, should be verified. Y + is the result calculated by solver and ubiquitous in turbulence problems. When meshing, the wall-adjacent cells are generally placed within the range established by the turbulent region outside the viscous sublayer, i.e. 30 < Y + < 300 [55], which can guarantee that the near-wall mesh resolution is acceptable. The Y + can be derived by the following equation [56], where y is the distance from the wall to the cell center, µ is the dynamic viscosity, ρ is the density of the fluid, and τ W is the wall shear stress. The simulation results indicate that for the squid wall-adjacent cells, Y + is about 90 and for most of these regions it does not drop significantly below 30 and satisfies the requirement.

Multiphase flow model
In the Euler-Euler method of the multiphase flow model in Fluent, the different phases are processed into interpenetrating continuums and thus introducing the concept of phasic volume fraction or relative volume fraction. The VOF (volume of fluid) multiphase flow model in the Euler-Euler method simulates two or more incompatible fluids by solving a same momentum equations and tracking the volume fraction of each fluid passing through the computational domain, which has a good applicability to the free surface flow [57,58]. Therefore, the VOF is adopted as the multiphase flow model. The volume fraction of the qth phase fluid in the VOF model can be expressed as a q (x, y , z, t), a continuous function of time and space. Therefore, when a q = 0, the control volume does not contain the phase q; when a q = 1, the control volume only contains the qth phase; when 0 < a q < 1, it indicates that the control volume contains both the qth phase fluid and other fluids. If the qth fluid is taken as an example, the continuity equation is, and the continuity equation of volume fraction can be expressed as, where t is time, V represents the velocity vector, and the sum of the volume fractions is 1 in the same control volume, i.e.
The multi-modal locomotion of flying squid only contains two fluids, i.e. water and air. Using a w represents the volume fraction of water, the density ρ and dynamic viscosity µ of each control volume of can be expressed as, where ρ a is density of air, and µ a denotes dynamic viscosity of air. Substituting ρ and µ into equations (26)- (28) to express multiphase flow motion. Furthermore, the geometric reconstruction method is used for tracking the free surface, and the direction and position of the interface are determined by geometric operations according to the value and gradient of the VOF factor in the control volume [59].

z-axis vorticity distribution and Q criterion
Vorticity criterion is widely used to recognize vortex regions [60]. The derivation of vorticity criterion is presented in the appendix. Figure 7 shows the z-axis vorticity w z contours of the flow field in the symmetrical slice at t = 0.22 s and the launch angle γ =45° when the squid is launching out of water. The slice coincides with the xOy plane of the ground coordinate system and passes through the body axis. The calculation of the vorticity is in accordance with the vorticity criterion in equation (A.2). The direction of the vortex follows the right-hand rule, with the positive direction in the vertical z-axis and the negative direction in the inward direction. The red part shown in the figure 7 is the positive vorticity, on the contrary, the blue present negative vorticity. In addition, the change in color expresses the change in vortex intensity (−5 to 5 s −1 ), with the red and blue portions being the largest, and the green portion being the fluid without spin.
According to the simulation results, the vorticity of the flow field is mainly composed of two parts, i.e. the part caused by the relative movement between the squid body and ambient fluid, and the part caused by the jet around the tail. The vorticity around the squid body is mainly caused by the viscous force of the fluid. When the fluid flows through the squid body surface (the solid side wall), a flow thin layer with a clear flow velocity gradient formed near the wall surface is called the boundary layer [61]. Since the fluid must satisfy the condition of no slip on the flow surface, the tangential velocity of the fluid on the wall must be zero, but as the distance from the surface increases, the velocity of the fluid must increase and gradually reach the speed of free flow. Hence, the vorticity around the body is largest and evenly distributed from the front to the rear of the squid body surface. The vorticity decreases as the normal distance to the surface increases. On the other hand, due to the influence of the jet, the vorticity direction behind the tail of the squid is exactly opposite to the part around body surface. At the front of the jet flow, the positive and negative vorticities are close to each other. However, at the rear of the jet flow, the positive vortices and negative vortices are separated and form a vortex ring. It should be noted that the magnitude of the vorticity does not absolutely characterize the intensity of the vortex. For instance, in the boundary layer, the dominant one is shear strain rather than vorticity magnitude, but the vorticity cannot identify them. Therefore, the area with larger vorticity in figure 7 is not necessarily a strong vortex area, but a comprehensive performance of shear strain and vorticity magnitude.
Q represents the balance between the shear strain and vorticity magnitude, i.e. the region where the vorticity magnitude is larger than the rate of strain [62,63]. Therefore, the vortex region can be determined by Q-criterion. The region of which second invariant Q of the velocity gradient tensor ∇ V in the flow field has a positive value can be defined as a vortex, i.e. Q > 0 [63]. Also, it requires that the pressure in the vortex region is lower than the ambient pressure. The 3D vortex identified by the Q-criterion is shown in figure 8(a). A pinch-off vortex ring is formed following the trailing jet. The isosurface when Q = 0.01 s −2 is colored by the velocity magnitude, which shows that velocity of the vortex ring is larger near the jet axis. The figure 8(b), the 2D vortex region, is colored by the value of Q. It can be seen that the whole part of trailing jet exists strong vorticity magnitude, and for the pinch-off vortex ring the Q varies from 0.4/s 2 to 0 /s 2 from the center of vortex ring to boundary. Both the trailing jet and pinch-off vortex ring help to generate thrust. In addition, the squid wakes have two different modes, which are characterized in terms of 'formation number' F, i.e. single vortex ring (F < 3, low in overall time-averaged force production but high in propulsive efficiency) and trailing jet with pinch-off vortex ring (F > 3, higher time-averaged force production but lower propulsive efficiency) [33]. F means the ratio of jet length (based on the vorticity extent (L ω )) to jet diameter (based on peak vorticity locations (D ω )). As presented in figures 7 and 8, the formation number L ω /D ω is 5.22, indicating that the squid jet strategy during launching phase is to produce greater time-averaged thrust rather than higher propulsion efficiency.

Kinematics
The trajectory, velocity, acceleration, and thrust of the squid movement from 0 to 1 s at different launch angles are shown in figure 9. By comparing the kinematic parameters when γ = 30°, 45°, and 60°, it can be found that launch angles have a significant impact on the motion of the squid launching phase. All kinematic parameters are obtained based on the model mass center. The motion of the flying squid equals the translation of the body coordinate system C-uvw relative to the ground coordinate system O-xyz. As shown in figure 9(a), in the simulation period, the locomotion of the squid contains two parts, underwater and in air. The squid flies out of the water with the ground coordinates (in meters) of (0.31, 0.25), which indicates that the underwater part is very short and the most of the trajectories are concentrated in the aerial part. Figures 9(b) and (c) present variation curves of velocity and acceleration during the launching phase, respectively. Through the time axis, it can be found that when the time is about 0.32 s, the squid firstly launched out of the water surface when γ = 60°, followed by 45° and 30°. During underwater stage, the velocity and acceleration are almost the same at different launch angles, but after jumping out of water, the lower launch angle has greater velocity and acceleration at the same time and the maximum acceleration reaches a peak of 26.53 m s −2 (approximate to the maximum acceleration of 30 m s −2 of natural squid [64]) at t = 0.51 s. This is because when moving underwater, there is buoyancy to balance gravity, so the work generated by the jet mainly contributes to speed increase. But when flying in the air, in a higher launch angle condition the gravity does more negative work reducing the flying velocity. In addition, the maximum flying velocity is 13 m s −1 (t = 0.94 s), 11.2 m s −1 (t = 0.90 s) and 10 m s −1 (t = 0.86 s) when γ = 30°, 45° and 60°, respectively. The results show that the maximum flying speed negatively correlates with the launch angle, indicating that a lower launch angle could result in a larger flying speed for the flying squid to escape.
Furthermore, the thrust T j shown in figure 9(d) at different launch angles is identical because they share the same mass transfer calculated by jet velocity shown in figure 4(a). It can be seen that the over-pressure in the mantle cavity leads to a large start-up propulsion, after that, the thrust is mainly controlled by the jet flow velocity. The maximum T j is 6.38 N when t = 0.48 s, which is 5 times gravity (1.27 N).
For the natural squid flying behavior, the video provided by Dr Ohizumi is believed to be the only reference-worthy record. By analyzing the frame sequence of the video and calibrating the feature points of the squid, the pixel distance of each frame image is obtained, and finally the data is calibrated and converted to obtain the speed and acceleration of the squid flying [26]. The results show that the entire jetting phase of S. oualaniensis lasted 0.61 s, and the flying distance was 4.12 m, and the average speed was 6.76 m s −1 . The compariso n between the CFD simulation results and data of the biological observations is shown in table 3, which shows that our simulation results are in good agreement with squid flying motion data from biological observations. Table 3 also indicates that the simulation results are valid and more comprehensive. Figure 10 presents the fluid forces exerted on the squid when the squid launches out of water, such as buoyancy and fluid drag, etc, in x-axis and y -axis., the results of which show that fluid forces are significantly different between water and air. This helps to explain the tradeoff strategies of multi-modal locomotion of the squid.

Fluid forces in different media
As shown in figures 10(a)-(c), during the underwater part, the water resistance of the squid increases along the negative direction of the x-axis as the squid speed increases. On the other hand, since the buoyancy in direction y + is greater than the water resistance in direction y − , the fluid force in the y -axis direction is positive during the underwater movement. For instance, when γ = 30°, the maximum fluid force is 1.5 N in the y -axis direction at 0.26 s, while the maximum value of fluid force in direction x-is 0.49 N at t = 0.31 s. Besides, around the water exit point on the fluid force curve in the x-axis, air gradually replaces water during the process of crossing the water surface, and the air resistance coefficient is much lower than that of water, so the absolute value of the fluid force is decreases continuously. For the fluid force in the y -axis, due to the rapid decrease of buoyancy, the fluid force in the positive direction of the y -axis also decreases rapidly until the buoyancy disappears, then the resistance of air becomes the main component of the fluid force value of which becomes negative. Figure 10(d) shows the resultant forces at different launch angles. It can be seen that during the whole underwater movement, the resultant force of a lower launch angle is greater than that of a higher launch angle. For example, the fluid resultant forces of 30°, 45°, and 60° are 1.51 N, 1.11 N, and 0.88 N around t = 0.25 s, respectively. However, when the squid leaves the water, the curves of different launch angles approach to similar small values which are the resistance forces generated by air. For instance, the fluid resultant forces of 30°, 45°, and 60° are all around 0.05 N after t = 0.52 s, which are nearly 30 times, 22    times, 17 times smaller than their corresponding maximum underwater fluid resultant forces, respectively. Moreover, further investigation on the effect of resistance force on squid was conducted, taking γ = 30° as an example. Compare the fluid force and moving speed at t = 0.31 s (underwater) and t = 0.6 s (in the air), the fluid force of the former in x-axis is 10 times larger than that of the latter (0.49 N to 0.05 N), while the velocity of the former is 4.4 times smaller than that of the latter (1.92 m s −1 to 8.39 m s −1 ), as shown in figures 9(b) and 10(a), the result of which might explain the tradeoff strategies of squid multimodal locomotion in water and air, i.e. the natural flying squid keeps its fins and arms folded during underwater launching phase to reduce the drag force, and spreads fins and arms in the air to generate the lift force while not inducing too much drag force. The revelation of the multi-modal locomotion trade-off mechanism of squid will also provide reference for the design and implementation of the aquatic-aerial vehicles.

Flow field velocity distribution during launching process
The velocity nephogram sequence of the squid in xOy plane under the conditions of different launch angles, i.e. γ = 30°, 45°, and 60°, is presented in figure 11. Each sequence includes three pictures which correspond to three different time instants at 0.18 s (underwater), 0.32 s (water-air interface), and 0.41 s (air), respectively. In order to show the variation of the velocity field clearly, the color bar is constrained to a velocity value range from 0 to 11 m s −1 . It should be pointed out that the flow field region of each nephogram cropped is the same in the ground coordinate O-xyz, and, the motion of the squid in the xOy plane is also presented in figure 11. The starting launching position of each set of nephogram sequences can be found in figure 5.
It can be seen that with the increase of time, the speed of the squid was also increasing. When t = 0.18 s, for all three launch angles the squid was underwater, and the body velocity of the three cases was approximated as 1.02 m s −1 and had no significant difference. The squid pierced the water surface at different angles when t is around 0.32 s, and the close-up views of water surface-crossing process when γ = 45° are shown in figure 12. At this time, the speed had some slight differences, which was 2.06, 2.03 and 2.00 m s −1 when γ = 30°, 45°, 60°, respectively. The results indicate that the speed of the squid crossing the water surface is about 2 m s −1 . For all three cases the squid was in the air when t = 0.41 s, and the velocities of different launch angles at 30°, 45°, and 60° are 3.58, 3.43, and 3.24 m s −1 , which shows that as the launch angle increases, the flying velocity decreases during aerial motion in the launching phase. In addition, as shown in figure 11, the displacement of the squid in the x direction when γ = 30° was always leading when t = 0.18, 0.32 and 0.41 s. But the squid first reached a higher altitude from the water surface when γ = 60°.
The velocity field near the body and the jet area of the squid was significantly affected due to interaction between the air, the water, the jetting and the body of the squid. The velocity magnitude of flow field far away from the squid body was almost zero due to weak influence from the squid movement. As seen from figure 11, the high speed region in the velocity field was concentrated around the jet flow area. According to figure 4(a), when t = 0.18 s, 0.32 s and 0.41 s, the jet velocity at the orifice of the squid was 15.29 m s −1 , 24.05 m s −1 , and 28.75 m s −1 , respectively. During the underwater motion part, it is obvious that the high-speed region of the jet flow (red and yellow portions, v jet (x, y , z, t) > 8 m s −1 ) is very small, only distributing near the orifice of the squid, and the velocity of the jet flow decreases rapidly along the backward jet direction. However, when the squid flies into the air, the high-speed area of the jet flow field increases significantly, and the interaction area of the jet flow on the surrounding flow field is also larger. Since the air is more easily disturbed by the motion of the squid than water, the area of the air flow field velocity around the squid body is significantly larger than that of the underwater condition. In addition, as the jet velocity increases, the area influenced by the jet flow also increases in the underwater motion part, as shown in figure 11 from v jet = 15.29 m s −1 when t = 0.18 s to v jet = 24.05 m s −1 when t = 0.32 s.

Conclusion
This paper investigated the motion of flying squid using an alternative approach, the computational fluid dynamics, which is distinguished from conventional approaches such as biological observation. The use of CFD can quantitatively investigate the squid kinematics and dynamics information. To make the results more accurate and reliable, a 3D biological prototype and initial parameters were obtained based on the investigation of natural squid. The 3D model was meshed and the boundary conditions were specified by GAMBIT, the integrated pre-process environment of Fluent. RANS (Reynolds averaged Navier-Stokes) model was used to calculate fluid motion and capture vortex feature, and VOF (volume of fluid) model was adopted as multiphase flow model. Movement of flying squid was given by the 6-DOF solver based on FVM in Fluent 16, provided with dynamic meshing technique.
By quantitative CFD simulation, several discoveries are concluded as follows: (1) The thrust-priority rather than efficiency-priority strategy works during flying squid's water-air locomotor transition. Both trailing jet and pinch-off vortex rings are formed to generate launching thrust. The formation number L ω /D ω is 5.22, indicating that the squid jet strategy during launching phase is to produce greater time-averaged thrust rather than higher propulsion efficiency; (2) The complete kinematic process of flying squid's water-air launching is reproduced, analyzed and revealed quantitatively for the first time. The maximum take-off acceleration of 26.53 m s −2 is discovered during squid's launching; (3) The maximum flying speed negatively correlates with the launch angle, indicating that a lower launch angle could result in a larger flying speed for the flying squid to escape.
To explore the amazing flying squid locomotion, further research on other motion phases of flying squid will be conducted, especially the gliding phase and diving phase. Also, the new concept of soft robotics method shows good potential application prospects [65][66][67], therefore, a soft bionic squid-like aquatic-aerial vehicle with water-air transition ability will be put into practice.