Sloshing waves in a heated viscoelastic fluid layer in an excited rectangular tank

Sloshing waves in a heated viscoelastic fluid layer in an excited rectangular tank

Physics Letters A 378 (2014) 3289–3300 Contents lists available at ScienceDirect Physics Letters A Sloshing waves in a ...

995KB Sizes 0 Downloads 57 Views

Physics Letters A 378 (2014) 3289–3300

Contents lists available at ScienceDirect

Physics Letters A

Sloshing waves in a heated viscoelastic fluid layer in an excited rectangular tank Magdy A. Sirwah Department of Mathematics, Faculty of Science, Tanta University, Tanta, Egypt

a r t i c l e

i n f o

Article history: Received 8 April 2014 Received in revised form 8 September 2014 Accepted 9 September 2014 Available online 26 September 2014 Communicated by C.R. Doering Keywords: Sloshing waves Moving tank Viscoelastic fluid Heat transfer Laplace transform

a b s t r a c t In this paper, we have investigated the motion of a heated viscoelastic fluid layer in a rectangular tank that is subjected to a horizontal periodic oscillation. The mathematical model of the current problem is communicated with the linearized Navier–Stokes equation of the viscoelastic fluid and heat equation together with the boundary conditions that are solved by means of Laplace transform. Time domain solutions are consequently computed by using Durbin’s numerical inverse Laplace transform scheme. Various numerical results are provided and thereby illustrated graphically to show the effects of the physical parameters on the free-surface elevation time histories and heat distribution. The numerical applications revealed that increasing the Reynolds number as well as the relaxation time parameter leads to a wider range of variation of the free-surface elevation, especially for the short time history. © 2014 Elsevier B.V. All rights reserved.

1. Introduction Usually, fluid in a container is influenced by the motion of the container. Such fluid assumes the form of waves, a phenomenon that is referred as sloshing. Liquid sloshing is such a complicated phenomenon. Depending on the type of disturbance and container shape, the free liquid surface can experience different types of motion including simple planar, non-planer, rotational, irregular beating, symmetric, asymmetric, quasi-periodic and chaos [1]. The sloshing of fluid in a tank is a well-studied problem that has applications in a number of practical situations. As indicated by Virella et al. [2], the problem is relevant to the safety of transporting fluids in tankers; Hermann and Timokha [3] also stress its relevance to the automotive, aerospace and shipbuilding industries. Fluid sloshing in road tankers may result in overturning of the vehicle, and resonant movement of fluid within ship cargo tanks is also of concern. These practical considerations are discussed further in the review article by Ibrahim and Pillipchuk [1]. They are reinforced by Frandsen [4], who also discusses the use of fluidfilled tanks to act as dampers on the motion of city buildings in high winds. If sloshing is prolonged or accentuated, the wall of the container is damaged by the body force of the fluid, which in turn, results in unstable behavior of the whole mechanical system. Therefore, we need to clearly understand the sloshing-behavior characteristics of fluids. Published research into the sloshing phenomenon of general fluids theoretically analyzes the characteristics of hydrodynamic behavior through numerical analysis. It also ad 0375-9601/© 2014 Elsevier B.V. All rights reserved.

dresses the design of apparatus and fluid–structure interaction, towards the reduction of sloshing. Liu and Lin, and Ikeda and Ibrahim analyzed sloshing using numerical analysis [5,6]. Abbas and Mansour investigated the damping of sloshing that arises from the use of baffles in fluid containers [7]. Also, Kim, Lee and Cho investigated optimization design technique for sloshing and free surface tracking for nonlinear liquid sloshing [8,9]. However, sloshing of general fluid is hard to control by itself; for reducing sloshing, additional structural elements are required in the fluid container [10–12]. Surface waves are primarily dominated by gravity and in the absence of solid boundaries viscosity usually has very little effect on the flow over a short period in time or over a short distance in space [13]. In other words, the effects of viscosity may become important only after many wave periods or after many wavelengths. It is quite common, therefore, that free surface flows of inviscid liquids are analyzed by the velocity potential theory. The problem becomes somewhat different when a wave encounters a body in its path, because of the sheared flow created by the body surface. But even in that case it is usually a common practice to deal with the free-surface effects and viscous effects separately. A typical example is linear wave interaction with an offshore structure. The interaction between the wave and the body can be analyzed by either wave-diffraction theory without viscosity or by viscous-flow theory without the free-surface effect, depending on the ratios of the characteristic dimension of the body to the wavelength and the wave amplitude [14]. Cases, however, do arise where the combined


M.A. Sirwah / Physics Letters A 378 (2014) 3289–3300

effect of the free surface and viscosity is important, some of which have been highlighted by Yeung and Yu [15]. A significant case is the flow near the liquid line of a floating body, i.e., at the intersection of the body surface and the free surface. Another example is the free surface flow of highly viscous fluid. In studying the behavior of the free surface interacting with the viscosity within a bounded region, it is assumed initially, common in potential flow analysis, that disturbance of the fluid is small and the flow is governed by the linearized Navier–Stokes (NS) equations. The justification and limitations of such an approximation have been discussed by Mei [16]. The analysis there is, however, based on the framework of boundary-layer theory. Such a case is easier than the case where the viscous effect is taken into account in the entire fluid domain. One difficulty in combined free-surface and viscous analysis lies at the intersection of the body surface and the free surface. Take a fixed body as an example: the no-slip condition in the NS equations suggests that the fluid particle there should remain stationary, but it can be observed experimentally that fluid moves up and down along the body surface. This difficulty was in fact mentioned by Lamb [17], and resolving it requires extensive experimental study such as that undertaken in [18]. Our current understanding of the flow structure near the intersection is still limited and many methods used to deal with the intersection are not entirely satisfactory. A commonly used scheme for water wave/structure interaction is based on two steps: (1) the NS equations are solved with the no-slip condition imposed on the body surface and (2) the motion of the intersection point is tracked through interpolation from the points on the free surface and near the intersection. This procedure has some clear defects. To ensure the result from the interpolation is accurate enough, the points used must be as close as possible to the intersection. However, if these points are sufficiently close to the body surface (for example when an extremely fine mesh is used in this region), the result will be the same as that based on the no-slip condition. Because of this difficulty, the analysis in this work remains modest. The no-slip condition on the body is replaced by a no-shear-force condition. The intention here is to show how the free surface will interact with the viscous flow, and results based on this model would be useful for this purpose. Indeed it was argued in [19] that the condition on the side walls may have little effect on the wave when analyzing Faraday’s instability. Furthermore, the equivalent of a zero shear force condition on the side walls was also used by Loh and Rasmussen [20], who solved the full Navier–Stokes equations for this problem based on the finite-difference method. When the temperature at the interface between two immiscible fluids is not uniform, a flow may be induced due to the temperature dependence of surface tension. This flow is usually called thermocapillary convection or Marangoni convection and can arise in a liquid–gas or a liquid–liquid interface subjected to a temperature gradient. As a consequence, the dependence of the surface tension upon the local temperature will create a shear stress on the liquid surface, which by viscous traction results in a Marangoni convection in the bulk of the liquid. In recent years, thermocapillary convection in systems with free boundaries have attracted much attention, mainly owing to their relevance in many processes of technological interest. Experimental study of thermocapillary convection is complicated by the presence of strong gravity convection under earthbound conditions and by residual mass–force accelerations and vibrations under microgravity conditions on board rocket probes and spacecraft. The simplest convective flow appears when a free surface of a single extended liquid layer is subject to a horizontal temperature gradient. As soon as a lateral heating is applied a basic flow settles down in the liquid layer and the resulting temperature profile is highly nonlinear. This flow destabilizes through an oscillatory instability that induces wave motions in simple fluids called hydrothermal waves (Smith and Davis [21]). Smith and

Davis [21] showed that the instability mechanism even in this simple case is quite complex, resulting from the interplay between the basic flow and thermal or velocity disturbances. At low Prandtl numbers (Pr), hydrothermal waves propagate in a direction perpendicular to the horizontal temperature gradient, while at high Pr, they advance parallel to the temperature gradient. At intermediate Pr, the waves form an angle with the streamwise direction (Smith [22]). Later on, some authors (Parmentier et al. [23], Mercier and Normand [24]) extended these calculations by taking into account buoyancy effects and thermal transfer properties at the interface. Daviaud and Vince [25] reported the first observation of hydrothermal waves in a shallow layer of 0.65 cSt silicon oil. These observations were complemented by other authors considering different liquids and geometries (Mukolobwiez et al. [26] and Pelacho et al. [27]). Although there is a numerous number of works has dealt with the study of the liquid sloshing, there remains considerable points are not still understood. This is mainly because sloshing is driven by free surface motion and related to many problems such as fluid–ship interaction, bubble formation, and probabilistic aspects. However, the reader may be referred to a rich series of papers on resonances and liquid sloshing, including the two and three dimensions, by Bridges et al. (see, e.g., Refs. [28–31]). In these works, they have discussed the free oscillation in a nearly square container for the case of standing waves. Also, they are interested in sloshing in shallow water in vessels that are undergoing a general rigid body motion. The motivation of the present work is to explore analytically the non-Newtonian influence together with the heat transfer on the surface wave profile formatted at the free interface of a finite viscoelastic liquid layer within a horizontally oscillated rectangular container. From this point of view, the current paper represents an extension to some previous papers such as [32], that deals with the effect of the viscosity on the sloshing of the free surface of a Newtonian isothermal liquid layer [33], where the problem concerning with the nonlinear sloshing dynamics of an isothermal ideal liquid layer within a horizontally oscillating rectangular tank, and Her´ and Weidman [34] who have presented an analytical and czynski experimental study on the periodic oscillation of free containers driven by an inviscid liquid sloshing where the rectilinear horizontal motion of the containers is assumed to be frictionless. The problem is formulated mathematically in Section 2. The basic state solution of the problem is obtained in Section 3. We then proceed to obtain the solution of the perturbed linearized problem by means of Laplace transforms. Some important limiting cases are investigated at the end of this section. In Section 4, we have presented the procedure of the numerical approach utilized to assign the perturbed quantities in terms of the time parameter and then some numerical computations are presented to demonstrate the effects of the various parameters on the behavior of the system. The concluding remarks are summarized in Section 5. 2. Problem formulation We consider the problem of viscoelastic nature fluid flow described by the Maxwell constitutive relation in a rectangular tank of length l, which undergoes a horizontal oscillation motion. A cartesian coordinate system O –xz is defined so that its origin is located at the center of the mean free surface and z-axis points upwards. The vertical sidewalls are insulated and the base of the tank is maintained at constant temperature T b while, the fluid free surface is subjected to Newtonian cooling with external ambient temperature T a and convective heat transfer coefficient h g . In addition, we will suppose that interfacial tension varies linearly with temperature, which can be expressed as follows [35]

σ = σ0 − Γ ( T − T a ),


M.A. Sirwah / Physics Letters A 378 (2014) 3289–3300

where σ0 is the surface tension at T b which is taken as the reference temperature and Γ = − ∂∂ σT | T = T b is the rate of change of surface tension with temperature which is a positive constant for most common fluids. The assumption of linear variation of surface tension with temperature is very much compatible with the experimental data. To account for the temperature-dependence of interfacial tension in the interfacial stress conditions, the temperature distributions at the interfaces must be determined. This requires the solution of the energy equation, which will be done later. The two-dimensional fundamental equations that govern and describe the behavior of Maxwell fluid consist of the continuity equation (for an incompressible fluid with density ρ ), momentum equation and energy equation for the temperature distribution (neglecting viscous dissipation effects) are respectively written as [36]:

∇. u = 0,  du ρ = −∇ pˆ + ∇.,


ρC p

dT dt

The tangential stresses are also balanced across the interface:

..t = ∇ σ ( T ).t , n

(9)  ∂η ∂η where t = (1, ∂ x )/ ( ∂ x )2 + 1 is the corresponding unit tangent vector. Newton’s law of cooling at the free surface that expresses the continuity of heat flux:

−kn.∇ T = h g ( T − T a ),

(ii) On the boundaries of the tank On the solid impermeable boundaries, we adopt the condition which assumes that the normal velocity components of the fluid particle and the body surface are the same and the frictional force is zero [38]. Consequently, we have

∂u ∂ v + = 0, at x = ±l/2, ∂z ∂x ∂u ∂ v + = 0, at z = −h, ∂z ∂x

u = U (t ),


= α∇ T ,


 = (u , v ) is the fluid velocity vector, where u

d dt

∂ ∂t

+ u .∇ stands


where k is the thermal conductivity.




v = 0,

(11) (12)

for the convective derivative, ∂∂t denotes the partial derivative with

respect to the time, ∇ ≡ and pˆ = p + ρ gz is the total hydrostatic pressure.  is the extra stress tensor of an incompressible Maxwell fluid given by the constitutive relationship

where U (t ) is the horizontal velocity of the tank defined by U (t ) = U 0 sin ωt. In this work, we assume that the motion starts smoothly, and then U (0) = 0. The thermal boundary conditions are adiabatic sidewalls, isothermal tank base, viz.

∂ = μΥ˙ , (5) ∂t where μ is the dynamical viscosity coefficient, λ denotes the re-

∂T = 0, at x = ±l/2, ∂x T = T b , at z = −h.

( ∂∂x , ∂∂z ) refers to the gradient operator,


 + (∇ u ) is the rate of strain laxation time parameter, and Υ˙ = ∇ u tensor, with Tr refers to the transpose. Then, Eq. (5) will give Tr

   ∂u j ∂ ∂ ui Πi j = μ , 1+λ + ∂t ∂xj ∂ xi


where Πi j denote the components of the extra stress tensor. In Eq. (4), T , α , and C p denote respectively the absolute temperature, thermal conductivity and the specific heat at constant pressure of the fluid. In order to close the system of the previous governing equations, the boundary conditions at the tank base z = −h, side walls x = ±l/2, and at the fluid free interface z = η(x, t ) are considered, where η is the height of the disturbed interface measured from the initial level z = 0, defined in the next section. The appropriate boundary conditions of the current problem are [36,37]: (i) At the free surface z = η(x, t ) The kinematic condition expresses the fact that the interface always comprises the same fluid particles (i.e. there is no cavitation in the interface between the fluids), and therefore the function η(x, t ) must satisfy the condition

∂η ∂η v= +u . ∂t ∂x


(8) ∂η

 = (− ∂ x , 1)/ ( ∂ x )2 + 1 is the where I is the identity tensor, n exterior pointing normal unit vector to the interface and 2H = ∂2η /(( ∂∂ηx )2 ∂ x2


x∗ , z∗ , h∗ = (x, z, h)/l,

+ 1)3/2 is the mean surface curvature.

u ∗ , v ∗ = (u , v )/ lg ,

Πi∗j = Πi j /ρ gl,

T∗ =

t t∗ =  , l/ g p p∗ = , ρ gl T − Ta Tb − Ta



and then we omit the star for simplicity where it is desirable to use the forms of the same physical parameters. Upon using these scales, the dimensionless momentum and heat equations, are thereby expressed as follows:

  ∂ pˆ ∂Πzx 1 ∂Πxx , + + dt ∂ x Re ∂ x ∂z   ∂ pˆ ∂Πzz 1 ∂Πxz dv , =− + + dt ∂z Re ∂x ∂z



.(− p I + ). n n = 2H σ ( T ),


In what follows, we will nondimensionalize the governing equations and the boundary conditions. In fact, we will often use dimensionless variables to provide improved insight into the physics and in order to understand hydrodynamic stability better. We assume that the stared variables are dimensionless


The balance of the normal stress (the dynamical boundary condition) yields:



dT dt

(16) (17)

= ∇2 T ,

where Re = ρ

(18) gl3 /μ is the gravitational Reynolds number, Pe = μC

Re Pr is the Peclet number, and Pr = α p denotes the Prandtl number. It represents the ratio of viscous diffusivity to thermal diffusivity. The Prandtl number is strictly a function of the properties of the fluid and not on flow properties and plays a strong role in determining the ratio of the thermal boundary layer to the momentum boundary layer.


M.A. Sirwah / Physics Letters A 378 (2014) 3289–3300

When the interfacial displacement is small, the motion may be described by the linearized equations together with the corresponding appropriate linearized boundary conditions. Consequently, the boundary conditions on perturbation interfacial variables need to be evaluated at the equilibrium position rather than at the interface. Therefore, it is necessary to express all the physical quantities involved in terms of Maclaurin series about z = 0. 3. Solution of the problem

In the undisturbed state, the fluid flow as well as the temperature field is assumed to be fully determined where the velocity field and the heat distribution are decoupled from each other. The base state solutions are obtained by assuming a quiescent initial state, therefore the base state velocity in the fluid domain is zero in which the flow is steady and fully developed. Thereby, the unknown are the basic heat distribution in the bulk of the fluid layer and pressure, which are varying linearly with the normal z. The equation for the non-dimensional temperature T (energy equation) in the undisturbed case reads:


= 0,


associated with the thermal boundary conditions

dT 0

= −B T 0,

dz T 0 = 1,

at z = 0,


at z = −h,


where B = h g l/k is the Biot number. Eq. (19) along with the conditions (20) and (21) gives the basic temperature field

T 0 ( z) =

1 − Bz 1 + Bh

∂ψ ∂η =− , at z = 0. ∂t ∂x ∂ T˜ = 0, at z = 0, Π˜ zx + Ma ∂x



and the hydrostatic pressure is given by

p 0 = − z.


(28) (29)

where Ma = Γ ( T b − T a )/μ lg is the Marangoni number.

∂ T˜ dT 0 + B T˜ + B η = 0, ∂z dz

3.1. Solution of the basic state

d2 T 0 ( z )

along with the boundary equations

˜ zz + p˜ − η − Π

1 B0

at z = 0.

(1 − CaΛ)


∂ 2η = 0, ∂ x2

at z = 0,


where Λ = 1+1Bh , B 0 = ρ g ρ l2 /σ0 is the Bond number, and Ca = Γ ( T b − T a )/σ0 is the capillary number.

T˜ = 0,

at z = −h.


∂ 2ψ ∂ 2ψ ∂ψ = U (t ), − = 0, ∂z ∂ z2 ∂ x2 ∂ T˜ = 0, at x = ±1/2. ∂x ∂ 2ψ ∂ 2ψ ∂ψ = 0, − = 0, at z = −h. ∂x ∂ z2 ∂ x2

(33) (34)

We now introduce Laplace transform with respect to time of the general form


  ¯f (x, z, s) = L f (x, z, t ) =

f (x, z, t )e −st dt ,



with s being the transform parameter. Then, we have

∂ 2 ψ¯ , (1 + λs) ∂ x∂ z  2  ∂ 2 ψ¯ ∂ ψ¯ 1 − 2 . Π¯˜ xz = Π¯˜ zx = 2 (1 + λs) ∂ z ∂x

Π¯˜ xx =

∂ 2 ψ¯ , (1 + λs) ∂ x∂ z 2

Π¯˜ zz = −


Accordingly, the governing equations transformed into 3.2. Solution of the perturbed flow

˜ , p˜ , Π˜ i j , T˜ ) for velocity, presIntroducing small disturbances (u sure, and temperature fields

 v˜ (x, z, t ) ,

u = u˜ (x, z, t ),

Πi j = Π˜ i j (x, z, t ),

p = p 0 ( z) + p˜ (x, z, t ),

T = T 0 ( z) + T˜ (x, z, t ).

∂ψ , ∂z

v˜ = −

∂ψ . ∂x



The equations of motion together with the relevant boundary conditions will be solved for these perturbations under the assumption that the perturbations are small, that is, all equations of motion and boundary conditions will be linearized in the perturbed quantities. Such quantities are governed by the equations

˜ xx ∂ 2 Π˜ zx ∂ 2 Π˜ xz ∂ 2  1 ∂ 2Π ∇ ψ = + − ∂t Re ∂ x∂ z ∂ z2 ∂ x2  ˜  dT 0 ∂ψ ∂T Pe = ∇ 2 T˜ , − ∂t dz ∂ x

 ∂ 2 Π˜ zz , − ∂ x∂ z

1 Re

¯ x, z, s), ∇ 4 ψ(

Pe s T¯ 1 (x, z, s) = ∇ 2 T¯ 1 (x, z, s),

(36) (37)

associated with the boundary conditions

The solution of the above system of governing equations and boundary conditions can be facilitated by defining a stream function, ψ(x, z, t ), which automatically satisfies Eq. (2) and is related to the velocity perturbations by

u˜ =

¯ x, z, s) = s(1 + λs)∇ 2 ψ(

(26) (27)

 ∂ 2 ψ¯ ∂ T¯ 1 ∂ 2 ψ¯ + Ma(1 + λs) − = 0, at z = 0. 2 2 ∂x ∂z ∂x ∂ψ B2 ∂ T¯ 1 1 + B T¯ 1 = − , at z = 0. ∂z s (1 + Bh) ∂ x   ∂ 3 ψ¯ ∂ ψ¯ 1 ∂ 3 ψ¯ 1 ∂ 2 ψ¯ −s + 3 2 + 3 + ∂z Re(1 + λs) s ∂ x2 ∂x ∂z ∂z 4 ¯ ∂ ψ 1 + (CaΛ − 1) 4 = 0, at z = 0. B0s ∂x ¯T 1 = 0, at z = −h. ∂ T¯ 1 ∂ ψ¯ = U¯ (s), = 0, at x = ±1/2. ∂z ∂x ∂ 2 ψ¯ ∂ 2 ψ¯ ∂ ψ¯ = 0, − = 0, at z = −h. ∂x ∂ z2 ∂ x2

(38) (39)

(40) (41) (42) (43)

A separable solution satisfying sidewall boundary conditions (42) is

M.A. Sirwah / Physics Letters A 378 (2014) 3289–3300

T¯ 1 = g¯ ( z, s) cos αn (x + 1/2),

In this case, we have the solution

ψ¯ = ¯f ( z, s) sin αn (x + 1/2) + U¯ z,


where αn = nπ and n is an integer. Consequently the functions ¯f , g¯ , are governed by the following equations

d ¯f 4


   d ¯f  − 2αn2 + Re s(1 + λs) + αn2 αn2 + Re s(1 + λs) ¯f = 0, 2 2


(45) Pe(s g¯ − αn ¯f ) =

d2 g¯ dz2

2 ¯. ng


ψ¯ =

∞ 2(−1 + (−1)n ) B 0 U¯ s(−2αn2 e γn z + (αn2 + γn2 )e αn z )

αn3 ( B 0 + αn2 )Re(1 + λs)

n =1

× sin αn (x + 1/2) + U¯ z.

ψ(x, z, t ) =

where γn = αn2 + Re s(1 + λs), δn = αn2 + Pe s. Thereby, Eqs. (36) and (37) will have the solutions

¯ x, z, s) = ψ(

T¯ (x, z, s) =



The coefficients c 1 , c 2 , ..., c 6 are given in Appendix A. In what follows, it is of interest to investigate some important limiting cases: (i) The case of infinite depth and isothermal flow We consider the simpler case of an infinitely deep container when the fluid is assumed to be isothermal and there is no heat transfer between the phase fluid and the ambient media. The equations of motion together with the appropriate boundary conditions are given by:

¯ x, z, s) = s(1 + λs)∇ 2 ψ(

1 Re

¯ x, z, s), ∇ 4 ψ(

2(−1 + (−1)n ) B 0 U 0 ωe αn z

αn3 ( B 0 + αn2 )Re(1 + λ2 ω2 )


(αn2 − δn2 )(δn2 − γn2 ) n =1   × αn2 − δn2 (c 1 cosh γn z + c 2 sinh γn z)   + γn2 − δn2 (c 3 cosh αn z + c 4 sinh αn z) + c 5 cosh δn z  1 − Bz + c 6 sinh δn z cos αn (x + 1/2) + , (50) 1 + Bh ∞ αn (c1 + c3 ) η¯ (x, s) = − (51) cos αn (x + 1/2). n =1


   × Re 1 + λ2 ω2 δ(t ) + 2αn2 cos ωt − e −t /λ    − ω Re − 2αn2 λ + Re λ2 ω2 sin ωt +

4(−1 + (−1)n ) B 0 ω U 0 e αn z 2αn2 Re ( B 0 + αn2 )(1 + λ2 ω2 )

  × −2αn e −t /λ + Re δ(t ) 1 + λ2 ω2 z + 2αn cos ωt    + 2αn λω − Re ω z − Re λ2 ω3 z sin ωt 

× sin αn (x + 1/2) + U z, (59) ∞ n 2(−1 + (−1) ) B 0 U 0 ω η(x, t ) = 2 α ( B 0 + αn2 )Re(1 + λ2 ω2 ) n =1 n   × 2αn2 cos ωt + λω sin ωt − e −t /λ     + ω Re − 2αn2 λ + Re λ2 ω2 sin ωt − Re 1 + λ2 ω2 δ(t )   (60) + 2αn2 e −t /λ − cos ωt cos αn (x + 1/2), where δ(t ) denotes the Dirac delta function. (ii) The case of infinite depth and inviscid liquid

¯ x, z, s) = 0, along with In such a case, we solve the equation ∇ 2 ψ( the boundary conditions (54)–(56) and then the problem will admit the following solutions for the stream function and surface elevation


∂ 2 ψ¯ ∂ 2 ψ¯ − = 0, at z = 0, ∂ z2 ∂ x2   ∂ 3 ψ¯ ∂ ψ¯ 1 ∂ 3 ψ¯ −s + 3 2 + 3 ∂z Re(1 + λs) ∂x ∂z ∂z 1 ∂ 2 ψ¯ 1 ∂ 4 ψ¯ − = 0, at z = 0, + s ∂ x2 B 0 s ∂ x4 ∂ ψ¯ = U¯ (s), at x = ±1/2, ∂z ψ¯ → 0, as z → −∞.

Re 1 + λ2 ω2 δ(t )

where L−1 stands for the inverse Laplace transform. When the Reynolds number is vanishingly small Re  1, then s(1+λs) z e γn z = e αn z (1 + 2α Re) + O (Re2 ). In this limit, the stream funcn tion together with the surface wave elevation are then given by

n =1

n =1

+ c 4 sinh αn z) sin αn (x + 1/2) + U¯ z, ∞ a Pe αn

αn3 ( B 0 + αn2 )Re(1 + λ2 ω2 )

  + 2αn2 cos ωt − e −t /λ    − ω Re − 2αn2 λ + Re λ2 ω2 sin ωt   4(−1 + (−1)n ) B 0 −1 se γn z − L (1 + λs)(s2 + ω2 ) αn ( B 0 + αn2 ) Re (58) × sin αn (x + 1/2) + U (t ) z,


∞ (c 1 cosh γn z + c 2 sinh γn z + c 3 cosh αn z

2(−1 + (−1)n ) B 0 U 0 ωe αn z   n =1


¯f = c 1 cosh γn z + c 2 sinh γn z + c 3 cosh αn z + c 4 sinh αn z, (47)   2 a Pe αn g¯ = αn − δn2 (c1 cosh γn z + c2 sinh γn z) (αn2 − δn2 )(δn2 − γn2 )   + γn2 − δn2 (c 3 cosh αn z + c 4 sinh αn z)c 5 cosh δn z  + c 6 sinh δn z , (48)


Consequently, the stream function will be expressed as a function of the time in the following form:

In the light of the boundary conditions (38)–(41) and (43), the functions ¯f , g¯ are given by




ψ(x, z, t ) =

∞ 2(−1 + (−1)n ) B 0 U 0 ωe αn z (δ(t ) − ω sin ωt ) n =1

αn3 ( B 0 + αn2 )

× sin αn (x + 1/2) + U (t ) z, (61) ∞ n 2(−1 + (−1) ) B 0 U 0 ω cos ωt η(x, t ) = − cos αn (x + 1/2). αn2 ( B 0 + αn2 ) n =1 (62)

(55) (56)

Such results are in well agreements with those obtained by Forbes ´ [33] and Herczynski and Weidman [34]. However, extensive results for inviscid flow are presented in various publications [39,40].


M.A. Sirwah / Physics Letters A 378 (2014) 3289–3300

4. Numerical procedure 4.1. Inversion of Laplace transforms It is too difficult to get the analytical inverse Laplace transform of the complicated obtained solutions for the stream function, the temperature field distribution, and the surface waves displacement, except for some limiting cases as depicted before. We then must resort to numerical calculations. The complex inversion formula for Laplace transforms can be written as [41]

f (x, z, t ) =

υ +i ∞

1 2π i

¯f (x, z, s)e st ds,


υ −i ∞

where ¯f (x, z, s) can be any of the transformed field quantities and υ is a positive free parameter that must be greater than the real parts of all singularities of ¯f (x, z, s). To obtain the unknown functions in the physical domain, we adopt the numerical inversion method of Durbin’s [42] based on a Fourier series expansion. In this method, the function f (x, z, t ) is approximated by the formula [43]

f (x, z, t ) =

eυ t

τ +

Fig. 1. Represents evolution of the free-surface profile with different values of K for a system having Ma = 1, Pe = 2, Ca = 0.05, B = 0.5, B 0 = 0.7, λ = 0.3, U 0 = 0.3, ω = 2, h = 0.3: Re = 5.

 1  − e ¯f (x, z, υ ) 2


   kπ kπ e ¯f x, z, υ + i cos t


k =0


    kπ kπ ¯ − m f x, z, υ + i sin t ,



0 < t < 2τ , (64)

where e {.} stands for the real part of {.}; while m{.} denotes the imaginary part. The function f refers to ψ , T , or η . K is a sufficiently large integer representing the number of terms in the truncated infinite Fourier series. K must be chosen such that


e υ t e ¯f

x, z, υ + i





    kπ kπ − m ¯f x, z, υ + i sin t ≤ ,




where  is a prescribed small positive number that corresponds to the degree of accuracy to be achieved. The optimal choice of υ was obtained according to the criteria described in [44]. The suggested value of υτ is between 5 and 10 for sufficient accuracy, with K ranging from 50 to 5000 [42]. Also, the solution interval is generally valid for τ ≥ 2t max where t max is the time up to which the results are to be calculated [45]. 4.2. Numerical results A number of numerical tests was achieved and it has been revealed that by choosing the following numerical values υτ = 9, K = 700, τ = 2t max = 200, Durbin inversion formula (64) provides stable and convergent results in all cases considered. The computations were performed on a cluster of core-i7 based personal computers, and the convergence of numerical solutions was established in a simple trial and error manner, by increasing the number of samples in Laplace domain (the truncation constant K ), while searching for stability in the numerical values of the computed solutions. The truncation constant K = 700 was found to led to uniform convergence for all configurations considered as inspected in Fig. 1. Fig. 2 shows the effect of increasing Reynolds Re, on the free surface time history, represented by Eq. (51). So, graphical

Fig. 2. Free-surface profile as a function of time corresponding to (51), for a system having Ma = 1, Pe = 2, Ca = 0.05, B = 0.5, B 0 = 0.7, λ = 0.3, U 0 = 0.3, ω = 2, h = 0.3: (a) Re = 5, (b) Re = 100, (c) Re = 500.

M.A. Sirwah / Physics Letters A 378 (2014) 3289–3300


Fig. 4. Influence of the relaxation time parameter on the wave elevation history for a system with parameters as those of Fig. 2, but with Re = 100 and for λ = 0.01, 0.5, 0.9.

Fig. 3. Three-dimensional spatial evolution of free-surface deformation for the same system as considered in Fig. 2.

representations of the numerical computations correspond to several values of Re are displayed, while maintaining all the remaining dimensionless parameters constant as given in the caption to Fig. 2. Fig. 2(a), where Re = 5, displays linear behavior of the free surface with irregular patterns for the short time history while the wave motions vary in a periodic manner with constant amplitude as the time goes on. Figs. 2(b), 2(c) depict situations in which all of the dimensionless parameters except Re, that has the values 100 and 500 respectively, are identical to those in Fig. 2(a). The results shown in Figs. 2(b), 2(c) clarify that there is no qualitative change in the behavior of the surface waves with the evolution in time at the different values of the Reynolds number, nevertheless a noticeable increase in values of η is observed. All one can say in

connection with Fig. 1 is that an increase in Re yields an overall destabilizing effect in the sense that it will lead to a wave motion disorder. This indicates that the liquid viscosity coefficient has a stabilizing influence for the selected values of input parameters due to the increased viscous dissipation. Such results are in a good agreement with those provided by Cerda and Tirapegui [46]. To illustrate a more clear depiction of the wave motion, the influence of increasing value of Reynolds number on the liquid free surface in three-dimensional space interface, correspond to the cases of Fig. 2, is demonstrated in the parts of Fig. 3. The numerical results displayed in Fig. 4 aim to investigate the effect of increasing relaxation time parameter λ on the surface wave profile with the evolution in time. Inspection of the graphs depicted in this figure clearly shows that the increase in the values of λ results in a significant increase in values of η , particularly when the value of t is less 20, while this change becomes much less at the large values of t (greater than 20).Consequently, the disturbing effect of the viscoelastic parameter λ on the motion of the free surface will diminish gradually with the evolution in time. On the other hand, it is observed at the same position inside the basin that the maximum and minimum values of the surface wave elevation η arise at different times dues to the variation of λ. This means that the motion of the viscoelastic liquid layer in a bounded region is less stable than the Newtonian one. However, viscoelastic forces try to destabilize the fluid motion. Fig. 5 displays the influence of the dimensionless thickness h of the fluid layer upon the free surface sloshing, where three values are assigned to h (= 0.1, 0.4, 0.8) for the sake of comparison. As the curves displayed in Fig. 5 are compared to each another, it is clarified that the higher the fluid thickness, the larger surface wave elevation, especially at the small values of t. Then, when we go towards the shallow water model, the flow pattern becomes more stable in bounded regions as reported in the lateral extent regions. Moreover, we observe that the positions of nodes are unchanged dues to the variation of h and that only the wave elevation changes up and down, which agrees with the main property of the linear standing waves. When the surface waves reflect off of the wall tank, it can interfere with its own reflection and then this creates a standing wave. Furthermore, the corresponding numerical wave profiles across the tank at different times are also shown in Fig. 6. The waveforms presented in Fig. 6 show wave propagation that has quasi-periodic waveform along the direction of propagation with a cross-tank variation from maxima on one wall and minima on the other, i.e. it has a time periodic half-wavelength variation across the tank. The


M.A. Sirwah / Physics Letters A 378 (2014) 3289–3300

Fig. 5. Represents the effect of the layer fluid thickness on the motion of the free surface, where Re = 100 and h = 0.1, 0.4, 0.8. The other parameters have the same values as considered in Fig. 2.

linear standing waves formed on the surface of the liquid moves up one side of the tank and down the other; then the up half-wave moves down and the down half-wave moves up, and so on. This figure clarifies that the wave has zero amplitude at t = 0, a positive peak at t = 10 at the left wall of the tank and a negative peak at the right wall. This is the fundamental antisymmetric wave. In Fig. 7, we examine the influence of the angular frequency of the periodic horizontal velocity of the tank on the long-time history of the surface waves deflection. It is to be noted that in-

creasing the angular frequency results in a notable growth in the wave amplitude and then the angular frequency has a destabilizing effect. This is due to increasing the kinetic energy transferred from the mechanical motion of the tank to wave motion and thereby the wave motion it goes towards the instability case. Fig. 8(a) summarizes the variation of the temperature distribution inside the basin with the vertical position z at various times and at the horizontal position x = −0.1. Values of the other dimensionless parameters are given in the caption to Fig. 8(a). We note that there is a continuous decrease in the values of the fluid temperature as well as we approach the basin free surface. In addition, such rate of decline in the temperature distribution be very large in the lower half of the liquid layer (−1 < z < −0.5), while decreasing significantly communicates to the positions in the upper half of the liquid layer (−0.5 < z < 0). Moreover, we observe that the evolution in time has a little influence on the heat transfer rate within the liquid layer at the positions close to the free surface. In Fig. 8(b), we repeat the same study carried out in Fig. 8(a), but at the position x = 0.1. We note that the behavior of the temperature distribution is similar to that investigated in Fig. 8(a). However, by comparing the results shown in Figs. 8(a) and 8(b), it is noteworthy that the values of the temperature at the position x = −0.1 is greater than its counterpart at the position x = 0.1, especially at large values of time. This means that the heat distribution at the same height varies from one position to another in spite of being regular in the basic case, which depends on the height only.

Fig. 6. Variation of the free-surface profile with the horizontal position for the same system as investigated in Fig. 2, with Re = 100, h = 0.4, corresponding to the times t = 0, 10, 20, 30, 40. Fig. 7. Evolution of the surface-wave profile with time at different values of the angular frequency ω , for the same system as considered in Fig. 2, with Re = 5.

Fig. 8. Temperature distribution inside the tank with the vertical position, z for a system having Ma = 1, Pe = 2, Ca = 0.05, B = 0.5, B 0 = 0.7, Re = 100, λ = 0.1, U 0 = 0.3, ω = 2, h = 0.7, corresponding to t = 1, 5, 10, 20: (a) x = −0.1, (b) x = 0.1.

M.A. Sirwah / Physics Letters A 378 (2014) 3289–3300

Fig. 9. Evolution of the temperature distribution with time for the same system as considered in Fig. 8, where x = 0.1 at different vertical locations, z = −0.9, −0.5, −0.3, −0.05.

Fig. 9 shows the variation of the thermal distribution within the liquid layer with time at different heights. Investigation of the numerical results presented in Fig. 9 reveals that the fluid temperature distribution at any position increases with time until it reaches the maximum value (the heat conduction inside the fluid particles is much faster than the heat convection away from it) and then starts decreasing gradually with increasing time (the heat conduction inside the body is much slower than the heat convection away from it). Also it is found that the maximum value of the fluid temperature increases gradually when the positions closest to the base of the container. This clarifies that the fluid becomes warmer downward and the thermocapillary force appearing at the interface acts in the opposite direction to the gravitational acceleration. The purpose of the numerical results shown in Fig. 10 is to explore the influence of the Biot number B, on the behavior of the temperature distribution with the evolution in time at various heights. Values of all the dimensionless parameters are given in the


caption of Fig. 9. In Fig. 10(a), where z = −0.9, we observe that increasing value of Biot number increases the maximum value of the fluid temperature. This indicates that the effect of the Biot number is to advance the heat transfer. Fig. 10(b) describes a situation in which all of the parameters are the same as those in Fig. 10(a) except z which equals −0.7. Fig. 10(b) shows that the parameter B performs the same role that was observed in Fig. 10(a) except that the maximum values of the temperature are lower than its counterpart in Fig. 10(a). As the Biot number increases, heat transfer occurs more rapidly between the fluid particles. Similar results for the influence of the Biot number are reported by many investigators [47,48]. As for the Fig. 10(c), we investigate the influence of B at z = −0.5 corresponds to the positions located at the midheight of the liquid layer. It should be noted that the parameter B has a very slight effect on the heat flow field. The results shown in Fig. 10(d) highlight the effect of the Biot number on the heat transfer at the position z = −0.2 located at the upper half of the liquid layer. In contrast to the results demonstrated in Figs. 10(a) and 10(b), it is observed that increasing B leads to a decrease in the maximum values of the fluid temperature at the positions placed above the initial mid-height, indicating that the effect of the Biot number is to retard the heat transfer throughout the upper half of fluid layer. By comparing the results displayed throughout the parts of Fig. 10, it is remarkable that the effect of increasing the Biot number on the fluid temperature flow depends strongly on the point position. To be better understand the flow behavior, we also present the streamlines in Fig. 11 and the corresponding isotherms in Fig. 12. The streamlines are effective tools to visualize a qualitative impression of the flow behavior during the motion. Part (a) depicts the numerical results of instantaneous streamlines of the stream function at Re = 3, while in parts (b) and (c), we take Re = 100 and 500, respectively, for the sake of comparison. We notice that a vortex appears whose center is close to the free surface (Fig. 11(a)) that will penetrate into deeper depth of the liq-

Fig. 10. Shows the influence of the Biot number on the time evolution of the temperature field inside the container as considered in Fig. 9: (a) z = −0.9, (b) z = −0.7, (c) z = −0.5, (d) z = −0.2.


M.A. Sirwah / Physics Letters A 378 (2014) 3289–3300

Fig. 11. Represents the streamlines contours for a system having the same parameters considered in Fig. 8: (a) Re = 3, (b) Re = 100, (c) Re = 500.

Fig. 12. Displays the isothermal contours for the same system as considered in Fig. 11.

uid with the increasing Reynolds number as shown in Figs. 11(a) and 11(c).

effect on the flow pattern and consequently the free surface motion be expected to be more perturbed and goes away from the regularity of the linear standing waves [46]. Such effects are also observed for the relaxation time parameter on the behavior of the free surface deformations, especially at the times close to the initial state, while this influence will lessen when the time goes on more. In other words, the change in the surface elevation is small as well as we go away from the initial state. The numerical results also showed that the rate of heat flow within the container is larger whenever the position is close to the base of the container, which represents the heat source in agreement with the results of Ezzat [49]. Further, we observe that monotonic increase in the values of the Biot number results in the augmenting the maximum value of the fluid temperature and this influence will strictly depend on the location inside the fluid layer. Consequently the increase in the Biot number trends to advance the hear transfer.

5. Conclusions A theoretical analysis of two-dimensional horizontal forcing of a rectangular basin containing a viscoelastic liquid layer when the base of such basin is maintain at a constant temperature has been presented. The problem is formulated in the frame of the linear theory of perturbation, where the perturbed quantities are assumed to be slightly disturbed away from its initial values. The solution of the problem is carried out using Laplace transform that is applied to the equations of motion as well as the relevant boundary equations. The inversion of such transforms is obtained utilizing Durbin’s numerical method. Solutions of special cases involving, infinite depth and isothermal flow, high viscous flow, and case of infinite depth and inviscid liquid are investigated. Numerical results are given and illustrated graphically to demonstrate the effects of the various physical parameters on the behavior of the surface waves pattern, the variation of the thermal distribution within the liquid layer, the streamlines, and isothermals. The results of our study showed that the variation of the Reynolds number has a substantial effect on the surface wave pattern. It was noted that growth of the Re results in a destabilization

Appendix A

c i = c˜ i /c , where

c˜ 1 = − −1 + (−1)n (1 + Bh)2 U s2 (1 + λs)

 2  3  × αn2 − γn2 αn2 − δn2 δn2 − γn2 sinh hγn

M.A. Sirwah / Physics Letters A 378 (2014) 3289–3300

× B Ma Pes(1 + λs)αn cosh hαn sinh hδn    + αn2 sinh hαn B 2s + B (Ma + 2hs + Ma sλ) sinh hδn  + 2(1 + Bh)sδn cosh hδn    + δn sinh hαn − B 2s + B (Ma + 2hs + Masλ) δn sinh hδn   − s cosh hδn B Ma Pe(1 + λs) + 2(1 + Bh)δn2 ,   c˜ 2 = −1 + (−1)n (1 + Bh)2 U s2 (1 + λs)  2  3  × αn2 − γn2 αn2 − δn2 γn2 − δn2 cosh hγn × B Ma Pe s(1 + λs)αn cosh hαn sinh hδn    + αn2 sinh hαn B 2s + B (Ma + 2hs + Masλ) sinh hδn  + 2(1 + β h)sδn cosh hδn    + δn sinh hαn − B 2s + B (Ma + 2hs + Masλ) δn sinh hδn   − s cosh hδn B Ma Pe(1 + λs) + 2(1 + Bh)δn2 ,  1  c˜ 3 = 2 −1 + (−1)n (1 + Bh)2 U s2 (1 + λs)


 3  2  × αn2 − γn2 αn2 − δn2 δn2 − γn2 sinh hαn  × −sδn sinh hγn cosh hδn − B Ma Pe αn2 (1 + λs)    + (1 + Bh) αn2 + γn2 γn2 − δn2  + B sinh hδn −Ma Pe s(1 + λs)αn2 γn cosh hγn   + sinh hγn s + B (Ma + hs + Ma sλ) αn2   + (1 + Bh)sγn2 δn2 − γn2 ,  1  c˜ 4 = − 2 −1 + (−1)n (1 + Bh)2 U s2 (1 + λs)



  2 3

 2 3

× αn2 − γn αn2 − δn δn2 − γn cosh hαn   × (1 + Bh)sγn2 sinh hγn ( B sinh hδn + δn cosh hδn ) γn2 − δn2  + αn2 B Ma Pe s(1 + λs)γn sinh hδn cosh hγn    + γn2 sinh hγn B s + B (Ma + hs + Ma sλ) sinh hδn  + (1 + Bh)sδn cosh hδn    − δn sinh hγn B s + B (Ma + hs + Ma sλ) δn sinh hδn   , + s cosh hδn B Ma Pe(1 + sλ) + (1 + Bh)δn2   1 −1 + (−1)n B (1 + Bh)U s2 (1 + λs) c˜ 5 = −



  2 2

 2 2

× αn2 − γn αn2 − δn δn2 − γn sinh hδn  × −αn2 sinh hαn B Ma Pe(1 + λs)  − 2(1 + Bh)αn2 + 2(1 + Bh)δn2 Pe sγn cosh hγn    + B sinh hγn γn2 − δn2 + Pe s + sinh hγn Pe sαn cosh hαn     + B sinh hαn αn2 − δn2 + Pe s −(1 + Bh)γn2 γn2 − δn2    , + αn2 B Ma Pe(1 + λs) − (1 + Bh) γn2 − δn2     1 −1 + (−1)n B (1 + Bh)U s2 (1 + λs) αn2 − γn2 c˜ 6 = −


2  2  × αn2 − δn2 δn2 − γn2 cosh hδn    × Pe sαn sinh hγn cosh hαn −(1 + Bh)γn2 γn2 − δn2    + αn2 B Ma Pe(1 + λs) + (1 + Bh) δn2 − γn2   + sinh hαn −Pe sαn2 γn cosh hγn B Ma Pe(1 + λs)   + 2(1 + Bh) δn2 − αn2


      + B αn2 − γn2 sinh hγn (1 + Bh) Pe s − δn2 γn2 − δn2    + αn2 P e 2s + B (Ma + 2hs + Ma sλ)   , + (1 + Bh) γn2 − δn2 c = |c i j |, c 11 = BPe αn

αn2 − δn2 cosh hγn ,

c 12 = − BPe αn


αn2 − δn2 sinh hγn , 

c 13 = − BPe αn δn2 − γn2 cosh hγn ,

c 14 = BPe αn δn2 − γn2 sinh hγn ,

c 15 = (1 + Bh) δn2 − γn2

c 16 = −(1 + Bh) δn2 − γn2 c 21 = s

αn2 − δn2 cosh hδn ,

αn2 − δn2 sinh hδn , 

αn2 − δn2 (1 + Bh)γn2 γn2 − δn2

   + αn2 B Ma Pe(1 + λs) + (1 + Bh) γn2 − δn2 ,

c 22 = 0, c 23 = αn2 s

γn2 − δn2 B Ma Pe(1 + λs) + 2(1 + Bh) αn2 − δn2 ,

c 24 = 0,

c 25 = −(1 + Bh) Ma αn s(1 + λs) Q 26 = 0,

c 31 = − B 2 αn (1 + λs)

αn2 − δn2 γn2 − δn2 ,

αn2 − δn2 Pe s + γn2 − δn2 , 

c 32 = − BPe αn γn s(1 + λs)

αn2 − δn2 ,

c 33 = − B 2 αn (1 + λs) Pe s + αn2 − δn2

c 34 = − BPe αn2 s(1 + λs)

γn2 − δn2 ,

γn2 − δn2 , 

c 35 = B (1 + Bh)s(1 + λs)

c 36 = (1 + Bh)s(1 + λs)δn

αn2 − δn2 γn2 − δn2 , αn2 − δn2 γn2 − δn2 , 

c 41 = Re αn3 B 0 + (1 − CaΛ)αn2 (1 + λs),

c 42 = B 0 αn γn s Re s(1 + λs) + 3αn2 − γn2 ,

c 43 = Re αn3 B 0 + (1 − CaΛ)αn2 (1 + λs),

c 44 = B 0 αn2 s Re s(1 + λs) + 2αn2 , c 45 = 0,

c 46 = 0,

c 51 = cosh hγn ,

Q 52 = − sinh hγn ,

c 54 = − sinh hαn , c 61 =

2 n

2 n

α +γ

cosh hγn ,

c 63 = 2αn2 cosh hαn , c 65 = 0,

c 55 = 0,

c 53 = cosh hαn ,

c 56 = 0,

c 62 = −

αn2 + γn2 sinh hγn ,

c 64 = −2αn2 sinh hαn ,

c 66 = 0.

References [1] R.A. Ibrahim, V.N. Pillipchuk, Recent advances in liquid sloshing dynamics, Appl. Mech. Rev. 54 (2) (2001) 133–199. [2] J.C. Virella, C.A. Prato, L.A. Godoy, Linear and nonlinear 2D finite element analysis of sloshing modes and pressures in rectangular tanks subject to horizontal harmonic motions, J. Sound Vib. 312 (2008) 442–460. [3] M. Hermann, A. Timokha, Modal modelling of the nonlinear resonant fluid sloshing in a rectangular tank I: a single-dominant model, Math. Models Methods Appl. Sci. 15 (2005) 1431–1458. [4] J.B. Frandsen, Sloshing motions in excited tanks, J. Comput. Phys. 196 (2004) 53–87. [5] D. Liu, P. Lin, A numerical study of three-dimensional liquid sloshing in tanks, J. Comput. Phys. 227 (8) (2008) 3921–3939. [6] T. Ikeda, R.A. Ibrahim, Nonlinear random responses of a structure parametrically coupled with liquid sloshing in a cylindrical tank, J. Sound Vib. 284 (1–2) (2005) 75–102.


M.A. Sirwah / Physics Letters A 378 (2014) 3289–3300

[7] M. Abbas, Z. Mansour, Sloshing damping in cylindrical liquid storage tanks with baffles, J. Sound Vib. 311 (1–2) (2008) 372–385. [8] H.S. Kim, Y.S. Lee, Optimization design technique for reduction of sloshing by evolutionary methods, J. Mech. Sci. Technol. 22 (2008) 25–33. [9] J.R. Cho, H.W. Lee, Free surface tracking for the accurate time response analysis of nonlinear liquid sloshing, J. Mech. Sci. Technol. 19 (7) (2005) 1517–1525. [10] S. Kamiyama, K. Ueno, Y. Yokota, Numerical analysis of unsteady gas–liquid two-phase flow of magnetic fluids, J. Magn. Magn. Mater. 201 (1–3) (1999) 271–275. [11] V.G. Bashtovoy, O.A. Lavrova, V.K. Polevikov, L. Tobiska, Computer modeling of the instability of a horizontal magnetic-fluids layer in a uniform magnetic field, J. Magn. Magn. Mater. 252 (2002) 299–301. [12] R. Badescu, D. Condurach, M. Ivanoiu, Ferrofluid with modified stabilisant, J. Magn. Magn. Mater. 202 (1) (1999) 197–202. [13] M.J. Lighthill, Waves in Fluids, CUP, Cambridge, 1978. [14] T. Sarpkaya, M. Issacson, Mechanics of Wave Forces on Offshore Structures, Van Nostrand, New York, 1981. [15] R.W. Yeung, X. Yu, Wave-structure interaction in a viscous fluid, in: Fourteenth Int. Conf. on Offshore Mech. and Arctic Eng., Copenhagen, Denmark, 1995. [16] C.C. Mei, The Applied Dynamics of Ocean Surface Waves, Wiley–Interscience, New York, 1983. [17] H. Lamb, Hydrodynamics, 6th edition, CUP, Cambridge, 1976. [18] Q. Chen, E. Rame, S. Garoff, The velocity field near a moving contact line, J. Fluid Mech. 337 (1997) 49–66. [19] E.A. Cerda, E.L. Tirapegui, Faraday’s instability in viscous fluid, J. Fluid Mech. 368 (1998) 195–228. [20] C.Y. Loh, H. Rasmussen, A numerical procedure for viscous free surface flow, Appl. Numer. Math. 3 (1987) 479–495. [21] M. Smith, S. Davis, Instabilities of dynamic thermocapillary liquid layers. Part 1. Convective instabilities, J. Fluid Mech. 132 (1983) 145–162. [22] M. Smith, Instability mechanisms in dynamic thermocapillary liquid layers, Phys. Fluids 29 (1986) 3182–3186. [23] P.M. Parmentier, V. Regnier, G. Lebon, Bouyant-thermocapillary instabilities in medium Pr number fluid layers subject to a horizontal temperature gradient, Int. J. Heat Mass Transf. 36 (1993) 2417–2427. [24] J. Mercier, C. Normand, Buoyant thermocapillary instabilities of differentially heated liquid layers, Phys. Fluids 8 (1996) 1433–1445. [25] F. Daviaud, J. Vince, Traveling waves in a fluid layer subjected to a horizontal temperature gradient, Phys. Rev. E 48 (1993) 4432–4436. [26] N. Mukolobwiez, A. Chiffaudel, F. Daviaud, Supercritical Eckhaus instability for surface-tension-driven hydrothermal waves, Phys. Rev. Lett. 80 (1998) 46–61. [27] M.A. Pelacho, A. Garcimartyn, J. Burguete, Local Marangoni number at the onset of hydrothermal waves, Phys. Rev. E 62 (2000) 477–483. [28] T.J. Bridges, Secondary bifurcation and change of type of three-dimensional standing waves in a finite depth, J. Fluid Mech. 179 (1987) 137–153. [29] H.A. Ardakani, T.J. Bridges, Shallow-water sloshing in rotating vessels: details of the numerical algorithm, Tech. Rep., Department of Mathematics University of Surrey, 2009.

[30] H.A. Ardakani, T.J. Bridges, Shallow-water sloshing in vessels undergoing prescribed rigid-body motion in two dimensions, Eur. J. Mech. B, Fluids 31 (2012) 30–43. [31] M.R. Turner, T.J. Bridges, Nonlinear energy transfer between fluid sloshing and vessel motion, J. Fluid Mech. 719 (2013) 606–636. [32] G.X. Wu, R.E. Taylor, D.M. Greaves, The effect of viscosity on the transient freesurface waves in a two-dimensional tank, J. Eng. Math. 40 (2001) 77–90. [33] L.K. Forbes, Sloshing of an ideal fluid in a horizontally forced rectangular tank, J. Eng. Math. 66 (2010) 395–412. ´ [34] A. Herczynski, P.D. Weidman, Experiments on the periodic oscillation of free containers driven by liquid sloshing, J. Fluid Mech. 693 (2012) 216–242. [35] A. Samanta, Stability of liquid film falling down a vertical non-uniformly heated wall, Physica D 237 (2008) 2587–2598. [36] K. Zakaria, M.A. Sirwah, S. Alkharashi, A two-layer model for superposed electrified Maxwell fluids in the presence of heat transfer, Commun. Theor. Phys. 55 (2011) 1077–1094. [37] M.A. Sirwah, Nonlinear dynamics and chaos of magnetized resonant surface waves in a rectangular container, Wave Motion 50 (2013) 596–618. [38] G.X. Wu, R.E. Taylor, D.M. Greaves, The effect of viscosity on the transient freesurface waves in a two-dimensional tank, J. Eng. Math. 40 (2001) 77–90. [39] G.X. Wu, Q.W. Ma, R.E. Taylor, Numerical simulation of sloshing waves in a 3D tank based on a finite element method, Appl. Ocean Res. 20 (1998) 337–355. [40] G.X. Wu, Second-order resonance of sloshing in a tank, Ocean Eng. 34 (2007) 2345–2349. [41] F.B. Hilderbrand, Methods of Applied Mathematics, Prentice Hall of India, 1992. [42] F. Durbin, Numerical inversion of Laplace transforms: an effective improvement of Dubner and Abate’s method, Comput. J. 17 (1973) 371–376. [43] S.C. Fan, S.M. Li, G.Y. Yu, Dynamic fluid–structure interaction analysis using boundary finite element method-finite element method, J. Appl. Mech. 72 (2005) 591–598. [44] G. Honig, U. Hirdes, A method for the numerical inversion of Laplace Transforms, J. Comput. Appl. Math. 10 (1984) 113–132. [45] Y.-C. Su, C.-C. Ma, Transient wave analysis of a cantilever Timoshenko beam subjected to impact loading by Laplace transform and normal mode methods, Int. J. Solids Struct. 49 (2012) 1158–1176. [46] E.A. Cerda, E.L. Tirapegui, Faraday’s instability in viscous fluid, J. Fluid Mech. 368 (1998) 195–228. [47] M.R. Sadiq, R. Usha, Linear instability in a thin viscoelastic liquid film on an inclined, non-uniformly heated wall, Int. J. Eng. Sci. 43 (2005) 1435–1449. [48] O.D. Makindea, T. Chinyoka, L. Rundoraa, Unsteady flow of a reactive variable viscosity non-Newtonian fluid through a porous saturated medium with asymmetric convective boundary conditions, Comput. Math. Appl. 62 (2011) 3343–3352. [49] M.A. Ezzat, Free convection effects on perfectly conducting fluid, Int. J. Eng. Sci. 39 (2001) 799.