- Email: [email protected]
Numerical study on liquid sloshing in baffled tank by nonlinear finite element method
Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598 www.elsevier.com/locate/cma
Numerical study on liquid sloshing in baffled tank by nonlinear finite element method J.R. Cho *, H.W. Lee C&AM Lab, School of Mechanical Engineering, Pusan National University, 30 Jangjeon-Dong, Kumjung-Ku, Pusan 609-735, South Korea Received 2 June 2003; received in revised form 4 December 2003; accepted 8 January 2004
Abstract This paper introduces a velocity–potential-based nonlinear finite element method for the accurate simulation of the large amplitude liquid sloshing in two-dimensional baffled tank subject to horizontal forced excitation. As well, the effects of baffle on the nonlinear liquid sloshing are parametrically examined. The free surface configuration is tracked by a direct time differentiation of the convective term in the kinematic boundary condition with the help of the four-step predictor–corrector methods. The flow velocity is interpolated from the velocity potential by a second-order least square method. And, the finite element mesh is adapted in a semi-Lagrangian approach combined with the free-surface height correction, in order to keep the mesh regularity and the total liquid volume unchanged. By comparing with the available reference solutions, we verify the numerical accuracy and stability of the nonlinear finite element method proposed. Through the numerical experiments performed by varying the installation height and the opening width of baffles, the hydrodynamic characteristics of the large-amplitude liquid sloshing are parametrically investigated. 2004 Elsevier B.V. All rights reserved. Keywords: Baffled tank; Nonlinear liquid sloshing; Direct time differentiation; Semi-Lagrangian free surface tracking; Baffle installation height; Opening width
1. Introduction Free surface fluctuation called liquid sloshing is the most prominent phenomenon of liquid motion in either stationary or moving tanks subject to forced external perturbations. Representative examples can be found from liquefied natural gas (LNG) carriers in shipping service and aboveground liquid-storage containers excited by earthquake wave [1]. When the liquid sloshing occurs, hydrodynamic forces acting on the container display the severe fluctuation in their magnitudes and distributions. Moreover, an excessive liquid sloshing may cause the structural failure or/and the manipulation loss, and which frequently leads to the tremendous loss of human, economic and environmental resources.
*
Corresponding author. Tel.: +82-51-510-2467; fax: +82-51-514-7640. E-mail address: [email protected] (J.R. Cho).
0045-7825/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2004.01.009
2582
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
In view of this critical effect of the liquid sloshing on liquid–structure interaction systems, very intensive research efforts have been made to suppress it from several decades ago [2,3], as a result, a number of passive devices for suppressing the liquid motion have been devised. Among them, the flexible baffle manufactured usually with metal is being widely used thanks to the practical convenience and the flexible adjustment in its installation to the container. However, the suppression efficiency of baffle is strongly affected by several design parameters, particularly the installation height and the opening width of baffle, as demonstrated through our previous works related [4,5]. Moreover, the suitable values of baffle parameters are dependent of the problem under consideration, so that the determination of appropriate baffle parameters needs to be made based upon the sufficient parametric investigation on the liquid sloshing in baffled tank. As is widely known, liquid sloshing, as a representative example of fluid–structure interactions, is an intrinsic nonlinear problem in which the free surface configuration is not known a priori. This nonlinearity becomes higher when baffles are introduced as well as when the excitation frequency is close to the natural sloshing frequencies or the excitation amplitude is large. Nevertheless, to the best of our literature survey, the early studies before 1970s had relied on the linear theory with the assumption of small amplitude sloshing [6,7]. However, after mid 1970s, active studies for developing the numerical techniques have been made to successfully implement the nonlinear fluid–structure interaction problems [8–15]. The main concern in the numerical implementation of fluid–structure interactions exhibiting a remarkable problem nonlinearity has focused on the insurance of the numerical stability and accuracy. Because the excessive distortion of fluid elements causes the numerical instability or even the computation incompleteness [16], and a small error in the time integration accumulates with time stage. Meanwhile, the choice of state variables (i.e. the underlying governing equations) and the continuity constraint enforcement for incompressible flows [15,17] characterize the numerical formulations and approximations. Besides, several another issues such as the coupling between solid and fluid have been also intensively studied for the reliable numerical simulation. The mesh distortion and the mesh adaptation are strongly related to the choice of kinematic description approach (e.g. Lagrangian, Eulerian, or ALE (arbitrary Lagrangian–Eulerian) [9]). Currently, this problem can be resolved, to a large extent, by employing the ALE method with the help of a suitable remeshing and smoothing algorithm [18]. Regarding this matter, the reader may refer to [4,9,12,13,16,17] for more detailed discussion. On the other hand, the fluid flow has been traditionally formulated by either the potential flow theory [5,8,11,14,19] or the full Navier–Stokes equations. In the former case, the free surface configuration is tracked by time-integrating the nonlinear kinematic and dynamic conditions, and the flow velocity and the dynamic pressure are usually interpolated from the velocity potential approximated. However, in the latter case in which both flow velocity and dynamic pressure are the primary state variables, the flow boundary identification varies to the choice of kinematic description approach. Referring to the papers by Soulaimani and Saad [20] and Cho and Lee [4], the boundary tracking in both Lagrangian and ALE approaches is straightforward because fluid mesh moves exactly with fluid particles. But, in the Eulerian approach the boundary tracking requires any special characteristic (volume fraction) function [21] with the mesh refinement along the flow boundary. The fluid mass conservation in incompressible flows is enforced in several ways. Among them are a penalty method or an iterative augmented Lagrangian method (see Oden et al. [15] and Childs and Reddy [22]), a correction of the free surface height [4,23], and an iterative split algorithm combined with the fractional step method [24]. Regarding the mesh interface between solid and fluid, which is beyond the scope of the current work, the reader may refer to the papers by Childs and Reddy [22] and Cruchaga et al. [25]. In the present study, we intend to introduce a nonlinear finite element method based on the fully nonlinear potential flow theory which accurately simulate the large amplitude sloshing flow in twodimensional baffled tank subject to horizontal forced excitation. In particular, the main concern focuses
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
2583
on the insurance of the numerical stability and accuracy by suppressing the numerical dissipation and the total liquid volume change, and the parametric investigation of the baffle effects on the nonlinear liquid sloshing. This kind of problem was researched by many investigators such as Ikegawa [8] and Faltinsen [11]. Restricting to the free surface tracking in the velocity-potential-based nonlinear sloshing analysis, Nakayama and Washizu [14] proposed a nonlinear boundary element method by introducing an error correction term into the dynamic boundary condition, in order to suppress the numerical instability and dissipation. Okamoto and Kawahara [26] introduced a velocity correction method into the Lagrangian finite element formulation, and they justified their numerical method by comparing the numerical solutions with independent experimental results. Chen et al. [27] employed the finite difference scheme based on the Crank–Nicolson time marching scheme in which a second-order numerical dissipation term was added to the kinematic boundary condition. Wu et al. [28] employed the open trapezoidal scheme for the stable time integration of nonlinear liquid sloshing in three-dimensional tank. On the other hand, we in the present study track the free surface configuration by directly time differentiating the convective term in the kinematic boundary condition with the help of the four-step predictor– corrector methods with the accuracy of OðDt4 Þ. As well, we introduce a semi-Lagrangian mesh adaptation and the liquid height correction algorithm, while assuming the container be a rigid body, in order to keep the mesh regularity and the total liquid volume unchanged. Furthermore, the flow velocity is interpolated from the velocity potential by the second-order least square method, and the dynamic pressure field is approximated with the same finite element framework constructed for the velocity potential approximation. The numerical stability and accuracy of the nonlinear finite element method introduced are verified through the comparison with the existing reference solutions. We also perform the parametric experiments by varying the installation height and the opening width of baffles, in order to examine the effects of baffle on the nonlinear liquid sloshing. Major quantities of interest are the maximum free-surface elevation, hydrodynamic force and moment resultants, and the variations of flow profile and hydrodynamic pressure. This paper is organized as follows. The problem formulation of nonlinear liquid sloshing in baffled tank is given in Section 2. A time-incremental nonlinear finite element approximation based on the direct time differentiation and the semi-Lagrangian approach and the evaluation of hydrodynamic forces and moments are presented respectively in Sections 3 and 4. Next, the verification and parametric numerical results are addressed in Section 5, and the final conclusion is made in Section 6.
2. Nonlinear liquid sloshing in baffled tank Referring to Fig. 1, we consider the sloshing flow of inviscid incompressible liquid in a two-dimensional baffled rigid tank of width d which is subjected to forced horizontal excitation. Where, a pair of rigid baffles with the opening width dB is installed at a height hB and liquid is filled up to a height hL in the stationary condition. Boundary of the liquid domain Xt 2 R2 is composed of the free surface oXtF and the liquid– structure interface oXtI , in which the superscript t indicates that both the liquid domain and the liquid boundary vary with time during the observation time interval ð0; t. For the sake of mathematical description, we use two Cartesian co-ordinate systems, fX0 g fixed in space and fX g moving with the tank, such that their origins are commonly positioned at the center of the stationary free surface and their axes are aligned parallel to each other. Then, a quantity Q can be defined in ~ y; tÞ. By denoting us ðtÞ ¼ fdxs =dt; dy s =dtg be either of two co-ordinate systems such that Qðx0 ; y0 ; tÞ ¼ Qðx; the rigid tank velocity, we have the relations between two systems: rQjfX0 g ¼ rQjfX g ;
ð1Þ
2584
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
Fig. 1. Baffled rigid tank under horizontal forced excitation.
oQ oQ ¼ us rQ: ot fX0 g ot fX g
ð2Þ
Hereafter, the subscript 0 refers to the quantities measured in the fixed Cartesian co-ordinates. When we denote /ðx0 ; y0 ; tÞ be the velocity potential, the continuity condition is described by the Laplace equation: r2 / ¼ 0;
in Xt :
ð3Þ
On the liquid–structure interface, the velocity potential should satisfy r/ n ¼ us n;
on oXtI
ð4Þ
with the outward unit vector n normal to the structure boundary. In the fixed Cartesian co-ordinates, the kinematic and dynamic conditions on the free surface are as follows: og0 o/ of0 o/ þ ¼ 0; ox0 ox0 oy0 ot
on oXtF ;
o/ 1 þ r/ r/ þ gf0 ¼ 0; ot 2
on oXtF ;
ð5Þ ð6Þ
in which f0 ðx0 ; tÞ denotes the free surface elevation at time t. Eq. (5) is derived from the kinematic condition: o/=oy0 ¼ Df0 =Dt on the free surface, while Eq. (6) is the Bernoulli equation when the atmospheric pressure is assumed to be zero on the free surface [29]. Since the free surface fluctuation is a relative motion to the container the moving co-ordinates is usually chosen. By denoting fðx; tÞ ¼ f0 ðx0 ; tÞ y s ðtÞ be the free surface elevation from the reference stationary level, together with Eqs. (1) and (2), above free surface conditions can be rewritten as of o/ of o/ þ ð7Þ usx usy ¼ 0; ot ox ox oy o/ 1 us r/ þ r/ r/ þ gðf þ y s Þ ¼ 0: ot 2
ð8Þ
In the moving co-ordinate system, the velocity potential /ðx; y; tÞ can be split into two parts, /h ðx; y; tÞ due to the internal sloshing flow and /p ðx; y; tÞ due to the tank motion, such that / ¼ /h þ /p ;
/p ¼ xusx þ yusy :
ð9Þ
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
2585
Substituting Eq. (9) into Eqs. (7) and (8), we respectively have of o/h of o/h þ ¼ 0; ot ox ox oy
ð10Þ
dusy 1 s 2 o/h 1 dus þ r/h r/h þ gf þ x x þ f ¼ ju j gy s : 2 2 ot dt dt
ð11Þ
It is worth to note that two terms in the RHS of Eq. (11) vanish because not only both are co-ordinateindependent quantities but also the velocity potential is a relative quantity in space. Then, the initial-boundary-value problem of nonlinear sloshing flow in the moving co-ordinate system is formulated as follows: r2 /h ¼ 0;
in Xt
ð12Þ
with initial conditions /h ¼ xusx ð0Þ yusy ð0Þ; f ¼ 0;
in X0 ;
ð13Þ
on oX0F ;
ð14Þ
and boundary conditions r/h n ¼ 0;
on oXtI ;
of o/ of o/h ¼ h þ ; ot ox ox oy
ð15Þ on oXtF ;
dusy o/h 1 dus ¼ r/h r/h gf x x f ; 2 ot dt dt
ð16Þ on oXtF :
ð17Þ
3. Nonlinear finite element approximation In order to solve the above nonlinear sloshing problem we divide the observation time interval t into a finite number of sub-intervals such that tn ¼ nDt (n ¼ 0; 1; 2; . . .). Starting from the initial solutions /0h and f0 , we successively solve the Laplace equation (12) defined in current flow domain Xn and boundary oXn to seek /nh with the boundary condition (15). From which we compute the flow velocity field vn and perform the free surface tracking. By time-integrating the kinematic and dynamic boundary conditions (16) and (17), together with the mass-conservative semi-Lagrangian remeshing, we identify fnþ1 , Xnþ1 and oXnþ1 . In this manner, the sloshing time response is to be numerically analyzed. 3.1. Interpolation of flow velocity field Referring to Fig. 2, we consider a time stage tn at which /nh and vn are to be sought with the geometry and field quantities determined at the previous time stage tn1 . According to the virtual work principle, we have the weak form of Eq. (12): find /nh such that Z rw r/nh dv ¼ 0 ð18Þ Xn
2586
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
Fig. 2. Flow status representation at time stage tn .
for every admissible velocity potential wðx; yÞ. The corresponding essential boundary condition on oXnF is /nh determined previously by the free surface tacking method which will be described later. With nine-node N n , we approximate the Lagrange-type finite element basis functions fNi gi¼1 and the nodal potential vector / h velocity potential n: /nh ¼ N / h
ð19Þ
Substituting Eq. (19) into Eq. (18) leads to n ¼ 0 K/
ð20Þ
with the matrix K defined by Z T K¼ ðrNÞ ðrNÞ dv:
ð21Þ
h
Xn
Once the velocity potential /nh is approximated, the intermediate velocity field ~vn is computed according to n: ~vn ¼ rN / h
ð22Þ
0
Owing to the C -continuity in the velocity potential field, the intermediate velocity field becomes discontinuous along the element interfaces. In order to secure the continuity in the velocity field, one may consider the linear interpolation of the intermediate velocity values at Gauss integration points. However, a small error in the simple velocity interpolation does not only affect the accuracy of the free surface tracking but also accumulate with the time stage. In order to suppress this crucial problem, we first interpolate elementwise velocity fields from ~vn with second-order polynomials according to the least square method. Let us denote ^vna (a ¼ x; y) be the element-wise velocity components to be interpolated: ^vna ¼ a1 þ a2 n þ a3 g þ a4 ng þ a5 n2 þ a6 g2 :
ð23Þ
And we define the element-wise error by Ea ¼
9 X
2
f^vna ðn‘ ; g‘ Þ ~vna ðn‘ ; g‘ Þg ;
a ¼ x; y;
ð24Þ
‘¼1
where ‘ stands for (3 · 3) Gauss integration points. Then, six coefficients involved in each velocity component interpolation (23) are to be determined from six simultaneous equations which are constructed from oEa ¼ 0; oak
k ¼ 1; 2; . . . ; 6:
ð25Þ
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
2587
After the element-wise velocity fields are interpolated, we compute the nodal velocities of all elements and take average of nodal values for the common nodes. By interpolating the nodal velocities vna with the same finite element basis functions used as for the velocity potential approximation (19), we finally have the global continuous velocity field: vna ¼ N vna ;
a ¼ x; y:
ð26Þ
3.2. Time marching of free surface Once both the velocity potential /nh and the flow velocity field vn are approximated, we integrate Eq. (16) to predict the free surface fnþ1 and Eq. (17) to determine the essential boundary condition /hnþ1 . As a result, the flow geometry and the boundary condition required for the next-step computation are to be established. A mass-conservative remeshing process follows the free surface tracking, while keeping the mesh regularity fairly, in order to update the liquid mesh. Here, the time integration of the kinematic boundary condition (16) may lead to the numerical divergence, the high frequency wiggle or the numerical dissipation, unless any special care is paid. This crucial problem is owing to the convection term in the RHS of Eq. (16), and its influence increases in proportional to the sloshing amplitude. In order to suppress the numerical instability, several kinds of correction techniques have been proposed by Nakayama and Washizu [14], Okamoto and Kawahara [26], and Chen et al. [27], as mentioned before. For the present study, however, we introduce a direct time differentiation scheme to track the free surface. Referring to Fig. 3, we predict the free surface configuration after the small time interval dt (dt ¼ aDt, 0 < a < 1) from the current time tn according to o/h o/h 0 0 xi ¼ xi þ dt ; yi ¼ yi þ dt ; ð27Þ ox oy ðxi ;yi Þ
ðxi ;yi Þ
where ðxi ; yi Þ denotes co-ordinates of node i located on the free surface at time tn . And, the elevation of node i during the time interval dt is calculated through f0i ¼ yj0 þ
0 yjþ1 yj0 ðxi x0j Þ: x0jþ1 x0j
ð28Þ
Then, we directly calculate the time derivative of the free surface elevation in Eq. (16) using the alternative relation given by of df f0 fn : ð29Þ ¼ ot n dt n dt
Fig. 3. Variation of the free surface during the small time interval dt.
2588
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
The parameter a influences the numerical dissipation and the incompressibility condition, so that its choice should be carefully made. As presented in Section 5.1, the former is getting improved as a becomes smaller, but the latter deteriorates either as a approaches zero or unity. In the current study, we choose a suitable a based upon its parametric effects on both the numerical dissipation and the liquid volume conservation. Next, we apply the predictor–corrector method composed of the explicit Adams–Bashforth four-step scheme [30] and the implicit Adams–Moulton three-step scheme [30] given respectively by ~fnþ1 ¼ fn þ Dt 55 of 59 of þ 37 of 9 of ; ð30Þ 24 ot ot ot ot n
Dt o~f n 9 ¼f þ 24 ot
n1
"
fnþ1
þ 19 nþ1
n2
of of of 5 þ ot n ot n1 ot n2
n3
# ð31Þ
in order to build up the free surface configuration oXFnþ1 at the next time stage tnþ1 . After that, we relocate the free surface nodes onto the new free surface by moving those only in the vertical direction, as depicted in Fig. 4(a). As well as, we calculate the total liquid area according to the trapezoidal rule Atot ¼
K X
Ai ;
Ai ¼ ðxiþ1 xi Þ
i¼1
fiþ1 þ fi : 2
ð32Þ
The mean free-surface elevation fmean can be calculated as follows: fmean ¼ Atot =d, and the nodal free surface heights finþ1 are to be corrected [4], referring to Fig. 4(b), such that fnþ1 ¼ fnþ1 fmean : i i
ð33Þ
Referring to Fig. 5, we next relocate all internal nodes such that the nodes located on the same vertical line should be in almost uniform spacing. In this manner, the free surface is tracked with the total liquid volume unchanged and the liquid mesh is adapted to keep the mesh regularity fairly. In the next step, we integrate Eq. (17) to specify the velocity potential /hnþ1 to the new free surface oXFnþ1 tracked. However, the time derivative of the free surface velocity potential in the semi-Lagrangian liquid mesh fXM g is related to one in the moving co-ordinate system fX g such that o/h o/h of o/h : ð34Þ ¼ ot oy ot ot fX g
fXM g
Thus, the time derivative of the free surface velocity potential with respect to the liquid mesh, at time stage tn , is rewritten as
Fig. 4. Volume-conservative free surface tracking: (a) relocation of free surface nodes; and (b) correction of the free surface height.
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
2589
Fig. 5. Liquid mesh adaptation by semi-Lagrangian remeshing scheme.
dusy o/h of o/nh 1 dusx n n n n r/ f ; ¼ r/ gf x h h ot n oy 2 ot n dt dt
on oXnF :
ð35Þ
After evaluating this time derivative, we apply the predictor–corrector methods (30) and (31) used for the free surface tracking once again, in order to determine the velocity potential /hnþ1 on the free surface oXnþ1 F . 4. Hydrodynamic forces and moment The total pressure pðx; y; tÞ in the liquid domain is expressed as p ¼ pd qgy;
ð36Þ
where q denotes the liquid density. Meanwhile, the dynamic pressure pd ðx; y; tÞ is computed according to dusy o/h 1 dusx þy pd ¼ q þ r/h r/h þ x : ð37Þ 2 ot dt dt In above equation, one can realize that an additional finite element approximation is required to obtain the time derivative of the velocity potential /h . Letting o/h =ot be v, we have the following Laplace equation from Eq. (12): r2 v ¼ 0;
in Xt
ð38Þ
with boundary conditions rv n ¼ 0;
on oXtI ;
ð39Þ
dusy 1 dus v ¼ r/h r/h gf x x f ; 2 dt dt
on oXtF :
ð40Þ
At any time stage tn , the essential boundary conditions (15) and (17) are to be time-integrated during the finite element approximation of the velocity potential /nh . So that, this boundary value problem can be easily solved with the same liquid mesh after the velocity potential approximation at time stage tn . Referring to Fig. 6, the force resultants acting on the baffled tank due to the liquid sloshing are computed through Z Z Fx ¼ p dy p dy; ð41Þ CR T
CL T
Z
Z
Fy ¼
CBT
p dx
Cþ B
p dx þ
Z p dx: C B
ð42Þ
2590
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
Fig. 6. Division of the liquid–structure interaction.
And the moment resultant about the origin O is obtained according to Z Z Z Z Z py dy py dy px dx px dx þ px dx: Mz ¼ CL T
CR T
Cþ B
CBT
ð43Þ
C B
Of course, the force and moment resultants are in function of time.
5. Numerical experiments According to the numerical formulae described so far we developed a test program coded in Fortran which is interfaced with MSC/Patran for visualization. The forced excitation of tank is taken as a horizontal harmonic movement xs ðtÞ ¼ a sin xt, but the amplitude and the frequency are variable depending on the simulation case. Referring to Fig. 1, the tank width, the liquid fill height, and the installation height and the opening width of baffles are also taken variable. Interior liquid is initially discretized with uniform mesh of nine-node quadratic elements and the time step size Dt is set by 5 · 103 s throughout the experiments. 5.1. Verification of numerical accuracy and stability We first perform the verification experiments with liquid-filled tanks without baffle, in order to validate the time-incremental nonlinear finite element method introduced and the test program developed. Two simulation cases are considered, as presented in Table 1, where case 1 was taken by Nakayama and Washizu [14] while case 2 by Chen et al. [27] and Okamoto and Kawahara [26]. In order to determine a suitable value of a, we first examine the effects of a on the numerical dissipation and the total liquid volume. This experiment was done with case 2 because the sloshing nonlinearity is higher in case 2 than case 1, as will be seen later. It is worthy noting that this preliminary simulation was performed without correcting the free-surface height according to Eq. (33) so that errors in both quantities accumulate with time. As shown in Fig. 7(a), the numerical dissipation is getting suppressed as a becomes smaller such that the time history responses of the free surface elevation between a ¼ 0:1 and 0.01 do not show a noticeable distinction. On the other hand, the total liquid volume, referring to Fig. 7(b), considTable 1 Simulation cases for the verification experiment Case 1 2
Geometry
Excitation
Tank width d (m)
Liquid height hL (m)
Amplitude a (m)
Frequency x (rad/s)
0.9 (16) 1.0 (24)
0.6 (12) 0.5 (12)
0.002 )0.0093
5.5 5.31123
Subdivision number to construct uniform liquid meshes.
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
2591
Fig. 7. The effects of a (case 2): (a) on the free surface elevation (at the right end x ¼ 0:5 m); (b) on the liquid volume conservation.
erably decreases when a is unity, and contrary it slightly increases when a 6 0:1. From these features of both quantities, we set a by 0.1 for the current study. Fig. 8 shows the comparison with the solution by Nakayama and Washizu [14] that has been frequently used as a reference solution in comparative numerical experiments. As mentioned before, they improved the numerical accuracy by introducing an error correction term into the dynamic boundary condition, and they justified the numerical accuracy by comparing with the FEM solution by Ikegawa [8]. The time history response of free surface elevation was measured at the right end x ¼ 0:45 m. In view of the fundamental sloshing frequency x0 ¼ 5:761 rad/s of case 1, the excitation frequency x is taken by x ¼ 0:9547x0 . One can observe that the numerical solution by the proposed method is in an excellent agreement with the reference solution. Comparative numerical results carried out with case 2 are given in Fig. 9 in which the free surface elevation is taken at the right end x ¼ 0:5 m. The sloshing response denoted by upwind scheme in Fig. 9(b) was obtained by our test program, but for which a five-point upwind scheme [31] was used to deal with the convection term in the kinematic boundary condition (16) instead of our direct differentiation scheme (29). The comparison with the Chen et al.Õs solution, for which the second-order numerical dissipation term was introduced into the kinematic boundary condition [27], shows that the present method provides the accurate sloshing response without causing any numerical instability. However, the case with the upwind scheme displays a crucial numerical dissipation so that the free surface elevation becomes considerably smaller after the peak elevation time, even though the numerical divergence or the high frequency wiggle is not observed. Fig. 10 represents the comparison of free surface configurations at six different time stages with the Okamoto and KawaharaÕs experiment. The time marching scheme introduced in the current study
Fig. 8. Comparison with Nakayama and Washizu [14]: time history response of the free surface elevation (case 1).
2592
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
Fig. 9. Time history response of the free surface elevation (case 2): (a) direct time differentiation scheme; (b) five-point upwind scheme.
Fig. 10. Comparison of free surface sloshing flows at different time steps (case 2).
accurately predicts the free surface fluctuation to the lapse of time that is in a good agreement with the experimental result up to the time of observation. Thus, the numerical accuracy and stability of the nonlinear finite element method introduced have been sufficiently verified. 5.2. Parametric experiments In order for the parametric analysis to the installation height and the opening width of baffles, we consider a tank with a pair of baffles. Referring to Fig. 1, the liquid height and the tank width are set by d ¼ 1:0 m and hL ¼ 1:0 m, and the liquid mesh is constructed with 16 · 24 uniform 9-node quadrilateral elements. Similar to the previous verification experiments, we apply the horizontal harmonic excitation with the amplitude a ¼ 0:005 m and the excitation frequency x ¼ 0:999x0 . Here, the fundamental sloshing frequency x0 varies with the two baffle parameters [5]. For the current parametric investigation, we consider the relative baffle heights hB =hL between 0.3 and 0.7 and the relative opening widths dB =d between 0.2 and 0.8, respectively, and we set the observation time duration by 40 s.
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
2593
Fig. 11 shows the variation of sloshing flow during one period for the case of hB =hL ¼ 0:7 and dB =d ¼ 0:5. We can observe two major features in the liquid sloshing when baffles are installed. One is the disappearance of sine-wave-like free-surface fluctuation shown in Fig. 10 that does commonly occur in the liquid sloshing in tanks without baffle. As a result, the free surface fluctuation of the baffle case becomes more complicated and the peak elevation does not always occur either at the right or left end. Another thing is the apparent difference in liquid flow such that the liquid motion is active above the baffle but almost dull below the baffle. This difference does also lead to the clear distinction between hydrodynamic pressure distributions above and below the baffle, as represented in Fig. 12. Where, we see the strict discontinuity in the dynamic pressure distribution across the baffle as well as the relative slowness in its vertical and temporal variations below the baffle. One can easily infer that these two features by baffle are strongly influenced by the baffle opening width. Next four figures starting from Fig. 13 represent the parametric results of the liquid sloshing with respect to the relative installation height hB =hL and the relative opening width dB =d of baffles: the peak elevation height, the maximum dynamic forces and moment acting on the tank during the observation period. Referring to Fig. 13, we first see that the peak elevation height fmax of free surface increases as a whole in proportional to the opening width reduction, regardless of the installation height. What is worse, it becomes larger than one of the no baffle case when dB =d is smaller than the specific critical values depending on the installation height. It is because the opening width reduction strengths the flow separation effect and which in turn decreases the fundamental sloshing frequency [5], while reminding that the free surface elevation is in the reverse relation with the fundamental sloshing frequency [1]. On the other hand, the baffle installation height shows a reciprocal parametric effect on the peak elevation height, starting from a critical opening width dB =d ¼ 0:5. In other words, the peak elevation height becomes larger with the installation height when dB =d < 0:5 and vice versa when dB =d > 0:5.
Fig. 11. Sloshing flow patterns during one period (hB =hL ¼ 0:7, dB =d ¼ 0:5).
2594
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
Fig. 12. Variation of hydrodynamic pressure distribution during one period (hB =hL ¼ 0:7, dB =d ¼ 0:5).
Fig. 13. Peak elevation height fmax of free surface versus the relative height and the relative opening width of baffles.
In the former case (i.e. when dB =d < 0:5), the flow separation effect by baffle is dominated, so that the liquid part participating in the sloshing motion becomes smaller as the baffle location becomes higher. Which in turn leads to the decrease of the fundamental sloshing frequency, as mentioned just above. Meanwhile, in the latter case (i.e. when dB =d > 0:5), the flow separation effect becomes more week with the opening width increase, while the deviation from the sine-wave-like free-surface configuration increases (except for dB =d reasonably close to unity). Moreover, the higher the baffle position is the more complex the free surface fluctuation is. It is not hard to realize that the peak elevation height decreases as this complexity increases. Which can be confirmed by comparing with the no baffle case shown in Fig. 13, where the minimum fmax is obtained when both dB =d and hB =hL are 0.7. As shown in Fig. 14(a), the maximum dynamic force Fx acting on the tank for the baffled case is smaller than one of the no baffle case, so that the dynamic damping effect of baffle is clearly observed. As a whole, the maximum x-force decreases in proportional to the installation height and the opening width reduction, except for the specific values of two baffle parameters. As observed from Figs. 11–13, the opening width reduction increases the flow separation effect that the liquid motion below the baffle becomes more inactive.
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
2595
Fig. 14. Maximum hydrodynamic force and moment versus the relative height and the relative opening width of baffles (dashed line: no baffle case): (a) Fx ; (b) Mz .
As well, the liquid portion showing inactive flow below the baffle becomes larger as the baffle is installed higher position. As shown in Fig. 14(a), the maximum dynamic force Fx acting on the tank for the baffled case is smaller than one of the no baffle case, so that the dynamic damping effect of baffle is clearly observed. As a whole, the maximum x-force decreases in proportional to the installation height and the opening width reduction, except for the specific values of two baffle parameters. As observed from Figs. 11–13, the opening width reduction increases the flow separation effect that the liquid motion below the baffle becomes more inactive. As well, the liquid portion showing inactive flow below the baffle becomes larger as the baffle is installed higher position. However, this trend is not shown any more when the baffle location exceeds a specific height and the opening width becomes smaller than a specific value. Referring to Figs. 12 and 13, the excessive flow separation by baffle installed higher position causes the extreme fluctuation in the free surface sloshing and the hydrodynamic pressure above the baffle, and which in opposition causes the increase of Fx . The maximum x-force reaches the lowest value when dB =d is 0.5 and hB =hL is 0.7. This parametric trend of Fx to the two baffle parameters gives rise to the similar variation of Mz , as represented in Fig. 14(b), where the maximum z-moment for the baffled case is shown to be smaller than one of the no baffle case regardless of the choice of dB =d and hB =hL . Meanwhile, Mz reaches the lowest level near dB =d ¼ 0:4 and hB =hL ¼ 0:65. As shown in Fig. 15(a), the maximum y-force increases consistently in proportional to the installation height and the opening width reduction. It is natural that the baffle opening width gives rise to this kind of parametric effect on the vertical hydrodynamic force acting on baffles. While, from the parametric effect of the baffle installation height, one can infer that the dynamic pressure jump across the baffle becomes larger as the baffle position becomes higher (referring to the vertical distribution of hydrodynamic pressure in Fig. 12). It is worth to note that the consistent variation of the dynamic force acting on baffles does not lead to the consistent improvement of the dynamic damping efficiency of baffle, as confirmed from Figs. 12–14. This can be also realized from the parametric variation of the maximum y-force acting on the tank bottom plate shown in Fig. 15(b). In other words, larger supporting of dynamic force by baffle does not result in the proportional decrease of dynamic force acting on the tank bottom plate. Instead, the maximum dynamic force acting on the bottom plate shows the parametric variation similar to one of the peak elevation height fmax shown in Fig. 13, so that the existence of baffle leads to the worse damping effect for
2596
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
Fig. 15. Maximum hydrodynamic force Fy versus the relative height and the relative opening width of baffles: (a) on baffles; (b) on the tank bottom plate (dashed line: no baffle case).
specific combinations of two baffle parameters. Regarding the opening width of baffle the appropriate range is shown to be 0.5–0.7, and the combination of dB =d ¼ 0:6 and hB =hL ¼ 0:7 produces the lowest level.
6. Conclusion Parametric investigation on the nonlinear liquid sloshing in baffled tank under horizontal forced excitation has been addressed. For the numerical analysis, a time-incremental nonlinear finite element method was introduced which is formulated based on the fully nonlinear potential flow theory in the semiLagrangian numerical approach. In order to secure the numerical accuracy and stability, the free surface was tracked by making use of the direct time differentiation and the predictor–corrector method, together with the volume-conservative mesh adaptation. The comparison with the existing reference solutions confirms that the present nonlinear method accurately predicts the sloshing time-history response without causing any numerical instability. Compared to the no baffle case, the liquid sloshing in baffled tank exhibited the prominent distinction between sloshing flows and hydrodynamic pressures above and below the baffle. While the liquid motion and the dynamic pressure variation above the baffle are active, those below the baffle are considerably dull. On the other hand, the free surface shows the severely complex fluctuation, while being remarkably influenced by the baffle opening width, so that traditional sine-wave-like profile does not appear any longer. Next, according to the parametric experiments to the installation height and the opening width of baffles, the following main observations are drawn. 1. The hydrodynamic force exerted on baffles exhibits the consistent parametric variation with respect to the two baffle parameters. However, larger supporting of dynamic force by baffle does not surely lead to the proportional reduction of other dynamic quantities associated with the liquid sloshing. 2. The peak free surface elevation and the maximum dynamic force acting on the tank bottom plate show the similar parametric variations to each other. As a whole, both dynamic quantities exhibit the reciprocal parametric dependence on the baffle installation height at the specific opening widths of baffle. Moreover, both become larger than those of the no baffle case for certain combinations of two baffle
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
2597
parameters. Restricting to the current model problem, the choice of 0.7 for dB =d and hB =hL leads to the lowest levels in both quantities. 3. The maximum dynamic force Fx and moment Mz do also display the quite similar parametric variations to the baffle parameters. However, differing from the above two dynamic quantities, these two are always smaller than those of the no baffle case. Both show the uniform reduction in proportional to the installation height and the opening width reduction, to some extent. But, those increase when the opening width is smaller than the specific value and the baffle location exceeds the specific height. Both dynamic quantities reach the lowest levels when dB =d is 0.4–0.5 and hB =hL is 0.65–0.7. As a result, the quantities of interest in the liquid sloshing display the strong dependence on the baffle design parameters, so that the baffle design for suppressing the liquid sloshing effect should be made through the comprehensive parametric investigation by a reliable numerical technique, as presented in the current study.
Acknowledgements The authors would like to thank the anonymous referees for valuable comments and suggestions. And the financial supports for this work by Korea Industrial Technology Foundation (June 2003–April 2006) and Pusan National University (2001.4–2004.3) are gratefully acknowledged.
References [1] J.R. Cho, J.M. Song, J.K. Lee, Finite element techniques for the free-vibration and seismic analysis of liquid-storage tanks, Finite Elem. Anal. Des. 37 (2001) 467–483. [2] A.M. Silvera, D.G. Stephens, H.W. Leonard, An experimental investigation of liquid oscillations in tanks with various baffles, NASA Technical Note D-715, 1961. [3] D.G. Stephens, Flexible baffles for slosh damping, J. Spacecraft Rockets 3 (5) (1966) 765–766. [4] J.R. Cho, S.Y. Lee, Dynamic analysis of baffled fuel-storage tanks using the ALE finite element method, Int. J. Numer. Methods Fluids 41 (2003) 185–208. [5] J.R. Cho, H.W. Lee, K.W. Kim, Free vibration analysis of baffled liquid-storage tanks by the structural-acoustic finite element formulation, J. Sound Vibr. 258 (5) (2003) 847–866. [6] G.W. Housner, Dynamic pressures on accelerated fluid containers, Bull. Seismol. Soc. Am. 47 (1957) 15–35. [7] S. Silverman, H.N. Abramson, Lateral sloshing in moving containers, NASA SP-106, 1966, pp. 13–78. [8] M. Ikegawa, Finite element analysis of fluid motion in a container, Finite Elem. Methods Flow Prob. (1974) 737–738. [9] C.W. Hirt, A.A. Amsden, J.L. Cook, An arbitrary Lagrangian Eulerian computing method for all flow speeds, J. Comput. Phys. 14 (1974) 227–253. [10] O.C. Zienkiewicz, P. Bettess, Fluid–structure dynamic interaction and wave forces. An introduction to numerical treatment, Int. J. Numer. Methods Engrg. 13 (1978) 1–16. [11] O.M. Faltinsen, A numerical non-linear method of sloshing in tanks with two-dimensional flow, J. Ship Res. 18 (4) (1978) 224– 241. [12] J.M. Kennedy, T. Belytschko, Theory and application of a finite element method for arbitrary Lagrangian–Eulerian fluids and structures, Nucl. Engrg. Des. 68 (1981) 129–146. [13] T.J.R. Hughes, W.K. Liu, T.K. Zimmermann, Lagrangian–Eulerian finite element formulation for incompressible viscous flows, Comput. Methods Appl. Mech. Engrg. 29 (1981) 329–349. [14] T. Nakayama, K. Washizu, The boundary element method applied to the analysis of two-dimensional nonlinear sloshing problems, Int. J. Numer. Methods Engrg. 17 (1981) 1631–1646. [15] J.T. Oden, N. Kikuchi, Y.J. Song, Penalty finite element method for the analysis of Stokesian flows, Comput. Methods Appl. Mech. Engrg. 31 (3) (1982) 297–330. [16] M. Souli, A. Ouahsine, L. Lewin, ALE formulation for fluid–structure interaction problems, Comput. Methods Appl. Mech. Engrg. 190 (2000) 659–675.
2598
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
[17] B. Ramaswamy, M. Kawahara, Arbitrary Lagrangian–Eulerian finite element method for unsteady, convective, incompressible viscous free surface fluid flow, Int. J. Numer. Methods Fluids 7 (1987) 1053–1075. [18] A.M. Winslow, Equipotential zoning of the interior of a three-dimensional mesh, Lawrence Radiation Laboratory, UCRL-7312, 1992. [19] Q.W. Ma, G.X. Wu, R.E. Taylor, Finite element simulation of fully non-linear interaction between vertical cylinders and steep waves. Part 1: Methodology and numerical procedure, Int. J. Numer. Methods Engrg. 36 (2001) 265–285. [20] A. Soulaimani, Y. Saad, An arbitrary Lagrangian–Eulerian finite element method for solving three-dimensional free surface flows, Comput. Methods Appl. Mech. Engrg. 162 (1998) 79–106. [21] C.W. Hirt, B.D. Nichols, Volume of fluids method for the dynamics of free boundaries, J. Comput. Phys. 39 (1981) 201–225. [22] S.J. Childs, B.D. Reddy, Finite element simulation of the motion of a rigid body in a fluid with free surface, Comput. Methods Appl. Mech. Engrg. 175 (1999) 99–120. [23] K. Yamamoto, M. Kawahara, Structural oscillation control using tuned liquid damper, Comput. Struct. 71 (1999) 435–446. [24] J. Sung, H.G. Choi, J.Y. Yoo, Time-accurate computation of unsteady free surface flows using an ALE-segregated equal-order FEM, Comput. Methods Appl. Mech. Engrg. 190 (2000) 1425–1440. [25] M. Cruchaga, D. Celentano, T. Tezduyar, A moving Lagrangian interface technique for flow computations over fixed meshes, Comput. Methods Appl. Mech. Engrg. 191 (2001) 525–543. [26] T. Okamoto, M. Kawahara, Two-dimensional sloshing analysis by Lagrangian finite element method, Int. J. Numer. Methods Fluids 11 (1990) 453–477. [27] W. Chen, M.A. Haroun, F. Liu, Large amplitude liquid sloshing in seismically excited tanks, Earthq. Engrg. Struct. Dyn. 25 (1996) 653–669. [28] G.X. Wu, Q.W. Ma, R.E. Taylor, Numerical simulation of sloshing waves in a 3D tank based on a finite element method, Appl. Ocean Res. 20 (1998) 337–355. [29] I.G. Currie, Fundamental Mechanics of Fluids, McGraw-Hill, New York, 1974. [30] D.G. Gill, M.R. Cullen, Advanced Engineering Mathematics, PWS Publishing Co., Boston, 1992. [31] P. Saucez, W.E. Schiesser, A.V. Wouver, Upwinding in the method of lines, Math. Comput. Simul. 56 (2001) 171–185.
Numerical study on liquid sloshing in baffled tank by nonlinear finite element method J.R. Cho *, H.W. Lee C&AM Lab, School of Mechanical Engineering, Pusan National University, 30 Jangjeon-Dong, Kumjung-Ku, Pusan 609-735, South Korea Received 2 June 2003; received in revised form 4 December 2003; accepted 8 January 2004
Abstract This paper introduces a velocity–potential-based nonlinear finite element method for the accurate simulation of the large amplitude liquid sloshing in two-dimensional baffled tank subject to horizontal forced excitation. As well, the effects of baffle on the nonlinear liquid sloshing are parametrically examined. The free surface configuration is tracked by a direct time differentiation of the convective term in the kinematic boundary condition with the help of the four-step predictor–corrector methods. The flow velocity is interpolated from the velocity potential by a second-order least square method. And, the finite element mesh is adapted in a semi-Lagrangian approach combined with the free-surface height correction, in order to keep the mesh regularity and the total liquid volume unchanged. By comparing with the available reference solutions, we verify the numerical accuracy and stability of the nonlinear finite element method proposed. Through the numerical experiments performed by varying the installation height and the opening width of baffles, the hydrodynamic characteristics of the large-amplitude liquid sloshing are parametrically investigated. 2004 Elsevier B.V. All rights reserved. Keywords: Baffled tank; Nonlinear liquid sloshing; Direct time differentiation; Semi-Lagrangian free surface tracking; Baffle installation height; Opening width
1. Introduction Free surface fluctuation called liquid sloshing is the most prominent phenomenon of liquid motion in either stationary or moving tanks subject to forced external perturbations. Representative examples can be found from liquefied natural gas (LNG) carriers in shipping service and aboveground liquid-storage containers excited by earthquake wave [1]. When the liquid sloshing occurs, hydrodynamic forces acting on the container display the severe fluctuation in their magnitudes and distributions. Moreover, an excessive liquid sloshing may cause the structural failure or/and the manipulation loss, and which frequently leads to the tremendous loss of human, economic and environmental resources.
*
Corresponding author. Tel.: +82-51-510-2467; fax: +82-51-514-7640. E-mail address: [email protected] (J.R. Cho).
0045-7825/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2004.01.009
2582
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
In view of this critical effect of the liquid sloshing on liquid–structure interaction systems, very intensive research efforts have been made to suppress it from several decades ago [2,3], as a result, a number of passive devices for suppressing the liquid motion have been devised. Among them, the flexible baffle manufactured usually with metal is being widely used thanks to the practical convenience and the flexible adjustment in its installation to the container. However, the suppression efficiency of baffle is strongly affected by several design parameters, particularly the installation height and the opening width of baffle, as demonstrated through our previous works related [4,5]. Moreover, the suitable values of baffle parameters are dependent of the problem under consideration, so that the determination of appropriate baffle parameters needs to be made based upon the sufficient parametric investigation on the liquid sloshing in baffled tank. As is widely known, liquid sloshing, as a representative example of fluid–structure interactions, is an intrinsic nonlinear problem in which the free surface configuration is not known a priori. This nonlinearity becomes higher when baffles are introduced as well as when the excitation frequency is close to the natural sloshing frequencies or the excitation amplitude is large. Nevertheless, to the best of our literature survey, the early studies before 1970s had relied on the linear theory with the assumption of small amplitude sloshing [6,7]. However, after mid 1970s, active studies for developing the numerical techniques have been made to successfully implement the nonlinear fluid–structure interaction problems [8–15]. The main concern in the numerical implementation of fluid–structure interactions exhibiting a remarkable problem nonlinearity has focused on the insurance of the numerical stability and accuracy. Because the excessive distortion of fluid elements causes the numerical instability or even the computation incompleteness [16], and a small error in the time integration accumulates with time stage. Meanwhile, the choice of state variables (i.e. the underlying governing equations) and the continuity constraint enforcement for incompressible flows [15,17] characterize the numerical formulations and approximations. Besides, several another issues such as the coupling between solid and fluid have been also intensively studied for the reliable numerical simulation. The mesh distortion and the mesh adaptation are strongly related to the choice of kinematic description approach (e.g. Lagrangian, Eulerian, or ALE (arbitrary Lagrangian–Eulerian) [9]). Currently, this problem can be resolved, to a large extent, by employing the ALE method with the help of a suitable remeshing and smoothing algorithm [18]. Regarding this matter, the reader may refer to [4,9,12,13,16,17] for more detailed discussion. On the other hand, the fluid flow has been traditionally formulated by either the potential flow theory [5,8,11,14,19] or the full Navier–Stokes equations. In the former case, the free surface configuration is tracked by time-integrating the nonlinear kinematic and dynamic conditions, and the flow velocity and the dynamic pressure are usually interpolated from the velocity potential approximated. However, in the latter case in which both flow velocity and dynamic pressure are the primary state variables, the flow boundary identification varies to the choice of kinematic description approach. Referring to the papers by Soulaimani and Saad [20] and Cho and Lee [4], the boundary tracking in both Lagrangian and ALE approaches is straightforward because fluid mesh moves exactly with fluid particles. But, in the Eulerian approach the boundary tracking requires any special characteristic (volume fraction) function [21] with the mesh refinement along the flow boundary. The fluid mass conservation in incompressible flows is enforced in several ways. Among them are a penalty method or an iterative augmented Lagrangian method (see Oden et al. [15] and Childs and Reddy [22]), a correction of the free surface height [4,23], and an iterative split algorithm combined with the fractional step method [24]. Regarding the mesh interface between solid and fluid, which is beyond the scope of the current work, the reader may refer to the papers by Childs and Reddy [22] and Cruchaga et al. [25]. In the present study, we intend to introduce a nonlinear finite element method based on the fully nonlinear potential flow theory which accurately simulate the large amplitude sloshing flow in twodimensional baffled tank subject to horizontal forced excitation. In particular, the main concern focuses
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
2583
on the insurance of the numerical stability and accuracy by suppressing the numerical dissipation and the total liquid volume change, and the parametric investigation of the baffle effects on the nonlinear liquid sloshing. This kind of problem was researched by many investigators such as Ikegawa [8] and Faltinsen [11]. Restricting to the free surface tracking in the velocity-potential-based nonlinear sloshing analysis, Nakayama and Washizu [14] proposed a nonlinear boundary element method by introducing an error correction term into the dynamic boundary condition, in order to suppress the numerical instability and dissipation. Okamoto and Kawahara [26] introduced a velocity correction method into the Lagrangian finite element formulation, and they justified their numerical method by comparing the numerical solutions with independent experimental results. Chen et al. [27] employed the finite difference scheme based on the Crank–Nicolson time marching scheme in which a second-order numerical dissipation term was added to the kinematic boundary condition. Wu et al. [28] employed the open trapezoidal scheme for the stable time integration of nonlinear liquid sloshing in three-dimensional tank. On the other hand, we in the present study track the free surface configuration by directly time differentiating the convective term in the kinematic boundary condition with the help of the four-step predictor– corrector methods with the accuracy of OðDt4 Þ. As well, we introduce a semi-Lagrangian mesh adaptation and the liquid height correction algorithm, while assuming the container be a rigid body, in order to keep the mesh regularity and the total liquid volume unchanged. Furthermore, the flow velocity is interpolated from the velocity potential by the second-order least square method, and the dynamic pressure field is approximated with the same finite element framework constructed for the velocity potential approximation. The numerical stability and accuracy of the nonlinear finite element method introduced are verified through the comparison with the existing reference solutions. We also perform the parametric experiments by varying the installation height and the opening width of baffles, in order to examine the effects of baffle on the nonlinear liquid sloshing. Major quantities of interest are the maximum free-surface elevation, hydrodynamic force and moment resultants, and the variations of flow profile and hydrodynamic pressure. This paper is organized as follows. The problem formulation of nonlinear liquid sloshing in baffled tank is given in Section 2. A time-incremental nonlinear finite element approximation based on the direct time differentiation and the semi-Lagrangian approach and the evaluation of hydrodynamic forces and moments are presented respectively in Sections 3 and 4. Next, the verification and parametric numerical results are addressed in Section 5, and the final conclusion is made in Section 6.
2. Nonlinear liquid sloshing in baffled tank Referring to Fig. 1, we consider the sloshing flow of inviscid incompressible liquid in a two-dimensional baffled rigid tank of width d which is subjected to forced horizontal excitation. Where, a pair of rigid baffles with the opening width dB is installed at a height hB and liquid is filled up to a height hL in the stationary condition. Boundary of the liquid domain Xt 2 R2 is composed of the free surface oXtF and the liquid– structure interface oXtI , in which the superscript t indicates that both the liquid domain and the liquid boundary vary with time during the observation time interval ð0; t. For the sake of mathematical description, we use two Cartesian co-ordinate systems, fX0 g fixed in space and fX g moving with the tank, such that their origins are commonly positioned at the center of the stationary free surface and their axes are aligned parallel to each other. Then, a quantity Q can be defined in ~ y; tÞ. By denoting us ðtÞ ¼ fdxs =dt; dy s =dtg be either of two co-ordinate systems such that Qðx0 ; y0 ; tÞ ¼ Qðx; the rigid tank velocity, we have the relations between two systems: rQjfX0 g ¼ rQjfX g ;
ð1Þ
2584
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
Fig. 1. Baffled rigid tank under horizontal forced excitation.
oQ oQ ¼ us rQ: ot fX0 g ot fX g
ð2Þ
Hereafter, the subscript 0 refers to the quantities measured in the fixed Cartesian co-ordinates. When we denote /ðx0 ; y0 ; tÞ be the velocity potential, the continuity condition is described by the Laplace equation: r2 / ¼ 0;
in Xt :
ð3Þ
On the liquid–structure interface, the velocity potential should satisfy r/ n ¼ us n;
on oXtI
ð4Þ
with the outward unit vector n normal to the structure boundary. In the fixed Cartesian co-ordinates, the kinematic and dynamic conditions on the free surface are as follows: og0 o/ of0 o/ þ ¼ 0; ox0 ox0 oy0 ot
on oXtF ;
o/ 1 þ r/ r/ þ gf0 ¼ 0; ot 2
on oXtF ;
ð5Þ ð6Þ
in which f0 ðx0 ; tÞ denotes the free surface elevation at time t. Eq. (5) is derived from the kinematic condition: o/=oy0 ¼ Df0 =Dt on the free surface, while Eq. (6) is the Bernoulli equation when the atmospheric pressure is assumed to be zero on the free surface [29]. Since the free surface fluctuation is a relative motion to the container the moving co-ordinates is usually chosen. By denoting fðx; tÞ ¼ f0 ðx0 ; tÞ y s ðtÞ be the free surface elevation from the reference stationary level, together with Eqs. (1) and (2), above free surface conditions can be rewritten as of o/ of o/ þ ð7Þ usx usy ¼ 0; ot ox ox oy o/ 1 us r/ þ r/ r/ þ gðf þ y s Þ ¼ 0: ot 2
ð8Þ
In the moving co-ordinate system, the velocity potential /ðx; y; tÞ can be split into two parts, /h ðx; y; tÞ due to the internal sloshing flow and /p ðx; y; tÞ due to the tank motion, such that / ¼ /h þ /p ;
/p ¼ xusx þ yusy :
ð9Þ
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
2585
Substituting Eq. (9) into Eqs. (7) and (8), we respectively have of o/h of o/h þ ¼ 0; ot ox ox oy
ð10Þ
dusy 1 s 2 o/h 1 dus þ r/h r/h þ gf þ x x þ f ¼ ju j gy s : 2 2 ot dt dt
ð11Þ
It is worth to note that two terms in the RHS of Eq. (11) vanish because not only both are co-ordinateindependent quantities but also the velocity potential is a relative quantity in space. Then, the initial-boundary-value problem of nonlinear sloshing flow in the moving co-ordinate system is formulated as follows: r2 /h ¼ 0;
in Xt
ð12Þ
with initial conditions /h ¼ xusx ð0Þ yusy ð0Þ; f ¼ 0;
in X0 ;
ð13Þ
on oX0F ;
ð14Þ
and boundary conditions r/h n ¼ 0;
on oXtI ;
of o/ of o/h ¼ h þ ; ot ox ox oy
ð15Þ on oXtF ;
dusy o/h 1 dus ¼ r/h r/h gf x x f ; 2 ot dt dt
ð16Þ on oXtF :
ð17Þ
3. Nonlinear finite element approximation In order to solve the above nonlinear sloshing problem we divide the observation time interval t into a finite number of sub-intervals such that tn ¼ nDt (n ¼ 0; 1; 2; . . .). Starting from the initial solutions /0h and f0 , we successively solve the Laplace equation (12) defined in current flow domain Xn and boundary oXn to seek /nh with the boundary condition (15). From which we compute the flow velocity field vn and perform the free surface tracking. By time-integrating the kinematic and dynamic boundary conditions (16) and (17), together with the mass-conservative semi-Lagrangian remeshing, we identify fnþ1 , Xnþ1 and oXnþ1 . In this manner, the sloshing time response is to be numerically analyzed. 3.1. Interpolation of flow velocity field Referring to Fig. 2, we consider a time stage tn at which /nh and vn are to be sought with the geometry and field quantities determined at the previous time stage tn1 . According to the virtual work principle, we have the weak form of Eq. (12): find /nh such that Z rw r/nh dv ¼ 0 ð18Þ Xn
2586
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
Fig. 2. Flow status representation at time stage tn .
for every admissible velocity potential wðx; yÞ. The corresponding essential boundary condition on oXnF is /nh determined previously by the free surface tacking method which will be described later. With nine-node N n , we approximate the Lagrange-type finite element basis functions fNi gi¼1 and the nodal potential vector / h velocity potential n: /nh ¼ N / h
ð19Þ
Substituting Eq. (19) into Eq. (18) leads to n ¼ 0 K/
ð20Þ
with the matrix K defined by Z T K¼ ðrNÞ ðrNÞ dv:
ð21Þ
h
Xn
Once the velocity potential /nh is approximated, the intermediate velocity field ~vn is computed according to n: ~vn ¼ rN / h
ð22Þ
0
Owing to the C -continuity in the velocity potential field, the intermediate velocity field becomes discontinuous along the element interfaces. In order to secure the continuity in the velocity field, one may consider the linear interpolation of the intermediate velocity values at Gauss integration points. However, a small error in the simple velocity interpolation does not only affect the accuracy of the free surface tracking but also accumulate with the time stage. In order to suppress this crucial problem, we first interpolate elementwise velocity fields from ~vn with second-order polynomials according to the least square method. Let us denote ^vna (a ¼ x; y) be the element-wise velocity components to be interpolated: ^vna ¼ a1 þ a2 n þ a3 g þ a4 ng þ a5 n2 þ a6 g2 :
ð23Þ
And we define the element-wise error by Ea ¼
9 X
2
f^vna ðn‘ ; g‘ Þ ~vna ðn‘ ; g‘ Þg ;
a ¼ x; y;
ð24Þ
‘¼1
where ‘ stands for (3 · 3) Gauss integration points. Then, six coefficients involved in each velocity component interpolation (23) are to be determined from six simultaneous equations which are constructed from oEa ¼ 0; oak
k ¼ 1; 2; . . . ; 6:
ð25Þ
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
2587
After the element-wise velocity fields are interpolated, we compute the nodal velocities of all elements and take average of nodal values for the common nodes. By interpolating the nodal velocities vna with the same finite element basis functions used as for the velocity potential approximation (19), we finally have the global continuous velocity field: vna ¼ N vna ;
a ¼ x; y:
ð26Þ
3.2. Time marching of free surface Once both the velocity potential /nh and the flow velocity field vn are approximated, we integrate Eq. (16) to predict the free surface fnþ1 and Eq. (17) to determine the essential boundary condition /hnþ1 . As a result, the flow geometry and the boundary condition required for the next-step computation are to be established. A mass-conservative remeshing process follows the free surface tracking, while keeping the mesh regularity fairly, in order to update the liquid mesh. Here, the time integration of the kinematic boundary condition (16) may lead to the numerical divergence, the high frequency wiggle or the numerical dissipation, unless any special care is paid. This crucial problem is owing to the convection term in the RHS of Eq. (16), and its influence increases in proportional to the sloshing amplitude. In order to suppress the numerical instability, several kinds of correction techniques have been proposed by Nakayama and Washizu [14], Okamoto and Kawahara [26], and Chen et al. [27], as mentioned before. For the present study, however, we introduce a direct time differentiation scheme to track the free surface. Referring to Fig. 3, we predict the free surface configuration after the small time interval dt (dt ¼ aDt, 0 < a < 1) from the current time tn according to o/h o/h 0 0 xi ¼ xi þ dt ; yi ¼ yi þ dt ; ð27Þ ox oy ðxi ;yi Þ
ðxi ;yi Þ
where ðxi ; yi Þ denotes co-ordinates of node i located on the free surface at time tn . And, the elevation of node i during the time interval dt is calculated through f0i ¼ yj0 þ
0 yjþ1 yj0 ðxi x0j Þ: x0jþ1 x0j
ð28Þ
Then, we directly calculate the time derivative of the free surface elevation in Eq. (16) using the alternative relation given by of df f0 fn : ð29Þ ¼ ot n dt n dt
Fig. 3. Variation of the free surface during the small time interval dt.
2588
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
The parameter a influences the numerical dissipation and the incompressibility condition, so that its choice should be carefully made. As presented in Section 5.1, the former is getting improved as a becomes smaller, but the latter deteriorates either as a approaches zero or unity. In the current study, we choose a suitable a based upon its parametric effects on both the numerical dissipation and the liquid volume conservation. Next, we apply the predictor–corrector method composed of the explicit Adams–Bashforth four-step scheme [30] and the implicit Adams–Moulton three-step scheme [30] given respectively by ~fnþ1 ¼ fn þ Dt 55 of 59 of þ 37 of 9 of ; ð30Þ 24 ot ot ot ot n
Dt o~f n 9 ¼f þ 24 ot
n1
"
fnþ1
þ 19 nþ1
n2
of of of 5 þ ot n ot n1 ot n2
n3
# ð31Þ
in order to build up the free surface configuration oXFnþ1 at the next time stage tnþ1 . After that, we relocate the free surface nodes onto the new free surface by moving those only in the vertical direction, as depicted in Fig. 4(a). As well as, we calculate the total liquid area according to the trapezoidal rule Atot ¼
K X
Ai ;
Ai ¼ ðxiþ1 xi Þ
i¼1
fiþ1 þ fi : 2
ð32Þ
The mean free-surface elevation fmean can be calculated as follows: fmean ¼ Atot =d, and the nodal free surface heights finþ1 are to be corrected [4], referring to Fig. 4(b), such that fnþ1 ¼ fnþ1 fmean : i i
ð33Þ
Referring to Fig. 5, we next relocate all internal nodes such that the nodes located on the same vertical line should be in almost uniform spacing. In this manner, the free surface is tracked with the total liquid volume unchanged and the liquid mesh is adapted to keep the mesh regularity fairly. In the next step, we integrate Eq. (17) to specify the velocity potential /hnþ1 to the new free surface oXFnþ1 tracked. However, the time derivative of the free surface velocity potential in the semi-Lagrangian liquid mesh fXM g is related to one in the moving co-ordinate system fX g such that o/h o/h of o/h : ð34Þ ¼ ot oy ot ot fX g
fXM g
Thus, the time derivative of the free surface velocity potential with respect to the liquid mesh, at time stage tn , is rewritten as
Fig. 4. Volume-conservative free surface tracking: (a) relocation of free surface nodes; and (b) correction of the free surface height.
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
2589
Fig. 5. Liquid mesh adaptation by semi-Lagrangian remeshing scheme.
dusy o/h of o/nh 1 dusx n n n n r/ f ; ¼ r/ gf x h h ot n oy 2 ot n dt dt
on oXnF :
ð35Þ
After evaluating this time derivative, we apply the predictor–corrector methods (30) and (31) used for the free surface tracking once again, in order to determine the velocity potential /hnþ1 on the free surface oXnþ1 F . 4. Hydrodynamic forces and moment The total pressure pðx; y; tÞ in the liquid domain is expressed as p ¼ pd qgy;
ð36Þ
where q denotes the liquid density. Meanwhile, the dynamic pressure pd ðx; y; tÞ is computed according to dusy o/h 1 dusx þy pd ¼ q þ r/h r/h þ x : ð37Þ 2 ot dt dt In above equation, one can realize that an additional finite element approximation is required to obtain the time derivative of the velocity potential /h . Letting o/h =ot be v, we have the following Laplace equation from Eq. (12): r2 v ¼ 0;
in Xt
ð38Þ
with boundary conditions rv n ¼ 0;
on oXtI ;
ð39Þ
dusy 1 dus v ¼ r/h r/h gf x x f ; 2 dt dt
on oXtF :
ð40Þ
At any time stage tn , the essential boundary conditions (15) and (17) are to be time-integrated during the finite element approximation of the velocity potential /nh . So that, this boundary value problem can be easily solved with the same liquid mesh after the velocity potential approximation at time stage tn . Referring to Fig. 6, the force resultants acting on the baffled tank due to the liquid sloshing are computed through Z Z Fx ¼ p dy p dy; ð41Þ CR T
CL T
Z
Z
Fy ¼
CBT
p dx
Cþ B
p dx þ
Z p dx: C B
ð42Þ
2590
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
Fig. 6. Division of the liquid–structure interaction.
And the moment resultant about the origin O is obtained according to Z Z Z Z Z py dy py dy px dx px dx þ px dx: Mz ¼ CL T
CR T
Cþ B
CBT
ð43Þ
C B
Of course, the force and moment resultants are in function of time.
5. Numerical experiments According to the numerical formulae described so far we developed a test program coded in Fortran which is interfaced with MSC/Patran for visualization. The forced excitation of tank is taken as a horizontal harmonic movement xs ðtÞ ¼ a sin xt, but the amplitude and the frequency are variable depending on the simulation case. Referring to Fig. 1, the tank width, the liquid fill height, and the installation height and the opening width of baffles are also taken variable. Interior liquid is initially discretized with uniform mesh of nine-node quadratic elements and the time step size Dt is set by 5 · 103 s throughout the experiments. 5.1. Verification of numerical accuracy and stability We first perform the verification experiments with liquid-filled tanks without baffle, in order to validate the time-incremental nonlinear finite element method introduced and the test program developed. Two simulation cases are considered, as presented in Table 1, where case 1 was taken by Nakayama and Washizu [14] while case 2 by Chen et al. [27] and Okamoto and Kawahara [26]. In order to determine a suitable value of a, we first examine the effects of a on the numerical dissipation and the total liquid volume. This experiment was done with case 2 because the sloshing nonlinearity is higher in case 2 than case 1, as will be seen later. It is worthy noting that this preliminary simulation was performed without correcting the free-surface height according to Eq. (33) so that errors in both quantities accumulate with time. As shown in Fig. 7(a), the numerical dissipation is getting suppressed as a becomes smaller such that the time history responses of the free surface elevation between a ¼ 0:1 and 0.01 do not show a noticeable distinction. On the other hand, the total liquid volume, referring to Fig. 7(b), considTable 1 Simulation cases for the verification experiment Case 1 2
Geometry
Excitation
Tank width d (m)
Liquid height hL (m)
Amplitude a (m)
Frequency x (rad/s)
0.9 (16) 1.0 (24)
0.6 (12) 0.5 (12)
0.002 )0.0093
5.5 5.31123
Subdivision number to construct uniform liquid meshes.
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
2591
Fig. 7. The effects of a (case 2): (a) on the free surface elevation (at the right end x ¼ 0:5 m); (b) on the liquid volume conservation.
erably decreases when a is unity, and contrary it slightly increases when a 6 0:1. From these features of both quantities, we set a by 0.1 for the current study. Fig. 8 shows the comparison with the solution by Nakayama and Washizu [14] that has been frequently used as a reference solution in comparative numerical experiments. As mentioned before, they improved the numerical accuracy by introducing an error correction term into the dynamic boundary condition, and they justified the numerical accuracy by comparing with the FEM solution by Ikegawa [8]. The time history response of free surface elevation was measured at the right end x ¼ 0:45 m. In view of the fundamental sloshing frequency x0 ¼ 5:761 rad/s of case 1, the excitation frequency x is taken by x ¼ 0:9547x0 . One can observe that the numerical solution by the proposed method is in an excellent agreement with the reference solution. Comparative numerical results carried out with case 2 are given in Fig. 9 in which the free surface elevation is taken at the right end x ¼ 0:5 m. The sloshing response denoted by upwind scheme in Fig. 9(b) was obtained by our test program, but for which a five-point upwind scheme [31] was used to deal with the convection term in the kinematic boundary condition (16) instead of our direct differentiation scheme (29). The comparison with the Chen et al.Õs solution, for which the second-order numerical dissipation term was introduced into the kinematic boundary condition [27], shows that the present method provides the accurate sloshing response without causing any numerical instability. However, the case with the upwind scheme displays a crucial numerical dissipation so that the free surface elevation becomes considerably smaller after the peak elevation time, even though the numerical divergence or the high frequency wiggle is not observed. Fig. 10 represents the comparison of free surface configurations at six different time stages with the Okamoto and KawaharaÕs experiment. The time marching scheme introduced in the current study
Fig. 8. Comparison with Nakayama and Washizu [14]: time history response of the free surface elevation (case 1).
2592
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
Fig. 9. Time history response of the free surface elevation (case 2): (a) direct time differentiation scheme; (b) five-point upwind scheme.
Fig. 10. Comparison of free surface sloshing flows at different time steps (case 2).
accurately predicts the free surface fluctuation to the lapse of time that is in a good agreement with the experimental result up to the time of observation. Thus, the numerical accuracy and stability of the nonlinear finite element method introduced have been sufficiently verified. 5.2. Parametric experiments In order for the parametric analysis to the installation height and the opening width of baffles, we consider a tank with a pair of baffles. Referring to Fig. 1, the liquid height and the tank width are set by d ¼ 1:0 m and hL ¼ 1:0 m, and the liquid mesh is constructed with 16 · 24 uniform 9-node quadrilateral elements. Similar to the previous verification experiments, we apply the horizontal harmonic excitation with the amplitude a ¼ 0:005 m and the excitation frequency x ¼ 0:999x0 . Here, the fundamental sloshing frequency x0 varies with the two baffle parameters [5]. For the current parametric investigation, we consider the relative baffle heights hB =hL between 0.3 and 0.7 and the relative opening widths dB =d between 0.2 and 0.8, respectively, and we set the observation time duration by 40 s.
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
2593
Fig. 11 shows the variation of sloshing flow during one period for the case of hB =hL ¼ 0:7 and dB =d ¼ 0:5. We can observe two major features in the liquid sloshing when baffles are installed. One is the disappearance of sine-wave-like free-surface fluctuation shown in Fig. 10 that does commonly occur in the liquid sloshing in tanks without baffle. As a result, the free surface fluctuation of the baffle case becomes more complicated and the peak elevation does not always occur either at the right or left end. Another thing is the apparent difference in liquid flow such that the liquid motion is active above the baffle but almost dull below the baffle. This difference does also lead to the clear distinction between hydrodynamic pressure distributions above and below the baffle, as represented in Fig. 12. Where, we see the strict discontinuity in the dynamic pressure distribution across the baffle as well as the relative slowness in its vertical and temporal variations below the baffle. One can easily infer that these two features by baffle are strongly influenced by the baffle opening width. Next four figures starting from Fig. 13 represent the parametric results of the liquid sloshing with respect to the relative installation height hB =hL and the relative opening width dB =d of baffles: the peak elevation height, the maximum dynamic forces and moment acting on the tank during the observation period. Referring to Fig. 13, we first see that the peak elevation height fmax of free surface increases as a whole in proportional to the opening width reduction, regardless of the installation height. What is worse, it becomes larger than one of the no baffle case when dB =d is smaller than the specific critical values depending on the installation height. It is because the opening width reduction strengths the flow separation effect and which in turn decreases the fundamental sloshing frequency [5], while reminding that the free surface elevation is in the reverse relation with the fundamental sloshing frequency [1]. On the other hand, the baffle installation height shows a reciprocal parametric effect on the peak elevation height, starting from a critical opening width dB =d ¼ 0:5. In other words, the peak elevation height becomes larger with the installation height when dB =d < 0:5 and vice versa when dB =d > 0:5.
Fig. 11. Sloshing flow patterns during one period (hB =hL ¼ 0:7, dB =d ¼ 0:5).
2594
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
Fig. 12. Variation of hydrodynamic pressure distribution during one period (hB =hL ¼ 0:7, dB =d ¼ 0:5).
Fig. 13. Peak elevation height fmax of free surface versus the relative height and the relative opening width of baffles.
In the former case (i.e. when dB =d < 0:5), the flow separation effect by baffle is dominated, so that the liquid part participating in the sloshing motion becomes smaller as the baffle location becomes higher. Which in turn leads to the decrease of the fundamental sloshing frequency, as mentioned just above. Meanwhile, in the latter case (i.e. when dB =d > 0:5), the flow separation effect becomes more week with the opening width increase, while the deviation from the sine-wave-like free-surface configuration increases (except for dB =d reasonably close to unity). Moreover, the higher the baffle position is the more complex the free surface fluctuation is. It is not hard to realize that the peak elevation height decreases as this complexity increases. Which can be confirmed by comparing with the no baffle case shown in Fig. 13, where the minimum fmax is obtained when both dB =d and hB =hL are 0.7. As shown in Fig. 14(a), the maximum dynamic force Fx acting on the tank for the baffled case is smaller than one of the no baffle case, so that the dynamic damping effect of baffle is clearly observed. As a whole, the maximum x-force decreases in proportional to the installation height and the opening width reduction, except for the specific values of two baffle parameters. As observed from Figs. 11–13, the opening width reduction increases the flow separation effect that the liquid motion below the baffle becomes more inactive.
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
2595
Fig. 14. Maximum hydrodynamic force and moment versus the relative height and the relative opening width of baffles (dashed line: no baffle case): (a) Fx ; (b) Mz .
As well, the liquid portion showing inactive flow below the baffle becomes larger as the baffle is installed higher position. As shown in Fig. 14(a), the maximum dynamic force Fx acting on the tank for the baffled case is smaller than one of the no baffle case, so that the dynamic damping effect of baffle is clearly observed. As a whole, the maximum x-force decreases in proportional to the installation height and the opening width reduction, except for the specific values of two baffle parameters. As observed from Figs. 11–13, the opening width reduction increases the flow separation effect that the liquid motion below the baffle becomes more inactive. As well, the liquid portion showing inactive flow below the baffle becomes larger as the baffle is installed higher position. However, this trend is not shown any more when the baffle location exceeds a specific height and the opening width becomes smaller than a specific value. Referring to Figs. 12 and 13, the excessive flow separation by baffle installed higher position causes the extreme fluctuation in the free surface sloshing and the hydrodynamic pressure above the baffle, and which in opposition causes the increase of Fx . The maximum x-force reaches the lowest value when dB =d is 0.5 and hB =hL is 0.7. This parametric trend of Fx to the two baffle parameters gives rise to the similar variation of Mz , as represented in Fig. 14(b), where the maximum z-moment for the baffled case is shown to be smaller than one of the no baffle case regardless of the choice of dB =d and hB =hL . Meanwhile, Mz reaches the lowest level near dB =d ¼ 0:4 and hB =hL ¼ 0:65. As shown in Fig. 15(a), the maximum y-force increases consistently in proportional to the installation height and the opening width reduction. It is natural that the baffle opening width gives rise to this kind of parametric effect on the vertical hydrodynamic force acting on baffles. While, from the parametric effect of the baffle installation height, one can infer that the dynamic pressure jump across the baffle becomes larger as the baffle position becomes higher (referring to the vertical distribution of hydrodynamic pressure in Fig. 12). It is worth to note that the consistent variation of the dynamic force acting on baffles does not lead to the consistent improvement of the dynamic damping efficiency of baffle, as confirmed from Figs. 12–14. This can be also realized from the parametric variation of the maximum y-force acting on the tank bottom plate shown in Fig. 15(b). In other words, larger supporting of dynamic force by baffle does not result in the proportional decrease of dynamic force acting on the tank bottom plate. Instead, the maximum dynamic force acting on the bottom plate shows the parametric variation similar to one of the peak elevation height fmax shown in Fig. 13, so that the existence of baffle leads to the worse damping effect for
2596
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
Fig. 15. Maximum hydrodynamic force Fy versus the relative height and the relative opening width of baffles: (a) on baffles; (b) on the tank bottom plate (dashed line: no baffle case).
specific combinations of two baffle parameters. Regarding the opening width of baffle the appropriate range is shown to be 0.5–0.7, and the combination of dB =d ¼ 0:6 and hB =hL ¼ 0:7 produces the lowest level.
6. Conclusion Parametric investigation on the nonlinear liquid sloshing in baffled tank under horizontal forced excitation has been addressed. For the numerical analysis, a time-incremental nonlinear finite element method was introduced which is formulated based on the fully nonlinear potential flow theory in the semiLagrangian numerical approach. In order to secure the numerical accuracy and stability, the free surface was tracked by making use of the direct time differentiation and the predictor–corrector method, together with the volume-conservative mesh adaptation. The comparison with the existing reference solutions confirms that the present nonlinear method accurately predicts the sloshing time-history response without causing any numerical instability. Compared to the no baffle case, the liquid sloshing in baffled tank exhibited the prominent distinction between sloshing flows and hydrodynamic pressures above and below the baffle. While the liquid motion and the dynamic pressure variation above the baffle are active, those below the baffle are considerably dull. On the other hand, the free surface shows the severely complex fluctuation, while being remarkably influenced by the baffle opening width, so that traditional sine-wave-like profile does not appear any longer. Next, according to the parametric experiments to the installation height and the opening width of baffles, the following main observations are drawn. 1. The hydrodynamic force exerted on baffles exhibits the consistent parametric variation with respect to the two baffle parameters. However, larger supporting of dynamic force by baffle does not surely lead to the proportional reduction of other dynamic quantities associated with the liquid sloshing. 2. The peak free surface elevation and the maximum dynamic force acting on the tank bottom plate show the similar parametric variations to each other. As a whole, both dynamic quantities exhibit the reciprocal parametric dependence on the baffle installation height at the specific opening widths of baffle. Moreover, both become larger than those of the no baffle case for certain combinations of two baffle
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
2597
parameters. Restricting to the current model problem, the choice of 0.7 for dB =d and hB =hL leads to the lowest levels in both quantities. 3. The maximum dynamic force Fx and moment Mz do also display the quite similar parametric variations to the baffle parameters. However, differing from the above two dynamic quantities, these two are always smaller than those of the no baffle case. Both show the uniform reduction in proportional to the installation height and the opening width reduction, to some extent. But, those increase when the opening width is smaller than the specific value and the baffle location exceeds the specific height. Both dynamic quantities reach the lowest levels when dB =d is 0.4–0.5 and hB =hL is 0.65–0.7. As a result, the quantities of interest in the liquid sloshing display the strong dependence on the baffle design parameters, so that the baffle design for suppressing the liquid sloshing effect should be made through the comprehensive parametric investigation by a reliable numerical technique, as presented in the current study.
Acknowledgements The authors would like to thank the anonymous referees for valuable comments and suggestions. And the financial supports for this work by Korea Industrial Technology Foundation (June 2003–April 2006) and Pusan National University (2001.4–2004.3) are gratefully acknowledged.
References [1] J.R. Cho, J.M. Song, J.K. Lee, Finite element techniques for the free-vibration and seismic analysis of liquid-storage tanks, Finite Elem. Anal. Des. 37 (2001) 467–483. [2] A.M. Silvera, D.G. Stephens, H.W. Leonard, An experimental investigation of liquid oscillations in tanks with various baffles, NASA Technical Note D-715, 1961. [3] D.G. Stephens, Flexible baffles for slosh damping, J. Spacecraft Rockets 3 (5) (1966) 765–766. [4] J.R. Cho, S.Y. Lee, Dynamic analysis of baffled fuel-storage tanks using the ALE finite element method, Int. J. Numer. Methods Fluids 41 (2003) 185–208. [5] J.R. Cho, H.W. Lee, K.W. Kim, Free vibration analysis of baffled liquid-storage tanks by the structural-acoustic finite element formulation, J. Sound Vibr. 258 (5) (2003) 847–866. [6] G.W. Housner, Dynamic pressures on accelerated fluid containers, Bull. Seismol. Soc. Am. 47 (1957) 15–35. [7] S. Silverman, H.N. Abramson, Lateral sloshing in moving containers, NASA SP-106, 1966, pp. 13–78. [8] M. Ikegawa, Finite element analysis of fluid motion in a container, Finite Elem. Methods Flow Prob. (1974) 737–738. [9] C.W. Hirt, A.A. Amsden, J.L. Cook, An arbitrary Lagrangian Eulerian computing method for all flow speeds, J. Comput. Phys. 14 (1974) 227–253. [10] O.C. Zienkiewicz, P. Bettess, Fluid–structure dynamic interaction and wave forces. An introduction to numerical treatment, Int. J. Numer. Methods Engrg. 13 (1978) 1–16. [11] O.M. Faltinsen, A numerical non-linear method of sloshing in tanks with two-dimensional flow, J. Ship Res. 18 (4) (1978) 224– 241. [12] J.M. Kennedy, T. Belytschko, Theory and application of a finite element method for arbitrary Lagrangian–Eulerian fluids and structures, Nucl. Engrg. Des. 68 (1981) 129–146. [13] T.J.R. Hughes, W.K. Liu, T.K. Zimmermann, Lagrangian–Eulerian finite element formulation for incompressible viscous flows, Comput. Methods Appl. Mech. Engrg. 29 (1981) 329–349. [14] T. Nakayama, K. Washizu, The boundary element method applied to the analysis of two-dimensional nonlinear sloshing problems, Int. J. Numer. Methods Engrg. 17 (1981) 1631–1646. [15] J.T. Oden, N. Kikuchi, Y.J. Song, Penalty finite element method for the analysis of Stokesian flows, Comput. Methods Appl. Mech. Engrg. 31 (3) (1982) 297–330. [16] M. Souli, A. Ouahsine, L. Lewin, ALE formulation for fluid–structure interaction problems, Comput. Methods Appl. Mech. Engrg. 190 (2000) 659–675.
2598
J.R. Cho, H.W. Lee / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2581–2598
[17] B. Ramaswamy, M. Kawahara, Arbitrary Lagrangian–Eulerian finite element method for unsteady, convective, incompressible viscous free surface fluid flow, Int. J. Numer. Methods Fluids 7 (1987) 1053–1075. [18] A.M. Winslow, Equipotential zoning of the interior of a three-dimensional mesh, Lawrence Radiation Laboratory, UCRL-7312, 1992. [19] Q.W. Ma, G.X. Wu, R.E. Taylor, Finite element simulation of fully non-linear interaction between vertical cylinders and steep waves. Part 1: Methodology and numerical procedure, Int. J. Numer. Methods Engrg. 36 (2001) 265–285. [20] A. Soulaimani, Y. Saad, An arbitrary Lagrangian–Eulerian finite element method for solving three-dimensional free surface flows, Comput. Methods Appl. Mech. Engrg. 162 (1998) 79–106. [21] C.W. Hirt, B.D. Nichols, Volume of fluids method for the dynamics of free boundaries, J. Comput. Phys. 39 (1981) 201–225. [22] S.J. Childs, B.D. Reddy, Finite element simulation of the motion of a rigid body in a fluid with free surface, Comput. Methods Appl. Mech. Engrg. 175 (1999) 99–120. [23] K. Yamamoto, M. Kawahara, Structural oscillation control using tuned liquid damper, Comput. Struct. 71 (1999) 435–446. [24] J. Sung, H.G. Choi, J.Y. Yoo, Time-accurate computation of unsteady free surface flows using an ALE-segregated equal-order FEM, Comput. Methods Appl. Mech. Engrg. 190 (2000) 1425–1440. [25] M. Cruchaga, D. Celentano, T. Tezduyar, A moving Lagrangian interface technique for flow computations over fixed meshes, Comput. Methods Appl. Mech. Engrg. 191 (2001) 525–543. [26] T. Okamoto, M. Kawahara, Two-dimensional sloshing analysis by Lagrangian finite element method, Int. J. Numer. Methods Fluids 11 (1990) 453–477. [27] W. Chen, M.A. Haroun, F. Liu, Large amplitude liquid sloshing in seismically excited tanks, Earthq. Engrg. Struct. Dyn. 25 (1996) 653–669. [28] G.X. Wu, Q.W. Ma, R.E. Taylor, Numerical simulation of sloshing waves in a 3D tank based on a finite element method, Appl. Ocean Res. 20 (1998) 337–355. [29] I.G. Currie, Fundamental Mechanics of Fluids, McGraw-Hill, New York, 1974. [30] D.G. Gill, M.R. Cullen, Advanced Engineering Mathematics, PWS Publishing Co., Boston, 1992. [31] P. Saucez, W.E. Schiesser, A.V. Wouver, Upwinding in the method of lines, Math. Comput. Simul. 56 (2001) 171–185.