Eswaran & Saha CFD Letters Vol.

4(4) 2012

[Link] Vol. 4 (4) – December 2012

Finite Difference Based Sigma - Transformation

Approach for Liquid Sloshing in a Rectangular
Tank under Regular Wave Excitation
Eswaran M1* and Ujjwal K. Saha2,*
1
Structural and Seismic Engineering Section, Reactor Safety Division
Bhaba Atomic Research Centre, Trombay, Mumbai – 400085, India
2
Department of Mechanical Engineering, Indian Institute of Technology Guwahati
Guwahati -781039, India.

Received: 20/12/2011 – Revised 26/11/2012 – Accepted 02/12/2012

Abstract

In this paper, a fully non-linear finite difference model has been developed based on the
inviscid flow equations, and a simple mapping function was used to remove the time-
dependence of the free surface in the fluid domain. The time-varying fluid surface is
mapped onto a rectangular domain by the Sigma (  )-transformation. This method is a
simple way to simulate non-breaking waves quickly and accurately especially that has a
low steepness. The fluid motion is solved in a unit square mesh in the transformed flow
domain (i.e., computational domain). Difference between the peaks and troughs of waves
are discussed for three different cases of horizontal (surge), vertical (sway) and combined
excitations of off and at resonance frequency of the tank. The spectrum analysis of
horizontally excited tank is presented. The stability and instability regions associated with
vertical and combined excitations conditions are discussed with the plots of free surface
elevation, phase-plane diagram and free surface profile.

Keywords: Sloshing; Surge and sway motion; Free surface; Mapping;  -transformation;
Finite difference method.

1. Introduction

The oscillation of the unrestrained free surface of the liquid in a partially filled container due
to external excitation is called as sloshing. These motions generate severe hydrodynamic loads that
can be dangerous for structural integrity in tanks and will raise the stability problems in rockets,
satellites, LNG ships, trucks and even stationary petroleum containers. Free surface of liquid in a
container attempts to attain the state of the equilibrium for the effective instantaneous acceleration
(gravitational, translational, etc.) felt by the fluid. However, the momentum of the fluid and external
forces on the fluid tank will prevent this state of equilibrium. The knowledge of liquid free surface
natural frequencies is important in the design of liquid containers subjected to different types of
excitation [1, 2]. The dynamic behavior of a free liquid surface depends on the type of excitation

*
Corresponding Authors:
1
Email: eswarm21@[Link] Telephone: +91 22-25591528
2
Email: saha@[Link] Telephone: +91 361 2582663 Fax: ++ 91 (361) 2690762
©2013 All rights reserved. ISSR Journals P II: S 21 80 -13 63 (12 )4 17 3 -X

173
 Finite Difference Based Sigma - Transformation Approach for Liquid Sloshing in a Rectangular Tank under Regular
Wave Excitation

and its frequency content. Civil engineers and seismologists have been studying liquid sloshing
effects on large dams, oil tanks and elevated water towers under ground motion.
In 1951 and 1952, Jacobsen and Ayre [3], and Graham and Rodriquez [4] performed some
basic studies relevant to this topic. Housner [5, 6] developed an analytical method for the
determination of hydrodynamic wall pressures under the assumption that the tank was a rigid
structure fixed at the base and only the fundamental sloshing mode was important. Applications in
the aerospace industry has been reviewed and discussed by Abramson [7] both analytically and
experimentally. Ink-trace experiments were conducted by Cole [8] to study the effect of baffle
thickness in the cylindrical tank. In 1970s and in early 1980s, sloshing phenomena in liquid oil
carrier vehicles have extensively been studied, since slosh-induced loads can cause serious damage
to cargo structure in marine engineering. Externally induced sloshing is studied through interface
location technique by Eswaran et al., [9]. A lot of researchers have done numerical simulation
either by using self made program or by using commercial CFD packages [10-12]. In general, the
fluid used in solving the nonlinear sloshing is assumed to be homogeneous, isotropic, viscous and
exhibits only limited compressibility. Various models and techniques have been used to solve the
problem.
The popular numerical methods like finite difference method, the boundary element method
and the finite element method [13] have been used for the sloshing analysis. Behr and Abraham
[14] used finite element method to solve the free surface between the inclined walls and explained
the difficulties associated with free-surface finite element flow simulations and overcoming
techniques. Recently, Sequentially-Coupled Arterial Fluid–Structure Interaction (SCAFSI)
technique was proposed by Tezduyar et al. [15] and applied to find the blood pressure profile of the
cardiac cycle. The various techniques used to handle the free surface behavior in the sloshing have
been the marker and cell (MAC), the volume of fluid (VOF), etc. [16-21]. Recently, Cruchaga et al.
[22] proposed the 3D remeshing algorithm to avoid the progressive distortion in the distribution of
markers in the domain. Recently, the coordinate transformation technique is also applied to this type
of problem etc. The particular type of coordinate transformation (stretching on vertical direction) is
called  -transformation. Since they all require complex computer programming in order to treat
the time varying free surface boundary and update the computational mesh, the  -transformation
gains popularity because of its simplicity. The  -transformation was applied to nonlinear steep
waves in fixed and base excited tanks by Chern et al. [23] and the waves in relatively deep water
was simulated by Turnbull et al. [24]. Frandsen [25] investigated numerically steep free surface
sloshing in fixed and base-excited rectangular tanks with a focus on moving liquid tank with
horizontal and vertical excitations. Recently, Chen and Nokes [26], and Dai and Xu [27] applied  -
transformation to predict the sloshing effects on horizontal cylindrical container and 2D rectangular
tank, respectively. The first few sloshing frequencies in a vertically accelerated container have been
reported by Eswaran and Saha [28, 29]. Originally,  - transformation was proposed for
meteorological forecasting by Phillips [30]. Later, Blumberg and Mellor [31] and Mellor and
Blumberg [32] applied in the context of oceanic and coastal flows.
In this paper, a fully non-linear model for idealized 2-D waves in a numerical wave tank has
been developed. The  -transformation technique is used to capture the liquid free surface which is
used to map the asymmetric liquid domain onto a rectangle, such that the moving free surface in the
physical plane becomes a fixed line in the computational mapped domain. The fourth order central
difference scheme and the Gauss–Seidel point successive over-relaxation iterative procedure are
used to capture the free surface wave profiles and free surface elevation plots of the fluid domain.
Liquid in a rectangular tank under different regular wave excitations (say horizontal, vertical and
combined excited conditions) is studied in detail. Spectrum analysis of horizontal excited tank; and
the stable and unstable conditions of vertical and combined motions are also discussed. Section 2 is
focused on the mathematical formulations of the present work. Section 3 presents the Mapping
procedure which transfers the physical domain to the computational domain. In section 4, the finite

174
 Eswaran & Saha CFD Letters Vol. 4(4) 2012

difference discritization of the computational domain is presented. The grid independence study and
the behavior of the liquid free surface in fixed tank is studied and discussed in section 5.1 – 5.2. The
free surface elevation of liquid, phase plane diagram, spectrum analysis for horizontal excitation is
discussed in section 5.3. The stability and instability regions of vertically excited and combined
motions (horizontally and vertically) of the tank which will create the interesting fluid free surface
behavior is discussed in section 5.4 and 5.5. Additionally, in section 5.5, the free surface profiles for
different time intervals and the 3 dimensional surface plots are also showed to explain the severity
of the combined motions.

2. Numerical details

A rectangular Cartesian coordinate system is first employed, with origin at the mean free-
surface at the left-hand side of the tank. A 2-D nonlinear wave problem is considered in this case, as
depicted in Figure 1, where  is the free-surface elevation above still water level, b is the length of
the tank, and hs is the still water depth. The fluid in the tank is assumed to be inviscid and
irrotational. On the above assumption that the fluid is governed by potential flow theory, the
velocity potential  satisfies the Laplace equation. The velocity components normal to the fixed
boundaries are zero by definition. The left, the right and the bottom boundary conditions are
indicated by L, R, and B, respectively (Figure 1).

y T1 T2
Mean or Still ζ
Water Level x

Standing Wave

L hs R

Figure 1. Physical domain

The free surfaces occur at the interface between two fluids. Such interfaces require two
boundary conditions to be applied, viz., (i) a kinematic condition that relates the motion of the free
interface to the fluid velocities at the free surface (i.e., T1) and (ii) a dynamic condition which is
concerned with the force balance at the free surface (i.e., T2). Tanks are studied satisfying the
condition that velocity at the bottom and at the side walls is zero. Zero pressure at the free surface
of the fluid is also considered in the analysis. In view of the fact mentioned above, kinematic and
dynamic conditions must be satisfied on the free surface. Therefore, the boundary conditions for the
physical domain are given by,

175
 Finite Difference Based Sigma - Transformation Approach for Liquid Sloshing in a Rectangular Tank under Regular
Wave Excitation

 2  2
  0 in the fluid domain (1)
x 2 y 2

 0 on the side walls (2)
n
   
  on the free surface (3)
t x x y
 1
  .  ( g  Yt )  0 on the free surface (4)
t 2
Here, Yt is the acceleration of the container in the vertical direction which can be neglected
from the free surface dynamic boundary condition for fixed container analysis.
The following quantities are introduced for generating dimensionless governing equations for
the present study,
x y  Y g 1
x   ; y   ;    ; Yt '  t ; t   t ;    (5)
b b A g b A bg
where, g is the acceleration due to gravity, A is the wave amplitude, Y is the acceleration of the
t

container and t is the time. Here x , y ,  , Yt, t  and   represent the dimensionless quantities. Using
Eq. (5) (hereafter, primes are omitted for simplification), the non-dimensional governing equation
and boundary conditions can be written as follows,
 2  2
 0 (6)
x 2 y 2

L: 0 on x=0 (7)
x

R: 0 on x=b (8)
x

B: 0 on y = - h s (9)
y
   
T1:  Eb  on y=  (10)
t x x y
 1
T2:  E b  .  (1  Yt )  0 on y=  (11)
t 2
where Eb is the amplitude-length ratio (=A/b). The Eqns. (6) through (11) form an initial boundary
value problem which is the Laplace equation with non-linear boundary conditions imposed on the
free surface. Here, the non-linearity is significant for two reasons. Primarily, the elevation of the
moving free surface is not known a priori at any given time instant and secondly, the boundary
conditions on the free surface [i.e., Eqns. (10) and (11)] contain second order differential terms.

3. Mapping procedure

The time-varying liquid free surface can be mapped onto a fixed plane surface by the proper
coordinate transformations, called the  -transformation, which prevents the need for free surface
smoothing for the cases considered herein. In this paper,  - transformation is applied in the
horizontal direction which stretches between the left and the right wall and in the vertical direction
which stretches between the moving liquid free surface and the bottom of the liquid container to
convert the moving free-surface physical domain onto a fixed square computational domain. The

176
 Eswaran & Saha CFD Letters Vol. 4(4) 2012

following section discusses transformation technique in 2-D and 3-D containers elaborately. During
transformation, the governing equation and boundary conditions will change appropriately.

3.1. Transformation of 2-D Rectangular Container

Initially, formulations are developed for the fixed container condition, so the horizontal
and vertical excitation terms are neglected from the Eq. 4 in the following formulation part.

3.2. Coordinate Transformation

The first transformation adopts the  -transformation technique to map the liquid
domain onto a rectangle, such that the moving free surface in the physical plane (Figure 2)
becomes a fixed horizontal line in the  -transformed domain (Figure 3). The mapping
function  ( x, t ) is defined as
y  hs
 , (12)
h
where h ( x , t )  E b  ( x, t )  h s , (13)

y
Initial wave profile 1
Still water level

 x

Figure 2. The 2-D Physical domain

(-1,1) y (1,1)

(-1,-1) (1,-1)
Figure 3. The 2-D Computational domain

177
Finite Difference Based Sigma - Transformation Approach for Liquid Sloshing in a Rectangular Tank under Regular
Wave Excitation

Here,  ( x, t ) is the stretching factor, which varies from 0 to 1. The value of  at the
bottom of the container is 0, while at the free surface is 1. The first-order derivatives of  can
be calculated as follows:
   y  hs      
     Eb and   Eb . (14)
t t  h  h t x h x
The potential function  ( x, y, t ) in the physical domain is transformed to the potential
function  ( X ,  , T ) in the  -transformed domain.
x  X , X  2 x  1 ; y  Y , Y  2  1 ; and t  T , T=t; (15)
Using the chain rule, the first set of derivatives of  with respect to x, y and t gets
transformed as
  4  
2  Eb (16)
x X h X Y
 2 
 (17)
y h Y
    
and   Eb (18)
t T h T 
The second set of derivatives of  with respect to x and y gets transformed as,
 2  2  2  2 
2
 4 2
 l1 2
 l 2  l3 (19)
x X Y XY Y
2 2
  4  
2
 2 (20)
y h Y 2
 4 h 
2
 8 h    4 h  2 8  2 h 
where, l1    ; l2    ; l3    h X   h X 2  .
 h X   h X     
Hence, by using the  -transformation, we can derive the new governing equation and
boundary conditions specified on the rectangular  -transformed domain. The governing
equation after the first coordinate transformation is given as:
 2   2  2
 L1  L 2  L3 0 (21)
X 2 Y XY Y 2
 4  h  2 2  2 h   h  1  2 h  2 
where, L1   2    ;
 2 L  2 ; L3   2    (22)
 h  X  h X 2  h X  h  h X  
Boundary conditions after the first coordinate transformation are given as:
 2 h 
L:  0 on X = -1 (23)
X h X Y
 2 h 
R:  0 on X = +1 (24)
X h X Y
 Y  1
B:  0, on (25)
Y
2
  2 8 Eb       
T1:       4 Eb on Y  1 (26)
T  h h  x   Y X X
2 2
       1   
T2:  M1  2 E b   M2      (1  Yt ) , on Y  1 (27)
T Y  X Y   h Y  

178
 Eswaran & Saha CFD Letters Vol. 4(4) 2012

2 h 2 h
where, M1  ; M2  . (28)
h T h X

4. Finite difference discretization in the computational plane

In the current study, the finite difference method is adopted. Assuming the transformed
domain to be rectangular and constructing on it a unit square mesh of uniform grid in  and 
directions respectively, the standard fourth order central difference approximation is given by

  i  2, j  16 i 1, j  30 i. j  16 i1, j   i  2, j

 2  2
 O( 4 ) (29)
12( )
1   i 1, j 1   i 1, j 1  i 1, j 1   i 1, j 1 
        O( 2 , 2 ) (30)
2  2(  ) 2(  ) 
2 2
where  (i, j ) denotes  ( i , j ) ;   ,   and   ,   are the first and second order central
difference operators along  and  directions respectively, and     is the mixed second order
central difference operator. Next to the boundary lines are discretized by second order
discretization. All the boundary equations are explicitly discretized. The discretized algebraic
governing and boundary equations are solved by Gauss–Seidel point successive over-relaxation
iterative procedure.

5. Results and discussion

In this problem, the initial wave profile is considered as  ( x, t )  0  A cos( K n x) and

 n  th
 ( , )   0  0 , where A is the initial wave amplitude, Kn is the wave number    for n mode
 b 
number (n = 0, 1, 2, …..), and x is the distance in the horizontal direction. The initial wave
2
steepness defined for fixed tank studies as E  A n / g , where gravity field g is usually considered
as 9.81 m/s2. Here, the ratio of h s to b is considered as 1:2. The linear natural sloshing frequencies
in the two dimensional rectangular tank are expressed by

 n  K n g tanh( K n hs ) , n=1, 2, 3… (31)

Now X D and YD are switched off for fixed tank condition in dynamic boundary condition. Two
quantities are usually kept in mind in the sloshing studies: the amplitude of the wave and the
excitation frequency. Here, the amplitude is measured by the wave steepness. The relation between
the acceleration and wave steepness is discussed above.

5.1. Independence Study

In order to display that the solution is grid and time independent, simulations have been
performed using different numbers of grid nodes and different values of t as shown in Figures
4 and 5.

The wave profiles along the tank at three different times for the first sloshing mode (n=2)
are plotted. Results for different grid resolutions are shown related to moderate a wave amplitude
(E=0.0338) for time steps of t  0.003 sec. Figure 4 shows the different grid size of 21  21,

179
Finite Difference Based Sigma - Transformation Approach for Liquid Sloshing in a Rectangular Tank under Regular
Wave Excitation

41  21, 41  41 and 41  61. Initially, grid size is increased in the horizontal direction from 21 to
41. As it is a fixed boundary, it is found that there is not much variation on horizontal direction,
and then, the grid points are increased from 21 to 41 and 61 in vertical direction. Increasing the
grid points in the vertical direction is found to be more effective in improving the accuracy than
increasing the grid points in the horizontal direction since it has the moving boundary at top. It
has been found that a grid size of 41  61 and 41  41 and a time step of 0.003 sec and 0.004 sec
provided sufficient accuracy to capture nonlinearities related to steep wave predictions (E >
0.02). There is not much variation in vertical direction above the grid size of 41  41 for moderate
wave steepness, and therefore, the grid size of 41  41 is sufficient for this problem.

0.6
21 x 21
41 x 21
0.4 41 x 41
41 x 61

0.2
/A

-0.2

-0.4

-0.6

0 0.5 1 1.5 2
x / hs

Figure 4. Grid independence study for E=0.033 at time 7.5 sec.

Time = 4.5 Sec, nt = 0.003
1 Time = 4.5 Sec, nt = 0.004
Time = 4.5 Sec, nt = 0.03
Time = 13.5 Sec, nt = 0.003
Time = 13.5 Sec, nt = 0.004
Time = 13.5 Sec, nt = 0.03
0.5
/A

-0.5

-1
0 0.5 1 1.5 2
x / hs

Figure 5. Time independence study for E=0.033 and grid size 41  61.

5.2. Effect of Wave Steepness in Fixed Tank

The wave characteristics include a crest at the top and a trough at the bottom. The
difference in elevation between the crests and trough is the wave height. The distance between
the adjacent crests or the troughs of wave is termed the wavelength. The ratio of wave height to
wavelength is the wave's steepness. While increasing the wave steepness the nonlinearity
increases. The free surface elevations at the left wall, in the middle and at right wall of the tank
are showed in Figure 6 (a) and 6 (c). When wave steepness increases, the considerable changes

180
 Eswaran & Saha CFD Letters Vol. 4(4) 2012

have been observed in the free surface wave profile. On the other hand, the wave phase plane
diagrams are shown in Figures 6 (b) and 6 (d). These profiles are a repeatable pattern which is
observed at left wall of the tank. The wave phase plane diagram for low steepness is almost a
perfect circle as seen from Figure 6(b), while at the same time, increasing the steepness makes
the phase plane diagram to become oval as depicted in 6 (d).

At tank left wall

1.5 At tank middle
At tank right wall
1
1

( / t ) / ( A x 2)
0.5
 /A

0
0

-0.5

-1

-1
-1.5

0 10 20 30 40 50 60
Time x 2 -1 0 1
/A
(a) (b)
2
At tank left wall
At tank middle
1.5 At tank right wall
1

1
( / t ) / ( A x 2)

0.5
 /A

0 0

-0.5

-1

-1
-1.5

-2
0 10 20 30 40 50 60
Time x 2 -1 0 1
/A
(c) (d)
Figure 6. Fixed tank elevation and phase plane diagram for n=2 with grid size 41  61 and
t  0.003 sec: (a) - (b) E=0.0338; (c) - (d) E=0.338.

5.3. Horizontally Excited Tanks

For the horizontally excited tank condition, YD was switched off from the dynamic
boundary condition at top wall. The horizontal excited acceleration is fixed as
X D  ( h Ah cos( h t )) . When the external horizontal forcing frequency is equal to the natural
sloshing frequency of the liquid, the resonance will occur. In this section, the free surface
motions are numerically examined off- and at resonance conditions. The initial wave profile is
considered as  ( , )  0  0 . For initial wave impulse is considered as  ( x, t )  0  A cos( K n x)
and for zero impulse condition is  ( x, t )  0  0 . And Kh is the measure of nonlinearity parameter
2
A
which is calculated from k h  h h , where Ah is the excitation amplitude. The sloshing motion
g
is more violent at the natural frequency of the container when the excitation frequency is equal to
the first mode rather than at the third mode, which is a well-known resonance phenomenon.

181
Finite Difference Based Sigma - Transformation Approach for Liquid Sloshing in a Rectangular Tank under Regular
Wave Excitation

5.3.1. Off and at resonance of horizontally excited tanks

  
The horizontal frequency ratio  x   h  is speckled as 0.7, 0.9 and 1.3 for off-
 n 
resonance and 1 for at resonance conditions. To validate the developed model, the present
data have been compared with the results of Frandsen [25]. The results have been found to be
in good agreement as seen from Figure 7(a). While increasing the excitation frequency, some
interesting free surface elevation shapes are observed (Figures 7a, 7c, 8a and 9a). When the
frequency ratio  x is small, standing waves are observed upto certain level. Then waves are
variations in waves are aroused slowly as shown in Figure 7(a) and while we come around
 x = 0.8 to 0.9 the waves are showing the complete different elevation patterns which has
the cunning points at certain intervals.
3
Present Work
Frandsen 2004 2
2

( / t ) / ( Ah x 1)
1
 / Ah

0 0

-1
-1

-2
-2

-3
0 20 40 60 80 100 -2 -1 0 1 2
Time x 1  / Ah
(a) (b)
10
At tank left wall
15 At tank middle
At tank right wall

10
5
( / t ) / ( Ah x 2)

5
 / Ah

0 0

-5

-10 -5

-15

0 30 60 90 120 -10
Time x 1 -10 -5 0 5 10
 / Ah
(c) (d)
Figure 7. Off-resonance horizontal excited tank elevation and phase plane diagram for n=1
with grid size 41  61 Kh=0.0034 and t  0.003 sec (a) - (b)  x = 0.7 and (c) - (d)  x = 0.9

182
 Eswaran & Saha CFD Letters Vol. 4(4) 2012

At tank left wall

40
At tank middle 40
At tank right wall

20
20

( / t ) / ( Ah x 1 )
/Ah

0
0

-20
-20

-40
-40
0 20 40 60 80
-40 -20 0 20 40
Time x 1  / Ah
(a) (b)
Figure 8. At resonance conditions horizontal excited tank elevation and phase plane
diagram for n=1 with grid size 41  61 Kh=0.0034 and t  0.003 sec.
6
4

3
2
( / t ) / ( Ah x 1 )
/A

0
0

-2
-3

-4

0 20 40 60 80 -6
-6 -3 0 3 6
Time x 1  / Ah
(a) (b)
Figure 9. Off-resonance horizontal excited tank elevation and phase plane diagram for n=1
with grid size 41  61, Kh=0.0034 and t  0.003 sec  x = 1.3

After reaching the  x =1 the free surface elevation are growing continuously with
constant increments (Figure 8). Since, at this point the excitation frequency is matched with
system natural frequency. Figure 8 (a) showed as the free surface elevation is observed 36
times to amplitude of excitation approximately at non-dimensional time 80. As discussed in
section 5.3, in Figure 8(a), free surface elevation is shown with two initial conditions viz. zero
impulse and initial impulse. After few seconds, both the results match very closely. It is
observed that the initial impulse has not undergone much variation from the original
behaviour of the liquid wave profiles during sloshing. However, at low frequencies, one can
expect very small deviations between the initial impulse and the zero impulse. Finally, the
frequency ratio was increased to  x = 1.3, where by the free surface elevation drastically
reduces to 3.4 times (from 36 times at resonance frequency) to the excitation amplitude
approximately at non-dimensional time 80. Figures 7(b), 7(d), 8(b) and 9(b) show the phase
plane diagram for their respective frequencies. In Figure 8(b), the phase plane diagram is
moving in a circular path and finally it attains the spiral shape, since it is drawn at the
resonance condition.

183
Finite Difference Based Sigma - Transformation Approach for Liquid Sloshing in a Rectangular Tank under Regular
Wave Excitation

5.3.2. Spectrum analysis for horizontally excited tanks

The Fast Fourier Transform (FFT) is extremely important in the area of frequency
(spectrum) analysis because it takes a discrete signal in the time domain and transforms that
signal into its discrete frequency domain representation. The FFT does not directly give the
spectrum of a signal. The Shift is required for visualizing the Fourier transform with the
zero-frequency component in the middle of the spectrum.
The spectra of a wave elevation are computed by the FFT, which is given by
N
( j 1)( k 1)
X ( K )   x( j ) N , (32)
j 1

where  n  e ( 2i ) / N , is a Nth root of unity. In order to avoid the effect of discontinuity at
the boundary, the hamming window as given by
 n
 (n)  0.54  0.46 cos 2 ,0  n  N (33)
 N
is used, where the window length is L=N+1. Figure 10 shows the spectra of wave elevations
at the tank left corner. It also caused peaks in the power spectra to become bigger and more
peaks appear at different frequencies. Figure 10 (a) shows that the maximum spectral peak
occurs at the excitation frequency when the excitation frequency is less than the first natural
frequency and a secondary peak occurs at the container natural frequency. While frequency
ratio  x approaches 0.9, the partial merging of first mode frequency and the excitation
frequency slightly increase the energy of the signal as shown in Figure 10 (b). In Figure 10
(c), a single peak of the dimensionless sloshing energy is observed at frequency ratio of 1,
which is irrespective of the magnitude of the excitation frequency. At the first modal
frequency, the response component is observed to increase as the frequency increases. When
the excitation frequency is greater than 1 , the domination of the first mode reduces with an
increase in the frequency ratio. When the excitation frequency is greater than the first mode
(Figure 10 d), the sloshing dominates at the first modal frequency up to second mode (  2 =
1.5 1 ) and till this frequency, the secondary peak is observed at the excitation frequency as
can be seen in Figure 10 (d).

5.3.3. Vertically excited tanks

The initial conditions for vertical excited tank are same as the sloshing motion
simulation of the fixed tank. It is difficult to simulate sloshing only with vertical excitation
by experiments. In order to have an initial perturbation in the free surface inside the
container, horizontal motions need to be excited before the vertical excitation. To avoid this
situation, the initial standing wave profile is assumed for this work. Initial wave impulse is
required for vertically excited condition and it is considered here as
 ( x, t )  0  A cos( K n x) . The vertical acceleration of tank is given by
YD  (  v Av cos( v t )) , where Av is the vertical forcing amplitude, t is the time,  v is the
angular frequency of forced vertical motion. The horizontally excited tank acceleration term
X D is switched off from the dynamic boundary condition at top wall for this analysis. The
initial velocity potential in the fluid domain is considered as  ( , )   0  0 .

184
 Eswaran & Saha CFD Letters Vol. 4(4) 2012
4 4
x 10 x 10
2.5 10

9 h,1
2 8


Power spectral density

Power spectral density
h
7

1.5 6

1 1 4

0.5 2
 1
2

0 0
0 0.5 1 1.5 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Frequency/N -3 Frequency/N -3
x 10
x 10
(a) (b)
5 4
x 10 x 10
3 3
1,h
2.5 2.5 1 h
Power spectral density

Power spectral density

2 2

1.5 1.5

1 1

Next lowest frequencies

0.5 0.5

0 0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Frequency/N -3
x 10 Frequency/N -3
x 10

(c) (d)
Figure 10. Spectrum analysis for horizontal excited tank off and at resonance frequencies.

The waves generated by the vertical excitation are called Faraday waves as explored
originally by Faraday [34] through his experiments. Faraday waves are the resonant waves
when the excitation frequency is twice the natural frequency for some initial perturbation in
the container. This resonance condition is called parametric resonance. The study dealing
with vertical excitation of liquids in a container is referred to as parametric sloshing. For the
2
vertically excited tank, the parameter k v  Av v / g is a measure of the importance of the
vertical forcing motion and E is the measure of nonlinearity. Frandsen [25] plotted the
instability map between  v   n /  v and k v and discussed results from stability and
instability regions. If any of the pairs of the parameters lie in the instability region, then the
corresponding mode grows exponentially with time. In this section, the profiles are given
with stability and instability region. Figures 11 (a) and (b) are show behaviour of the liquid
free surface and phase plane diagrams respectively from the stability region (  v = 1.38 and
Kv = 0.4). The unstable regions results from  v = 1 and Kv = 0.4 for elevation and phase
plane are plotted in Figures 11 (c) and (d).

185
Finite Difference Based Sigma - Transformation Approach for Liquid Sloshing in a Rectangular Tank under Regular
Wave Excitation

0.4 0.3
At tank left wall
At tank middle
At tank right wall
0.3
0.2

0.2

( / t ) / ( Av x 1)
0.1

0.1
 / Av

-0.1
-0.1

-0.2
-0.2

-0.3 -0.3
0 20 40 60 80 -0.3 -0.2 -0.1 0 0.1 0.2 0.3
Time x 1  / Av

(a) (b)
0.8 0.8

0.6 0.6

0.4 0.4

( / t ) / ( Av x 1 )
0.2 0.2
 / Av

0 0

-0.2 -0.2

-0.4 -0.4

-0.6 -0.6
0 50 100 150 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6
Time x 1  / Av

(c) (d)
Figure 11. Stable and unstable regions from vertically excited tank elevation and phase
plane diagram for n=1 with grid size 41  61 and t  0.003 sec (a) - (b) Stable solution at
 v = 1.38 and Kv = 0.4; (c) - (d) Unstable region at  v = 1 and Kv = 0.4.

5.3.4. Combined excitation study (Sway and surge motion)

The combined motion of horizontal and vertical excitation is discussed in this

section. The vertical excited container is used for an ideal earthquake. And the combined
motions of horizontal and vertical excitations are two major considerable motions during
earthquake. The water waves at resonance may create the impact pressure rise inside the
tanks. When the liquid is striking on the wall, the impact pressure which might have been
raised near the critical pressure range of the tank material will cause the structural damage in
the tanks. The violent sloshing of combined excitation creates localized high impact loads
on the tank roof and walls which may damage the tank. During the combined excitation of
tank, the flow behaviour becomes tumultuous which creates the intricate free surface shapes
during this type of excitation. Due to the vertical excitation being present with this combined
motion, the instability regions exist as discussed in section 5.4. The initial wave impulse is
required for vertically excited condition and it is considered here as
 ( x, t )  0  A cos( K n x) . The vertical acceleration of tank is YD  (  v Av cos( v t )) and the
horizontal excited acceleration is fixed as X D  ( h Ah cos( h t )) . The initial velocity
potential in the fluid domain is considered as  ( , )   0  0 . The stable and unstable
regions from horizontally and vertically excited tank elevation and phase plane diagram are
presented Figure 12 (a) through (d). Figure 12 (a) and (b) are stable solutions from n=1 with

186
Eswaran & Saha CFD Letters Vol. 4(4) 2012

grid size 41  61, t  0.003 sec,  v = 1.38,  x =0.7 and Kv = 0.4 while Figure 12 (c) and
(d) are from unstable region at n=1 with grid size 41  61, t  0.003 sec,  v = 0.5,  x =0.7
and Kv = 0.4. As discussed, if any of the pairs of the parameters lie in the instability region,
then the corresponding mode grows exponentially with time as shown in Figure 12(c) and
(d).
0.03 0.025

0.02
0.02 0.015

0.01

( / t ) / ( Av x 1)
0.01
0.005
 / Av

0
0
-0.005

-0.01 -0.01

-0.015

-0.02 -0.02

-0.025
-0.03
0 20 40 60 80 100 -0.03
Time x 1 -0.02 -0.01 0 0.01 0.02
 / Av

(a) Stable solution at  v = 1.38, (b) Stable solution at  v =

 x =0.7 and Kv = 0.4. 1.38,  x =0.7 and Kv = 0.4.

60

40
 / Ah

20

-20

-40
0 10 20 30 40 50
Time x 1

(c) Unstable region at  v = 0.5, (d) Unstable region

at  v = 0.5,  x =0.7 and
 x =0.7 and Kv = 0.4.
Kv = 0.4.
Figure 12. Stable and unstable regions from horizontally and vertically excited tank
elevation and phase plane diagram for n=1 with grid size 41  61 and t  0.003
sec.

The stable and unstable region surface plots are shown in Figure 13 (a) and (b). The
surface plot is drawn between time and tank width and free surface elevation which is
showed in Figure 13 (a) from the stable region (n=1 with grid size 41  61, t  0.003 sec,
 v = 1.38,  x =0.7 and Kv = 0.4) and 14(b) from unstable region (n=1 with grid size 41 
61, t  0.003 sec,  v = 0.5,  x =0.7 and Kv = 0.4). Figure 14 shows the free surface
elevation from time 2.1 sec to 5.4 sec with the time interval of 0.3 seconds. The free surface
elevation along the tank width is illustrated here. It is observed that the waves are moving up
and down, and it must not be a uniform wave as observed in a fixed tank. One can easily
find the considerable variation due to the combined excitation.

187
Finite Difference Based Sigma - Transformation Approach for Liquid Sloshing in a Rectangular Tank under Regular
Wave Excitation

(a) Stable region  v = 1.38 and Kv = 0.4

(b) Unstable region  v = 0.5 and Kv = 0.4

Figure 13. Stable and unstable regions from horizontally and vertically excited tank surface
plot for n=1 with grid size 41  61 and t  0.003 sec at  x =0.7.

188
 Eswaran & Saha CFD Letters Vol. 4(4) 2012
2

-1

-2
0 0.5 1 1.5 2
V1

(a) Time =2.1 sec (b) Time =2.4 sec (c) Time =2.7 sec (d) Time =3.0 sec

(e) Time =3.3 sec (f) Time =3.6 sec (g) Time =3.9 sec (h) Time =4.2 sec

(i) Time =4.5 sec (j) Time =4.8 sec (k) Time =5.1 sec (l) Time =5.4 sec
Figure 14. Free surface profile diagram for horizontally and vertically excited tank n=1 with
grid size 41  61, t  0.003 sec,  v = 1.38,  x =0.7 and Kv = 0.4.

6. Summary and Conclusions

This paper deals with the non-linear effects of standing wave motion of liquid in a 2-D
rectangular tank. Governing equation and boundary conditions have been developed based on the
potential flow theory. A fully non-linear finite difference model had been developed based on the
inviscid flow equations, and a simple mapping function is used to remove the time-dependence of
the free surface in the fluid domain. Mapped governing equations and boundary conditions are
solved by finite difference method. Results of liquid sloshing induced by horizontal, vertical and
combined (horizontal and vertical) base excitations have been presented for small to steep non-
breaking waves. Simulations are limited to water fill ratio (hs/b) of 0.5. However, this numerical
model is valid for any water depth except shallow and deep sloshing. For the shallow sloshing,
viscous effects would become important, which we considered as inviscid fluid in our numerical
model. Primarily, the model was validated for sloshing motions and the effect of the steepness of
wave for a fixed tank. The  -transformation is limited to small steep non overturning waves. The
numerical model captured the free surface displacement at the left wall, right wall and at the center
of the tank. A good agreement between our numerical model and previously published result [25]
had been obtained for low steeping waves. The grid independent test was conducted for grid sizes
of 21  21, 41  21, 41  41 and 41  61 and finally, 41  61 is selected as grid resolution. Increasing
grid points in the vertical direction was found to be more effective in improving the accuracy. The
wave profiles along the tank at three different times for the first sloshing mode were plotted. Phase-
plane plot shows the behaviour of the free surface with repeatable patterns for the peaks and troughs
in bounded orbits.

189
 Finite Difference Based Sigma - Transformation Approach for Liquid Sloshing in a Rectangular Tank under Regular
Wave Excitation

The numerical wave tank captured steep waves generated by horizontal, vertical and
combined forcing amplitudes. The horizontally excited tank free surface elevation and phase plane
diagram was discussed with off and at resonance frequencies. Behaviours show the effect of wave
excitation while matching with resonance frequency. The spectrum analysis of horizontally excited
tank is also presented. The vertical and excitation causes the instability associated with parametric
resonance of the combined motion for a certain set of frequencies and amplitude of the vertical
motion. The initial condition for the surface elevation was an important parameter as there should
be some initial perturbation in the system for the generation of waves due to vertical excitation.
These conditions were also applied to combined excitation. Here, the data were chosen from the
stability region as well as instability region as discussed in section 5.4 and 5.5. Early simulations of
the liquid sloshing problem have mostly been performed with waves of small steepness. But the
present work reported with low and moderate steepness of wave. The present work can be extended
to a 3D tank and can also be solved by some other higher order numerical methods like compact
scheme. Moreover,  -transformation can be applicable for polar coordinates geometries as well.

Nomenclature

A Wave amplitude, m
Av Vertical forcing amplitude, m
Ax Horizontal forcing amplitude, m
b Length of the tank, m
E Wave steepness
hs Still water depth, m
h Instant water height from tank bottom, m
Kn Wave number
Kv Nonlinearity parameter in vertical direction
Kh Nonlinearity parameter in horizontal direction
n Mode number
YD Vertical acceleration of the tank, m/s2
X D Horizontal acceleration of the tank, m/s2
Greek symbols
n Natural sloshing frequency , rad/s
v Frequency of vertical motion, rad/s
h Frequency of horizontal motion, rad/s
 Free-surface elevation, m
 Stretching Factor
 Velocity potential function at physical domain ( x, y , t )
 Velocity potential function at transformed domain ( X ,  , T )
 Velocity potential function at Computational domain ( , , )
v Frequency ratio in vertical direction (   n /  v )
x Frequency ratio in vertical direction   h  n 
 , 2 First and second order central difference operators
  ,  Mixed second order central difference operator.

190
 Eswaran & Saha CFD Letters Vol. 4(4) 2012

References

[1] Maleki, A., and Ziyaeifar, M., Sloshing damping in cylindrical liquid storage tanks with
baffles. Journal of Sound and Vibration, 2008. 311: p. 372–385.
[2] Ibrahim, R.A., Liquid Sloshing Dynamics: Theory and Applications. 1 st edition, Cambridge
University Press, 2005, New York.
[3] Jacobsen, L.S., and Ayre, R.S., Hydrodynamic experiments with rigid cylindrical tanks
subjected to transient motions. Bulletin of the Seismological Society of America, 1951. 41:
p. 313-346.
[4] Graham, E.W., and Rodriguez, A.M., Characteristics of fuel motion which affect airplane
dynamics, Journal of Applied Mechanics, 1952. 19: p. 381–388.
[5] Housner, G.W., Dynamic pressures on accelerated fluid containers. Bulletin of the
Seismological Society of America, 1957. 47: p. 15-35.
[6] Housner, G.W., The dynamic behavior of water tanks. Bulletin of the Seismological Society
of America, 1963. 53: p. 381-387.
[7] Abramson, H.N., The dynamics of liquids in moving containers. Report SP 106, 1966,
NASA.
[8] Cole, H.A., Baffle thickness effects in fuel sloshing experiments. TN D-3716, 1966. NASA.
[9] Eswaran, M., Singh, A., and Saha, U.K., Experimental measurement of the surface velocity
field in an externally induced sloshing tank, Proceedings of the Institution of Mechanical
Engineers, Part M, Journal of Engineering for the Maritime Environment, 2011. 225(2): p.
133-148.
[10] Modi, V.J., and Seto, M.L., Suppression of flow-induced oscillations using sloshing liquid
dampers: analysis and experiments. Journal of Wind Engineering and Industrial
Aerodynamics, 1997. 67 & 68: p. 611-625.
[11] Frandsen, J.B., Borthwick, A.G.L., Simulation of sloshing motions in fixed and vertically
excited containers using a 2-D inviscid  -transformed finite difference solver. Journal of
Fluids and Structures, 2003, 18: p. 197-214.
[12] Godderidge, B., Turnock, S., Tan, M., Earl, C., An investigation of multiphase CFD
modeling of a lateral sloshing tank. Computers and Fluids, 38: p. 183–193.
[13] Cho, J.R., Lee, H.W., Numerical study on liquid sloshing in baffled tank by nonlinear finite
element method. Computer Methods in Applied Mechanics and Engineering, 2004. p. 193:
2581–2598.
[14] Behr, M., Abraham, F., Free-surface flow simulations in the presence of inclined walls.
Computer Methods in Applied Mechanics and Engineering, 2002. 191: p. 5467–5483.
[15] Tezduyar, T.E., Schwaab, M., and Sathe, S., Sequentially-coupled arterial fluid–structure
interaction (SCAFSI) technique. Computer Methods in Applied Mechanics and Engineering,
2009. 198:p. 3524–3533.
[16] Popov, G., Sankar, S., Sankar, T.S., and Vatistas, G.H., Dynamics of liquid sloshing in
horizontal cylindrical road containers. Proceedings of the Institution of Mechanical
Engineers, Part C: Journal of Mechanical Engineering Science, 1993. 207-C6: p. 399-406.
[17] Akyildiz, H., and Unal, N., Sloshing in a three-dimensional rectangular tank: numerical
simulation and experimental validation. Ocean Engineering, 2006. 33: p. 2135–2149.
[18] Cho, J.R., Lee, H.W., and Ha, S.Y., Finite element analysis of resonant sloshing response in
2D baffled tank. Journal of Sound Vibration, 2005. 288: p. 829–845.
[19] Arafa, M., Finite element analysis of sloshing in rectangular liquid-filled tanks. Journal of
Vibration Control, 2007. 13: p. 883–903.
[20] Hirt, C.W., and Nichols, B.D., Volume of fluid (VOF) method for the dynamics of free
boundaries. Journal Computational Physics, 1981. 39: p. 201–205.

191
 Finite Difference Based Sigma - Transformation Approach for Liquid Sloshing in a Rectangular Tank under Regular
Wave Excitation

[21] Eswaran, M., Saha, U.K., and Maity, D., Effect of baffles on a partially filled cubic tank:
Numerical simulation and experimental validation. Computers and Structures, 2009, 87(3-
4): p. 198–205.
[22] Cruchaga, M., Celentano, D., Breitkopf, P., Villon, P., and Rassineux, A., A surface
remeshing technique for a lagrangian description of 3D two-fluid flow problems.
International Journal of Numerical Methods in Fluids, 2010. 63: p. 415-430.
[23] Chern, M.J., Borthwick, A.G.L., and Taylor, R.E., A pseudospectral  -transformation
model of 2-D nonlinear waves. Journal of Fluids and Structures, 1999. 13: p. 607- 630.
[24] Turnbull, M.S., Borthwick, A.G.L., and Taylor, R.E., Numerical wave tank based on a  -
transformed finite element inviscid flow solver. International Journal for Numerical Methods
in Fluids, 2003, 42: p. 641–663.
[25] Frandsen, J.B., Sloshing in excited tanks. Journal of Computational Physics, 2004. 196: p.
53-87.
[26] Chen, B.F., and Nokes, R., Time-independent finite difference analysis of fully non-linear
and viscous fluid sloshing in a rectangular tank. Journal of Computational Physics, 2005.
209: p. 47-81.
[27] Dai, L., and Xu, L., A numerical scheme for dynamic liquid sloshing in horizontal
cylindrical containers. Proceedings of the Institution of Mechanical Engineers, Part D,
Journal of Automobile Engineering, 2006. 20: p. 901-918.
[28] Eswaran, M., and Saha, U.K., Low steeping waves simulation in a vertical excited container
using -transformation. Paper No. OMAE2009–80248, 28th International Conference on
Ocean, Offshore and Arctic Engineering, 2009, Honolulu, Hawaii, USA.
[29] Eswaran, M., and Saha, U.K., 'Wave simulation in an excited cylindrical tank using sigma
transformation’, Paper No. IMECE2010-39752, ASME International Mechanical
Engineering Congress & Exposition, 2010, Vancouver, Canada.
[30] Phillips, N.A., A coordinate system having some special advantages for numerical
forecasting. Journal of Atmospheric Sciences, 1957, 14: p.184–185.
[31] Blumberg, A.F., and Mellor, G.L., A coastal ocean numerical model, Mathematical
Modeling of Estuarine Physics, 1980. P. 203-219.
[32] Mellor, G.L., and Blumberg, A.F., Modeling vertical and horizontal diffusivities with the
sigma coordinate system. Monthly Weather Review, 1985. 113: p. 1379-1383.
[33] Faltinsen, O.M., A nonlinear theory of sloshing in rectangular tanks. Journal of Ship
Research, 1974, 18: p. 224-241.
[34] Faraday, M., On a peculiar class of acoustical figures; and on certain forms assumed by
groups of particles upon vibration elastic surface. Philosophical Transactions of the Royal
Society, 1831. London 121: p. 299–340.

192