Numerical simulation of structure response outfitted with a tuned liquid damper

Computers and Structures 87 (2009) 1154–1165

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Numerical simulation of structure response outfitted with a tuned liquid damper M. Marivani, M.S. Hamed * Thermal Processing Laboratory, Department of Mechanical Engineering, McMaster University, 1280 Main Street West, Hamilton, Ontario, Canada L8S 4L7

a r t i c l e

i n f o

Article history: Received 17 November 2008 Accepted 17 May 2009 Available online 16 June 2009 Keywords: Fluid–structure interaction Tuned liquid dampers Finite difference Volume of fluid Sloshing motion Vibration damping

a b s t r a c t An integrated fluid–structure numerical model has been developed to simulate the response of a single degree of freedom (SDOF) structure outfitted with a Tuned Liquid Damper (TLD). The structure is exposed to random external excitations. A non-linear, two-dimensional, flow model has been developed using the finite-difference method. Unlike most existing flow models, the present model does not include any linearization assumptions; it rather solves the entire nonlinear, moving boundary, flow problem under conditions leading to large interfacial deformations. The free surface has been reconstructed using the volume of fluid method and the donor–acceptor algorithm. The Duhamel integral method has been used to determine the response of the structure. The effectiveness and accuracy of the flow model has been validated using a set of benchmark problems and experimental data. The numerical results of this model have been compared with results of an equivalent TMD model. The present fluid–structure model can be used as a valuable tool for performance evaluation and design of more effective tuned liquid dampers. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction Tuned Liquid Dampers (TLDs) are economical and effective dynamic vibration absorbers. They have been increasingly used to control the dynamic response and protect structures from damage due to external excitations. A Tuned Liquid Damper (TLD) is properly designed, partially filled, water tank used as a vibration absorber to decrease the dynamic motion of a structure. It is a passive damping device. Unlike active or semi-active dampers, passive dampers do not require any external supply of power to operate [1]. The main function of a passive damping device is to absorb portion of the input energy associated with external dynamic excitation acting on the structure. Examples for external excitations are wind and earthquakes. By doing so, the passive damping device minimizes or eliminates the possibility of structural damages. A TLD consists of one or multiple rigid tanks attached to the structure. The vibration of the structure under an external dynamic load causes sloshing motion of the contained liquid. By tuning the natural frequency of this sloshing motion to the natural frequency of the structure, the liquid motion imparts inertia forces approximately anti-phase to the dynamic forces exciting the structure, thereby reducing structural motion. As such, the attachment of TLDs modifies the frequency response of the structure in a way similar to increasing its effective damping. The use of TLDs has increased due to their many advantages over other conventional damping devices. They can be used for * Corresponding author. Tel.: +1 905 572 9140x26113; fax: +1 905 572 7944. E-mail address: [email protected] (M.S. Hamed). 0045-7949/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2009.05.010

both small (wind) and large amplitude (earthquake) vibrations. They are easy to install to existing structures and require very low maintenance and operating cost. TLDs were used to stabilize marine vessels against rocking and rolling motions [2,3] in offshore platforms [4,5] and in tall structures [6–10]. Although the construction of a TLD is simple, the liquid sloshing motion inside it has a very nonlinear and complex nature. Surface slopes can approach infinity and the liquid may encounter the top cover in case of enclosed tanks. The amplitude of the sloshing motion depends on: the nature (e.g. harmonic or random), amplitude, and frequency of the applied external excitation, the geometry of the tank, the depth of the liquid layer, and the properties of the contained liquid [11,12]. Because of the difficulty in modeling the non-linear behavior of the liquid inside the TLD, a semi-empirical equivalent tuned-mass damper (TMD) model was introduced. [13–16]. In this model, the TLD is substituted with an equivalent ideal tuned-mass damper (TMD) with certain effective mass, frequency, and damping. These effective parameters are determined using experimental results. The main advantage of the TMD analogy is its simplicity; however, this approach does not give good results under some operating conditions [13–16]. The sloshing behavior of liquid in tanks subjected to dynamic excitations has been studied extensively through experimental investigations [17–24]. Numerical modeling has also been used to study sloshing behavior of liquids. The numerical investigations have been carried out mostly by using simplified linear theories that are applicable only to small interfacial deformations. The potential flow theory, which considers inviscid, irrotational flows, has been widely used by many researchers [25–31]. Shallow

M. Marivani, M.S. Hamed / Computers and Structures 87 (2009) 1154–1165

1155

Nomenclature A gx gy L h p ps t tf u,v x,y F Xe d dc FUX FUY VISX FVX FVY VISY mw MS CS

amplitude of external dynamic excitation (m) horizontal acceleration (m/s2) gravitational acceleration (m/s2) length of tank (m) height of the initial flat free surface (m) pressure (Pa) surface pressure (Pa) time (s) final simulation time (s) velocity components in the horizontal and vertical directions, respectively Cartesian coordinates volume of fluid per unit volume of cell (I, J) external harmonic excitation (mm) distance between free surface and center of the neighbor fluid cell distance between the center of the free surface cell and the center of the neighbor fluid convective term of the u velocity in the x-direction convective term of the u velocity in the y-direction diffusion term of the u velocity in the x-direction convective term of the v velocity in the x-direction convective term of the v velocity in the y-direction diffusion term of the v velocity in the y-direction water mass in TLD (Kg) mass of structure (Kg) damping of structure (Kg/s)

water wave theory has also been applied by others [20,32] where the wave height or the amplitude of the interfacial deformation was assumed to be small compared to the mean depth of the liquid layer, and the horizontal velocity was assumed uniform throughout the liquid layer. A number of studies [33–36] combined the shallow water theory with the boundary layer theory, where the liquid viscosity was considered only in the boundary layer region. An improvement over predictions based on the potential flow theory was provided by Zang et al. [37] using a linearized form of the Navier–Stokes equations by neglecting the convective acceleration terms. They indicated that liquid viscosity has an important effect on sloshing motion near rigid walls, especially under excitation frequency near resonance. The linear theory is applicable only when the amplitude of interfacial deformations are expected to be small. This is valid in cases of small external amplitudes or at frequency of excitation away from the natural frequency of the TLD. Consequently, when the excitation of the oscillation is near the natural frequency of the TLD, or the excitation amplitude is high, the linear model becomes inadequate [38]. Ramaswamy et al. [39] have solved the full non-linear Navier– Stokes equations using a Lagrangian approach and the finite-element method. They implemented the velocity-correction method for a two-dimensional liquid sloshing problem. Only small amplitude oscillations were considered in this study, and therefore the behavior of the liquid sloshing was assumed linear. A TLD is usually designed to operate at or near the resonant frequency of the structure in order to maximize the absorbed and dissipated energy, and thus maximize the benefit of the TLD. To accurately predict the sloshing motion inside a TLD, a numerical model should: (1) consider all physical mechanisms (inertia and viscosity), and (2) be able to solve such moving boundary problem while allowing large interfacial deformations. In such a problem, not only the transport of the momentum and mass are coupled,

KS XS Fe FTLD F0 fw P

stiffness of structure (N/m) displacement of structure (mm) external excitation force (N) sloshing force (N) normalized sloshing force natural frequency of TLD total momentum of fluid (Kg m/s)

Greek symbols distance coefficient factor density (kg/m3) kinematics viscosity (m2/s) dynamic viscosity (kg/m/s) b pressure correction coefficient e convergence criteria a upwind fraction factor sij shear stress (N/m2) r surface tension (N/m) / phase angle (rad/s) xr relaxation parameter xn natural frequency of structure (rad/s) xD damped frequency of structure (rad/s) x frequency of excitation (rad/s) n damping ratio of structure j local free surface curvature excitation frequency ratio bw

g q t l

but also the formation, evolution, and dynamics of the interface play major roles in the defining the system behavior. Actually, the main difficulty in solving such a problem is in the determination of the location of the liquid free surface as an integrated part of the solution of the system of equations governing the fluid flow problem. Algorithms for moving boundary problems were reviewed by Floryan and Rasmussen [40]. The moving free surface can be described either in terms of fixed grids, adaptive grids, or by applying analytical mapping techniques. Thé et al. [41] solved the full form of the primitive variable equations of motion of a viscous, incompressible, free surface flow using the finite-volume method. They tracked the interface using an adaptive grid. However, they did not consider large deformations of the free surface; only cases with moderate surface deformations were considered Yamamoto and Kawahara [42] solved the Navier–Stokes equations in the form of arbitrary Lagrangain–Eulerian (ALE) formulation. In the ALE method, the remeshing technique has been used to maintain the computational stability. The remeshing means in the computation at every computational step, the new finite element mesh is reformed. To compute a large amplitude sloshing problem, computation often tends to be unstable in this technique because the unevenness of nodal points is formed on the free surface. To reduce this instability, in the modified ALE method, a smoothing technique was introduced which it needed to a smoothing factor. Finding a proper value for this factor is the major issue of this technique because it is very difficult to choose a unique constant for the entire of computation. Siddique et al. [43] applied an analytical mapping technique. In this approach, an analytical mapping function is used to transfer the problem from the irregular physical domain onto a rectangular computational domain. The mapping function is unknown and had to be determined as part of the solution procedure. Although this

1156

M. Marivani, M.S. Hamed / Computers and Structures 87 (2009) 1154–1165

approach was successful in dealing with fairly large deformations, the method could not be applied to cases were surface discontinuities exist, such as in case of wave breaking and flow through submerged screens. The main objective of this work is to develop an integrated flow-structure numerical model that can be used to simulate the sloshing motion inside a TLD without using any linearization assumptions and to investigate the effect of the TLD on the response of the structure. 2. Mathematical formulation The TLD considered in this study is a rigid rectangular tank of width L and height H filled with water to a stationary depth h, as shown in Fig. 1. The coordinate system is attached to the tank with the origin located at the bottom left corner. The tank is subjected to a general external harmonic excitation in the horizontal (x) direction, defined by a time dependent displacement Xe(t). By neglecting the end effects, the motion of the liquid inside the tank can be assumed two-dimensional in the x–y directions. 2.1. Governing equations and boundary conditions The two dimensional, incompressible, free surface, fluid flow, problem has been modeled using an Eulerian frame of reference. The fixed points ~ x in the domain (~ x ¼ x^i þ y^j) are described using Cartesian coordinates. The velocity field ~ V is function of space and time, i.e.,

~ V ¼ uðx; y; tÞ^i þ v ðx; y; tÞ^j: The governing equations of the incompressible, Newtonian, laminar flow in the Cartesian coordinate system shown in Fig. 1 are the following continuity and momentum equations:

@u @ v þ ¼ 0; @x @y @u @u @u 1 @p 1 @ sxx 1 @ sxy þu þv ¼ þg þ þ ; @t @x @y q @x x q @x q @x @v @v @v 1 @p 1 @ syy 1 @ sxy þg þ þu þv ¼ þ : @t @x @y q @y y q @y q @y

stress boundary condition at the free surface can be written in tensor form as [52]:

^ i ¼ ðdik  n ^i n ^k Þ ðrj  ps Þn

@r ^k ;  sik n @xk

ð5Þ

^ i is the unit force normal to the where r is the fluid surface tension, n surface (into the fluid) and j is the local free surface curvature. For ^ two dimensional flows, projecting Eq. (5) along the unit normal n and unit tangent ^t results in an equivalent set of scalar boundary conditions. These are the normal stress boundary condition given by:

ps  rj ¼ 2lnk

@uk @n

ð6Þ

and the tangential stress boundary condition given by:



l ti

@ui @uk þ nk @n @s

 ¼

@r ; @s

ð7Þ

@ @ ^  r is the normal where, @s ¼ ^t  r is the surface derivative and @n ¼n derivative. The viscous effects at the free surface have been neglected and the liquid is treated as an ideal fluid. The surface tension, r, is assumed to be constant. Since the curvature radius of the free surface is expected to be large under conditions of interest in this study, the surface pressure (rj) effects has been ignored at the free surface. Thus, the normal stress boundary condition becomes:

ps ¼ 0:

ð8Þ

Eq. (8) is the inviscid free surface normal stress condition. Our validation showed that it is an acceptable approximation for the conditions considered in this study. 2.2. Discretization of the governing equations

The time evolution of liquid region is computed solving the following equation [44–51],

Discretization of the governing equations has been carried out in the physical space using the finite-difference method. A forward staggered grid arrangement is used, in which the scalar quantities are located at the geometric center of the cell, while the velocity components u and v lie at the midpoints of the cell faces. The locations of the various dependent variables using this grid configuration are shown in Fig. 2. The discretization of Eqs. (2) and (3) leads to the following expressions in the x- and y-directions:

@F @F @F þu þv ¼ 0; @t @x @y

unþ1 iþ12;j

ð1Þ ð2Þ ð3Þ

ð4Þ

where F is the local volume fraction of liquid phase. In the cells occupied with the liquid phase, F is unity, and in the cells occupied with the gas phase, F is zero. For the cells containing the interface bounding the liquid and gas phases, F lies between zero and unity. Eqs. (1)–(3) are subject to the no-slip velocity boundary condition, i.e., zero tangential velocity at the wall. The normal component of the velocity at the wall is also set to be zero due to non-penetrating wall boundary condition. The exact surface

"

# nþ1 nþ1 1 piþ1;j  pi;j ¼ þ dt  þ g x  FUX  FUY þ VISX ; ð9Þ q xiþ1  xi " # nþ1 nþ1 1 pi;jþ1  pi;j n v nþ1 ¼ v þ dt  þ g  FVX  FVY þ VISY ; ð10Þ 1 1 y i;jþ2 i;jþ2 q yjþ1  yj uniþ1;j 2

where, FUX, FUY, FVX and FVY represent the convection terms of the momentum equations in the u and v directions, respectively. Also VISX and VISY are associated with the diffusion terms of the

i,j+1/2

p ui-1/2,j

i,j

Fi,j

ui+1/2,j

vi,j-1/2 Fig. 1. Dimensions of the TLD.

Fig. 2. Location of dependent variables.

1157

M. Marivani, M.S. Hamed / Computers and Structures 87 (2009) 1154–1165

momentum equations in the x and y directions. Since the upwind scheme usually introduces significant numerical damping and the central difference scheme generates numerical instability, a combination of these two schemes usually yield a more accurate algorithm. Thus, the general finite-difference formula for the advection of u in the x and y directions become:

"      @u 1 @u @u n FUX ¼ u  ¼ uiþ1;j dxi þ dxiþ1 2 Dxa @x iþ1;j @x iþ1;j @x i;j 2 !#       @u @u þa:sign uiþ1;j dxiþ1 ;  dxi 2 @x i;j @x iþ1;j "      @u 1 n  @u þ @u ¼ v D y þ D y 1 iþ2;j @yiþ1;j @y iþ1;jþ1 @y iþ1;j1 Dy a 2 2 2 !#2 2       þ @u  @u þa:sign v iþ1;j Dy ;  Dy 2 @y iþ1;j1 @y iþ1;j1

dp ¼ 

ð11Þ

FUY ¼ v

2

2

2

ð13Þ

  Dya ¼ Dyþ þ Dy þ a:sign v iþ1;j ðDyþ  Dy Þ;

ð14Þ

1 2

where dxi, Dy+ and Dy have been in the Fig. 3. Here, a, the donor-cell fraction and the convective differencing, is a linear combination of the upwind and the centered differencing. When a = 0, this formulation reduces to a second order accurate, centered, difference approximation. When a = 1, the first order upwind form is recovered. The value of a in this study is 1. The diffusion term, VISX, is expressed as:

 @ 2 u @ VISX ¼ t  @x2 

iþ12;j

 @ 2 u þ 2 @y 

1 A:

ð16Þ

iþ12;j

Similar expressions are used for the advection terms FVX, FVY and the diffusion term VISY.

2

2

ð19Þ

2

1 ; dxi ðdxiþ1 þ dxi Þ 1 ¼ ; dxi ðdxi1 þ dxi Þ 1 ¼ ; dyj ðdyjþ1 þ dyj Þ

kiþ1 ¼

ð20Þ

ki1

ð21Þ

2

h i 1 dxiþ1 v i;jþ1 þ dxiþ1 v i;j1 þ dxi v iþ1;jþ1 þ dxi v iþ1;j1 ; 2 2 2 2 2ðdxi þ dxiþ1 Þ ð15Þ

0

q

 ; dt kiþ1 þ ki1 þ njþ1 þ nj1

where

2

  Dxa ¼ dxi þ dxiþ1 þ a:sign uiþ1;j ðdxiþ1  dxi Þ;

2

ð18Þ

2

ð12Þ

2

2

S ; @S=@p

where S is the value of the left side of Eq. (17) evaluated with the most updated values of u, v and p. Eq. (18) is derived by substituting the right side of the Eqs. (25)–(28) into the continuity equation, Eq. (17). 1:0 The quantity b ¼ @S=@p is obtained from:

b¼

where,

v iþ ;j ¼

computational cell occupied by the fluid. The solution is obtained by using the following iterative procedure: Eqs. (9) and (10) are evaluated using quantities at the previous time step, i.e., at t = n producing a provisional velocity field used as an estimate of the advanced time velocities. In cells that contain fluid, a pressure correction is calculated from:

njþ1 2

nj1 ¼ 2

ð22Þ

1 : dyj ðdyj1 þ dyj Þ

ð23Þ

The kth iteration for the pressure would be:

pki;j ¼ pk1 i;j þ dp

ð24Þ

and the velocities at all four faces are updated using:

dp

ukiþ1;j ¼ uk1 þ dt iþ1;j

qðxiþ1  xi Þ

 dt uki1;j ¼ uk1 i1;j

qðxi  xi1 Þ

2

2

2

2

v ki;jþ

1 2

¼ v k1 þ dt i;jþ1

v ki;j

1 2

¼ v k1  dt i;j1

2

2

dp

;

ð25Þ

;

ð26Þ

;

ð27Þ

:

ð28Þ

dp

qðyjþ1  yj Þ dp

qðyj  yj1 Þ

2.3. The numerical algorithm

This procedure is modified for the cells that contain free surface. For these cells, the pressure, pi,j, is obtained using linear interpolation between the surface pressure (ps), and a neighboring pressure (pn) inside the fluid in a direction perpendicular to the free surface, i.e.:

The velocities calculated using Eqs. (9) and (10) will not satisfy the continuity equation, which, after discretization, takes the form:

pi;j ¼ ð1  gÞpn þ gps ;

unþ1  unþ1 iþ1;j i1;j 2

2

dxi

þ

v nþ1  v nþ1 i;jþ i;j 1 2

1 2

dyj

ð29Þ

where the distance coefficient factor g is defined by:

¼ 0:0:

ð17Þ

In order to satisfy the continuity equation simultaneously throughout the mesh, pressures and velocities must be adjusted in each

g¼

dc ; d

ð30Þ

d and dc are defined in Fig. 4, where dc is the distance between the two adjacent grids and d is calculated from the orientation of the

j+1 y+

Free Surface

i,j yj i-1

i+1

i,j x-

x+ y-

d

dc i,j-1

j-1 xi Fig. 3. The computational cell.

Fig. 4. Geometric parameters of a surface cell.

1158

M. Marivani, M.S. Hamed / Computers and Structures 87 (2009) 1154–1165

free surface at the previous time step. The neighboring pressure (pn) is at (i, j  1) in Fig. 4. The new pi,j is obtained iteratively. The iteration is continued until all residuum (S) are sufficiently smaller than a small number, e which is typically of the order of 104. The last values of the velocity and pressure fields are used as the new guess in the next time step. In most cases, convergence of the iterations may be accelerated by using successive-over-relaxation by multiplying the dp by an over relaxation parameter xr. A value of xr = 1.7 has been used in this study. After the velocity and pressure fields have been determined at the new time step, the volume of fluid function, F, is determined at the current time step. The new F-field is calculated by solving Eq. (4). Combining this equation with the continuity equation gives:

@F @ðFuÞ @ðF v Þ þ þ ¼ 0:0: @t @x @y

ð31Þ

lows. If the free surface is convected mostly normal to itself, the acceptor cell F-value is used, otherwise the donor cell F-value is used. If the acceptor cell is empty, or if the cell upstream of the donor cell is empty; the acceptor cell is used to determine the flux, regardless of the orientation of the free surface. Determination of whether the interface is mostly horizontal or vertical is based on the value of the local slop of the fluid–void interface. One must note that the interface can be represented either as a single-valued function of (x) i.e., Y(x) or as a single-valued function of (y), i.e., X(y), depending on its orientation, which is determined based on the value of F in the surface cell and its eight neighbors. A good approximation of Y(xi) is:

Yðxi Þ ¼ F i;j1 dyj1 þ F i;j dyj þ F i;jþ1 dyjþ1 then

" #   ðY iþ1  Y i Þdxi1 ðY i  Y i1 Þdxiþ1 dY 1 2 2  ¼ þ ; dx i dxi1 þ dxiþ1 dxiþ1 dxi1 2

This conservative form of the F-equation allows one to write a difference approximation that conserves the fluid volume. The discretized form of Eq. (31) is written as follows: " #  1   1  nþ1 n nþ1 n nþ1 n nþ1 n nþ1 n F i;j ¼ F i;j  dt u 1 F 1  ui1;j F i1;j þ v 1 F 1  v i;j1 F i:j12 : 2 2 2 dxi iþ2;j iþ2;j dyj i;jþ2 i;jþ2

ð32Þ The convection algorithm must resolve the sharp interface at the free surface. Also, it must guarantee that there is not any flux across any interfacial cell more than the available amount in that cell. A donor–acceptor method has been used in order to satisfy these conditions. At each boundary of each cell, one of the two interface cells is designated as a donor cell and the other is designated as an acceptor cell. Cell quantities are given the subscript D for donor and A for acceptor, respectively. Each computational cell will have four separate assignments of D or A corresponding to each of the cell boundaries, as shown in Fig. 5. This method is based on the calculation of the amount of F fluxed and advected through the right-hand face of the donor cell during a time step dt. The volume flux crossing this cell face per unit cross sectional area is:

V x ¼ unþ1  dt; iþ1;j

ð33Þ

2

where u is the normal velocity at the cell face. The sign of u determines the donor and the acceptor cells. For example, if u is positive, the upstream or left cell is the donor and the downstream or right cell is the acceptor. The amount of F advected across the cell face in one time step is dF times the face cross section area, where:

dF ¼ MINðF AD jV x j þ CF; F D dxD Þ;

ð34Þ

and; CF ¼ MAXfð1:0  F AD ÞjV x j  ð1:0  F D ÞdxD ; 0:0g:

ð35Þ

The double subscript, AD, refers to either the acceptor (A) or donor (D), depending on the orientation of the interface relative to the direction of flow. The rules of choosing AD = A or AD = D is as fol-

ð36Þ

where; dxiþ1 ¼ 2

2

2

dxi þ dxiþ1 : 2

A similar calculation of

ð37Þ

2

ð38Þ

  dX dy

j

can be made from,

Xðyj Þ ¼ F i1;j dxi1 þ F i;j dxi þ F iþ1;j dxiþ1

ð39Þ

" #   ðX jþ1  X j Þdyj1 ðX j  X j1 Þdyjþ1 dX 2 2 ¼ þ then dy j dyjþ1 dyj1 2

ð40Þ

dyj þ dyjþ1 : 2

ð41Þ

2

where; dyjþ1 ¼ 2

2

1  ; dyjþ1 þ dyi1 2

The general rule to determine the orientation of the free surface is as follows:   dX     >   and dY  < 0, then free surface is vertical with the If dY dx dy dx liquid  left side.  the dY  on dY    dX If  dx  >  dy  and  dx  > 0, then free surface is vertical with the liquid left side.   on the     dY    > dx and dX < 0, then free surface is horizontal with If dX dy  dy  the the bottom side.   on   liquid    dY    > dx and dX > 0, then free surface is horizontal with If dX dy  dy  the liquid on the top side. Once the slope of the interface and the side occupied by the liquid are determined, a line can be constructed in the cell with the correct amount of F placed on the liquid side. This line is used as an approximation of the actual interface and provides the information necessary to calculate d and dC, and hence, the distance coefficient factor, g, for the application of the free surface pressure condition, Eq. (29). 2.4. Conditions for numerical stability Once a mesh size is chosen, the value of the time step must be determined based on some stability conditions. In this work, two conditions have been applied. First, the fluid cannot move through more than one cell in one time step. Therefore, the time step must satisfy the following condition:



Dt  min

Fig. 5. Donor–acceptor cell configuration.

 Dx Dy : ; jui;j j jv i;j j

ð42Þ

The minimum of Dt is calculated for every cell in the entire mesh. Second, momentum must not diffuse more than one cell in one time step. Linear stability analysis shows that this limitation implies that:

1159

M. Marivani, M.S. Hamed / Computers and Structures 87 (2009) 1154–1165

Dt 

1 Dx2 Dy2 : 2t Dx2 þ Dy2

ð43Þ

2.5. The fluid–structure interaction model The TLD is coupled with a SDOF structure as shown in Fig. 6. In this fluid–structure system, the generated force by the water sloshing motion inside the TLD, FTLD, acts as a resisting or damping force to the external force, Fe. The equation of motion of this fluid–structure system can be written in the following form [23]:

€ S þ C S X_ S þ K S X S ¼ F e þ F TLD ; MS X

mij ¼ q  Dxi Dyj :

ð45Þ

The total momentum of the fluid inside each control volume is: nj ni X X i¼1

mij uij :

ð46Þ

j¼1

The damping force, FTLD, can then be determined from the rate of change of the fluid momentum using the following equation:

F TLD ¼

1 ðPðtÞ  Pðt þ DtÞÞ: Dt

ð47Þ

The Duhamel integral method [53] has been used to solve the equation of motion of the structure. This integral method can be used to determine the response of the SDOF system under any type of external excitations (harmonic and random). The total displacement of a damped SDOF system subjected to an arbitrary external force can be expressed as:

  u0 þ X 0 nxn X S ðtÞ ¼ enxn t X 0 cos xD t þ sin xD t :

xD

X S ðtÞ ¼

1 M S xD

t

enxn t fAD ðtÞ sin xD t  BD ðtÞ cos xD t g; M S xD

ð50Þ

where, AD and BD can be evaluated from:

AD ðti Þ ¼ AD ðt i1 Þ þ BD ðt i Þ ¼ BD ðt i1 Þ þ

Z Z

ti

F t ðsÞenxn s cosðxD sÞds;

ð51Þ

F t ðsÞenxn s sinðxD sÞds:

ð52Þ

ti1 ti

t i1

Considering a linear piecewise loading function, the forcing function Ft(s) can be approximated by:

F t ðsÞ ¼ F t ðti1 Þ þ

DF i ðs  t i1 Þ; Dt i

ti1  s  t i ;

ð53Þ

where,

DF i ¼ F t ðti Þ  F t ðti1 Þ;

ð54Þ

F t ðt i Þ ¼ F e ðt i Þ þ F TLD ðti Þ; Dti ¼ t i  t i1 :

ð55Þ ð56Þ

In Eqs. (51), (52), AD (ti) and BD(ti) can be evaluated from:

  DF i DF i I1 þ AD ðti Þ ¼ AD ðt i1 Þ þ F t ðt i1 Þ  ti1 I4 ; Dt i Dt i   DF i DF i BD ðt i Þ ¼ BD ðt i1 Þ þ F t ðti1 Þ  ti1 I3 ; I2 þ Dt i Dt i

ð57Þ ð58Þ

where the integrals I1, I2, I3 and I4, are evaluated as follows:

I1 ¼

Z

ti

enxn s cos xD sds

t i1

ð48Þ

Setting X0 = 0 and u0 = Ft(s)ds/MS and substituting t for t  s in Eq. (48) and integration of this equation over the entire loading interval results in:

Z

X S ðtÞ ¼

ð44Þ

where, MS, CS and KS are the structure mass, generalized damping, and stiffness, respectively. XS is the structure displacement. The damping force due to the sloshing motion has been calculated by applying the momentum theory. The mass of each fluid control volume is calculated suing the following equation:

P¼

using Duhamel’s integral. The displacement is calculated by numerical integration of Eq. (49). Using the trigonometric identity sin x(t  s) = sin xt cos xs  cos xt sin xs the displacement can be calculated from:

¼

ti   ðn x cos x s þ x sin x s Þ  ; n D D D 2 2  ðnxn Þ þ xD t enxn s

ð59Þ

i1

I2 ¼

Z

ti

enxn s sin xD sds

t i1

F t ðsÞenxn ðtsÞ sin xD ðt  sÞds;

ð49Þ

0

where, Ft is the sum of the external excitation force and the TLD sloshing force. XS(t) in Eq. (49) is the response of a damped system

ti   ¼ ðnxn sin xD s þ xD cos xD sÞ ; 2  ðnxn Þ þ x2D t enxn s

ð60Þ

i1

I3 ¼

Z

ti

s:enxn s sin xD s ds

t i1

ti   ¼s I2 þ I1  ; 2 2 2 2 ðnxn Þ þ xD ðnxn Þ þ xD t n xn

w

FTLD

Z

ti

s:enxn s cos xD sds

t i1

¼

KS

CS

s

n xn ðnxn Þ2 þ x2D

! I1 

ti   I : 2 2 2 ðnxn Þ þ xD t

xD

ð62Þ

i1

The substitution of Eqs. (54) and (55) into Eq. (50) gives the displacement at time ti as:

X S ðt i Þ ¼

Fig. 6. A SDOF structure outfitted with a TLD.

ð61Þ

i1

I4 ¼

MS

Fe

xD

enxn ti fAD ðt i Þ sin xD t i  BD ðti Þ cos xD ti g: M S xD

ð63Þ

The flowchart of the entire fluid–structure numerical algorithm is shown in Fig. 7.

1160

M. Marivani, M.S. Hamed / Computers and Structures 87 (2009) 1154–1165

Fig. 7. Flow chart of the present fluid–structure algorithm.

3. Model validation The algorithm has been validated using the following benchmark and test problems. 3.1. The broken dam problem The first benchmark problem is the broken dam problem. A rectangular column of water, originally under hydrostatic equilibrium, is confined between two vertical walls. The appearance of both vertical and horizontal free surfaces in this problem provides a rigorous test of the capability of the algorithm to handle the free surfaces that are not single valued with respect to x or y. The water column is 1.0 unit wide and 2.0 units high, as shown in Fig. 8. Gravity is acting downwards with unit magnitude. At the beginning of the calculation, the right wall of the dam is removed suddenly and water is allowed to flow out along a dry horizontal floor. Experimental data of this problem have been reported in [54] providing the position of the leading edge of the water surface as function of time as water flows to the right. Figs. 9 and 10 show the predicted horizontal and vertical positions of the free surface, respectively, as function of time compared with the experimental data. Fig. 11 shows the numerically predicted shape of the free surface using the present algorithm. These results show good agreement with the experimental data.

Fig. 9. Horizontal position of free surface as function of time versus experimental data in [52].

3.2. Sloshing motion of a wave under the effect of gravity The second test case is the problem reported in Raad et al. [55], related to the sloshing motion of a low-amplitude surface wave

Fig. 10. Vertical position of the free surface as function of time versus experimental data reported in [52].

Fig. 8. The layout of the broken dam problem.

under the influence of gravity as shown in Fig. 12. Initially the quiescent fluid has an average depth of 0.05 m, and its surface is initially defined by: y(x) = 0.05 + 0.005 cos (px/L). The computational

M. Marivani, M.S. Hamed / Computers and Structures 87 (2009) 1154–1165

1161

Fig. 11. Free surface profile at various times.

where k is the wave number and h is the average fluid depth. Fig. 13 shows a reasonable agreement between the theoretical predictions of the height of the liquid free surface at the left boundary reported in [53] and the numerical results obtained using the present algorithm. Fig. 14 shows the evolution of the free surface during one period (T) of the sloshing motion predicted by the present algorithm. The system is initially at rest. After t = T/4, the potential energy of the system transferred into kinetic energy and the fluid velocity reaches its maximum. After t = T/2, all the kinetic energy is transferred back into potential energy and the fluid velocity goes back to zero. 3.3. Sloshing of a shallow water layer in a TLD

Fig. 12. Initial profile of the gravity induced sloshing motion of the low-amplitude surface wave reported in [53].

Fig. 13. Position of the interface at the left boundary as a function of time.

domain is a 0.1 m wide by 0.065 m high rectangle. The fluid begins to slosh solely under the effect of gravity. The theoretical period of the first mode this sloshing motion is:

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi T ¼ 2p gk tanhðkhÞ ¼ 0:3739;

ð64Þ

The third problem is sloshing of a shallow water layer in a TLD subjected to a sudden sinusoidal external excitation in the horizontal direction. The TLD is partially filled with water and is initially at rest. The external excitation imposes a horizontal displacement given by Xe = Asin(xt + /), where A, x, and / are the amplitude, the frequency, and the phase angle of the forced oscillation, respectively. The results presented here are for the amplitude, A, the period of oscillation, T, and the phase angle, /, of 0.259 cm, 1.681 s, and 4.0°, respectively. The initial depth of the liquid layer, h, is 0.119 m, and the tank width, L, is 0.966 m. These parameters have been selected according to the experiments reported in [23]. The numerical simulations have been carried out using a 200  100 uniform mesh. Grid dependence has been checked using 120  60 and 160  80 meshes. Table 1 shows the maximum difference in the free surface deflection using these three different mesh sizes. Fig. 15 shows the variation of the height of the free surface as function of time at x = 0.05  L. Numerical results of the present algorithm are in good agreement with the experimental data reported in [23]. It is worth noting that the interfacial deformations depicted in Fig. 15 are about 20% of the original height of the liquid layer, which confirms the capability of the present algorithm in capturing large interfacial deformations. There is a slight phase shift between the numerical and the experimental results shown in Fig. 15 for time less than 30 s. This phase shift is due to a slight phase shift between the actual excitation used in the experiments and the one used in the numerical computations, see Fig. 16. In agreement with the results presented in Figs. 15 and 16 shows that this phase shift disappeared after about 30 s.

1162

M. Marivani, M.S. Hamed / Computers and Structures 87 (2009) 1154–1165

Fig. 14. Development of the free surface during the first period of the sloshing motion.

Table 1 Grid dependence test. Mesh size

Maximum deviation

200100 16080 12060

– (Selected mesh) 1.8% 11%

3.4. The TLD–TMD analogy The performance of a TLD is often treated by analogy with a TMD. The TMD parameters are usually determined by matching

the energy dissipation or the shear force produced by the TLD with that produced by an equivalent single degree of freedom TMD. The energy dissipation or the shear force produced by the TLD is usually determined experimentally. Numerical results of the present algorithm have been compared with results of an equivalent amplitude dependent TMD model. Parameters of the equivalent TMD model were obtained experimentally [23]. The TMD parameters were determined by matching the energy dissipated by the TLD and by an equivalent single degree of freedom TMD. The TLD was exposed to a sinusoidal excitation with amplitudes equal to 2.5 mm and 12.5 mm. Figs. 17 and 18 show the normalized sloshing force as function of the excitation frequency for the two amplitudes. The sloshing force has been normalized by the inertia force Numerical Results Experimental Results [23]

Deflection of free surface (mm)

20 15 10 5 0 -5 -10 -15

0

10

20

30

40

50

Time(s) Fig. 15. Surface deformations at a location 5% of tank length from the left wall (i.e., x = 0.0483 m) compared with experimental data reported in [23].

1163

M. Marivani, M.S. Hamed / Computers and Structures 87 (2009) 1154–1165

Excitation amplitude (mm)

10

Real excitation [21] Sinousoidal excitation

5

0

-5

-10

10

20

30

40

50

Time(s) Fig. 16. Comparison of the real forced excitation reported in [21] and the excitation used for code validation.

1 fw ¼ 2p

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ffi pg ph : tanh L L

ð65Þ

Substituting L and h into the above equation, leads to value of fw  0.545 Hz. Figs. 17 and 18 show a comparison of the normalized sloshing force calculated by the present numerical model, the measured sloshing force reported in [23], and the estimated base shear force obtained using the equivalent TMD model. In both small and large amplitude cases, the agreement between the present model and the experiment is excellent. However, results of the equivalent TMD model, especially in the case of large amplitude of excitation, did not produce satisfactory results. Figs. 17 and 18 clearly indicate that the present numerical model can be used to carry out an extensive investigation of the performance of any TLD using the TMD analogy. 4. Response of a SDOF structure outfitted with a TLD Fig. 17. Variation of the normalized sloshing force F, as function of excitation frequency ratio bw. Comparison between results of experiments [23] equivalent TMD analogy, and present numerical model ay amplitude of excitation of 2.5 mm.

The effect of TLD on the response of a SDOF structure has been investigated using the present algorithm. The TLD–SDOF system considered in this study is shown in Fig. 6. The characteristics of this system are listed in Tables 2 and 3. These values are the same as those used in the experimental study carried out by Tait [23]. In that study, the response of the SDOF structure to a random excitation has been tested. The actual random excitation used in [23] is shown in Fig. 19. Tait measured the structure response with and without a TLD. Numerical simulations of the two cases have been carried out using the present model. Fig. 20 shows the predicted response of the structure with and without the TLD. Our numerical results indicate that the TLD reduced the structure response by more than 70%, which is in agreement with the trend reported in [23]. 5. Summary and conclusions A numerical algorithm has been developed to investigate the response of a SDOF structure outfitted with a TLD. In order to predict Table 2 Properties of the SDOF structure.

Fig. 18. Variation of the normalized slohing force F, as function of excitation frequency ratio bw. Comparison between results of experiments [23] equivalent TMD analogy, and present numerical model ay amplitude of excitation of 12.5 mm.

of the liquid treated as a solid mass, mw(2pf)2A. The excitation frequency ratio, bw, is the ratio of the frequency of excitation to the natural frequency of the TLD calculated from the linear wave theory according to the following equation:

MS (Kg)

KS (N/m)

f (%)

4480

55 100

0.1

Table 3 Properties of the TLD. L (m)

H (m)

H (m)

0.966

0.3

0.119

1164

M. Marivani, M.S. Hamed / Computers and Structures 87 (2009) 1154–1165

Fig. 19. External random excitation used in [23].

Fig. 20. Structure responses with and without the TLD.

the effect of the sloshing motion inside the TLD on the structure response, a non-linear, two-dimensional, flow model has been developed. The flow model solves the complete set of Navier–Stokes equations using primitive variables and the finite-difference method. Tacking of the moving free surface has been carried out using the volume of fluid method and the donor–acceptor algorithm. The developed flow model has been validated using three sloshing flow problems. Two of which have experimental data, and the third has a theoretical solution. The results show that the present non-linear flow model is capable of predicting non-linear sloshing motions under conditions leading to large interfacial deformations. The fluid–structure model has been tested. Results show that the model is capable of capturing the right expected damping effect of the TLD on the structure response. The present numerical model can be utilized to carry out extensive investigations of the performance of any TLD using the TMD analogy. Acknowledgments The authors would like to acknowledge the financial support received from the National Sciences and Engineering Research Counsel of Canada (NSERC). We would also like to acknowledge the Shared Hierarchical Academic Research Computing Network (SHARCNET). This work was made possible utilizing their computational facilities. References [1] Soong TT, Dargush GF. Passive energy dissipation systems in structural engineering. John Wiley; 1997. [2] Watanabe S. Methods of vibration reduction. Proc Jpn Naval Archit Soc Symp 1969:156–89.

[3] Matsuuara Y, Matsumoto K, Mizuuki M, Arima K, Jouuchi H, Hayashi S. On a mean to reduce excited vibration with the sloshing in a tank. J Soc Naval Archit Jpn 1986;160:424–32. [4] Vandiver JK, Mitone S. Effect of liquid storage tanks on the dynamic response of offshore platforms. Appl Ocean Res 1978;1:67–74. [5] Lee SC, Reddy DV. Frequency tuning of offshore platforms by liquid sloshing. Appl Ocean Res 1982;4:226–31. [6] Kareem A, Sun WJ. Stochastic response of structures with fluid containing appendages. J Sound Vib 1987;21:389–408. [7] Tamura Y, Fujii K, Sato T, Wakahara T, Kosugi M. Wind-induced vibration of tall towers and practical applications of tuned sloshing damper. In: Proceedings of the workshop on serviceability of buildings, Ottawa, Ontario, 1988; p. 228–41. [8] Noji T, Yoshida H, Tatsumi E, Kosaka H, Hagiuda H. Study on vibration control damper utilizing sloshing of water. J Wind Eng 1988;37:557–66. [9] Fujii K, Tamura Y, Sato T, Wakahara T. Wind-induced vibration of tower and practical applications of tuned liquid damper. J Wind Eng Ind Aerodynam 1990;33:263–72. [10] Wakahara T. Wind-induced response of TLD-structure coupled system considering non-linearity of liquid motion. Shimizu Tech Res Bull 1993;12:41–52. [11] Koh CG, Mahatma S, Wang CM. Theoretical and experimental studies on rectangular liquid dampers under arbitrary excitations. Earthquake Eng Struct Dynam 1994;23:7–31. [12] Celebi MS, Akyildiz H. Nonlinear modeling of liquid sloshing in a moving rectangular tank. Ocean Eng 2002;29:1527–53. [13] Casciati F, De Stefano A, Matta E. Simulation a conical tuned liquid damper. Simulat Model Practice Theory 2003;11:353–70. [14] Sun LM, Fujino Y, Chaiseri P, Pacheco PM. The properties of tuned liquid dampers using a TMD analogy. Earthquake Eng Struct Dynam 1995;24:625–36. [15] YU JK. Nonlinear characteristics of tuned liquid damper. PhD thesis, University of Washington; 1997. [16] Yu JK, Wakahara T, Reed DA. A non-linear numerical model of the tuned liquid damper. Earthquake Eng Struct Dynam 1999;28:671–86. [17] Chester W. Resonant oscillations of water waves. II. Experiment. Proc Royal Soc London 1968;306(A):23–39. [18] Miles JW. Nonlinear surface waves in closed basins. J Fluid Mech 1976;75:419–48. [19] Clough RW, Niwa A, Clough DP. Experimental seismic study of cylindrical tanks. J Struct Div ASCE 1978;105:2565–90. [20] Reed D, Jinkyu Y, Harry Y. Investigation of tuned liquid dampers under large amplitude excitation. J Eng Mech ASCE 1988;124(4):405–13.

M. Marivani, M.S. Hamed / Computers and Structures 87 (2009) 1154–1165 [21] El Damatty AA. Studies on the application of tuned liquid dampers (TLD) to upgrade the seismic resistance of structures. Report submitted to Institute of Catastrophic Loss Reduction; December 2001. [22] Ju YK, Yoon SW, Kim SD. Experimental evaluation of a tuned liquid damper system. Struct Build 2004;157(4):251–62. [23] Tait M. The performance of 1-D and 2-D tuned liquid dampers. PhD thesis. University of Western Ontario, London, Canada; 2004. [24] Hamelin J. The effect of screen geometry on the performance of a tuned liquid damper. MSc. thesis, McMaster University; 2007. [25] Nakayama T, Washizu K. The boundary element method applied to the analysis of two-dimensional nonlinear sloshing problems. Int J Numer Methods Eng 1981;17:1631–46. [26] Nakayama T. Boundary element analysis of nonlinear water wave problems. Int J Numer Methods Eng 1983;19:953. [27] Ohyama T, Fujii K. A boundary element method analysis for two-dimensional nonlinear sloshing problem. JSCE J Struct Eng 1989;35(A):575–84. [28] Tosaka N, Sugino R. Boundary element analysis of potential flow with nonlinear free surface. In: Proceedings of the fourth China–Japan symposium on boundary element methods; 1991. p. 137–44. [29] Chen W, Harou MA, Liu F. Large amplitude liquid sloshing in seismically excited tanks. Earthquake Eng Struct Dynam 1996;25:653–69. [30] Warnitchai P, Pinkaew T. Modeling of liquid sloshing in rectangular tanks with flow-dampening devices. Eng Struct 1998;20(7):593–600. [31] Dutta S, Laha MK. Analysis of the small amplitude sloshing of a liquid in a rigid container of arbitrary shape using a low-order boundary element method. Int J Numer Methods Eng 2000;47:1633–48. [32] Shimizu T, Hayama S. Nonlinear response of sloshing based on the shallow water wave theory. JSME Int J 1987;30(263):806–13. [33] Sun LM, Fujino Y, Pacheco BM, Isabe M. Nonlinear waves and dynamic pressures in rectangular TLD-simulation and experimental verification. J Struct Eng Earthquake Eng 1989;6:251–62. [34] Fujino Y, Sun LM, Pacheco BM, Chaiseri P. Tuned liquid damper for suppressing horizontal motion of structures. ASCE J Eng Mech 1992;118:2017–30. [35] Sun LM, Fujino Y. A semi-analytical model for tuned liquid damper (TLD) with wave breaking. J Fluids Struct 1994;8:471–88. [36] Sun LM, Fujino Y, Koga K. A model of tuned liquid damper for suppressing pitching motions of structure. Earthquake Eng Struct Dynam 1995;24:625–36. [37] Zang Y, Xue S, Kurita S. A boundary element method and spectral analysis model for small-amplitude viscous fluid sloshing in couple with structural vibrations. Int J Numer Methods Fluids 2000;32:79–96.

1165

[38] Lepelletier TG, Raichlen F. Nonlinear oscillations in rectangular tanks. J Eng Mech 1988;114:1–23. [39] Ramaswamy B, Kawahara M, Nakayama T. Lagrangian finite element method for the analysis of two-dimensional sloshing problems. Int J Numer Methods Fluids 1986;6:659–70. [40] Floryan JM, Rasmussen H. Numerical methods for viscous flows with moving boundaries. Appl Mech Rev 1989;42:323. [41] Thé JL, Raithby GD, Stubley GD. Surface adaptive final volume method for solving free surface flows. Numer Heat Transfer 1994;26(B):367–80. [42] Yamamoto K, Kawahara M. Structural oscillation control using tuned liquid damper. Comput Struct 1999;71:435–46. [43] Siddique MR, Hamed MS, El Damatty AA. A nonlinear numerical model for sloshing motion in tuned liquid dampers. Int J Numer Methods Heat Fluid Flow 2004;15(3):306–24. [44] Hirt CW, Nichols BD. Volume of fluid (VOF) method for the dynamics of free boundaries. J Comp Phys 1981;399:201. [45] Torrey MD, Cloutman LD, Mjolsness RC, Hirt CW. NASA-VOF2D: a computer program for incompressible flows with free surfaces. LA-10612-MS, Los Alamos National Lab; 1985. [46] Lemos CM. Wave breaking: a numerical study. Lecture notes in engineering. Springer-Verlag; 1992. [47] Kothe DB. Comments on modeling interfacial flows with volume of fluid methods. Los Alamos National Laboratory, 65P05; 1995. [48] Kothe DB, Rider WJ. Volume tracking of interface having surface tension in two and three dimensions. AIAA J 1996. 96-0859. [49] Rider WJ. Reconstructing volume tracking. J Comp Phys 1998;141:112–52. [50] Pengzhi Lin. A fixed grid model for simulation of moving body in free surface flows. Comput Fluids 2007;36:549–61. [51] Chen S, Johnson DB, Raad PE. Velocity boundary conditions for the simulation of free surface fluid flow. J Comput Phys 1995;116:262–76. [52] Philip LF, Lin Pengzhi. A numerical model for breaking waves: the volume of fluid method. Technical report, University of Delaware, Center for applied coastal research; 1997. [53] Paz M. Structural dynamics theory and computation. 3rd ed. USA: Chapman and Hall; 1980. [54] Ubbink O. Numerical prediction of two fluid systems with sharp interface. PhD thesis, Imperial College, London. UK; 1997. [55] Raad PE, Johnson DB, Chen S. Simulation of impacts of fluid free surfaces with solid boundaries. Int J Numer Methods Fluids 1994;19:153–76.