Robust Adaptive Control of Underwater Vehicles with Biologically Inspired Vortex Ring Thrusters

Robust Adaptive Control of Underwater Vehicles with Biologically Inspired Vortex Ring Thrusters Y. Xu*. K. Mohseni*, ⋆, †. *Department of Mechanical and Aerospace Engineering, University of Florida, Gainesville, FL 32611-6250 USA (e-mail: {yxu, mohseni}@ufl.edu). ⋆Department of Electrical and Computer Engineering, University of Florida, Gainesville, FL 32611-6250 USA. †Institute for Networked Autonomous Systems, University of Florida, Gainesville, FL 32611-6250 USA. Abstract: A novel underwater thruster was designed to provide quantized level of thrusts for underwater vehicle maneuvering with a minimal effect on the drag profile of the vehicle. The design was inspired by the efficient and agile locomotion mechanism in squid and jellyfish, which is characterized by formation of vortex rings from pulsatile jets. There have been several successful implementations of the device on a series of autonomous underwater vehicles in our group. In this paper, a robust adaptive control algorithm is formulated for future vehicles with new, more compact thrusters of this type. Control strategy for motion in 6 degrees of freedom, together with corresponding force allocation method is developed and is tested in simulation, where an asymptotic position and velocity tracking result is obtained. Keywords: Autonomous underwater vehicle, vortex ring, propulsion, robust control, adaptive control. 1. INTRODUCTION Due to increased applications of underwater vehicles, performance and safety in position and velocity control systems is receiving more focus (Fossen and Blanke, 2000). Accurate underwater vehicle control needs to overcome poorly modelled hydrodynamic forces that are dependent on relative vehicle velocities, and to stabilize against all kinds of environmental disturbances. Current underwater vehicles may fall into two categories. Torpedo shaped underwater vehicles are characterised by long streamlined bodies for optimized forward motion, generating steering forces with the help of control surfaces. This type of vehicle is efficient in long-distance high-speed travels, but lacks agility and positioning accuracy at lower speeds since the control surfaces cannot provide deflective forces without a steady relative velocity of flow. The other type of underwater vehicles has a box shaped body that utilizes multiple thrusters at different locations to reinforce low-speed maneuverability. However, this type of design impedes the vehicle’s ability to travel at high speeds, because thrusters at multiple locations may produce drag and other undesirable effects as the vehicle traveling forward. Even for low-speed manoeuvres, traditional propellers have poor performance for small corrections that are on the level of a single shaft rotation. (Yoerger et al., 1990; Mclean, 1991.) Alvarez et al. (2009) consider an energy-efficient design that combines oceanographic gliders with steering jet-pumps, which allows for maneuvering even at low speed, but the lack of control surfaces also prevents navigation without propulsion. Research studies have shown that locomotion mechanism in cephalopods and jellyfishes (as in Fig. 1 from Krieg and Mohseni, 2008) allows them to perform efficient and accurate manoeuvres in cluttered marine environment (Sahin and Mohseni, 2009; Sahin, Mohseni, and Colins, 2009; Lipinski

and Mohseni, 2009). A novel scheme for compact pulsatile jet actuators was proposed by Mohseni (2004). The device was later referred as a vortex ring thruster (VRT), which is capable of producing quantized propulsive force without significantly affecting the drag profile of the vehicle. The VRT (as illustrated in Fig. 2) consists of a large cavity that is similar to the squid mantle. Inside, there is a fluid manipulator which can change the volume of cavity and force fluid in and out through an opening. The thruster can create an array of high-momentum vortex rings by successive ingestion and expulsion of fluid. Both theoretical and experimental studies were carried out in understanding the mechanisms of vortex formation for development of VRTs, which suggests that VRT can provide adjustable thrusts for accurate slow-speed maneuvering and positioning, and for flow control and drag reduction at cruising speeds (Mohseni et al., 1998; Thomas et al., 2003; and Mohseni, 2006). Since the VRT has no protruding components that increase drag, it does not affect forward travel ability of the vehicle in comparison to other current propellers.

Fig. 1. Conceptual diagram of locomotion for squid A number of VRTs have been successfully implemented on a series of automatic underwater vehicles (AUV), namely, Calmar-E, Kraken, and Hydra. As new VRT designs are more efficient and compact, they are to be implemented and tested on a new generation of AUVs (as shown in Fig. 3). In this paper, a robust adaptive control scheme is developed aiming at accurate control of AUVs with VRTs.

Transformation between the coordinate systems is given by

Flexible cavity

(2)

� � J (� )� , � � J �1 (� )� .

Ejected fluid

(a)

It is assumed that the dynamic equation of motion can be written in the following form. (Definitions and properties are found in Fossen, 1994) (b)

M� � C (� )� � D(� )� � g (� ) � � ;

Fig. 2. (a) Schematic diagram and (b) actual fabrication of the VRT unit ③ ④

⑤ ①

② 0.3 m (1 ft.)

⑥

⑦

Fig. 3. Current design of underwater vehicle in our group. ①Buoyancy device ②Propeller ③VRT ④RF antenna plug ⑤Wetmate plug ⑥Camera ⑦Acoustic transmitter/receiver. The paper is arranged as follows. In Section 2, a general model for underwater vehicles is presented, together with models for the VRTs and other actuating units. In Section 3, control design and force allocation algorithm is formulated. Simulation results of vehicle control are presented in Section 4. Conclusions and discussions are in Section 5. 2. VEHICLE MODEL AND PROPERTIES

To describe motion of the underwater vehicle model, two reference frames are introduced. The coordinates frame fixed on the Earth is considered to be inertial. A body-fixed frame is conveniently defined on the underwater vehicle. The 6 degrees of freedom motion of the vehicle can be described by position and orientation vector � � 6 in the earth-fixed frame, together with linear and angular velocity vector � � 6 in the body-fixed frame. (1)

� � [ x, y, z, � ,� ,� ] , � � [u, v, w, p, q, r ] . T

T

Components in the vectors are defined according to Table 1.

Table 1. Notation used for modelling Mode

Descriptions

1 2 3 4 5 6

Motion in x direction Motion in y direction Motion in z direction Rotation about x axis Rotation about y axis Rotation about z axis

Linear and Position angular and Euler velocities angles u v w p q

r

where M � 6�6 denotes the inertia matrix that is symmetric and positive-definite. It includes inertia of rigid body and added mass. One should note that changes in buoyancy may also affect the rigid body inertia. However, since the changes are part of the control input, which is known, the inertia matrix is considered as a known function with uncertain constant parameters. Rewrite the inertia matrix as

M � [M1T

x y

M 2T ]T , M1 , M 2 �

3�6

(4)

.

Then the Coriolis and centripetal matrix C (� ) � formulated as � 0 C (� ) � � 3�3 � � �[ M1� ]

6�6

can be

�[ M1� ]� � �; �[ M 2� ]� �

(5)

where [�]� is defined such that a � b � [a]� b for all a, b � that is � 0 [a] � �� a3 �� �a2 �

2.1 Body-fixed Vehicle Model

(3)

�a3 0 a1

a2 � �a1 �� , �a � [a1 , a2 , a3 ]T � 0 ��

3

3

,

(6)

.

D(� ) � 6�6 denotes hydrodynamic damping terms. For an underwater vehicle moving at low speed, it can be assumed that the damping effects are non-coupled, and that offdiagonal terms and terms higher than second order are negligible. Therefore the damping matrix is considered to be diagonal with only linear and quadratic damping terms. D(� ) � D1 D� � D2 , D� � diag � u , v , w , p , q , r

�;

where diag{} � denotes a diagonal matrix, and D1 , D2 � are constant diagonal matrices.

(7) 6�6

Vector g (� ) � 6 represents restoring forces and moments that results from gravity and displacement of fluid, and � � 6 represents propulsion forces and moments. Note that effect from the buoyancy units on the vehicle is considered as propulsion forces, therefore it is excluded in g (� ) . 2.2 Vehicle Model Properties

z �

The dynamic model described by (2) and (3) has several properties that are useful in subsequent development and analysis (Fossen, 1994).

� �

a) There exist positive constants cm1 , cm 2 �

cm1 �

2

2

� � T M � � cm2 � , �� �

�

such that

6

;

(8)

where � stands for standard Euclidean norm. b) Time derivative M� (� ) of the inertia matrix and the centripetal-Coriolis matrix C (� ) satisfies

� T [ 12 M � C(� )]� � 0, �� ,� �

6

(9)

.

c) It is assumed that there are positive constants such that the following relationships are true:

C (� a ) � cc1 � a , C (� a ) � C (� b ) � cc 2 � a �� b , D(� a ) � cd 1 � a � cd 2 ,

�� ,� a ,� b �

6

;

(10)

D(� a ) � D(� b ) � cd 3 � a �� b , �

� f � (13) TBi � � Bi � , i � 1, 2. � r f Bi � � Bi b) Propeller: It is assumed that the propeller can provide desired force through an independent control unit, thus its dynamic model is not considered here. c) VRT: Dynamics for the VRT, as shown in Fig. 5 and Fig. 6, can be characterized by combination of a first-order damper with settling time t s and a sinusoidal mode in pulsation frequency f (Krieg and Mohseni, 2010).

TVRT (t ) � cth f 2 {[1 � exp(� ft ts )] � cw sin(2� ft )}.

Q(� ,� a ) � Q(� ,� b ) � cq � a �� b , where cc1 , cc 2 , cd1 , cd 2 , cd 3 , cq �

For i th buoyancy unit, the force can be expressed in vector of buoyancy center rBi � 3 (with respect to the origin) as

. Q(� ,� ) �

6�6

is defined to

This model can be simplified as

TVRT (t ) � cth f 2 .

be time derivative of the inverse transformation matrix: (11)

Q(� ,� ) � J �1 (� ).

In Appendix A, a possible way to obtain the constants in (10) will be illustrated, especially for the afore-described vehicle model. Hence, these inequalities may be considered as properties more than assumptions.

(14)

(15)

Forces and moments generated by the actuators depend on the geometry of installation. The allocation for control forces and moments will be further discussed in 3.3.

Models of the actuators are necessary for the control design so that the control forces can be properly allocated to each of the actuating units. Current design of the underwater vehicle has 7 actuators: a rear-facing propeller for forward movement, four VRTs on both sides for accurate maneuvers, and two buoyancy units on both ends of the vehicle (see Fig. 4). The vehicle is so designed that its x-z plane always stays upright, that is � � 0 for X-Z-Y rotations. Origin of the body-fixed frame is selected as midpoint of the buoyancy units.

y

Propeller x

3

Buoyancy units

1

2

1

4

2

Vortex ring thrusters

Time (s) Fig. 5. Average thrust response of VRT from ignition though several actuation cycles (from Krieg and Mohseni, 2008)

Thrust (N)

x

Thrust (N)

2.3 Actuator Models

z

Fig. 4. Schematic diagram of actuating units on the vehicle a) Buoyancy: Forces from buoyancy units can be transformed into the body-fixed reference frame, and are defined as

f Bi � J Bi (� ) Bi ,

f Bi , J Bi (� ) �

3

, i � 1, 2;

(12)

where Bi is the magnitude of force generated by each of the buoyancy units in the earth-fixed frame.

Frequency (Hz) Fig. 6. Average actuator thrust from VRT versus frequency (from Krieg and Mohseni, 2008)

3. CONTROL DEVELOPMENT

regression matrix of known functions in the position vector � , the virtual control vector � b , and vectors of the desired trajectory � d , � d , and � d .

3.1 Control Design An adaptive control scheme is integrated to robust control strategy, such that uncertainties that are linear in the errors can be eliminated, and that an asymptotic tracking, rather than a uniformly ultimately bounded result, can be obtained (Khalil, 2001). The control objective is to track a continuous bounded desired trajectory �d ,�d ,�d � 6 with measurement of the position and orientation vector � and the velocity vector � (assumed as measurable). Vectors of the desired trajectory are bounded by

�d � c�1 ,

�d � c� 2 , c�1 , c� 2 �

�

.

(16)

Performance of the controller is represented by position and orientation tracking error e � 6 . e � � � �d .

(17)

� � �Y�ˆ � ke e � � k� � � � k� � � � k f � � J T (� )e;

are control gains, and �ˆ � 6 represents the adaptive vector for the system parameters, whose adaptation law is designed to be where ke , k� , k� , k f �

�

(26)

�ˆ � Y T ��. ��

6�6

is a constant diagonal matrix.

Consequently, closed-loop dynamics for the back stepping error can be written as

M � � �C (� ) � � [C (� ) � C (� b )]� b � [ D(� ) � D(� b )]� � D(� b ) � � kb M [Q(� ,� ) � Q(� ,� b )](� � � d )

(27)

� kb M � � M [Q(� ,� ) � Q(� ,� b )]� d � Y�

(18)

where � �

6

represents the mismatch of the adaptive vector: (28)

ˆ � � � � ��

Define virtual control signal as

� b � �kb J �1 (� )e � J �1 (� )�d , kb �

�

.

(19)

3.2 Stability Analysis

(20)

Performance of the control design can be evaluated using Lyapunov analysis (see Khalil, 2001).

Define back stepping error as

� � � �� b .

Closed-loop dynamics for tracking error can be obtained by combining (19) and (20) into (18). e � J (� )� b � J (� )� ��d � �kb e � J (� )�.

(25)

� ke e � � k� � � � k� � � � J T (� )e;

Time derivative of the tracking error can be written as e � J (� )� � �d .

The control force is designed to be

Define Lyapunov candidate function to be (29)

V � 12 eT e � 12 � T M � � 12 � T � �1� ;

(21)

which is positive definite, radially unbounded, and decrescent.

Substituting (17) and (18) into time derivative of (19) yields dynamics for the virtual control signal.

Substitute (21), (26), (27), and (28) into the time derivative of (29). The following inequality can be obtained using (10).

� b � �kb Q(� ,� )e � kb J �1 (� )e � Q(� ,� )�d � J �1 (� )�d

2

(22)

� �kb Q(� ,� )(� � �d ) � kb� � kb J �1 (� )�d � Q(� ,� )�d � J �1 (� )�d .

Open-loop dynamics for the back stepping error can be derived from (20) using (3) and (22).

M � � M� � M� b � �C (� )� � D(� )� � g (� ) � � � kb MQ(� ,� )(� � �d ) � kb M� � kb MJ (� )� d �1

(23)

� MQ(� ,� )�d � MJ (� )� d . �1

2

3

(30)

2

� (k� � �� ) � � � ( k f � � f ) � ;

where �e , � � , �� �

�

are known constants.

�e � kbcm 2cq , � � � cc 2 � cd 1, �� � cc 2 � cd1 � cd 3 , � f � cd 2 � kbcm 2 � cm 2cqc�1.

(31)

If the control gains are chosen such that

ke � �e , k� � �� , k� � �� , k f � �� � ka �

�

;

(32)

then (30) can be written as

The dynamics in (23) can be linearly parameterized into

2

V � �kb e � ka �

Y� � �C (� b )� b � D(� b )� b � g (� ) � kb MQ(� ,� b )(� � �d ) � kb M� b � kb MJ �1 (� )�d

2

V � �kb e � (ke � �e ) e � � (k� � � � ) �

(24)

� MQ(� ,� b )�d � MJ �1 (� )�d ; where � � p denotes a vector containing constant system parameters of the dynamic system, and Y � 6� p denotes the

2

� � min{kb , ka } eT

2

�T .

(33)

The function V (e(t ), � (t ),� (t )) is negative semi-definite, and therefore it is bounded, that is, V (e(t ), � (t ),� (t )) � L� ; hence,

e(t ), � (t ),� (t ) � L� . By using (17), (19), (20), and (28) it can be proved that � (t ),� (t ),� (t ),�ˆ (t ) � L . From (18) and (25), b

�

one may find that � (t ), e(t ) � L� . As a result, the control force � is defined and bounded, and the position tracking error e is uniformly continuous. It can be shown from (33) that e(t ) � L2 . By Barbalat’s lemma, it can be concluded that

lim e(t ) � 0. t ��

(34)

Consequently, a globally asymptotically tracking result can be obtained with this control design. 3.3 Control Allocation The desired control forces and moments from the aforedescribed control design are allocated to the actuators.

ke � 1530, k� � k� � 105, k f � 1670.

(41)

The desired trajectory is defined as movements in x, y and z directions, and rotation about z axis, as described with hyperbolic tangent functions. xd � yd � zd � � d � 0.1tanh(t � 5) � 0.1tanh(5), �d � 0

(42)

Simulation results are shown in Fig. 7, Fig. 8, and Fig. 9, which only have errors whose magnitudes are small in comparison with that of the desire trajectory (0.2). Position error (m)

Position error (rad)

Define the desired forces and moments as

� � [�1 ,� 2 ,� 3 ,� 4 ,� 5 ,� 6 ]T .

(35)

As mentioned in 2.3, the vehicle is passively stable at � � 0 (X-Z-Y). Note that Mode 6 should not be activated without � � 0 otherwise this stability might be affected by VRTs. a) Buoyancy: The buoyancy units can only affect Mode 1, 3 and 5 (see Table 1). Mode 1 is also affected by the propeller. The controlled buoyancy can be found by solving (36).

J B 2 (� ){3} � � B1 � �� 3 � � J B1 (� ){3} �[r � J (� )]{2} [r � J (� )]{2}� � B � � �� � ; � B1 B1 �� 2� � 5� B2 B2

Time (s)

Time (s)

Fig. 7. Position error from the desired trajectory versus time Velocity error (m/s)

Velocity error (rad/s)

(36)

where {} � denotes the index of a component in a vector. b) Propeller: Force from propeller can be determined from

TP � �1 � J B1 (� ){1}B1 � J B 2 (� ){1}B2 .

c) VRT: Due to the location and orientation of VRTs, it is convenient to divide them into two groups, one at either end, so that the force allocation is similar to that of buoyancy. TVRTa � TVRT 1 � TVRT 2 , TVRTb � TVRT 3 � TVRT 4 .

Time (s)

Time (s)

(37) Fig. 8. Velocity error for desired trajectory versus time Control force (N)

Control force (N)

(38)

Thus, forces of the VRTs can be determined by

� 1 1 � �TVRTa � �� 2 � � (39) �l l � �T � � �� � , la � �lb � l � ; � a b � � VRTb � � 6 � where l is longitudinal distance of a VRT to center of mass. 4. SIMULATION TESTS Simulation tests are carried out to verify the convergence of the control design. Dynamic parameters for the preceding vehicles are used, and are defined in (40). 0 0 0 1.81 0 � �37.02 � 0 � 74.96 0 1.81 0 0.49 � � � � � 0 0 74.96 0 0.49 0 � M �� �, 0 0.41 0 0 � �1.81 � 0 � 1.81 0 0.49 0 9.01 0 � � � � 0 0.49 0 0 0 9.01 � � D � diag{6.79, 133.42, 133.42, 0.001, 7.11, 7.11}, rG � [0 0 0.5]T , rB1 � � rB 2 � [0.34 0 0]T , l � 0.324.

The control gains are chosen accordingly as

(40)

Time (s)

Time (s)

Fig. 9. Control force for desired trajectory versus time 5. CONCLUSIONS In summary, a robust adaptive control scheme is designed and tested in simulation to operate an underwater vehicle with VRTs. Implementation of VRTs is expected to reduce the need for solving coupled relationship between actuators as in vehicles with propellers and fins. Its compact structure and ability to generate quantized level of thrusts enables the vehicle to perform accurate maneuvers at low speed, without sacrificing the streamlined profile for efficient high speed traveling. Further development for the control algorithm is needed before implementation since sensor noises, actuator capacities and other practical aspects should be considered.

REFERENCES

Appendix A. DERIVATION OF UPPER BOUNDS

Alvarez, A., Caffaz, A., Caiti, A., Casalino, G., Gualdesi L., Turetta, A., and Viviani, R. (2009). Fòlaga: A low-cost autonomous underwater vehicle combining glider and AUV capabilities, Ocean Engineering, vol. 36, pp. 24-38. Fossen, T. I. (1994). Guidance and Control of Ocean Vehicles. pp. 4-52, John Wiley & Sons Ltd. Fossen, T. I., and Blanke, M. (2000). Nonlinear Output Feedback Control of Underwater Vehicle Propellers using Feedback from Estimated Axial Flow Velocity. IEEE J. Ocean. Eng., vol. 25, no. 2, pp. 241-255. Horn, R. (1990). Matrix Analysis, 3rd edn, Chap. 1, 5, and 7, Cambridge University Press, NY. Khalil, H. (2001). Nonlinear Systems, 3rd edn, Chap. 3-5, and 14, Prentice Hall, NJ. Krieg, M., and Mohseni, K. (2008). Thrust characterization of pulsatile vortex ring generators for locomotion of underwater robots. IEEE J. Ocean. Eng., vol. 33, no. 2, pp. 123-132. Krieg, M., and Mohseni, K. (2010). Dynamic modeling and control of biologically inspired vortex ring thrusters for underwater robot locomotion. IEEE Trans. Robotics, vol. 26, no. 3, pp. 1212-1242. Lipinski, D., and Mohseni, K. (2009). A numerical investigation of flow structures and fluid transport with applications to feeding for the hydromedusae Aequorea victoria and Sarsia tubulosa. J. Exp. Biol., vol. 212, pp. 2436-2447. Mclean, M. (1991). Dynamic performance of small diameter tunnel thrusters. Ph.D. dissertation, Naval Postgrad. School, Monterey, CA. Mohseni, K., and Gharib, M. (1998). A model for universal time scale of vortex ring formation. Phys. Fluids., vol. 10, no. 10, pp. 2436-2438. Mohseni, K. (2004). Zero-mass pulsatile jets for unmanned underwater maneuvering. Paper presented at 3rd AIAA Unmanned Unlimited Tech. Conf., Workshop Exhib., Chicago, IL, Sep. 20-23, AIAA Paper 2004-6386. Mohseni, K. (2006). Pulsatile vortex generators for lowspeed maneuvering of small underwater vehicles. Ocean. Eng., vol. 33, no. 16, pp. 2209-2223. Sahin, M., Mohseni, K., and Colins, S. (2009). The numerical comparison of flow patterns and propulsive performances for the hydromedusae Sarsia tubulosa and Aequorea victoria. J. Exp. Biol., vol. 212, pp. 2656-2667. Sahin, M., and Mohseni, K. (2009). An arbitrary LagrangianEulerian formulation for the numerical simulation of flow patterns generated by the hydromedusa Aequorea victoria. J. Comput. Phys., vol. 228, pp. 4588-4605. Thomas, T., Krieg, M., and Mohseni, K (2003). A Thrust characterization for bioinspired pulsatile vortex ring thrusters with variable exit nozzle diameter. Paper presented at the ASME Int. Mech. Eng. Congr. Expo., Buena Vista, FL, Nov. 13-19. Yoerger, D., Cooke, J., and Slotine, J.-J. (1990). The influence of thruster dynamics on underwater vehicle behavior and their incorporation into control system design. IEEE J. Ocean. Eng., vol. 15, no. 3, pp. 167-178.

In this appendix, a possible method to derive the constants in the inequalities in (10) is presented. It should be noted that this is not intended to provide a rigorous proof for the inequalities. Useful properties can be found in Horn (1990). From (9), norm of the inertia matrix can be upper bounded by M � cm 2 .

(41)

To derive the upper bounds for C (� ) , it is useful to examine the equation for the Euclidean norm of a matrix A � 6�6 , which is its largest singular value. A � max � A� : � �

Define unit vector x �

For any vector � �

6

x2T ]T , x1, x2 �

(42)

3

(43)

.

, it is true from (5) that 2

�( M 1� ) � x2

2

C (� ) x �

with � � 1 �.

as

6

x � [ x1T

6

�( M 1� ) � x1 � ( M 2� ) � x2

� M 1� � M�

2 2

2

x2 � M 1�

�x

1

2

� x2

2

2

2

x1 � M 2�

�� M

which yields

C (� ) � max � C(� ) x

For all vectors � a ,� b �

2

�

��c

m2

2

2

x2

(44)

2 2

� cm2 2 � ;

� .

C(� a ) � C(� b ) � C(� a �� b ).

(46)

Substitute into (10) gives C(� a ) � C(� b ) � C(� a �� b ) � cm 2 � a �� b . For A �

m�n

A�

(45)

6

(47)

, the following inequalities hold true. n � A � m A �,

A � n A 1.

(48)

Of course (48) also applies for vectors since they may be considered as matrices as well. Therefore D� � 6 diag�u, v, w, p, q, r� � � 6 �

�

� 6� .

(49)

According to (7), D(� ) � 6 D1 � � D2 .

For all vectors � a ,� b �

(50)

6

D� a � D� b � 6 diag � ua � ub , va � vb , wa � wb , pa � pb , qa � qb , ra � rb

�

�

� 6 � a �� b .

Similarly, it can be deducted that D(� a ) � D(� b ) � 6 D1 � a �� b . Q(� ,� ) can be arranged as a multiplication of Q� (� ) �

(51)

(52) 6�6

that consists of trigonometric functions in the Euler angles, and matrix Q� (� ) � 6�6 in the velocity vector � . Q(� ,� ) � Q� (� )Q� (� ).

(53)

The constants for Q(� ,� ) in (11) can be obtained accordingly.