Dynamic Stability of Plane Free Surface of Liquid in Axisymmetric Tanks

Abstract

When liquid filled containers are excited vertically, it is known that, for some combinations of frequency and amplitude, the free surface undergoes unbounded motion leading to instability, called parametric instability or parametric resonance, while for other combinations the free surface remains plane. In this paper, the stability of the plane free surface is investigated theoretically when the vessel is a vertical axisymmetric container. The effect of coupled horizontal excitation on the stability is examined. The dynamics of sloshing flows under specified excitations are simulated numerically using fully nonlinear finite element method based on non-linear potential flow theory. A mixed Eulerian-Lagrangian technique combined with 4th-order Runge-Kutta method is employed to advance the solution in time. A regridding technique based on cubic spline is applied to the free surface for every finite time step to avoid possible numerical instabilities.

1. Introduction

The motion of the unrestrained free surface of the liquid due to external excitation in the liquid filled containers is known as sloshing. Sloshing is likely to be seen whenever we have a liquid with a free surface in the presence of gravity. At equilibrium the free surface of the liquid is static, when the container is perturbed; an oscillation is set up in the free surface. The phenomenon of liquid sloshing occurs in a variety of engineering applications such as sloshing in liquid propellant launch vehicles, liquid oscillation in large storage tanks by earthquake, sloshing in the nuclear reactors of pool type, nuclear fuel storage tanks under earthquake, and the water flow on the deck of ship. Such liquid motion is potentially dangerous problem to engineering structures and environment leading to failure of engineering structures and unexpected instability. Thus, understanding the dynamic behaviour of liquid free surface is essential. As a result the problem of sloshing has attracted many researchers and engineers motivating to understand the complex behaviour of sloshing and to design the structures to withstand its effects.

Liquid sloshing can be stimulated by a variety of container excitations. The container excitation can be horizontal, vertical, or rotational. Under horizontal excitations the liquid free surface experiences normal sloshing; the sloshing frequency will be equal to excitation frequency. When the external excitation frequency is equal to fundamental slosh frequency, the free surface undergoes resonance. Extensive research has been done on sloshing response under pure horizontal excitations. When the liquid filled container is subjected to vertical excitations, for some combinations of amplitude and frequency of the external excitation, the free surface undergoes unbounded motion leading to instability called parametric instability, and for few other combinations the free surface undergoes bounded motion. In case of parametric instability the rate of growth of amplitude is exponential, and this exponentially growing response is potentially dangerous to the system.

The problem of liquid response under vertical excitations was first studied experimentally by Faraday [1], reporting that the frequency of the liquid vibrations on free surface is half of the external excitation frequency. The sloshing waves generated under vertical excitation are sometimes referred to as Faraday waves. Rayleigh [2, 3] analyzed Faraday’s study and confirmed Faraday’s observations. Matthiessen [4] conducted experiments and reported that the fluid free surface vibrations are synchronous to the external excitation. The discrepancy between Faraday’s observations and Matthiessen’s observations was explained mathematically by Benjamin and Ursell [5]. Benjamin and Ursell [5] investigated the problem theoretically. They considered linearized inviscid potential flow model with surface tension. They concluded that the response of the plane free surface of fluid under vertical excitation is governed by Mathieu equation. The solution of Mathieu equation [6] may be stable, periodic or unstable depending on the system parameters. The stability and instability of the Mathieu equation are shown in the form of plots of amplitude versus frequency of external excitation which gives regions of stability and instability. Benjamin and Ursell concluded that, if the response of the free-surface is unstable, the resulting motion can have a frequency equal to (1/2)𝑛Ω, where 𝑛 is an integer and Ω is the natural sloshing frequency. The problem of sloshing under vertical excitation in various geometric shapes of the container has been studied by various researchers. For example, Dodge et al. [7] and Miles [8] have studied the liquid surface oscillations under vertical excitation in cylindrical tanks. Khandelwal and Nigam [9] have studied the parametric instability in rectangular tanks. The problem of parametric oscillations was discussed by Miles and Henderson [10], in their review paper. Most of the studies available in the old literature on the problem of parametric sloshing under vertical excitation are both experimental or analytical involving complex equations and heavy derivations. In order to understand the complex behaviour of sloshing under vertical excitations including nonlinear terms, numerical simulations are advantageous compared to analytical solutions. Analytical solutions get complicated when the shape of container is not regular. Numerical methods like finite element method, finite volume method, finite difference method, boundary element method, and so forth are available for numerical modeling of sloshing waves. The studies on parametric sloshing using numerical modeling are very few. Frandsen [11] analysed the problem numerically and theoretically considering fully nonlinear inviscid potential flow equations. Frandsen applied finite difference method for sloshing response. Frandsen study was on 2-dimensional rectangular tanks. Wu et al. [12] applied finite element method for solving sloshing 2D and 3D sloshing problems. Wu et al. discussed the sloshing response under vertical excitations in rectangular tanks. Ning et al. [13] applied boundary element method to study the liquid sloshing in rectangular containers under coupled horizontal and vertical excitations. In the present paper the sloshing response under vertical excitations in axisymmetric tanks, that is, upright cylindrical and annular cylindrical tanks, is taken up. First, the stability of plane free surface of liquid in axisymmetric tanks is obtained from the linearized equations and the sloshing response of fluid is numerically simulated for various frequencies and amplitude of the external excitation using finite element method under pure vertical excitations. Then the numerical work is extended to explore the sloshing behaviour in axisymmetric tanks under coupled horizontal and vertical excitations. According to the authors’ knowledge there is no study available in the literature on parametric sloshing using finite element method in axisymmetric tanks.

2. Governing Equations

Consider that an axisymmetric tank fixed in cylindrical coordinate system 𝑂𝑟𝜃𝑧 is moving with respect to inertial cylindrical coordinate system 𝑂₀𝑟₀𝜃₀𝑧₀. The origins of these systems are at center of the tank at the free surface and pointing upwards in 𝑧-direction. These two coordinate systems coincide when the tank is at rest. Figure 1 shows the axisymmetric tanks in the moving cylindrical coordinate system 𝑂𝑟𝜃𝑧 along with the prescribed boundary conditions.

(a)

(b)

Figure 1 

Axisymmetric tanks and their boundary conditions.

Let the tank be displaced horizontally along 𝑟 by 𝑟_𝑡(𝑡) and vertically along 𝑧 by 𝑧_𝑡(𝑡). The displacements of the tank are governed by 𝑋=𝑟_𝑡(𝑡),𝑧_𝑡(𝑡). (1) Let the tank be filled with fluid to a depth ℎ. Fluid is assumed to be inviscid and incompressible, and motion is irrotational. Therefore the fluid motion is governed by Laplace’s equation with the unknown as velocity potential 𝜙: ∇²∅=0,where∇²=𝜕²𝜕𝑟²+1𝑟𝜕+1𝜕𝑟𝑟²𝜕²𝜕𝜃²+𝜕²𝜕𝑧². (2)(3) Equation (2) is an elliptic partial differential equation and needs Neumann and Dirichlet boundary conditions to solve. The fluid obeys Neumann boundary conditions on the walls of the container and Dirichlet boundary conditions on the liquid free surface. In the moving coordinate system the velocity component of the fluid normal to the walls is zero. Hence, on the bottom and on the walls of the tank (Γ_𝐵) we have 𝜕∅𝜕𝑟_{𝑟=𝑅,𝑅in}=0;𝜕∅𝜕𝑧_𝑧=−ℎ=0. (4) On the free surface (Γ_𝑆), liquid obeys dynamic and kinematic boundary conditions, respectively, given by 𝜕∅𝜕𝑡_𝑧=𝜉+12∇∅⋅∇∅+𝑔+𝑧_𝑡^′′𝜉+𝑟𝑟_𝑡^′′=0,𝜕𝜉+𝜕𝑡𝜕∅𝜕𝑟𝜕𝜉−𝜕𝑟𝜕∅𝜕𝑧=0, (5)(6) where 𝜉 is the free surface elevation measured vertically above the still water level, 𝑟_𝑡^′′ and 𝑧_𝑡^′′ are the horizontal and vertical accelerations of the tank, and 𝑔 is acceleration due to gravity.

3. Governing Equation for Dynamic Stability of Free Surface

The possible general solution of the Laplace equation, which satisfies the previously given boundary conditions on the rigid surfaces for axisymmetric tanks, can be represented as [14] ∅=^∞∞_𝑚=0𝑛=1cosh𝜌_𝑚𝑛/𝑅(𝑧+ℎ)cosh𝜌_𝑚𝑛/𝑅ℎcos𝑚𝜃𝑀_𝑚𝜌_𝑚𝑛𝑅𝑟𝐹_𝑚𝑛(𝑡),𝜉=^∞∞_𝑚=0𝑛=1cos𝑚𝜃𝑀_𝑚𝜌_𝑚𝑛𝑅𝑟𝑧_𝑚𝑛(𝑡), (7) where 𝑚 is the circumferential mode number and 𝑛 is the radial wave number. 𝐹_𝑚𝑛(𝑡),𝑧_𝑚𝑛(𝑡) are the time evolution functions and will be determined from the initial boundary conditions. The function 𝑀_𝑚(𝜌_𝑚𝑛𝑟) is different for cylindrical and annular tanks and is given by [14] 𝑀_𝑚𝜌_𝑚𝑛𝑅𝑟=𝐽_𝑚𝜌_𝑚𝑛𝑅𝑟𝐽forcylindricaltanks,_𝑚𝜌_𝑚𝑛𝑅𝑟𝑌^′_𝑚𝜌_𝑚𝑛𝑅−𝐽^′_𝑚𝜌_𝑚𝑛𝑅𝑌_𝑚𝜌_𝑚𝑛𝑅𝑟forannulartanks, (8) where in (9), 𝐽_𝑚 is the Bessel function of the first kind of order 𝑚 and 𝑌_𝑚 is the Bessel function of the second kind of order 𝑚. For cylindrical tank, 𝜌_𝑚𝑛 are the roots of 𝜕𝐽_𝑚𝜌_𝑚𝑛/𝑅𝑟𝜕𝑟_𝑟=𝑅=0. (9) And for annular cylindrical tank, 𝜌_𝑚𝑛 are the roots of 𝐽^′_𝑚𝜌_𝑚𝑛𝑌^′_𝑚𝑘𝜌_𝑚𝑛−𝑌^′_𝑚𝜌_𝑚𝑛𝐽^′_𝑚𝑘𝜌_𝑚𝑛=0, (10) where 𝑘=𝑅_in/𝑅 is the ratio of inner radius to outer radius of the tank.

To obtain the equation which governs the free surface stability, the free surface boundary conditions (5) and (6) are linearized and the solution of Laplace equation is substituted in the linearized boundary conditions. The linearized free surface boundary conditions obtained considering pure vertical excitation exists are 𝜕∅𝜕𝑡_𝑧=𝜉+𝑔+𝑧_𝑡^′′𝜉=0,𝜕𝜉−𝜕𝑡𝜕∅𝜕𝑧=0. (11) Substituting (7) into (11), we get 𝑑𝐹_𝑚𝑛(𝑡)+𝑑𝑡𝑔+𝑧_𝑡^′′𝑧_𝑚𝑛(𝑡)=0,𝑑𝑧_𝑚𝑛(𝑡)𝑑𝑡−𝜆_𝑚𝑛tanh𝜆_𝑚𝑛ℎ𝐹_𝑚𝑛(𝑡)=0. (12)(13) Differentiating (12) w.r.t. time and substituting (13), we get 𝑑²𝑧_𝑚𝑛(𝑡)𝑑𝑡²+𝜔²_𝑚𝑛𝑧1+^𝑡′′𝑔𝑧_𝑚𝑛(𝑡)=0, (14) where 𝜔_𝑚𝑛=𝑔𝜌_𝑚𝑛𝑅𝜌tanh^𝑚𝑛ℎ𝑅. (15) Equation (15) gives the slosh frequencies of the free surface for cylindrical and annular tanks. The values of 𝜌_𝑚𝑛 are different for cylindrical and annular tanks and are calculated from (9), (10), respectively. The tank is assumed to be excited periodically with 𝑧_𝑡(𝑡)=𝑎_𝑣cos(𝜔_𝑣𝑡), where 𝑎_𝑣 is vertical forcing amplitude, 𝜔_𝑣 is vertical forcing angular frequency, and 𝑡 is time. Under harmonic vertical excitation (14) reduces to 𝑑²𝑧_𝑚𝑛(𝑇)𝑑𝑇²+Ω²_𝑚𝑛1−𝑘_𝑣cos(2𝑇)𝑧_𝑚𝑛(𝑇)=0, (16) where 𝑇=(1/2)𝜔_𝑣𝑡, Ω_𝑚𝑛=𝜔_𝑚𝑛/𝜔_𝑣, and 𝑘_𝑣=𝑎_𝑣𝜔²_𝑣/𝑔. Equation (16) is Mathieu’s equation. The stability and instability of the free surface under vertical excitation are guided by (16).

Mathieu equation is a second-order differential equation with periodic coefficients. The solution of Mathieu equation may be bounded or unbounded, that is, stable or unstable. If the solution is bounded, it may be periodic, or nonperiodic. The theory of Mathieu equation is well documented by McLachlan [6], and Mathieu equation applications concerning parametric instability are documented by Bolotin [15]. The form of the solution for Mathieu equation can be obtained by Floquet’s theory [16]. According to Floquet theory the solution of Mathieu equation can be expressed as a linear combination of two linearly independent Floquet solutions 𝐹_𝑟(𝑧) and 𝐹_𝑟(−𝑧). The Floquet solution can be represented in the form 𝐹_𝑟(𝑧)=𝑒^𝑖𝑟𝑧𝑝(𝑧), (17) where 𝑟 is called characteristic exponent and depends on Ω_𝑚𝑛,𝑘_𝑣 and  𝑝(𝑧) is a periodic function with period 𝜋. The behaviour of the solution can be obtained from (17). Solution of Mathieu equation is bounded when the value of 𝑟 is real and unbounded when the value of 𝑟 is complex. In stable solution, if 𝑟 is irrational then solution is not periodic, if 𝑟 is rational then solution is periodic but not with period 𝜋 or 2𝜋, and if 𝑟 is an integer then solution is periodic with period 𝜋 or 2𝜋.

The behaviour of the solution can be obtained from the stability chart, a plot of system parameters Ω_𝑚𝑛, 𝑘_𝑣. The plots have regions of stability and instability, from which behaviour of solution can be sought. The stable and instable regions are separated by boundary curves, on which the solution is periodic with period 𝜋 or 2𝜋. To plot the stability chart it is needed to plot only the boundary curves. The periodic solutions on the boundary curves can be expanded as a Fourier series [15]. The periodic solution with a period 2𝜋 can be written in the form 𝑧(𝑡)=^∞_𝑘=1,3,5𝑎_𝑘sin𝑘𝑡2+𝑏_𝑘cos𝑘𝑡2. (18) Substituting the series (18) into (16) and equating the coefficients of identical sine and cosine terms lead to the following system of linear homogenous algebraic equations: 𝑘1+^𝑣2−14Ω²_𝑚𝑛𝑎₁−𝑘_𝑣2𝑎₃𝑘=0,1−²4Ω²_𝑚𝑛𝑎_𝑘−𝑘_𝑣2𝑎_𝑘−2+𝑎_𝑘+2𝑘=0𝑘=3,5,7,…,1−^𝑣2−14Ω²_𝑚𝑛𝑏₁−𝑘_𝑣2𝑏₃𝑘=0,1−²4Ω²_𝑚𝑛𝑏_𝑘−𝑘_𝑣2𝑏_𝑘−2+𝑏_𝑘+2=0𝑘=3,5,7,…. (19) The periodic solution with period 𝜋 can be expressed in Fourier series as 𝑧(𝑡)=𝑏₀+^∞_𝑘=2,4,6𝑎_𝑘sin𝑘𝑡2+𝑏_𝑘cos𝑘𝑡2. (20) Substituting the series (20) into (16) and equating the coefficients of identical sine and cosine terms lead to the following system of linear homogenous algebraic equations: 11−Ω²_𝑚𝑛𝑎₂−𝑘_𝑣2𝑎₄𝑘=0,1−²4Ω²_𝑚𝑛𝑎_𝑘−𝑘_𝑣2𝑎_𝑘−2+𝑎_𝑘+2𝑏=0𝑘=4,6,8,…,₀−𝑘_𝑣2𝑏₂1=0,1−Ω²_𝑚𝑛𝑏₂−𝑘_𝑣2𝑏₀+𝑏₄𝑘=0,1−²4Ω²_𝑚𝑛𝑏_𝑘−𝑘_𝑣2𝑏_𝑘−2+𝑏_𝑘+2=0𝑘=4,6,8,…. (21) The system of linear homogenous equations (19) and (21) has a nontrivial solution when the determinant composed of the coefficients is zero. The determinants are written as 𝑘1±^𝑣2−14Ω²_𝑛−𝑘_𝑣2−𝑘0⋯_𝑣291−4Ω²_𝑛−𝑘_𝑣2⋯𝑘0−^𝑣21−254Ω²_𝑛⋯⋮⋮⋮⋱=0. (22) Equation (22) gives determinant obtained from both the conditions of (19) combined under the ± sign: 11−Ω²_𝑛−𝑘_𝑣2−𝑘0⋯_𝑣241−Ω²_𝑛−𝑘_𝑣2⋯𝑘0−^𝑣21−16Ω²_𝑛⋯𝑘⋮⋮⋮⋱=0,1−^𝑣200⋯−𝑘_𝑣11−Ω²_𝑛−𝑘_𝑣2𝑘0⋯0−^𝑣241−Ω²_𝑛−𝑘_𝑣2⋯𝑘00−^𝑣21−16Ω²_𝑛⋯⋮⋮⋮⋮⋱=0. (23) By solving the previously determinants, a stability map can be plotted. The stability chart for (16) is shown in Figure 2. From Figure 2, we can predict the stability of the free surface. If the amplitude and frequency of the external excitation lie inside the curves (grey colored region), the free surface is unstable, and the amplitude grows unbounded exponentially. If the parameters of external excitation are outside the curves the free surface is stable. Thus from the stability charts, the stability of the plane free surface of fluid can be predicted.

Figure 2 

Stability chart for dynamic stability of free surface under vertical excitations.

4. Numerical Modeling of Nonlinear Sloshing

In the present section, numerical modeling of nonlinear sloshing response using finite element method is discussed. The external excitation will excite only the antisymmetric modes of sloshing in an axisymmetric tank; one can expand the velocity potential 𝜙 in Fourier series of cos𝑚𝜃-type modes. Thus the velocity potential takes the form 𝜙(𝑟,𝑧,𝜃,𝑡)=𝜙(𝑟,𝑧,𝑡)cos𝑚𝜃 [17]. By substituting this into (2), we get the Laplace equation and the Neumann boundary condition for axisymmetric tanks as ∇²∅=0,where∇²=𝜕²𝜕𝑟²+1𝑟𝜕+𝜕𝜕𝑟𝜕𝑧²²−𝑚²𝑟²,𝜕∅𝜕𝑟_{𝑟=𝑅,𝑅in}=0;𝜕∅𝜕𝑧_𝑧=−ℎ=0, (24)(25)(26) where the angular variable 𝜃 is separated. The governing equation is independent of the variable 𝜃 now; similarly the boundary conditions become independent of the variable 𝜃. Thus the 3-dimensional problem of sloshing in axisymmetric tanks can be modeled numerically as a 2-dimensional problem in the plane 𝑂𝑟𝑧. In the present study, non-linear sloshing analysis is carried out completely using 2-dimensional finite element method. The finite element model based on mixed Eulerian-Lagrangian scheme is adopted. The free surface nodes behave like Lagrangian particles and interior nodes behave like Eulerian particles. For this formulation, the free-surface kinematic and boundary conditions (5) and (6), respectively, are modified and written in Lagrangian form as [18] 𝜕∅𝜕𝑡_𝑧=𝜉=12∇∅⋅∇∅−𝑔+𝑧_𝑡^′′𝜉−𝑟𝑟_𝑡^′′=0,𝑑𝑟=𝑑𝑡𝜕∅;𝜕𝑟𝑑𝑧=𝑑𝑡𝜕∅.𝜕𝑧 (27)(28) The problem of sloshing is non-linear because the free-surface position is not known a priori and the boundary conditions have non-linear terms. The sloshing problem is evaluated as an initial boundary value problem; to solve the problem, fluid is assumed to be at rest with some initial perturbation on the free surface. The initial conditions for the free surface in the moving Cartesian system at 𝑡=0 and 𝑧=0 are assumed as ∅(𝑟,0,0)=−𝑟𝑑𝑟_𝑡(𝑡)𝑑𝑡−𝑧𝑑𝑧_𝑡(𝑡),𝜉()𝜉()𝑑𝑡𝑟,0=0forhorizontalexcitation,𝑟,0=𝜉₀forverticalexcitation, (29)(30)(31) where 𝜉₀ is the initial elevation of the free surface. Initial condition in (31) is used in case of pure vertical excitations only because in case of vertical excitation some initial perturbation is needed on free surface and without this initial perturbation, there will not be any oscillation in the fluid free surface. Equations (24)–(31) give complete behaviour of nonlinear sloshing flow under coupled horizontal and vertical excitations. In order to solve this nonlinear sloshing problem, time interval 𝑡 is divided into a finite number of time steps, 𝑡_𝑛=𝑛Δ𝑡 (𝑛=0,1,2,3,…), at a particular time step (𝑛=0), the initial boundary conditions (29)–(31) are known, and using these initial conditions along with the boundary condition (26), Laplace equation (24) is solved to get velocity potential 𝜙, with which velocity 𝑣 is evaluated; with these evaluated velocities the kinematic and dynamic free surface boundary conditions, (27)-(28) are time integrated and the position of free surface is updated to get the free surface position for the next time step (𝑛=1). In this manner the sloshing response is numerically simulated. The authors have developed a finite element numerical formulation for nonlinear sloshing response for sloshing in 2D rectangular tanks [19]. The same numerical formulation is followed here; the present simulation is an extension of [19] to axisymmetric containers. In the present paper the numerical formulation is explained briefly; for complete in-depth explanation and algorithm readers can refer to [19].

4.1. Finite Element Formulation

Four noded axisymmetric isoparametric quadrilateral ring elements are used to discretize the fluid domain. By introducing the finite element shape functions the liquid velocity potential can be approximated as 𝜙(𝑥,𝑧)=^𝑛_𝑗=1𝑁_𝑗(𝑥,𝑧)𝜙_𝑗, (32) where 𝑁_𝑗 is finite element shape function, 𝑛 is the number of nodes in the element, and 𝜙_𝑗 is the nodal velocity potential. On applying Galerkin residual method to Laplace equation, we get 𝐾𝜙=0 (33) with matrix 𝐾 defined by 𝐾=_Ω(∇𝑁)^𝑇(∇𝑁)𝑑Ω, (34) where Ω is complete liquid domain. Equation (33) is used to calculate velocity potential 𝜙, for the interior nodes using the known free surface velocity potential.

4.2. Velocity Recovery

To get the free surface position for the next time step it is required to compute velocities. In order to derive a smoothed, accurate, and continuous velocity, patch recovery technique [20] is applied. In patch recovery technique, the continuous velocity field is obtained by considering the linear interpolation of the velocities at the Gauss integration points, 𝑣̂=𝑎₁+𝑎₂𝜉+𝑎₃𝜂+𝑎₄𝜉𝜂, (35) where 𝑣̂ is any velocity component (𝑣_𝑥̂ or 𝑣_𝑦̂), 𝜉,𝜂 are the Gauss locations, and 𝑎₁, 𝑎₂, 𝑎₃, 𝑎₄ are unknowns which need to be evaluated. To evaluate these unknowns, a least square fit is considered between 𝑣̂ and 𝑣, 𝐹(𝑎)=⁴_𝑖=1𝑣̂𝜉_𝑖,𝜂_𝑖−𝑣𝜉_𝑖,𝜂_𝑖², (36) where 𝑖 is 2×2 order Gauss integration points. Then, the four unknown coefficients are determined from four simultaneous equations obtained from 𝜕𝐹(𝑎)𝜕𝑎_𝑘=0;𝑘=1,2,3,4. (37) Substituting the obtained 𝑎_𝑘’s in (35) gives the velocity values for individual elements and these are averaged for the common nodes. Finally, a smoothed velocity field which is interelement continuous is constructed by interpolating the finite element shape functions used in (32) and nodal averaged velocities. The global continuous velocity field 𝑣(𝑣_𝑥or𝑣_𝑦) is given as 𝑣=𝑁⋅𝑣.̂ (38)

4.3. Numerical Time Integration and Free Surface Updating

After calculating the velocity at a time step 𝑡, we need to calculate the position of free surface from (28) and determine the potential on the free surface using (27) for the next time step 𝑡+Δ𝑡. As a result, the liquid mesh and the boundary condition required for the next time step are established. The fourth-order Runge-Kutta method using explicit time integration is employed in the present paper for time marching [19].

4.4. Regridding Algorithm

At the beginning of the numerical simulation, the free surface nodes are uniformly distributed along the free surface. As the time proceeds the free surface nodes are spaced unequally and cluster into a steep gradient leading to numerical instability. This problem occurs for a long time simulation; to avoid this instability an automatic regridding condition using cubic spline is employed when the movement of the nodes is 75% more or less than the initial grid spacing [19].

5. Numerical Results and Discussion

A code is developed following the previous numerical formulation for computing sloshing response under vertical excitations and coupled excitation in axisymmetric tanks. A rigid cylindrical tank of radius 1.115 m filled with water to a depth of half of radius and a rigid annular tank of outer radius 1.115 m, inner radius of 0.2230 m (0.2 times outer radius, 𝑘=𝑅_in/𝑅=0.2) filled with water to a depth of half of net radius, are considered for the numerical simulation of liquid sloshing response. The finite element mesh for the axisymmetric tanks considered is shown in Figure 3. For different values of 𝑚=0,1,2… different circumferential modes are obtained. The dominant mode in most of the applications is the first mode, 𝑚=1, because this mode exerts most fluid forces onto the tank wall [17]. In the present numerical simulation 𝑚=1 is considered. From now onwards, the suffix 𝑚 from all the equations is dropped and it should be noted that 𝑚=1.

(a)

(b)

Figure 3 

A schematic mesh for axisymmetric tanks.

5.1. Free Vibration Analysis

A free vibration problem is solved first to validate the code for stiffness matrix formulation. A mass matrix 𝑀 (as given in (39)) for the fluid free surface is computed: 1𝑀=𝑔_Γ𝑠𝑁^𝑇𝑁𝑑Γ_𝑠. (39) If 𝜔_𝑛(𝑚=1) denotes the 𝑛th natural slosh frequency of the fluid and {𝜓_𝑛} the corresponding mode shape, the free vibration problem to be solved is 𝐾−𝜔²_𝑛𝑀Ψ_𝑛=0. (40) The natural slosh frequencies obtained by solving (40) are compared with theoretical slosh frequencies obtained from (15). The roots 𝜌_𝑛 for various modes of cylindrical tanks and for different 𝑘 values of annular tanks can be obtained from [21]. For the present problem, the root values for the first five modes for cylindrical tanks are 𝜌_𝑛=1.841, 5.335, 8.535, 11.205, 14.850 and for annular tank (𝑘=0.2) are 𝜌_𝑛=1.7051, 4.9608, 8.4331, 12.1650, 15.9932. Table 1 shows the first four slosh frequencies in rad/s obtained for the present tanks using finite element method and the pervious analytical formula. Both the results are in good match.

Table 1 

Slosh frequencies (rad/s) in axisymmetric tanks.

5.2. Vertically Excited Tanks

The free surface slosh response is evaluated using the finite element numerical formulation discussed previous for different vertical harmonic excitation amplitudes and frequencies lying inside and outside the regions of parametric instability. In the numerical simulations, 20 nodes along the 𝑥-direction and 20 nodes along the 𝑧-direction are taken, and a time step of 0.01 s is adopted. The different cases considered in simulation are marked on the stability chart as shown in Figure 4. The excitation parameters for the cases shown in Figure 4 are given in Table 2. The test cases considered are similar to the cases considered by Dodge et al. [7].

Table 2 

Excitation parameters for the test cases shown in Figure 4.

Figure 4 

Stability chart for sloshing response under vertical excitations with test cases marked.

To simulate the vertical slosh response the initial boundary condition on the free surface is assumed to be 𝜉=𝜉₀𝑅cos𝜋𝑟𝑅, (41) where 𝜉₀=0.002.

Figures 5(a) and 5(b) show the free-surface elevation at 𝑟=𝑅 for case 1 for cylindrical and annular tanks, respectively, and associated phase plane plots. As case 1 lies in stable region, the slosh response is stable. The time histories of the free surface elevation are nondimensionalised. It can be observed from the phase plane plot that the solution is bounded.

(a)

(b)

(c)

(d)

Figure 5 

Slosh response and associated phase-plane plot for case 1.

Figures 6(a) and 6(b) show the free-surface elevation at 𝑟=𝑅 for case 2 for cylindrical and annular tanks respectively, and associated phase plane plots. As case 2 lies in unstable region, the slosh response is unstable. The free surface undergoes instability; exponential growth in the amplitude of free surface can be seen. It can be observed from the phase plane plot that the solution is unbounded.

(a)

(b)

(c)

(d)

Figure 6 

Slosh response and associated phase-plane plot for case 2.

Figures 7(a) and 7(b) show the free-surface elevation at 𝑟=𝑅 for case 3 for cylindrical and annular tanks, respectively. Case 3 also lies in unstable region, but in second instability region. This case corresponds to instability in the sloshing mode lying in the second instability region. According to theory, the effect of parametric resonance gradually reduces as we move to higher regions of instability. As expected, the amplitudes do not grow rapidly in this instability region compared to first instability region response shown in Figures 6(a) and 6(b). In case of annular tank, first the amplitude of the slosh response started growing exponentially in a resonance mode and then after certain time the response reduced gradually. As the amplitude increases the natural frequency of the system changes and creates low frequency amplitude oscillations leading to decrease in amplitudes of response. This behaviour is called detuning effect; under parametric excitation of frequency close to twice the natural frequency of a certain mode, the free surface oscillates exhibiting the shape of that mode. As the excitation amplitude increases, the natural frequency changes and the input energy can excite the other neighbour modes. If the excited neighbour nodes are stable, the increase in the amplitude will be suppressed leading to detuning effect. This detuning effect can be captured only in non-linear systems. In case of linear systems [5], the response will be always increasing; this detuning effect cannot be captured. The present finite element non-linear numerical model can capture this detuning effect effectively. No detuning effect was found in cylindrical tank for the present case, the amplitude increases exponentially, but the magnitude is lower compared to Figure 6(a), and time taken to build up instability is large compared to Figure 6(a). Figures 7(c) and 7(d) show the respective phase-plane plots. Phase-plane plot for cylinder shows a very slow buildup of instability.

(a)

(b)

(c)

(d)

Figure 7 

Slosh response and associated phase-plane plot for case 3.

Figures 8(a) and 8(b) show the slosh response at 𝑟=𝑅 for case 4. Case 4 corresponds to instability lying in first instability region with lower excitation amplitude compared to case 2. As the excitation amplitude is low, the rate of increase of amplitude and magnitude of amplitude will be less than the response compared to case 2.

(a)

(b)

Figure 8 

Slosh response for case 4.

Figures 9(a) and 9(b) show the slosh response for case 5. Case 5 lies in stable region just outside the instability region. As expected, a stable response is obtained.

(a)

(b)

Figure 9 

Slosh response for case 5.

Figures 10(a) and 10(b) show the slosh response for case 6. The parameters for case 6 are close to the parameters of case 5, but it lies just inside the instability region. Although the parameters of case 5 and case 6 are close, as expected the response of the free surface is unstable. For this case the excitation amplitude is high compared to cases 2 and 4 and thus a rapid increase in the slosh response can be found compared to Figures 6 and 8.

(a)

(b)

Figure 10 

Slosh response for case 6.

5.3. Horizontally and Vertically Excited Tanks

In this section the sloshing response when the tank is subjected to combined horizontal and vertical excitations is considered. The tank is assumed to be excited periodically with 𝑧_𝑡(𝑡)=𝑎_𝑣cos(𝜔_𝑣𝑡), along vertical direction, and 𝑟_𝑡(𝑡)=𝑎_ℎcos(𝜔_ℎ𝑡), along the horizontal direction. The initial conditions required for this simulation are 𝜙()𝜉()𝑟,0,0=0,𝑟,0=0. (42) No initial perturbation is required as needed in pure vertical excitation; the horizontal harmonic excitation creates the perturbation needed for slosh response under vertical excitations. The equation for the free surface obtained from (16) for pure vertical excitations differs in the case of combined excitations; a forcing term appears on the right side due to horizontal excitation. Equation (16) under combined excitation of tank can be written as 𝑑²𝑧_𝑛(𝑇)𝑑𝑇²+Ω²_𝑛1−𝑘_𝑣cos(2𝑇)𝑧_𝑛(𝑇)=𝑑²𝑟_𝑡𝑑𝑡². (43) Equation (43) is a nonhomogenous Mathieu Hill equation. It should be noted that the stability chart of the Mathieu- Hill equation is independent of the term on right-hand side. The stability chart shown in Figure 2 is still valid, but the response of the free surface is affected by the presence of horizontal loading. This horizontal term can produce resonance, which is recognized by the linear growth of amplitude in time. It is known that, under horizontal excitation alone, when external excitation frequency is equal to fundamental sloshing frequency the free surface undergoes resonance. To verify this, response of the free surface is simulated under horizontal excitation alone with external horizontal frequency equal to fundamental sloshing frequency. Figures 11(a) and 11(b) show free-surface response for cylindrical and annular container, respectively, under pure horizontal excitation when excitation frequency is equal to fundamental frequency. As expected the response of the free surface shows resonance phenomenon.

(a)

(b)

Figure 11 

Slosh response for horizontal excitation alone; 𝜔ℎ=𝜔𝑛(1), 𝑎ℎ=0.002 m.

In case of pure horizontal motion the system has only one resonance frequency; but under combined motion, system has infinite resonance frequencies. When the horizontal excitation frequency 𝜔_ℎ is close to fundamental slosh frequency and when sum or difference of horizontal and vertical frequencies 𝜔_ℎ,𝜔_𝑣 is closer to fundamental slosh frequency system undergoes resonance. This resonance is characterized by linear growth in the amplitude if the vertical excitation parameters are in stable region. If the vertical excitation parameters are in unstable region, system grows exponentially in time. Figures 12–14 show the slosh response of free surface in cylindrical and annular containers for three main resonant frequencies of the horizontal motion (𝜔_ℎ/𝜔₁=0.98, 0.18, 1.78) and under vertical excitation parameters lying in a stable region test for case 1 as shown in Figure 4.

(a)

(b)

Figure 12 

Slosh response for the main resonance with horizontal excitation 𝜔ℎ/𝜔1=0.98, 𝑎ℎ=0.002 m with vertical excitation parameters of case 1.

(a)

(b)

Figure 13 

Slosh response for the first resonance with horizontal excitation 𝜔ℎ/𝜔1=0.18, 𝑎ℎ=0.002 m with vertical excitation parameters of case 1.

(a)

(b)

Figure 14 

Slosh response for the second resonance with horizontal excitation 𝜔ℎ/𝜔1=1.78, 𝑎ℎ=0.002 m with vertical excitation parameters of case 1.

Figure 12 shows the slosh response for the strongest of the resonant frequencies in cylindrical and annular tanks. The horizontal forcing frequency is closer to the first slosh natural frequency (𝜔_ℎ/𝜔₁=0.98). Figures 13 and 14 show the slosh response for coupled frequencies (𝜔_ℎ±𝜔_𝑣) closer to first slosh natural frequency. It can be observed from the figures that the sloshing response is high in main resonant frequencies compared to secondary resonances. Influence of vertical excitation with horizontal excitation on slosh response is that if the vertical excitations parameters are in stable region and the horizontal frequency or coupled frequencies are closer to sloshing natural frequency resonance takes place which is characterized by linear increase in the response.

Figures 15–18 show the slosh response in the unstable regions in the presence of horizontal forcing excitation. The vertical excitation parameters lie in unstable regions, and a horizontal forcing amplitude 𝑎_ℎ=0.002 m with a forcing frequency 𝜔_ℎ/𝜔₁=0.5 is considered. The difference between the present simulation and simulations carried out in Section 5.2 is that in the present case tank is also excited horizontally. This horizontal forcing creates necessary initial perturbation required for vertical slosh response. Figures 15(a) and 15(b) show the slosh response in cylindrical and annular tanks with vertical parameters of case 2. This corresponds to the instability of the first mode in the first instability region. The sloshing response of the combined forced tank motion, compared with slosh response in case of pure vertical excitation shown in Figure 6, illustrates that the horizontal forcing excitation delays the occurrence of instability in the free surface. Figures 16(a) and 16(b) show the slosh response in cylindrical and annular tanks with vertical excitation parameters of case 3. This corresponds to the instability of the first mode in the second instability region and thus shows a less strong parametric resonance compared with response in Figure 15. On comparing this combined response with pure vertical excitation response shown in Figure 7, one can infer that detuning effect takes place in cylindrical tank, whereas the present detuning effect is lost in annular tank displaying instability. Thus the horizontal excitation with vertical parameters in unstable region can make free-surface undergo detuning and can cancel off the existed detuning effect. Figures 17 and 18 show the slosh response in cylindrical and annular tanks with vertical excitation parameters of case 4 and case 5, respectively. Figures 8 and 9 show slosh response in case 4 and case 5 for pure vertical excitation. From the figures we can infer that, the effect of horizontal excitation is to delay the occurrence of instability.

(a)

(b)

Figure 15 

Slosh response under combined excitation with vertical excitation parameters of case 2.

(a)

(b)

Figure 16 

Slosh response under combined excitation with vertical excitation parameters of case 3.

(a)

(b)

Figure 17 

Slosh response under combined excitation with vertical excitation parameters of case 4.

(a)

(b)

Figure 18 

Slosh response under combined excitation with vertical excitation parameters of case 6.

6. Conclusion

Sloshing response in axisymmetric tanks under pure vertical excitations and combined horizontal and vertical excitations is analyzed. In these cases, free surface boundary condition reduces to a Mathieu equation; free surface undergoes parametric instability for some combinations of vertical excitation frequencies and amplitudes. Stability chart is plotted for the dynamic stability analysis from linear equations. A fully non-linear finite element numerical model has been developed based on the potential flow theory for simulating the slosh response. The slosh response is simulated for different cases lying in the stability chart. In the stable regions, the free surface response is always bounded. In the unstable regions, the free-surface undergoes parametric resonance characterized by an unbounded response of the free-surface. In the unstable regions, even small excitations can cause the growth of small initial perturbations, if the tank is excited for a sufficiently long time. The slosh response obtained is in exact agreement with the theoretical predictions of the stability chart. The present numerical model can capture detuning effects as well. Detuning is due to change of frequencies during the amplitude growth.

Sloshing response under combined horizontal and vertical base excitations is carried out. Slosh response for combined excitation is unstable if the vertical excitation parameters lie in unstable region. For vertical excitation parameters lying in the stable region, it is found that system undergoes resonance, when coupled frequencies are equal to slosh frequencies. In addition to the resonant frequency of the pure horizontal excitation, there exist an infinite number of additional resonance frequencies due to the combined motion of the tank. The horizontal excitation, when the vertical excitation parameters lie in unstable region, delays the instability and can introduce detuning or cancel detuning leading to slow buildup of instability.