A Finite-Element Computation For The Sloshing Motion in LNG Tank
A Finite-Element Computation For The Sloshing Motion in LNG Tank
A Finite-Element Computation For The Sloshing Motion in LNG Tank
Proceedings of The Twelfth (2002) International Offshore and Polar Engineering Conference
Kitakyushu, Japan, May 2631, 2002
Copyright 2002 by The International Society of Offshore and Polar Engineers
ISBN 1-880653-58-3 (Set); ISSN 1098-6189 (Set)
2
Seoul National University
San 56-1, Shinlim Dong, Kwanak Ku, Seoul 151-742, Korea
tower is a jacket-type structure that supports discharge Mikelis (1984), Arai (1994), and Kim (2000). More
and loading pipes. The large surge motion of the liquid at thorough review on the recent theoretical and numerical
the low filling level will induce significant drag force work based on the finite-difference scheme can be found
concentrated at the lower part of the tower. Mild impact in Faltinsen & Rognebakke (2000) and Faltinsen &
due to the steep wave front is also anticipated. Timokha (2001). Another popular numerical method is
Reinforcement of the lower part and the mounting system finite-element method. To name a few, Nakayama and
may be necessary to take care of the sloshing load on the Washizu (1990), Okamoto (1990) and Wu et al. (1998)
lower part of the pump tower. For the assessment of the applied the finite-element method to sloshing problems.
all these structural issues, correct estimation of the
hydrodynamic pressure and the flow field is very The theoretical approach has not shown successful
important. estimation of the impact pressure at the low filling yet.
The earlier work of Verhagen, J.H.G, and van
Besides the hydrodynamic impact, sloshing in LNG tank Wijngaarden (1965), based on shallow-water wave
involves more complicated physical phenomena. To theory, provides a good estimation of the total sloshing
name a few, the effect of compressibility and phase force at the low filling level. The numerical methods
transition during the impact, air entrapment, cushioning based on the finite-difference method also have shown
effect due to corrugation, elastic effect of the tank wall, inaccuracies in the conservation of mass and numerical
etc. Also important is the dynamic response of the instability, which made the long-time simulation difficult
structure to the impact load. The dynamic response of the (Faltinsen & Timokha, 2001). Some numerical method
structural system is known to be sensitive to the ratio also could not capture the impact pressure at the low
between the duration of the impulse and the longest filling level because of their artificial numerical damping
natural period of the structure (Abramson et al., 1974). and/or compressibility. The numerical damping usually
makes the wave front smoother than the reality, which
Due to this complexity, analysis of the sloshing load results in the loss or weakening of the impact pressure.
relies on model tests. Experiment has usually been done The artificial compressibility also introduces numerical
with the rigid container filled with water, where the error in pressure and mass conservation when severe
physical phenomena listed above are ignored. To account impact occurs.
for the neglected phenomena, a simple calibration has
been done by comparing the measured peak pressure in In this study, we have developed a finite-element code,
the model test and the estimated collapse pressure from SLOFE, to analyze the sloshing pressure in two
the previous damage in the LNG carriers (Tanaka et al., dimensions. Emphasis is made on the correct simulation
1984). The calibrated pressure is much lower than the of the impact pressure at the lower filling level, where
pressure from similarity law based on the Froude and the large impact load due to the hydraulic jump has been
Euler scaling. However, the calibration should be done observed. The existing finite-difference codes based on
very carefully when one has to consider a new SOLA-SURF and SOLA-VOF scheme have not been so
operational condition. It is doubtful whether one can use successful in predicting this impact pressure at the lower
the same calibration factor that has been used at the high filling level. The new finite-element method presented
filling conditions for the sloshing pressure at the low here is based on variational principle in stream-function
filling level. The answer should be made after formulation. This method satisfies the continuity
investigating the difference between the impact pressure equation and the kinematic boundary conditions exactly.
at these different filling levels. Previous theoretical and As a result, the conservation of mass and mechanical
experimental study showed that the flow pattern at the energy is exactly satisfied before we introduce numerical
high filling level is that of the standing wave, which is damping to consider the viscous damping empirically.
quite different from the hydraulic jump at the low filling We also developed a unique scheme to obtain the
level case. It is important to find the difference in the impulsive part of the pressure separately without any
sloshing pressure pattern at these two different situations. numerical differentiation that can introduce numerical
error for short-duration impact pressure. As a result, we
Many theoretical and numerical methods have been could obtain accurate description of the impact pressure
proposed to analyze the flow and pressure pattern in in spatial and temporal space.
sloshing problem. Faltinsen (1978), Faltinsen &
Rognebakke (2000), Faltinsen et al (2000) and Faltinsen A couple of numerical technique is introduced to handle
& Timokha (2001) used the semi-analytic analysis based the presence of tank top and to allow the dry bottom at
on nonlinear modal analysis and perturbation method. the tank bottom when sloshing motion is severe. The
The finite-difference method has been one of the most kinematic condition on the tank top, which can be stated
popular numerical method for sloshing problem since the that the fluid particle cannot not move above the tank
development of SOLA, VOF and SURF schemes (see top, is replaced by a dynamic condition by introducing
Chan and Street, 1970 and Hirt et al, 1975). Recent penalty pressure to suppress the fluid particle from
practical application of these scheme can be found in moving above the top. The penalty pressure is also
exerted when the thickness of the liquid layer is smaller Then the governing equation in the fluid domain is given
than a threshold thickness to prevent the numerical by
scheme from being unstable due to the dry bottom or the
zero thickness of the liquid layer. 2 = 2 in D . (2)
To validate the present numerical scheme, the model test On the tank wall, SW, no mass flux is allowed
performed by Abramson et al. (1974) is simulated. The
impact pressure at the knuckle points at the upper and
= 0 on SW (3)
lower chamfer agreed well with the model test. The
spatial and temporal pattern of the impact pressure at the
four different filling levels is compared. To measure the On the free surface, SF, the evolution of two canonical
effective area of the sloshing pressure, we introduced the variables (Miles, 1977), the free surface (x, t ) and the
idea of panel pressure. It has been found that the impact free-surface potential (x, t ) , is governed by
pressure at the low filling level has more effective area
than the one at the high filling level. The consequence of
this difference in the high- and low filling level is also = on S F (4)
discussed. t x
)( )
MATHEMATICAL FORMULATION 1
+ 2x 2x + 2 x x x
( 2
t 2 1 + x
(5)
2 pf
y xx+gx+ Ax = on S F .
x 2
2 1 2
u= , v= (1) p(x, y, t ) pi (x, y, t ) g x + A x x x +
y x 2 2
(7)
Then the impulsive pressure p i (x, y, t ) is governed by coefficients and will also be determined empirically.
2 1 2
Here bottom * = h + and will be determined
p i (x, y, t ) = p f + g x + A x x x + on S F
2 2 empirically.
(10) FINITE-ELEMENT METHOD
The governing equations for the impact pressure, Eqs. Variational Formulation
(7)~(10), can easily be derived by integrating the Euler
equations of the liquid. The initial/boundary value problem given in the previous
section can be equivalently formulated as a variational
MODELING OF TANK TOP AND BOTTOM problem to make the following functional stationary.
We assumed that the tank surface is rigid. When the
surface elevation is high enough, it can hit the tank top
L = +
2
(
+ 4 + 2 r dxdy
t x 2
2
)
but it cannot go over the tank top. The constraint on the
tank top can be given as the following inequality: (
g + Ay 2 )
Ax x + x +
2 x
top (11)
(14)
The above kinematic condition is not a convenient form
for numerical modeling. Instead, we introduce an Here, Ax and Ay denote the x and y component of the tank
equivalent dynamic condition. When the free surface hit acceleration A. The variational formulation given here
the ceiling, an impact pressure is exerted on the fluid can be obtained from the on given in Kim et al. (2001),
surface. This pressure can be interpreted as the constraint by generalizing the formulation to the rotational
force to satisfy the inequality (11). As a result, we can coordinate system.
satisfy the inequality by introducing a additional pressure
(or Lagrange multiplier) on the free surface. This Finite-Element Discretization
approach, however, also introduces numerical difficulty
with additional unknowns. One of the practical ways to The stream function (x, y, t ) is interpolated by a
enforce this kind of constraint in structural analysis is the bilinear finite-element interpolation function. The surface
use of contact spring, which is a stiff spring only elevation (x, t ) and the velocity potential (x, t ) are
triggered on when the free surface is close to the tank
interpolated by linear interpolation functions:
wall. The rigid ceiling is replaced by a spring-damper
system, with sufficiently high stiffness to allow Ne
deformation within an admissible range. We also
introduce damping to suppress the short numerical waves
( x, y , t ) = (t )N (x, y )
i =1
i i (15)
due to the high stiffness of the ceiling. Then the
additional penalty pressure is given by Nx
(x, t ) = i (t )M i (x ) (16)
( *
top +
pf =
)
t
, > top *
(12)
i =1
0, < top * Nx
(x, t ) = i (t )M i (x ) (17)
i =1
number of nodes. Substituting the Eqs. (15)~(17) into Eq. Table 1 Input parameters for sloshing simulations for
(14), and taking variations with respect to the variables various filling levels
i (t ), i (t ) and i (t ) , we can obtain ordinary differential
equations for i (t ), i (t ) and algebraic equation for Case A Case B Case C Case D
i (t ) . The Runge-Kutta 4th order method is used to
d/b 0.12 0.20 0.40 0.71
integrate the ordinary differential equations. The details
of the procedure are quite similar to the one given in Kim d/h 0.16 0.27 0.54 0.96
& Bai (2000).
A/b 0.10 0.10 0.10 0.10
NUMERICAL RESULTS
/ g 82 82 82 82
A two-dimensional prismatic tank that has been used by
Abramson et al. (1974) for their model test is used for the
numerical simulation. The sketch of the geometry is / gh 2.1 2.1 2.1 2.1
given in Fig. 2. The dimensions of the tank are given by
b = 1.38m, h = 1.02m, hu = 0.295m, hl = 0.14m / gb2 3.14 x 10-4 3.14 x 10-4 3.14 x 10-4 2.00 x 10-5
and u = b = 135o.
/h 2.35 x 10-4 2.35 x 10-4 2.35 x 10-4 2.35 x 10-4
P1
hl b Fig. 3 shows the comparison of the impact pressure with
model test results for Case A. Abramson et al. (1974)
Fig. 2 Definition sketch of a two-dimensional LNG tank reported that the peak values of the impact pressure are
random. The experimental values in Fig. 3 are the 10%
exceedance limit of the peak values they provided.
The tank is excited by a regular harmonic sway motion:
However, the calculated pressure shows periodic pattern
after a few transient periods, as can be seen in Fig. 4.
2t
r (t ) = A sin ,0 , This is presumably due to the three-dimensional effect in
T the experiment. The impact pressure at the knuckle point
in intermediate filling level (Case C) is given in Fig. 5.
where A is the sway amplitude and T is the period of the The maximum pressure in the SLOFE result occurred at
motion. slightly shorter excitation period than the model test. The
theoretical result of Faltinsen and Rognebakke (2000)
Numerical simulation has been performed for the four shows significantly lower values than the experiment and
filling levels given in Table 1. The sway amplitude is SLOFE result. Fig. 6 shows the time history of impact
given as 10% of the breadth of the tank, which pressure when T = 1.6 sec. In general, the present
corresponds to a realistic extreme sway motion of a ship numerical scheme predicts the magnitude of the impact
(Abramson et al., 1974). The other input parameters used pressure quite well.
in the numerical simulation are also given in Table 1.
The parameters , and are determined by calibrating
the impact pressure from the model tests.
4
2
3
p/gb
p/gb
2
1
0
0
0 2 4 6 8 10
0 2 4 6 8 10
t/T
Fig. 4 Time history of impact pressure at P1. t/T
d/b = 0.12, A/b = 0.1, T = 1.85 sec. Fig. 6 Time history of impact pressure at P2.
d/b = 0.12, A/b = 0.1, T = 1.6 sec.
Fig. 8 Pressure distribution in the fluid domain and on the tank Fig.10 Pressure distribution in the fluid domain and on the tank
walls for Case B. Pmax denotes the maximum value of p/gb on walls for Case D. Pmax denotes the maximum value of p/gb on
the tank wall. the tank wall
Fig. 9 Pressure distribution in the fluid domain and on the tank We compare two pressures to show the spatial
walls for Case C. Pmax denotes the maximum value of p/gb on distribution of the impact load. Fig. 12 shows the time
the tank wall history of the maximum and average panel pressure at
the four different filling level. As we expected, the ratio
the high filling level. The ratio is higher than 0.5 for Case 6
B, where the impact pressure on the bulkhead is most
p/gb
significant. At high filling level, Case D, the ratio is 4
about 0.2. This is because the peak of the impact load is
highly localized at this filling level as can be seen in Fig.
2
10. This clearly shows that the impact pressure at the low
filling level has wider effective area than the one at the
high filling level. 0
0 2 4 6 8 10
As mentioned in the introduction, the calibration of the t/T
pressure from the model test or simulation has been made (a) Case A, T = 1.85 sec
based on collapse pressure of the insulation box at the
high filling level. The collapse of the insulation box 8 Max. panel pressure
Avg. panel pressure
should depend on the total load on a finite area of
insulation system. The calibration factor at the high 6
filling level may not be valid for the low filling level case
p/gb
where the impact load has wider effective area as shown 4
in the present numerical simulation. In addition, the angle
between the wave front of and the tank wall was usually
2
greater than 15 degree at the low filling level (see Fig. 7
and 8). This indicates that the air entrapment is less likely
0
to happen at the low filling level, which may provide
another reason for the calibration factor to be different 0 2 4 6 8 10
t/T
from the one at high filling level.
(b) Case B, T = 1.6 sec
8
Max. panel pressure
Avg. panel pressure
6
p/gb
0 2 4 6 8 10
t/T
Fig. 11 Panels to evaluate panel pressure
(c) Case C, T = 1.6 sec
8
Max. panel pressure
Avg. panel pressure
6
p/gb
0 2 4 6 8 10
t/T
(d) Case D, T = 0.87 sec
Chan, K.C. and Street, R.L. 1970 SUMMAC - A Wu, G. X., Ma, Q. W. AND Taylor, R. E 1998
numerical model for water waves. Stanford University, Numerical simulation of sloshing waves in a 3D tank
Report 135 based on a finite element method. Applied Ocean
Research, vol. 20
Faltinsen, O.M. 1990 Sea Loads on Ships and Offshore
Structures, Cambridge University Press. Zhao, R., and Faltinsen, O.M. 1993 Water Entry of Two-
Dimensional Bodies, J. Fluid Mech., Vol. 246, pp. 593-
Faltinsen, O. M. 1978 A numerical non-linear method of 612.
sloshing in tanks with two-dimensional flow: Journal of
Ship Research, vol. 18.