Spiraling motion of underwater gliders: Modeling, analysis, and experimental results

Ocean Engineering 60 (2013) 1–13

Contents lists available at SciVerse ScienceDirect

Ocean Engineering journal homepage: www.elsevier.com/locate/oceaneng

Spiraling motion of underwater gliders: Modeling, analysis, and experimental results Shaowei Zhang a,b, Jiancheng Yu a, Aiqun Zhang a, Fumin Zhang c,n a

State Key Laboratory of Robotics, Shenyang Institute of Automation, Chinese Academy of Sciences, No. 114, Nanta Street, Shenyang, Liaoning 110016, China Graduate School of Chinese Academy of Sciences, Beijing 100049, China c School of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta, GA 30332, USA b

a r t i c l e i n f o

abstract

Article history: Received 4 January 2012 Accepted 1 December 2012 Available online 16 January 2013

This paper presents a thorough approach characterizing the spiraling motion of underwater gliders. The dynamic model for underwater gliders, steered by a single internal movable and rotatable mass, is established. Spiraling motions are equilibria of the dynamics, for which equations are derived and then solved by a recursive algorithm with fast convergence. This theoretical method is applied to the Seawing underwater glider whose hydrodynamic coefficients are computed using computational fluid dynamics (CFD) software packages. In a recent experiment in the South China Sea, the Seawing glider produced a spiraling motion against strong ocean current, agreeing with theoretical predictions. Hence the recursive algorithm may be used to compute control input to achieve desired spiraling motion for underwater gliders in practice. & 2012 Elsevier Ltd. All rights reserved.

Keywords: Underwater glider Dynamics and control Hydrodynamic modeling

1. Introduction Underwater gliders are a class of underwater vehicles not using external active propulsion systems. Commercially available gliders such as the Slocum (Webb et al., 2001), the Spray (Sherman et al., 2001) and the Seaglider (Eriksen et al., 2001), are capable of long range missions with low energy consumption, low cost, and great endurance. They have found broad application in physical and biological oceanography. By changing its buoyancy with an internal pumping system, an underwater glider can move up and down in a water column. The forward gliding motion is generated by the hydrodynamic lift forces exerted on a pair of wings attached to the glider hull. Hence a glider does not require extra energy consumption to move horizontally, other than adjusting the buoyancy. As a result, the glider usually follows a sawtooth motion pattern in the vertical plane, and progresses along a straight line in the horizontal plane. Thus, a glider can be controlled to follow piecewise linear paths. This has enabled various mission designs in real world applications, where either a single or a group of gliders are deployed (Zhang et al., 2007; Paley et al., 2008; Leonard et al., 2010; Smith et al., 2009, 2011). The sawtooth pattern can be viewed as a steady state motion of the glider dynamics. The dynamic models of gliders steered by a linearly movable internal mass were established in Graver and Leonard (2001) and Bhatta and Leonard (2008). In Bender et al.

n

Corresponding author. Tel.: þ1 404 385 2751; fax: þ1 404 894 4641. E-mail addresses: [email protected] (S. Zhang), [email protected] (J. Yu), [email protected] (A. Zhang), [email protected] (F. Zhang). 0029-8018/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.oceaneng.2012.12.023

(2008) and Leonard and Graver (2001), the stability of the sawtooth gliding motion was analyzed, and then a model-based feedback control method was developed to stabilize the motion. In addition, a dynamic model was derived for USM’s underwater glider in Hussain et al. (2011) and Isa and Arshad (2011). Other than the sawtooth motion, an underwater glider can perform another steady state motion, which has an underwater trajectory as a spiral. Similar to the sawtooth motion, the spiraling motion is a result of hydrodynamic forces, and is also energy efficient. This spiraling motion can be used to change the direction of the glider movement underwater. It can be adjusted by a rotating internal mass, e.g., the battery pack, or a rudder attached to the tail section of the glider. Even though spiraling motions have been observed in most glider systems, the underlying mechanism is less understood than the mechanism for the linear sawtooth motion. This is because the complex dynamic model and the coupled hydrodynamic effects make it very challenging to compute the relationship between the spiraling motion and the control inputs. In Bhatta and Leonard (2008), an attempt was made to find numerical solutions for the equations that characterize the spiraling motion. Some simulation results about spiraling motion were presented in Kan et al. (2008). In Mahmoudian et al. (2010) and Mahmoudian and Woolsey (2008), an approximate analytical solution for steady spiraling motion was derived by applying perturbation theory. Our previous conference paper Zhang et al. (2011) simplified the equations of the glider dynamics, which led to a recursive algorithm that produced fast solutions for spiraling motion. In order to apply the theoretical methods to real gliders, a complete characterization of the hydrodynamic coefficients in the dynamic models is necessary. Compared to towing tank experiments,

2

S. Zhang et al. / Ocean Engineering 60 (2013) 1–13

Nomenclature

g ¼ ½PT , PT T the generalized momentum of the glider in the

e0 : ðx,y,zÞ the body frame E0 : ði,j,kÞ the inertial frame p0 : ðp1 , p2 , p3 Þ the flow frame V ¼ ½V 1 ,V 2 ,V 3 T translational velocity in the body frame X ¼ ½p,q,rT angular velocity in the body frame b ¼ ½x,y,zT glider position in the inertial frame h ¼ ½f, y, cT glider attitude in the inertial frame ms mass of the static block rs position of the static block in the body frame Is inertia matrix of the static block in the body frame mr mass of the semi-cylindrical movable block rr position of the semi-cylindrical movable block in the body frame Rr eccentric offset of the semi-cylindrical movable block g rotation angle of the semi-cylindrical movable block around the x-axis Ir ðgÞ inertia matrix of the semi-cylindrical movable block in the body frame mb mass of the adjustable net buoyancy rb position of the net buoyancy in the body frame m mass of the displaced fluid by a glider ub rate of change of the mass of the net buoyancy qsE , qbE , qrE positions of the center of mass of the static block, the net buoyancy, and the movable block in the inertial frame, respectively f ext total external force generated by the wings in the inertial frame sext total external moment generated by the wings in the inertial frame MA added mass terms in the body frame IA added inertia terms in the body frame CA added coupling terms in the body frame Tf total kinetic energy of surrounding fluid Lglider glider hull length Dglider glider hull diameter p translational momentum of the glider in the inertial frame p angular momentum of the glider in the inertial frame P translational momentum of the glider in the body frame P angular momentum of the glider in the body frame

m ¼ ½VT , XT T the generalized velocity of the glider in the

body frame

computational fluid dynamics (CFD) software offers a less expensive alternative to estimating hydrodynamic coefficients for underwater vehicles (Racine and Paterson, 2005; Boger and Dreyer, 2006). In Geisbert (2007), the added mass and added inertia terms of underwater vehicles were computed using the software package USAERO. In Morgansen et al. (2007), the drag and lift forces generated by the wings of a fin-actuated underwater vehicle were computed. In Tang et al. (2009) and Jagadeesh et al. (2009), the hydrodynamic coefficients of underwater vehicles were verified by comparing CFD results with experiments. Glider shape design had a significant impact on the performance. The hydrodynamic characteristic of a submarine launched underwater gliders was presented in Rodgers and Wharington (2010), and the glider shape was optimized to obtain a higher lift to drag ratio over a large range of attack angle. The effect of using actuated wings on a glider to enhance maneuverability were verified by experiments in a tank in Arima et al. (2008, 2009). This paper establishes a thorough approach characterizing the spiraling motion of underwater gliders: including the use of modeling, analysis, and experiments. We extend the theocratical effort

body frame the generalized inertia matrix translational and angular velocity of the static block kinetic energy of the static block translational speed of the movable block along x with respect to e0 in the body frame orx angular speed of the movable block around x with respect to e0 in the body frame V r , Xr translational and angular velocity of the movable block in the body frame Tr kinetic energy of the movable block T the total energy of the glider system Fh hydrodynamic forces in the flow frame Th hydrodynamic moments in the flow frame F hydrodynamic forces in the body frame T ¼ ½T 1 T 2 T 3  hydrodynamic moments in the body frame a the attack angle b the flip angle Re the Reynolds number of the glider Vc velocity of the buoyancy center in CFX simulation Rc turning radius of the buoyancy center in CFX simulation V inlet velocity of a fluid particle in the inlet surface in CFX simulation Rinlet turning radius of a fluid particle in the inlet surface in CFX simulation o the angular speed of fluid in CFX simulation g the gravitational constant m the fluid viscosity coefficient r the fluid density K D0 , K D coefficients of drag force D with respect to V 2 , a2 V 2 Kb coefficients of side force SF with respect to bV 2 K L0 , K a coefficients of lift force L with respect to V 2 , aV 2 K MR , K p coefficients of roll moment T DL1 with respect to bV 2 , pV 2 K M0 , K M , K q coefficients of pith moment T DL2 with respect to V 2 , aV 2 , qV 2 K MY , K r coefficients of yaw moment T DL3 with respect to bV 2 , rV 2

M Vs , Xs Ts V rx

initiated in Zhang et al. (2011) to establish an analytical model for spiraling motion of underwater gliders steered by an internal movable mass block. This theoretical model is then applied to the Seawing glider developed by the Shenyang Institute of Automation, a subdivision of the Chinese Academy of Sciences. The Seawing glider, illustrated in Fig. 1, was prototyped in 2003 and has gone through several iterations of improvements. It successfully accomplished a series of ocean tests in 2011. We then compute all the necessary hydrodynamic coefficients using CFD software, which allows us to make theoretical predictions on the spiraling motion given a control input, e.g., position of the battery pack. Furthermore, we present results from an ocean experiment performed at the South China Sea, where the glider produced a spiraling motion against strong ocean current. The experimental results agree with our theoretical predictions. To our knowledge, the experimental results on spiraling motion in the ocean have not been previously presented or analyzed in the literature. The study of the spiraling motion is relevant to new vehicle designs, for example, hybrid underwater gliders that combine the

S. Zhang et al. / Ocean Engineering 60 (2013) 1–13

3

Table 1 Glider mechanical property. List of symbol

Value

Static block ms Position of the static block rs Inertia of the static block Is Movable block mr Position of the movable block rr

54.28 kg [  0.0814 0 0.0032] m diag([0.60 15.27 15.32]) kg m2 11 kg p p rgr ; Rr ¼ 0:014 m, 2 2 0:3516 m o r rx o 0:4516 m diag([0.02 10.16 0.17]) kg m2

Inertia of the movable block I0r Net buoyancy mb Position of the red net buoyancy rb Displaced fluid mass m Glider hull diameter Dglider Glider hull length Lglider

0:5 kgr mb r 0:5 kg [0 0 0] 65.28 kg 0.22 m 1.99 m

Fig. 1. Seawing underwater glider.

long endurance of an underwater glider with the high maneuverability of a propeller or jet driven vehicle (Alvarez et al., 2008; Wang et al., 2007, 2011). The spiraling motion will offer a low energy turning behavior, serving as an alternative to fast turning behaviors that involve active propulsion. Path planning methods developed in Isern-Gonza´lez et al. (2011) and Pereira et al. (2011) may also benefit from an understanding of the spiraling motion so that a glider can be temporarily parked at a spot against a strong ocean current. The paper is organized as follows. Section 2 derives the dynamic model for underwater gliders steered by internal movable mass. Section 3 describes the procedure to compute the hydrodynamic coefficients used by the dynamic model. Section 4 presents a recursive algorithm to produce fast solutions for the nonlinear equations that establish the relationship between the spiraling motion and the control input. Section 5 presents the experimental results collected during an experiment at the South China Sea. Section 6 provides some discussions and conclusions.

2. Motion model for an underwater glider The Seawing glider in Yu et al. (2012) is driven by an internal pumping system and is steered by a movable and rotatable battery pack. The two tail wings aligned in the vertical plane are utilized to steady the vertical gliding motion. Detailed specifications of the Seawing gliders are given in Table 1. In this section, we establish the dynamic model for a Seawing glider steered by a sliding and rotating movable mass block. The glider has the following external structures (Fig. 2): the prolate ellipsoid as the rigid hull, the CTD sensor module mounted on top of the hull, two main wings located in the horizontal symmetry plane e0 2xy near the middle of the hull, the tail wings in the vertical symmetry plane e0 2yz, and the tail mast serving as the satellite communication antenna. 2.1. Coordinate frames Fig. 2 shows the three coordinate frames, inertial frame, body frame, and flow frame, established to describe the motion of an underwater glider. The body frame e0 : ðx,y,zÞ is established at the buoyancy center of the glider. The x-axis coincides with the longitudinal axis of the glider. The y-axis lies in the wing plane, pointing to the right when viewed along the direction of x. The z-axis is selected as x  y, as shown in Fig. 2. The inertial frame is

Fig. 2. A 3D illustration of the Seawing underwater glider.

described by E0 : ði,j,kÞ. In the body frame, the translational velocity and the angular velocity of the underwater glider are defined as V ¼ ½V 1 ,V 2 ,V 3 T and X ¼ ½p,q,rT , respectively. In the inertial frame, the position and the attitude of the underwater glider are described by b ¼ ½x,y,zT and h ¼ ½f, y, cT , respectively. A rotation matrix REB maps V in the body frame to the rate of change of b in the inertial frame as b_ ¼ REB V

ð1Þ

Using the simplified notation c  ¼ cosðÞ and s  ¼ sinðÞ, REB has the form of 2 3 cfsycc þ sfsc cycc sfsycccfsc 6 7 REB ¼ 4 cysc cfcc þ sfsysc sfcc þcfsysc 5 ð2Þ sy sfcy cfcy According to Fossen (2011), the relationship between X and h_ is given by 2 3 1 sf tan y cf tan y cf sf 7 h_ ¼ 6 ð3Þ 40 5X 0 sf sec y cf sec y Other than the inertial and the body frames, it is usually more convenient to compute and analyze hydrodynamics in a flow frame depicted as p0 : ðp1 , p2 , p3 Þ in Fig. 2. In order to define the flow frame, we first define the attack angle a and the slip angle b as tan a ¼

V3 V1

4

S. Zhang et al. / Ocean Engineering 60 (2013) 1–13

V2 sin b ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi V 21 þ V 32 þ V 33

The flow frame is defined relative to the body frame as follows. First we rotate the body frame around the y-axis by the angle a, as a result the z-axis is rotated to a new position which is now defined as axis p3 . We then rotate the new frame around p3 by the angle b. Thus the x becomes the flow frame axis p1 , and the y becomes the flow frame axis p2 . 2.2. Model of mechanics Fig. 3 shows a model of mechanics as a system of mass blocks: the static block with mass ms , the sliding and rotating movable block with mass mr , and the adjustable net buoyancy mb , which represents the difference between the total buoyancy and the total mass of the glider. We model the movable block as a semicylinder with eccentric offset Rr . Its center of mass is located at position r rx along x, and is rotated with an angle g around x. The vectors qrE , qbE , and qsE denote the positions of the centers of mass for the movable block, the net buoyancy, and the static block in the inertial frame, respectively. When mb ¼ 0, the displaced fluid mass of the glider is equal to the total mass of the movable block and the static block, then the glider is neutrally buoyant. If the movable mass block mr stays at the initial position, then the glider floats horizontally without motion relative to the surrounding fluid. When mb 4 0, the movable mass block mr slides forward to generate a negative pitch angle to drive the glider to dive down; and when mb o0, the movable mass block mr slides backward to generate a positive pitch angle to drive the glider to climb up. The sawtooth motion is shown in Fig. 4(a).

j

i x

qrE

Rr

mr

y

q sE

rr

rs

ms

e0

E0

b qbE

rb

k

mb

z

Fig. 3. Mass distribution in the body frame and the inertial frame.

Other than the linear motion, the movable mass block mr can rotate around the x-axis to achieve turning motion shown in Fig. 4(b): when mr rotates for an angle g, the glider hull generates a roll angle in the opposite direction of g, creating a rotational offset for the static block. As a result, the lift force generated by the pair of wings attached to the hull yields a force component along j in the inertial frame serving as the centripetal force for the turning motion. 2.3. Mathematical model of dynamics We define p as the translational momentum and p as the angular momentum of the glider in the inertial frame. We define P and P as the translational momentum and angular momentum relative to the glider body frame. The related translational velocity and angular velocity of the movable mass block mr relative to the body frame, caused by the driving motor, are assumed to satisfy V rx ¼ 0 and orx ¼ 0 for all time. These assumptions are reasonable in practice since the movable mass block moves slowly. As in Leonard and Graver (2001), according to Newton’s laws, we know that p_ ¼ f ext þ ðms þmr þ mb mÞgk

p_ ¼ qrE  mr gk þqbE  mb gk þ qsE  ms gk þ sext

ð4Þ

where k is the unit vector pointing to the direction of gravity, f ext is the external force generated by the wings attached to the glider, and sext is the external moment. We map the generalized momentum ½pT pT T from the inertial frame to the generalized momentum g ¼ ½PT PT T in the body frame as p ¼ REB P

p ¼ REB P þ b  p

ð5Þ

Then we combine the generalized translational velocity V and the angular velocity X of the glider in the body frame as m ¼ ½VT XT T . The relationship between g and m is

g ¼ Mm

ð6Þ

where M is the generalized inertia matrix of the glider system. By differentiating Eq. (6) with respect to time, we get the relationship among the generalized force, velocity and acceleration as _ m þ Mm_ g_ ¼ M

ð7Þ

Then we need to derive the total kinetic energy of the glider _ m , and m_ . system to find equations for M, M, We now derive the total kinetic energy of the glider system. Since the center of the movable block has a constant offset from the center of the body frame e0 , we express the generalized translational velocity and angular velocity of the static block in the body frame as Vs ¼ Vr^ s X Xs ¼ X

Fig. 4. Mechanisms for the linear sawtooth motion and the spiraling motion. (a) Generate linear sawtooth motion and (b) Generate turning spiraling motion.

S. Zhang et al. / Ocean Engineering 60 (2013) 1–13

We define the operator ^ as the skew-symmetric matrix constructed from a vector. For the three-dimensional column ^ ¼ x  y. Then the kinetic energy of the vectors x and y, we have xy static block is T s ¼ 12 ms JVs J2 þ 12Xs  Is Xs ¼ 12m T Ms m

ð8Þ

The net buoyancy is controlled by a pump system that changes the volume of a bladder by pumping oil to and from an oil chamber. The oil chamber is located at the center of buoyancy of the glider. Hence the net buoyancy does not contribute to the inertial terms or to the kinetic energy of glider dynamic systems. The velocities and kinetic energy of mr vary as r rx and g change. The position of the movable block rr is rr ¼ r rx x þRr ðcosðg þ p=2Þy þ sinðg þ p=2ÞzÞ

ð9Þ

We can express the translational velocity and the angular velocity of the movable block in the body frame as ð10Þ

So the kinetic energy of the movable block is expressed as 1 mr JVr J2 þ 12XTr  Ir ðgÞXr 2 1 ¼ m T Mr m 2

Tr ¼

ð11Þ

The inertia matrix of the movable block is a function of g as Ir ðgÞ ¼ RTx ðgÞI0r Rx ðgÞ 1

6 Rx ¼ 4 0 0

0

ð12Þ

0

cos g

^ PÞ p_ ¼ REB ðP_ þ X ^ Þ þREB V  p þ b  p_ _ þ XP p_ ¼ REB ðP

ð17Þ

Substituting Eq. (17) into Eq. (4) will give P_ ¼ P  X þ mb gðRTEB kÞ þRTEB f ext

sin g

cos g

1 T f ¼ m T Mf m 2

ð13Þ

CA

CTA

IA

ð18Þ

Then substituting Eq. (18) into Eq. (16), we get the dynamic model as in Eq. (19), where F ¼ REB f ext , and T ¼ REB sext are the hydrodynamic forces and moments in the body frame: ( " # PX _ mþ m_ ¼ M1 M P  X þP  V " #  ) mb gðRTEB kÞ F þ þ T ðmr rr þms rs þmb rb Þg  ðRTEB kÞ ð19Þ

2.4. Hydrodynamic forces

When a glider accelerates in the flow, the surrounding fluid is also accelerated. This results in an apparent increase in mass and inertia. These forms of acceleration-dependent hydrodynamic effects are represented by the added mass, added inertia, and cross terms. Following Leonard and Graver (2001), we express the kinetic energy generated by the added mass MA , the added inertia matrix IA and the cross term CA due to the fluid effect as

MA

ðRTEB kÞ þRTEB sext

sin g 7 5

and represents the principal inertia matrix of mr computed in the stationary state with g ¼ 0. Eqs. (9) and (12) model the effect of the movable/rotatable mass in both the sawtooth motion and the spiraling motion. When the movable mass rotates, the roll angle changes. From Eq. (12), we know that the inertia of Ir ðgÞ changes as g varies.

Mf ¼

Differentiating Eq. (5) on both sides and applying the kinematic matrix in Eqs. (1) and (3), we get the relationship among g, m , m_ , and g_ as

_ b ¼ ub m

3

I0r

where "

where Mt ¼ ðmr þms ÞI3 þ MA , Ct ¼ CA ms r^ s mr r^ r , and It ¼ Is þ Ir ðgÞ þIA mr r^ r r^ r ms r^ s r^ s . Applying the Legendre transform to derive the equations for linear momentum and angular momentum gives Eq. (6). Then we take the time derivative of g and obtain " # V_ _ mÞ m_ ¼ _ ¼ M1 ðg_ M ð16Þ X

_ ¼ P  X þ P  V þðms rs þ mr rr þ m r Þg P b b

Vr ¼ Vr^ r X Xr ¼ X

where 2

5

# ð14Þ

In the flow frame, the hydrodynamic force Fh ¼ ½D SF LT and hydrodynamic moment Th ¼ ½T DL1 T DL2 T DL3 T are usually expressed as D ¼ ðK D0 þK D a2 ÞV 2 SF ¼ K b bV 2 L ¼ ðK L0 þ K a aÞV 2 T DL1 ¼ ðK MR b þK p pÞV 2 T DL2 ¼ ðK M0 þK M a þ K q qÞV 2 T DL3 ¼ ðK MY b þK r rÞV 2 A rotation matrix RBC is defined as 2 3 cos a cos b cos a sin b sin a 6 sin b cos b 0 7 RBC ¼ 4 5 cosa sin a cos b sin a sin b

ð20Þ

ð21Þ

RBC maps the hydrodynamic force and moment from the flow frame to the body frame as F ¼ RBC Fh ,

T ¼ RBC Th

ð22Þ

Finally, we get the total energy of the glider/fluid system as T ¼ T r þT s þ T f 1 ¼ m T Mm 2 The generalized inertia matrix is then " # Mt Ct M¼ CTt It

2.5. Added mass and added inertia terms ð15Þ The added terms in Eq. (14) are computed by applying strip theory. The glider hull is approximated as a slender hull. We consider it as a prolated ellipsoid, and apply the strip theory (Fossen, 2011) to compute the added mass and inertia terms numerically. Since the glider has two symmetry planes as e0 2xz

6

S. Zhang et al. / Ocean Engineering 60 (2013) 1–13

3.1. Turbulence models The Reynolds number Re is the ratio between the inertial effects and the viscous effect in the fluid, and the solution of viscous flow is affected by the Reynolds number. The Reynolds number is computed as

rVLglider ð23Þ m 3 where r ¼ 1025 kg=m is the water density, V is the glider velocity, Lglider is the total length (not including the tail), and m is fluid viscosity

Re ¼ Fig. 5. Parameters for the external geometry of the glider.

and e0 2xy, the added mass matrix is approximately diagonal 2 3 X u_ 0 0 6 0 Y 0 7 MA ¼ 4 5 v_ 0 0 Z w_ where X u_ , Y v_ and Z w_ are the acceleration fluid mass terms around the glider generated from the force along the x, y, z axes, respectively. Also, due to the symmetry of the glider, the added inertia matrix IA has a diagonal form as 2 3 K p_ 0 0 6 0 M 0 7 IA ¼ 4 5 q_ 0 0 Nr_ where K p_ , M q_ , Nr_ are the acceleration fluid inertia terms generated from the moments around x, y, z axes, respectively. The matrix CA of cross terms contains the coriolis and centripetal terms. Considering the symmetric properties of the glider, we utilize the strip theory to estimate the following four parameters in matrix CA : N v_ , M w_ , Y r_ , Z q_ , where M w_ and Nv_ are the pitch moment with respect to the acceleration in the z direction, and yaw moment with respect to acceleration in the y direction, respectively, and Y r_ and Z q_ are the side force due to yaw acceleration and the lift force due to pitch acceleration, respectively. So CA has the form 2

0 60 CA ¼ 4 0

0 0 Nv_

3

2 0 7 M w_ 5 ¼ 6 40 0 0 0

3T

0

0

0

Y r_ 7 5

Z q_

0

The principle of the strip theory includes the following steps. First, we partition the glider into a number of strips. Then we compute the two-dimensional hydrodynamic derivatives of the added terms for each strip. Finally we integrate the two-dimensional hydrodynamic derivatives to get the three-dimensional added terms. Detailed dimension and property of the glider are shown in Fig. 5 and Table 1. Using strip theory, we have numerically computed the added hydrodynamic terms as X u_ ¼ 1:48 kg, Y v_ ¼ 49:58 kg, Z w_ ¼ 65:92 kg, K p_ ¼ 0:53 kg m2 , Mq_ ¼ 7:88 kg m2 , N r_ ¼ 10:18 kg m2 , N v_ ¼ 2:57 kg m, Z q_ ¼ 3:61 kg m.

coefficient. The glider velocity ranges from 0.257 m/s (0.5 knot) to 0.514 m/s (1.0 knot). Hence the Reynolds number ranges from 4:32  105 ðV ¼ 0:257 m=sÞ to 8:65  105 ðV ¼ 0:514 m=sÞ. When 1  105 oRe o 1  106 , the low Reynolds turbulence models give good estimation for underwater vehicle hydrodynamics (Jagadeesh et al., 2009). Two popular low Reynolds turbulence models (Jagadeesh et al., 2009; Yoon et al., 2007; Jagadeesh and Murali, 2005) are used in our work: (1) the standard k2e turbulence model and (2) the standard ko turbulence model. In this paper, we use the k2e model to simulate the surrounding flow when the glider is performing linear motion, and use the ko turbulence model to simulate the two-dimensional turbulent flow when the glider is performing spiraling motion. 3.2. Simulation setup For the k2e simulation, the fluid volume is selected as 5Lglider  18Dglider  18Dglider . The glider buoyancy center is 9Dglider from the ceiling and the floor of the fluid volume, and 9Dglider from the left side and the right side of the fluid volume, respectively. The distance between the glider buoyancy center and the inlet of the fluid volume is 1:5Lglider and the distance between the glider buoyancy center and the outlet of the fluid volume is 3:5Lglider . The fluid volume in the ko simulation is ring-shaped to simulate the turning motion. The dimension of the volume is 5Lglider  18Dglider  18Dglider . The total arc length of the centerline of the ring is 5Lglider , the arc length along the centerline from the inlet to the buoyancy center is 1:5Lglider . Other parameters of the fluid volume are the same as in the k2e simulation. We use Gridgen as the pre-processing tool to create a mesh for the glider. The fluid volume is meshed as two unstructured grid fields: a coarse mesh field with grid size 200 mm and a refined mesh field with grid size 80 mm. The refined mesh field is confined to the fluid volume close to the glider hull to enhance the mesh topology performance near the glider surface. The coarse mesh is utilized for the fluid field further away from the glider hull to expedite the computing process. After the fluid volume is defined, we specify the fluid boundary conditions for the fluid volumes. The glider surface boundary is set as a no-slip wall. In the k2e simulation, the upstream inlet is set as a velocity-inlet boundary with V inlet ¼ 0:257 m=s, which is

3. Hydrodynamic modeling When fluid flows past underwater vehicles, viscous turbulent flow happens at the fluid vehicle boundary layer. The Reynolds Averaged Navier–Stokes (RANS) equations method is a popular method for hydrodynamic simulations of underwater vehicles. In this study, we use a computational fluid dynamics (CFD) software package, the CFX, that implements RANS to compute the hydrodynamic coefficients. Proper setup of CFD simulation requires knowledge of the Reynolds number, proper selection of a turbulence model, and proper settings for grid generation and boundary conditions.

2

,

1 3

Fig. 6. A section of the ring shaped water volume.

S. Zhang et al. / Ocean Engineering 60 (2013) 1–13

perpendicular to the inlet flow surface. The downstream outlet satisfies the zero static pressure condition, and the other four boundaries of the fluid volume are set as free-slip walls. The configuration for the ko simulation is shown in Fig. 6. The linear speed at the glider buoyancy center is fixed at V c ¼ 0:257 m=s. The angular speed of the fluid volume o is specified by changing the curvature of the ring. Let the radius Rc be the distance between the curvature center of the ring and the glider buoyancy center, then the angular speed of the fluid volume is

o¼

Vc Rc

ð24Þ

The speed of each point on the inlet surface varies as a linear function of the radius of curvature Rinlet . However, the angular speed o of the fluid volume is kept constant. The speed distribution at the inlet surface can be generated by the angular speed and the distance between any point on the inlet surface and the curvature center of the ring. Let the position of the curvature center of the ring in the flow frame be Rc ¼ ½Rx Ry Rz T , and let the position of a point on the inlet surface in the flow frame be r0 ¼ ½x0 y0 z0 T . The radius of curvature for a point on the inlet qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi surface Rinlet is obtained as Rinlet ¼ ðx0 Rx Þ2 þðy0 Ry Þ2 þ ðz0 Rz Þ2 . Therefore, the speed of flow at that point is V inlet ¼ oRinlet . Other boundary conditions in the ko simulation are identical to those used in the k2e simulation. In the k2e simulation, the angular speed of the fluid is zero. The glider buoyancy center is fixed at the origin of the flow frame. We rotate the fluid volume to achieve different values for the attack angle a and the side slip angle b, then the CFX software computes the hydrodynamic force in each case. In the ko simulation, we rotate the ring-shaped fluid volume with attack angle a and flip angle b, and change the radius of the ring. The CFX software computes the hydrodynamic moments in each case, as shown in Fig. 7. Table 2 lists the combinations of the values of a, b, and Rc for all CFD simulations performed.

3.3. Data analysis and curve fitting From the boundary conditions V c , V inlet , a, b, we can determine the glider velocities V 1 , V 2 , V 3 , p, q, r. Then from the measured forces and moments, the least mean square method is used to find out the hydrodynamic parameters based on Eq. (20). The glider velocity V ¼ ½V 1 ,V 2 ,V 3 T is computed by converting the inlet velocity from the flow frame to the body frame as 2 3 3 2 3 2 V inlet cos a cos b V1 V inlet 6V 7 7 6 7 6 V ð25Þ 4 2 5 ¼ RBC  4 0 5 ¼ 4 5 inlet sin b V inlet cos b sin a V3 0

7

When the glider glides with this angular velocity aligned with p3 in the flow frame, we map the angular velocity to the body frame to compute p, q, r as 2 3 2 3 2 3 p 0 o sin a 6 7 6 7 6 7 0 ð26Þ 5 4 q 5 ¼ RBC 4 0 5 ¼ 4 r

o

o cos a

On the other hand, when the glider glides with its angular velocity aligned with p2 in the flow frame, we get the angular velocity of the glider in the body frame as 3 2 3 2 2 3 o cos a sin b 0 p 7 6 7 6 6 7 o cos b ð27Þ 5 4 q 5 ¼ RBC 4 o 5 ¼ 4 o sin a sin b 0 r Fig. 7(a) gives a snapshot of the ko simulation when the glider turns left around p3 -axis, and Fig. 7(b) gives a snapshot when the glider turns down around p2 -axis. Fig. 8 illustrates the curve fitting procedure for estimating the hydrodynamic coefficients. Data from the k2e simulations is used to estimate the hydrodynamic force coefficients based on Eq. (20). The drag force D and lift force L are functions of a only. Hence K D0 , K D , K L0 , and K a can be determined from the data collected by letting b ¼ 0 and a vary from  121 to 121, as shown in Fig. 8(b) and (c). The side force SF is a function of the slip angle b. Hence K b is determined from the data collected by letting a ¼ 0 and b vary from  121 to 121, as shown in Fig. 8(f). The moment T DL2 around the y-axis in the k2e simulation is a function of a when q¼0. Hence we can determine K M and K M0 from the data collected by letting a vary from 121 to 121, as shown in Fig. 8(d). Similarly, data from the ko simulation is used to find the coefficients for the hydrodynamic moments based on Eq. (20). We know that the moment T DL3 is a function of b and r. Hence K MY and K r can be determined from data collected when glider turns around the p3 -axis, as shown in Fig. 8(g). The roll moment T DL1 in Eq. (20) is a function of b and p. Hence K MR and K p are determined from data collected when glider turns around the p2 -axis, as shown in Fig. 8(h) and (i). T DL2 is a function of a and q. Since K M and K M0 are already determined from the k2e simulations, the coefficient K q can be determined from the data collected when Table 2 Combinations of parameter values for CFD simulation. Simulation

V inlet , V c , a, b

k2e

V inlet ¼ 0:257 m=s, (b ¼ 0, a ¼ 12,11, . . . 12), ða ¼ 0, b ¼ 12,11, . . . 12) V c ¼ 0:257 m=s, Rc ¼ 5,10,20,30,40 m; ða ¼ 0, b ¼ 3,0,3,6,9,12; b ¼ 0, a ¼ 3,3,6,9Þ

k2o in horizontal plane and vertical plane, separately

Fig. 7. Fluid speed distribution in the ko simulations. (a) Turning around p3 and (b) Turning around p2 .

8

S. Zhang et al. / Ocean Engineering 60 (2013) 1–13

Fig. 8. CFD simulation results. (a) Ratio of lift and drag. (b) Drag force D vs a. (c) Lift force L vs a. (d) Pitch moment T DL2 vs a. (e) Pitch moment T DL2 vs q. (f) Side force SF vs b. (g) Steer moment T DL3 vs r. (h) Roll moment T DL1 vs b and p (descending and. (i) Roll moment T DL1 vs b and p (ascending).

In order to enhance glider performance, we need to determine the angle of attack that generates the maximum lift–drag ratio. From Fig. 8(a), we see that the optimal angle of attack is approximately 771.

Table 3 Hydrodynamic coefficients of the Seawing glider. List of symbol

Value

Added mass terms MA Added inertia terms IA Added coupling terms CA Coefficients of drag force D

diag([1.48 49.58 65.92]) kg diag([0.53 7.88 10.18]) kg m2 N v_ ¼ 2:57 kg m, Z q_ ¼ 3:61 kg m K D0 ¼ 7:19 kg=m, 2

Coefficients of side force SF Coefficients of lift force L Coefficients of T DL1 Coefficients of T DL2

Coefficients of T DL3

K D ¼ 386:29 kg=m=rad K b ¼ 115:65 kg=m=rad K L0 ¼ 0:36 kg=m, K a ¼ 440:99 kg=m=rad K MR ¼ 58:27 kg=rad, K p ¼ 19:83 kg s=rad K M0 ¼ 0:28 kg, K M ¼ 65:84 kg=rad, K q ¼ 205:64 kg s=rad K MY ¼ 34:10 kg=rad, K r ¼ 389:30 kg s=rad

2

2

glider turns around the p2 -axis, as shown in Fig. 8(e). We list all the hydrodynamic coefficients, added mass and added inertia terms in Table 3.

4. Characteristics and analysis of spiraling motion Since we have determined all the coefficients in the dynamic model of the underwater glider, we can fully characterize its steady state motion. The steady state motion can be characterized as a steady linear motion and a steady spiraling motion. When the movable block slides from its initial position along x, the gravitational center shifts away from the buoyancy center e0 along x, and the moment caused by this offset drives glider up and down in the vertical plane. When the movable block rotates around x, the gravitational center deviates from buoyancy center e0 along y. The moment generated by the offset of the movable block causes the glider to roll around x until it is balanced by T1 and the static block. Since glider wing-forces are no longer collinear with z, the lift force generates a vertical component to balance glider net

S. Zhang et al. / Ocean Engineering 60 (2013) 1–13

buoyancy and a horizonal component as the centripetal force, causing the glider to enter a spiraling motion. With the assumption that the movable block is fixed at rr and _ b ¼ 0, that there is a constant net buoyancy, which implies that m the glider equilibria equations in three dimensions are 0¼P

X þ mb gðRTEB kÞ þF

0 ¼ P  X þ P  V þðmr rr þ ms rs Þg  ðRTEB kÞ þ T

ð28Þ ð29Þ

And P and P are simplified as P ¼ Mt Vðms r^ s þmr r^ r ÞX

ð30Þ

P ¼ ðms r^ s þmr r^ r ÞV þ It X

ð31Þ

The steady state linear motion has been investigated in Leonard and Graver (2001). In this paper, we analyze the steady state spiraling motion where yaw angle changes at a constant rate while the pitch and roll angles are constant. This implies that RTEB k is constant. By taking the time derivative of RTEB k, we get

X  ðRTEB kÞ ¼ 0

ð32Þ

From Eq. (32), we know that the glider moves with constant speed along a circular helix which is aligned with gravity, hence the angular velocity is X ¼ RTEB ko3 . The spiraling motion can be projected into a rotational motion around k and a linear motion qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi along k. We project the total velocity V ¼ V 21 þV 22 þ V 23 into the turning direction and the vertical direction as V cosðyaÞ ¼ o3 R

ð33Þ

V sinðyaÞ ¼ V vertical

ð34Þ

From Eq. (34), we get the vertical velocity in the inertial frame when the glider glides in 3D spiraling motion. Using the fact that X ¼ RTEB ko3 , we expand the force equation along x and y in Eq. (28), and moment equation around y in Eq. (29). This gives the following equations: 0 ¼ V o3 ðmt2 sin b cos f cosymt3 sin a cos b sin f cos yÞ   sin 2y cosðf þ gÞ mb g sin y þ mr o23 r rx cos2 y þ Rr 2   sin 2 y cos fr s3 þ ms o23 cos2 yr s1 þ 2 D cos a cos bSF cos a sin b þ L sin a

þcos gðmr gRr sin y þmr Rr r rx o23 ðcos2 f cos2 ysin2 yÞÞ

ð37Þ

where mt1 , mt2 , mt3 are the diagonal terms of Mt , and mA1 , mA2 , mA3 , IA1 , IA2 , IA3 are the diagonal terms of MA and IA , respectively. The steady spiraling motion can be determined by 10 parameters. V, a, and b describe the velocity of the glider. o3 , f, and y describe the angular velocity. r rx and g describe the position of the movable mass block. The other two parameters are the net buoyancy mb and the turning radius R. By combing Eqs. (28), (29) and (33), we obtain seven equations of 3D spiraling motion. Since these equations may not be independent, we need to specify at least three parameters out of 10 to solve for the other seven parameters. Solving seven coupled nonlinear equations is computationally challenging for the low power embedded computer on-board the glider. Hence we derive a recursive algorithm to obtain fast solutions. We consider the situation where V, a, and b are known, from which we know the expressions of T, F, and RBC . Then we solve for the other seven parameters. In order to simplify Eqs. (28) and (29), we take the inner products with respect to X on both sides to get 0¼

mb g

o3

X  X þF  X

0 ¼ ðP  V þ TÞ  X

ð38Þ ð39Þ

For the Seawing glider, the vector ms rs þmr rr is approximately located in the plane formed by the vector rr and RTEB k. Therefore, we take the inner products on both sides of Eq. (29) with respect to rr , and we get the following equation: 0 ¼ ms rs g  ðRTEB kÞ  rr þ ðIt X  XÞ  rr þ ðMt V  VÞ  rr þ T  rr

ð40Þ

From Eq. (38), we can solve for mb as mb ¼

F  ðRTEB kÞ g

o3 ¼  ð35Þ

ð41Þ

ð36Þ

sin 2g sin f sin 2y 4 sin 2a sin 2f cos2 y sin g 2 2 mr Rr r rx o23 þ ðmA3 mA1 ÞV cos b 2 2

mr V o3 r rx ðcos a cos b sin f cos y þ sin y sin bÞ ! Isx Isz þ Irx þ IA1 IA3 sin 2y þ o23 cos f 2 2 Iry sin g þ ðmr Rr Irz Þcos2 gmr r 2rx 2 mr gr rx cos f cos y þ T DL1 sin b þ T DL2 cos b   sin 2y þr s3 r s1 ðcos2 f cos2 ysin2 yÞ þ ms o23 ðr 2s3 r 2s1 Þcos f 2

ðF þmb gRTEB kÞ  V T  ðRTEB kÞ

ð42Þ

Note that T is also a function of o3 . Hence Eq. (42) is a quadratic equation for o3 , which can be solved analytically if y, f are known. Therefore, supposing solutions of y, f are found, we can find mb and o3 from Eqs. (41) and (42). After that, we can solve for the turning radius from Eq. (33) as R¼

0 ¼ ðIry Irz mr R2r Þo23

þ ms o3 V sin bðcos f cos yr s3 sin yr s1 Þ ms gðr s3 sin y þr s1 cos f cos yÞ

ms V o3 sin f cos y cos bðr s3 sin ar s1 cos aÞ þmr VRr o3 cos g cos yðsin b cos fsin f sin a cos bÞ

Using Eqs. (39) and (28), we can get the angular velocity o3 as

0 ¼ V o3 ðmt3 sin a cos b sin ymt1 cos a cos b cos f cos yÞ þ mb g sin f cos yD sin b þSF cos b   sin f sin 2y sin 2f cos2 y r s1  r s3 þ ms o23 2 2   r sin 2 y sin f R sin 2f cos2 y cos g rx r  þ mr o23 2 2 þ mr o23 ðRr sin gðcos2 f cos2 y þ sin2 yÞÞ

9

V cosðyaÞ

o3

ð43Þ

With the assumption that we know V, a, and b, we can solve for mb , o3 and R as functions of y and f. This process has reduced the unknown parameters from seven to four. We still need to solve for y, f, g, and r rx . When o3 is small, Eqs. (35) and (36) can be used to solve for y and f. Then Eqs. (37) and (40) can be used to solve for g and r rx . When o3 is not small, Eqs. (35)–(37), (40)–(42) have to be solved altogether. We define a recursive vector as D ¼ ½y,R, f,mb , g, o3 ,r rx , and establish a recursive equation as

Dk ¼ FðDk1 Þ

ð44Þ

to solve for the recursive vector D. Equations for mb , o3 , and R have been derived as Eqs. (41)–(43). We now derive the equations for y, f, g, and r rx .

10

S. Zhang et al. / Ocean Engineering 60 (2013) 1–13

From Eq. (35), we solve for y as f

y ffi ly y ¼ arcsin qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

ð45Þ

2 2 f y1 þf y2

where f y1 ¼ V o3 ðmt2 sin b cos fmt3 sin a cos b sin fÞ f y2 ¼ mb g f y ¼ D cos a cos b þ SF cos a sin bL sin a   sin 2y cosðf þ gÞ mr o23 r rx cos2 y þ Rr 2   sin 2 y cos fr s3 ms o23 cos2 yr s1 þ 2 f y1 f y2 ffi sin ly ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi , cos ly ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 2 f y1 þ f y2 f y1 þf y2

When o23 is small, terms that contains r rx and g disappear from Eqs. (45) and (46). Therefore after we replace mb by Eq. (41), Eqs. (42), (45) and (46) can be solved recursively for y and f. We next derive the recursive equations from r rx and g. From Eq. (37), we express g as f

g ffi lg g ¼ arcsin qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 2 where sin lg ¼ f g1 = f g1 þ f g2 , cos lg ¼ f g2 = f g1 þ f g2 , and f g1 , f g2 , and f g are defined in Eq. (48). From Eq. (40), we get r rx as Eq. (49). Therefore, we have shown that the solution of steady spiraling motion can be derived with a recursive algorithm given by Eq. (44): f g1 ¼ mr VRr o3 cos yðsin f sin a cos bsin b cos fÞ þ mr gRr sin y

Similarly, from Eq. (36), we express f as

mr Rr r rx o23 ðcos2 f cos2 ysin2 yÞ

ff ffi lf f ¼ arcsin qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 f f1 þ f f2

ð46Þ

f g2 ¼ mr Rr r rx o23 cos f cos2 y sin f

where f g ¼ ðIry Irz mr R2r Þo23

sin 2y sin 2y r s1 mr o23 r rx f f1 ¼ mb g cos yms o 2 2 f f2 ¼ V o3 mt1 cos a cos b cos y 2 3

sin 2g sin f sin 2y 4

þðmA3 mA1 ÞV 2 cos2 b

f f ¼ SF cos bD sin b sin 2f cos2 y r s3 V o3 mt 3 sin a cos b sin yms o23 2   2 Rr sin 2f cos y cos g 2 Rr sin gðcos2 f cos2 y þ sin yÞ þ mr o23  2 f f2 sin lf ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi , 2 2 f f1 þ f f2

ð47Þ

2

f g1 þf g1

mr r rx o3 Vðcos a cos b sin f cos y þ sin y sin bÞ

þ

f f1 cos lf ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 f f1 þ f f2

sin 2a 2

Isx Isz þ Irx þ IA1 IA3 Iry sin2 g þ ðmr R2r Irz Þ cos2 gmr r 2rx

!

o23 cos f

sin 2y 2

  sin 2y þ r s3 r s1 ðcos2 f cos2 ysin2 yÞ þms o23 ðr 2s3 r 2s1 Þcos f 2

Table 4 Signs for initial values for the recursive vector. Gliding motion

V

a

b

mb

r rx

g

o3

y

f

DL DR AR AL

þ þ þ þ

þ þ  

þ  þ 

þ þ  

þ þ  

 þ þ 

 þ þ 

þ   þ

  þ þ

þ ms o3 V sin bðcos f cos yr s3 sin yr s1 Þ þ T DL1 sin b þ T DL2 cos b

ms V o3 sin f cos y cos bðr s3 sin ar s1 cos aÞ ms gðr s3 sin y þ r s1 cos f cos yÞmr gr rx cos f cos y

ð48Þ

Fig. 9. Left: Simulated trajectory of a glider spiraling downwards. Right: The nonlinear relationship between the turning radius R and the angle of the battery pack g. (a) Spiraling motion and (b) Relationship between R and g.

S. Zhang et al. / Ocean Engineering 60 (2013) 1–13

 1 V 2 ðmt1 mt3 ÞRr sin22a cos2 b sin g þ ðmt1 mt2 ÞRr sin22b cos a cos g C 1 0 C sin g cos fðIsx þ IA1 Isz IA3 Þ C C B þsin f cos gðI þI I I Þ C C sx sy A1 A2 C B C C B C B þIrx sinðg þ fÞIrz sin 2g cosðgfÞ C C C B C 2 sin 2y B C C o3 Rr 2 B þIry cos 2g sinðgfÞ C C C B C B þ ms Rr sin go2 sin 2y cos fðr 2 r2 Þ C C 3 s1 s3 C 2 B C  A @ C 2 2 2 C þr s3 r s1 ðsin ycos f cos yÞ C C C C T DL2 Rr sin g þT DL3 Rr cos g C   2 C yr s3 sin 2y sin fr s1 C ms Rr o23 cos grs1 sin 2f cos  2 2 A þms gððr s3 sin y þ rs1 cos f cos yÞRr sin g þr s1 Rr cos g sin f cos yÞ   1 0 o23 cos2 y ðIsz þIA3 Isy IA2 Þ sin22f þ ðIrz Iry Þ sin 2ð2fgÞ C B C B C B T DL1 þ ðmt3 mt2 ÞV 2 sin a sin22b C B   A @ sin 2f 2 sin 2y 2 þms o3 r s3 2 sin fr s1  2 cos yrs3 þms gr s3 sin f cos y

0

1032

B B B B B B B B B B B B B B B B B B B B B B B @ rrx ¼

1030 1028 1026 1024 1022 1020 0

ð49Þ

200

400 600 Depth(m)

800

1000

Fig. 10. Measured sea water density versus depth.

0.8 V V

0.6 Max Velocity(m/s)

In order to get a rapid convergence in this recursive algorithm to the desired spiraling motion, proper initialization of the recursive algorithm is necessary. As we know, the three-dimensional spiraling motion consists of four gliding situations: (DL) glider descends and turns left; (DR) glider descends and turns right; (AL) glider ascends and turns left; (AR) glider ascends and turns right. The signs of the parameters in the recursive vector are different with respect to different gliding state. This is shown in Table 4. We demonstrate a glider 3D spiraling motion equilibrium simulated by MATLAB in Fig. 9. The parameters are: g ¼ 451, r rx ¼ 0:4216 m, mb ¼ 0:3 kg, V¼0.490 m/s, a ¼ 1:2671, b ¼ 1:2831, o3 ¼ 0:0039 rad=s, f ¼ 13:7031, y ¼ 35:6411. The radius of the spiraling motion is R  100:83. We have also computed the relationship between the battery angle g and the glider turning radius R as shown in Fig. 9. Other groups of parameters are also selected for simulation. The results given by the recursive algorithm are within 5% of the simulation results based on the full glider dynamics.

11

V

East North Vertical

0.4

0.2

0

−0.2 0

200

400 depth(m)

600

800

5. Experimental results In July 2011, the Seawing glider was tested in Western Pacific Ocean (longitude: 13011:260 , latitude:3111:440 ) near Mindanao, Philippines. A series of gliding experiments have been performed. We present results collected for a downward spiraling motion during a 1.5 h time window from 12:28 BJT to 14:10 BJT on July 20th, 2011. The nominal net buoyancy was set to neutral, e.g., mb ¼ 0, and the nominal position of the weight block was at r rx ¼ 0:4016 m and g ¼ 0. When the glider dove, the net buoyancy was set to be mb ¼ 0:5 kg, and the movable block was controlled near the nominal position with an offset 9Dr rx 9  0:02 m. Our main goal was to verify whether the theoretical model developed for the Seawing glider agreed with the experimental observations. If the model was consistent with real experiments, then the recursive algorithm we developed can be used to control the control input for the Seawing glider in practice. During the experiments, the depth of the glider was inferred from the CTD (the Conductivity-Temperature-Depth profiler), and the attitude angles were measured by a TCM3 digital compass. The sampling periods of TCM3 and CTD were approximately 6 s. Seawater density depends on temperature, depth, and salinity. The density at different depths was estimated from the CTD data, as shown in Fig. 10. The rate of change for density became smaller as the depth increases. Seawater density increases by about 0.7% from the surface to 800 m depth. The glider reached the maximum gliding depth, approximately 800 m, before it became

Fig. 11. Measured maximum velocity of ocean current versus depth. V East is the strength of the latitudinal current with positive direction eastward. V North is the strength of the longitudinal current with positive direction northward. V Vertical is the strength of the up and down current with positive direction aligned with gravity.

neutrally buoyant. It then ascended to the surface by increasing the buoyancy. At the test site, there exists an ocean current with velocity Vcurrent , an advanced doppler current profiler (ADCP) (not installed on the glider) was utilized to measure the ocean current Vcurrent with an 8 m depth interval from the ocean surface to the depth of 800 m. Fig. 11 plots the maximum current measured in 24 h during the experiment day. The current measured by the ADCP at the depth where the glider has penetrated can then be used to calculate the velocity of the glider relative to water as Vr ¼ VR1 EB V current

ð50Þ

In the dynamic model (Eq. (19)) for the glider, the term V in Eq. (19) should now be replaced by Vr to generate the prediction of the motion of the glider based on computer simulation. The experimental data collected and the simulation results are compared in Fig. 12. The simulation results were generated by a MATLAB simulation with the same parameters and added ocean currents as measured in the experiments. Fig. 12(e) and (f) plots the r rx and g for the movable block. The position of the movable block r rx was adjusted to control the pitch angle y. The angle

12

S. Zhang et al. / Ocean Engineering 60 (2013) 1–13 Pitch Angle θ(°)

Depth 0

Roll Angle φ (°)

0

100

Experiment Simulation

−5

Experiment Simulation

40

−10

30

200

−15 300

−20

400

−25

20 φ (°)

Glider Depth(m)

Experiment Simulation

10

−30

500

0

−35 600

−40

−10

−45

700 1000

2000

3000

4000

5000

6000

500

Yaw Angle ψ ( ) 60

1000

1500

500

rrx of mr

0.43

1500

2000

γ of m r

10 Experiment Simulation

0.425

20

1000 times(s)

Experiment Simulation

40

2000

time(s)

time(s)

Experiment Simulation

0 −10 −20

0

0.42

−20

−30

°

−40

−40

−60

0.415

−50

−80

−60

−100

0.41 1000

2000

3000

4000

5000

6000

−70 500

time(s)

1000

1500

2000

time(s)

500

1000

1500

2000

time(s)

Fig. 12. Comparing experimental results with simulation results. (a) Glider depth. (b) Glider pitch angle y. (c) Glider roll angle f. (d) Glider yaw angle c. (e) Position of mr along x and (f) Rotational angle g of mr .

g stayed at  601 most of the time, but were moved to 01 intermittently. The purpose for this intermittent switching was to counteract the strong ocean currents. Since at g ¼ 601, the vertical distance between the center of buoyancy and the center of mass of the glider decreased compared to g ¼ 01 that made the attitude of the glider less stable. By intermittently switching the angle g from  601 to 01, we increased the tolerance to strong current. However, when g varied, the mass distribution of the glider changed, which further led to variations of roll angle c and small variations in pitch angle y to balance the change of g. As shown in Fig. 12(a), the glider reached a depth of 800 m. The pitch angle was approximately  401 as shown in Fig. 12(b). The roll angle varied around 281 as shown in Fig. 12(c). The yaw angle changed from  501 to 501 at a relatively constant rate as shown in Fig. 12(d). These results imply that a downward spiraling motion has been successfully produced against the ocean current. From Fig. 12(b) and (c), we do see small oscillations caused by the intermittent rotation of the movable mass. We did not test uncontrolled spiraling motion because the strong ocean current often exceeded glider speed, which would have carried the glider out of the detection range of the supporting ship. There are other factors that may have contributed to the difference between experimental data and simulation results. These include water density variation, change of mass distribution, and the compression of the glider hull at greater depth. In the glider dynamic model, the static mass and movable mass are treated as rigid mass blocks. However, in the actual system, the net buoyancy is controlled by pumping oil between different chambers inside the glider. Hence the mass distribution is difficult to determine. This uncertainty affects the moment around y-axis, which may account for the disagreements between the simulation and the experiment in Fig. 12(b). The Seawing glider has two main wings aligned in the horizontal plane and two tail wings aligned in the vertical plane of the body frame. During the spiraling motion, the tail wings produce lift forces that

are orthogonal to the lift forces produced by the main wings. The combined force provides the centrifugal force for turning. The effect of the tail wing is omitted during the theoretical analysis in this paper since the size of tail wings is much smaller compared to the main wings. However, we suspect that the effect of the tail wings on turning may not be neglected in certain circumstances. The optimal ratio between the size of the tail wings and the size of the main wings is a topic for future research. But nevertheless, in these figures, we see a consistent match between the simulations and the experimental results, which provides evidence that the motion model captures the glider dynamics to a satisfactory accuracy.

6. Conclusions In summary, we have demonstrated that the dynamic model of underwater gliders can be simplified to solve for the parameters for the spiraling motion using a recursive algorithm. All hydrodynamic coefficients in the dynamic model can be computed using CFD softwares. The theoretical glider model is consistent with experimental results. We therefore conclude that the proposed recursive algorithm based on the theoretical model may be used to compute the desired control input to produce a circular helical motion in practice.

Acknowledgments This work is supported by the State Key Laboratory of Robotics, Grant no. 2009-Z05, and Knowledge Innovation Program of Chinese Academy of Sciences, Grant no. KZCX2-YW-JS205. The authors would like to thank Zhier Chen, Wenming Jin, and Yan Huang for the cooperating work in the glider experiment during July 2011. And thank Zhiqiang Hu and Haitao Gu for their help on CFD simulation.

S. Zhang et al. / Ocean Engineering 60 (2013) 1–13

References Alvarez, A., Caffaz, A., Caiti, A., Casalino, G., Gualdesi, L., Turetta, A., Viviani, R.,  2008. Folaga: a low-cost autonomous underwater vehicle combining glider and AUV capabilities. Ocean Eng. 36 (1), 24–38. Arima, M., Ichihashi, N., Ikebuchi, T., 2008. Motion characteristics of an underwater glider with independently controllable main wings. In: IEEE OCEANS 2008, Kobe, Japan, pp. 1–7. Arima, M., Ichihashi, N., Miwa, Y., 2009. Modelling and motion simulation of an underwater glider with independently controllable main wings. In: IEEE OCEANS 2009, Bremen, Germany, May, pp. 1–6. Bender, A., Steinberg, D.M., Friedman, A.L., Williams, S.B., 2008. Analysis of an autonomous underwater glider. In: Australasian Conference on Robotics and Automation, Canberra, Australia, pp. 1–10. Bhatta, P., Leonard, N.E., 2008. Nonlinear gliding stability and control for vehicles with hydrodynamic forcing. Automatica 44 (5), 1240–1250. Boger, D.A., Dreyer, J.J., 2006. Prediction of hydrodynamic forces and moments for underwater vehicles using overset grids. In: Proceedings of 44th AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, USA, pp. 1–13. Eriksen, C.C., Osse, T.J., Light, R.D., Wen, T., Lehman, T.W., Sabin, P.L., Ballard, J.W., Chiodi, A.M., 2001. Seaglider: a long range autonomous underwater vehicle for oceanographic research. IEEE J. Oceanic Eng. 26 (4), 424–436. Fossen, T.I., 2011. Handbook of Marine Craft Hydrodynamics and Motion Control. Wiley, UK. Geisbert, J.S., 2007. Hydrodynamic Modeling for Autonomous Underwater Vehicles using Computational and Semi-empirical Methods. Master’s Thesis, Virginia Polytechnic Institute and State University, Virginia. Graver, J., Leonard, N.E., 2001. Underwater glider dynamics and control. In: 12th International Symposium on Unmanned Untethered Submersible Technology, Durham, pp. 1–14. Hussain, N.A.A., Arshad, M.R., Mohd-Mokhtar, R., 2011. Underwater glider modelling and analysis for net buoyancy, depth and pitch angle control. Ocean Eng. 38 (16), 1782–1791. Isa, K., Arshad, M.R., 2011. Dynamic modeling and characteristics estimation for USM underwater glider. In: Proceedings of the 2011 IEEE Control and System Graduate Research Colloquium, Shah Alam, pp. 12–17. Isern-Gonza´lez, J., Sosa, D.H., Ferna´ndez-Perdomo, E., Cabrera-Ga´mez, J., Domı´nguez-Brito, A.C., 2011. Path planning for underwater gliders using iterative optimization. In: Proceedings of 2011 IEEE International Conference on Robotics and Automation, Shanghai, China, pp. 1538–1543. Jagadeesh, P., Murali, K., 2005. Application of low-Re turbulence models for flow simulations past underwater vehicle hull forms. J. Nav. Archit. Mar. Eng. 1 (2), 41–54. Jagadeesh, P., Murali, K., Idichandy, V.G., 2009. Experimental investigation of hydrodynamic force coefficients over AUV hull form. Ocean Eng. 36 (4), 113–118. Kan, L., Zhang, Y., Fan, H., Yang, W., Chen, Z., 2008. MATLAB-based simulation of buoyancy-driven underwater glider motion. J. Ocean Univ. China (English ed.) 7, 113–118. Leonard, N.E., Graver, J., 2001. Model-based feedback control of autonomous underwater gliders. IEEE J. Oceanic Eng. 26 (4), 633–645. Leonard, N.E., Paley, D.A., Davis, R.E., Fratantoni, D.M., Lekien, F., Zhang, F., 2010. Coordinated control of an underwater glider fleet in an adaptive ocean sampling field experiment in Monterey bay. J. Field Robot. 27 (6), 718–740.

13

Mahmoudian, N., Geisbert, J., Woolsey, C., 2010. Approximate analytical turning conditions for underwater gliders: implications for motion control and path planning. IEEE J. Oceanic Eng. 35 (1), 131–143. Mahmoudian, N., Woolsey, C., 2008. Underwater glider motion control. In: Proceedings of the 47th IEEE Conference on Decision and Control, Cancun, Mexico, December, pp. 552–557. Morgansen, K.A., Triplett, B.I., Klein, D.J., 2007. Geometric methods for modeling and control of free-swimming fin-actuated underwater vehicles. IEEE Trans. Robot. 23 (6), 1184–1199. Paley, D.A., Zhang, F., Fratantoni, D.M., Leonard, N.E., 2008. Cooperative control for ocean sampling: the glider coordinated control system. IEEE Trans. Control Syst. Technol. 16 (4), 735–744. Pereira, A.A., Binney, J., Jones, B.H., Ragan, M., Sukhatme, G.S., 2011. Toward risk aware mission planning for autonomous underwater vehicles. In: Proceedings of 2011 IEEE/RSJ International Conference on Intelligent Robots and Systems, San Francisco, CA, USA, pp. 3147–3153. Racine, B.J., Paterson, E.G., 2005. CFD-based method for simulation of marinevehicle maneuvering. In: Proceedings of 35th AIAA Fluid Dynamics Conference and Exhibit, Toronto, Ontario, Canada, June, pp. 1–22. Rodgers, J.D., Wharington, J.M., 2010. Hydrodynamic implications for submarine launched underwater gliders. In: IEEE OCEANS 2010, Sydney, May, pp. 1–8. Sherman, J., Davis, R.E., Owens, W.B., Valdes, J., 2001. The autonomous underwater glider ‘‘spray’’. IEEE J. Oceanic Eng. 26 (4), 437–446. Smith, R.N., Chao, Y., Jones, B.H., Caron, D.A., Li, P.P., Sukhatme, G.S., 2009. Trajectory design for autonomous underwater vehicles based on ocean model predictions for feature tracking. In: The 7th International Conference on Field and Service Robots (FSR 2009), MIT Campus, Cambridge, MA, pp. 263–273. Smith, R.N., Schwager, M., Smith, S.L., Rus, D., Sukhatme, G.S., 2011. Persistent ocean monitoring with underwater gliders: adapting sampling resolution. J. Field Robot. 28 (5), 714–741. Tang, S., Ura, T., Nakatani, T., Thornton, B., Jiang, T., 2009. Estimation of the hydrodynamic coefficients of the complex-shaped autonomous underwater vehicle TUNA-SAND. J. Mar. Sci. Technol. 14 (3), 373–386. Wang, S., Sun, X., Wang, Y., Wu, J., Wang, X., 2011. Dynamic modeling and motion simulation for a winged hybrid-driven underwater glider. China Ocean Eng. 25 (1), 97–112. Wang, S., Xie, C., Wang, Y., Zhang, L., Jie, W., Hu, S.J., 2007. Harvesting of PEM fuel cell heat energy for a thermal engine in an underwater glider. J. Power Sources 169 (2), 338–346. Webb, D.C., Simonetti, P.J., Jones, C.P., 2001. SLOCUM, an underwater glider propelled by environmental energy. IEEE J. Oceanic Eng. 26 (4), 447–452. Yoon, H.K., Son, N.S., Lee, G.J., 2007. Estimation of the roll hydrodynamic moment model of a ship by using the system identification method and a free running model test. Ocean Eng. 32 (4), 798–806. Yu, J., Zhang, F., Jin, W., Tian, Y., 2012. Motion parameter optimization and sensor scheduling for the sea-wing underwater glider. IEEE J. Oceanic Eng., http://dx. doi.org/10.1109/JOE.2012.2227551. Zhang, F., Fratantoni, D.M., Paley, D., Lund, J., Leonard, N.E., 2007. Control of coordinated patterns for ocean sampling. Int. J. Control 80 (7), 1186–1199. Zhang, S., Yu, J., Zhang, A., Zhang, F., 2011. Steady three dimensional gliding motion of an underwater glider. In: Proceedings of the 2011 IEEE Conference on Robotics and Automation, Shanghai, China, pp. 2392–2397.