High performence of sloshing problem in cylindrical tank with various barrels by isogeometric boundary element method

Engineering Analysis with Boundary Elements 114 (2020) 148–165

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

High performence of sloshing problem in cylindrical tank with various barrels by isogeometric boundary element method Jun Liu a,b, Quansheng Zang a,b,∗, Wenbin Ye a,b, Gao Lin a,b a b

State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China School of Hydraulic Engineering, Faculty of Infrastructure Engineering, Dalian University of Technology, China

a r t i c l e

i n f o

Keywords: Isogeometric Boundary element method NURBS Sloshing Barrel Multiply connected domain

a b s t r a c t The NURBS basis is introduced into the boundary element method (BEM) for solving the liquid sloshing problem in the multiply connected domain. An ideal liquid assumed to be inviscid, incompressible is used to study the sloshing characteristics in the cylindrical tanks with barrels of the arbitrary shapes, dimensions and number. Considering the boundary conditions at the free surface and the tank bottom as well as the constant section of the tank, the Laplace equation for potential flow problem can be decoupled into 2D Helmholtz and modified Helmholtz equations. The NURBS widely used in the CAD industry is adopted in this paper to approximate both the boundary geometry and the field variables. The proposed model is verified against the available results from some existing literatures. Meanwhile, the effects of geometric parameters, number and arrangement of internal barrels on liquid sloshing characteristics (such as the sloshing force, surface elevation and natural frequencies) in cylindrical tanks are also investigated in detail. This paper is expected to provide useful guidance for design of cylindrical tanks equipped with several internal tubular structures (such as the fast reactor with pumps and heat exchangers).

1. Introduction Liquid sloshing is a common phenomenon in the partially-filled container subjected to the external excitation, such as the fast reactor under the seismic load and the liquid storage tanks on vehicles. It may result in violent surface and even bulk turbulence, which will further generate additional dynamic loads on the tank wall and the structures inside the tank. And in some extreme cases, the tremendous dynamic loads may threaten the structural integrity of the whole system. It is noteworthy that many kinds of tanks equipped with internal structures (barrels, poles, etc.) widely used in many fields, such as storage tanks in chemical plant and pharmaceutical factory, fast reactor with pumps and heat exchangers in nuclear power station, shafts of ocean platforms even fuel rocket tanks. The internal structures of this tanks may take different forms, but their dynamic effects on liquid contained in the tanks and the entire tank system are the same, that is the barrels (poles or something others) can change the natural sloshing frequencies and obtain extra repression owing to the flow separation. Therefore, a comprehensive understanding of the influence factors on liquid sloshing in tanks with internal barrels will provide important guidance for the structural design of tank system. Studies on liquid sloshing have attracted attentions of many researchers during the past few decades. The first analysis of liquid slosh-

∗

ing was found in a NASA report by Abramson [1], which aimed to investigate the small amplitude sloshing in a fuel container of the aircraft with the linear potential flow theory. Thereafter, many investigations with regard to sloshing in different kinds of liquid containers emerged, such as sloshing motions in 2D, 3D rectangular tanks, and circular or annular cylindrical tanks were also analytically or experimentally investigated by Fujita et al. [21], Jeon et al. [33], Su et al. [53], Grotle et al. [27], Radnić et al. [47], Yue et al. [64], Raynovskyy and Timokha [48], Moradi et al. [44], Bazazzadeh et al. [4] and Kabiri et al. [35]. Frosina et al. [20] adopted a 3D Computational Fluid Dynamics (CFD) method combined with the experimental analysis to study the sloshing behavior in an actual fuel tank with complex geometry. Meanwhile, some commercial software such as ADINA [37] and ANSYS [39] were also used to simulate the sloshing response in practical tanks. As the rapid development of computer technology during recent decades, more and more numerical methods were adopted by an increasing number of researchers to investigate the sloshing characteristics in liquid tanks. For instance, the finite volume method (FVM) method for 3D rectangular tank ([7,13]), the FVM method for cylindrical tank [10], the moving particle semi-implicit (MPS) method for 2D rectangular tank [32], and the 3D BEM for cylindrical and rectangular tanks ([12,65]). However, most of the aforementioned references were investigations with regard to liquid sloshing in tanks without any internal structures,

Corresponding author at: State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China. E-mail address: [email protected] (Q. Zang).

https://doi.org/10.1016/j.enganabound.2020.02.014 Received 24 November 2019; Received in revised form 26 February 2020; Accepted 27 February 2020 Available online 19 March 2020 0955-7997/© 2020 Elsevier Ltd. All rights reserved.

J. Liu, Q. Zang and W. Ye et al.

Engineering Analysis with Boundary Elements 114 (2020) 148–165

and the dynamic force actually caused by the liquid sloshing in such tank is always extremely large, which may induce structural failure of the tank system. Thus, some special devices were developed to suppress the sloshing effect, among which the baffle might be one of the most popular equipments, and a great number of researchers have devoted themselves into the study of the suppressing effect of baffle on liquid sloshing. For instance, Zheng et al. [66] used the finite-volume method (FVM) to estimate the effects of the transverse design of baffle on suppressing the liquid sloshing in tank vehicles. Demirel and Aral [15] provide a novel slot-baffle to study the suppression effects on liquid sloshing in 2D rectangular tank. By combining the finite element method (FEM) with arbitrary Lagrangian–Eulerian (ALE) method, Kargbo et al. [36] investigated the multiphase sloshing in two-dimensional rectangular tank with a baffle or a submersed block. Liquid sloshing in vertical cylindrical tanks with annular baffles were studied by Gavrilyuk et al. [22], Wang et al. [60] , Wang et al. [59] with the analytical method, and Gedikli and Ergüven [23], Ebrahimian et al. [17], Gnitko et al. [24] with the VBEM or BEM method, as well as Wang et al. [61] with the SBFEM method. Moreover, a more extensive survey of some other literatures ([43]; Faltinsen et al.; Saoudi et al., 2013; Molin and Remy, 2013; Cho and Kim, 2016; Cho et al., 2017) evidently elucidated that the porous structure with a reasonable porosity could possess more effective repression influence on the liquid sloshing than the solid baffle. It is obvious that most of studies mentioned above can be achieved by solving the sloshing problem within convex fluid domains. Nevertheless, relatively few researchers have investigated the liquid sloshing within tanks with multiply connected domains. For examples, Jeong [34] adopted the finite Fourier transformation together with the Graf’s additional theorem and Beltrami’s theorem to investigate the influences of annular flow gap and the eccentricity of the circular cylindrical shell on the liquid natural frequencies. Drake [16] used a BEM to investigate the effect of arbitrary array of pipes on the liquid sloshing frequencies in circular cylindrical tank, and the same examples were later solved by Timokha [58] with the Trefftz method. VOF together with the CFD methods were applied by Lu et al. [41] to study the sloshing response in annular tank mounted with several circular cylinders, and experimental data were also measured by them for the verification purpose. FirouzAbadi et al. [19] developed a reduced order BEM model for solutions of sloshing problem in an off-center annular cylindrical tank, and the surface elevation and natural frequencies were analyzed. Based on the Galerkin method, Takahara et al. [56] analyzed the frequency response of the nonlinear liquid sloshing in a cylindrical tank with a coaxial barrel. Hasheminejad and Soleimani [29] developed a 3D hydrodynamic analysis for sloshing in doubly connected fluid domain of an upright elliptical cylindrical tank with an internal elliptical core barrel based on the linear water wave theory and separation of variables. Sanapala et al. [49] performed an experiment to study the interaction of an internal vertical pole and the dynamics of liquid sloshing in a 3D rectangular tank. Generally, all the aforementioned studies were about the sloshing in cylindrical or rectangular tanks with simply circular or elliptical cylindrical barrels (pipe or pole), and most of them only analyzed the natural frequencies and mode shapes. Being different from the above descriptions, one main intention of this paper is to perform a comprehensive study on the effects of geometric parameters, number and arrangement of barrels on liquid sloshing characteristics (involving the sloshing force, surface elevation and natural frequencies) in the cylindrical tanks. Meanwhile, as far as authors know, there is no report on modeling the liquid sloshing problem by using the isogeometric boundary element method (IGABEM) yet, which is a novel numerical method based on the concept of isogeometric analysis developed firstly by Hughes et al. [31], i.e. the numerical analysis approach (originally, the isogeometric finite element method) can exactly describe the geometry of problem domain generated by CAD software without any change. That is because both the numerical method and the CAD system use the same bases, and there is no approximation exists in the conversion from CAD geometry to the numerical model. The non-uniform rational B-spline (NURBS) bases were

incorporated into the FEM by Hughes et al. [31] for the reason NURBS bases were widely used in the CAD industry. Being different from the traditional polynomial bases, NURBS can efficiently handle the geometry with very few control points, which allows the isogeometric analysis to obtain more accurate solutions even with less computational effort. In view of the advantages of NURBS, some years later the isogeometric concept was integrated into the boundary element method (BEM) by Takahashi and Matsumoto [57], Simpson et al. [51] and Simpson et al. [50], and they termed it as the isogeometric boundary element method (IGABEM). Since then, the IGABEM has been widely developed by many researchers for various problems. For examples, Beer et al. [5] and Marussig et al. [42] used the IGABEM to simulate the underground excavations. Simpson et al. [52] extended their work in the simulation of acoustic phenomena by using the T-splines enhanced IGABEM. Heltai et al. [30] proposed a non-singular formulation of IGABEM for three dimensional Stokes flows. Bai et al. [3] employed the IGABEM to evaluate the effective elastic properties as well as the stress states of doubly periodic complex shaped inclusions. Wang et al. [62] proposed a nonsingular IGABEM for 3D elastostatic problems by coupling trimmed NURBS with the nonsingular BEM. Lian et al. [38] adopted the NURBS-based IGABEM for the 2D gradient-based shape optimization. Gong and Dong [25], Zhou et al. [67] and Wang et al. [59] investigated the 3D and 2D potential problems with the NURBS enhanced IGABEM. Coox et al. [14] applied an indirect IGABEM to solve acoustic problems within open-boundary domains. Peng et al. [45] and Sun et al. [54] used the IGABEM for the static fracture and crack propagation problems. Qu et al. [46] and Gong et al. [26] studied the steady state thermal conduction in heterogeneities with the IGABEM. Falini et al. [18] proposed an adaptive IGABEM, which was used to solve the crack problem on a slit. Recently, Chen et al. [11] and Sun et al. [54] extended the IGABEM to the three-dimensional acoustic problems. Meanwhile, many of the aforementioned researchers (such as [2,45,50,51,67]) pointed out that the IGABEM possessed more excellent performance in accuracy and convergence than the traditional BEM. This paper extends the IGABEM to the solutions of the Helmholtz and modified Helmholtz equations for the investigation of sloshing phenomena in cylindrical tanks with multiply connected domains, and it can be organized as follows. Section 2 gives the detailed derivations of the isogeometric boundary integral equations for liquid sloshing problem based on 2D Helmholtz and modified Helmholtz equations. In Section 3, the proposed IGABEM model is verified against the available results from some other researchers, and then effects of geometric parameters, number and arrangement of internal barrels on liquid sloshing characteristics (sloshing force, surface elevation and natural frequencies) in cylindrical tanks are studied in detail. Finally, some main conclusions are summarized in Section 4. 2. Basic formulations of liquid sloshing problem 2.1. Description of the basic theory Fig. 1 provides the schematic diagram of a partially-filled cylindrical tank subjected to a lateral excitation u(x), the origin of the Cartesian coordinate system oxyz is located at the center of the liquid free surface, and the internal barrel is considered to move together with the tank system. The height from the bottom of the tank to the liquid free surface is defined as H, and the radius of the tank’s cross-section is r. An ideal liquid without compressibility and viscosity together with the linear boundary condition at the liquid free surface is used in this paper. Under these definitions and assumptions, the solution of sloshing can be described as a function of velocity potential, which is given as [43] Φ(𝑥, 𝑦, 𝑧, 𝑡) = 𝜑(𝑥, 𝑦, 𝑧)𝑒−𝑖𝜔𝑡

(1)

where 𝜑√ denotes the velocity potential of the flow field, the imaginary unit 𝑖 = −1, and t as well as 𝜔 are the time and excitation frequency, respectively. It is clear from Eq. (1) that the item 𝜑 is expressed as the 149

J. Liu, Q. Zang and W. Ye et al.

Engineering Analysis with Boundary Elements 114 (2020) 148–165

𝜕 2 𝜑𝑚 𝜕 𝑥2

+

𝜕 2 𝜑𝑚 𝜕 𝑦2

− 𝑘2𝑚 𝜑𝑚 = 0

(𝑚 = 1, 2, … , ∞)

(9)

and then, the velocity potential 𝜑(x, y, z) can be found via Eq. (5) after 𝜑0 and 𝜑m being solved from Eqs. (8) and (9). However, the boundary condition at the tank wall is also needed in order to solve Eqs. (8) and (9), which is given as 𝜑(𝑥, 𝑦, 𝑧),𝒏 = −𝑖𝜔𝐴 cos 𝜃

on 𝑆𝑟

(10)

where Sr indicates the rigid boundaries involving the walls of the tank and the internal barrels, 𝜑(x, y, z), n is the normal derivative of 𝜑(x, y, z), 𝜃 denotes the angle from the normal vector n to the positive x-axis, and A represents the amplitude of the lateral excitation. One can note that there are no terms of 𝜑0 and 𝜑m in Eq. (10), which implies that it cannot be applied directly to solve Eqs. (8) and (9). Therefore, some transformations are proposed to be made to Eq. (10). It is clear that 𝜑(x, y, z), n in Eq. (10) can be expanded as 𝜕𝜑(𝑥, 𝑦, 𝑧) 𝜕𝜑(𝑥, 𝑦, 𝑧) 𝜕𝜑(𝑥, 𝑦, 𝑧) 𝑛𝑥 + 𝑛𝑦 + 𝑛𝑧 𝜕𝑥 𝜕𝑦 𝜕𝑧

𝜑(𝑥, 𝑦, 𝑧),𝒏 =

(11)

in which nx , ny and nz are components of n along the x, y and z-axis directions. For the cylindrical tank mounted with vertical barrel (barrels) as shown in Fig. 1, the value of nz is equal to zero at Sr , hence substituting Eq. (5) into Eq. (11) yields 𝜕𝜑(𝑥, 𝑦, 𝑧) 𝜕𝜑(𝑥, 𝑦, 𝑧) 𝑛𝑥 + 𝑛𝑦 𝜕𝑥 𝜕𝑦 ( ) cosh 𝑘0 (𝑧 + 𝐻) 𝜕 𝜑0 (𝑥, 𝑦) 𝜕 𝜑0 (𝑥, 𝑦) = 𝑛𝑥 + 𝑛𝑦 cosh 𝑘0 𝐻 𝜕𝑥 𝜕𝑦 ( ) ∞ ∑ 𝜕 𝜑𝑚 (𝑥, 𝑦) 𝜕 𝜑𝑚 (𝑥, 𝑦) + cos 𝑘𝑚 (𝑧 + 𝐻) 𝑛𝑥 + 𝑛𝑦 𝜕𝑥 𝜕𝑦 𝑚=1

𝜑(𝑥, 𝑦, 𝑧),𝒏 =

Fig. 1. Cylindrical tank with a barrel and the Cartesian coordinate system.

function of space coordinate only, which can be explained by the threedimensional Laplace equation ∇2 𝜑 =

𝜕2 𝜑 𝜕2 𝜑 𝜕2 𝜑 + + =0 𝜕 𝑥2 𝜕 𝑦2 𝜕 𝑧2

in Ω

𝜕 𝜑 (𝑥, 𝑦) cosh 𝑘0 (𝑧 + 𝐻) 𝜕 𝜑0 (𝑥, 𝑦) ∑ + cos 𝑘𝑚 (𝑧 + 𝐻) 𝑚 cosh 𝑘0 𝐻 𝜕 𝒏(𝑥, 𝑦) 𝜕 𝒏(𝑥, 𝑦) 𝑚=1 ∞

=

= −𝑖𝜔𝐴 cos 𝜃

(2)

(12) cosh 𝑘0 (𝑧+𝐻) cosh 𝑘0 𝐻

Due to the orthogonality relation between

and

in which Ω is the entire flow domain within the cylindrical tank. Meanwhile, the kinematic condition at the free surface based on the linear potential theory and the Bernoulli equation is explained as

cos km (z + H) for the eigenvalue problems, Eq. (12) can further be decoupled as follow

𝜕𝜑 𝜔2 − 𝜑=0 𝜕𝑧 𝑔

𝜕 𝜑0 (𝑥, 𝑦) = −𝑖𝜔𝐴𝜆0 cos 𝜃 𝜕 𝒏(𝑥, 𝑦)

on 𝑧 = 0

(3)

with g denoting the gravity acceleration. Additionally, the impermeable boundary condition at the tank bottom is given as 𝜕𝜑 =0 𝜕𝑧

on 𝑧 = −𝐻

𝜕 𝜑𝑚 (𝑥, 𝑦) = −𝑖𝜔𝐴𝜆𝑚 cos 𝜃 𝜕 𝒏(𝑥, 𝑦)

(4)

cosh 𝑘0 (𝑧 + 𝐻) ∑ + 𝜑𝑚 (𝑥, 𝑦) cos 𝑘𝑚 (𝑧 + 𝐻) cosh 𝑘0 𝐻 𝑚=1

(

) ( ) 0 cosh2 𝑘 (𝑧 + 𝐻) cosh 𝑘0 (𝑧 + 𝐻) 0 𝑑𝑧 ( ) 𝑑𝑧 ∕ ( ) ∫−𝐻 cosh 𝑘0 𝐻 ∫−𝐻 cosh2 𝑘 𝐻 0 ( ) 2 sinh 2𝑘0 𝐻 = ( ) sinh 2𝑘0 𝐻 + 2𝑘0 𝐻

(5)

(

in which m denotes the order of evanescent mode, and k0 and km are the wave numbers of the propagating and evanescent modes, respectively, which can be determined by solving the dispersion equations as follows

𝜔2 = −𝑔 𝑘𝑚 tan 𝑘𝑚 𝐻

𝜆𝑚 = =

(6) (𝑚 = 1, 2, … , ∞)

𝜕 𝑥2

+

𝜕 2 𝜑0 𝜕 𝑦2

+ 𝑘20 𝜑0 = 0

0

∫−𝐻

) ( cos𝑘𝑚 (𝑧 + 𝐻)𝑑𝑧 ∕ (

)

4 sin 𝑘𝑚 𝐻 ( ) sin 2𝑘𝑚 𝐻 + 2𝑘𝑚 𝐻

0

∫−𝐻

(15)

) cos 𝑘𝑚 (𝑧 + 𝐻)𝑑𝑧 2

(𝑚 = 1, 2, … , ∞)

(16)

Thereafter, 𝜑0 and 𝜑m can be determined by Eqs. (8) and (9), and then the velocity potential at any position within the flow domain can be calculated by Eq. (5), after which the liquid dynamic pressure p as well as the elevation of the liquid free surface 𝜂 can be acquired as

(7)

Meanwhile, 𝜑0 and 𝜑m in Eq. (5) are the propagating mode and evanescent mode potentials, which satisfy the two-dimensional Helmholtz and modified Helmholtz equations, respectively, which can be further written as [63] 𝜕 2 𝜑0

(14)

0

𝜆0 =

∞

𝜔2 = 𝑔 𝑘0 tanh 𝑘0 𝐻

(𝑚 = 1, 2, … , ∞)

with

By considering the boundary conditions mentioned above and introducing the corresponding eigen-functions with respect to the z coordinate, the term of 𝜑(x, y, z) in Eq. (1) can be rewritten as [40] 𝜑(𝑥, 𝑦, 𝑧) = 𝜑0 (𝑥, 𝑦)

(13)

𝑝(𝑥, 𝑦, 𝑧) = 𝑖𝜔𝜌𝜑 𝜂=

(8)

150

𝑖𝜔 𝜑 𝑔

(17) (18)

J. Liu, Q. Zang and W. Ye et al.

Engineering Analysis with Boundary Elements 114 (2020) 148–165

𝜑∗0 (𝑃 , 𝑄) =

𝑖 𝑌 (𝑘 𝑅) 4 0 0

𝜑∗𝑚 (𝑃 , 𝑄) =

1 ′ 𝑌 (𝑘 𝑅) 2𝜋 0 𝑚

(26) (𝑚 = 1, 2, … , ∞)

(27)

in which R indicates the distance from point P to Q as shown in Fig. 2. 𝑌0 together with Y′0 are described as 𝑌0 (𝑘0 𝑅) = 𝐼0 (𝑘0 𝑅) + 𝑖𝐵0 (𝑘0 𝑅)

(28)

( ) ( )2𝑙 ∞ ∑ 𝑘 𝑅 𝜓(𝑙) 𝑘𝑚 𝑅 𝑌0′ (𝑘𝑚 𝑅) = − ln 𝑚 + ℎ 𝐼0′ (𝑘𝑚 𝑅) + (𝑚 = 1, 2, … , ∞) 2 2 2 𝑙=1 (𝑙!) (29) with Fig. 2. Multiply connected domain with the source and field points.

𝜓(𝑗) = with 𝜌 denoting the liquid density. Finally, the total sloshing force on the whole tank system in the positive x-direction can be obtained as ( ) 𝑁 0 ∑ 𝐹𝑥 = 𝑓𝑤 + 𝑓𝑏𝑗 𝑑𝑧 (19) ∫−𝐻 𝑗

𝑓𝑏 =

2𝜋

∫0 ∫Γ𝑏

𝑝𝑤 𝑟 cos 𝜃𝑑𝜃

𝑞0∗ (𝑃 , 𝑄)=

(20)

(21)

|𝐹𝑥 |

2.2. Boundary integral equations for Helmholtz and modified Helmholtz equations Based on the weighted residual method and the Green’s third identity boundary integral equations for the Helmholtz and modified Helmholtz equations can be provided as (see [8,6]) ( ) 𝐶 𝑃 𝜑𝑃0 + 𝑞0∗ (𝑃 , 𝑄) − 𝜑∗0 (𝑃 , 𝑄)𝑞0𝑄 𝑑 Γ𝑟 = 0 (23) 𝜑𝑄 0 ∫Γ𝑟 (

) ∗ ∗ 𝑄 𝜑𝑄 𝑚 𝑞 𝑚 ( 𝑃 , 𝑄 ) − 𝜑 𝑚 ( 𝑃 , 𝑄 ) 𝑞 𝑚 𝑑 Γ𝑟 = 0

(𝑚 = 1, 2, … , ∞)

=−

] 𝜕𝑅 𝑘0 𝑖 𝑘 𝑖[ 𝜕𝑅 𝑌 (𝑘 𝑅) = − 0 𝐼1 (𝑘0 𝑅) + 𝑖𝐵1 (𝑘0 𝑅) 4 1 0 𝜕𝒏 4 𝜕𝒏 (31)

𝜕𝜑∗𝑚 (𝑃 , 𝑄) 𝜕𝒏 [

(𝑚 = 1, 2, … , ∞)

where I1 and I′1 are the first order Bessel function and modified Bessel function of the first kind, and B1 is the first order Bessel function of the second kind. When the distance from point P to Q tends to zero i.e. R → 0, the terms of 𝑌0 (𝑘𝑅) in Eq. (26) and Y′0 (kR) in Eq. (27) tend to (2/𝜋)ln R and − ln R, both of which show the logarithmic singularity, thus for the case P and Q are located at the same element, a particular Gauss quadrature method proposed by Guiggiani and Casalini [28] is adopted to treat the singularity in this paper.

with V denoting the volume of the liquid in the tank.

∫Γ𝑟

(30)

(32)

(22)

𝜌𝐴𝜔2 𝑉

𝐶 𝑃 𝜑𝑃𝑚 +

𝜕𝒏

𝑘 𝜕𝑅 = − 𝑚 𝑌 ′ 1 (𝑘𝑚 𝑅) 2𝜋 𝜕𝒏 ( ) 𝑘𝑚 𝑘𝑚 𝑅 1 =− + ln + 𝐶 𝐼 ′ 1 (𝑘𝑚 𝑅) 2𝜋 𝑘𝑚 𝑅 2 ( )2𝑛+1 ] ∞ 1 ∑ 𝜓(𝑙 + 1) + 𝜓(𝑙) 𝑘𝑚 𝑅 𝜕𝑅 − 2 𝑙=0 𝑙!(𝑙 + 1)! 2 𝜕𝒏

where Γb denotes the boundary of one internal barrel on the xy plane, and pw together with pb are the liquid dynamic pressures on the tank wall and internal barrel, respectively. Furthermore, a dimensionless factor of 𝜌A𝜔2 V is used in this paper to normalize the value of Fx as |𝐹𝑥 | =

𝜕𝜑∗0 (𝑃 , 𝑄)

𝑞𝑚∗ (𝑃 , 𝑄) = 𝑝𝑏 cos 𝜃𝑑 Γ𝑏

𝑗 = 1, 2, …

h in Eq. (29) stands for the Euler’s constant, I0 as well as 𝐼0′ in Eqs. (28) and (29) are the zero order Bessel function and modified Bessel function of the first kind, respectively, and B0 in Eq. (28) represents the zero order Bessel function of the second kind. The derivatives of 𝜑∗0 and 𝜑∗𝑚 with respective to n are defined by 𝑞0∗ and 𝑞𝑚∗ as described in Eqs. (23) and (24), which can be depicted as

in which N represents the total number of internal barrels, and fw as well as 𝑓𝑏𝑗 denote the dynamic sloshing forces acting on the unit vertical length of the tank wall and the jth internal barrel, respectively, which can be expressed as follow 𝑓𝑤 =

𝑗 ∑ 1 , 𝑙 𝑙=1

2.3. Isogeometric BEM model for liquid sloshing The philosophy of IGA-BEM is to approximate the field variables by the same basis functions used for representing the boundary geometry of the problem domain, and it is known that the NURBS is a widespread and thoroughly developed technology in CAD industry, therefore the NURBS is adopted in this paper for the simulation of the liquid sloshing in cylindrical tanks with multiply connected domain. The NURBS curve can be constructed as

(24)

where 𝑞0 and 𝑞𝑚 are the derivatives of 𝜑0 and 𝜑𝑚 , which can be called propagating mode and evanescent mode flux densities, respectively. P and Q are the source and field points as shown in Fig. 2, and Γr denotes the rigid boundaries involving the tank wall and the internal barrels. CP is a coefficient associated with the internal angle at point P, and for the two-dimensional problem it can be described as 𝛼 𝐶𝑃 = (25) 2𝜋

𝑪 (𝑢) =

𝑛 ∑ 𝑖=1

𝑅𝑖,𝜅 (𝑢)𝑷 𝑖

(33)

in which Pi is the coordinate vector of the ith control point, u indicates the knot defined in the parameter space, and the knot vector 𝚵 = {u1 ,u2 ,…, un + 𝜅 + 1 } is also an essential factor for defying the NURBS, n denotes the number of NURBS basis functions, and Ri,𝜅 (u) means the ith NURBS basis function with 𝜅 denoting the degree of NURBS curve, which is defined as

where 𝛼 is the internal angle at point P. In particular, for the cases P is outside or inside the flow domain 𝛼 = 0 or 𝛼 = 2𝜋. Meanwhile, 𝜑∗0 and 𝜑∗𝑚 in Eqs. (23) and (24) are the fundamental solutions of the two-dimensional Helmholtz and modified Helmholtz equations, which are expressed as

𝑅𝑖,𝜅 (𝑢) = ∑𝑛

𝑤𝑖 𝑁𝑖,𝜅 (𝑢)

𝑗=1

151

𝑤𝑗 𝑁𝑗,𝜅 (𝑢)

(34)

J. Liu, Q. Zang and W. Ye et al.

Engineering Analysis with Boundary Elements 114 (2020) 148–165

Fig. 3. NURBS curve and the associated NURBS basis functions. (a) NURBS curve, and (b) NURBS basis functions.

where 𝜑0 and 𝑞0 are the propagating mode potential and flux density, while 𝜑𝑚 and 𝑞𝑚 are the evanescent mode potential and flux density. The superscript e implies the index of the element, i is the index of the control point in the eth element. In order to perform the numerical integration conveniently along the boundary of the problem domain, the Gauss–Legendre rule is commonly used by many methods, and this paper is not an exception, which utilizes the local coordinate 𝜉 ∈ [ − 1, 1] as the parent element, and a Jacobian is also adopted to transform the aforementioned integral over a specified element into the parent element. Eq. (40) gives the formula of the Jacobian (see Simpson et al. [50] for the detailed derivations), which involves the transformation from physical space to parameter space and that from parameter space to parent space

where wi is the weight associated with the ith NURBS basis function, which is used to control the curvature of the NURBS curve, and Ni,𝜅 (u) indicates the ith B-spline basis function. Detailed definition of NURBS and its properties can be found in the existing literature published by Hughes et al. [31]. In order to make the concept of NURBS to be more intuitive, an example of a quadratic NURBS curve defined upon the knot vector of 𝚵={0, 0, 0, 1, 2, 3, 3, 4, 5, 6, 6, 6} and the associated NURBS basis functions are proposed in Fig. 3. It can be seen that, 9 control points correspond to 9 basis functions, and the values of the first and last knots 0 and 6 are repeated three times in 𝚵={0, 0, 0, 1, 2, 3, 3, 4, 5, 6, 6, 6}, while the value of knot 3 is repeated two times, these make the NURBS basis functions interpolatory at the repeated knots (see R1,2 , R5,2 , R9,2 in Fig. 3(b)), and accordingly the NURBS curve shown in Fig. 3(a) show geometric corners at control points 1, 5, 9 (1 coincides with 9). The partition of unity property is one of the main features of the NURBS basis. According to this property the parameter space can be divided into several small knot spans, in any of which there are at most 𝜅 + 1 non-zero NURBS basis functions, and the elements for IGABEM are defined on these knot spans. One can note that the definition of IGABEM element is quite different from that of the conventional BEM element. After the elements are defined, based on the partition of unity property of the NURBS basis functions, the geometry as well as the field variations of the domain boundary can be readily estimated as follows 𝒙( 𝑢 ) =

𝜅+1 ∑ 𝑗=1

𝑅𝑒𝑖,𝜅 (𝑢)𝑷 𝑖,𝑒

𝑖=1

𝜑𝑚 (𝑢) =

𝑞0 (𝑢) =

𝑞𝑚 (𝑢) =

𝑖=1

𝜅+1 ∑ 𝑖=1 𝜅+1 ∑ 𝑖=1

and 𝑖

𝜅+1 ∑ 𝑗=1

(36) =

𝜅+1 ∑

′ 𝑅𝑒𝑗,𝜅

(𝑚 = 1, 2, … , ∞)

𝑅𝑒𝑖,𝜅 (𝑢)𝑞0𝑖,𝑒

𝑅𝑒𝑖,𝜅 (𝑢)𝑞𝑚𝑖,𝑒

(40)

[

𝑒=1 𝑗=1

1

∫−1

[

1

∫−1

] 𝑞𝑚∗ (𝑖, 𝑄(𝜉))𝑅𝑒𝑗,𝜅 (𝜉)𝐽 𝑒 (𝜉)𝑑𝜉

𝜑𝑗,𝑒 𝑚

]

𝜑∗𝑚 (𝑖, 𝑄(𝜉))𝑅𝑒𝑗,𝜅 (𝜉)𝐽 𝑒 (𝜉)𝑑𝜉 𝑞𝑚𝑗,𝑒

(𝑚 = 1, 2, … , ∞)

(42)

in which e′ denotes the element where the source point i is located, 𝜉′ is the local coordinate with respect to point i, and the boundary of problem domain is divided into NE elements. Moreover, Je′ and Je are the Jacobian according to elements e′ and e, respectively. Furthermore, it should be noted that the boundaries of the problem domain are defined on different knot vectors, thus the loops over every boundary have to be introduced to the implementation of IGABEM for the multiply connected domain. For instance, the discretized equations of the multiply

(37)

(38)

(𝑚 = 1, 2, … , ∞)

𝑁𝐸 𝜅+1 ∑ ( ′ ) 𝑗,𝑒′ ∑ 𝜉 𝜑𝑚 +

𝑁𝐸 𝜅+1 ∑ ∑ 𝑒=1 𝑗=1

𝑅𝑒𝑖,𝜅 (𝑢)𝜑𝑖,𝑒 𝑚

𝜉 ∈ [−1, 1]

𝑒=1 𝑗=1

(35)

𝑅𝑒𝑖,𝜅 (𝑢)𝜑𝑖,𝑒 0

𝑑 Γ𝑟 𝑑 Γ𝑟 𝑑𝑢 = 𝑑𝜉 𝑑𝑢 𝑑𝜉

Considering that the coordinate in parameter space can be expressed as a function of 𝜉 as u(𝜉), the discretized equations can be achieved by substituting -(Eqs. (35)–(40) into Eqs. (23) and (24) as ] [ 𝜅+1 𝑁𝐸 𝜅+1 1 ∑ ′ ( ) 𝑗,𝑒′ ∑ ∑ 𝐶𝑖 𝑅𝑒𝑗,𝜅 𝜉 ′ 𝜑0 + 𝑞0∗ (𝑖, 𝑄(𝜉))𝑅𝑒𝑗,𝜅 (𝜉)𝐽 𝑒 (𝜉)𝑑𝜉 𝜑𝑗,𝑒 0 ∫−1 𝑗=1 𝑒=1 𝑗=1 ] [ 𝑁𝐸 𝜅+1 1 ∑ ∑ = 𝜑∗ (𝑖, 𝑄(𝜉))𝑅𝑒𝑗,𝜅 (𝜉)𝐽 𝑒 (𝜉)𝑑𝜉 𝑞0𝑗,𝑒 (41) ∫−1 0

𝐶

𝜅+1 ∑

𝜑0 (𝑢) =

𝐽 (𝜉) =

(39) 152

J. Liu, Q. Zang and W. Ye et al.

Engineering Analysis with Boundary Elements 114 (2020) 148–165

By performing Eqs. (43) and (44) at all the mentioned collocation points, the discretised equations in IGABEM for sloshing can be expressed as 𝑯 0 𝝋0 = 𝑩 0 𝒒 0 𝑯 𝑚 𝝋𝑚 = 𝑩 𝑚 𝒒 𝑚

connected domain with N boundaries can be deduced from Eqs. (41) and (42) as 𝐶

𝜅+1 ∑ 𝑗=1

=

′ 𝑅𝑒𝑗,𝜅

𝑁 ∑ 𝑁𝐸 𝜅+1 ∑ ( ′ ) 𝑗,𝑒′ ∑ 𝜉 𝜑0 +

[ 𝑁 ∑ 𝑁𝐸 𝜅+1 ∑ ∑ 𝑙=1 𝑒=1 𝑗=1

𝑙=1 𝑒=1 𝑗=1

1

∫−1

[

1

∫−1

𝑞0∗

] ( ) 𝑒,𝑙 𝑙 𝑒,𝑙 𝑖, 𝑄 (𝜉) 𝑅𝑗,𝜅 (𝜉)𝐽 (𝜉)𝑑𝜉 𝜑𝑗,𝑒,𝑙 0

] ( ) 𝑒,𝑙 𝜑∗0 𝑖, 𝑄𝑙 (𝜉) 𝑅𝑒,𝑙 ( 𝜉) 𝐽 ( 𝜉) 𝑑𝜉 𝑞0𝑗,𝑒,𝑙 𝑗,𝜅

(𝑚 = 1, 2, … , ∞)

(47)

in which H0 , B0 , Hm and Bm are the coefficient matrices, and 𝜑0 and q0 are the velocity potentials and flux densities of the propagating mode at all the control points. Meanwhile, 𝜑m together with qm are the velocity potentials and flux densities of the evanescent mode at all the control points, respectively. Finally, 𝜑0 and 𝜑m can be determined by solving Eqs. (46) and (47), and then the liquid dynamic sloshing force as well as the elevation of the liquid free surface can be achieved. It should be noted that m = 1, 2, …, 4 for Eq. (48) is enough to obtain high precision results, thus, m = 1, 2, …, 4 is adopted for all the numerical examples in this paper. A flowchart of the multiply connected IGABEM code is shown in Fig. 5 to explain the computer program for implementing the present methodology. It can be seen from Fig. 5 that the preprocessing processes here are much different from those of the standard BEM, the initial information for model reconstruction is read directly from the CAD software, and there is almost no error generated in this process. Meanwhile, connectives between the domain boundary, patches and elements should be built carefully, which directly affect the correctness and calculation efficiency of the present code. Moreover, the blocks involving “Hgsubmatrix_PGQ_0″, “HGsubmatrix_PGQ_m”, “HGsubmatrices_GLQ_0″ and “HGsubmatrices_GLQ_m” illustrate the core subroutines of the proposed algorithm, and “Hgsubmatrix_PGQ_0″ (for singular integral) as well as “HGsubmatrices_GLQ_0″ (for regular integral) are used for calculating the coefficient matrixes in Eq. (46), while “HGsubmatrix_PGQ_m” together with “HGsubmatrices_GLQ_m” are applied to evaluate the coef′ ′ ficient matrixes in Eq. (47). And 𝑯 𝑒0 , 𝑯 𝑒0 , 𝑩 𝑒0 , 𝑩 𝑒0 in Fig. 5 denote the

Fig. 4. Collocation points of BURBS curve.

𝑖

(46)

element coefficient matrixes according to Eq. (46), while 𝑯 𝑒𝑚 , 𝑯 𝑒𝑚 , 𝑩 𝑒𝑚 , 𝑩 𝑒𝑚 are those according to Eq. (47). ′

(43)

′

and 𝐶

𝑖

𝜅+1 ∑ 𝑗=1

=

′ 𝑅𝑒𝑗,𝜅

𝑁 ∑ 𝑁𝐸 𝜅+1 ∑ ( ′ ) 𝑗,𝑒′ ∑ 𝜉 𝜑𝑚 +

[ 𝑁 ∑ 𝑁𝐸 𝜅+1 ∑ ∑ 𝑙=1 𝑒=1 𝑗=1

𝑙=1 𝑒=1 𝑗=1

1

∫−1

[

1

∫−1

𝑞𝑚∗

] ( ) 𝑒,𝑙 𝑙 𝑒,𝑙 𝑖, 𝑄 (𝜉) 𝑅𝑗,𝜅 (𝜉)𝐽 (𝜉)𝑑𝜉 𝜑𝑗,𝑒,𝑙 𝑚

3. Numerical examples In this section, the accuracy of the proposed IGABEM model is firstly verified against the available results achieved by several other researchers, and then some numerical examples are carried out for the studies of liquid sloshing in circular cylindrical tanks with different types of barrels. The effects of geometric parameters, number, and arrangement of the internal barrels on the sloshing forces, natural frequencies, and free surface elevations are discussed in detail. In addition, the amplitude of the lateral excitation, liquid density and gravitational acceleration are fixed at A = 1, 𝜌 = 1 and g = 9.81 in this paper for simplicity.

]

( ) 𝑒,𝑙 𝑗,𝑒,𝑙 𝜑∗𝑚 𝑖, 𝑄𝑙 (𝜉) 𝑅𝑒,𝑙 𝑗,𝜅 (𝜉)𝐽 (𝜉)𝑑𝜉 𝑞𝑚

(𝑚 = 1, 2, … , ∞) (44)

Another obvious difference between IGABEM and BEM is the definition of the collocation point, and it is clearly found from Fig. 3(a) that most of the control points are not located on the boundary, which can not be directly applied to the collocation of IGABEM formulation as the conventional BEM does. Hence the additional definition for the collocation points is needed, and the Greville abscissae definition [28,55] is used in this paper, which can be expressed as 𝑢̂ 𝑖 =

𝑢𝑖+1 + 𝑢𝑖+2 + … 𝑢𝑖+𝜅 , 𝑖 = 1, 2, … , 𝑛 − 1 𝜅

3.1. Verification examples The existing data obtained by Cao [9] and Yue et al. [64] are used for the verification purpose of the present IGA-BEM model in this section. Fig. 6 sketches the cross section of the cylindrical tank with a coaxial barrel, in which the radii of the tank wall and the internal barrel are denoted as r (the same definition is used in the rest of this paper) and r1 , respectively. Meanwhile, the liquid depth is defined as H, which is not shown here, and it can be seen in Fig. 1 (no special explanation will be given in the rest of this paper). In Addition, 40 elements with 48 control points (the same number of collocation points) are adopted for the proposed verification examples, and Fig. 7 illustrates the initial IGABEM meshes of the flow domain in cylindrical tank with a coaxial barrel, for which the information vectors of the control points at the

(45)

where 𝑢̂ 𝑖 indicates the knot value of the ith collocation point, after the knots of collocation points are defined, and the collocation points in physical space can be completely determined by Eq. (35). Fig. 4 shows the collocation points of the NURBS curve described in Fig. 3(a), and one can find that the all the collocation points are located on the boundary. As mentioned above, the field variations of the domain boundary are approximated based on the control points, while the given variations are always the boundary conditions at the collocation points, thus the inverse calculations are needed to transform the boundary conditions to those at the control points. 153

J. Liu, Q. Zang and W. Ye et al.

Engineering Analysis with Boundary Elements 114 (2020) 148–165

Fig. 5. Flowchart of the multiply connected IGABEM program.

inner and outer boundaries are

and

⎡−𝑟1 ⎢ ⎢−𝑟1 ⎢ 0 ⎢ ⎢ 𝑟1 ⎢𝑟 ⎢ 1 ⎢ 𝑟1 ⎢ 0 ⎢ ⎢−𝑟1 ⎢−𝑟 ⎣ 1

⎡−𝑟 ⎢ ⎢−𝑟 ⎢0 ⎢ ⎢𝑟 ⎢𝑟 ⎢ ⎢𝑟 ⎢0 ⎢ ⎢−𝑟 ⎢−𝑟 ⎣

0 𝑟1 𝑟1 𝑟1 0 −𝑟1 −𝑟1 −𝑟1 0

1 ⎤ ⎥ 2∕2⎥ 1 ⎥ √ ⎥ 2∕2⎥ 1 ⎥ √ ⎥ 2∕2⎥ 1 ⎥ √ ⎥ 2∕2⎥ 1 ⎥⎦ √

(48)

154

0 −𝑟 −𝑟 −𝑟 0 𝑟 𝑟 𝑟 0

1 ⎤ ⎥ 2∕2⎥ 1 ⎥ √ ⎥ 2∕2⎥ 1 ⎥ √ ⎥ 2∕2⎥ 1 ⎥ √ ⎥ 2∕2⎥ 1 ⎥⎦ √

(49)

J. Liu, Q. Zang and W. Ye et al.

Engineering Analysis with Boundary Elements 114 (2020) 148–165

Fig. 6. Cross section of the cylindrical tank with a coaxial barrel. Fig. 8. Comparison of the present surface elevation at the side wall with the analytic solution.

Fig. 7. Coarse meshes of the flow domain with IGABEM.

Fig. 9. Comparison of the second mode natural sloshing frequencies in cylindrical tank with a coaxial barrel with the analytical solutions.

respectively, with the first two columns are the coordinates of the control points while components of the last column indicate the weights. Furthermore, the knot vectors for both the inner and outer boundaries are given as 𝚵 = {0, 0, 0, 1∕4, 1∕4, 1∕2, 1∕2, 3∕4, 3∕4, 1, 1, 1}

3.2. Cylindrical tank with an eccentric elliptical barrel The cross section dimension and the coordinate systems of the cylindrical tank with an eccentric elliptical barrel are depicted in Fig. 10, in which a and b denote the lengths of major and minor half axes of the eccentric barrel, respectively. The Cartesian coordinate system o′x′y′ is independent to oxy, the origin of o′x′y′ is located at the center of the internal barrel, and the positive x′ and y′ axes are associated with the major and minor half axes of the internal barrel as shown in Fig. 10. Moreover, e in Fig. 10 represents the distance between the two origins o and o′, and 𝛼 1 together with 𝛼 2 are the angles from o′x′ and oo′ to the ox direction, respectively. The radius of the barrel and the liquid depth are fixed at r = 1 and H = 2 for this example. Furthermore, the eccentricity of the barrel is defined as 𝑒̄ = 𝑒∕(𝑟 − 𝑎). Fig. 11 shows the coarse meshes of the flow domain within the present tank, and a total of 40 elements together with 48 control points are used for this model, the information vector for the outer boundary is the same as Eq. (50) and that for the

(50)

Fig. 8 plots the comparison of the present surface elevation at the side wall of the tank with the analytic solution given by Cao [9]. In which the geometric parameters are fixed at r1 = 0.0001, r2 = 1 and H = 2 according to the existing literature. The surface elevation is expressed as a function of non-dimensional frequency 𝜔/(g/H)1/2 . The curve denotes the results of Cao [9], and the dots are the present results. Excellent agreement can be readily observed from Fig. 8. The comparison of second mode natural sloshing frequencies Ω (defined as 𝜔/2𝜋) against the radius ratio r1 /r with the fixed liquid depth H = 1.5 is depicted in Fig. 9, where the lines are the results of the present model, and the dots indicate the solutions of Yue et al. [64] with the simplified formulas. Once again, good agreements for different size of the cylindrical tank verify the proposed model. 155

J. Liu, Q. Zang and W. Ye et al.

Engineering Analysis with Boundary Elements 114 (2020) 148–165

Fig. 10. Cross section of the cylindrical tank with an eccentric elliptical barrel. Fig. 12. Variation of the non-dimensional sloshing force versus 𝜒.

and knot vectors for both the inner and outer boundaries are the same as Eq. (50). By defining the ratios of b/r and a/r at b/r = a/r = 𝜒, the effect of the radius of a circular barrel on sloshing force on the tank system is investigated as shown in Fig. 12 with 𝑒̄ = 0, kr = 2 and 𝛼 1 = 𝛼 2 = 0. It can be observed from Fig. 12 that the non-dimensional liquid sloshing forces |Fx | acting on both the whole system and the tank wall firstly decrease rapidly to their minimums, respectively, and then increase monotonously. When the value of 𝜒 is equal to about 0.38, |Fx | achieves its minimum for the tank wall, and 𝜒 = 0.3 seems to be the optimum choice for the circular barrel to suppress the sloshing force acting on the whole system. By contrast, the sloshing force on the inner barrel first increases gradually to a peak value at about 𝜒 = 0.5 and then decreases until 𝜒 equals to about 0.9. Finally, when 𝜒 ≥ 0.9 a sudden increase of |Fx | on the inner barrel and the tank wall can be clearly identified from Fig. 12, this is probably because when the radius of the internal barrel extends to close to the tank wall, and the annular region between the barrel and the tank wall becomes extremely narrow, which leads to more complex interaction between the internal barrel and the liquid, as well as the liquid and the tank wall. Fig. 13 exhibits the variation of the non-dimensional sloshing force against the kr at different eccentricity 𝑒̄ with the radius ratios b/r = a/r = 0.2, 𝛼 2 = 0° and 90° are selected for this example. When 𝛼 2 = 0°, an evident increase in the peak value of the sloshing force with the increase of 𝑒̄ can be found from Fig. 13(a). Fig. 13(a) also illustrates that with increasing 𝑒̄, the first natural frequency increases monotonously. However, changing the value of 𝑒̄ has almost negligible effect on either the sloshing force or the natural sloshing force in the case of 𝛼 2 = 90° as shown in Fig. 13(b). The above discussion reveals that changing the position of the internal barrel along the excitation direction has a great influence on the sloshing frequency of liquid in the container. By fixing the geometric parameters at a/r = 0.5, 𝑒̄ = 0, 𝛼 1 = 𝛼 2 = 0°, the non-dimensional sloshing force |Fx | is plotted in Fig. 14 as a function of kr with different ratios of b/a. One can note from Fig. 14 that the maximum of the sloshing force decreases generally with the increase in the ratio of b/a except for the two cases of b/a= 0.3 and b/a = 0.7. Meanwhile, increasing the ratio of b/a yields significant decrease in the first natural frequency. Similar to Fig. 14, Fig. 15 illustrates the variation of the nondimensional sloshing force |Fx | versus the kr at different ratios of a/b with the fixed b/r = 0.5, 𝑒̄ = 0 and 𝛼 1 = 𝛼 2 = 0°. It can be seen from Fig. 15 that with increasing the ratio of a/b, the first natural sloshing frequency shows a slight decrease, and the peak value of |Fx | first

Fig. 11. Coarse meshes of the flow domain in the cylindrical tank with an eccentric barrel.

inner ellipse can be determined as ⎡−𝑟1 ⎢ ⎢−𝑟1 ⎢ 0 ⎢ ⎢ 𝑟1 ⎢𝑟 ⎢ 1 ⎢ 𝑟1 ⎢ 0 ⎢ ⎢−𝑟1 ⎢−𝑟 ⎣ 1 ⎡1 ⎢1 ⎢ ⎢1 ⎢1 ⎢ + ⎢1 ⎢1 ⎢ ⎢1 ⎢1 ⎢ ⎣1

0 𝑟1 𝑟1 𝑟1 0 −𝑟1 −𝑟1 −𝑟1 0 1 1 1 1 1 1 1 1 1

1 ⎤ ⎥ 2∕2⎥ 1 ⎥ √ ⎥ 2∕2⎥ ⎡ cos 𝛼1 1 ⎥ ⎢− sin 𝛼1 √ ⎥⎢ 2∕2⎥ ⎣ 0 1 ⎥ √ ⎥ 2∕2⎥ 1 ⎥⎦ √

1⎤ 1⎥⎥ 1⎥ 1⎥ ⎡𝑒 cos 𝛼2 ⎥ 1⎥ ⎢ 0 ⎢ 1⎥ ⎣ 0 ⎥ 1⎥ 1⎥ ⎥ 1⎦

sin 𝛼1 cos 𝛼1 0

0 𝑒 sin 𝛼2 0

0⎤ 0⎥ ⎥ 0⎦

0⎤ 0⎥ ⎥ 1⎦

(51)

156

J. Liu, Q. Zang and W. Ye et al.

Engineering Analysis with Boundary Elements 114 (2020) 148–165

Fig. 13. Variation of the non-dimensional sloshing force against the kr at different eccentricity 𝑒̄. (a) 𝛼 2 = 0, and (b) 𝛼 2 = 90°.

Fig. 14. Variation of the non-dimensional sloshing force |Fx | vs. the kr at different b/a.

Fig. 16. Variation of the non-dimensional sloshing force |Fx | versus eccentricity 𝑒̄ at different a/b.

decreases to its minimum and then increases with the further increase in a/b, which implies that there exists an optimal value of a/b for the internal barrel to evoke the smallest sloshing force. Furthermore, based on the comparison between Figs. 14 and 15, one can clearly discover that the first natural frequency is more sensitive to the size change of internal barrel along the oy direction (the orthogonal direction of external excitation). As expected, the barrel can be regard as a baffle when the ratio of a/r is small enough. Based on this hypothesis, Fig. 16 describes the variation of the non-dimensional total sloshing force |Fx | versus eccentricity 𝑒̄ at a/r = 0.001, 𝛼 1 = 𝛼 2 = 0° and kr = 1.9, the ratios of b/r = {0.4, 0.6, 0.8} are selected to investigate the influence of baffle width on sloshing force. Due to the symmetry of the tank system and the antisymmetry of the external excitation, all the curves in Fig. 16 are symmetric. When the barrel (baffle) is located at the center of the cylindrical tank (𝑒̄ = 0.0), the non-dimensional sloshing force reaches to its minimum. Meanwhile, the sloshing force decreases monotonously with the increase in b/r for the case |𝑒̄| < 0.5.

Fig. 15. Variation of the non-dimensional sloshing force |Fx | versus the kr at different a/b. 157

J. Liu, Q. Zang and W. Ye et al.

Engineering Analysis with Boundary Elements 114 (2020) 148–165

Fig. 17. Cross sections of the cylindrical tanks with different number of circler cylindrical barrels. (a) Three barrels, (b) five barrels, (c) seven barrels, and (d) nine barrels.

Fig. 18. Variation of non-dimensional sloshing force acting on different parts of the tank system against kr with different number of barrels. (a) Total sloshing force on the whole system, and (b) sloshing force on the tank wall.

3.3. Cylindrical tanks with multiple barrels

numbered from 1 to n (n = 3, 5, 7, 9) with n denoting the total number of barrels, and the radii of these barrels are defied as 𝑟𝑖 , 𝑖 = 1, 2, … , 𝑛. As shown in Fig. 17, the barrels 1 are located at the center of the tanks, and barrels 2, 3 in Fig. 17(a), barrels 3, 5 in Fig. 17(b), barrels 3, 6 in Fig. 17(c), barrels 3, 8 in Fig. 17(d) are located on the y axis of these tanks, other barrels are evenly distributed between them along the circumferential direction. In order to ensure the cross sections of barrels in each tank having the same total area, the radii of barrels satisfy the

In this section, cylindrical tanks with multiple barrels are considered to investigate how the number and arrangement of barrels on liquid sloshing force, natural frequency and surface elevation. Fig. 17 depicts the cross sections of the cylindrical tanks mounted with different number of barrels, and the distribution and dimensions of the barrels are clearly described in these figures, in which the barrels are 158

J. Liu, Q. Zang and W. Ye et al.

Engineering Analysis with Boundary Elements 114 (2020) 148–165

Fig. 19. Surface elevation along the wall of the cylindrical tank with different number of barrels at kr = 1.2.

Fig. 21. Variation of non-dimensional sloshing force acting on different parts of the tank system against 𝛽. .

Fig. 20. Cross section of the cylindrical tank with three circler barrels. Fig. 22. Dimensionless surface elevation along the wall of the cylindrical tank with three circler cylindrical barrels at different 𝛽.

following relationship √ 𝑟𝑖 = 𝑟∕2 𝑛, 𝑖 = 1, 2, … , 𝑛

(52) with xi and yi denoting the coordinates of the center of the ith circle. While knot vectors for all the inner and outer boundaries are the same as Eq. (50). Fig. 18 exhibits the variation of the non-dimensional sloshing forces acting on the whole tank system and the tank wall against kr with different number of barrels. It is obvious that the peak values of sloshing forces acting on both the whole tank system and the tank wall increase with the increase in n. Meanwhile, as the number of barrels increases, the natural frequency increases gradually, and the increase rate of natural frequency becomes much small when n ≥ 5, from which it can be inferred that when the number of barrels increases to a certain one, and the frequency of the structure will not be affected by the further increase of n. The non-dimensional surface profiles along the tank wall are depicted in Fig. 19 with the fixed normalized wave number kr = 1.2. All the lines are of the obviously symmetric character due to the symmetry geometries of the tanks as shown in Fig. 17. It is clearly found that the curve shape for tank with three internal barrels is much different from those for other cases. This is because the phase variation occurs when the value of wave number kr is near about 1.2 for the case n = 3 as shown in Fig. 18(a). Fig. 19 further shows that as the number of barrels

Moreover, the radii of the tank wall and liquid depth are the same as those in Section 3.1. The distance from the origin o to the i-th (i ≠ 1) barrel is defined as e, whose value is equal to r/3 as described in Fig. 17. 𝛼 is the angle between the vector of a certain point at the tank wall and the ox direction. 20 elements comprised 24 control points for each circular boundary (involving all the internal barrels and the tank wall) are used for this model, the information vector for the outer boundary is the same as Eq. (50) and those for the ith inner circular barrel can be determined as ⎡− 𝑟 𝑖 ⎢ ⎢− 𝑟 𝑖 ⎢ 0 ⎢ ⎢ 𝑟𝑖 ⎢𝑟 ⎢ 𝑖 ⎢ 𝑟𝑖 ⎢ 0 ⎢ ⎢− 𝑟 𝑖 ⎢− 𝑟 ⎣ 𝑖

0 𝑟𝑖 𝑟𝑖 𝑟𝑖 0 −𝑟𝑖 −𝑟𝑖 −𝑟𝑖 0

1 ⎤ ⎡1 ⎥ 2∕2⎥ ⎢1 ⎢ 1 ⎥ ⎢1 √ ⎥ ⎢ 1 2∕2⎥ ⎢ 1 ⎥ + ⎢1 √ ⎥ 2∕2⎥ ⎢1 ⎢ 1 ⎥ ⎢1 √ ⎥ ⎢ 2∕2⎥ ⎢1 1 ⎥⎦ ⎣1 √

1 1 1 1 1 1 1 1 1

1⎤ 1⎥⎥ 1⎥ 1⎥ ⎡ 𝑥 𝑖 ⎥ 1⎥ ⎢ 0 ⎢ 1⎥ ⎣ 0 ⎥ 1⎥ 1⎥ ⎥ 1⎦

0 𝑦𝑖 0

0⎤ 0⎥ ⎥ 0⎦

(53)

159

J. Liu, Q. Zang and W. Ye et al.

Engineering Analysis with Boundary Elements 114 (2020) 148–165

Fig. 23. Contours of the dimensionless free surface elevation 𝜂/Ain tanks with three barrels at different 𝛽. (a) 𝛽 = 90°, (b) 𝛽 = 180°, and (c) 𝛽 = 270°. Fig. 24. Cross sections of the cylindrical tanks with internal barrels in different shapes. (a) Elliptical barrel, (b) equilateral triangular barrel, (c) square barrel, and (d) circular barrel.

increases, the sloshing wave tends to more stable. Based on Figs. 14, 16, 18 and 19 taking the cost of manufacture into account, three circular cylindrical barrels are enough to suppress the sloshing force and surface elevation. A cylindrical tank with three circler barrels is plotted in Fig. 20, in which the radii of these internal barrels are ri = r/5, i = 1, 2, 3, o2 together with o3 are centers of barrel 2 and barrel 3, o2 is located at the negative direction of ox axis. And the distances from the origin o to both o2 and o3 is e = 0.6. The definition of 𝛼 and the value of r are the same as those in Fig. 17. The liquid depth is still H = 2, and 𝛽 in Fig. 20 is the angle from vector oo2 to vector oo3 . 20 elements with 24 control points for each circular boundary are employed for this example.

The results of sloshing forces acting on different parts of the tank system versus 𝛽 is depicted in Fig. 21 at kr = 1.2. From Fig. 21, one can see that the sloshing force acting on the tank wall is much larger than that on the inner barrels, and the curve for the whole system is between the two mentioned above, which implies that the direction of the sloshing force on the tank wall is opposite to that on the inner barrels. Meanwhile, for all the three cases, when the values of 𝛽 is equal to 90° or 270°, the sloshing forces reach their maximum, and 𝛽 = 180° seems to be the optimum choice for the arrangement of the three barrels. Fig. 22 shows the curves of the dimensionless surface elevation along the wall of the cylindrical tank with three circler cylindrical barrels at kr = 1.2, and three different values of 𝛽 (90°, 180° and 270°) is 160

J. Liu, Q. Zang and W. Ye et al.

Engineering Analysis with Boundary Elements 114 (2020) 148–165

Fig. 25. Coarse meshes of the flow domains in the cylindrical tanks with internal barrels in different shapes. (a) Elliptical barrel, (b) triangular barrel, (c) square barrel, and (d) circular barrel.

√ √ parameters are fixed at 𝑎 = 0.8 2∕𝜋, 𝑏 = 0.4 2∕𝜋, L = 0.8 for square √ √ 4 barrel, 𝐿 = 1.6∕ 3 for triangular barrel and 𝑟1 = 0.8∕ 𝜋 to ensure the cross section of each barrel has the same area, which ensures that the volume of liquids in all tanks considered in this section is equal. Moreover, it should be noted that all the centers of the mentioned barrels are located at the centers of the tanks and the minor axis of the elliptical barrel, and one midline of the triangular barrel and one midpoint of the square are located at the ox axes. The definition of 𝛼 and the value of r and H are the same as those in Section 3.2. Meanwhile, Fig. 25 shows the initial meshes of the flow domain within the proposed tanks, and 20 elements together with 24 control points are used for each boundary except for the triangular one (15 elements together with 18 control points). The information vectors of the control points for all the outer boundaries considered in this example are the same as Eq. (49), and information vectors as well as the knot vectors for the internal boundaries can be expressed as follows

considered. As expected, the surface profile for the case of 𝛽 = 180° is the gentlest one compared with the other two cases, which implies that the three barrels arranged along the direction of excitation may be the optimal barrels arrangement for suppressing the sloshing behaviors. Furthermore, from the comparison between the curves for cases of 𝛽 = 90° and 𝛽 = 270°, one can find out that the curve shifts to the right when 𝛽 = 90° and shifts to the opposite direction when 𝛽 = 270°, which suggests that the internal barrel performs significant influence on the sloshing response of the adjacent liquids. According to Fig. 22, the contours of the dimensionless free surface elevation 𝜂/A are plotted in Fig. 23, which can further validate the above viewpoint. Furthermore, an interesting phenomenon can be observed from Fig. 23 that the liquid flow between two barrels tends to be more rapid if the barrels are arranged orthogonal to the direction of external excitation, whereas the contrary trend can be seen if the barrels are arranged along the direction of external excitation. 3.4. Cylindrical tanks with internal barrels in different shapes

⎡− 𝑏 ⎢ ⎢− 𝑏 ⎢0 ⎢ ⎢𝑏 ⎢𝑏 ⎢ ⎢𝑏 ⎢0 ⎢ ⎢− 𝑏 ⎢− 𝑏 ⎣

In this section, four barrels with different geometries (i.e. elliptical barrel, equilateral triangular barrel, square barrel and circular barrel) are utilized to study the effects of barrel shapes on sloshing response in cylindrical tanks. Cross sections of the cylindrical tanks with different internal barrels are plotted in Fig. 24, in which a and b are the lengths of major and minor half axes of the elliptical barrel as shown in Fig. 24(a), L denotes the side length of the square barrel or the triangular barrel as shown in Fig. 24(b) and (c), and r1 is the radius of the circular barrel as shown in Fig. 24(d). For comparison purpose, the aforementioned 161

0 𝑎 𝑎 𝑎 0 −𝑎 −𝑎 −𝑎 0

1 ⎤ ⎥ 2∕2⎥ 1 ⎥ √ ⎥ 2∕2⎥ 1 ⎥ √ ⎥ 2∕2⎥ 1 ⎥ √ ⎥ 2∕2⎥ 1 ⎥⎦ √

(54)

J. Liu, Q. Zang and W. Ye et al.

Engineering Analysis with Boundary Elements 114 (2020) 148–165

and

for the triangular barrel, and

𝚵 = {0, 0, 0, 1∕4, 1∕4, 1∕2, 1∕2, 3∕4, 3∕4, 1, 1, 1}

⎡−𝐿∕2 ⎢−𝐿∕2 ⎢ ⎢−𝐿∕2 ⎢ 0 ⎢ ⎢ 𝐿∕2 ⎢ 𝐿∕2 ⎢ ⎢ 𝐿∕2 ⎢ 0 ⎢ ⎣−𝐿∕2

(55)

for the elliptical barrel, √ ⎡ − 3𝐿∕3 ⎢ √ ⎢− 3𝐿∕12 ⎢ √ 3𝐿∕6 ⎢ ⎢ √ 3𝐿∕6 ⎢ ⎢ √ 3𝐿∕6 ⎢ ⎢ √ ⎢− 3𝐿∕12 ⎢ √ ⎣ − 3𝐿∕3

0 𝐿∕4 𝐿∕2 0 −𝐿∕2 −𝐿∕4 0

1⎤ ⎥ 1⎥ ⎥ 1⎥ ⎥ 1⎥ ⎥ 1⎥ ⎥ 1⎥ ⎥ 1⎦

1⎤ 1⎥⎥ 1⎥ 1⎥ ⎥ 1⎥ 1⎥ ⎥ 1⎥ 1⎥ ⎥ 1⎦

(58)

and (56)

𝚵 = {0, 0, 0, 1∕4, 1∕4, 1∕2, 1∕2, 3∕4, 3∕4, 1, 1, 1}

(59)

for the square barrel, and those for the circular barrel are the same as Eqs. (48) and (49). It can be seen from -(Eqs. (56)–(59) that, all the weights of control points for the triangular and square boundaries are 1, and the values of knots in 𝚵 are repeated three (the first and the last knots) or two (middle members of 𝚵) times, these factors together with the coordinates of the control points determine whether the corners of geometric shapes can be displayed. Actually, this initial information can be directly obtained from the CAD software, and there

and 𝚵 = {0, 0, 0, 1∕3, 1∕3, 2∕3, 2∕3, 1, 1, 1}

−𝐿∕2 0 𝐿∕2 𝐿∕2 𝐿∕2 0 −𝐿∕2 −𝐿∕2 −𝐿∕2

(57)

Fig. 26. Variation of non-dimensional sloshing force acting on different parts of the tank system against kr for cylindrical tanks with internal barrels in different shapes. (a) Total sloshing force on the whole system, (b) sloshing force on the tank wall, and (c) sloshing force on the internal barrel. 162

J. Liu, Q. Zang and W. Ye et al.

Engineering Analysis with Boundary Elements 114 (2020) 148–165

is no need to distinguish the geometries of the problem domains during the numerical analysis. Fig. 26 exhibits the variation of the non-dimensional sloshing forces acting on the whole tank system, the tank wall and the inner barrel against kr. It can be seen for the whole tank system and the tank wall the sloshing forces for different tanks show similar trends, i.e. the maximum of sloshing force for the tank with triangular barrel is slightly smaller than those for other tanks. Little discrepancies can be observed between the non-dimensional sloshing forces for square and circular barrels. The sloshing natural frequency of liquid in tank with elliptical barrel is smallest. The second mode resonance phenomenon for the tank with triangular barrel appears earlier than other cases. From Fig. 26(c), one can note that the peak value of sloshing force acting on the elliptical barrel is larger than those on the square and circular barrels, and the sloshing force on the triangular barrel is the least one, from which it can be inferred that the magnitude of liquid sloshing may be related to the frontal area of the internal barrel. Considering that the external excitation frequency encountered in practical engineering is usually smaller than the first-order natural sloshing frequency, thus it can be seen from Fig. 26 that the tank with circular barrel will be subjected smaller dynamic force than those equipped with other kinds of barrels. Fig. 27 shows the dimensionless surface elevations 𝜂/A along the wall of the cylindrical tanks with different barrels at kr = 1.0. It is clear from Fig. 27 that the motion of liquid in the tank with a circular barrel

Fig. 27. Dimensionless surface elevations along the wall of tanks with internal barrels in different shapes.

Fig. 28. Contours of the dimensionless free surface elevation 𝜂/A in tanks with internal barrels in different shapes at kr = 1.0. (a) Elliptical barrel, (b) triangular barrel, (c) square barrel, and (d) circular barrel.

163

J. Liu, Q. Zang and W. Ye et al.

Engineering Analysis with Boundary Elements 114 (2020) 148–165

References

is the gentlest one, and the maximum of the surface elevation in the tank with an elliptical barrel is relatively large by contrast. This confirms again the tank with circular barrel being the optimal choice for suppressing liquid sloshing. Finally, the contours of the dimensionless free surface elevation 𝜂/A in the tanks mentioned above are illustrated in Fig. 28 with the normalized wave number fixed at kr = 1.0, which further demonstrates the conclusions obtained from Figs. 26 and 27. In addition, these figures also show that liquid motions at both the windward and leeward sides of the internal barrels are much gentle than those in other regions.

[1] Abramson HN. The dynamic behavior of liquids in moving containers with applications to space vehicle technology. NASA SP; 1966. p. 106. [2] An Z, Yu T, Bui TQ, Wang C, Trinh NA. Implementation of isogeometric boundary element method for 2-D steady heat transfer analysis. Adv Eng Softw 2018;116:36–49. [3] Bai Y, Dong CY, Liu ZY. Effective elastic properties and stress states of doubly periodic array of inclusions with complex shapes by isogeometric boundary element method. Compos Struct 2015;128:54–69. [4] Bazazzadeh S, Shojaei A, Zaccariotto M, Galvanetto U. Application of the peridynamic differential operator to the solution of sloshing problems in tanks. Eng Comput (Swansea) 2019;36(1):45–83. [5] Beer G, Marussig B, Duenser C. Isogeometric boundary element method for the simulation of underground excavations. Géotechnique Lett 2013;3(3):108–11. [6] Blyth MG, Pozrikidis C. A comparative study of the boundary and finite element methods for the Helmholtz equation in two dimensions. Eng Anal Bound Elements 2007:35–49. [7] Bouabidi A, Driss Z, Abid MS. Study of hydrostatic pump created under liquid sloshing in a rectangular tank subjected to external excitation. Int J Appl Mech 2016;08(01):1650004. [8] Brebbia CA, Dominguez J. Boundary elements: an introductory course[M]. WIT press; 1994. [9] Cao F. The application of scaled boundary finite element method in potential flow theory Doctor thesis. Dalian, Liaoning, China: Dalian University of Technology; 2009. [10] Caron PA, Cruchaga MA, Larreteguy AE. Study of 3D sloshing in a vertical cylindrical tank. Phys Fluids 2018;30(8):082112. [11] Chen LL, Lian H, Liu Z, Chen HB. Structural shape optimization of three dimensional acoustic problems with isogeometric boundary element methods. Comput Methods Appl Mech Eng 2019;355:926–51. [12] Chen Y, Hwang W, Ko C. Sloshing behaviours of rectangular and cylindrical liquid tanks subjected to harmonic and seismic excitations. Earthquake Engng Struct 2007(36):1701–17. [13] Chen Y, Xue M. Numerical simulation of liquid sloshing with different filling levels using Openfoam and experimental validation. Water (Basel) 2018;10(12):1752. [14] Coox L, Atak O, Vandepitte D, Desmet W. An isogeometric indirect boundary element method for solving acoustic problems in open-boundary domains. Comput Methods Appl Mech Eng 2017;316:186–208. [15] Demirel E, Aral M. Liquid sloshing damping in an accelerated tank using a novel slot-baffle design. Water (Basel) 2018;10(11):1565. [16] Drake KR. The effect of internal pipes on the fundamental frequency of liquid sloshing in a circular tank. Appl Ocean Res 1999;21(3):133–43. [17] Ebrahimian M, Noorian MA, Haddadpour H. A successive boundary element model for investigation of sloshing frequencies in axisymmetric multi baffled containers. Eng Anal Bound Elem 2013;37(2):383–92. [18] Falini A, Giannelli C, Kanduč T, Sampoli ML, Sestini A. An adaptive iga-bem with hierarchical B-splines based on quasi-interpolation quadrature schemes. Int J Numer Methods Eng 2018;117(10):1038–58. [19] Firouz-Abadi RD, Haddadpour H, Ghasemi M. Reduced order modeling of liquid sloshing in 3D tanks using boundary element method. Eng Anal Bound Elem 2009;33(6):750–61. [20] Frosina E, Senatore A, Andreozzi A, Fortunato F, Giliberti P. Experimental and numerical analyses of the sloshing in a fuel tank. Energies 2018;11(3):682. [21] Fujita K, Ito T, Okada K. "Seismic response of liquid sloshing in the annular region formed by coaxial circular cylinders. Eng Comput (Swansea) 1985;2(4):299–306. [22] Gavrilyuk I, Lukovsky I, Trotsenko Y, Timokha A. Sloshing in a vertical circular cylindrical tank with an annular baffle. part 1. linear fundamental solutions. J Eng Math 2006;54(1):71–88. [23] Gedikli A, Ergüven ME. Seismic analysis of a liquid storage tank with a baffle. J Sound Vib 1999;223(1):141–55. [24] Gnitko V, Naumemko Y, Strelnikova E. Low frequency sloshing analysis of cylindrical containers with flat and conical baffles. Int J Appl Mech Eng 2017;22(4):867–81. [25] Gong YP, Dong CY. An isogeometric boundary element method using adaptive integral method for 3D potential problems. J Comput Appl Math 2017;319:141–58. [26] Gong YP, Dong CY, Qu XY. An adaptive isogeometric boundary element method for predicting the effective thermal conductivity of steady state heterogeneity. Adv Eng Softw 2018;119:103–15. [27] Grotle ELG, Bihs H, æsøy V, Pedersen E. Computational fluid dynamics simulations of nonlinear sloshing in a rotating rectangular tank using the level set method. J Offshore Mech Arctic Eng 2018;140(12). [28] Guiggiani M, Casalini P. Direct computation of cauchy principal value integrals in advanced boundary elements. Int J Numer Methods Eng 1987;24:1711–20. [29] Hasheminejad SM, H S. Linear solution for liquid sloshing in an upright elliptical cylindrical container with an eccentric core barrel. J Eng Mech 2017;143(9). [30] Heltai L, Arroyo M, Desimone A. Nonsingular isogeometric boundary element method for Stokes flows in 3D. Comput Methods Appl Mech Eng 2014;268:514–39. [31] Hughes TJR, Cottrell JA, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Comput Methods Appl Mech Eng 2005;194:4135–95. [32] Jena D, Chandra BK. Violent sloshing and wave impact in a seismically excited liquid-filled tank: meshfree particle approach. J Eng Mech 2017;144(3). [33] Jeon SH, Seo MW, Cho YU, Park WG, Jeong WB. Sloshing characteristics of an annular cylindrical tuned liquid damper for spar-type floating offshore wind turbine. Struct Eng Mech 2013;47(3):331–43. [34] Jeong K. Dynamics of a concentrically or eccentrically submerged circular cylindrical shell in a fluid-filled container. J Sound Vib 1999;224(4):709–32.

4. Conclusion This paper developed an NURBS based isogeometric boundary element method (IGABEM) for the liquid sloshing problem in multiply connected domain. The hypothesis of inviscid, incompressible liquid together with the linear boundary condition at the free surface are employed for solving the present problem. Meanwhile, the liquid sloshing in three-dimensional tank is reduced to 2D problem by taking the constant cross section of the tank into account. The proposed IGABEM model is verified against the available results from some existing literatures, good agreement is achieved, and then the effects of geometric parameters, number and arrangement of internal barrels on the sloshing force, surface elevation and natural frequencies of liquid in cylindrical tanks are discussed in detail. In the light of the analysis in Section 3, the tank with three circular barrels arranged along the excitation direction is expected to be the optimal barrel arrangement for suppressing liquid sloshing. Certainly, if some specified parts of the tank need carefully protection, appropriate design can be made according to the following conclusions: (1) Circular barrel with a reasonable radius can significantly suppress the sloshing force acting on the whole tank system. (2) Changing the position of the internal barrel along the excitation direction has great effect on the natural frequency of liquid in cylindrical tanks. (3) The first natural frequency is more sensitive to the size change of internal barrel along the orthogonal direction of the external excitation. (4) The peak values of sloshing forces acting on both the whole tank system and the tank wall increase with the increase in number of barrels. (5) As the number of barrels increases, the natural frequency increases gradually, and the increase rate of natural frequency becomes much small when the number of barrels increases to a relatively large one. (6) Increasing the number of barrels can obviously suppress the sloshing elevation. (7) The internal barrel performs evident influence on the sloshing response of the adjacent liquids, and the liquid flow between two barrels tends to be more rapid if the barrels are arranged orthogonal to the direction of external excitation. (8) Liquid sloshing motions at both the windward and leeward sides of the internal barrels are much gentle than those in other regions. (9) Triangular barrel has obvious effects on reducing the liquid sloshing force, and the circular barrel can make more contribution to suppress the free surface elevation. Acknowledgments This research was supported by Grant 51779033 from the National Natural Science Substrate of China, and Grant DUT18LK16 from the fundamental research funds for the Central Universities for which the authors are grateful. 164

J. Liu, Q. Zang and W. Ye et al.

Engineering Analysis with Boundary Elements 114 (2020) 148–165

[35] Kabiri MM, Nikoomanesh MR, Danesh PN, Goudarzi MA. Numerical and experimental evaluation of sloshing wave force caused by dynamic loads in liquid tanks. J Fluids Eng 2019;141(11). [36] Kargbo O, Xue M, Zheng J. Multiphase sloshing and interfacial wave interaction with a baffle and a submersed block. J Fluids Eng Trans Asme 2019;141(7). [37] Li X, Song C, Zhou G, Wei C, Lu M. Experimental and numerical studies on sloshing dynamics of PCS water tank of nuclear island building. Sci Technol Nuclear Installat 2018;2018:1–13. [38] Lian H, Kerfriden P, Bordas SPA. Implementation of regularized isogeometric boundary element methods for gradient-based shape optimization in two-dimensional linear elasticity. Int J Numer Meth Engng 2016;106(12):972–1017. [39] Liu Z, Feng Y, Lei G, Li Y. Fluid sloshing dynamic performance in a liquid hydrogen tank. Int J Hydrogen Energy 2019;44(26):13885–94. [40] Lin G, Liu J, Li J, Hu Z. A scaled boundary finite element approach for sloshing analysis of liquid storage tanks. Eng Anal Bound Elem 2015;56:70–80. [41] Lu D, Yin T, Liu H. Experimental and numerical investigation of sloshing behavior in annular region separated by several cylinders related to fast reactor design. Nuclear Eng Design 2018;339:235–43. [42] Marussig B, Beer G, Duenser C. Isogeometric boundary element method for the simulation in tunneling. AMM 2014;553:495–500. [43] Miao G. Analytical solutions for the sloshing loading on circular cylindrical liquid tanks with interior semi-porous barriers. J Hydrodyn 2001;32-:39. [44] Moradi R, Behnamfar F, Hashemi S. Mechanical model for cylindrical flexible concrete tanks undergoing lateral excitation. Soil Dyn Earthquake Eng 2018;106:148–62. [45] Peng X, Atroshchenko E, Kerfriden P, Bordas SPA. Isogeometric boundary element methods for three dimensional static fracture and fatigue crack growth. Comput Methods Appl Mech Eng 2017;316:151–85. [46] Qu XY, Dong CY, Bai Y, Gong YP. Isogeometric boundary element method for calculating effective property of steady state thermal conduction in 2D heterogeneities with a homogeneous interphase. J Comput Appl Math 2018;343:124–38. [47] Radnić J, Grgić N, Sunara Kusić M, Harapin A. Shake table testing of an open rectangular water tank with water sloshing. J Fluids Struct 2018;81:97–115. [48] Raynovskyy IA, Timokha AN. Damped steady-state resonant sloshing in a circular base container. Fluid Dyn Res 2018;50(4):45502. [49] Sanapala VS, Sajish SD, Velusamy K, Ravisankar A, Patnaik BSV. An experimental investigation on the dynamics of liquid sloshing in a rectangular tank and its interaction with an internal vertical pole. J Sound Vib 2019;449:43–63. [50] Simpson RN, Bordas SPA, Lian H, Trevelyan J. An isogeometric boundary element method for elastostatic analysis: 2D implementation aspects. Comput Struct 2013;118:2–12.

[51] Simpson RN, Bordas SPA, Trevelyan J, Rabczuk T. A two-dimensional isogeometric boundary element method for elastostatic analysis. Comput Methods Appl Mech Eng 2012;209-212:87–100. [52] Simpson RN, Scott MA, Taus M, Thomas DC, Lian H. Acoustic isogeometric boundary element analysis. Comput Methods Appl Mech Eng 2014;269:265–90. [53] Su Y, Liu Z, Gao Z. Shallow-water sloshing motions in rectangular tank in general motions based on boussinesq-type equations. J Hydrodyn 2018;30(5):958–61. [54] Sun FL, Dong CY, Yang HS. Isogeometric boundary element method for crack propagation based on Bézier extraction of Nurbs. Eng Anal Bound Elem 2019;99:76–88. [55] Greville TNE. Numerical procedures for interpolation by spline functions. J Soc Ind Appl Math Ser B Numer Anal 1964;1(1):53–68. [56] Takahara H, Kimura K. Frequency response of sloshing in an annular cylindrical tank subjected to pitching excitation. J Sound Vib 2012;331(13):3199–212. [57] Takahashi T, Matsumoto T. An application of fast multipole method to isogeometric boundary element method for Laplace equation in two dimensions. Eng Anal Bound Elem 2012;36(12):1766–75. [58] Timokha A. Analytically approximate natural sloshing modes and frequencies for an upright circular container with poles. J Eng Math 2016;101(1):47–54. [59] Wang J, Wang C, Liu J. Sloshing reduction in a pitching circular cylindrical container by multiple rigid annular baffles. Ocean Eng 2019;171:241–9. [60] Wang JD, Zhou D, Liu WQ. Response of liquid in cylindrical tank with rigid annual baffle considering damping effect. AMR 2011;255-260:3687–91. [61] Wang W, Peng Y, Zhou Y, Zhang Q. Liquid sloshing in partly-filled laterally-excited cylindrical tanks equipped with multi baffles. Appl Ocean Res 2016;59:543–63. [62] Wang Y, Benson DJ, Nagy AP. A multi-patch nonsingular isogeometric boundary element method using trimmed elements. Comput Mech 2015;56(1):173–91. [63] Ye W, Liu J, Lin G, Zhou Y, Yu L. High performance analysis of lateral sloshing response in vertical cylinders with dual circular or arc-shaped porous structures. Appl Ocean Res 2018;81:47–71. [64] Yue H, Chen J, Xu Q. Sloshing characteristics of annular tuned liquid damper (ATLD) for applications in composite bushings. Struct Control Health Monit 2018;25(8):e2184. [65] Zhang C, Ning D, Teng B. Secondary resonance of liquid sloshing in square-base tanks undergoing the circular orbit motion. Eur J Mech B/Fluids 2018;72:235–50. [66] Zheng X, Li X, Ren Y, Wang Y, Ma J. Effects of transverse baffle design on reducing liquid sloshing in partially filled tank vehicles. Math Probl Eng 2013;2013:1–13. [67] Zhou W, Liu B, Wang Q, Cheng Y, Ma G, Chang X, Chen X. NURBS-enhanced boundary element method based on independent geometry and field approximation for 2D potential problems. Eng Anal Bound Elem 2017;83:158–66.

165