Computation of the Added Masses of an Unconventional Airship

This paper presents a modelling of an unmanned airship.We are studying a quadrotor flyingwing. The modelling of this airship includes an aerodynamic study. A special focus is done on the computation of the added masses. Considering that the velocity potential of the air surrounding the airship obeys the Laplace’s equation, the added masses matrix will be determined by means of the velocity potential flow theory. Typically, when the shape of the careen is quite different from that of an ellipsoid, designers in preprocessing prefer to avoid complications arising frommathematical analysis of the velocity potential. They use either complete numerical studies, or geometric approximation methods, although these methods can give relatively large differences compared to experimental measurements performed on the airship at the time of its completion. We tried to develop here as far as possible the mathematical analysis of the velocity potential flow of this unconventional shape using certain assumptions. The shape of the careen is assumed to be an elliptic cone. To retrieve the velocity potential shapes, we use the spheroconal coordinates. This leads to the Lamé’s equations. The whole system of equations governing the interaction air-structure, including the boundary conditions, is solved in an analytical setting.


Introduction
The rapid expansion of the capabilities of airships in the last decade and the growing of the range of missions they designed to support lead to the design of complex shapes of the careens.
Exhaustive studies in that topic were presented by 1, 2 .Traditionally, ellipsoidal shapes are used for airships 3-5 .However, and in order to optimize their performances, different original shapes have been tested in the last years.This is due to the advances in aerodynamics, control theory, and embedded electronics.The airship studied here departs with the traditional shapes.The MC500 is a flying wing Figure 1 , developed by the French network DIRISOFT.
The MC500 is an experimental prototype for a set of great innovating airships.A precise dynamics model is needed for this kind of unmanned airships including the airstructure interaction.This will enable the elaboration of convenient algorithms of control, stabilization, or navigation of these flying objects.
The aerodynamic forces have a large influence on the dynamic behaviour of the airships.Among the aerodynamic solicitations, the added masses phenomenon is one of the most important.In fact, when hovering or in case of low speed displacement, the lift and drag forces and torques could be neglected.
The added masses phenomenon is well known for airships and similarly for submarines.When an airship moves in an incompressible and infinite inviscid fluid, the kinetic energy of the fluid produces an effect equivalent to an important increase of the mass and the inertia moments of the body.This contribution may be of the same magnitude as the terms of mass or inertia of the airship.Apart from ellipsoidal shapes 6 or thin rectangular plates 7 where the theory is well established for many decades, the analysis of this problem for unconventional shapes in preprocessing is usually done by approximate methods.We can quote the geometric method, consisting of an extension of a 2D analysis and requiring the introduction of empirical parameters 8, 9 , or the fully numerical methods consisting of modeling both the airship and the surrounding air 10, 11 .As opposed to other treatments of this problem, the derivation proposed here is based on the velocity potential flow theory 6, 7, 12 .We tried to pursue the analytical study to an advanced stage with some assumptions.The shape of the careen of the MC500 is considered as a cone with elliptic section.When we consider that the velocity potential of the air obeys the Laplace's equation, the added masses matrix could be determined by solving this equation, using the spheroconal coordinates.Solving the Laplace's equation for such configuration leads to the Lamé's equations.We have dedicated a section of this paper for the determination of the ellipsoidal harmonic series issued from these equations with the specific boundary conditions.Such series developed for the first time by Lamé 13 were improved particularly by Liouville and Sturm in their famous theory and by Hermite 14 and were applied in different fields.Byerly 15 and later Hobson 16 gave an important comprehensive literature about this topic.
Significant recent works, such as the works of Boersma and Jansen 17 in electromagnetic field, and those of Garmier and Barriot 18 in astrophysics could be quoted.In this work, we are trying to apply this theory for the air-structure interaction.

Kinematics
The airship MC500 is assumed to be a rigid flying object.
We use two reference frames.First an earth-fixed frame R 0 O, X 0, Y 0 , Z 0 assumed to be Galilean.Then a local reference frame R m G, X m, Y m , Z m fixed to the airship.Its axes are selected as follows: X m and Y m are the transverse axis of the airship, Z m : the normal axis directed downwards.
To describe the orientation of the airship, we use a parameterization in yaw, pitch, and roll.The configuration of the object is described by means of three rotations defined by three angles of orientation, that is, the yaw ψ, pitch θ, and roll φ.
The whole transformation between the local reference frame R m and the global frame is given by Goldstein 19 : We denote by cθ cos θ; sφ sin φ.Using the rotation matrix J 1 η 2 , the expression of the linear speed in the reference frame R 0 is given by The angular speed is expressed as follows: It is noticed that the parameterization by the Euler angles has a singularity in Θ π/2 kπ.This parameterization is acceptable because it is not possible for an airship to reach this singular orientation of 90 degrees pitching angle.

Dynamics
As mentioned previously, a specificity of the lighter than air vehicles is illustrated by the added masses phenomenon.When a big and light object moves in the air, the kinetic energy of the particles of air produces an effect equivalent to an important growing of the mass and inertia of the body.As the airship displays a very large volume, its added masses and inertias become significant.
The basis of the analysis of the motion of a rigid body immersed in a perfect fluid is established by Lamb 6 .He proves that the kinetic energy T of the body with the surrounding fluid can be expressed as a quadratic function of the six components of the translation and rotation velocity as follows: where M a is the added mass matrix due to the motion of the surrounding air and M b is the mass matrix of the body.
For an airship fully immersed in the air, the 6 × 6 added mass matrix M a is assumed to be positive and definite.As the added kinetic energy T a , it will be discussed in the next section.The whole mass matrix M is assumed symmetric block-diagonal matrix: For the computation of the whole dynamics model, we choose to use the Kirchoff's equation 20 : where ∧ is the vectorial product, τ 1 and τ 2 are, respectively, the external forces and torques, including the rotors effects, the weight mg , the buoyancy B, and the aerodynamic lift F L and drag F D .
The dynamical system of the airship becomes 21-23

Description of the Rotors
The MC500 has four electric engines driving rotors.Each rotor has two parallel contrarotating propellers to avoid any aerodynamic torque Figure 2 .The rotor can swivel in two directions.A rotation of angle β i around the Y axis −180 • ≤ β i ≤ 180 • and a rotation of angle γ i around an axis Z iR normal to Y and initially coinciding with the Z axis −30 • ≤ γ i ≤ 30 • .A fictive axis X iR completes the rotor frame.
Let us denote P i by the position of the rotor i.We can then define a rotation matrix J i 3 between the frame P i , X iR , Y, Z iR and the local frame R m such as 2.9 If we suppose that the intensity of the thrust force of the rotor i is F i , this force could be introduced in the second member of the dynamic equation as where e X m is a unitary vector along the X m axis.
The torque produced by this rotor in the centre of inertia G is F i ∧ P i G.

Weight and Buoyancy
An important characteristic of the airships is the buoyancy B u .This force represents a natural static lift, corresponding roughly to 1 Kg for each m 3 of helium involved in the careen.We suppose here that this force is applied in the centre of buoyancy B different from the centre of inertia G: where V is the volume of the careen, ρ air is the density of the air, and g the gravity.
Let us note F WB and M WB the force and the moment due to the weight and buoyancy.We have

2.12
Aerodynamic Forces F A Such as other flying objects, the airships are subjected to aerodynamic forces.The resultant of these forces could be divided into two component forces, one parallel to the direction of the relative wind and opposite to the motion, called Drag, and the other perpendicular to the relative wind, called Lift.The MC500 is designed with an original shape oriented to a best optimization of the ratio lift upon drag forces.However, and as first study, we try to evaluate the behaviour of the airship in the case of low velocity or while hovering.In these cases, the effect of these forces could be neglected.
The global dynamic system could be expressed in a compact form as follows 21, 24 : 2.15

Spheroconal Coordinates
Much of the current airships developed and studied in the literature are considered as ellipsoidal.This particular shape has a wide popularity in this field due to the existence of a large and reliable knowledge and experimentations for both ellipsoidal airships and their alter ego submarines.However for our flying wing airship, we should use a more convenient approximation for the aerodynamic study.The MC500 airship is a collection of aerodynamic profiles, with symmetry around the x-z axis.Each transverse section of the careen parallel to the plan y-z gives roughly an ellipse.This motivated us to model as first assumption the airship as an elliptic cone Figure 3 .
To take into account the interaction of the airship with the surrounding fluid, a model of the flow is needed.
A variation of this study may be performed by calculating the fluid pressure around the airship through the Bernoulli equation 25 .
Here, we use the potential flow theory with the following assumptions.
a The air can be considered as an ideal fluid with irrotational flow and uniform density ρ f , that is, an incompressible fluid.
b A velocity potential Φ f exists and satisfies the Laplace's equation throughout the fluid domain: and satisfies the nonlinear free surface condition, body boundary condition, and initial conditions.
Finally we suppose that the velocity of the air is null far from the airship: The velocity of the fluid obeys the equation: In his study 6 , Lamb proves that the kinetic energy T a of the fluid surrounding a moving rigid body can be expressed as a quadratic function of the six components of the translation and rotation velocity ν u, v, w, p, q, r T and can also be expressed as function of the velocity potential of this fluid by the following relation: n O is an outward normal vector.Φ f can be splitted on where Φ appears as a collection of six spatial potential Φ i functions of x, y, z that verify the Laplace's equation.To extract the terms of the added mass matrix, we should define the spatial potential Φ of the moving fluid.In some cases it could be determined using geometric means 8, 26 .Here we choose to use the potential flow theory.We can then express the components of the added mass matrix M a in function of the velocity potential flow of the air surrounding the airship Φ such as For solving the Laplace's equation 3.1 and according to the configuration of the careen, we use the spheroconal coordinates 16 .This assumption seems acceptable in that stage.
A first parameterisation by the coordinates ρ, Θ, ϕ is used, for example, by Boersma and Jansen 17 , Kraus and Levine 27 Figure 4 , where ρ is the distance of a point to the origin, ϕ is the azimuthal angle 0 ≤ ϕ ≤ 2π , and Θ is the longitudinal angle.Θ Θ 0 represents the external surface of the elliptic cone.
This parameterisation has the advantage to be physically significant.However, in our case, it leads to intractable calculations.We prefer to use the parameterisation ρ, μ, ζ given by Hobson 16 .The equation of the surface of an elliptic cone is given by G. A. Korn and T. M. Korn 28 such as a, b, and μ 0 are parameters that define the elliptic cone.With ρ is always the distance of a point to the origin.We define then the Cartesian coordinates x, y, z as 3.9 By analogy with the azimuthal angle in the previous parameterisation, we can introduce an angle ϕ defined by The surface of the elliptic cone is now defined by μ μ 0 .By varying this parameter, one can see its influence on the shape of the cone in Figure 5.
To express the components of the fluid velocity V f in the conical elliptic reference frame, we define three unitary vectors e ρ , e ϕ , e μ .With h i are the scale factors.One can easily verify that these three angles form an orthonormal basis, and the velocity of the fluid could be expressed as

3.12
The surface of the elliptic cone is defined by μ μ 0 ; e μ appears as the normal vector to the surface in a given point.

Laplace's Equation
With the spheroconal coordinates, the Laplace's equation 3.1 becomes 16

3.14
To solve the Laplace's equation 3.13 , we use the traditional separation of variables: Replacing 3.16 in 3.12 gives This equation can be splitted on three separated differential equations:

3.17
And by replacing α, β, and γ by their expressions in 3.14 , one can obtain

3.18
The solutions of these equations are ellipsoidal harmonics or Lamé's functions.The Lamé's equations are usually written in a general form as For a given n integer, we can find z a 2 b 2 • s such that the particular solution of 3.19 can be written as 18

3.20
Those solutions are the Lamé's functions of first kind with a j x 2j is Lamés s polynomial, 3.21 where m depends on the integer k, as shown in Table 1: , if n is odd.

3.22
ψ n x is called the principal product.There are four different Lamé's functions that differ from their characteristics.We present them in Table 1.
z i is an eigenvalue.We remark that z i corresponds to a value of the parameter z in 3.19 , which determines as many equations as there are values of z i .
For a given value of n, the unknowns are the coefficients a j 3.23 and the parameters z i .We consider as a first step the functions K z n x ; the other functions could be defined similarly.
If we introduce the expressions 3.20 and 3.21 in 3.19 , we obtain the following recurrence's relation: or λ j a j−1 j − z a j σ j a j 1 0 for j 0, . . ., k.

3.24
Knowing that a k 1 0, the iteration will stop at rank k.And if we introduce the conditions: the whole recurrence's relations 3.28 can be written in a matrix form such as or Ω K is a square matrix of dimension k 1 .The vector Λ is an eigenvector associated to the eigenvalue z.
There are basically 2k 1 eigenvalues and eigenvectors.In fact we show that it is possible to find a diagonal matrix D and a symmetric matrix S K such as , where

3.29
The matrix S k is diagonalizable such as S k R T • Δ • R, where Δ is diagonal and R is an orthogonal matrix.Accordingly, Ω K is diagonalizable and admits 2k 1 separate eigenvalues z i associated to k 1 eigenvectors Λ i , so k 1 functions K z i n x for i 0, . . ., k. Ω K has the same eigenvalues than S k .Knowing that S k is symmetric, its eigenvalues z i and eigenvectors Λ s are obtained conventionally using the algorithm QR.Then the eigenvectors of Ω K are given by Λ D −1 Λ s .
The computation of the three other functions L z i n x , M z i n x , and N z i n x is in the same manner.
Each of the functions K z i n x , L z i n x , M z i n x , and N z i n x admits exactly i zeros in the interval a, b .
For i ≥ k/2, the value of z i corresponding to K z i n x is equal to z i−1 corresponding to x , and the values of z i corresponding to L z i n x and M z i n x are identical.
However, the computation shows a significant numerical instability if the value of n increases beyond 10.This is due to the computation of the coefficients a j in the expression 3.21 .
A technique was proposed by Dobner and Ritter 29 to stabilize such computation.They proposed to use another expression of the polynomial P n such as

3.30
This variation gives a more stable computation.
The product E r • E μ • E ν satisfies the potential equation 3.19 within the elliptic cone.However for the external space other kinds of solutions are needed which vanish at infinity.
There is another kind of function F z n x solution of 3.19 which vanishes at infinity.These functions are called Lamé's functions of second kind.They are defined by Hobson 16 : For each value of E z i n x we have a function F z i n x .We have then 2k 1 Lamé's functions of second kind.
Recalling that surface of the elliptic cone is given by the relation μ μ 0 , the spatial velocity potential flow Φ, solution of 3.13 , is then 16 where A n,i are constants to be defined using the boundary conditions.

Boundary Condition
In addition to the relations 3.1 -3.3 the velocity potential flow Φ verifies a kinematical boundary condition on the surface of the elliptic cone such as

3.36
Let us denote k a/b.And suppose that 0 ≺ α ≺ K and 0 ≺ γ ≺ 4K with K is the complete elliptic integral F k, π/2 .For a function of α, β, and γ on the boundary surface μ μ 0 we will have

3.37
We can then solve the problem for the potential Φ in a given point in the space with the boundary conditions defined on the surface of the airship μ μ 0 .The relations 3.32 and 3.35 give Using 3.37 and 3.38 we obtain A n,i and then Φ and by the way M aij .
We can now deduce the different components of the mass matrix M. Finally, the vector of gyroscopic forces Q G can then be expressed in a developed form as 3.39

Journal of Applied Mathematics
This leads to the developed dynamic model: 3.40b

Results and Discussion
In this section we present the computation results of the added mass matrix of the MC500.Since any comparison with a similarly shaped airship is impossible at present, we conducted the calculation by the geometric method that we present here briefly in order to compare the results of our method.Geometric method is well known and discussed extensively elsewhere 6, 8, 12, 20 .In a 2D analysis the planar coefficients m ij are established for the standard shapes rectangle, where m 11 is a 2D added mass coefficient for the forward motion.
According to the large difference of size between the diagonal and off-diagonal terms, we will neglect these last terms, keeping only the term M a46 .

Results
We present here some characteristics of the geometry of the airship: z G 0.5 m; a 2.5 m; b 6.5 m; c 2 m; b 1 5.4 m; b 3 6.5 m; volume V 500 m 3 .Numerical results are presented in Table 2.

Discussion
The first results described here show that thanks to the application of the velocity potential flow theory to this unconventional airship it is possible to obtain reasonable values of the added masses matrix.To our knowledge, this is the first attempt to compute the added masses using this technique for an elliptic cone-shaped airship.Some differences can be observed between the two methods to certain terms of the added masses matrix.Experimental studies on the prototype nearing completion will confirm the accuracy of our method.
Although some geometric assumptions are made, it nonetheless demonstrates the capability of this method to compute interesting values of the added masses matrix.
Validation of our technique will allow a good estimation of the added masses matrix in preprocessing phase for this kind of airship.The development of an airship is usually done by iterative techniques model, stability, propulsion, sizing of the control, rudders, etc. ; obtaining a fairly accurate aerodynamic model early design allows better refinement of the final model the airship.

Conclusion
In this paper a dynamic model of an unconventional airship is presented.The original shape of the careen induces a necessary reformulation of the dynamic and aerodynamic study of these flying objects.A special focus is put on the computation of the added masses.According with the original shape, an assumption of elliptic cone was made to define the careen.The ellipsoidal harmonics and the Lamé's equations are revisited for the analysis of the velocity potential flow, according to the constraints of fluid-structure interaction.The first results seem promising.

Nomenclature η 1
x 0 , y 0 , z 0 T : Vector position of the origin expressed in the fixed reference frame R 0 η 2 φ, θ, ψ T : Vector orientation of the pointer R m in regards to R and given by the Euler angles η η 1 , η 2 T : Vector attitude compared to R 0 η : Velocity vector compared to R 0 expressed in R 0 ν 1 u, v, w T : Velocity vector expressed in R m ν 2 p, q, r T : Vector of angular velocities expressed in R m ν ν 1 , ν 2 T : The6× 1 velocity vector m : The mass of the airship I 3 : The identity matrix 3 × 3.

Figure 2 :
Figure 2: Position of the rotors.

Figure 4 :
Figure 4: Representation of the spheroconal coordinates in the elliptic cone.

Figure 5 :
Figure 5: Influence of the parameter μ on the shape of the cone.

Table 1 :
Characteristics of lamé's functions of first kind.

Table 2 :
Comparison between the two methods of computation., etc. .Following the exhaustive study ofBrenner 8, we model our flying wing as a truncated cone T with elliptic section.The airship is divided into a dozen cross-sections to optimize the inclusion of changes in transverse dimensions in the 3D calculation.The computing of the terms of the added masses matrix can be seen, for example, as