Finite Element Simulation of Dynamic Stability of Plane Free-Surface of a Liquid under Vertical Excitation

Abstract

When partially filled liquid containers are excited vertically, the plane free-surface of the liquid can be stable or unstable depending on the amplitude and frequency of the external excitation. For some combinations of amplitude and frequency, the free-surface undergoes unbounded motion leading to instability called parametric instability or parametric resonance, and, for few other combinations, the free-surface undergoes bounded stable motion. In parametric resonance, a small initial perturbation on the free-surface can build up unboundedly even for small external excitation, if the excitation acts on the tank for sufficiently long time. In this paper, the stability of the plane free-surface is investigated by numerical simulation. Stability chart for the governing Mathieu equation is plotted analytically using linear equations. Applying fully nonlinear finite element method based on nonlinear potential theory, the response of the plane free-surface is simulated for various cases.

1. Introduction

In many engineering applications, it is of great interest to study the oscillatory motion of a mechanical system under external excitations. In all of the vibrating systems, it is essential to know the response of the system and the nature of the oscillations. It is very important to get information about the different resonance cases of the system to avoid the harmful ones in its design. Under external periodic excitations, the system may undergo two kinds of oscillations: forced oscillations and parametric oscillations. Forced oscillations correspond to the oscillatory response of the system in the direction of external excitation, and the system undergoes resonance when external excitation is equal to natural frequency of the system. In resonance, systems amplitude increases linearly. Parametric oscillations refer to an oscillatory motion in a mechanical system due to time-varying changes in the system parameters caused due to external excitation. The system parameters can be inertia or stiffness. The response of the system is orthogonal to the direction of external excitation. System undergoes parametric resonance when the external excitation is equal to integral multiple of natural frequency of the system. In parametric resonance, systems amplitude increases exponentially and may grow without limit. This exponential unlimited increase of amplitude is potentially dangerous to the system. Parametric resonance is also known as parametric instability or dynamic instability. Although parametric instability is secondary, the system may undergo failure near the critical frequencies of parametric instability. Parametric instability is of concern in structures, aircrafts, and marine crafts. Structural systems when subjected to vertical ground motion may undergo parametric resonance leading to large motion or deformation of structure due to which structure may damage or collapse. The problem of parametric instability in various structures is studied and well documented by Bolotin [1]. Parametric resonance also occurs in case of liquid sloshing which is the present interest of the paper.

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 water flow on the deck of the 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 behavior 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 sloshing 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 parametric instability, and, for few other combinations, the free surface undergoes bounded motion. In the present paper, the response of the fluid-filled containers under vertical excitation is undertaken.

The problem of liquid response under vertical excitations was first studied experimentally by Faraday [2] 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 as to Faraday waves. Matthiessen [3] conducted experiments and reported that the fluid free surface vibrations are synchronous to the external excitation. The Faraday study has been analyzed by Rayleigh [4, 5], and he confirmed Faraday’s observations. The discrepancy between Faraday’s observations and Matthiessen’s observations was explained mathematically by Benjamin and Ursell [6]. Benjamin and Ursell [6] 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 the Mathieu equation. The solution for the Mathieu equation [7] 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 [6] 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. Dodge et al. [8] analyzed the sloshing response in propellant tanks of space vehicles theoretically and experimentally under vertical excitations. Mitchell [9] extended the studies on stability of free-surface of liquid under vertical random excitations. Khandelwal and Nigam [10] have studied the parametric instability in rectangular tanks analytically. Miles [11] studied the parametric resonance considering nonlinear terms in cylindrical tanks theoretically. The problem of parametric oscillations was discussed by Miles and Henderson [12] in their review paper. The study of Faraday waves is also reported in works by Perlin and Schultz [13] and Jiang et al. [14]. Most of the studies available in the old literature on the problem of parametric sloshing under vertical excitation are experimental or theoretical 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. Wu et al. [15] applied finite element method for solving sloshing problem. Wu et al. also discussed the sloshing response under vertical excitations in rectangular tanks using a 3D model. Frandsen [16] analyzed the problem numerically and theoretically considering fully nonlinear inviscid potential flow equations. Frandsen applied modified sigma-transformed finite difference method for sloshing response. Frandsen study was on 2-dimensional rectangular tanks. Ning et al. [17] applied boundary element method to study the liquid sloshing in rectangular containers under coupled horizontal and vertical excitation. In the present paper, the sloshing response under vertical excitations in rectangular tank is taken up. First, the stability of plane free surface of liquid in rectangular tank is obtained from the linearized equations, and, then, the obtained stability chart is validated numerically applying finite element method. The objective of the present paper is to study the stability of free-surface of liquid in rectangular tank under vertical excitations numerically. The tank is assumed to be rigid with aspect ratio (ℎ/𝐿), of 0.5; ℎ is depth of fluid, and 𝐿 is length of the tank. First, the stability of the free-surface of the liquid is analyzed theoretically through the governing linearized equations. Stability chart for the liquid free-surface is plotted from these linearized equations using harmonic balance method. The plotted stability chart is validated numerically. A finite element numerical formulation based on nonlinear potential theory is developed for simulating sloshing response.

2. Governing Equations

Consider a rectangular tank fixed in the Cartesian coordinate system 𝑂𝑥𝑧, which is moving with respect to the inertial Cartesian coordinate system 𝑂₀𝑥₀𝑧₀ vertically along the 𝑧-direction. The origins of these coordinates systems are at the left end of the tank wall at the free surface and pointing upwards to the 𝑧-direction. These two Cartesian coordinate systems coincide when the tank is at rest. Figure 1 shows the tank in the moving Cartesian coordinate system 𝑂𝑥𝑧 along with the prescribed boundary conditions.

Figure 1 

Sloshing tank and its boundary condition.

Let the displacement of the tank be governed by the directions of axes as 𝑧=𝑧_𝑡(𝑡). (1) Fluid is assumed to be inviscid, incompressible, and irrotational. Therefore, the fluid motion is governed by Laplace’s equation with the unknown velocity potential Φ as follows: ∇²𝜙=0. (2) Fluid obeys the Neumann boundary conditions on the walls of the container and the Dirichlet boundary condition at the liquid free surface. In the moving coordinate system, the velocity component of the fluid normal to the walls is zero. Hence, at the bottom and on the walls of the tank (Γ_𝐵), we have that 𝜕𝜙𝜕𝑛_𝑥=0,𝐿=0,𝜕𝜙𝜕𝑛_𝑧=−ℎ=0. (3) On the free surface (Γ_𝑆), dynamic and kinematic boundary conditions hold, and they are given by 𝜕𝜙𝜕𝑡_𝑧=𝜁+12∇Φ⋅∇Φ+𝑔+𝑧_𝑡^′′𝜁=0,𝜕𝜁+𝜕𝑡𝜕𝜙𝜕𝑥𝜕𝜁−𝜕𝑥𝜕𝜙𝜕𝑧=0, (4)(5) where 𝜁 is the free surface elevation, measured vertically above still water level, 𝑧_𝑡^′′ is the vertical acceleration of the tank, and 𝑔 is the acceleration due to gravity.

Equations (1)–(5) give the complete behaviour of nonlinear sloshing in fluids under vertical base excitation of the tank. The position of the fluid free surface is not known a priori; to solve the problem, the fluid is assumed to be at rest with some initial perturbation on the free surface. Thus, the initial conditions for the free surface in the moving Cartesian system at 𝑡=0 and 𝑧=0 can be written as Φ(𝑥,0,0)=−𝑧𝑑𝑧_𝑡(𝑡),𝜁()𝑑𝑡𝑥,0=𝜁₀, (6)(7) where 𝜁₀ is the initial elevation of the free surface. It should be noted that it is not possible to attain the initial boundary condition (7) maintaining (6); in reality, it is a nonphysical condition. This condition is used in case of vertical excitations alone, because, in this excitation, some initial perturbation is needed without which there will not be any oscillation in the fluid free surface.

3. Governing Equation for Dynamic Stability of Free-Surface

In this section, the governing equation for dynamic stability of free-surface of liquid under vertical excitation is derived. The general solution for the Laplace equation in the rectangular domain satisfying the rigid boundary conditions (5) can be written as 𝜙=^∞_𝑛=0cosh𝑘_𝑛(𝑧+ℎ)cosh𝑘_𝑛ℎcos𝑘_𝑛𝑥𝐹_𝑛(𝑡),𝜁=^∞_𝑛=0cos𝑘_𝑛𝑥𝑧_𝑛(𝑡), (8) where 𝑘_𝑛=𝑛𝜋/𝐿 is the wave number for mode number 𝑛. 𝐹_𝑛(𝑡), 𝑍_𝑛(𝑡) are the time evolution functions of the respective 𝑛th mode and can be calculated by substituting the general solution (8) in the linear free-surface boundary conditions obtained from (4) and (5). The linearized free-surface boundary conditions are 𝜕𝜙𝜕𝑡_𝑧=𝜁+𝑔+𝑧_𝑡^′′𝜁=0,𝜕𝜁−𝜕𝑡𝜕𝜙𝜕𝑧=0. (9) Substituting (8) into (9) leads to 𝑑𝐹_𝑛(𝑡)+𝑑𝑡𝑔+𝑧_𝑡^′′𝑧_𝑛(𝑡)=0,𝑑𝑧_𝑛(𝑡)𝑑𝑡−𝑘_𝑛tanh𝑘_𝑛ℎ𝐹_𝑛(𝑡)=0. (10)(11) Substituting (10) into (11) and arranging terms, 𝑑²𝑧_𝑛(𝑡)𝑑𝑡²+𝜔²_𝑛𝑧1+^𝑡′′𝑔𝑧_𝑛(𝑡)=0, (12) where 𝜔_𝑛=𝑔𝑘_𝑛tanh𝑘_𝑛ℎ. (13) Equation (13) gives the linear sloshing frequencies. The tank is assumed to be subjected to harmonic vertical excitation given as 𝑧_𝑡(𝑡)=𝑎_𝑣cos𝜔_𝑣𝑡. (14) Equation (11) reduces to 𝑑²𝑧_𝑛(𝑇)𝑑𝑇²+Ω²_𝑛1−𝑘_𝑣cos(2𝑇)𝑧_𝑛(𝑇)=0, (15) where 𝑇=1/2𝜔_𝑣𝑡, Ω_𝑛=𝜔_𝑛/𝜔_𝑣, and 𝑘_𝑣=𝑎_𝑣𝜔²_𝑣/𝑔. Equation (15) is a Mathieu equation. The stability and instability of the free-surface are guided by (15).

The Mathieu equation is a second-order differential equation with periodic coefficients. The solution for the 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 the Mathieu equation is well documented by Mc Lachlan [7], the only reference in the literature which covers exclusively the Mathieu equation. The form of the solution for the Mathieu equation can be obtained by Floquet’s theory [18]. According to the Floquet theory, the solution for the 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 𝐹_𝑟(𝑧)=𝑒^𝑖𝑟𝑧𝑝(𝑧), (16) where 𝑟 is called the characteristic exponent and depends on Ω_𝑛, 𝑘_𝑣 and 𝑝(𝑧) is a periodic function with period 𝜋. The behaviour of the solution can be obtained from (16). Solution for the 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 unstable 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 [1]. The periodic solution with a period 2𝜋 can be written in the form 𝑧(𝑡)=^∞_𝑘=1,3,5𝑎_𝑘sin𝑘𝑡2+𝑏_𝑘cos𝑘𝑡2. (17) Substituting the series (17) into (15) 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,…. (18)(19)(20)(21) The periodic solution with period 𝜋 can be expressed in the Fourier series [1] as 𝑧(𝑡)=𝑏₀+^∞_𝑘=2,4,6𝑎_𝑘sin𝑘𝑡2+𝑏_𝑘cos𝑘𝑡2. (22) Substituting the series (22) into (15) 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,…. (23)(24)(25)(26)

The system of linear homogenous algebraic equation (19) and (21) and (24) and (26) 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. (27) Equation (27) gives the determinant obtained from both (19)–(21) combined under the ± sign as follows: 11−Ω²_𝑛−𝑘_𝑣2−𝑘0⋯_𝑣241−Ω²_𝑛−𝑘_𝑣2⋯𝑘0−^𝑣21−16Ω²_𝑛⋯𝑘⋮⋮⋮⋱=0,1−^𝑣200⋯−𝑘_𝑣11−Ω²_𝑛−𝑘_𝑣2𝑘0⋯0−^𝑣241−Ω²_𝑛−𝑘_𝑣2⋯𝑘00−^𝑣21−16Ω²_𝑛⋯⋮⋮⋮⋮⋱=0. (28)

The determinants given above are called the Hill determinants, and they are of infinite order. The Hill determinants are tridiagonal in nature. It is clear from these determinants that Ω_𝑛 always appear on the diagonal of the matrices; one can invoke the analogy with the eigenvalue problem and refer to Ω_𝑛 as an eigenvalue. Then, for given values of 𝑘_𝑣, it is possible to calculate values of Ω_𝑛 corresponding to periodic solutions of periods 𝜋 and 2𝜋. By solving the above eigenvalue problems, a stability chart is plotted. The stability chart for (15) 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 (white region), the free-surface is stable. Thus, from the stability chart, the stability of the free-surface can be predicted.

Figure 2 

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

4. Numerical Procedure

In the present section, the response of the free-surface under vertical excitations of the tank is simulated numerically using finite element method. The finite element model based on the mixed Eulerian-Lagrangian scheme is adopted. The free surface nodes behave like the Lagrangian particles, and the interior nodes behave like the Eulerian particles. For this formulation, the free-surface kinematic and boundary conditions (4) and (5), respectively, are modified and written in the Lagrangian form as [19] 𝑑Φ𝑑𝑡_𝑧=𝜁=12∇Φ⋅∇Φ−𝑔+𝑧_𝑡^′′𝜁,𝑑𝑥=𝑑𝑡𝜕Φ,𝜕𝑥𝑑𝑧=𝑑𝑡𝜕Φ.𝜕𝑧 (29)(30)

The problem of sloshing is nonlinear because the free-surface position is not known a priori and the boundary conditions have nonlinear terms. The sloshing problem is evaluated as an initial boundary value problem. The initial boundary is taken as Φ(𝑥,0,0)=−𝑧𝑑𝑧_𝑡(𝑡)𝜁()𝑑𝑡=0,𝑥,0=𝜁₀=𝑎cos𝑘_𝑛𝑥, (31) where 𝑎 is the amplitude of the initial perturbation. In order to solve this nonlinear sloshing problem, a finite element numerical approach based on the mixed Eulerian-Lagrangian scheme is employed. The free surface nodes behave like the Lagrangian particles, and interior nodes behave like the Eulerian particles. A four-node isoparametric element is used in the analysis. Time interval 𝑡 is divided into finite number of time steps, 𝑡_𝑛=𝑛Δ𝑡 (𝑛=0,1,2,3,…), at a particular time step (𝑛=0); the initial boundary conditions (31) are known; using these initial conditions along with the boundary condition (3), the Laplace equation, (2), is solved to get velocity potential Φ, with which velocity 𝑣 is evaluated; with these evaluated velocities, the kinematic and dynamic free surface boundary conditions, (29)-(30), 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 velocities from the potential are calculated using patch recovery [20]. The fourth-order Runge-Kutta method is employed to advance the solution in time. As the time proceeds in the simulation due to the Lagrangian behaviour of the free surface nodes, these nodes move closer and develop a steep gradient leading to numerical instability. To get rid of this problem, a cubic spline interpolation is used to regrid the free-surface uniformly. The authors have developed a finite element numerical formulation for nonlinear sloshing response for sloshing in 2D rectangular tanks and vertical axisymmetric tanks under horizontal excitations and vertical coupled excitations [21, 22]. The same numerical formulation is followed here; the present simulation is an extension of [21, 22] sloshing response under vertical excitations. In the present paper, the numerical formulation is explained briefly; for complete in depth explanation and algorithm, readers can refer to [21, 22].

4.1. Finite Element Formulation

The solution for the nonlinear sloshing boundary value problem is obtained using finite elementmethod. The entire liquid domain is discretized by using four-node isoparametric quadrilateral elements. A typical mesh structure assuming some initial perturbation on the free-surface is shown in Figure 3. By introducing the finite element shape functions, the liquid velocity potential can be approximated as 𝜙(𝑥,𝑧)=^𝑛_𝑗=1𝑁_𝑗(𝑥,𝑧)𝜙_𝑗, (32) where 𝑁_𝑗 is the shape function, 𝑛 is the number of nodes in the element, and 𝜙_𝑗 is the nodal velocity potential.

Figure 3 

Finite element discretization of liquid domain.

The potential on the free surface at a particular time step is obtained from the free surface boundary condition (like at 𝑡=0, (6)). It is needed to calculate the velocity potential for the interior nodes. Applying the Galerkin residual method to the Laplace equation along with the Neumann boundary conditions and taking the free surface nodes, where the potential is known to the right-hand side, will lead to the following system of finite element equations: _Ω∇𝑁_𝑖^𝑚_𝑗=1Φ_𝑗∇𝑁_𝑗𝑑Ω|_{𝑗,𝑖∉Γ𝑠}=−_Ω∇𝑁_𝑖^𝑚_𝑗=1𝜙_𝑗∇𝑁_𝑗𝑑Ω|_{𝑗∈Γ𝑠,𝑖∉Γ𝑠}, (33) where 𝑚 is the total number of nodes in the liquid domain.

4.2. Velocity Recovery

To track the free surface, (29) and (30) need velocities which are computed from the calculated potential field using 𝑣=∇𝑁_𝑗⋅Φ_𝑗. (34) The velocities calculated using (34) are the velocities at the Gauss integration points, and they do not possess interelement continuity and have a low accuracy at nodes and element boundaries. Utmost care should be taken to calculate the velocities; a small error in the velocity recovery will affect the accuracy of free surface updating or tracking, get accumulated with time, and lead to underestimation of the solution. In order to derive a smoothed 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 as follows: 𝑣̂=𝑎₁+𝑎₂𝜉+𝑎₃𝜂+𝑎₄𝜉𝜂, (35) where 𝑣̂ is any velocity component (𝑣_𝑥̂ or 𝑣_𝑦̂), 𝜉, 𝜂 are the Gauss locations, and 𝑎₁, 𝑎₂, 𝑎₃, and 𝑎₄ are unknowns which need to be evaluated. To evaluate these unknowns, a least square fit is considered between 𝑣̂ and 𝑣 as follows: 𝐹(𝑎)=⁴_𝑖=1𝑣̂𝜉_𝑖,𝜂_𝑖−𝑣𝜉_𝑖,𝜂_𝑖², (36) where 𝑖 is the 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 these calculated 𝑎_𝑘’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 the nodal averaged velocities. The global continuous velocity field is given as 𝑣=𝑁⋅𝑣,̂ (38) where 𝑣 is a velocity component (𝑣_𝑥 or 𝑣_𝑦).

4.3. Numerical Time Integration and Free Surface Tracking

After calculating the velocity at a time step 𝑡, it is required to calculate the position of free surface from (30) and determine the potential on the free surface using (29) for the next time step 𝑡+Δ𝑡. As a result, the liquid mesh and the boundary condition required for the next time step are established. This is done using a finite difference numerical procedure. The numerical time integration scheme plays a major role in any time marching problem. The fourth-order Runge-Kutta explicit time integration method is employed in the present paper. The nodal coordinates of the free surface and the associated velocity potential at a current time step 𝑖 are known and can be represented in a single variable as 𝑠_𝑖=𝑥_𝑖,𝑧_𝑖,Φ_𝑖, (39) where 𝑥_𝑖=𝑥₁,𝑥₂,…,𝑥_{𝑁𝑋+1𝑖},𝑧_𝑖=𝑧₁,𝑧₂,…,𝑧_{𝑁𝑋+1𝑖},Φ_𝑖=Φ₁,Φ₂,…,Φ_{𝑁𝑋+1𝑖}, (40) where 𝑁𝑋 is number of segments along the free surface. Similarly, the time derivative can be written as 𝐷𝑠_𝑖𝑡𝐷𝑡=𝐹_𝑖,𝑠_𝑖=𝐹_𝑖. (41) The free surface position and the associated velocity potential at the next time step 𝑖+1 can be expressed as 𝑠_𝑖+1=𝑠_𝑖+𝑠₁6+𝑠₂3+𝑠₃3+𝑠₄6, (42) where 𝑠₁=∇𝑡𝐹𝑡_𝑖,𝑠_𝑖,𝑠₂=∇𝑡𝐹𝑡_𝑖+∇𝑡2,𝑠_𝑖+𝑠₁2,𝑠₃=∇𝑡𝐹𝑡_𝑖+∇𝑡2,𝑠_𝑖+𝑠₂2,𝑠₄=∇𝑡𝐹𝑡_𝑖+∇𝑡,𝑠_𝑖+𝑠₃. (43) After obtaining the new positions and potential of the free surface, the liquid domain is remeshed based on these obtained new coordinate positions.

4.4. Regridding Algorithm

At the beginning of the numerical simulation, the free surface nodes are uniformly distributed along the 𝑥-direction. As the time proceeds, the free surface nodes are spaced unequally, and they 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 B-spline interpolation is employed when the movement of the nodes is 75% more or less than the initial grid spacing. For the regridding, first, the free surface length 𝐿_𝑓 is obtained. Then, the free surface is divided into 𝑁𝑋 segments with the identical arc length. The coordinates of node are denoted as (𝑥_𝑖,𝑦_𝑖) (𝑖=1,2,…,𝑁𝑋+1), and let the arc length between two successive points 𝑖 and 𝑖+1 be 𝑆_𝑖. Being a uniform regridding, 𝑆_𝑖 can be expressed as 𝑆_𝑖=𝑖𝐿_𝑓.𝑁𝑋 (44) The coordinates of the nodes (𝑥_𝑖,𝑦_𝑖) are a function of the arc length 𝑆_𝑖: 𝑥_𝑖,𝑧_𝑖=𝑓𝑆_𝑖. (45)

The cubic B-spline interpolation is used to calculate the coordinates (𝑥_𝑖,𝑧_𝑖), and the velocity potential on the new uniform free surface is also obtained in a similar fashion.

4.5. Complete Algorithm for Sloshing Response

Including all of the steps above, the algorithm for numerical procedure of nonlinear sloshing is as follows.(1)Assume the initial free surface elevation and the associated velocity potential from (31). (2)The liquid domain is discretized using finite elements, and the nodal connectivity data and nodal positions data are obtained.(3)Velocity potential for the interior nodes is calculated using (33). (4)Recover the velocity from the velocity potential using patch recovery technique.(5)Update the free surface nodes using (29) and (30) to a new position, and find the new associated velocity potential using the fourth-order Runge-Kutta method. (6)Check for the requirement of regridding. (7)After finding the new free surface position and the associated velocity potential, repeat step 2 until the final time step is reached.

5. Numerical Results and Discussion

A code is developed following the pervious numerical formulation for simulating the dynamic stability of free-surface sloshing response under various external vertical excitation amplitudes and frequencies. The liquid sloshing response inside a 2.0 m wide rectangular rigid container with liquid filled to a depth of 1 m is simulated. Aspect ratio of 0.5 is maintained. For the present tank, the order of sloshing frequencies from (13) is 𝜔₁=3.7594rad/s and 𝜔₃=6.7986rad/s. The tank is assumed to be subjected to the vertical harmonic excitation given by (14). In the numerical simulations, 40 nodes along the 𝑥-direction and 20 nodes along the 𝑧-direction are taken, and a time step of 0.01 s is adopted. The sloshing response is simulated for different cases with increasing wave steepness, inside and outside the regions of parametric instability or parametric resonance. The steepness parameter depends on the initial conditions considered: 𝜀=𝑎−𝜔²_𝑛/𝑔 [16]. 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 1. The test cases considered are similar to the cases considered by Frandsen [16].

Table 1 

Excitation parameters for the test cases shown in Figure 4.

Figure 4 

Stability chart of sloshing response under vertical excitations (test cases are also shown).

The first test case is in stable region, with frequency ratio Ω₁=1.253 and forcing amplitude 𝑘_𝑣=0.5: test case 1 as shown in Figure 4. The sloshing responses on the left wall of the tank for low (𝜀=0.0014) and high (𝜀=0.288) wave steepness are shown in Figures 5(a) and 5(b), respectively. Figures 5(c) and 5(d) show the respective phase-plane plots for the small and steep wave cases. The time histories of the free-surface elevation are nondimensionalized. The sloshing response obtained with the present simulation is compared with numerical results of Frandsen [16]. Both results are in excellent agreement. The sloshing response for low steep waves is symmetric; that is, amplitudes of peaks and troughs are equal, whereas, for high steep waves, the response is asymmetric showing different amplitudes for peaks and troughs. This is an indication of nonlinear response. This nonlinear behaviour can be noticed from the respective phase-plane plots. The phase-plane plot of the sloshing response with low wave steepness shown in Figure 5(c) has a closed orbit displaying a linear behaviour, whereas the phase-plane plot of the sloshing response with high wave steepness has nonrepeatable and nonclosed orbits displaying a nonlinear characteristic.

(a)

(b)

(c)

(d)

Figure 5 

Free-surface elevation on the left wall in the stable region (test case 1): Ω1=1.253, 𝑘𝑣=0.5 for (a) 𝜀=0.0014, (b) 𝜀=0.288, and the respective phase-plane plots.

The second test case lies in the unstable region, with frequency ratio Ω₁=0.5 and forcing amplitude 𝑘_𝑣=0.3: test case 2 as shown in Figure 4. Figure 6(a) shows the free surface elevation on the left wall of the tank for a low wave steepness of 𝜀=0.0014. As the excitation parameters lie in unstable region, an unbounded response is expected; the sloshing response plot displays the expected behaviour. Figure 6(b) shows the corresponding phase-plane diagram. The phase-plane plot and the response plot clearly show that the nonlinear effects are predominant. A moving mesh generated at different time steps in this unstable region for parametric resonance is shown in Figure 7. The animation of simulation in this case can be seen in supplementary file Dynamic Instability.gif; see the supplementary file in Supplementary Material available online at http://dx.doi.org/10.1155/2013/252760.

(a)

(b)

Figure 6 

Free-surface elevation on the left wall in the unstable region (test case 2) for Ω1=0.5, 𝑘𝑣=0.5, 𝜀=0.0014, and the respective phase-plane plot.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 7 

Moving mesh generated at different type steps in parametric resonance.

Figure 8 also shows the sloshing response time histories in unstable regions. A low wave steepness parameter of 𝜀=0.0014 is considered. Figure 8(a) shows the sloshing response time history for frequency ratio Ω₁=1.0 and forcing amplitude 𝑘_𝑣=0.5: test case 3 as shown in Figure 4. This case corresponds to instability in the first sloshing mode lying in the second instability region. According to the theory, the effect of parametric resonance gradually reduces as we move to the higher regions of instability. As expected, the amplitudes do not grow rapidly in this instability region compared to the first instability region response shown in Figure 6(a). First, the amplitude of the sloshing response starts growing exponentially in a resonance mode, and, then, after certain time, the response reduces 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 neighbor nodes. If the excited neighbor nodes are stable, the increase in the amplitude will be suppressed leading to detuning effect. This detuning effect can be captured only in nonlinear systems. In case of linear systems [6], the response will be always increasing; this detuning effect cannot be captured. The present finite element nonlinear numerical model can capture this detuning effect effectively. Figure 8(c) shows the respective phase-plane plot. Figure 8(b) shows the sloshing response on left wall of tank for frequency ratio Ω₃=0.5 and forcing amplitude 𝑘_𝑣=0.2: test case 4 as shown in Figure 4. This case corresponds to instability in the second sloshing mode lying in the first instability region. As the instability is in the second mode, the amplitudes do not grow rapidly when compared to instability in the first mode, Figure 6(a). After certain time, the amplitude comes down showing the detuning effect. In this case, the free-surface oscillates exhibiting the third sloshing mode, and, as the amplitude increases, the input parametric excitation excites the first sloshing mode, which is stable, and the amplitudes fall down. Similar free-surface behaviour is reported by Frandsen [16]. Figure 8(d) shows the respective phase-plane plot. The phase-plane plots for the responses display a linear behaviour for the system. Figure 9 shows the moving mesh generated for the response shown in Figure 8(b) at various time steps. The animation of simulation in this case can be seen in supplementary file Detuning Effect.gif.

(a)

(b)

(c)

(d)

Figure 8 

Free-surface elevation on the left wall in the unstable regions with 𝜀=0.0014: for (a) Ω1=1.0, 𝑘𝑣=0.5 (test case 3) (b) Ω3=0.5, 𝑘𝑣=0.4 (test case 4), and the respective phase-plane plots.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 9 

Moving mesh generated at different time steps in unstable region displaying detuning effect.

Figure 10 shows the sloshing response for test case 5 as shown in Figure 4, with frequency ratio Ω₁=0.6 and forcing amplitude 𝑘_𝑣=0.5. This point lies in the stable region but very close to instability region. Figure 10(a) shows the sloshing response for low wave steepness parameter, 𝜀=0.0014, and Figure 10(b) shows the sloshing response for high wave steepness parameter, 𝜀=0.288. As expected, the point is in the stable region, and the sloshing response is stable.

(a)

(b)

Figure 10 

Free-surface elevation on the left wall in the stable region (test case 5) for Ω1=0.55, 𝑘𝑣=0.5: (a) 𝜀=0.0014, and (b) 𝜀=0.288.

Figure 11 shows the sloshing response for test case 6 with frequency ratio Ω₁=0.55 and forcing amplitude 𝑘_𝑣=0.5 lying in the unstable region as shown in Figure 4. A low steepness parameter 𝜀=0.0014 is taken. As the point lies in the unstable region, the response is also unstable as expected. The low steepness response is sufficient to be demonstrated to show the rapid increase in the amplitudes.

Figure 11 

Free-surface elevation on the left wall in unstable region (test case 6) for Ω1=0.55, 𝑘𝑣=0.5, and 𝜀=0.0014.

6. Conclusion

The free-surface of fluid under vertical excitation of the tank reduces to the Mathieu equation. Free-surface undergoes parametric resonance for some combinations of vertical excitation frequencies and amplitudes; this resonance is characterized by exponentially unbounded growth of amplitude. Stability chart is plotted for the dynamic stability analysis from linear equations. A fully nonlinear finite element numerical model has been developed based on the potential flow theory for simulating the sloshing response. The sloshing 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 sloshing 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.

Supplementary Materials

Dynamic Instability.gif shows the animation of sloshing fluid for the second test case which lies in the unstable region, with frequency ratio Ω1 = 0.5 and forcing amplitude kV =0.3: test case 2 as shown in Figure 4. Figure 6 shows the free-surface elevation and Figure 7 shows the moving mesh generated in this case.

Detuning Effect.gif shows the animation of sloshing fluid for the fourth test case, with frequency ratio Ω3 = 0.5 and forcing amplitudeV = 0.2: test case 4 as shown in Figure 4. This case corresponds to instability in the second sloshing mode lying in the first instability region. Figure 8(b) shows the free-surface elevation and Figure 9 shows the moving mesh generated in this case.

1. Supplementary Material