A new fluid-structure coupling model: Case of a pipe under internal pressure

Abstract

This work presents the eigenfrequencies of a study on the fluid-structure interaction for a pipe with a circular section, under pressure constraints in laminar, incompressible and irrotational flow, followed by a comparison between stiffness coupling and bar coupling with different geometric ratio Radius/Length; Radius/Thickness, as a frequency geometric identity of the pipe in the industry. The realization of this study is to have better coupling, where we use the relations of mechanical behaviour of the solid on displacement-stress and for the fluid, we use the Navier-Stokes equations in cylindrical coordinates under the speed-stress form transformed into displacement by the Galerkin-Temp theory. Either the force of the fluid will be distributed on the internal surface of the cylinder by the coupling (interface) between fluid and solid. To obtain the results of the eigenfrequencies we used the hierarchical finite element method (HFEM) for both mechanical and fluid entities. For the natural frequencies, we apply the hierarchical finite elements of the Legendre polynomial which is defined in a well-known interval, thanks to this mathematical calculation we arrive at the general equations of motion matrix form of the solid, the interface and the fluid. The calculation was carried out by a MATLAB programme which determines the eigenfrequencies of the evolved system by different geometric and physical parameters mentioned. We validated this work by a table of comparison between the experimental values and the values of the programme.

Introduction

The fluid-structure interaction is the behaviour between the mechanics of fluids and that of the solid by action and reactions to the stresses applied to each other by natural or industrial influences on the flow inside or/and outside of the structure and is the phenomenon in it by all the sides of our industrial life, on the vibration which is propagated continuously and arbitrarily, causes the structural destruction in the movement of the fluid in the pipe.

There are several modes of mechanical contact between two media of different or identical natures. The patterns of this rubbing contact can be considered penetrating or non-penetrating. The particularity of the contact problem requires the knowledge of information on the geometry at the level of the boundaries: distance between two nodes, tangential velocity, fluid pressure and the existence of friction or the nature of the contact under study (Figure 1). ¹

Figure 1.

Representation of normal contact.

We will establish the continuity of the work of Hicham et al., ² by a new coupling method with the search for less frequent possibilities at different values of the ratios between radius / length and radius / thickness of interaction fluid structure by the finite element method hierarchical as the method of numerical resolution of a cylindrical structure with circular section conveys a non-turbulent laminar fluid, incompressible cylindrical coordinated irrotational, which depends on the polynomial of modelling used. Quite simply it is the h version instead of the p version which does not require a lot of storage space for the convergence of less convergence time, and also minimizes the chosen finite element refinement followed by the degree of freedom of the polynomial of legend which is presented by Nguyen et al., ³ The convergence of our work has been validated by the work of Zhang et al., ⁴ for then passing to the study of the method used which gives good compatibility between the experimental values and the hierarchical method, which now also uses the speed polynomial in this domain with different possible boundary conditions (single-single, integrated-integrated, free-free, simple-integrated, integrated-free, integrated-simple) to arrive at the minimum natural frequency of the system. We arrive at the final matrix form of the solid and fluid movement by the temporary Lagrangian approach of the uniform fluid by the principle of virtual works with the Eulerian approach of the solid. We will be interested in a geometrically deformable body which is cylindrical, adapted to the chosen boundary conditions with the physical properties modelled by a stress-displacement relation of the pipe, or the Navier-Stokes equation is modelled in stress-velocity and integrated by Newmark’s method which gives us homogeneity to the equation of fluid motion by the penalty coupling in the form of a bar or the influence of the negligible mass which is linked to the cylinder inside. In all this, we develop a programme under MATLAB, validated by several studied examples, to solve the problem of fluid-pipe interaction with different variables of geomatic parameters concerning natural frequencies.

The majority of researchers who are published on the fluid-structure interaction apply the finite element method like Causin et al., ⁵ to solve the differential equations in simulation with the different approaches of the possible methods, which brings us to plan a lot of time to arrive at the convergence or the resolution of the problem by refinement but on the other hand, one finds the method of the hierarchical finite elements by Houmat, ⁶ which succeeds the method of the finite elements or one adds a polynomial of different degrees of freedom instead of refinement (HFEM) by Ouissi and Houmat, ⁷ it is economical from the point of view of storage (minimization of finite elements) and resolution (minimization of the size of the matrices of stiffness and mass) a simple algorithm, or most works do not practice this method on the domain of fluid-structure interaction.

We will highlight some work that has been carried out in the field of fluid-structure interaction with a quotation of the coupling used in these articles on the vibration side or the eigenfrequencies obtained for cylindrical pipes carrying fluid. Among these works we find in the study of Alshabatat, ⁸ which seeks to optimize the natural frequencies of cylindrical shells with thin walls for materials with variable gradient in the longitudinal direction by the trigonometric law of volume fractions which uses the finite element method to modify the stiffness and mass of the structure by the theory of shear deformation formulated by Hamilton’s principle. This validation by the p-version of the element to have compatible values that take a lot of convergence time to refine 36 × 40 elements. Duan et al., ⁹ use the finite element method to study the modal analysis of tubes at high temperature and high pressure to apply the fluid-structure interaction (coupling) of a way where the force of the fluid equals the force of the solid in the form of an inertia matrix to arrive at the dynamic equation of the system which gives a decrease in the eigenfrequencies if there is greater displacement, also the more the pressure increases it there is an increase in the stresses which also gives a low natural frequency for the tube because the higher the temperature, the lower the natural frequency. This time Mi and Zhou, ¹⁰ applied Galerkin’s method for a cylindrical pipe carrying an internally flowing fluid composed of oil-gas-water at different viscosity values to influence natural frequencies with critical velocities found with the coupling uses added mass types in the vibrational kinetic equation of motions, also Wu et al., ¹¹ studied the underwater vibrational causes of a shaft pipeline submarines where they used ANSYS software for geometry then reported in Fluent for vibration analysis and pressure distribution with subsea pipeline velocity. They applied the Navier-Stokes equation following the standard k−ε turbulent model for the state of the moving fluid in the pipeline in this work which is based on the finite element method FEM, on the other hand, Selmi and Hassis, ¹² studied the natural frequencies of the three modes in the vibration analysis on a pipe in graded material which transports fluid. The material properties are assumed to be graded in the thickness direction after solving the problem and obtained by buckling deformation with various boundary conditions. Whether the influence of the speed and the density of the fluid in free vibration is sought, also we have the work of Hyvarinen et al., ¹³ use the experimental software ANSYS by the finite element method FEM for the solid and BEM for the fluid where they will solve their problem by investigating the failure of the hydraulic hose which is dynamic with the presence of the phenomenon of fluid-structure interaction (FSI). They discovered that fatigue in hydraulic hoses comes from modular vibrations or natural frequencies about the speed and pressure of the hydraulics. The works of Wu et al., ¹⁴ are also based on the k−ε model of the turbulence with the pressure force added through the convergence time to have the exact solution of the equation of motion to be able to find the natural frequencies of the pipe with different modes from 1 to 5, with the comparison of these experimental and numerical results of the coupling proposed by coherent with the calculation method used to have the effect of the pressure on the vibration of the pipeline. We also find Quan et al., ¹⁵ studied the axial vibration responses of the fluid-structure interaction (FSI) for the hydraulic lines of the C919 aircraft where they used the coupling in a model of friction like the contact interface between fluid-structure applying high speeds and high pressures, with of course the presence of the Reynolds number in their calculation. Using MATLAB to find the frequencies of the different characteristic poses in this article is to forget the comparison between numerical and experimental simulation results. Fuster et al., ¹⁶ their work applied the Lagrangian method and the Weymouth and Yue (WY) method Euler-implicit method for conservation of mass and momentum of fluids to improve stability and to couple its applied density formulas which vary with incompressible Navier-stokes flow. On the other hand, we find different couplings in previous works such as Guo et al., ¹⁷ applied the coupling in the form of a stress-strain relation by strain-displacement contact with the mass coefficients of the solid and the fluid according to a dimensional summation in inverse time which is equal to zero, to ensure the stability of the fluid-piping-flexible interaction for the two methods used FSI-FVM, the different effects of the free supports, fixed and elastic. Also, Cui et al., ¹⁸ use a multi-field mathematical model of the fluid-solid-magnetic coupling, based on two bidirectional fluid-solid couplings and a unidirectional magnetic structure caused by the fluid pressure deformation of the weld were obtained by the permutation of the fluid penetrated as a force applied the opening of the internal crack which causes the growth of this crack in the weld of the pipe with a discussion on the field of magnetic in the process of internal crack growth. The works of Liang et al., ¹⁹ investigate natural frequencies of fluid velocity, fluid density with the presence of functional gradient materials (FGMs) the volume variation of the internal material by its power gradient for a pipe composed of metal and ceramics the transport of fluid and the proposed coupling it is coupling by the mass added in the form of density multiplied by the volume. Kesimli et al. ²⁰ studied the linear vibration of a fluid transport pipe for single supports of variable position between two supports clamped at the end, or they used the coupling in the form Lagrangian in kinetic energies for fluid and potential energies for the structure of the fluid-structure interaction system. Using Hamilton’s principle of the integral in time which is zero, or the results obtained in graphical form. Also found in work by Xu et al., ²¹ a Fluid-Structure Interaction (FSI) study applies the coupling that fluid force equals solid force, by models of finite elements of a bent pipe filled with fluid in free vibration at high frequency which meets the kinematic conditions on the interface of the ship by ADINA software. Shows that the results of FSI seriously affect the motor modal frequencies where the liquid reduces the motor frequencies and that the fluid-structure interaction hardly affects the liquid modal frequencies, if the fluid velocity drops towards stability then one has the vibrational continuity of driving and vibration disperses as time increases. Catinaccio, ²² studied the general behaviour of thin isotropic pipes, thick isotropic pipes and straight uniform orthotropic thin pipes with recessed ends, they considered that pipes under internal pressure had a comparison between Timoshenko’s formula ‘Theory of elastic stability’ and The longitudinal force on the pipe, therefore the overall compression/tension force P are approximations based on the small deformations and the incompressibility of the chosen material with the load lateral q. Use of finite elements in ANSYS. In this work, we find another notion of the coupling of the fluid-structure interaction where Küttler et al., ²³ gives an idea of the coupling where the incompressible fluid contains internal stresses defined by Navier-Stokes equation where the considered fluid is assumed to be Newtonian, its coupling depends on the distance between fluid and solid by the Dirichlet-Neumann iteration or it can be found in the article by Klöppel et al., ²⁴ which is based on Cartesian coordinates as a limit condition of the contract to obtain the possibility of applying the balloon deformed by time.

Mathematical formulation

Structure

Consideration is given to the fluid-structure interaction of a cylindrical pipe of circular section which transports fluid of diameter D, of internal radius R1, external R2 of length L, thickness h of a flow at internal pressure P over the entire length of the tube as shown in Figure 2.

Figure 2.

Cylindrical to internal pressure.

By relating the stresses to the deformations, the anisotropic, orthotropic and isotropic materials will be considered, except in the case of one-dimensional elements (axial, torsor and beam) where only the isotropic case is treated. The state of strain in an elastic body is defined by the strain components.

$$\left[\epsilon \right]=\left[\begin{array}{ccc}{\epsilon }_{rr}& {\gamma }_{r\theta }& {\gamma }_{rz}\\ {\gamma }_{\theta r}& {\epsilon }_{\theta \theta }& {\gamma }_{\theta z}\\ {\gamma }_{zr}& {\gamma }_{z\theta }& {\epsilon }_{zz}\end{array}\right]$$

(1)

If the displacement components in the directions of the axes are denoted by $(u_r,u_{\theta },u_z)$ , then the strain-displacement relations are ²⁵ :

$$\begin{array}{l}{\epsilon }_{rr}=\frac{\partial u_r}{\partial r},{\epsilon }_{\theta \theta }=\frac{1}{r}\left(\frac{\partial u_{\theta }}{\partial \theta }+u_r\right),{\epsilon }_{zz}=\frac{\partial u_z}{\partial z}\\ {\epsilon }_{r\theta }=\left(\frac{\partial u_{\theta }}{\partial r}-\frac{u_{\theta }}{r}+\frac{1}{r}\frac{\partial u_r}{\partial \theta }\right){\epsilon }_{rz}=\left(\frac{\partial u_r}{\partial z}+\frac{\partial u_z}{\partial r}\right)\\ {\epsilon }_{\theta z}=\left(\frac{\partial u_{\theta }}{\partial z}+\frac{1}{r}\frac{\partial u_z}{\partial \theta }\right)\end{array}$$

(2)

Consider a transversely isotropic elastic solid cylinder of constant section. The isotropic plane of the material is perpendicular to the axis of the cylinder, which is considered to be the same as the z-axis of the cylindrical coordinate system (r, θ, z), see Figure 3.

Figure 3.

Constraint in cylindrical coordinates.

The stress ( ${\sigma }_{\mathrm{ij}}$ )-displacement (u_i) constitutive relations for materials with transverse isotropy are. ²⁶

$$\begin{array}{c}\multicolumn{1}{c}{{\sigma }_{\mathrm{rr}}=C_{11}\frac{\partial u_r}{\partial r}+C_{12}(\frac{1}{r}\frac{\partial u_{\theta }}{\partial \theta }+\frac{u_r}{r})+C_{13}\frac{\partial u_z}{\partial z}}\\ \multicolumn{1}{c}{{\sigma }_{\theta \theta }=C_{12}\frac{\partial u_r}{\partial r}+C_{11}(\frac{1}{r}\frac{\partial u_{\theta }}{\partial \theta }+\frac{u_r}{r})+C_{13}\frac{\partial u_z}{\partial z}}\\ \multicolumn{1}{c}{{\sigma }_{\mathrm{zz}}=C_{13}\frac{\partial u_r}{\partial r}+C_{13}(\frac{1}{r}\frac{\partial u_{\theta }}{\partial \theta }+\frac{u_r}{r})+C_{33}\frac{\partial u_z}{\partial z}}\\ \multicolumn{1}{c}{{\sigma }_{\mathrm{rz}}=C_{44}(\frac{\partial u_z}{\partial r}+\frac{\partial u_r}{\partial z})}\\ \multicolumn{1}{c}{{\sigma }_{\theta z}=C_{44}(\frac{\partial u_z}{\partial r}+\frac{\partial u_r}{\partial z})}\\ \multicolumn{1}{c}{{\sigma }_{r\theta }=C_{66}(\frac{1}{r}\frac{\partial u_r}{\partial \theta }+\frac{\partial u_r}{\partial \theta }-\frac{u_r}{r})}\end{array}$$

(3)

where C₆₆ = (C₁₁−C₁₂)/2 and $C_{\mathrm{ij}}$ are the elastic constants. The equations of motion in cylindrical coordinates are:

$$\begin{array}{c}\frac{\partial {\sigma }_{\mathrm{rr}}}{\partial r}+\frac{1}{r}\frac{\partial {\sigma }_{r\theta }}{\partial \theta }+\frac{\partial {\sigma }_{\mathrm{rz}}}{\partial z}+\frac{1}{r}({\sigma }_{\mathrm{rr}}-{\sigma }_{\theta \theta })+F_r={\rho }_s\frac{{\partial }^2{u}_r}{\partial t^2}\hfill \\ \frac{\partial {\sigma }_{r\theta }}{\partial r}+\frac{1}{r}\frac{\partial {\sigma }_{\theta \theta }}{\partial \theta }+\frac{\partial {\sigma }_{\theta z}}{\partial z}+\frac{2}{r}{\sigma }_{r\theta }+F_{\theta }={\rho }_s\frac{{\partial }^2{u}_{\theta }}{\partial t^2}\hfill \\ \frac{\partial {\sigma }_{\mathrm{rr}}}{\partial r}+\frac{1}{r}\frac{\partial {\sigma }_{\theta z}}{\partial \theta }+\frac{\partial {\sigma }_{\mathrm{zz}}}{\partial z}+\frac{1}{r}{\sigma }_{\mathrm{rz}}+F_z={\rho }_s\frac{{\partial }^2{u}_z}{\partial t^2}\hfill \end{array}$$

(4)

Where ρ is the density and t are the time variable. Unlike the conventional displacement method ²⁷ which eliminates the six stress variables σ_ij from equations (3) and (4) to obtain three second-order coupled partial differential equations around the three displacement components u_i, the state space approach generally keeps three stress components and three displacement components and transforms the governing equations into a set of first-order ordinary differential equations concerning a particular coordinate variable. The state-space approach has particular advantages over the displacement method in solving certain kinds of mechanical and engineering problems, and the interested reader is referred to Chen et al. ²⁸ for more details and references. Where [D] is the matrix that links the displacements and deformations. The stress-strain relationship takes the Hooke form:

$$\left\{\sigma \right\}=\left[D\right]\left\{\epsilon \right\}$$

(5)

Where [D] is a symmetric matrix. An anisotropic material, it contains 21 independent material constants. In the case of an isotropic material, it is:

$$\begin{array}{l}\left[D\right]=\frac{E}{\left(1+\nu \right)\left(1-2\nu \right)}\\ \left[\begin{array}{cccccc}\left(1-\nu \right)& 0& 0& 0& 0& 0\\ 0& \left(1-\nu \right)& 0& 0& 0& 0\\ 0& 0& \left(1-\nu \right)& 0& 0& 0\\ 0& 0& 0& \frac{1}{2}\left(1-2\nu \right)& 0& 0\\ 0& 0& 0& 0& \frac{1}{2}\left(1-2\nu \right)& 0\\ 0& 0& 0& 0& 0& \frac{1}{2}\left(1-2\nu \right)\end{array}\right]\end{array}$$

(6)

Fluid

The basic equations of fluid dynamics are called Navier-Stokes ²⁹ equations. In the case of an isothermal flow, that is to say, a flow at a constant temperature, they represent two physical laws of conservation: the conservation of mass and the conservation of linear momentum. There are different ways to derive these equations. The flow will be described with the variables which are supposed to be sufficiently smooth functions in the time interval [0, T] and the domain Ω⊂R³.

• ρ(t,x): volumic mass [kg/m³],

• v(t,x) : velocity [m/s],

• P(t,x) : pressure [Pa = N/m³],

The Navier-Stokes (N-S) equations ³⁰ are the widely applied mathematical model to examine changes in these properties during dynamic and/or thermal interactions. The equations are adjustable according to the problems to be simulated and are expressed based on the principles of conservation of mass, momentum and energy:

• ➢ Conservation of mass: continuity equation.

• ➢ Conservation of momentum: Newton’s second law.

We consider an incompressible isothermal Newtonian ³¹ flow (density ρ = constant, viscosity μ = constant), with a velocity field $\overrightarrow{V}=(v_r,v_{\theta },v_z)$ , which gives us the Navier-Stokes equation as follows:

$$\rho (\frac{\partial v_r}{\partial t}+v\frac{\partial v_r}{\partial r}+\frac{v_{\theta }}{r}\frac{\partial v_r}{\partial \theta }-\frac{v_{\theta }^2}{r}+v_r\frac{\partial v_r}{\partial z})=-\frac{\partial p}{\partial r}+\rho g_r$$

(7)

$$\rho (\frac{\partial v_{\theta }}{\partial t}+v_r\frac{\partial v_{\theta }}{\partial r}+\frac{v_{\theta }}{r}\frac{\partial v_{\theta }}{\partial \theta }+\frac{v_r{v}_{\theta }}{r}+v_z\frac{\partial v_{\theta }}{\partial z})=-\frac{1}{r}\frac{\partial p}{\partial \theta }+\rho g_{\theta }$$

(8)

$$\rho (\frac{\partial v_z}{\partial t}+v_r\frac{\partial v_z}{\partial r}+\frac{v_{\theta }}{r}\frac{\partial v_z}{\partial \theta }+v_z\frac{\partial v_z}{\partial z})=-\frac{\partial p}{\partial z}+\rho g_z$$

(9)

It was also called the mixed finite element method. To present how to deal with this difficulty, let us recall that the pressure acts as a Lagrangian multiplier of the incompressibility constraint. Then, the algebraic system for the nodal values of velocity and pressure in a Galerkin formulation ³² will be governed by a partitioned matrix with a zero submatrix on the diagonal. The resolution of the algebraic system depends on an appropriate choice of finite element spaces for velocity and pressure. To better understand how to deal with the difficulty of the incompressibility condition, the procedure of the finite element method based on the formulation of primitive variables will be presented step by step. We will start with the discretization of the incompressible Navier-Stokes equation in weak form using the weak formulation. Next, the selection of finite element spaces for velocity and pressure will be discussed.

The application of the Galerkin method (PTV). ³³

$$w={\int }_{\mathrm{\Omega }}\left(\mathrm{\Psi }f\left(v\right)+F\right)d\mathrm{\Omega }=0$$

(10)

$$\mathrm{w}={\sum }_i^e\mathrm{avec}w=w_{N-S}+w_{\mathrm{continuit}é}$$

(11)

Discretization

The eight-node hexahedron element is a generalization of the bilinear quadrilateral element. Under a tri-linear transformation, the reference element (unit cube; Figure 4) defined in the space η, ζ, ξ is the image of the domain Ω_e defined in the space x, y and z. The basis of the vector space generated by polynomials of degree 1 is {1, ξ, η, ζ, ξη, ηζ, ξζ, ξηζ}. By using the properties of the shape functions, we deduce the expressions of the shape functions.

Figure 4.

Tri-linear hexahedron element.

From the shape, we deduce the expressions of the shape functions. Spatial functions for the HFEM approximation of the velocity component become ³⁴ :

$$\begin{array}{c}{N}_1=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$4$}\right.\left(1-\xi \right)\left(1-\eta \right)\\ \begin{array}{c}{N}_2=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$4$}\right.\left(1+\xi \right)\left(1-\eta \right)\\ N_3=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$4$}\right.\left(1+\xi \right)\left(1+\eta \right)\\ \begin{array}{c}{N}_4=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$4$}\right.\left(1-\xi \right)\left(1+\eta \right)\\ N_5=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\left(1-\xi \right)\times \sum _{s=3}^p{g}_s\left(\eta \right)\\ \begin{array}{c}{N}_6=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\sum _{s=3}^p{g}_s\left(\xi \right)\times \left(1-\eta \right)\\ \begin{array}{c}{N}_7=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\left(1+\xi \right)\times \sum _{s=3}^p{g}_s\left(\eta \right)\\ \begin{array}{c}{N}_8=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\sum _{s=3}^p{g}_s\left(\xi \right)\times \left(1+\eta \right)\\ \begin{array}{l}{N}_9=\sum _{s=3}^p{g}_s\left(\xi \right)\times \sum _{s=3}^p{g}_s\left(\eta \right)\\ (continuityofthevariationoftheelementinthecenter)\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}$$

(12)

Spatial functions for pressure approximation:

$$\begin{array}{c}{M}_1=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$4$}\right.\left(1-\xi \right)\left(1-\eta \right)\\ \begin{array}{c}{M}_2=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$4$}\right.\left(1+\xi \right)\left(1-\eta \right)\\ \begin{array}{c}{M}_3=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$4$}\right.\left(1+\xi \right)\left(1+\eta \right)\\ M_4=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$4$}\right.\left(1-\xi \right)\left(1+\eta \right)\end{array}\end{array}\end{array}$$

(13)

Cylindrical vibration by the HFEM method

From the conditions of the choice the finite element, whatever is in 2D or 3D, applies the hierarchical finite element (HFEM) of it replaces the nodes 5, 6, 7 and 8 of the element of Taylor-Hood by polynomial bar following an interpolation to be established by the legend function ³⁵ which is defined as follows:

$$f_m\left(\xi \right)=\frac{1}{m!{(-2)}^m}\frac{d^m}{d{\xi }^m}\left\{{(1-{\xi }^2)}^m\right\}-1\le \xi \le +1$$

(14)

Which satisfies the orthogonality relation ${\int }_{-1}^1{\mathrm{f}}_{\mathrm{i}}\left(\xi \right){\mathrm{f}}_{\mathrm{j}}\left(\xi \right)\mathrm{d}\mathrm{\xi }=\frac{2}{2\mathrm{i}+1}{\delta }_{\mathrm{ij}}$ where δij is the Kronecker delta function. The most commonly used hierarchical functions come from the polynomials of Legendre Nguyen et al. ³ namely:

$$f_m^r\left(\xi \right)=\sum _{n=0}^{m/2}\frac{{\left(-1\right)}^n}{n!2^n}\frac{\left(2m-2n-2s-1\right)!!}{\left(m-2n\right)!}{\xi }^{m-2n}$$

(15)

With

$$m!!=m(m-2)\cdots (2\mathrm{ou}1),0!!=1(-1)=1!!$$

The orthogonality of the polynomial of Legendre Groen et al. ³⁶ in the interval [−1 1], is made by integration, and all their lower order derivatives disappear at ξ = ±1, we obtain the following result numerical formula Babuška et al. ³⁷

The H-FEM can be used to calculate the vibration characteristics, point displacements of a cylindrical shell or the tangential components v_r, v_θ, v_z, of the normal component. The routes are given by:

$$|\begin{array}{c}{v}_r(\xi ,\eta ,\zeta )=\sum _{r=1}^q\sum _{s=1}^p\sum _{t=1}^l{g}_r\left(\xi \right)g_s\left(\eta \right)g_l\left(\eta \right)U_{\mathrm{rst}}\hfill \\ \multicolumn{1}{c}{v_{\theta }(\xi ,\eta ,\zeta )=\sum _{r=1}^q\sum _{s=1}^p\sum _{t=1}^l{g}_r\left(\xi \right)g_s\left(\eta \right)g_l\left(\eta \right)V_{\mathrm{rst}}}\\ \multicolumn{1}{c}{v_z(\xi ,\eta ,\zeta )=\sum _{r=1}^q\sum _{s=1}^p\sum _{t=1}^l{g}_r\left(\xi \right)g_s\left(\eta \right)g_l\left(\eta \right)W_{\mathrm{rst}}}\end{array},$$

(16)

Couplings between fluid and structure FEM-FEM penalty method

The non-zero components of matrix B are displacements of the gradients ξ or η, in which the directions of the nodes of displacement its resultant in the tractions and displacements relating to the central point of integration. In this class of interface elements, displacements relative to the n pairs of nodes are used instead of a field of relative displacement interpolated in integration points. The element stiffness matrix is now expanded as:

$$\mathrm{K}=\sum _{i=1}^{n_i}{A}_i{\left({\mathrm{B}}_{\mathrm{d}}^{\mathrm{T}}\right)}_i{\left({\mathrm{D}}_{\mathrm{d}}\right)}_i{\left({\mathrm{B}}_{\mathrm{d}}\right)}_i$$

(17)

where the matrix ${\left({\mathrm{B}}_{\mathrm{d}}\right)}_i$ for the element of the interface to the set of nodes i is written

$${\left({\mathrm{B}}_{\mathrm{d}}\right)}_i=\left[\begin{array}{ccc}\begin{array}{cc}-1\hfill & +1\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\hfill & \begin{array}{cc}0\hfill & 0\hfill \\ \multicolumn{1}{c}{-1}& +1\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\hfill & \begin{array}{cc}0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \\ \multicolumn{1}{c}{-1}& +1\hfill \end{array}\hfill \end{array}\right]$$

(18)

and A_i, the i-node surface. The sub-matrix K_sou for a linear interface element with four pairs of nodes and a surface A becomes:

$${\mathrm{K}}_{sou}=\frac{{\mathrm{Ad}}_{sou}}{4}\left[\begin{array}{cc}\begin{array}{c}\begin{array}{cc}1& 0\\ 0& 1\end{array}\\ \begin{array}{c}\begin{array}{cc}0& 0\\ 0& 0\end{array}\\ \begin{array}{c}\begin{array}{cc}-1& 0\\ 0& -1\end{array}\\ \begin{array}{cc}0& 0\\ 0& 0\end{array}\end{array}\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{cc}0& 0\\ 0& 0\end{array}\\ \begin{array}{c}\begin{array}{cc}1& 0\\ 0& 1\end{array}\\ \begin{array}{c}\begin{array}{cc}0& 0\\ 0& 0\end{array}\\ \begin{array}{cc}-1& 0\\ 0& -1\end{array}\end{array}\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{cc}-1& 0\\ 0& -1\end{array}\\ \begin{array}{c}\begin{array}{cc}0& 0\\ 0& 0\end{array}\\ \begin{array}{c}\begin{array}{cc}1& 0\\ 0& 1\end{array}\\ \begin{array}{cc}0& 0\\ 0& 0\end{array}\end{array}\end{array}\end{array}& \begin{array}{c}\begin{array}{cc}0& 0\\ 0& 0\end{array}\\ \begin{array}{c}\begin{array}{cc}-1& 0\\ 0& -1\end{array}\\ \begin{array}{c}\begin{array}{cc}0& 0\\ 0& 0\end{array}\\ \begin{array}{cc}1& 0\\ 0& 1\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right]$$

(19)

Bien évidemment, aucun couplage n’existe entre les degrés de liberté des différentes paires de nœuds. On note que les éléments d’interface à intégration localisée sont proches des éléments d’interface dits nodaux ou ponctuels (Ngo et Scordelis 1967), qui peuvent, au sens mécanique, être interprétés comme des interfaces reliées par des ressorts aux nœuds.

The FEM-FEM penalty method ³⁸ is characterized, on the one hand, by the link that exists between one the structure element and one fluid element, and, on the other, by the stiffness at the interface; the radial force at the interface can therefore be determined. In the present work, the stiffness is represented by a bar in order to take into account the axial force at the interface, in addition to the radial force, as clearly shown in Figure 5.

Figure 5.

FEM- FEM spring coupling and bar coupling.

The method used consists in arranging dummy bars, at the interface, capable of generating radial and axial forces penetrating the contact surface, for all the nodes. The equilibrium position of these bars corresponds to a slave node positioned on the master segment. The forces are determined by:

$${{\mathrm{F}}_r}_{(\mathrm{m}->\mathrm{s})}=-{\mathrm{k}}_r({{\mathrm{u}}_r}_{\mathrm{s}}-{{\mathrm{u}}_z}_{\mathrm{f}})$$

(20)

$${\mathrm{F}}_{z\left(\mathrm{m}->\mathrm{s}\right)}=-{\mathrm{k}}_z\left({\mathrm{u}}_{r\mathrm{s}}-{\mathrm{u}}_{z\mathrm{f}}\right)$$

(21)

This approach makes it possible to determine the equation of motion as follows:

$$\begin{array}{c}\multicolumn{1}{c}{\{\begin{array}{c}\left[M_{\mathrm{sf}}\right]\left\{\stackrel{··}{\mathrm{u}}\right\}+\left[K_{\mathrm{sf}}\right]{\delta }_{\mathrm{n}}=\left\{f\right\}\mathrm{si}{\delta }_{\mathrm{n}}<0\hfill \\ \multicolumn{1}{c}{\left[M_{\mathrm{sf}}\right]\{\stackrel{··}{\mathrm{u}}\}=\left\{f\right\}\mathrm{si}{\delta }_{\mathrm{n}}\ge 0}\end{array}}\end{array}$$

(22)

The boundary conditions are defined with respect to the solid and the fluid. The boundary conditions in the case of the solid are:

$$-\left[{\sigma }_{\mathrm{s}}\right]\left\{\mathrm{n}\right\}+\left\{{\mathrm{f}}_{\mathrm{s}}\right\}=0;\mathrm{sur}{\mathrm{S}}_{\mathrm{s}}\iff \left[{\sigma }_{\mathrm{s}}\right]\left\{\mathrm{n}\right\}=\left\{{\mathrm{f}}_{\mathrm{s}}\right\}$$

(23)

The boundary conditions in the case of the fluid are:

$$-\left[{\sigma }_{\mathrm{f}}\right]\left\{\mathrm{n}\right\}+\left\{{\mathrm{f}}_{\mathrm{f}}\right\}=0;\mathrm{sur}{\mathrm{S}}_{\mathrm{f}}\iff \left[{\sigma }_{\mathrm{f}}\right]\left\{\mathrm{n}\right\}=\left\{{\mathrm{f}}_{\mathrm{f}}\right\}$$

(24)

Since the interaction occurs between the fluid and the structure, the action-reaction law imposes the equality of the two expressions (23) and (24).

Therefore, one can write ³⁹ :

$$\left[{\sigma }_{\mathrm{s}}\right]\left\{\mathrm{n}\right\}=\left[{\sigma }_{\mathrm{f}}\right]\left\{\mathrm{n}\right\}$$

(25)

And

$$\left\{{\mathrm{f}}_{\mathrm{s}}\right\}=\left\{{\mathrm{f}}_{\mathrm{f}}\right\}$$

(26)

The contact forces in the penalty case are determined by:

$$\left\{{\mathrm{f}}_{\mathrm{s}}\right\}={\mathrm{\alpha }}_{\mathrm{P}}\left[{\mathrm{g}}_{\mathrm{s}}\right]\mathrm{et}\left\{{\mathrm{f}}_{\mathrm{f}}\right\}={\mathrm{\alpha }}_{\mathrm{P}}\left[{\mathrm{g}}_{\mathrm{f}}\right]$$

(27)

α_p is the penalty factor.

To determine the quantities $\left[{\mathrm{g}}_{\mathrm{s}}\right]$ and $\left[{\mathrm{g}}_{\mathrm{f}}\right]$ , one must first find the stiffness matrix for the bar; this can be modeled by the hierarchical finite element method (HFEM). The shape functions for displacements ${{\mathrm{u}}_r}_1,{{\mathrm{u}}_r}_2,{{\mathrm{u}}_z}_1,{{\mathrm{u}}_z}_2$ and rotations ${\theta }_1,{\theta }_2$ at the nodes of the bar are given by the following relation

$$\begin{array}{c}\multicolumn{1}{c}{\{\begin{array}{c}\begin{array}{c}{f}_1\left(\xi \right)=\frac{1+\xi }{2}\hfill \\ \multicolumn{1}{c}{f_2\left(\xi \right)=\frac{1}{2}-\frac{3}{4}\xi +\frac{1}{4}{\xi }^3}\end{array}\hfill \\ \multicolumn{1}{c}{\begin{array}{c}{f}_3\left(\xi \right)=(\frac{1}{8}-\frac{1}{8}\xi -\frac{1}{8}{\xi }^2+\frac{1}{8}{\xi }^3)l_e\hfill \\ \multicolumn{1}{c}{f_4\left(\xi \right)=\frac{1-\xi }{2}}\end{array}}\\ \multicolumn{1}{c}{\begin{array}{c}{f}_5\left(\xi \right)=\frac{1}{2}+\frac{3}{4}\xi -\frac{1}{4}{\xi }^3\hfill \\ \multicolumn{1}{c}{f_6\left(\xi \right)=(-\frac{1}{8}-\frac{1}{8}\xi +\frac{1}{8}{\xi }^2+\frac{1}{8}{\xi }^3)l_e}\end{array}}\end{array}}\end{array}$$

(28)

To check the penalty coupling condition, ⁴⁰ it is necessary that the ratio $\frac{\mathrm{EA}}{L}\approx 1$ for the radial movement, and $\frac{8}{L^3}\mathrm{EI}\approx 1$ for the axial movement. Note that this approach makes it possible to determine the coupling matrix^41,42 as follows:

$$\begin{array}{c}{\left[\mathrm{K}\right]}_{\mathrm{T}}={\int }_0^L\mathrm{EA}{\left[B\right]}^T\left[B\right]\mathrm{dx}=\frac{1}{L}{\int }_{-1}^1\mathrm{EA}{\left[\begin{array}{c}-1\hfill \\ \multicolumn{1}{c}{1}\end{array}\right]}^T\left[\begin{array}{cc}-1\hfill & 1\hfill \end{array}\right]d\xi \frac{\mathrm{EA}}{L}\hfill \\ \left[\begin{array}{cc}1\hfill & -1\hfill \\ \multicolumn{1}{c}{-1}& 1\hfill \end{array}\right]={\mathrm{\alpha }}_{\mathrm{T}}\left[{\mathrm{G}}_{\mathrm{Te}}\right]={\mathrm{\alpha }}_{\mathrm{P}}\left[{\mathrm{g}}_{\mathrm{Ts}}\right]\hfill \\ {\left[\mathbf{K}\right]}_{\mathbf{F}}={\int }_0^{\mathit{L}}{\left[\mathit{B}\right]}^{\mathit{T}}\left[\mathit{D}\right]\left[\mathit{B}\right]\mathit{dx}={\int }_{-1}^1\frac{4}{{\mathit{L}}^2}{\left[\mathit{B}\right]}^{\mathit{T}}\left[\mathit{D}\right]\frac{4}{{\mathit{L}}^2}\left[\mathit{B}\right]\frac{\mathit{L}}{2}\mathit{d}\xi \hfill \\ {\left[\mathbf{K}\right]}_{\mathbf{F}}=\frac{10\mathit{L}}{12}\times \frac{8}{{\mathit{L}}^3}\mathit{EI}{\int }_{-1}^1\hfill \\ \left[\begin{array}{cc}\begin{array}{cc}{\left({\mathit{N}}_{1,\xi \xi }\right)}^2& \frac{\mathit{L}}{2}{\mathit{N}}_{1,\xi \xi }{\mathit{N}}_{2,\xi \xi }\\ & \frac{{\mathit{L}}^2}{4}{\left({\mathit{N}}_{2,\xi \xi }\right)}^2\end{array}& \begin{array}{cc}{\mathit{N}}_{1,\xi \xi }{\mathit{N}}_{3,\xi \xi }& \frac{\mathit{L}}{2}{\mathit{N}}_{1,\xi \xi }{\mathit{N}}_{4,\xi \xi }\\ \frac{\mathit{L}}{2}{\mathit{N}}_{2,\xi \xi }{\mathit{N}}_{3,\xi \xi }& \frac{{\mathit{L}}^2}{4}{\mathit{N}}_{2,\xi \xi }{\mathit{N}}_{4,\xi \xi }\end{array}\\ \mathit{S}\mathit{Y}\mathit{M}& \begin{array}{cc}{\left({\mathit{N}}_{3,\xi \xi }\right)}^2& \frac{\mathit{L}}{2}{\mathit{N}}_{3,\xi \xi }{\mathit{N}}_{4,\xi \xi }\\ & \frac{{\mathit{L}}^2}{4}{\left({\mathit{N}}_{4,\xi \xi }\right)}^2\end{array}\end{array}\right]\mathit{d}\xi \hfill \\ {\left[\mathrm{K}\right]}_{\mathrm{F}}=\frac{10L}{12}\times \frac{8}{L^3}\mathrm{EI}\left[\begin{array}{cc}1\hfill & -1\hfill \\ \multicolumn{1}{c}{-1}& 1\hfill \end{array}\right]={\mathrm{\alpha }}_{\mathrm{F}}\left[{\mathrm{G}}_{\mathrm{Fe}}\right]={\mathrm{\alpha }}_{\mathrm{P}}\left[{\mathrm{g}}_{\mathrm{Fs}}\right]\hfill \end{array}$$

(29)

The multiplication of the matrices ${\left[\mathrm{K}\right]}_{\mathrm{T}}\mathrm{by}\frac{\mathrm{EA}}{L}\mathrm{and}{\left[\mathrm{K}\right]}_{\mathrm{F}}$ by $\frac{10L}{12}\times \frac{8}{L^3}\mathrm{EI}$ makes it possible to obtain the stiffness matrix of the structure-fluid assembly. The coupling can then be represented as:

$$\left[{\mathrm{K}}_{\mathrm{SYS}}\right]=\left[{\mathrm{K}}_{\mathrm{FS}}\right]+{\mathrm{\alpha }}_{\mathrm{P}}\underset{{\mathrm{G}}_{\mathrm{sf}}}{\underbrace{\left[\underset{g_{Ts}}{\underbrace{{\mathrm{\alpha }}_{\mathrm{T}}\left[{\mathrm{G}}_{Te}\right]}}\right]+{\mathrm{\alpha }}_{\mathrm{P}}\left[\underset{g_{Ts}}{\underbrace{{\mathrm{\alpha }}_{\mathrm{F}}\left[{\mathrm{G}}_{Fe}\right]}}\right]}}$$

(30)

whith:

$$\left[\mathrm{G}\right]={\mathrm{\alpha }}_{\mathrm{P}}\left[{\mathrm{G}}_{\mathrm{sf}}\right]$$

(31)

The matrix $\left[{\mathrm{g}}_{\mathrm{sf}}\right]$ is then multiplied by its transposed matrix to obtain a square matrix.

$$\left[\mathrm{G}\right]={\mathrm{\alpha }}_{\mathrm{P}}\times [\left[{\mathrm{g}}_{\mathrm{sf}}\right]\times {\left[{\mathrm{g}}_{\mathrm{sf}}\right]}^{\mathrm{T}}]={\left[\mathrm{g}\right]}^2$$

(32)

It should be noted that α_p must be between 10⁴ and 10⁹; $\left[{\mathrm{K}}_{\mathrm{FS}}\right]$ is the matrix of global stiffness of the fluid-structure assembly.

Equation of motion

The potential energy of the structure is given by:

$$\delta {\mathrm{U}}_{\mathrm{S}}=\frac{1}{2}\int {\left\{\delta {\epsilon }_{\mathrm{S}}\right\}}^{\mathrm{T}}\left\{{\sigma }_{\mathrm{S}}\right\}\mathrm{dV}$$

(33)

The expression (33) for the cylindrical shell, in cylindrical coordinates, may be expressed as:

$$\delta {\mathrm{U}}_{\mathrm{S}}=\frac{1}{2}{\int }_{{\mathrm{r}}_1}^{{\mathrm{r}}_2}{\int }_0^{\mathrm{H}}{\int }_0^{2\pi }\left(\delta {\left\{\begin{array}{c}\begin{array}{c}{\epsilon }_{\mathrm{rr}}\hfill \\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{zz}}}\\ \multicolumn{1}{c}{{\epsilon }_{\theta \theta }}\end{array}\hfill \\ \multicolumn{1}{c}{\begin{array}{c}{\gamma }_{\mathrm{rz}}\hfill \\ \multicolumn{1}{c}{{\gamma }_{\mathrm{r}\mathrm{\theta }}}\\ \multicolumn{1}{c}{{\gamma }_{\mathrm{z}\mathrm{\theta }}}\end{array}}\end{array}\right\}}^{\mathrm{T}}\left[\mathrm{D}\right]\left\{\begin{array}{c}\begin{array}{c}{\epsilon }_{\mathrm{rr}}\hfill \\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{zz}}}\\ \multicolumn{1}{c}{{\epsilon }_{\theta \theta }}\end{array}\hfill \\ \multicolumn{1}{c}{\begin{array}{c}{\gamma }_{\mathrm{rz}}\hfill \\ \multicolumn{1}{c}{{\gamma }_{\mathrm{r}\mathrm{\theta }}}\\ \multicolumn{1}{c}{{\gamma }_{\mathrm{z}\mathrm{\theta }}}\end{array}}\end{array}\right\}\right)\mathrm{rdrdzd}\mathrm{\theta }$$

(34)

For a solid to be in equilibrium, the work of the external forces must be equal to the kinetic energy:

$$\mathrm{\delta }\mathrm{T}=\mathrm{\delta }\mathrm{W}$$

(35)

The elementary kinetic energy is given by:

$$\mathrm{\delta }\mathrm{T}=\delta \left({\int }_{{\mathrm{r}}_1}^{{\mathrm{r}}_2}{\int }_0^{\mathrm{H}}{\int }_0^{2\pi }({\left(\frac{\mathrm{du}}{\mathrm{dt}}\right)}^2+{\left(\frac{\mathrm{dv}}{\mathrm{dt}}\right)}^2+{\left(\frac{\mathrm{dw}}{\mathrm{dt}}\right)}^2)\right)$$

(36)

The application of the principle of Hamilton ⁴² allows writing:

$${\int }_{{\mathrm{t}}_1}^{{\mathrm{t}}_2}(\mathrm{\delta }\mathrm{T}-\mathrm{\delta }\mathrm{U})\mathrm{dt}=0$$

(37)

This then makes it possible to obtain the equation of movement of the structure.

$$({\left[\mathrm{M}\right]}_{\mathrm{s}}^{\mathrm{e}}\left\{\stackrel{··}{\mathrm{U}}\right\}+{\left[\mathrm{K}\right]}_{\mathrm{s}}^{\mathrm{e}}\left\{\mathrm{U}\right\}+{\left\{\mathrm{F}\right\}}_{\mathrm{s}}^{\mathrm{e}})=0$$

(38)

The stresses due to the displacement of the fluid, ⁴³ as determined by the Navier-Stokes equation, are given by the following relation:

$$\begin{array}{c}\left\{\sigma \left(M\right)\right\}=2\mu \hfill \\ \left\{\begin{array}{c}\begin{array}{c}\frac{\partial {\stackrel{·}{u}}_r}{\partial r}-\frac{p}{2\mu }+\frac{1}{2}(\frac{1}{r}\frac{\partial {\stackrel{·}{u}}_{\theta }}{\partial \theta }-\frac{{\stackrel{·}{u}}_{\theta }}{r}+\frac{1}{r}\frac{\partial {\stackrel{·}{u}}_r}{\partial \theta })+\frac{1}{2}(\frac{\partial {\stackrel{·}{u}}_z}{\partial z}+\frac{\partial {\stackrel{·}{u}}_r}{\partial r})\hfill \\ \multicolumn{1}{c}{\frac{1}{2}(\frac{1}{r}\frac{\partial {\stackrel{·}{u}}_{\theta }}{\partial \theta }-\frac{{\stackrel{·}{u}}_{\theta }}{r}+\frac{1}{r}\frac{\partial {\stackrel{·}{u}}_r}{\partial \theta })+(\frac{1}{r}\frac{\partial {\stackrel{·}{u}}_{\theta }}{\partial \theta }+\frac{{\stackrel{·}{u}}_r}{r})-\frac{p}{2\mu }+\frac{1}{2}(\frac{1}{r}\frac{\partial {\stackrel{·}{u}}_z}{\partial \theta }+\frac{\partial {\stackrel{·}{u}}_{\theta }}{\partial z})}\\ \multicolumn{1}{c}{\frac{1}{2}(\frac{\partial {\stackrel{·}{u}}_r}{\partial z}+\frac{\partial {\stackrel{·}{u}}_z}{\partial r})+\frac{1}{2}(\frac{1}{r}\frac{\partial {\stackrel{·}{u}}_z}{\partial \theta }+\frac{\partial {\stackrel{·}{u}}_{\theta }}{\partial z})+\frac{\partial {\stackrel{·}{u}}_z}{\partial z}-\frac{p}{2\mu }}\end{array}\hfill \\ \multicolumn{1}{c}{\begin{array}{c}\frac{1}{2}(\frac{1}{r}\frac{\partial {\stackrel{·}{u}}_{\theta }}{\partial \theta }-\frac{{\stackrel{·}{u}}_{\theta }}{r}+\frac{1}{r}\frac{\partial {\stackrel{·}{u}}_r}{\partial \theta })\hfill \\ \multicolumn{1}{c}{\frac{1}{2}(\frac{\partial {\stackrel{·}{u}}_r}{\partial z}+\frac{\partial {\stackrel{·}{u}}_z}{\partial r})}\\ \multicolumn{1}{c}{\frac{1}{2}(\frac{1}{r}\frac{\partial {\stackrel{·}{u}}_z}{\partial \theta }+\frac{\partial {\stackrel{·}{u}}_{\theta }}{\partial z})}\end{array}}\end{array}\right\}\hfill \end{array}$$

(39)

The equation of motion of the fluid is given by ⁴⁴ :

$${\left[\mathrm{M}\right]}_{\mathrm{f}}^{\mathrm{e}}\left\{\stackrel{·}{\mathrm{U}}\right\}+{\left[\mathrm{K}\right]}_{\mathrm{f}}^{\mathrm{e}}\left\{\mathrm{U}\right\}=0$$

(40)

The Navier-Stokes equation (38) contains strong variational terms. The Galerkin method is applied in order to weaken the formulation for the numerical discretization. In addition, both equations (37)and (39) are used to find the equation of motion for the coupled system in the following matrix form:

$${\left[\begin{array}{ccc}{M}_s\hfill & M_s\hfill & 0\hfill \\ \multicolumn{1}{c}{M_s}& M_{\mathrm{sf}}\hfill & M_f\hfill \\ \multicolumn{1}{c}{0}& M_f\hfill & M_f\hfill \end{array}\right]}_G{\left\{\stackrel{··}{U}\right\}}_G+{\left[\begin{array}{ccc}{K}_s\hfill & K_s\hfill & 0\hfill \\ \multicolumn{1}{c}{K_s}& K_{\mathrm{sf}}\hfill & K_f\hfill \\ \multicolumn{1}{c}{0}& K_f\hfill & K_f\hfill \end{array}\right]}_G{\left\{U\right\}}_G=0$$

(41)

The equation of motion of the entire system is determined in terms of the radial, circumferential and axial displacements in the cylindrical coordinate system (r, θ, z). The free vibrations of the components $u_r$ , $u_{\theta }$ and $u_z$ may be written in the form^45,46:

$$U=\{\begin{array}{c}{u}_r(\mathrm{r},\theta ,z,\mathrm{t})=\sum _1^m\sum _1^n{u}_r\cos \left(\mathrm{m}\mathrm{\theta }\right)·\cos \left(\lambda z\right)·\cos \mathrm{wt}\hfill \\ \multicolumn{1}{c}{u_{\theta }(\mathrm{r},\mathrm{\theta },z,\mathrm{t})=\sum _1^m\sum _1^n{u}_{\theta }\sin \left(\mathrm{m}\mathrm{\theta }\right)·\sin \left(\lambda \mathrm{rx}\right)·\cos \mathrm{wt}}\\ \multicolumn{1}{c}{u_z(\mathrm{r},\mathrm{\theta },z,\mathrm{t})=\sum _1^m\sum _1^n{u}_z\cos \left(\mathrm{m}\mathrm{\theta }\right)·\sin \left(\lambda x\right)·\cos \mathrm{wt}}\end{array}$$

Where $u_r$ , $u_{\theta }$ and $u_z$ are the displacements along the three cylindrical directions.

The coefficient λ is determined from the boundary conditions:

• - Cylinder supported at both ends $\lambda =n\pi \frac{R}{L}$

• - Cylinder whose ends are embedded $\lambda =(n+0.3)\pi \frac{R}{L}$

Validation

A programme, which allows obtaining the natural frequencies of the pressurized cylinder with a fluid flow, was developed. It was validated by studying its convergence and the results obtained were compared with those reported by other researchers.

Convergence

The study of convergence is clearly shown in Figures 6 and 7.

Figure 6.

Convergence for different degrees of the polynomial: Case of the solid.

Figure 7.

Convergence for different degrees of the polynomial: Case of the fluid.

Comparison of results

The eigenfrequencies and the eigenmodes of vibration of the pipes containing the fluid are determined. The physical and geometrical parameters of the pipe are given in the following Table 1 of the system clamped–clamped boundary.

Table 1.

Physical and geometric parameters of the pipe.

E	$\nu$	L/R	R/h	${\rho }_s$	${\rho }_f/{\rho }_s$
2.0e11 N/m2	0.32	6.0	666	7916 kg/m3	$0.00136$

The eigen frequencies of the pipe are determined as a function of the pressure variation for the first longitudinal mode and for different circumferential modes. The results obtained are compared with the experimental findings ⁴⁷ and those found numerically by the FEM ⁴⁸ and Zhang et al. ⁴⁹ ; they are given in Table 2 and illustrated in Figures 8 to 10. The percentages of errors $\epsilon {\%}_{ij}$ are defined in the tables by $\epsilon {\%}_{ij}=\left|\left(100\times \left(f_i-f_j\right)/f_i\right)\right|$ and Figure 11.

Table 2.

Validation of the model.

Mode		P 0 (KPa)	Frequency (Hz)					Relative diffirence (%)
m	n		Experimental 47	Provided FEM 48	Provided FEM 49	Present HFEM Spring	Present HFEM Bar	Provided FEM 48	Provided FEM 49	Present HFEM Spring	Present HFEM Bar
1	7	20.68	441	457.3	471.3	447.82	420.3	3.70	6.87	1.55	4.69
		27.57	497	525.9	542.4	518.67	485.01	5.81	9.13	4.36	2.41
		34.47	551	585.5	603.9	577.72	543.08	6.26	9.60	4.85	1.44
		41.36	606	640.7	660.9	633.26	590.2	5.73	9.06	4.50	2.61
		48.26	652	690.6	712.5	671.51	634.77	5.92	9.28	2.99	2.64
		62.05	740	780.4	805.3	763.65	729.65	5.46	8.82	3.20	1.40
		68.94	775	821.26	847.58	807.2	760.7	5.97	9.37	4.15	1.85
1	8	13.78	419	432.4	442.6	420.96	416.05	3.20	5.63	0.47	0.70
		20.68	510	526.3	539	520.89	508.54	3.20	5.69	2.14	0.29
		27.57	564	605.8	620.3	583.62	560.08	7.41	9.98	3.48	0.70
		34.47	633	675.6	691.9	662.75	624.33	6.73	9.30	4.70	1.37
		41.36	696	738.6	756.6	724.63	686.35	6.12	8.71	4.11	1.39
		48.26	751	796.4	815.8	783.07	740.91	6.05	8.63	4.27	1.34
		62.05	847	900.4	922.4	888.66	837.42	6.30	8.90	4.92	1.13
		68.94	892	947.8	971	934.23	881.68	6.26	8.86	4.73	1.16
1	9	20.68	569	596.4	607.6	583.08	574.59	4.82	6.78	2.47	0.98
		34.47	717	764.3	780	750.51	724.74	6.60	8.79	4.67	1.08
		41.36	793	836.8	852.9	827.65	805.74	5.52	7.55	4.37	1.61
		62.05	971	1020.5	1040.3	1009.9	980.2	5.10	7.14	4.01	0.95
		68.94	1017	1074.3	1095.2	1060.4	1025.7	5.63	7.69	4.27	0.86

Figure 8.

Variation of frequency versus flow pressure, for m = 1and n = 7.

Figure 9.

Variation of frequency versus flow pressure, for m = 1 and n = 8.

Figure 10.

Variation of frequency versus flow pressure, for m = 1 and n = 9.

Figure 11.

Relative error $\mathit{\epsilon }{\%}_{\mathit{i}\mathit{j}}$ of Table 2.

Figures 8 to 10 clearly show the similarity between the results obtained from a bar coupling model and the experimental ones.

Case study

After the validation of our work in being able to establish a graphic study of different geometrical parameters according to the applied coupling (bar and spring) with a frequency graphic presentation of Figures 11 to 16, we use different circumferential and longitudinal modes on the structure and fluid as we show on the theory part of the equation (42). This study allows us to have a general idea about the behaviour of the coupling between the natural frequencies as a function of the length and the radius of the pipe, as it is followed by a graphic appearance between the natural frequencies of interface (coupling) stiffness and bar or finds that the natural frequencies of stiffness rise than the natural frequencies of the bar with increasing pipe length, but when changing radius, the reverse occurs where the natural frequencies of bar rise than the natural frequencies of cylinder stiffness. We allow saying that the coupling by bar gives results close to the reality of the natural frequencies that coupling by stiffness which absorbs the applied forces of the fluid which causes damping for frequential (minimizes the frequencies to its real value) therefore we can say that the interface of the bar and perfect for coupling between solid and frequential fluid. As shown in Figures 12 to 17.

Figure 12.

Variation of the eigenfrequencies of the bar coupling as a function of the length for the different longitudinal modes or the constant circumferential modes m = 1.

Figure 13.

Variation of the eigenfrequencies of the stiffness coupling as a function of the length for the different longitudinal modes or the constant circumferential modes m = 1.

Figure 14.

Comparison of frequency variations versus pipe length for the cases spring coupling and bar coupling with m = 1.

Figure 15.

Variation of the eigen frequencies coupling stiffness versus pipe radius with m = 1.

Figure 16.

Variation of the eigen frequencies with bar coupling versus pipe radius m = 1.

Figure 17.

Comparison of eigen frequencies by versus pipe radius for spring coupling and bar coupling (mode m = 1, n = 1).

Results interpretation

When we found the opposite behaviour from the results obtained previously between the natural frequencies of the coupling stiffness and bar share ratio the variation of the length and the radius made it possible to make a study on the eigenfrequencies with the different possible geometric ratio between Radius /Length, Radius/Thickness (R/L, R/e) to have a frequential analytical explanation on this behaviour of the bar coupling, then we carried out a study for which we left all the physical parameters of the validation in such a way facilitates us to conclude the results that we obtain.

This is the reason why we change the ratios from R/L = 1 to 7 and from R/e = 0.002 to 0.14 where we found that the influence of the geometric ratios on the dynamic behaviour of the system of fluid-structure interaction from the natural frequency contains significant because there is a sudden decrease between 0.1 and 1 m then follows the increase of these natural frequencies up to the 7 m length of the pipe, also the frequencies decrease when the ratio between Radius/Length and Radius/Thickness increase as shown in Figures 18 to 20. For different longitudinal modes at constant circumferential mode. This presentation allowed us to conclude that the less frequent part varies between 0.01 and 1 m these or can choose the less frequent pipe for the industrialist.

Figure 18.

Variation of the frequencies of bar coupling versus pipe length with mode m = 1, n = 1.

Figure 19.

Variation of the frequencies of bar coupling versus pipe length with mode m = 2, n = 1.

Figure 20.

Variation of the frequencies of bar coupling versus pipe length with mode m = 9, n = 1.

Case study

Several examples have been considered to study the influence of geometric parameters in the two cases of stiffness and member coupling. Showing us some of the results obtained are shown in Figures 18 to 20.

Conclusion

In this work, a new model of penalty coupling is proposed where the spring element of the interface is replaced by a bar. This model makes it possible to take into account the effect of the radial contact and that of the axial forces. In this study, the model of the case of an internal pressure cylindrical pipe with incompressible laminar flow is applied, of an internal pressure cylindrical pipe with incompressible laminar flow; this was possible by modeling the system using the hierarchical finite element method.

The behavioural equations were determined using the principle of virtual work for the solid, and the Navier-Stokes equations for the fluid. Modeling of the structure-fluid system was done by the hierarchical finite element method using a finite element of the rectangular four-node type with three degrees of freedom per node. A MATLAB programme was developed to calculate the natural frequencies of the pipe carrying the fluid.

The comparison between the results obtained by our programme and the experimental ones indicates that the bar coupling is better than the spring coupling.

In addition, some case studies were carried out and the results obtained are illustrated in a number of graphs which represent the influence of the different geometrical and physical parameters on the eigenfrequencies, in the case of the two types of coupling, i.e. bar coupling and spring coupling, for the fluid-structure interaction phenomenon.

The following conclusions may therefore be drawn:

• - The results obtained using the bar coupling model are very close to those obtained experimentally.

• - Increasing the diameter of the pipe, with constant thickness, increases the stiffness of the structure. Also, as the thickness of the pipe increases, the stiffness of the structure goes down. In the presence of the fluid, the mass of the structure rises.

• - Increasing the pipe length leads to lower structural stiffness, for the different circumferential modes.

• - Finally, increasing the pipe length results in a decrease in the structural stiffness, for the different circumferential modes.