Computation of the Added Masses of an Unconventional Airship

Abstract

This paper presents a modelling of an unmanned airship. We are studying a quadrotor flying wing. 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 from mathematical 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.

1. 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.

Figure 1 

The flying wing Airship MC500.

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 air-structure 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.

2. Dynamic Model

2.1. Kinematics

The airship MC500 is assumed to be a rigid flying object.

We use two reference frames. First an earth-fixed frame 𝑅₀=(𝑂,𝑋_0,𝑌₀,𝑍₀) assumed to be Galilean. Then a local reference frame 𝑅_𝑚=(𝐺,𝑋_𝑚,𝑌_𝑚,𝑍_𝑚) fixed to the airship. Its axes are selected as follows:𝑋_𝑚 and 𝑌_𝑚 are the transverse axis of the airship, 𝑍_𝑚: 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 𝑅_𝑚 and the global frame is given by Goldstein [19]: 𝐽^𝑇₁𝜂₂=𝑐𝜃⋅𝑐𝜓𝑐𝜃⋅𝑠𝜓−𝑠𝜃𝑠𝜙⋅𝑠𝜃⋅𝑐𝜓−𝑐𝜙⋅𝑠𝜓𝑠𝜙⋅𝑠𝜃⋅𝑠𝜓+𝑐𝜙⋅𝑐𝜓𝑠𝜙⋅𝑐𝜃𝑐𝜙⋅𝑠𝜃⋅𝑐𝜓+𝑠𝜙⋅𝑠𝜓𝑐𝜙⋅𝑠𝜃⋅𝑠𝜓−𝑠𝜙⋅𝑐𝜓𝑐𝜙⋅𝑐𝜃, (2.1) such as 𝐽^𝑇₁(𝜂₂)⋅𝐽₁(𝜂₂)=𝐽₁(𝜂₂)⋅𝐽^𝑇₁(𝜂₂)=𝐼₃. We denote by 𝑐𝜃=cos𝜃; 𝑠𝜙=sin𝜙.

Using the rotation matrix 𝐽₁(𝜂₂), the expression of the linear speed in the reference frame 𝑅₀ is given by 𝜂̇₁=𝐽₁𝜂₂⋅𝜈₁. (2.2) The angular speed is expressed as follows: 𝜂̇₂=𝐽₂𝜂₂⋅𝜈₂ (2.3) with 𝐽₂𝜂₂=01𝑠𝜙tan𝜃𝑐𝜙tan𝜃0𝑐𝜙−𝑠𝜙𝑠𝜙𝑐𝜃𝑐𝜙𝑐𝜃. (2.4) It is noticed that the parameterization by the Euler angles has a singularity in Θ=𝜋/2+𝑘𝜋.

This parameterization is acceptable because it is not possible for an airship to reach this singular orientation of 90 degrees pitching angle.

2.2. 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 𝑇 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: 1𝑇=2𝜈^𝑇𝑀_𝑏+𝑀_𝑎𝑀𝜈=𝑇_𝑏+𝑇_𝑎, (2.5) where 𝑀_𝑎 is the added mass matrix due to the motion of the surrounding air and 𝑀_𝑏 is the mass matrix of the body.

For an airship fully immersed in the air, the 6×6 added mass matrix 𝑀_𝑎 is assumed to be positive and definite. As the added kinetic energy 𝑇_𝑎, it will be discussed in the next section. The whole mass matrix 𝑀 is assumed symmetric block-diagonal matrix: 𝑀𝑀=^𝑇𝑇00𝑀_𝑅𝑅. (2.6) For the computation of the whole dynamics model, we choose to use the Kirchoff’s equation [20]: 𝑑𝑑𝑡𝜕𝑇𝜕𝜈₁+𝜈₂∧𝜕𝑇𝜕𝜈₁=𝜏₁,𝑑𝑑𝑡𝜕𝑇𝜕𝜈₂+𝜈₂∧𝜕𝑇𝜕𝜈₂+𝜈₁∧𝜕𝑇𝜕𝜈₁=𝜏₂, (2.7) where ∧ is the vectorial product, 𝜏₁ and 𝜏₂ are, respectively, the external forces and torques, including the rotors effects, the weight (mg), the buoyancy 𝐵, and the aerodynamic lift (𝐹_𝐿) and drag (𝐹_𝐷).

The dynamical system of the airship becomes [21–23] 𝑀_𝑇𝑇00𝑀_𝑅𝑅𝜈̇₁𝜈̇₂=𝜏₁−𝜈₂∧𝑀_𝑇𝑇𝜈₁𝜏₂−𝜈₂∧𝑀_𝑅𝑅𝜈₂−𝜈₁∧𝑀_𝑇𝑇𝜈₁. (2.8)

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 𝛽_𝑖 around the 𝑌 axis (−180^∘≤𝛽_𝑖≤180^∘) and a rotation of angle 𝛾_𝑖 around an axis 𝑍_𝑖𝑅 normal to 𝑌 and initially coinciding with the 𝑍 axis (−30^∘≤𝛾_𝑖≤30^∘). A fictive axis 𝑋_𝑖𝑅 completes the rotor frame.

Figure 2 

Position of the rotors.

Let us denote 𝑃_𝑖 by the position of the rotor 𝑖. We can then define a rotation matrix 𝐽^𝑖₃ between the frame (𝑃_𝑖,𝑋_𝑖𝑅,𝑌,𝑍_𝑖𝑅) and the local frame 𝑅_𝑚 such as 𝐽^𝑖₃=𝑐𝛾_𝑖𝑐𝛽_𝑖−𝑠𝛾_𝑖𝑐𝛽_𝑖𝑠𝛽_𝑖𝑠𝛾_𝑖𝑐𝛾_𝑖0−𝑐𝛾_𝑖𝑠𝛽_𝑖𝑠𝛾_𝑖𝑠𝛽_𝑖𝑐𝛽_𝑖. (2.9) If we suppose that the intensity of the thrust force of the rotor 𝑖 is ‖𝐅_𝑖‖, this force could be introduced in the second member of the dynamic equation as 𝐅_𝑖=𝐽^𝑖₃𝐹_𝑖⋅𝑒_𝑋𝑚, (2.10) where 𝑒_𝑋𝑚 is a unitary vector along the 𝑋_𝑚 axis.

The torque produced by this rotor in the centre of inertia 𝐺 is 𝐅_𝑖∧𝑃_𝑖𝐺.

Weight and Buoyancy
An important characteristic of the airships is the buoyancy 𝐵_𝑢. This force represents a natural static lift, corresponding roughly to 1 Kg for each m³ of helium involved in the careen. We suppose here that this force is applied in the centre of buoyancy 𝐵 different from the centre of inertia 𝐺: 𝐵_𝑢=𝜌_air⋅𝑉⋅𝑔, (2.11) where 𝑉 is the volume of the careen, 𝜌_air is the density of the air, and 𝑔 the gravity.

Let us note 𝐹_𝑊𝐵 and 𝑀_𝑊𝐵 the force and the moment due to the weight and buoyancy.

We have 𝐹_𝑊𝐵=𝑚𝑔−𝐵_𝑢⋅𝐽^𝑇₁𝑒_𝑍,𝑀_𝑊𝐵=𝐵_𝑢⋅𝐽^𝑇₁⋅𝑒_𝑍∧𝐵𝐺. (2.12)

Aerodynamic Forces 𝐹_𝐴
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.13) with 𝜏𝜏=¹𝜏₂, (2.14) such as 𝜏₁=⁴_𝑖∑₌₁𝐹_𝑖⋅𝑐𝛾_𝑖⋅𝑐𝛽_𝑖−𝑚𝑔−𝐵_𝑢⋅𝑠𝜃_4𝑖∑₌₁𝐹_𝑖𝑠𝛾_𝑖+𝑚𝑔−𝐵_𝑢−⋅𝑠𝜙⋅𝑐𝜃_4𝑖∑₌₁𝐹_𝑖𝑐𝛾_𝑖⋅𝑠𝛽_𝑖+𝑚𝑔−𝐵_𝑢,𝜏𝑐𝜙⋅𝑐𝜃₂𝑐=−⁴_𝑖∑₌₁𝐹_𝑖⋅𝑠𝛾_𝑖+𝑏₁𝐹₁𝑐𝛾₁⋅𝑠𝛽₁−𝐹₂𝑐𝛾₂⋅𝑠𝛽₂+𝑏₃𝐹₃𝑐𝛾₃⋅𝑠𝛽₃−𝐹₄𝑐𝛾₄⋅𝑠𝛽₄+𝐵_𝑢𝑧_𝐵⋅𝑠𝜙⋅𝑐𝜃−𝑐∑⁴_𝑖=1𝐹_𝑖𝑐𝛾_𝑖⋅𝑐𝛽_𝑖+𝑎𝐹₄𝑐𝛾₄⋅𝑠𝛽₄+𝐹₃𝑐𝛾₃⋅𝑠𝛽₃−𝐹₁𝑐𝛾₁⋅𝑠𝛽₁−𝐹₂𝑐𝛾₂⋅𝑠𝛽₂+𝐵_𝑢𝑧_𝐵𝑏⋅𝑠𝜃₁𝐹₁𝑐𝛾₁⋅𝑐𝛽₁−𝐹₂𝑐𝛾₂⋅𝑐𝛽₂+𝑏₃𝐹₃𝑐𝛾₃⋅𝑐𝛽₃−𝐹₄𝑐𝛾₄⋅𝑐𝛽₄+𝑎𝐹₄𝑠𝛾₄+𝐹₃𝑠𝛾₃−𝐹₁𝑠𝛾₁−𝐹₂𝑠𝛾₂,Q_𝐺=−𝜈₂∧𝑀_𝑇𝑇𝜈₁−𝜈₂∧𝑀_𝑅𝑅𝜈₂−𝜈₁∧𝑀_𝑇𝑇𝜈₁. (2.15)

3. Computation of the Added Masses

3.1. 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 𝑥-𝑧 axis. Each transverse section of the careen parallel to the plan 𝑦-𝑧 gives roughly an ellipse. This motivated us to model as first assumption the airship as an elliptic cone (Figure 3).

Figure 3 

Transverse sections of the MC500.

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 𝜌_𝑓, that is, an incompressible fluid.(b)A velocity potential Φ_𝑓 exists and satisfies the Laplace’s equation throughout the fluid domain: ∇²Φ_𝑓=0, (3.1) 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: Φ_𝑓∞→0. (3.2)(c) The velocity of the fluid obeys the equation: 𝑉_𝑓=∇Φ_𝑓. (3.3)

In his study [6], Lamb proves that the kinetic energy 𝑇_𝑎 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 𝜈=(𝑢,𝑣,𝑤,𝑝,𝑞,𝑟)^𝑇 and can also be expressed as function of the velocity potential of this fluid by the following relation: 𝑇_𝑎=12𝜈^𝑇𝑀_𝑎1𝜈=−2𝜌_𝑓𝜕𝑆Φ_𝑓𝜕Φ_𝑓𝜕𝑛_𝑂𝑑𝑆. (3.4)𝑛_𝑂 is an outward normal vector. Φ_𝑓 can be splitted on Φ_𝑓=𝜈^𝑇Φ=𝑢Φ₁+𝑣Φ₂+𝑤Φ₃+𝑝Φ₄+𝑞Φ₅+𝑟Φ₆, (3.5) where Φ appears as a collection of six spatial potential Φ_i functions of 𝑥, 𝑦, 𝑧 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 𝑀_𝑎 in function of the velocity potential flow of the air surrounding the airship Φ such as 𝑀_𝑎𝑖𝑗=−𝜌_𝑓𝑆Φ_𝑖𝜕Φ_𝑗𝜕𝑛_𝑂𝑑𝑆. (3.6) 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. Θ=Θ₀ represents the external surface of the elliptic cone.

Figure 4 

Representation of the spheroconal coordinates in 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 𝑥²𝜇²+𝑦²𝜇²−𝑎²+𝑧²𝜇²−𝑏²=0. (3.7)𝑎, 𝑏, and 𝜇₀ are parameters that define the elliptic cone.

With 𝑏²≻𝜁²≻𝑎²≻𝜇²≻𝜇²₀, (3.8)𝜌 is always the distance of a point to the origin. We define then the Cartesian coordinates (𝑥, 𝑦, 𝑧) as 𝑥=𝜌𝜇𝜁,𝜌𝑎𝑏𝑦=𝑎𝜇²−𝑎²𝜁²−𝑎²𝑎²−𝑏²,𝜌𝑧=𝑏𝜇²−𝑏²𝜁²−𝑏²𝑏²−𝑎². (3.9) By analogy with the azimuthal angle in the previous parameterisation, we can introduce an angle 𝜑 defined by cos²𝑏𝜑=²−𝜁²𝑏²−𝑎²,sin²𝜁𝜑=²−𝑎²𝑏²−𝑎²,suchthat𝜁=𝑎²cos²𝜑+𝑏²sin²𝜑, (3.10) with (0≤𝜑≤2𝜋).

The surface of the elliptic cone is now defined by 𝜇=𝜇₀. By varying this parameter, one can see its influence on the shape of the cone in Figure 5.

Figure 5 

Influence of the parameter 𝜇 on the shape of the cone.

To express the components of the fluid velocity 𝑉_𝑓 in the conical elliptic reference frame, we define three unitary vectors ^→𝑒_𝜌, ^→𝑒_𝜑, ^→𝑒_𝜇. With ^→𝑒_𝜌=1ℎ_𝜌𝜕𝑥𝜕𝜌𝜕𝑦𝜕𝜌𝜕𝑧,𝜕𝜌^→𝑒_𝜑=1ℎ_𝜑𝜕𝑥𝜕𝜑𝜕𝑦𝜕𝜑𝜕𝑧,𝜕𝜑^→𝑒_𝜇=1ℎ_𝜇𝜕𝑥𝜕𝜇𝜕𝑦𝜕𝜇𝜕𝑧𝜕𝜇, (3.11)ℎ_𝑖 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 𝜇=𝜇₀; ^→𝑒_𝜇 appears as the normal vector to the surface in a given point.

3.2. Laplace’s Equation

With the spheroconal coordinates, the Laplace’s equation (3.1) becomes [16] 𝜇²−𝜁²𝜕²Φ𝜕𝛼²+𝜌²−𝜁²𝜕²Φ𝜕𝛽²+𝜌²−𝜇²𝜕²Φ𝜕𝛾²=0 (3.13) with 𝛼=𝑏^𝜌_𝑏𝑑𝑡𝑡²−𝑎²𝑡²−𝑏²,𝛽=𝑏^𝜇_𝑎𝑑𝑡𝑏²−𝑡²𝑡²−𝑎²,𝛾=𝑏^𝜁₀𝑑𝑡𝑎²−𝑡²𝑏²−𝑡². (3.14) To solve the Laplace’s equation (3.13), we use the traditional separation of variables: Φ𝜌,𝜇,𝜁=𝑅𝜌⋅𝑍𝜇⋅𝑌(𝜁). (3.15) Replacing (3.16) in (3.12) gives 𝜇²−𝜁²𝑅𝑑²𝑅𝑑𝛼²+𝜌²−𝜁²𝑍𝑑²𝑍𝑑𝛽²+𝜌²−𝜇²𝑌𝑑²𝑌𝑑𝛾²=0. (3.16) This equation can be splitted on three separated differential equations: 𝑑²𝑅𝑑𝛼²−𝑛(𝑛+1)𝜌²−𝑎²+𝑏²𝑠𝑑𝑅=0,₂𝑍𝑑𝛽²+𝑛(𝑛+1)𝜇²−𝑎²+𝑏²𝑠𝑑𝑍=0,₂𝑌𝑑𝛾²−𝑛(𝑛+1)𝜁²−𝑎²+𝑏²𝑠𝑌=0. (3.17) And by replacing 𝛼, 𝛽, and 𝛾 by their expressions in , one can obtain 𝜌²−𝑎²𝜌²−𝑏²𝑑²𝑅𝑑𝜌²+𝑟2𝜌²−𝑎²−𝑏²𝑑𝑅−𝑛()𝜌𝑑𝜌𝑛+1²−𝑎²+𝑏²𝑠𝜇𝑅=0,₂−𝑎²𝜇²−𝑏²𝑑²𝑍𝑑𝜇²+𝜇2𝜇²−𝑎²−𝑏²𝑑𝑍−𝑛()𝜇𝑑𝜇𝑛+1²−𝑎²+𝑏²𝑠𝜁𝑍=0,₂−𝑎²𝜁²−𝑏²𝑑²𝑌𝑑𝜁²+𝜁2𝜁²−𝑎²−𝑏²𝑑𝑌−𝑛()𝜁𝑑𝜁𝑛+1²−𝑎²+𝑏²𝑠𝑌=0. (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 𝑥²−𝑎²𝑥²−𝑏²𝑑²𝐸_𝑛𝑑𝑟²+𝑥2𝑥²−𝑎²−𝑏²𝑑𝐸_𝑛−𝑛()𝑥𝑑𝑥𝑛+1²−𝑎²+𝑏²𝑠𝐸_𝑛=0. (3.19) For a given 𝑛 integer, we can find 𝑧=(𝑎²+𝑏²)⋅𝑠 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 𝑃^𝑧_𝑛(𝑥)=^𝑚_𝑗=0𝑎_𝑗𝑥^2𝑗isLamés′spolynomial, (3.21) where 𝑚 depends on the integer 𝑘, as shown in Table 1: 𝑛𝑘=2,if𝑛iseven,𝑛−12,if𝑛isodd. (3.22)𝜓_𝑛(𝑥) is called the principal product.

Table 1 

Characteristics of lam’s functions of first kind.

There are four different Lamé’s functions that differ from their characteristics. We present them in Table 1.

𝑧_𝑖 is an eigenvalue. We remark that 𝑧_𝑖 corresponds to a value of the parameter 𝑧 in (3.19), which determines as many equations as there are values of 𝑧_𝑖.

For a given value of 𝑛, the unknowns are the coefficients 𝑎_𝑗 (3.23) and the parameters 𝑧_𝑖. We consider as a first step the functions 𝐾^𝑧_𝑛(𝑥); 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: 2𝑘−𝑗+12𝑛−2𝑘+2𝑗−1_𝜆𝑗𝑎_𝑗−1+𝑎²+𝑏²𝑛−2𝑘+2𝑗𝑛−2𝑘+2𝑗−𝑧_𝜛𝑗𝑎_𝑗−𝑎²𝑏²𝑛−2𝑘+2𝑗+2𝑛−2𝑘+2𝑗+1_𝜎𝑗𝑎_𝑗+1=0, (3.23) or 𝜆_𝑗𝑎_𝑗−1+𝜛_𝑗−𝑧𝑎_𝑗+𝜎_𝑗𝑎_𝑗+1=0for𝑗=0,…,𝑘. (3.24) Knowing that 𝑎_𝑘+1=0, the iteration will stop at rank 𝑘. And if we introduce the conditions: 𝑎₋₁=𝑎_𝑘+1=0, (3.25) the whole recurrence’s relations (3.28) can be written in a matrix form such as 𝜛₀𝜎₀𝜆000⋯0₁𝜛₁𝜎₁00⋯00𝜆₂𝜛₂𝜎₂0⋯000𝜆₃𝜛₃𝜎₃⋮⋮⋮⋱000⋯0𝜆_𝑘−1𝜛_𝑘−1𝜎_𝑘−100⋯00𝜆_𝑘𝜛_𝑘Ω𝐾⋅𝑎₀𝑎₁⋮⋮⋮𝑎_𝑘−1𝑎_𝑘Λ𝑎=𝑧⋅⁰𝑎₁⋮⋮⋮𝑎_𝑘−1𝑎_𝑘Λ, (3.26) or Ω_𝐾Λ=𝑧⋅Λ. (3.27)Ω_𝐾 is a square matrix of dimension (𝑘+1). The vector Λ is an eigenvector associated to the eigenvalue 𝑧.

There are basically (2𝑘+1) eigenvalues and eigenvectors. In fact we show that it is possible to find a diagonal matrix 𝐷 and a symmetric matrix 𝑆_𝐾 such as 𝑆_𝐾=𝐷⋅Ω_𝐾⋅𝐷⁻¹, (3.28) with 𝑐𝐷=⁰00⋯00𝑐₁0⋯⋮00𝑐₂⋯⋮⋱000⋯0𝑐_𝑘𝑐,where⁰𝑐=1,_𝑗=𝜎_𝑗−1𝜆_𝑗⋅𝑐_𝑗−1for𝑗=1,…,𝑘. (3.29)

The matrix 𝑆_𝑘 is diagonalizable such as 𝑆_𝑘=𝑅^𝑇⋅Δ⋅𝑅, where Δ is diagonal and 𝑅 is an orthogonal matrix. Accordingly, Ω_𝐾 is diagonalizable and admits 2𝑘+1 separate eigenvalues 𝑧_𝑖 associated to 𝑘+1 eigenvectors Λ_𝑖, so 𝑘+1 functions 𝐾^𝑧𝑖_𝑛(𝑥) for 𝑖=0,…,𝑘.

Ω_𝐾 has the same eigenvalues than 𝑆_𝑘. Knowing that 𝑆_𝑘 is symmetric, its eigenvalues 𝑧_𝑖 and eigenvectors Λ_𝑠 are obtained conventionally using the algorithm QR. Then the eigenvectors of Ω_𝐾 are given by Λ=𝐷⁻¹Λ_𝑠.

The computation of the three other functions 𝐿^𝑧𝑖_𝑛(𝑥), 𝑀^𝑧𝑖_𝑛(𝑥), and 𝑁^𝑧𝑖_𝑛(𝑥) is in the same manner.

Each of the functions 𝐾^𝑧𝑖_𝑛(𝑥), 𝐿^𝑧𝑖_𝑛(𝑥), 𝑀^𝑧𝑖_𝑛(𝑥), and 𝑁^𝑧𝑖_𝑛(𝑥) admits exactly 𝑖 zeros in the interval ]𝑎,𝑏[.

For 𝑖≥𝑘/2, the value of 𝑧_𝑖 corresponding to 𝐾^𝑧𝑖_𝑛(𝑥) is equal to 𝑧_𝑖−1 corresponding to 𝑁^𝑧𝑖−1_𝑛(𝑥), and the values of 𝑧_𝑖 corresponding to 𝐿^𝑧𝑖_𝑛(𝑥) and 𝑀^𝑧𝑖_𝑛(𝑥) are identical.

However, the computation shows a significant numerical instability if the value of 𝑛 increases beyond 10. This is due to the computation of the coefficients 𝑎_𝑗 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 𝑃_𝑛 such as 𝑃^𝑧_𝑛(𝑥)=^𝑚_𝑗=0𝑎_𝑗𝑥1−²𝑎². (3.30) This variation gives a more stable computation.

The product 𝐸(𝑟)⋅𝐸(𝜇)⋅𝐸(𝜈) 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 𝐹^𝑧_𝑛(𝑥) solution of (3.19) which vanishes at infinity. These functions are called Lamé’s functions of second kind. They are defined by Hobson [16]: 𝐹_𝑛(𝑥)=(2𝑛+1)𝐸_𝑛(𝑥)⋅^∞_𝑥𝑑𝑡𝐸_𝑛(𝑥)²√𝑡²−𝑎²𝑡²−𝑏². (3.31) For each value of 𝐸^𝑧𝑖_𝑛(𝑥) we have a function 𝐹^𝑧𝑖_𝑛(𝑥). We have then (2𝑘+1) Lamé’s functions of second kind.

Recalling that surface of the elliptic cone is given by the relation 𝜇=𝜇₀, the spatial velocity potential flow Φ, solution of (3.13), is then [16] Φ=^∞_𝑛=0^2𝑛+1_𝑖=1𝐴_𝑛,𝑖𝐹^𝑧𝑖_𝑛𝜇𝐹^𝑧𝑖_𝑛𝜇₀𝐸^𝑧𝑖_𝑛𝜌⋅𝐸^𝑧𝑖_𝑛(𝜁), (3.32) where 𝐴_𝑛,𝑖 are constants to be defined using the boundary conditions.

3.3. 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.33)𝑛_𝑂 is an outward normal vector to the surface: 𝑛_𝑜=𝜕𝑥,𝜕𝜇𝜕𝑦,𝜕𝜇𝜕𝑧𝜕𝜇^𝑇. (3.34) The relation (3.33) can then be written as 𝑢𝜕𝑥𝜕𝜇+𝑣𝜕𝑦𝜕𝜇+𝑤𝜕𝑧(𝜁)=𝜕𝜇=𝜌⋅𝑔𝜕Φ_𝑓𝜕𝜇 (3.35) with 𝜕𝑥=𝜕𝜇𝜌𝜁,𝑎𝑏𝜕𝑦=𝜕𝜇𝜌𝜇𝑎𝜁²−𝑎²𝑎²−𝑏²𝜇²−𝑎²,𝜕𝑧=𝜕𝜇𝜌𝜇𝑏𝜁²−𝑏²𝑏²−𝑎²𝜇²−𝑏²,𝑔(𝜁)=𝑢⋅𝜁+𝑎𝑏𝑣⋅𝜇₀𝑎𝜁²−𝑎²𝑎²−𝑏²𝜇²₀−𝑎²+𝑤⋅𝜇₀𝑏𝜁²−𝑏²𝑏²−𝑎²𝜇²₀−𝑏². (3.36) Let us denote 𝑘=𝑎/𝑏. And suppose that 0≺𝛼≺𝐾 and 0≺𝛾≺4𝐾 with 𝐾 is the complete elliptic integral 𝐹(𝑘,𝜋/2).

For a function of 𝛼, 𝛽, and 𝛾 on the boundary surface 𝜇=𝜇₀ we will have 𝑓1𝛼,𝛽,𝛾𝑑𝑆=𝑏^2𝐾0𝑑𝛼₀^4𝐾𝑓𝛼,𝛽₀,𝛾𝜌²−𝜁²𝜌²−𝜇²𝜁²−𝜇²𝑑𝛾. (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 𝜇=𝜇₀.

The relations (3.32) and (3.35) give 𝜕Φ𝜕𝜇_𝜇=𝜇0=^∞_𝑛=0^2𝑛+1_𝑖=1𝐴_𝑛,𝑖𝐸^𝑧𝑖_𝑛𝜌𝐸^𝑧𝑖_𝑛(𝜁)𝐹̇^𝑧𝑖_𝑛𝜇₀𝐹^𝑧𝑖_𝑛𝜇₀=𝜌⋅𝑔(𝜁). (3.38) Using (3.37) and (3.38) we obtain 𝐴_𝑛,𝑖 and then Φ and by the way 𝑀_𝑎𝑖𝑗.

We can now deduce the different components of the mass matrix 𝑀.

Finally, the vector of gyroscopic forces 𝑄_𝐺 can then be expressed in a developed form as 𝑄_𝐺=𝑀₂₂𝑣𝑟−𝑀₃₃𝑀𝑞𝑤₃₃𝑝𝑤−𝑀₁₁𝑀𝑢𝑟₁₁𝑢𝑞−𝑀₂₂𝑣𝑝−𝑀₄₆𝑝𝑞+𝑀₅₅−𝑀₆₆𝑀𝑞𝑟₄₆𝑝²+𝑀₆₆−𝑀₄₄𝑝𝑟−𝑀₄₆𝑟²𝑀₄₄−𝑀₅₅𝑝𝑞+𝑀₄₆𝑞𝑟. (3.39) This leads to the developed dynamic model:𝑀₁₁⋅𝑢̇=⁴_𝑖=1𝐹_𝑖𝑐𝛾_𝑖⋅𝑐𝛽_𝑖−𝑚𝑔−𝐵_𝑢𝑠𝜃−𝑀₃₃𝑞⋅𝑤+𝑀₂₂𝑀𝑟⋅𝑣,₂₂⋅𝑣̇=⁴_𝑖=1𝐹_𝑖𝑠𝛾_𝑖+𝑚𝑔−𝐵_𝑢𝑠𝜙⋅𝑐𝜃+𝑀₃₃𝑝⋅𝑤−𝑀₁₁𝑀𝑟⋅𝑢,₃₃⋅𝑤̇=−⁴_𝑖=1𝐹_𝑖𝑐𝛾_𝑖⋅𝑠𝛽_𝑖+𝑚𝑔−𝐵_𝑢𝑐𝜙⋅𝑐𝜃+𝑀₁₁𝑢⋅𝑞−𝑀₂₂̇1𝑣⋅𝑝,𝑝=𝑀₄₄𝑀₆₆−𝑀²₄₆×−𝑀₆₆𝑐⁴_𝑖=1𝐹_𝑖𝑠𝛾_𝑖+𝑀₄₆−𝑀₆₆𝑏₁𝐹₁𝑐𝛾₁⋅𝑠𝛽₁−𝐹₂𝑐𝛾₂⋅𝑠𝛽₂+𝑀₄₆−𝑀₆₆𝑏₃𝐹₃𝑐𝛾₃⋅𝑠𝛽₃−𝐹₄𝑐𝛾₄⋅𝑠𝛽₄+𝑀₄₆𝑎𝐹₄𝑠𝛾₄+𝐹₃𝑠𝛾₃−𝐹₁𝑠𝛾₁−𝐹₂𝑠𝛾₂−𝑀₆₆𝐵_𝑢𝑧_𝐺𝑠𝜙⋅𝑐𝜃−𝑀₄₆𝑀₄₄−𝑀₅₅+𝑀₆₆𝑝𝑞+𝑀₅₅𝑀₆₆−𝑀²₄₆−𝑀²₆₆,𝑀𝑞𝑟₅₅⋅𝑞̇=−𝑐⁴_𝑖=1𝐹_𝑖𝑐𝛾_𝑖⋅𝑐𝛽_𝑖+𝑎𝐹₄𝑐𝛾₄⋅𝑠𝛽₄+𝐹₃𝑐𝛾₃⋅𝑠𝛽₃−𝐹₁𝑐𝛾₁⋅𝑠𝛽₁−𝐹₂𝑐𝛾₂⋅𝑠𝛽₂−𝐵_𝑢𝑧_𝐺⋅𝑠𝜃+𝑀₄₆𝑝²+𝑀₆₆−𝑀₄₄𝑝𝑟+𝑀₄₆𝑟²,̇1𝑟=𝑀₄₄𝑀₆₆−𝑀²₄₆×𝑀₄₆⋅𝑐⁴_𝑖=1𝐹_𝑖𝑠𝛾_𝑖+𝑀₄₆−𝑀₄₄𝑏₁𝐹₁𝑐𝛾₁⋅𝑠𝛽₁−𝐹₂𝑐𝛾₂⋅𝑠𝛽₂+𝑀₄₆−𝑀₄₄𝑏₃𝐹₃𝑐𝛾₃⋅𝑠𝛽₃−𝐹₄𝑐𝛾₄⋅𝑠𝛽₄−𝑀₄₄𝑎𝐹₄𝑠𝛾₄+𝐹₃𝑠𝛾₃−𝐹₁𝑠𝛾₁−𝐹₂𝑠𝛾₂+𝑀₄₆𝐵_𝑢𝑧_𝐺+𝑀𝑠𝜙⋅𝑐𝜃₂₄₄+𝑀²₄₆−𝑀₄₄𝑀₅₅𝑝𝑞+𝑀₄₆𝑀₄₄−𝑀₅₅+𝑀₆₆𝑞𝑟. (3.40a)(3.40b)

4. 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 𝑚_𝑖𝑗 are established for the standard shapes (rectangle, circle, ellipse, etc.). Following the exhaustive study of Brenner [8], we model our flying wing as a truncated cone (𝑇) 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 𝑀_𝑎11=_(𝑇)𝑚₁₁𝑦,𝑧𝑑𝑥, (4.1) where 𝑚₁₁ 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 𝑀_𝑎46.

4.1. Results

We present here some characteristics of the geometry of the airship: 𝑧_𝐺=0.5 m; 𝑎=2.5 m; 𝑏=6.5 m; 𝑐=2 m; 𝑏₁=5.4 m; 𝑏₃=6.5 m; volume 𝑉=500 m³. Numerical results are presented in Table 2.

Table 2 

Comparison between the two methods of computation.

4.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.

5. 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=[𝑥0,𝑦0,𝑧0]𝑇:	Vector position of the origin expressed in the fixed reference frame  𝑅0
𝜂2=[𝜙,𝜃,𝜓]𝑇:	Vector orientation of the pointer 𝑅𝑚 in regards to 𝑅 and given by the Euler angles
𝜂=[𝜂1,𝜂2]𝑇:	Vector attitude compared to 𝑅0
𝜂̇∶	Velocity vector compared to 𝑅0 expressed in  𝑅0
𝜈1=[𝑢,𝑣,𝑤]𝑇:	Velocity vector expressed in  𝑅𝑚
𝜈2=[𝑝,𝑞,𝑟]𝑇:	Vector of angular velocities expressed in  𝑅𝑚
𝜈=[𝜈1,𝜈2]𝑇:	The 6×1 velocity vector
𝑚∶	The mass of the airship
𝐼3∶	The identity matrix 3×3.