Igor Mozolevski
email: igor.mozolevski@ufsc.br

Edson Luiz Valmorbida
email: edsonvalmorbida@utfpr.edu.br

Efficient Equilibrated Flux Reconstruction in High Order Raviart -Thomas Space for Discontinuous Galerkin methods

We develop an efficient and computationally cheap method of equilibrated fluxes reconstruction for high -order dG solutions to elliptic problems using a specific computational basis in high order Raviart -Thomas space. The computational basis is designed in such a way that coordinates of equilibrated fluxes with respect to this basis can be easy calculated from the moments of the numerical fluxes of dG method. Some applications of this method in implementation of a posteriori error estimators for elliptic boundary value problems are considered.

Introduction

Equilibrated fluxes reconstruction in Raviart -Thomas space is used in finite element methods for development of fully computable (not involving unknown constants), efficient and reliable a posteriori error error estimates for elliptic, convectiondiffusion and parabolic problems, see e.g [START_REF] Kim | A posteriori error estimators for locally conservative methods of nonlinear elliptic problems[END_REF][START_REF] Ern | An accurate H(div) flux reconstruction for discontinuous Galerkin approximations of elliptic problems[END_REF][START_REF] Ern | Flux reconstruction and a posteriori error estimation for discontinuous Galerkin methods on general nonmatching grids[END_REF][START_REF] Creusé | A posteriori error estimator based on gradient recovery by averaging for discontinuous Galerkin methods[END_REF][START_REF] Creusé | A posteriori error estimator based on gradient recovery by averaging for convection-diffusion-reaction problems approximated by discontinuous Galerkin methods[END_REF][START_REF] Braess | An equilibrated a posteriori error estimator for the interior penalty discontinuous Galerkin method[END_REF][START_REF] Mozolevski | Goal-oriented error estimation based on equilibrated-flux reconstruction for finite element approximations of elliptic problems[END_REF][START_REF] Ern | Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations[END_REF]. As an another important application the equilibrated velocity recuperation from a discontinuous Galerkin solution to the Darcy equation in the multiphase flow in heterogeneous porous media should be mentioned, see [START_REF] Ern | Accurate velocity reconstruction for discontinuous Galerkin approximations of two-phase porous media flows[END_REF][START_REF] Erm | Discontinuous Galerkin method for two-component liquid gas porous media flows[END_REF]. One of the attractive properties of the equilibrated fluxes a posteriori error error estimates is the robustness in respect to the order of polinomial approximation (cf. [START_REF] Braess | An equilibrated a posteriori error estimator for the interior penalty discontinuous Galerkin method[END_REF], [START_REF] Ern | Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations[END_REF]), whereas the efficiency of residual type estimates can decrease with the degree (cf. [START_REF] Melenk | On residual-based a posteriori error estimation in hp-FEM[END_REF], [START_REF] Braess | Equilibrated residual error estimates are p-robust[END_REF]).

Owing to the local conservation properties, the discontinuous Galerkin (dG) finite element methods allow easy flux reconstruction in Raviart-Thomas space by locally prescription of the numerical flux moments as the degrees of freedom, cf. [START_REF] Di Pietro | Mathematical Aspects of Discontinuous Galerkin Methods[END_REF]. Such an approach offers cheap and efficient computational algorithm for implementation of flux reconstruction in lowest order Raviart-Thomas space.

Nevertheless, reconstruction of equilibrated fluxes from the prescribed moments of discrete numerical fluxes in higher order Raviart-Thomas space can be computationally involved procedure. In this article we introduce a specific modal basis in high order Raviart-Thomas space such that calculation of reconstructed flux coefficients in respect to this from the prescribed moments is extremely easy owing to orthogonal properties of edge elements of the basis. Using this tool we develop an efficient and computationally cheap method of equilibrated fluxes reconstruction from high-order dG solutions to elliptic problems. To demonstrate the potential of the method we consider an application to adaptive mesh refinement, where the method is used for equilibrated fluxes calculation needed for a posteriori error estimator.

Modal basis in high order Raviart-Thomas space

Let Ω be a polygonal domain in R 2 . Let us denote by H k (Ω ) the Sobolev space of order k ∈ N 0 . The space of vector functions u ∈ [L 2 (Ω )] 2 with weak divergence ∇ • u in L 2 (Ω ) is denoted by H(div , Ω ) . The reader is referred to [START_REF] Adams | Sobolev Spaces[END_REF] and [START_REF] Boffi | Mixed finite element methods and applications[END_REF] among others, where standard properties of the Sobolev and H(div ) spaces are exposed.

For discrete approximation of H(div ) spaces let us define in Ω a shape regular family T h of triangular meshes (see e.g. [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF]), where h = max T ∈ T h h(T ) denotes the mesh size and h(T ) is the diameter of the mesh element T . We denote the set of all mesh edges as E and decompose it in the set E i of all interior edges (interfaces between adjacent mesh elements) and the set of all boundary faces E ∂ . Next we define vector field n E : E → R 2 of edge normals, where n E (E) = n E is the fixed unit vector orthogonal to E which coincides with the external normal to ∂ Ω on the boundary edges. We also denote as n T the external normal to ∂ T for any mesh element T . For E ∈ E we denote by T E = {T ∈ T h : E ⊂ ∂ T } the set of all mesh elements sharing the edge E.

For any triangle T ∈ T h the local Raviart-Thomas space is defined by

RT k (T ) = [P k (T )] 2 + x y P k (T ), (1)
where P k (T ) denotes the space of polynomials in T of degree less on or equal to k ∈ N 0 . For u ∈ RT k (T ) the degrees of freedom are given by

∂ T (u • n T )p, ∀p ∈ P k (∂ T ); (2)
T u • q, ∀q ∈ [P k-1 (T )] 2 if k ≥ 1. (3)
Associated with the triangulation T h the global Raviart-Thomas finite element space is defined as

RT k (T h ) = {u h ∈ H(div , Ω ) | u h | T ∈ RT k (T ), ∀T ∈ T h }. (4)
The computational implementation of RT 0 elements is typically included in finite element software packages and it was carefully discussed in [START_REF] Bahriawati | Three matlab implementations of the lowest-order Raviart-Thomas MFEM with a posteriori error control[END_REF]. Here we aim at introducing a computational basis in high order RT k (T h ) such that the coordinates of an element in this basis can be easily calculated from its degrees of freedom ( 2) -(3). We start with the definition of the basis in the master element and then extend the definition to any T ∈ T h using Piola transformation.

Let us consider the reference triangle T with vertexes v1 = (-1, -1) , v2 = (1, -1) , v3 = (-1, 1) ;

T = {v 1 , v2 , v3 } = {(r, s) | r, s ≥ -1; r + s ≤ 0}.
For any

T ∈ T h , T = {v 1 , v 2 , v 3 }, v 1 = (x 1 , y 1 ) , v 2 = (x 2 , y 2 ) , v 3 = (x 3
, y 3 ) we fix the canonical affine application π T : T → T as:

π T (r, s) = - r + s 2 v 1 + r + 1 2 v 2 + s + 1 2 v 3 = (x(r, s), y(r, s)) . (5)
The Piola transformation corresponding to π T is defined for û ∈ [L 2 ( T )] 2 by

P T û(x, y) = 1 | det J T | J T û(π T (r, s)), (6)
where J T denotes the Jacobian matrix of π T .

Lemma 1 (Properties of Piola transformation, see e.g. [START_REF] Boffi | Mixed finite element methods and applications[END_REF]) For any u ∈ H(div , T ) and v ∈ H 1 (T ) we have

T (∇ • u)v = T ( ∇ • û) v; T u • ∇v = T û • ∇ v; ∂ T u • n T v = ∂ T û • n T v,
where û = P -1 T u and v = π -1 T v. Now we are ready to formulate the theorem that provides a construction of the basis.

Theorem 1 In the local Raviart-Thomas space RT k (T ), k ∈ N 0 , T ∈ T h , T = {v 1 , v 2 , v 3 } there exists a basis {Φ ∂ T i,l ,Ψ T m }, i = 1, 2, 3, l = 1, . . . , k + 1, m = 1, . . . , 2M, M = k(k+1) 2 such that (B1) Φ ∂ T i,l • n T i | E i = δ i,i L i l , i, i ∈ {1, 2, 3}, l ∈ {1, . . . , k + 1},
where E i denotes the triangle's edge opposite to the vertex v i , n T i is the unit normal to the edge E i external to T and

{L i l } k l=0 is orthonormal system of Legendre polynomials in L 2 (E i ) . (B2) Ψ T m , m = 1, . . . , 2M form basis in L 2 ([P k-1 ( T )] 2 )
and

Ψ T m • n i | E i = 0, i ∈ {1, 2, 3}, m ∈ {1, . . . , 2M};
Proof. Following [START_REF] Ervin | Computational bases for RT k and BDM k on triangles[END_REF] let us consider in the master element T vector functions

e 1 (r, s) = 1 2 r + 1 s + 1 , e 2 (r, s) = 1 2 r -1 s + 1 , e 3 (r, s) = 1 2 r + 1 s -1 , (7)
t 1 (r, s) = s + 1 2 r + 1 s -1 , t 2 (r, s) = r + 1 2 r -1 s + 1 . ( 8
)
We define

Φi,l (r, s) = L l-1 (s)e i (r, s), i = 1, 2, Φ3l (r, s) = L l-1 (r)e 3 (r, s), l ∈ {1, . . . , k + 1},
where {L n } k n=0 are (normalized) Legendre polynomials that form an orthonormal system in L 2 ([-1, 1]) . Since Φi,l satisfies the property

e 1 • n j | Ê j = √ 2 2 δ 1 j , e i • n j | Ê j = δ i j , i ∈ {2, 3}, j ∈ {1, 2, 3}
we obtain

Φ∂ T i,l • n T i | Êi = δ i,i L i l .

Note that

Êi

( Φ∂ T i,l • ni )( Φ∂ T i,l • ni ) = δ ll , l, l ∈ {1, . . . , k + 1}, i ∈ {1, 2, 3}. (9)
Next we define Ψm = pm (r, s)t 1 (r, s), m = 1, . . . , M and Ψm = pm-M (r, s)t 2 (r, s), m = M + 1, . . . , 2M, where polynomials pm form the orthonormal Dubiner basis in P k-1 ( T ), cf. [START_REF] Dubiner | Spectral methods on triangles and other domains[END_REF]. Since

t j (r, s) • n i | Êi = 0, j ∈ {1, 2}, i ∈ {1, 2, 3}, we have Ψm • n i | Êi = 0. Using the Piola transformation the respective basis in RT k (T ), T ∈ T h is defined as Φ ∂ T i,l = P • Φ il • π -1 T , Ψ T m = P • Ψm • π -1
T and the required properties (B1) -(B2) follow directly from Lemma 1.

Next we will demonstrate how to recuperate the coefficients in respect to this basis from the degrees of freedom of a finite element in RT k space.

Let us consider in

[P k-1 (T )] 2 a basis P T m 2M m=1 , P T m = p m 0 , m = 1, . . . , M; P T m = 0 p M-m j = M + 1, . . . , 2M;
where polynomials p m form the Dubiner basis in P k-1 (T ).

Lemma 2 Assume that for

u h ∈ RT k (T h ) u T = u h | T = ∑ i ,l c ∂ T i ,l Φ ∂ T i ,l + ∑ m c T m Ψ T m ( 10
)
be local representation in respect to the basis

{Φ ∂ T i,l ,Ψ T m } in RT k (T ). Let µ ∂ T i,l (u T ) = E i (u T • n i )L i l , i ∈ {1, 2, 3}, l ∈ 1, . . . , k + 1; (11) µ T m (u T ) = T u T • P m , m ∈ {1, . . . , 2M}, (12)
be the degrees of freedom of u T . Then

c ∂ T i,l = µ ∂ T i,l (u T ), i ∈ {1, 2, 3}, l ∈ 1, . . . , k + 1 and c T = G -1 T F T , (13)
where

c T = (c T 1 , . . . , c T 2M ) , G T = T Ψ T i • P T j 2M×2M , F T = µ T m (u) -∑ i,l µ ∂ T i,l T Φ ∂ T i,l • P T m 1×2M
.

Proof. From the edge moments ( 11) of [START_REF] Dubiner | Spectral methods on triangles and other domains[END_REF] we have:

µ ∂ T i,l (u T ) = ∑ i ,l c ∂ T i ,l E i (Φ ∂ T i ,l • n i )L i l + ∑ m c T m E i (Ψ T m • n i )L i l = c ∂ T i,l E L l L l = c ∂ T i,l
owing to (B1), (B2) and orthogonality of the Legendre polynomials [START_REF] Di Pietro | Mathematical Aspects of Discontinuous Galerkin Methods[END_REF]. Similarly from the element's moments [START_REF] Ern | Accurate velocity reconstruction for discontinuous Galerkin approximations of two-phase porous media flows[END_REF] we obtain:

µ T m (u T ) = ∑ i ,l c ∂ T i ,l T Φ ∂ T i ,l • P T m + ∑ m c T m T Ψ T m • P T m , that is G T c T = F T .
Note 1 Let us note that Lemma 2 provides extremely cheap method for the flux recuperation from the moments: in each element of the mesh we only need to solve a small linear system.

Note 2 Since T Ψ T m •Ψ T m = |J T | T P T Ψ T m • P T Ψ T m = T (J T J T )/|J T | Ψm • Ψ T m , we im- mediately obtain c T m = µ T m (u T ) -∑ i,l µ ∂ T i,l (u T ) T Φ ∂ T i,l •Ψ T m , m ∈ 1, ..., 2M
for triangles where (J T J T )/|J T | = Id. This situation occurs for rectangular equilateral elements for example, so for such structured triangular meshes the flux reconstruction can be calculated directly from the moments and does not require a solution of the local systems.

Equilibrated flux reconstruction for discontinuous Galerkin method

Let us present an application of the computational basis introduced in previous section to equilibrated fluxes reconstruction from discrete gradient of discontinuous Galerkin approximation to a solution of elliptic boundary value problem.

We consider in Ω next model problem:

-∇ • (D∇u) = f in Ω , (14) u
= g on ∂ Ω .
Here the diffusion coefficient D > 0 is supposed to be constant in Ω , f ∈ L 2 (Ω ) and g ∈ H 3/2 (∂ Ω ).

For a shape regular family T h of triangular meshes in Ω we introduce the (discontinuous) finite element spaces V k h as:

V k h := {v h ∈ L 2 (Ω ) : v h | T ∈ P k (T ), ∀T ∈ T h }.
Symmetric version of the interior penalty dG method is formulated as:

find u h ∈ V k h such that B h (u h , v h ) = F(v h ), ∀v h ∈ V k h , (15)
where

B h (u h , v h ) = T h D∇ h u h • ∇ h v h - E {{n E • D∇ h u h }}[[v h ]] + E (-{{n E • D∇ h v h }} + γ E [[v h ]]) [[u h ]], F(v) = Ω f v + E ∂ (-{{n • D∇ h v h }} + γ E [[v h ]]) g.
Here we are using standard definition (see e.g. [START_REF] Di Pietro | Mathematical Aspects of Discontinuous Galerkin Methods[END_REF]) for discrete gradient, mean value and jump at the edges; the penalty parameter

γ E | E = 2.5D(k + 1) 2 h -1
E is considered to be reasonable ( [START_REF] Hesthaven | Nodal discontinuous Galerkin methods[END_REF]) for the stabilization.

For v ∈ H 1 (T h ) + V k h we define the energy norm associated with the dG method by

|||v h ||| dG =   D 1 2 ∇ h v h 2 L 2 (Ω ) + E γ E [[v h ]] 2   1 2 . ( 16
)
For the energy norm of the error in dG approximation the following a priori estimate holds true, see e.g. [START_REF] Houston | Discontinuous hp-finite element methods for advectiondiffusion problems[END_REF].

Theorem 2 Let u ∈ H k+1 (Ω ) be a weak solution to (14) and u h ∈ V k h be the dG finite element approximation of u. Then the estimate

|||u -u h ||| dG ≤ Ch k u H k+1 (Ω ) (17)
holds with a constant C > 0 independent of h.

We will term the potential the weak solution u ∈ H 1 (Ω ) to the problem ( 14) and the flux σ (u) = -D∇u ∈ H(div , Ω ). Similarly, the dG solution u h ∈ V k h to ( 15) is said to be discrete dG potential and the discrete dG flux is defined as

σ h (u) = -D∇ h u h . Definition 1. A vector field s h ∈ H(div , Ω ) is called equilibrated up to order l ∈ N 0 if div s h -π l h ( f ) = 0, where π l h : L 2 (Ω ) → V l h denotes the orthogonal projection operator.
Following [START_REF] Ern | An accurate H(div) flux reconstruction for discontinuous Galerkin approximations of elliptic problems[END_REF], let us consider a flux t k-1 h (u h ) ∈ RT k-1 (T h ) with the degrees of freedom of ( 11) and ( 12) locally prescribed by :

k ≥ 1 : µ ∂ T i,l (t k-1 h (u h )) = E i -{{n • D∇ h u h }} + γ E [[u h ]] g L i l , E i ∈ ∂ T, i ∈ {1, 2, 3}, l ∈ 1, . . . , k; (18)
k ≥ 2 : µ T m (t k-1 h (u h )) = - T D∇ h u h • P m (19)
+ ∑ i E i χ e (E i )D(P m • n i )[[u h ]] g , m = 1, . . . , 2M,
where

[[u h ]] g = [[u h ]]| E , for E ∈ E i , u h | E -g, for E ∈ E D .
Note that this definition is correct since the numerical fluxes of dG method are uniquely defined at the edges of the mesh. Moreover, we will refer to the following lemma.

Lemma 3 Let u be a weak solution to (14) and u h ∈ V k h be the dG finite element approximation of u. The flux t k-1 h (u h ), defined by (18) - [START_REF] Kim | A posteriori error estimators for locally conservative methods of nonlinear elliptic problems[END_REF], is equilibrated up to order k -1 and there exists a constant C > 0, independent of h, such that

D 1 2 ∇u -D -1 2 t k-1 h (u h ) L 2 (Ω ) ≤ C|||u -u h ||| dG . (20)
See [START_REF] Ern | An accurate H(div) flux reconstruction for discontinuous Galerkin approximations of elliptic problems[END_REF] for the proof and more details. For computational reconstruction of t k-1 h (u h ) from the prescribed moments ( 18) and ( 19) we use the algorithm presented in Lemma 2.

Numerical examples

To demonstrate a potential of the suggested algorithm we consider an application to adaptive mesh refinement in dG approximation of elliptic boundary value problems, where the equilibrated fluxes are used in a posteriori error estimator in the energy norm [START_REF] Ervin | Computational bases for RT k and BDM k on triangles[END_REF].

We consider the error estimator introduced in [15]

η 2 (u h ) = ∑ T ∈T h η 2 (T ) = ∑ T ∈T h ((η O (T ) + η ∇ (T )) 2 + η 2 H (T ))
where

η O (T ) = f -∇•t h (u h ) L 2 (T ) is the oscillations term, η ∇ (T ) = D -1 2 (σ h (u h )- t h (u h )) L 2 (T ) measures the deviation of the discrete flux σ h (u h ) from H(div, Ω ) and η H (T ) = D 1 2 (∇ h (u h ) -∇ h (u O h )) L 2 (T ) measures the deviation of u h from H 1 (Ω ). Here u O h ∈ H 1 (Ω )
is Oswald interpolator of u h and the equilibrated flux t k-1 h (u h ) is reconstructed by prescription from ( 18) - [START_REF] Kim | A posteriori error estimators for locally conservative methods of nonlinear elliptic problems[END_REF] in the computational basis

{Φ ∂ T i,l ,Ψ T m } in RT k (T )
. This type of estimator has proven to be reliable, efficient and robust with respect to polynomial order of approximation space, see [START_REF] Braess | An equilibrated a posteriori error estimator for the interior penalty discontinuous Galerkin method[END_REF], [START_REF] Ern | Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations[END_REF]. The quality of the error estimator η is assessed in terms of the effectivity index

I η = η(u h ) D 1 2 (∇u-∇ h u h ) L 2 (Ω )
evaluated on sequences of uniformly and adaptively refined meshes.

Test case 1: uniform mesh refinement

Firstly we confirm the order of approximation of the exact flux by the reconstructed flux for a smooth solution to elliptic problem. So let us consider the model problem [START_REF] Ern | Flux reconstruction and a posteriori error estimation for discontinuous Galerkin methods on general nonmatching grids[END_REF] in Ω = (0, 1) 2 with D = 1, the right-hand side and the homogeneous Dirichlet boundary condition corresponding to the exact solution u(x, y) = sin(πx) sin(πy). The model problem was solved numerically using the dG method of order k = 1, 2, 3, 4 on a sequence of meshes obtained by successive uniform bisection from the initial unstructured mesh of 48 elements.

Let us denote by e(u

h ) = D -1 2 (σ (u) -σ h (u h )) L 2 (Ω ) the error in dG approx- imation, by e(t h ) = D -1 2 (σ (u) -t h (u h )) L 2
(Ω the error in approximation of the exact flux by the reconstructed Raviart-Thomas flux and by e(∇

• t h ) = f -∇ • t h (u h ) L 2 (Ω )
the error in equilibration of the reconstructed flux. In Table 1 we show the errors and the convergence rates of the dG method and of the equilibrated flux, reconstructed from discrete solution using suggested computational basis. We observe that the dG method accurately approximate the exact solution to the problem and exhibits the optimal order of convergence k predicted in the Theorem 2. The error in approximation of the exact flux by reconstructed equilibrated flux is also of order k, which is the optimal order of approximation of the exact flux by dG method and the optimal order of the projection of the exact flux on the Raviart-Thomas space RT k-1 . Finally, the order of convergence of Inspired by propagation of saturation front in a two-fase flow in porous media, we consider here the model problem [START_REF] Ern | Flux reconstruction and a posteriori error estimation for discontinuous Galerkin methods on general nonmatching grids[END_REF] in Ω = (0, 1) 2 with the solution u(x, y) = x(x -1)y(y -1) arctan 60 (x -5/4) 2 + (y + 1/4) 2 ) -1 , that exhibits a steep front in the interior of the domain. This example is commonly used for testing adaptive refinement algorithms for elliptic equations, c.f. [START_REF] Mitchell | A collection of 2D elliptic problems for testing adaptive grid refinement algorithms[END_REF]. Here we also consider D = 1, and the right-hand side and the homogeneous Dirichlet boundary condition correspond to the exact solution.

∇ • t h → f is k + 1, that is the flux t h is equilibrated up to order k -1 .
Meshes are adapted from the same initial mesh, using a refinement strategy based on the method proposed by Dörfler, whereby the elements in a minimal set M ⊂ T h , such that ∑ T ∈M η(T ) ≥ θ ∑ T ∈T h η(T ), are refined. Elements are refined using the longest edge bisection technique and additional refinements of the mesh are considered in order to eliminate hanging nodes.

Figure 1 displays in the first column the energy norm of the error and the convergence order calculated for uniform and adaptive mesh refinements with θ = 0.25, θ = 0.5 and θ = 0.75 in the Dörfler marking as a function of DOF on a logarithmic scale for order of dG method k = 1 to k = 4. The second column shows the respective effectivity indices and the third column presents the adaptively refined meshes corresponding to the error ≈ 0.01 for θ = 0.25. We can see that for given order of dG method the energy norm of the errors are very close for all values of the parameter in the Dörfler marking and asymptotically exhibit the optimal convergence rates. The effectivity indices remain above 1 and are asymptotically close to 1 even the order of approximation increases. We observe also that the number of DOF necessary to achieve the same global approximation error decreases with increasing polynomial degree k.

Conclusions

A specific modal computational basis, in which the coordinates of the equilibrated fluxes can be easy calculated from the numerical fluxes of dG method, is designed for high order Raviart-Thomas space. Optimal convergence of equilibrated fluxes and the robustness of the reconstruction method in respect to the order of dG method are confirmed numerically. Adaptive mesh refinement, guided by the equilibrated error estimator calculated in this basis, exhibits robust effectivity index and provides final meshes with less DOF for the same error tolerance for higher orders of the dG method.

Fig. 1

 1 Fig. 1 The energy norm of error, the convergence order (left column) and the effectivity index (middle column) as a function of the degrees of freedom (DOF) on a logarithmic scale for various θ in the Dörfler marking and for various orders of dG method. Right column: adaptively refined meshes, corresponding to the error ≈ 0.01 in the energy norm, with DOF number N DOF = 104991 for k = 1, N DOF = 7692 for k = 2, N DOF = 4270 for k = 3 and N DOF = 3930 for k = 4.

Table 1

 1 Convergence rates o(e(u h )) and o(t h (u h )) for different orders of approximation calculated for the refinement of order N r . 10 -1 2.951 3.664 × 10 -1 2.605 6.013 × 10 -2 5.544 2 1.394 × 10 -1 0.995 1.879 × 10 -1 0.964 1.521 × 10 -21.983 3 6.984 × 10 -2 0.997 9.425 × 10 -2 0.995 3.813 × 10 -31.996 4 3.495 × 10 -2 0.999 4.711 × 10 -2 1.000 9.540 × 10 -41.999 5 1.749 × 10 -2 0.999 2.354 × 10 -2 1.001 2.385 × 10 -42.000 2 1 2.480 × 10 -2 5.058 3.184 × 10 -2 4.275 3.725 × 10 -3 8.005 2 6.227 × 10 -3 1.994 7.819 × 10 -3 2.026 4.685 × 10 -4 2.991 3 1.561 × 10 -3 1.996 1.944 × 10 -3 2.008 5.864 × 10 -5

	k N r	e(u h )	o(e(u h ))	e(t h )	o(e(t h ))	e(∇ • t h )	o(e(∇ • t h ))
	1	1 2.779 × 2.998
		4 3.910 × 10 -4 1.998 4.852 × 10 -4 2.003 7.332 × 10 -6	3.000
		5 9.783 × 10 -5 1.999 1.212 × 10 -4 2.001 9.165 × 10 -7	3.000
		1 1.011 × 10 -3 7.457 1.294 × 10 -3 6.806 2.743 × 10 -4	9.262
	3	2 1.294 × 10 -4 2.965 1.648 × 10 -4 2.973 1.697 × 10 -5 3 1.623 × 10 -5 2.995 2.050 × 10 -5 3.007 1.059 × 10 -6	4.015 4.003
		4 2.029 × 10 -6 3.000 2.547 × 10 -6 3.008 6.613 × 10 -8	4.001
		5 2.536 × 10 -7 3.001 3.172 × 10 -7 3.005 4.133 × 10 -9	4.000
		1 6.619 × 10 -5 8.177 8.910 × 10 -5 7.072 8.448 × 10 -6 11.581
	4	2 4.042 × 10 -6 4.033 5.414 × 10 -6 4.041 2.801 × 10 -7 3 2.514 × 10 -7 4.007 3.371 × 10 -7 4.005 8.856 × 10 -9	4.915 4.983
		4 1.570 × 10 -8 4.001 2.109 × 10 -8 3.999 2.775 × 10 -10	4.996
		5 9.812 × 10 -10 4.000 1.319 × 10 -9 3.999 1.161 × 10 -11	4.579
	4.2 Test case 2: adaptive mesh refinement

Acknowledgements The authors gratefully acknowledge the support by CNPq, Brazil, grant 477935/2013-3.


J C Thomas

Z Jiang

C Wielgosz Gem

Continuous and Finite Element Methods for the Vibrations of Inflatable Beams

Keywords: Continuous element, mass matrix, stiffness matrix, follower forces, inflatable tubes, natural modes

Inflatable structures are under increasing development in various domains. Their study is often carried out by using 3D membrane finite elements and for static loads. There is a lack of knowledge in dynamic conditions, especially for simple and accurate solutions for inflatable beams. This paper deals with the research on the natural frequencies of inflatable Timoshenko beams by an exact method: the continuous element method (CEM), and by the classical finite element method (FEM). The dynamic stiffness matrix D(ω) is here established for an inflatable beam; it depends on the natural frequency and also on the inflation pressure. The stiffness and mass matrixes used in the FEM are deduced from D(ω). Natural frequencies and natural modes of a simply supported beam are computed, and the accuracy of the CEM is checked by comparisons with the finite element method and also with experimental results.

INTRODUCTION

Inflatable structures are increasingly used nowadays. They find applications in well-known structures (tires, boats) and in more recent applications, mainly in the spatial industry (inflatable antennas, shields, houses) or in civil engineering (buildings, temporary hospitals, inflatable scales). Other industrial sectors are developing inflatable structures. Inflatable wings are aeronautical applications. The naval industry is interested in many different structures: boats, inflatable laths to rigidify the sails and inflatable landing stages. Environmental protection is concerned with inflatable dams. Medical inflatable prosthesis are being increasingly used.

The growing interest for inflatable structures comes from the large amount of advantages that they offer: they are light, easy to fold, quite easy to built, easy to deploy. In the spatial industry, the focus is put on their ability to attain large sizes for a small storage volume and to present low specific mass and an interesting reliability. For these applications, the inflation pressure is quite low. Earthly structures are mainly used for temporary applications, but durable ones are also built nowadays: autonomous or in association with metallic reinforcements (tensairity systems [START_REF] Pedretti | European Congress on Computationnal Methods in Applied Sciences and Engineering[END_REF]). Generally, they need greater pressure than for spatial applications and may become very strong structures if the pressure is large enough. All these developments require engineers to be able to calculate their mechanical behaviour.

The type of material used, the geometrical characteristics of the beams, and the internal pressure lead to the strength and stiffness of the inflatable structure. Several types of fabrics are used. They are often manufactured from high modulus fibres such as Kevlar, Vectran. For instance, Vectran was used for the Mars Pathfinder [START_REF] Cadogan | Development and evaluation of the mars pathfinder inflatable airbag landing system[END_REF]. Polyester fabrics with a PVC coating are more often used for inflatable structures in civil engineering: Ferrari pre-stressed fabric is one of them. Materials are robust and can be packed and deployed many times with minimal degradation in performance.

Inflatable structures are built from elementary parts like tubes, panels, cones or tori. Due to the geometry of these elements, the strength of materials is well adapted to the study of these structures in order to predict their behaviour. A beam approach gives useful tools for easy model assembly of longitudinal elements. Several authors have worked on static analysis. Comer and Levy [START_REF] Comer | Deflections of an inflated circular-cylindrical cantilever beam[END_REF] have proposed an Euler Bernoulli formulation for an inflatable cantilever beam. Main and al. [START_REF] Main | Load deflection behaviour of space-based inflatable Beams[END_REF] have realized experiments and compared their results to Comer and Levy's theory. In fact, the effect of the pressure isn't obtained explicitly by using the usual beam theory. Clearly, the deflections of fabric tubes are dependent on the internal pressure, on the fabric's constitutive law, on the initial geometry and on the external loads. Fichter [START_REF] Fichter | A theory for inflated thin-walled cylindrical beams[END_REF] used the minimization of the global potential energy with Timoshenko's beam assumptions and was the first to put forward a theory for the deflections including the pressure effects. Thomas [START_REF] Thomas | Deflections of highly inflated fabric tubes[END_REF] has constructed a model for the tubes, which lead the way to quite complicated formulas and to a nonsymmetrical stiffness matrix. Even if the comparison between theoretical and experimental results have given evidence of a good accuracy of these formulas, the choice has been made here to use the recent Le Van's approach [START_REF] Van | Bending and buckling of inflatable beams: some new theoretical results[END_REF] who obtained simpler static beam formulations, which leads to bending formulas where the inflation pressure appears also explicitly in the rigidities of the beam. In fact, Le Van's study is an improvement on Fichter's. Equations are written by the use of the virtual work principle in lagrangian variables with finite rotations. This approach has been extended in this paper by introducing the dynamic terms in order to predict the influence of the pressure effects on the natural frequencies.

There are only a few papers which deal with theoretical studies on the dynamics of inflatable beams with a strength of materials approach. Membrane or shell analysis is generally used. For instance, the vibrations of inflatable dams (rubber dams) have lead to several studies where the authors mainly use displacement finite element computations [START_REF] Mysore | Dynamic analysis of single-anchor inflatable dams[END_REF]. Recently, Jha and al. [START_REF] Jha | Importance of geometric nonlinearity and follower pressure load in the dynamic analysis of a gossamer structure[END_REF] have suggested a study focused on the vibrations of an inflatable torus made of kapton dedicated to the spatial industry, and the model used was based on a shell theory. Concerning the beam approach, Main and al [START_REF] Main | Dynamic analysis of space-based inflated beam structures[END_REF] put forward work dealing principally with the damping mechanisms; an Euler-Bernoulli theory was used again. In fact the experimental results clearly show the dependency of the frequencies on the pressure. The beam dynamic theory that will be established in this paper follows Le Van's work on static analysis. We have chosen to use the continuous element method to analytically obtain the dynamic stiffness matrix because it allows an exact vibration analysis of the problem. The most popular computational method used in structural dynamic is the finite element method. In this work, the stiffness matrix and the mass matrix used in the FEM have also been calculated.

The physical intuition is that increasing the pressure in a structure leads to an increase of the material properties. In order to take into account the effects of the inflation pressure and the shear deformation, a Timoshenko beam theory has been used in this paper. The internal pressure is the pre-stress effect: it allows pre-stressing the walls and inducts the initial geometry and stiffness. The second special feature comes from the fact that the force inducted locally follows the walls and will be denoted by the follower force effect. The equations of motion therefore have to be written on the current configuration, which explains the choice of the lagrangian approach. This means that the follower load effect due to the pressure can be taken into account properly, which is of importance for these kinds of structures.

The framework of this paper is the following: in the first part, the differential equations of motion are obtained by using the virtual work principle written in lagrangian variables following [START_REF] Van | Bending and buckling of inflatable beams: some new theoretical results[END_REF]. In the second part, the dynamic stiffness matrix of an inflatable beam is built by connecting the values of the loads and the displacements at the two ends of the beam. FEM stiffness and mass matrixes are deduced from the dynamic stiffness matrix. In the fourth part, the stiffness and mass matrix are directly derived from the FEM and the coherence with the continuous element analysis is verified. Then the natural frequencies are calculated for an inflatable simply, supported beam by using the CEM. Comparisons are finally made with experimental results and show the accuracy of the theory.

THE DIFFERENTIAL EQUATIONS OF MOTION

Inflatable structures are made of pressurized fabrics. The inflation pressure leads to a pre-stress state in the walls of the structure, which ensures its stiffness. The shear deformations of these thin-walled structures must not be overlooked and one must use Timoshenko's beam kinematics because the crosssection doesn't remain orthogonal to the neutral fibre (see Fig 2 in [6]). Fig. 1 shows a straight uniform Timoshenko beam model.

The use of the virtual work principle written in lagrangian variables leads to differential equations of motion where the effect of the inflation pressure appears. This principle is written as follows:

(1) U denotes the displacement field and U* the virtual displacement field. Π is the first Piola-Kirchhoff stress tensor. ρ 0 is the density in the reference configuration, Ω 0 is the volume of the element in the reference configuration and δΩ 0 its boundary. N is the outward normal unit in the reference configuration, and f 0 the body force per unit reference volume. In the following, details are expressed for the virtual work of the inertia terms. The principle of calculus and the corresponding results for the other virtual works (internal virtual stress work, external virtual dead loads work, external virtual work due to the pressure) have just been condensed. Details can be found in [START_REF] Van | Bending and buckling of inflatable beams: some new theoretical results[END_REF].

The virtual work of the inertia terms

The cross-section is assumed to have a rigid body displacement. In a right-handed rectangular Cartesian coordinate system, Fig. 1 shows an inflatable tube of length ᐉ 0 in the reference configuration and a current cross section on the reference and in its current configuration. The displacement of any point P 0 of the cross section comes from:

x y z G 0 P 0 Y Z . . . P U(X,t) V(X,t) θ(x,t) G X 0 curent configuration reference configuration
Ω ρ Ω Ω Ω Ω δ Ω Ω U U N Π U U f U Π U 0 ⋅ = ⋅ ⋅ + ⋅ + - ∈ ∀ ∫ ∫ ∫ ∫ (2)
Here, U(G 0 , t) = (U(X, t), V(X, t), 0) is the displacement of the centre of the current cross-section, and R corresponds to the rotation of the section.

The virtual displacement field is chosen as follows:

(

where U*(G 0 ) = (U*(X), V*(X), 0), Ω* = (0, 0, θ*(X)) and where the particularity comes from GP which corresponds to the final configuration. This gives in lagrangian variables:

The virtual work of the inertia terms is:

(

where [START_REF] Thomas | Deflections of highly inflated fabric tubes[END_REF] This leads to: [START_REF] Van | Bending and buckling of inflatable beams: some new theoretical results[END_REF] By introducing the second moment of area of the beam cross section about the z axis I, and considering that the first moment of area is nil, the virtual work of the inertia terms is:

(8) ∫ ∂ ∂ + ∂ ∂ + ∂ ∂ = 0 0 2 2 * 2 2 * 2 2 * 0 0 * acc dX t I t V V t U U S W θ θ ρ             ∫ ∫ ∂ ∂ + ∂ ∂ - ∂ ∂ - ∂ ∂ - ∂ ∂ + ∂ ∂ - ∂ ∂ - + ∂ ∂ + ∂ ∂ = 0 0 0 S * 2 2 2 * 2 * 2 2 * 2 2 * 2 * 2 2 * 2 2 2 2 * 2 2 * 0 * acc dX dS t Y V cos t V sin t sin t V U sin t U cos t cos t U Y t V V t      0 sin Y V cos Y U . 0 cos t Y sin t Y t V sin t Y cos t Y t U ) P ( t ) t , P ( θ θ θ - - ∂ ∂ - ∂ ∂ - ∂ ∂ ∂ ∂ + ∂ ∂ - ∂ ∂ = ∂ ∂ U U             0 * 2 2 0 * ac d t W 0 Ω ρ Ω U U ⋅ ∂ ∂ = ∫ Z cos Y sin Y Z Y 0 1 0 0 0 cos sin 0 sin cos : θ θ θ θ θ θ - = - = = 0 0 P G R GP         GP Ω U U ∧ + = * 0 * 0 * ) G ( ) P ( ( ) 0 ) cos 1 ( Y V sin Y U Z Y 0 0 0 0 0 1 cos sin 0 sin 1 cos 0 V U . ) t , G ( ) t , P ( θ θ θ - - - = - - - + = - + = + + = = 0 0 P G R 0 0 0 0 P G I R U GP G G G P P P U 0 0                  

Strain tensor, stress tensor, constitutive law

The Green strain tensor is given by E = 1 _ 2 (H + H T + H T • H), where H = gradU is the displacement gradient tensor. One obtains: (9) [START_REF] Main | Dynamic analysis of space-based inflated beam structures[END_REF] The matrix of the second Piola-Kirchhoff stress tensor Σ is assumed to be a "beam" matrix: [START_REF] Wittrick | A general algorithm for computing natural frequencies of elastic structures[END_REF] The following generalized stresses can therefore be introduced: [START_REF] Banerjee | Exact dynamic stiffness matrix for composite Timoshenko beams with applications[END_REF] which are respectively considered as the normal force, the bending momentum and the shear force acting on the cross-section in the reference configuration.

Different kinds of materials are commonly used to build inflatable structures. Thin isotropic membrane may be used (example of the kapton in the spatial industry). Composites made of fabric embedded in a matrix material are also used (temporary buildings are often made of polyester fibres and rubber coating). Even if this kind of membrane follows an orthotropic behaviour, the material is considered here as an hyperelastic isotropic Saint-Venant Kirchhoff material and the constitutive laws are: [START_REF] Banerjee | Development of dynamic stiffness matrix for free vibration analysis of spinning beams[END_REF] ∑ 0 XX and ∑ 0 XY are the initial stresses. They are induced by the initial inflation with inducts pre-stress effects in the walls of the structures. E and G are respectively the Young modulus and the shear modulus.

The virtual work of the internal forces

The virtual work of the internal forces is: [START_REF] Richards | An accurate method in structural vibration analysis[END_REF] Π is calculated by using the deformation gradient F: Π = FΣ On the assumption that the material is linear and elastic and by taking into account the existence of prestressing, we have: [START_REF] Wielgosz | An inflatable fabric beam finite element[END_REF] [

dX dX d X K X KE 2 1 S I N sin X V cos X U 1 M sin X V cos X U 1 T X cos X V X sin X U 1 M dX dV cos T X sin M X V N dX dU sin T X cos M X U 1 N W * 0 2 0 0 * * 0 * * int 0 ∂ ∂ + ∂ ∂ + + ∂ ∂ + ∂ ∂ + + ∂ ∂ + ∂ ∂ + - ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂ + - + + ∂ ∂ + ∂ ∂ + - ∂ ∂ + ∂ θ θ θ                                                                                  0 * T * int d : ) ( W 0 Ω Ω gradU Π ∫ - = XY 0 XY XY XX 0 XX XX E . G 2 E . E + = + = Σ Σ Σ Σ ∫ = 0 S 0 XX dS N Σ , ∫ - = 0 S 0 XX dS y M Σ and ∫ = 0 S 0 XY dS T Σ [ ]= 0 0 0 0 0 0 XY XY XX Σ Σ Σ ∑           ∂ ∂ + - ∂ ∂ = θ θ sin X U 1 cos X V 2 1 E XY and 0 E YY =             ∂ ∂ ∂ ∂ - ∂ ∂ ∂ ∂ - ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ - ∂ ∂ = X sin X V Y 2 X cos X U Y 2 X Y X V X U 2 1 X cos Y X U E θ θ θ                        
Here, and γ 0 comes from the development of the initial stresses. The strain expression allows writing ( 16)

The virtual work of the external loads

Two kinds of load can be considered for these structures: dead loads (loads which always go in the same direction) and follower loads due to the inflation pressure (which inducts local forces acting normally to the surface, and following this surface from the reference configuration to the current configuration). The inflation pressure is considered here as an external force even if it is acting inside the structure.

The virtual work of the external loads (except the pressure) is given by: (

Where T = Π.N and p x , p y , µ are respectively the in-plane components of the dead load and momentum per unit length; X(.), Y(.), Γ(.) are the components of the resultant force and the resultant momentum, which will lead to the boundary conditions.

Concerning the external loads due to the pressure, it is supposed that the beam is subjected to a uniform pressure p on the entire internal surface. The surface is broken down into two different parts: the lateral walls and the surfaces at the end sections. In the following, the force due to the pressure acting on the end section in the initial configuration will be noted: P = pS 0 . Under the assumption that reference volume is a cylinder, the virtual work of the surface density p n is calculated.

(18) = P S * pre dS . W n * p U ∫ ) ( ) ( ) ( V ) ( Y ) ( U ) ( X ) 0 ( ) 0 ( ) 0 ( V ) 0 ( Y ) 0 ( U ) 0 ( X dX ) V p U p ( dS ) ) G ( )( t , P ( d ) ) G ( ).( t , P ( dS ) P ( ). t , P ( d ) P ( ). t , P ( W 0 * 0 0 * 0 0 * 0 * * * 0 * * Y * X S 0 0 0 0 0 0 0 0 0 0 0 0 * ext 0 0 0 0 0 θ Γ θ Γ µθ Ω Ω Ω Ω Ω + + + + + + + + = ∧ + + ∧ + = + = ∂ GP U T GP U W W f U T U f * * * * 0 * * 0 ∫ ∫ ∫ ∫ ∫ 2 XX Y Y γ β α Σ + + = and 2 0 0 0 0 XX Y Y γ β α Σ + + = . - = 0 S

Initial configuration

S p corresponds to the lateral external surfaces and the two bases of the tube. After development (see [START_REF] Van | Bending and buckling of inflatable beams: some new theoretical results[END_REF]), the work can be written as follows:

(19)

Equations of motion and boundary conditions:

The use of the virtual work principle written in lagrangian variables then leads to the following equations where the dynamic terms appear:

(20)

After part integration and considering the terms U*, V*, θ*, we can obtain the equations of motion:

(21) (22) θ θ θ ρ = ∂ ∂ + - ∂ ∂ - ∂ ∂ + ∂ ∂ + - ∂ ∂ + ∂ ∂ + + ∂ ∂ ∂ ∂ - ∂ ∂ + ∂ ∂ - ∂ ∂ )) X U 1 ( sin X V (cos P X K X EK 2 1 S I N T ) sin X V cos ) X U 1 (( sin ) X V M ( X cos )) X U 1 ( M ( X t I X , 0 2 0 2 2 0                   Y 2 2 0 0 p X cos P ) cos T ( X ) X sin M ( X ) X V N ( X t V S = ∂ ∂ + ∂ ∂ - ∂ ∂ ∂ ∂ - ∂ ∂ ∂ ∂ - ∂ ∂ θ θ θ θ θ ρ X 2 2 0 0 p X sin P ) sin T ( X ) X cos M ( X )) X U 1 ( N ( X t U S = ∂ ∂ - ∂ ∂ + ∂ ∂ ∂ ∂ - ∂ ∂ + ∂ ∂ - ∂ ∂ θ θ θ θ θ ρ * int * pre * ext * acc W W W W + + = 0 0 0 * * * 0 * * * pre sin V cos U P dX )) X U 1 ( sin X V (cos X cos V X sin θ θ θ + + ∂ ∂ + - ∂ ∂ + ∂ ∂ - ∂ ∂ = ∫           ) 0 ( X ) 0 ( cos P ) 0 ( sin ) 0 ( T ) 0 ( X ) 0 ( cos ) 0 ( M )) 0 ( X U 1 )( 0 ( N - = - - ∂ ∂ + ∂ ∂ + θ θ θ θ ) ( X ) ( cos P ) ( sin ) ( T ) ( X ) ( cos ) ( M )) ( X U 1 )( ( N 0 0 0 0 0 0 0 0 0 = - + ∂ ∂ + ∂ ∂ + θ θ θ θ θ ) 0 ( Y ) 0 ( sin P ) 0 ( cos ) 0 ( T ) 0 ( X ) 0 ( sin ) 0 ( M ) 0 ( X V ) 0 ( N - = - + ∂ ∂ + ∂ ∂ θ θ θ ) ( Y ) ( sin P ) ( cos ) ( T ) ( X ) ( sin ) ( M ) ( X V ) ( N 0 0 0 0 0 0 0 0 0 = - + ∂ ∂ + ∂ ∂ θ θ θ θ ) 0 ( ) 0 ( X K ) 0 ( X EK 2 1 S I ) 0 ( N ) 0 ( X V ) 0 ( sin ) 0 ( M ) 0 ( cos )) 0 ( X U 1 )( 0 ( M 0 2 0 0 Γ θ γ θ θ θ - = ∂ ∂ + ∂ ∂ + + ∂ ∂ + ∂ ∂ +                         ( ) ) ( X K ) ( X EK 2 1 S I ) ( N ) ( X V ) ( sin ) ( M ) ( cos )) ( X U 1 )( ( M 0 0 0 0 2 0 0 0 0 0 0 0 0 0 Γ θ γ θ θ θ = ∂ ∂ + ∂ ∂ + + ∂ ∂ + ∂ ∂ +
In order to calculate analytical solutions, the same assumptions than those made in [START_REF] Van | Bending and buckling of inflatable beams: some new theoretical results[END_REF] are made. After linearization, the reference and spatial configurations can be considered as the same approximately. The computational viewpoint remains Lagrangian; however, use will be made of small characters (i.e. Eulerian notations) in order to designate the variable. The constitutive laws become: (25) here, k is the correction shear coefficient. For circular thin tubes, k =0.5.

The equations of motion can therefore be deduced:

(26) (27) (28)
The boundary conditions are written for the generalized stresses:

(29)

(30) (31) (32) (33)
(34)

CONSTRUCTION OF THE ELEMENTARY DYNAMIC STIFFNESS MATRIX 3.1 Solutions for the differential equation of motion

In the following, we will consider an inflated tube submitted to bending local external loads and torques at its two ends, and to its internal pressure. Hence, we consider that p x , p y , µ are nil. Eqns (25), ( 29) and (30) give:

(35) P ) ( N ) 0 ( N 0 0 0 = = P ) ( N ) 0 ( N 0 0 0 = = ) 0 ( ) 0 ( x ) S I ) 0 ( N EI ( 0 0 0 0 Γ θ - = ∂ ∂ + ) ( Y ) ( ) kGS P ( ) ( x V ) kGS ) ( N ( 0 0 0 0 0 0 0 = + - ∂ ∂ + θ ) 0 ( Y ) 0 ( ) kGS P ( ) 0 ( x V ) kGS ) 0 ( N ( 0 0 0 - = + - ∂ ∂ + θ ) ( X P ) ( N 0 0 0 = - ) 0 ( X P ) 0 ( N 0 - = - µ θ θ θ ρ = - ∂ ∂ + - ∂ ∂ + - ∂ ∂ ) x V )( kGS P ( x ) S I N EI ( t I 0 2 2 0 0 0 0 2 2 0 0 y 0 2 2 0 0 2 2 0 0 p x ) kGS P ( x V ) kGS N ( t V S = ∂ ∂ + + ∂ ∂ + - ∂ ∂ θ ρ x x , 0 2 2 0 0 p N t U S = - ∂ ∂ ρ 0 N N = ) x V ( kGS T 0 θ - ∂ ∂ = x EI M 0 ∂ ∂ = θ Figure 3.

nodal forces and torques

By neglecting the effect of rotation inertia, the equation of motion for free vibrations can be obtained:

(36) (37)
which allows to write, after some manipulations:

(38)

One assumes that the solution takes the shape:

(39)

The differential equation can be written as follows:

(40)

and leads to the usual result:

(41)

= ⋅ + ⋅ ⋅ + - ⋅ ⋅ + ρ ρ ) t ( g ) x ( v ) t , x ( V ⋅ = 0 x V t x V kGS P S t V S PI EI S 4 4 2 2 4 0 0 0 2 2 0 0 0 0 0 = ∂ ∂ + ∂ ∂ ∂ ⋅ + - ∂ ∂ ⋅ + ρ ρ 0 ) x V )( kGS P ( x ) S I P EI ( 0 2 2 0 0 0 = - ∂ ∂ + - ∂ ∂ + - θ θ 0 x x V ) kGS P ( t V S 2 2 0 2 2 0 0 = ∂ ∂ - ∂ ∂ + - ∂ ∂ θ ρ       F 1 F 2 Γ 1 Γ 2 v(0) θ(0 ( ) 0 v ( ) 0 A B x y θ
Harmonic variation of g with circular frequency ω is assumed then g(t) = cos(ωt -ϕ); then:

(42)

Let us consider a solution of eqn (42) in the exponential form:

(43)

Introducing solution eqn (43) into eqn (42) yields to the characteristic equation:

(44)

The roots can now be obtained:

+ + + + + = ρ ω ρ ω ρ ω λ 1 1 i s λ ± = ; 2 2 s λ ± =

Expression of displacements in the matrix form

The expression of the transversal displacement and the cross-section rotation can be written as:

= {v 1 , θ 1 , v 2 , θ 2 } [ ] - = ) ( ch ) ( sh ) cos( ) sin( ( sh ) ) ( ch ) sin( ) cos( 0 0 0 1 0 1 ) ( 0 2 2 2 1 0 2 2 2 1 0 1 1 2 2 0 1 1 2 2 0 2 0 2 0 1 0 1 2 2 1 1 2 2 λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ ω W                   { } [ ] ) t cos( D C B A ) ( ϕ ω ω - = W U               ) t cos( D ) ( ch C ) ( sh B ) cos( A ) sin( 0 2 0 2 0 1 0 1 2 ϕ ω λ Θ λ Θ λ Λ λ Λ θ - ⋅ ⋅ + ⋅ + ⋅ - ⋅ =     ) t cos( ) ( Dsh ) ( Cch ) sin( B ) cos( A V 0 2 0 2 0 1 0 1 2 ϕ ω λ λ λ λ - ⋅ + + + =     ) t cos( ) D B ( 1 ϕ ω Θ Λ θ - ⋅ ⋅ + ⋅ - = ) t cos( ) C A ( V 1 ϕ ω - ⋅ + = 2 2 1 λ λ Λ - = ⋅ ; 2 1 2 λ λ Θ = ⋅ ) t cos( D ) x ( ch C ) x ( sh B ) x cos( A ) x sin( 2 2 1 1 ϕ ω λ Θ λ Θ λ Λ λ Λ θ - ⋅ ⋅ + ⋅ + ⋅ - ⋅ =     [ ] [ ] [ ] [ ] ) t cos( ) x ( Dsh ) x ( Cch ) x sin( B ) x cos( A ) t , x ( V 2 2 1 1 ϕ ω λ λ λ λ - • + + + = ) ( ) x ( Dch ) x ( Csh ) x cos( B ) x sin( A ) x ( ~) x ( Dch ) x ( Csh ) x cos( B ) x sin( A ) x ( ~) x ( Dch ) x ( Csh ) x cos( B ) x sin( A ) x ( ~) t cos( ) x ( kGS P S ) x ( kGS P S PI EI ) x ( ) t , x ( 2 2 2 2 2 1 1 1 1 3 2 3 2 2 3 2 1 3 1 1 3 1 2 2 2 2 2 1 1 1 1 1 3 0 0 0 2 0 0 0 0 1 ω λ λ λ λ λ λ λ λ θ λ λ λ λ λ λ λ λ θ λ λ λ λ λ λ λ λ θ ϕ ω θ ρ θ θ θ - • + + + - = + + - = + + + - = - • + - + + + =                

Expression of the external loads in the matrix form

In the same way, we will obtain the expression of the external loads at the ends of the beam by using the boundary conditions (31-34):

-For the transversal forces F 1 and F 2 :

(59)

-For the torques Γ 1 and Γ 2:

(60)

Substituting into ∂V ∂x , ∂θ ∂x , θ and simplifying gives:

(61) With (62) (63)
The values of the shear loads and bending momenta at the ends of beam are finally: At the left-end x = 0:

(64) (65)
At the right-end x = ᐉ 0 :

(66)

The external loads are also written in matrix form:

(68) { } [ ] ) t cos( D C B A ) ( ϕ ω ω - = Z F               ) t cos( D ) ( sh C ) ( ch B ) sin( A ) cos( ) EI ( 0 2 2 1 0 2 2 1 0 1 2 2 0 1 2 2 p 2 ϕ ω λ λ λ λ λ λ λ λ Γ - ⋅ ⋅ + ⋅ + ⋅ - ⋅ - =     [ ] ) t cos( D ) ( ch ) ( C ) ( sh ) ( B ) cos( A ) sin( ) ( F 0 2 1 0 2 1 0 1 2 0 1 2 2 ϕ ω λ λ λ λ λ λ λ λ ϑ - ⋅ ⋅ - + ⋅ - + ⋅ + ⋅ - = ) t cos( ) C A ( ) EI ( 2 1 2 2 p 1 ϕ ω λ λ Γ - ⋅ ⋅ + ⋅ - = - ) t cos( ) D B ( F 1 2 1 ϕ ω λ λ ϑ - ⋅ ⋅ - ⋅ = - 0 0 0 P S PI EI ) EI ( + = and
) kGS P ( ) kGS ( 0 p + = p 0 0 2 ) EI ( S ρ ω ϑ = 1 0 0 0 ) 0 ( x ) S I P EI ( Γ θ - = ∂ ∂ + 2 0 0 0 0 ) ( x ) S I P EI ( Γ θ = ∂ ∂ + 1 0 F ) 0 ( ) 0 ( x V ) kGS P ( - = - ∂ ∂ + θ 2 0 0 0 F ) ( ) ( x V ) kGS P ( = - ∂ ∂ + θ             [ ] ) t cos( D ) x ( ch ) ( C ) x ( sh ) ( B ) x cos( A ) x sin( ) ( ) x V ( ) kGS ( 2 1 2 1 1 2 1 2 p ϕ ω λ λ λ λ λ λ λ λ ϑ θ - ⋅ ⋅ - + ⋅ - + ⋅ + ⋅ - = - ∂ ∂ ) t cos( D ) x ( sh C ) x ( ch B ) x sin( A ) x cos( ) EI ( x ) EI ( 2 2 1 2 2 1 1 2 2 1 2 2 p p ϕ ω λ λ λ λ λ λ λ λ θ - ⋅ ⋅ + ⋅ + ⋅ - ⋅ - = ∂ ∂    
where (69)

and with

{F} T = {F 1 , Γ 1 , F 2 , Γ 2 }

The dynamic stiffness matrix

The dynamic stiffness matrix [D(ω)] is given by:

(70) and [D(ω)] = [Z(ω)][W(ω)] -1 (71)
The dynamic stiffness matrix is symmetrical, and its components are:

( ω D             ( ) [ ] = 2 2 1 1 2 2 1 1 V V w F F θ θ Γ Γ D                             [ ] - - - - - - - = ) ( sh ) EI ( ( ch ) EI ( ) ) sin( ) EI ( ) cos( ) EI ( ) ( ch ) ( sh ) cos( ) sin( 0 ) EI ( 0 ) EI ( 0 0 ) ( 0 2 2 1 p 0 2 2 1 p 0 1 2 2 p 0 1 2 2 p 0 2 1 0 2 1 0 1 2 0 1 2 2 1 p 2 2 p 1 2 λ λ λ λ λ λ λ λ λ ϑλ λ ϑλ λ ϑλ λ ϑλ λ λ ϑλ ϑλ ω Z               [ ] ) ( sh ) sin( ) ( d D D 0 2 3 1 0 1 3 2 2 2 2 1 2 1 31 13 λ λ λ λ λ λ λ ϑλ + + - = = [ ] ) ( ch ) cos( ) ( d D D 0 2 0 1 2 2 2 1 2 2 2 1 41 14 λ λ λ λ λ ϑλ - + - = = [ ] ) ( sh ) cos( ) ( ch ) sin( ) ( d ) EI ( D 0 2 0 1 3 2 0 2 0 1 3 1 2 2 2 1 2 1 p 22 λ λ λ λ λ λ λ λ λ λ - + = [ ] ) ( ch ) cos( ) ( d D D 0 2 0 1 2 2 2 1 2 2 2 1 32 23 λ λ λ λ λ ϑλ - + = = ) ( ch ) sin( ) ( sh ) cos( ) ( d D 0 2 0 1 3 2 0 2 0 1 3 1 2 2 2 1 2 1 11 λ λ λ λ λ λ λ λ λ ϑλ + + = ) 1 ) ( ch ) )(cos( ( ) ( sh ) sin( d D D 0 2 0 1 3 2 1 2 3 1 0 2 0 1 4 2 4 1 2 1 21 12 - - + + = = λ λ λ λ λ λ λ λ λ λ λ ϑλ            
where

The inflation pressure appears in the equations and each term of the matrix D is a transcendental function of ω.

Simple changes of coordinates give the matrixes needed in order to compute the bending natural frequencies and mode shapes of inflatable tubes. This leads to the non-linear eigenvalue problem:

(73) Natural frequencies can be calculated by the use of the well-known algorithm defined by Wittrick and Williams [START_REF] Wittrick | A general algorithm for computing natural frequencies of elastic structures[END_REF], which enables the calculation of the number of natural frequencies of the structure that exist below an arbitrary chosen trial frequency. It allows us to calculate the natural frequencies with a chosen accuracy. Numerous papers deal with the determination of natural frequencies and natural modes of beams with the CEM with accuracy ( [START_REF] Banerjee | Exact dynamic stiffness matrix for composite Timoshenko beams with applications[END_REF], [START_REF] Banerjee | Development of dynamic stiffness matrix for free vibration analysis of spinning beams[END_REF]). In this paper, we have preferred to focus our attention on two points: the first is the coherence of the results with the conventional finite element method; the second point is concerned with the influence of the inflation pressure in the results which is of great importance for inflatable structures.

F.E. STIFFNESS MATRIX AND MASS MATRIX

The aim of this section is to deduce FEM stiffness and mass matrixes dedicated to inflatable tubes from the dynamic stiffness matrix.

From CEM to FEM:

Free vibrations may be calculated by the use of the dynamic stiffness matrix, which leads to the non-linear eigenvalue problem:

(74) This problem may be rewritten as follows:

(75)

The mass and stiffness matrixes are obtained with the help of Leung's theory [START_REF] Richards | An accurate method in structural vibration analysis[END_REF]:

(76) then (77) ( ) ( ) ( ) ω ω ω ω M D K 2 + = ( ) ( ) 2 ω ω ω ∂ ∂ - = D M ( ) ( ) ( ) 0 = U M - K 2 ω ω ω ( ) 0 U D = ω ( ) F U D = ω ) 1 ) ( ch ) (cos( 2 ) ( sh ) sin( ) ( d 0 2 0 1 3 2 3 1 0 2 0 1 6 2 6 1 - - - = λ λ λ λ λ λ λ λ [ ] ) ( sh ) sin( ) ( d ) EI ( D D 0 2 3 2 0 1 3 1 2 2 2 1 2 1 p 42 24 λ λ λ λ λ λ λ λ - + - = = [ ] ) ( ch ) sin( ) ( sh ) cos( ) ( d D 0 2 0 1 3 2 0 2 0 1 3 1 2 2 2 1 2 1 33 λ λ λ λ λ λ λ λ λ ϑλ + + = ( ) [ ] ) 1 ) ( ch ) )(cos( ( ) ( sh ) sin( d D D 0 2 0 1 3 2 1 2 3 1 0 2 0 1 4 2 4 1 2 1 43 34 - - + + - = = λ λ λ λ λ λ λ λ λ λ λ ϑλ [ ] ) ( sh ) cos( ) ( ch ) sin( ) ( d ) EI ( D 0 2 0 1 3 2 0 2 0 1 3 1 2 2 2 1 2 1 p 44 λ λ λ λ λ λ λ λ λ λ - + =
Free vibrations may also be studied by the use of the finite element method. The displacement field is discretized:

(78)
where N is the shape function matrix, and V is the vector of generalized co-ordinates. The FEM leads to the linear eigenvalues problem:

(79)

The FEM stiffness and mass matrixes will then be calculated with:

(80) (81) The pressure P appears explicitly in the stiffness matrix. If the pressure tends to 0, one can recognize the Timoshenko beam stiffness matrix. This stiffness matrix can be deduced from the results derived by Wielgosz [START_REF] Wielgosz | An inflatable fabric beam finite element[END_REF], by replacing E*I* by EI o +PI o /S o and pbh by P+kGS 0 .

The stiffness and mass matrix obtained with the dynamic stiffness matrix

The limit of M(ω) with ω→0, gives the mass matrix M:

(84)

φ φ φ φ K               ( ) ( ) 2 0 0 lim lim ω ω ω ω ω ∂ ∂ - = = → → D M M ( ) ω ω D K 0 lim → = ( ) 0 = U M - K 2 ω [ ]{ } ) t ( ) x ( ) t , x ( v V N = [ ] + + + + - + + + + - + + + + + + - + + + + + + + =
φ φ φ φ φ φ φ φ φ φ φ φ φ φ φ φ φ φ ρ M                                                    
The inflation pressure also appears in the mass matrix. If P tends to 0, the mass matrix also becomes the Timoshenko beam mass matrix.

The stiffness matrix by the left quotient operation:

The foundations of this algebraic operation can be found in [START_REF] Nayroles | Opérations algébriques en mécanique des structures[END_REF]. The diagram given figure 4 can draw the left quotient operation. The continuous element is described by a "mechanical element": the two vector spaces V (displacement space) and Φ (loads space) are dual for a bilinear form <.;.> 1 , which represents the work of a load {ϕ} ∈ Φ in a displacement {v} ∈ V. The load law L (application from V into Φ) summarizes the equilibrium equations and the boundary conditions. This has already been written for an inflatable panel [START_REF] Wielgosz | Equilibrium finite elements for inflatable beams[END_REF] element submitted to nodal concentrated loads. The vector spaces are:

(85)

The static load law and the load vector are:

(86) with {ϕ} T = {F 1 , Γ 1 , F 2 , Γ 2 }
The principle of the left quotient operation is to describe the equilibrium element by another "mechanical element". Φ 1 is the discreet subspace of the loads space vector and contains the nodal loads (here Φ 1 = R 4 ). If j is the canonical injection from Φ 1 into Φ, then V/Φ 1 0 is the quotient space (i.e. the subspace of the nodal displacements, V/Φ 1 0 = R 4 ). {ϕ ~} ∈ Φ 1 represents the nodal loads (F 1 , Γ 1 , F 2 , Γ 2 ) and {v ~} ∈ V/Φ 1 0 represents the nodal displacements (v 1 , θ 1 , v 2 , θ 2 ). These two vector spaces are dual for the bilinear form <.;.> 2 , work of the nodal loads in the nodal displacements:

(87)

The principle of the left quotient operation is to solve the following problem: if the nodal loads {ϕ ~} ∈ Φ 1 are given, find the nodal displacements {v ~} ∈ V/Φ 1 0 . The displacement {v ~} ∈ V is obtained by integration of the equilibrium equations. The nodal displacements are therefore given by {v ~} = j t ({v}). The result gives the compliance of the element {v ~} = (j t o L -1 o j)({ϕ ~}). If the nodal displacements {v ~} ∈ V/Φ 1 0 are given, find the nodal loads {ϕ ~} ∈ Φ 1 . The result now gives the stiffness matrix of the element: {ϕ ~} = (j t o L -1 o j) -1 ({ϕ ~})

We have to inverse the compliance given by the first problem, but it always leads to a singular system (because

{ }{ } { }{ } { } { } { } { }) v ( j v ) ( j , v , v ) , / V ( ) , v ( 1 t 1 2 1 0 1 - ∈ ∀ > =< > < → ∈ ϕ ϕ Φ Φ ϕ ( ) ( ) , ) 0 ( ) 0 ( dx dv PS kGS , dx d dx v d , dx dv ) kGS P ( dx d EI ) , v ( L 0 0 2 2 0 2 2 0 - + - - - + + = θ θ θ θ θ ( ) ( ) ( ) - + - ) ( dx d EI , ) ( ) ( dx dv pS kGS ), 0 ( dx d EI 0 0 0 0 0 0 0 θ θ θ                         [ ] ( ) R ; , 0 C V 0 4 = [ ] ( ) 4 0 R R ; , 0 C × = Φ V V/ j <.,.> 1 Φ T ↓ ↑ j Φ 1 o <.,.> 2 Φ 1
of the rigid body displacements) and the building of the stiffness matrix needs other equations, which come from the equilibrium equations. Following the process detailed in [START_REF] Wielgosz | Equilibrium finite elements for inflatable beams[END_REF], one obtains:

(88)

Concerning the inflatable tubes, the load law is obtained with the static equilibrium equations obtained by elimination of the dynamic terms in eqns (36-37-59-60):

(89) Hence, an obvious analogy allows us to obtain the exact stiffness matrix of the inflatable tube. Only the bending rigidity changes (just replace EI 0 by

). This proves that the stiffness matrix obtained by calculating the limit of [D(ω)] is the exact stiffness matrix.

The FEM stiffness and mass matrix

The left quotient operation doesn't allow the calculation of the mass matrix. Once the coherence between FEM and CEM is verified for the stiffness matrix, the usual method is used to propose shape functions and calculate the stiffness and mass matrix by the use of the virtual work principle. After linearization, by neglecting the distributed loads from eqn [START_REF] Wielgosz | Equilibrium finite elements for inflatable beams[END_REF], the following expression for the flexion of an inflatable beam in dynamics is obtained:

(90)
The use of the constitutive equations allows us to write the virtual work principle in terms of displacements, and after some manipulation, the equation may be written under a simple and practical form:

(91) + - - ∂ ∂ + = + + + ∫ ) dx dV )( x V )( kGS P ( ) ( ) ( V F ) 0 ( ) 0 ( V F * 0 * 0 0 * 2 0 * 2 * 1 * 1 0 θ θ θ Γ θ Γ dx t I V t V S dx dx d x S PI EI 0 0 * 2 2 0 0 * 2 2 0 0 * 0 0 0 ∫ ∂ ∂ + ∂ ∂ + ∂ ∂ + θ θ ρ ρ θ θ                   = + - ∂ ∂ + ∂ ∂ - + + + + ∫ 0 0 0 0 * * 0 * 0 * 2 0 * 2 * 1 * 1 V P dx ) x V ( x V P ) ( ) ( V F ) 0 ( ) 0 ( V F θ θ θ θ θ Γ θ Γ dx dx d x S I P M T dx dV ) T x V P ( dx t I t V V S * 0 0 * 0 * 0 * 2 2 0 0 2 2 * 0 0 0 θ θ θ θ θ ρ ρ ∂ ∂ + + - + + ∂ ∂ + ∂ ∂ + ∂ ∂ ∫ ∫                       0 0 0 S PI EI + ( ) , ) 0 ( ) 0 ( dx dv pS kGS , dx d dx v d , dx dv ) kGS pS ( dx d S PI EI ) , v ( L 0 0 2 2 0 0 2 2 0 0 0 - + - - - + + + = θ θ θ θ θ ( ) + - + + - ) ( dx d S PI EI , ) ( ) ( dx dv pS kGS ), 0 ( dx d S PI EI 0 0 0 0 0 0 0 0 0 0 0 θ θ θ                                           [ ] + - - - + - + = ) 4 ( sym 6 12 ) 2 ( 6 ) 4 ( 6 12 6 12 ) 1 ( ) EI ( p 2 0 0 p 2 0 0 p 2 0 0 0 p 3 0 0 φ φ φ φ K with 2 0 0 0 P ) kGS P ( ) EI ( 12 + = φ              
The virtual work principle for static analysis can then be written, which is used here to obtain the stiffness matrix:

(92)

The equations of motion are easily obtained by elimination of the dynamic terms in eqns (36-37). For an Euler Bernoulli beam, the displacement is written as follows:

(93)

where the shape functions are the usual cubic expansions:

(94)

The relation is valuable for an Euler Bernoulli beam. It allows us to obtain the shape functions for the section's rotation. For inflatable beams, which are Timoshenko's beams, the cross-section does not remain orthogonal to the neutral fibre. Hence dv/dx ≠ θ. The preceding functions interpolate displacements and the derivative of displacements, which is now written as:

(95) where the shape functions

N v ′ 1 (x) and N v ′ 2 (x) are: (96)
In this case, the following boundary conditions are, of course, checked: The boundary conditions give the linear system:

(97) ' 2 0 2 0 ' 1 1 v ) ( dx dV v ) ( V v ) 0 ( dx dV v ) 0 ( V = = = = - = 0 0 2 v x 1 x ) x ( N ' 1 - - = 0 0 2 v x 1 x ) x ( N ' 2             [ ] [ ] {} V N v ) x ( v v v v ) x ( N ) x ( N ) x ( N ) x ( N ) x ( V ' 2 2 ' 1 1 v v v v ' 2 2 ' 1 1 = =               + - = 0 2 0 1 v x 2 1 x 1 ) x ( N - = 0 2 0 2 v x 2 3 x ) x ( N - = 0 0 2 1 x 1 x ) x ( N θ - - = 0 0 2 2 x 1 x ) x ( N θ                                     [ ] [ ]{ } V N = = ) t ( ) t ( v ) t ( ) t ( v ) x ( N ) x ( N ) x ( N ) x ( N ) x ( V 2 2 1 1 v v 2 2 1 1 θ θ θ θ               dx dx d dx d S PI EI ) dx dV )( dx dV )( kGS P ( ) ( ) ( V F ) 0 ( ) 0 ( V F * 0 0 0 * 0 * 0 0 * 2 0 * 2 θ Γ θ Γ + + - - + = + + + ∫       By changing the unknowns {V ~} = {v 1 , v ′ 1 , v 2 , v ′ 2 } to the unknowns {V} = {v 1 , θ 1 , v 2 , θ 2 },
(100)

The following relation is then easily obtained:

(102)

And the shape functions can be calculated by:

(103)

The shape functions depends on the inflation pressure, the Young modulus and the shear modulus:

) The virtual work of the internal loads can then be written as:

P 2 0 0 P 0 2 P 2 0 0 P 0 0 1 P 3 0 2 0 P 0 2 2 v P 3 0 2 0 2 0 P 0 2 0 1 v 1 x 2 x x 2 1 ) x ( N 1 2 x 2 x x 2 1 ) x ( N 1 x 3 x 2 x ) x ( N 1 x x 2 x ) x ( N φ φ φ φ φ φ φ φ θ θ + + - = + - - - = + - - - = + - - - - = [ ][ ]{ } V A N = ) x ( V [ ]{ } V A = + - - - + - + =
v v 2 1 2 0 1 0 0 2 2 1 0 0 0 1 1 1 v v v v θ θ φ φ φ φ φ φ φ φ φ                                           + - - - + - + =
P P ' 2 ' 1 v v 2 1 2 2 2 1 1 1 v v θ θ φ φ φ φ φ φ φ φ φ                                 [ ] = ' 2 2 ' 1 1 42 2 1 v v v v M θ θ                    
(106) And by use of the shape functions:

(107)

The stiffness matrix of the inflatable tube is then calculated, which is nothing else than the stiffness matrix calculated from the dynamic stiffness matrix:

The use of the shape functions now allows the calculation of the mass matrix:

(108) (109)
The term gives this mass matrix. The calculation of the different terms of the matrix is simplified by the use of formal calculus software. Finally, the mass matrix and the stiffness matrix obtained here are exactly the same as those given by eqns (82-84), which shows the coherence of these developments.

SOLUTIONS FOR A SIMPLY SUPPORTED INFLATABLE BEAM

The theory is now applied to a simply supported inflatable beam. The boundary conditions are:

(110) 0 ) t , 0 ( V = 0 ) t , ( V 0 = 0 ) t , 0 ( x = ∂ ∂θ 0 ) t , ( x 0 = ∂ ∂θ [ ] [ ]dx ) x ( N ) x ( N S 0 0 T 0 0 ∫ ρ { } [ ] [ ] [ ] [ ] = ∂ ∂ ∫ ∫ 2 2 0 T T T * 0 0 0 2 2 * dt V d A dx ) x ( N ) x ( N A V S 0 0 S d x t V V 0 0 ρ ρ             [ ][ ]{ } ) t ( V A ) x ( N ) t , x ( V = [ ] + - - - + - + + = ) 4 ( sym
φ φ φ φ K               { }[ ] { } V dx dx dN dx dN N dx dN N dx dN ) 12 ( A V S PI EI W T v 0 T v P 2 0 * 0 0 0 * 0 + - - + = ∫ θ θ θ θ φ δ                               dx dx d dx d ) dx dV )( dx
- - + = + + - - + ∫ ∫             + + + - + - = - 2 2 1 1 p p 0 p p p p 0 p p v v 1 2 1 1 2 1 dx dV θ θ φ φ φ φ φ φ φ φ θ                    
The natural frequencies of vibration of an inflatable simply supported beam are solutions of:

One gets:

Hence, with d ≠ 0, one can obtain:

(114)

We find:

The natural frequency can finally be calculated by the use of eqn (48):

(116)

The inflation pressure appears explicitly in the formulation of the frequency, which increases with the pressure.

+ + + = π ρ ρ π 0 1 nπ λ ± = or 0 2 = λ 0 ) sin( 0 1 = λ or 0 ) ( sh 0 2 = λ 0 ) 1 ) ( ch ) (cos( 2 ) ( sh ) sin( ) ( ) ( sh ) sin( 0 2 0 1 3 2 3 1 0 2 0 1 6 2 6 1 0 2 0 1 = - - - λ λ λ λ λ λ λ λ λ λ     0 ) ( sh ) sin( ( ) ( d EI ) ( sh ) cos( ) ( ch ) sin( ) ( d EI 2 0 2 3 2 0 1 3 1 2 2 2 1 2 1 0 2 0 2 0 1 3 2 0 2 0 1 3 1 2 2 2 1 2 1 0 = - + - - - + λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ                     0 D D D D 44 42 24 22
= their influence on the solution. Increasing the pressure inducts an increase of the natural frequencies, but these growths are quite small and depend on the natural frequency. Multiplying by three the pressure leads in this case to an increase of the frequencies by 2.7%, 6% and 9% respectively for n=1, 2, 3 in the most varying cases. In fact, the first frequency is mainly dependent on the Young modulus. For a given Young's modulus, the shear coefficient G has a great influence on the interval between the natural frequencies.

COMPARISONS WITH THE FINITE ELEMENT METHOD RESULTS

The first three natural frequencies are calculated with the two methods. The CEM values are simply obtained with eqn (116). The three first frequencies are calculated for the simply supported inflatable tube, which is used for the experiments (see section 7). Geometry and material are such that: ᐉ 0 = 1,858m, R 0 = 0,0831m, Ee = 179000 Pa.m, Ge = 20000 Pa.m, ρ 0 e = 0,3759kg.m -2 . Figure 6. Evolution of frequency according to the number of element reference Figure [START_REF] Thomas | Deflections of highly inflated fabric tubes[END_REF] shows that the results of FEM converge with the results of CEM with the increase of the number of elements. If the CEM gives the exact frequency with a chosen accuracy, it is not the case of the FEM, However, the results of FEM are always greater than the ones of the CEM, as expected.

COMPARISONS WITH THE EXPERIMENTAL RESULTS

Experiments have been carried out to measure the natural frequencies of an simply inflatable supported beam. Sewing two layers of fabric manufactures the beam used for these experiments. This fabric is a woven network, for which the fibbers are initially orthogonal and imply an orthotropic behaviour. Because we focus the studies on beam models, it will be considered here as an hyperelastic isotropic Saint-Venant Kirchhoff material. A thin PVC membrane incorporated inside the beam ensures the air tightness. The vibration of the beam comes from an electromagnet. The induction coil has been separated form the electromagnet. It is very light and just sticked on the beam. In order to compare experiments and theoretical results, the Young modulus has to be identified. This is done with the variation of the initial geometry, which clearly shows the importance of the pres-tress state. It is easy to calculate the circumferential and longitudinal stresses in the membrane due to the internal pressure of the beam in the initial configuration:

(117) e corresponds to the thickness of the fabric. Since we consider an elastic isotropic material, the length and the radius may be modelled with:

(118) Here, ᐉ φ and R φ are the length and the radius for a very low internal pressure, and can be obtained with the curves Fig. 8. b with the linear interpolations of the two curves length versus pressure and radius versus pressure. The Young modulus is also determined with these linear interpolations. The average value of the product of the Young modulus with the thickness e is estimated to Ee=178973 Pam. The Young modulus is the main influent parameter for the first natural frequency. The product Ge is then estimated in order to ensure the interval between the other natural frequencies (here Ge=20000 Pam). The comparison between theoretical and experimental results is shown fig 10 a. The analytical results are close to the experimental results for the first pressure (50 kPa) but the slopes of the theoretical curves are smaller than the experimental ones. This can be corrected by the use of the influence of the pressure on the radius and on the length on the natural frequencies. The initial state is then correctly defined.

CONCLUSIONS

Inflatable structures are increasingly used nowadays. They are often made of beam elements. In this paper, two methods to calculate the eigenvalues and eigenmodes of structures made of inflatable beams have been propounded. The equations of motion have been written in lagrangian variables for a Timoshenko beam in order to take into account the following forces due to the internal pressure and the shear behaviour of this kind of thinwalled fabric structure. We first have dealt with the continuous element method to establish the dynamic stiffness matrix, which allows us to get exact results. The finite element method stiffness and mass matrixes have been derived from the CEM and also directly from the virtual work principle. In each case, the inflation pressure appears in the matrixes. The natural frequencies have been calculated for a simply supported inflatable beam with the dynamic stiffness matrix, and these results have been compared to FEM and to experimental ones. The influence of the material parameters has been pointed out. The first frequency is mainly dependent on the Young modulus and for a given Young's modulus, the shear coefficient has a great influence on the interval between the natural frequencies. The comparison between the experimental results and the theoretical results has shown the importance of the pre-stress effect, which increases the radius and the length of the initial configuration.

Figure 1

 1 Figure 1. reference and current configuration

Figure 2 .

 2 Figure 2. following forces

λ 1

 1 and λ 2 are functions of the circular frequency ω and of the internal pressure.

-

  eqn (50), A, B, C and D are undetermined constants. After reduction, one can obtain: the displacements and rotations at the ends of the beams are: -At the left-end x = 0: At the right-end x = ᐉ 0 : (55) (56)We can now write the displacements in a matrix form:

  By calculating the limit of [D(ω)] when ω→0, one obtains [K]. (82) where Φ P has been introduced for convenience: (83)

Figure 4 .

 4 Figure 4. the left quotient operation

Figure 5 .

 5 Figure 5. natural frequencies of a simply supported inflatable tube: influence of the internal pressure

Figure 7 .

 7 Figure 7. Simply supported beam experiment

Figure 9 .

 9 Figure 9. pre-stress effect

  Fig 10. b presents the comparison between the experimental frequencies and the theoretical frequencies taking into account the variation of the length and radius versus pressure, and the results show a correct accuracy between theory and experiment. a) without pre-stress effect b) with pre-stress effect Figure10. Comparison between theoretical and experimental results