Analytical Solution of Sloshing in a Cylindrical Tank with an Elastic Cover

: Coupled ﬂuid-structure is signiﬁcant in many aspects of engineering applications such as aerospace fuel tanks, the seismic safety of storage tanks and tuned liquid dampers. Numerical investigation of the effects of thin plate cover over a cylindrical rigid fuel tank ﬁlled by an inviscid, irrotational, and incompressible ﬂuid is investigated. Governing equations of ﬂuid motion coupled by plate vibration are solved analytically. A parameter study on the natural frequency of coupled ﬂuid-structure interaction is performed. The results show the non-dimensional natural frequency of coupled ﬂuid-structure is a function of mass ratio, plate elasticity number and aspect ratio. This function is derived numerically for high aspect ratios which in companion with a semi-analytical could be used in the engineering design of liquid tanks with a cover plate.


Introduction
The coupling problem of structure and fluid interaction is a classic problem interested in various engineering applications such as vehicle dynamics, aircraft dynamics, and storage tanks [1][2][3][4][5]. The storage tanks usually built from strong materials such as steels (for drums), stainless steels, concrete, aluminum alloys, fiber-reinforced plastics (for fiber drum), and Timber which can bear large elastic vibrational energy. These tanks are exposed to dead and live loads such as snow, wind (vortex oscillation), seismic, buckling, and earth pressures. When the structure is assumed to have the many mass vibration system, modal analysis is a useful tool for the design of yield shear force of the convective mass vibration, damage to the structure (cumulative plastic strain energy), and water pressure imposed on the tank. That analysis needs the natural frequency of the maximum number of vibration modes that influence seismic responses.
One of the first analytical methods proposed to solve that problem is the Rayleigh-Ritz method [6]. The application of the method of Rayleigh-Ritz on the structure and fluid interaction with weak coupling [5] and strong coupling [7] was investigated by Amabili [7]. Various methods are applied to find the solution of the coupled fluid-structure interaction such as artificial spring method [8], finite element formulation [9], Galerkin method [10], time-marching technique [11], collocation method [12], Rayleigh quotients method [13], boundary elements method [14], boundary integral equation method [15]. Various numerical methods such as finite element, collocation, Rayleigh, boundary elements, and boundary integral method are the approximation to the exact solution.
In the Rayleigh-Ritz method, the approximations to eigenvalue equations are found [5]. The finite element method analytical solution of boundary value problems is used to get a system of algebraic equations. In the Galerkin method [10] which is the popular method of finite element formulation [9], the differential equation is converted from a continuous operator problem to a discrete problem. In the collocation method, the idea is to choose a finite-dimensional space of proposed polynomial solutions and several domain points that satisfy the given equation at the collocation points [12]. In the boundary element method [14] the integral equations are solved by Green's function elements connecting pairs of source and field patches defined by the mesh form a matrix. The boundary element method [15] is often more efficient than other methods [2] when there is a small surface to volume ratio. Conceptually, it works by constructing the surface mesh. However, for many problems boundary element methods are significantly less efficient than volume-discretization methods (finite element method, finite difference method, finite volume method). The added complexity of the boundary integral equation method could be overwhelmed by Compression techniques. Computational time and the storage requirements by boundary element formulations grow according to the square of the problem size while the finite element matrices are typically grown linearly with the problem size.
The fluid-structure coupling problem considered for various fluids such as dense fluids [16], the frictionless liquid in zero gravity [17], and compressible fluid [18], filled with non-viscous liquid [19] and various filling ratios for fluid-filled [19] and partially filled [20] were discussed in the literature. One of the most interesting geometries is a circular cylinder container. A drum (also called barrels in common usage) used in the shipping of dangerous powders and liquids, with 208 liters volume, has an 880 millimeters tall and 610 millimeters diameter (H/a ≈ 2.9). Similar geometry is used for a barrel of crude oil (72 cm × 43 cm).
On the other hand, the types of solid structures connected with fluid were composed of cylindrical baffles [21], annular cylindrical tanks [22], ring-stiffeners and flexible bottom [23], annular baffle [24], vertical baffle [25], submerged components [26], internal bodies [27,28], two identical plates [29], immersed plate in a container [30]. As well the elastic structure bonding the fluids were initially bent cylindrical shells [31] and composite shells with orthogonal stiffeners [32]. Based on the size the containers are flexible and rigid. The design of a fluid-structure coupling problem in a cylindrical tank is an interesting research theme with possible engineering practical applications on stiffeners. Drums have chimes or rims at ends (chines). Most steel drums have reinforcing rolling hoops or rings of thickened metal or plastic. Amabili [23] investigated the effect of ring-stiffeners and flexible bottom on sloshing modes in cylindrical containers and drums.
Amabili [33] addresses a research theme in the area of fluid-solid interaction, focusing the analytical solution of an interesting problem of coupled vibration of a fluid-filled cylindrical rigid tank with an elastic top cover. The two kinds of drums are barrels (welded top) and the open top. Vibrations of a circular flexible membrane or elastic plates resting on a sloshing liquid is a class of problem which first studied by Amabili [33], for annular plate cover by Kim and Lee [34] and fully covered case by Bauer [35]. Bauer [35] makes use of the Fourier-Bessel series expansion to analytically derive the solution of the addressed problem, and perform a parametric study on the natural frequency of the coupled fluid-structure interaction in clamped boundary condition and for varying characteristic values of the system.
Considering all of the above, there is no research performed on the analytical solution of rigid tank coupling with elastic cover with different support conditions. The analytical solution is not easy to use by engineers but is the exact solution of the problem and has the benefits of a minimum number of parameters and accuracy over the other methods such as finite element, collocation, Rayleigh, boundary elements, and boundary integral method. The aim of the current study is finding the analytical solution of a fluid-filled cylindrical container with rigid sides and bottom and elastic circular plate over the top surface. In the next section mathematical model of the system in the fluid and solid parts is presented.

Fluid Model
The schematic of the problem is shown in Figure 1. If the fluid field equations are linearized for small displacements based on the linearized theory of water the flow potential is found from the Laplace equation. By assuming a simple harmonic motion with a radian frequency ω, the velocity potential is In order to calculate the flow potential distribution, φ(r, θ, z, t), in the above equation ( i.e., Equation (1)) based on a typical power series, the fluid domain is defined mathematically as The modal expansion method could be used for the fluid velocity potential as where φ is time-independent velocity potentials, the W ns is time-dependent coefficient of n-th mode, and φ ns is n-th sloshing position dependent eigenfunctions. The governing equation in the fluid region for φ is which satisfies the following boundary conditions: on the symmetric line (for −H < z < 0), on the rigid wall (for −H < z < 0), and on the container bottom (for 0 < r < a). The three-dimensional view of the problem with fluid-solid interaction surface and the side boundary condition is shown in Figure 2. The solution of above equations are (see Equation (1.40 where α mn and β mn are time-dependent coefficients could be found from the initial conditions of free-surface or pressure distribution, λ mn is mode coefficient in fluid domain, J m (λ mn r) is order m of first kind Bessel function. The above definition should satisfy the linearized (neglecting from circumferential and radial displacement) fluid-surface Kinematic boundary condition coupled by solid for 0 < r < a where in the case of fluid without a cover the eigenvalue of the system is found from the linearized boundary condition at free-surface for 0 < r < a in which, g denotes the gravitational acceleration. The use of the analytical solution in the above equation leads to (see Equation (1.43a) in [1]) where ω mn is the frequency of mode mn in fluid domain. As well for the sake of simplicity, the axi-symmetric modes just considered then the solution in the fluid domain is and by assuming the axi-symmetric modes and wall boundary condition at the outer wall (See Equation (6)) then in the fluid domain, the solution is the roots of Equation (13) are presented in Table 1.

Solid Model
The governing equation of the solid motion with the variable of w for the transverse displacement of a bent thin plate is presented by the following partial differential equation where the flexural rigidity (D) is defined by and ∇ 4 is the biharmonic operator that can be expanded as follows where As the pressure could be presented in the form of a harmonic function in time p(t, r, θ, z) = e iωt p(r, θ, z) and displacement could be presented in the form of a harmonic function in time the classical differential equation of solid displacement from Equations (14) can be rewritten as The pressure is made of two components of fluid pressure, p f , and the pressure caused by surface elevation, p z , As the fluid is considered to be stuck to the solid at to boundary, the surface elevation pressure is found by [35] p z = ρ f g w(r, θ, z, t) Using Bernoulli equation for the fluid pressure, p f , and neglecting the nonlinear terms, the pressure of the fluid could be related to the time derivative of velocity potential of the fluid ( φ) and fluid density (ρ f ) [33] Finally the simplified problem of solid motion from Equations (1), (20), (21), (22) and (23) which should be solved is with the free surface condition of

Analytic Solution of Solid Deflection with Fluid Load
The general idea, for the solution of an inhomogeneous linear PDE with Equation (24) is to split its solution into two parts The first term, w h , is the homogeneous solution of Equation (24) and next term, w p , is the in-homogeneous solution of the Equation (24) which has the same boundary condition (plus the initial conditions as the time is a variable) of the full problem.
For the solution of the homogeneous equation, one can assume the Fourier components in θ as, where W nm (r) = a nm J n (k nm r) + b nm Y n (k nm r) + c nm I n (k nm r) + d nm K n (k nm r), and respectively, where Y n and J n are the second and first kind Bessel functions, respectively. As well, the K n , and I n , are second and first kind modified Bessel functions. The mode coefficients (a nm , b nm , c nm , and d nm ) determine the significance of the mode shape and could be found from the initial condition and the boundary conditions. Thus, the general solution of a solid governing equation in polar coordinates is (See Equation (1.18) in [5]) When the polar coordinate origin is considered the same as the circular plate center (without internal hole), the terms of Equation (1.18) in [5] involving Y n (kr) and K n (kr) neglected as the value of the deflection at singular point (r = 0) has a finite value. When these simplifications are employed, the modal displacement will be reduced to the following equation for a typical mode The particular solution of Equation (24) with the aid of Equation (23) and the solution of fluid potential (see Equation (12)) is The solid deflection is cos(nθ) a n J n (kr) + b n I n (kr) + where difference between the above equation (Equation (34)) and Equation (30) is that here the k is a function of ω and could be different from the natural frequency of the plate.

Clamped Plate
The boundary condition of zero deflection at the perimeter is used both for clamped plate and simply supported case. For the simply supported case and clamped plate, the plate displacement is zero on the perimeter of the plate (r = a), By substitution of a instead of r in Equation (33) w(a, θ) = ∞ ∑ n=0 cos(nθ) a n J n (ka) which could be expanded to Since the boundary condition of Equation (35) for constant terms (i.e., n = 0) leads to and for each cosine term (i.e., cos(nθ) where n = 0) leads to a n J n (ka) The second boundary condition which is devoted to just clamped plate case is The first derivative of deflection with respect to the r at r = a is cos(nθ) a n kJ n (ka) + b n kI n (ka)+ By considering the wall boundary condition for the fluid (See Equation (13)) the boundary condition of Equation (41) for constant terms ( i.e., n = 0 ) leads to and for each cosine term (i.e., cos(nθ) where n = 0) leads to (i.e., m = 0 and n = 1, ..., ∞) a n J n (ka) + b n I n (ka) = 0.

Simply Supported Plate
Simply supported plate is a combination of zero deflection at perimeter (i.e., Equations (38) and (39)) and zero radial twisting moment at perimeter. As well the radial twisting moment where at plate edge must be zero. i.e., The first deriviate of deflection respect to the r at r = a is retrieved from Equation (41) and the second deriviate of deflection respect to the r at r = a is The second derivative of deflection with respect to θ at r = a is cos(nθ) a n J n (ka) + b n I n (ka) + Since cos(nθ) a n k 2 J n (ka) + b n k 2 I n (ka)+ The boundary condition of Equation (45) for constant terms (i.e., n = 0) leads to or in a simplified way and for each cosine term (i.e., cos(nθ) where n = 0) leads to a n k 2 J n (ka) + b n k 2 I n (ka)+ rearranged in the form a n k 2 J n (ka) + νk a J n (ka) − n 2 ν a 2 J n (ka) + b n k 2 I n (ka) + νk a I n (ka) − n 2 ν a 2 I n (ka) + mn λ 2 mn J n (λ mn a) + νλ mn a J n (λ mn a) − n 2 ν a 2 J n (λ mn a) = 0,

Kinematic Boundary Condition at Free Surface
Finally by considering the Kinematic boundary condition (See Equation (9)) of the fluid-solid interaction for the region of the top surface of the fluid under the circular plate with the time dependency relations (See Equations (1) and (19)) for all values of r, θ, t on the plate position which is simplified by the aid of Equations (1) and (19) to The equation (Equation (54)) can be explained by the system of equations with fluid motion modes expansion.
cos(nθ) a n J n (k nm r) + b n I n (k nm r) + The boundary condition at top surface (Equation (54)) a n cos(nθ) a n J n (kr) + b n I n (kr) + The boundary condition of Equation (54) for constant terms (i.e., n = 0) leads to and for each cosine term (i.e., cos(nθ) where n = 0) leads to a n J n (kr) + b n I n (kr)+ mn − tanh(λ mn (H))λ mn J n (λ mn r) = 0, Multiplying the above equation (Equation (56)) by J m (λ mn r)r and integrating over the whole plate domain (i.e., × a 0 J m (λ mn r)rdr) for m = 0 to m = ∞ the system of equations are closed.
If the variables seen in a matrix form The equations can be re-written in the form matrix as  Table 3. Matrix of equations and variables for simply supported case.

Equation Number
Number of Equations α 00 α m0 a 0 b 0 a n b n α mn Table 4. Matrix of equations and variables for clamped case.

Equation Number
Number of Equations α 00 α m0 a 0 b 0 a n b n α mn

Simply Supported Plate
In simply supported case, the rows of A are derived from Equations (38), (48), and (56) respectively.
The coefficient matrix would be obtained as where • A 11 is a 1 × 1 matrix as • A 12 is a 1 × M matrix as where m is from 1 to M. • A 23 is a 1 × 1 matrix as • A 24 is a 1 × 1 matrix as were l is from 0 to M and in calculation of A 31 (l + 1, 1) instead of J 0 (λ 00 r) the uniform unit function ( f (r) = 1) should be considered.i.e., where k is from 0 to M and m is from 1 to M. As well in calculation of A 32 (k + 1, m) at k =0 instead of J 0 (λ 00 r) the uniform unit function ( f (r) = 1) should be considered.
where k is from 0 to M. In calculation of A 33 (k + 1, 1) at k =0 instead of J 0 (λ 00 r) the uniform unit function ( f (r) = 1) should be considered.
• A 45 is a diagonal N × N matrix where the element of A 45 (n, n) (n is from 1 to N) where defined as • A 46 is a diagonal N × N matrix where the element of A 46 (n, n) (n is from 1 to N) which defined as • A 47 is a N × M matrix as A 47 = ρ f ω 2 λ 2 mn J n (λ mn a) + νλ mn a J n (λ mn a) − n 2 ν a 2 J n (λ mn a) n is from 1 to N and m is from 1 to M. • A 55 is a diagonal N × N matrix where the element of A 55 (n, n) (n is from 1 to N) where defined as • A 56 is a diagonal N × N matrix where the element of A 56 (n, n) (n is from 1 to N) which defined as • A 57 is a diagonal N × N matrix where the element of A 57 (n, n) (n is from 1 to N) wheredefined as

Clamped Plate
The matrix A is defined by The coefficient matrix would be obtained as The coefficients of A in above equation (Equation (86)) are same as simply supported case except: • A 23 is a 1 × 1 matrix as • A 24 is a 1 × 1 matrix as • A 55 is a diagonal N × N matrix where the element of A 55 (n, n) (n is from 1 to N) are defined as • A 56 is a diagonal N × N matrix where the element of A 56 (n, n) (n is from 1 to N) which defined as

Results and Discussion
The analytical solutions which are obtained in the previous section, are investigated numerically for the practical case. The numerical results in the current section have been derived by developing code inside the MATLAB software. The natural frequencies are presented in Hz unit where f = ω 2π . Usual material properties of the water tank storage are given in Table 5. The values are not fixed and are averaged from various references. For example, the stone Poisson ratio is reported between 0.2 and 0.3 in various references but the value of 0.25 is chosen for the table. The last two rows of the Table 5 calculated for the special case of h = 1 mm and a = 1 m. As shown for that values the first parameter  Table 6. Relative errors of computing axi-symmetric modes of the clamped cover plates are presented in Table 7. Table 7 presents the relative error ( f f exact − 1) versus mode number and aspect ratios where the solution with ten expansion terms is considered as the exact frequency (M = 10). The performance of the method is comparable with the Table 4 in Amabili [33] where the convergence of the solution was discussed. As shown the use of the current method causes the precise results in most cases by three terms in fluid expansion series with three terms. Table 4 in Amabili [33] used 10 terms for the solid motion (in Table 7 just one term is used for solid deflection) and 50 terms for liquid potential (As seen in Table 7 Tables 1 and 2 of Amabili [33] are presented in Table 8. As shown a good agreement seen among the current study results and those of Bauer ([35]) and Amabili ([33]) for axi-symmetric modes of the clamped cover plates (for n = 0 and n = 1 modes). Results were found by one-term expansion of the solid plate deflection and 4 terms for the liquid potential.     The non-dimensional fluid-structure frequencies for simply supported boundary condition (presented by the dashed line in Figure 3) and clamped boundary condition (presented by solid line in Figure 3) at ρ p h p ρ f a = 10 −3 , ρ f ga 4 D = 500, and axi-symmetric mode n = 0 is plotted in Figure 3. As shown, the agreement of results of the current study with non-dimensional fluid-structure frequencies for clamped plate obtained by Bauer ([35]) is perfect. As exposed in all cases the clamped boundary condition caused higher frequencies than the simply supported plate.

Clamped Plate
Assessment of the present results and the Bauer ( [35]) results for other figures in plate cover case has been performed in Figure 4. The dominant modes of the structure i.e., the first three mode shapes are illustrated in Figure 4. The non-dimensional fluid-structure frequencies for clamped plate = 500, and axi-symmetric mode n = 2 is plotted in Figure 4. As shown, the agreement of results of the current study with non-dimensional fluid-structure frequencies for clamped plate obtained by Bauer ([35]) are matching. As for the common storage tank the aspect ratio ( h a ) is above the unity the parameter study of the natural frequencies as the function of the elastic parameter (

Simply Supported Plate
The modal analysis is required to find the inherent dynamic properties of the any domain in terms of its natural frequencies. The non-dimensional fluid-structure frequencies for simply supported plate = 500, and axi-symmetric mode n = 2 is plotted in Figure 8. The dominant modes of the structure i.e., the first three mode shapes are illustrated in Figure 8. As illustrated in Figure 8, the natural frequencies are increased by an increase of aspect ratio. The comparison of the Figure 8 and Figure 4 is not simple such as the results of zero modes Figure 3 for clamped plate and simply supported case. As exposed in some cases the clamped boundary condition caused higher frequencies than the simply supported plate and in some cases it was vice versa. The non-dimensional fluid-structure frequencies for clamped plate versus ρ f ga 4 D and axi-symmetric mode n = 0 for various ρ p h p ρ f a is plotted in Figure 9. As shown the parameter study of the natural frequencies as the function of the elastic parameter ( ρ f ga 4 D ) and mass ratio ( ρ p h p ρ f a ) presents a linear relation between the parameters in comparison with the nonlinear case was shown for clamped case in Figure 5. As shown by the increase of elastic parameter ( ρ f ga 4 D ) and mass ratio ( ρ p h p ρ f a ) the non-dimensional circular frequency (ω * = ωa 2 µ D ) is increased. The sensitivity analysis of the results show that even the mass ratio ( ρ p h p ρ f a ) is the most important parameter on the detect of natural frequency, the elastic parameter ( ρ f ga 4 D ) is more important than the nonlinear case was shown for clamped case in Figure 5. The non-dimensional fluid-structure frequencies for clamped plate versus ρ f ga 4 D and axi-symmetric mode (n = 1) for various ρ p h p ρ f a is plotted in Figure 10. In this case in contract with linear behavior of the clamped case in Figure 5 a parabola shape is seen. It means that there are some optimal values of elastic parameter ( ρ f ga 4 D = 500) where the natural angular velocity of the axi-symmetric mode (n = 1) is minimized. This frequency is important in horizontal movement of the sloshing tanks. The non-dimensional fluid-structure frequencies for clamped plate versus ρ f ga 4 D and axi-symmetric mode n = 2 for various ρ p h p ρ f a is plotted in Figure 11. The fundamental angular velocity of the axi-symmetric mode (n = 2) shows a parabola shape same as the fundamental angular velocity of the axi-symmetric mode (n = 1).

Conclusions
The analytical investigation of the effects of thin plate cover over a cylindrical rigid fuel tank filled by an inviscid, irrotational, and incompressible fluid was investigated in the current study. Governing equations of fluid motion coupled by plate vibration are found by the application of the Fourier-Bessel series expansion. A parameter study on the fundamental angular velocity of coupled fluid-structure motion is performed. The results show the non-dimensional fundamental angular velocity of a coupled fluid-structure system is a function of mass ratio, plate elasticity number and aspect ratio. This function is derived numerically for high aspect ratios which in companion with a semi-analytical could be used in the engineering design of liquid tanks with a cover plate. As illustrated by an increase of the elastic parameter (