ARTICLE IN PRESS JID: CAF [m5G;August 14, 2018;18:50] Computers and Fluids 0 0 0 (2018) 1–21 Contents lists available at ScienceDirect Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid Long duration SPH simulations of sloshing in tanks with a low ﬁll ratio and high stretching Mashy D. Green∗, Joaquim Peiró∗ Department of Aeronautics, Imperial College London, United Kingdom a r t i c l e i n f o Article history: Received 23 May 2017 Revised 10 July 2018 Accepted 11 July 2018 Available online xxx Keywords: SPH Sloshing GPU a b s t r a c t Achieving stable and accurate simulations of long duration sloshing with low ﬁll ratios using smooth particle hydrodynamics (SPH) is a challenging problem. Its solution requires a judicious choice of SPH formulation to minimize the effect of errors introduced into the simulation by boundary conditions, dissipation terms or computer arithmetic. We assess the diﬃculties and common pitfalls of such simulations and propose a SPH method to deal with them effectively. The formulation combines an eﬃcient use of double precision, wall boundary conditions using ghost particles, and a δ -SPH scheme for stability with minimal energy dissipation. The proposed SPH formulation accurately matches the experimental data, both in terms of surface elevations and forces on the tank, of the study of tuned liquid dampers by Reed et al. 1998, later reproduced by ESA/ESTEC, over a wide range of frequencies. © 2018 Elsevier Ltd. All rights reserved. 1. Introduction SPH is a fully Lagrangian particle method for the solution of partial differential equations. Their discretisation uses particles which follow the ﬂuid ﬂow and the relevant properties are calculated at the location of each particle. The origin of the SPH method can be traced back to references [1,2]. It was developed for modelling astrophysical problems and it has since been adapted for many other types of ﬂuid ﬂow, e.g. free-surface and multiphase ﬂows, and for structural modelling. Two recent general reviews of the state-of-the-art SPH are [3,4]. The ﬁrst application of SPH to free-surface ﬂows is described in [5] and the modelling of sloshing using SPH is reviewed in reference [6]. Our aim is to produce SPH simulations that reproduce an experiment ﬁrst described in reference [7] studying tuned liquid dampers, and later reproduced by ESA/ESTEC [8,9]. The experiment consists of a tank of length 0.56 m, width 0.322 m and height 0.25 m, which is subjected to sinusoidal lateral acceleration in the x-axis direction over extended periods of time. The experiment sweeps through multiple frequencies, dwelling on each for prolonged durations while gradually increasing the maximum displacement until reaching a steady state. The SPH method should in principle be perfectly suited for simulating sloshing in partially ﬁlled tanks because of its ability ∗ Corresponding author. E-mail addresses: mdgreen@imperial.ac.uk (M.D. Green), j.peiro@imperial.ac.uk (J. Peiró). to naturally capture highly non-linear behaviour and discontinuities of the free surface without needing any special treatment [10], thus greatly reducing the complexity of the simulation compared to alternative CFD methods. The simulation of violent sloshing using either weakly-compressible or incompressible formulations of SPH has been reported, for instance, in references [11] and [12], respectively, which include further references therein. However, in trying to simulate long duration sloshing with low ﬁll ratios, several diﬃculties are encountered which arise speciﬁcally from both the duration of the simulation and the stretched geometry. Consistent SPH formulations simulate advection exactly in the absence of boundaries, and hence, any error introduced into the simulation lingers over time [13]. This is a serious concern when simulations are run for long durations and millions of time-steps. Wall boundary conditions and numerical dissipation methods play a vital role in the SPH algorithm when concerning such cases, and careful consideration of each is required to achieve a successful, or even stable SPH simulations. The main challenge is to devise a SPH scheme that balances spatial and temporal accuracy with an eﬃcient computer implementation to produce accurate simulations at a reasonable computing cost. A review of the literature indicated that wall boundary treatments such as the ﬁxed boundary particle method by Marrone et al. [14] and by Adami et al. [15], and the semi-analytical boundary condition [16,17] are suitable candidates for the simulations of interest here since they address the loss of consistency of the smoothing kernel approximation in the proximity https://doi.org/10.1016/j.compﬂuid.2018.07.006 0045-7930/© 2018 Elsevier Ltd. All rights reserved. Please cite this article as: M.D. Green, J. Peiró, Long duration SPH simulations of sloshing in tanks with a low ﬁll ratio and high stretching, Computers and Fluids (2018), https://doi.org/10.1016/j.compﬂuid.2018.07.006 JID: CAF 2 ARTICLE IN PRESS [m5G;August 14, 2018;18:50] M.D. Green, J. Peiró / Computers and Fluids 000 (2018) 1–21 of the wall and limit the errors induced by the wall boundary condition. The method by Adami et al. [15] is the least complicated to implement and is also computationally cheap, which is a signiﬁcant consideration when dealing with long duration simulations. SPH schemes generally require some form of artiﬁcial dissipation to stabilise the simulation and overcome spurious effects resulting from propagation and growth of numerical errors from the initial and boundary conditions and other sources. Various dissipative terms have been proposed in the SPH literature. The most common is the artiﬁcial viscosity [18], based on the Neumann–Richtmyer artiﬁcial viscosity, which is added to the momentum conservation equation. It is effective in reducing oscillations appearing from pseudo-sound waves generated by the ﬂuid particles not being initialised in a zero-energy state [13] and helps stabilise the scheme. Due to its similarity to the physical viscosity, it has also been used to simulate it [19,20]. Spurious spatial ﬂuctuations in the density (and hence also pressure) ﬁeld might arise due to the presence of spurious modes, or chequerboarding, where a zero gradient at a particle results from contributions of non-constant values at its neighbours [21,22]. To address this issue, several artiﬁcial dissipation schemes applied to the density ﬁeld have been proposed [16,22–26]. Reference [27] reviews a variety of such schemes and shows the method of Antuono et al. [26] to be suitable for long duration simulations. Hence it is the artiﬁcial dissipation chosen in this work. It is worth mentioning that the use of a Riemann solver in the particle interactions is an alternative way of incorporating the required dissipation to dump the spurious oscillations, e.g. [28]. However, such formulation was not explored in this work. In addition to the above, the stretched nature of the ﬂuid domain requires a fairly high resolution. The reasons for this are that we require a sizeable number of particles in the vertical direction to capture the signiﬁcant changes in free-surface position and ﬂuid properties along its height, and also a regular distribution of particles to ensure SPH accuracy. This results in a very large number of ﬂuid particles because of the large stretching ratio of the computational domain under consideration here. In addition, the large number of particles along the height often leads to very small particle spacing, and hence small smoothing lengths which reduce the size of the time-step. This in turn can result in thousands or even hundreds of thousands of steps per second. The duration of the simulation together with the time-step restriction lead to considerable computational demands for SPH codes. Realistic simulations require a massively parallel code. Here the open-source code DualSPHysics [29] has been modiﬁed to include the desired features and the simulations run on GPUs. The proposed combination of the ﬁxed ghost particle method [15] and the corrected δ -SPH scheme [26], and their implementation in the eﬃcient open-source GPU code DualSPHysics is veriﬁed through hydrostatic and linear sloshing simulations. These simulations show the ability of the method to remain accurate over prolonged durations. The scheme is then validated against experimental data of lateral sloshing in a low ﬁll-ratio rectangular tank [8,9] where 2D and 3D simulations demonstrate its ability to accurately capture both wave heights and sloshing forces. 2. Numerical scheme This section describes the SPH algorithm with an emphasis on our choices of wall boundary condition treatment, numerical dissipation, and time-stepping schemes, and discuss the reasons why such choices are deem to be suitable for long duration simulation of water sloshing in low ﬁll-ratio tanks. 2.1. Governing equations The Lagrangian form of the compressible isentropic Navier– Stokes equations is Dρ = −ρ∇ · v = M Dt Dv 1 = ∇σ + fe = F Dt ρ Dx =v Dt (1) where D/Dt denotes the material derivative, ρ is the density, x is the position vector, v is the velocity, and fe is an external force per unit mass. The stress tensor σ is given by 2 σ = −pI + μ − (∇ · v ) I + ∇ v + ∇ vT 3 and the pressure, p, is calculated using a suitable equation of state. The weakly-compressible SPH formulation uses Tait’s constitutive equation p( ρ ) = p 0 ρ γ ρ0 −1 (2) ρ c2 with p0 = 0γ 0 , where ρ 0 is a reference density, c0 is the reference speed of sound and γ is the polytropic constant. For water, γ = 7 and c0 is chosen to have a nominal low Mach number, for instance Ma < 0.1. 2.2. Basic SPH The semi-discrete SPH form of the governing equations is Dρa = −ρa ∇ · va = Mab = Ma Dt (3) Dva 1 = ∇ (−pa ) + fe = Fab + fe = Fa Dt ρa (4) Dx a = va Dt (5) b b with the discrete equation of state to close the system p a = p( ρa ) = p 0 ρ γ a ρ0 −1 (6) The terms of the right-hand side of the continuity and the momentum equations can be expressed in Lagrangian conservative form, see e.g. [30], as Ma = ρa Vb (va − vb ) · ∇aWab + Da (7) b Fa = − 1 ρa Vb ( pa + pb )∇aWab + a + fe (8) b where Vb = mb /ρb is the volume associated to particle b, Wab = W (xa − xb , h ) is the SPH kernel, and ∇ a Wab is the spatial gradient of kernel evaluated at the location of particle a. For the simulations described here we have selected the Wendland kernel [31] given by W ( x − x , h ) = αd hd (2q + 1 )(1 − and its derivative ∇ W ( x − x , h ) = αd hd+1 q 4 ) ; 2 −5q 1 − 0≤q≤2 q 2 (9) 3 x − x ; 0≤q≤2 x − x Please cite this article as: M.D. Green, J. Peiró, Long duration SPH simulations of sloshing in tanks with a low ﬁll ratio and high stretching, Computers and Fluids (2018), https://doi.org/10.1016/j.compﬂuid.2018.07.006 ARTICLE IN PRESS JID: CAF [m5G;August 14, 2018;18:50] M.D. Green, J. Peiró / Computers and Fluids 000 (2018) 1–21 (10) where q = x − x , d is the number of dimensions of the problem, h is the smoothing length, and the values of the coeﬃcient α d are 47π and 1621π in two and three dimensions, respectively. The variables Da and a represent diffusion terms. The ﬁrst term is referred to as δ -SPH which is used to numerically dissipate ﬂuctuations in the density ﬁeld. The second term is a viscous term, acting as either a laminar viscous stress term for laminar ﬂows or a numerical dissipation to stabilise the scheme for inviscid ﬂows. The forms adopted for these terms are described in more detail and their use justiﬁed in Section 2.4. 1 h 2.3. Solid wall boundary treatment Solid wall boundaries should prevent particles penetrating the wall, i.e. v · n = 0 at the wall, where n is the normal to the solid wall. The wall boundary conditions should also apply a no-slip or free-slip condition depending whether the ﬂuid is modelled as viscous or inviscid respectively. However, the truncation of the kernel support by the solid boundary leads to a kernel inconsistency, namely ∗ W ( x − x , h ) d x = 1 (11) where ∗ is the truncated kernel support. To address this issue, a ﬁxed ghost particle method proposed by Adami et al. [15] is employed. In this method ghost particles are placed outside the ﬂuid domain and their properties obtained by interpolating a force balancing term between the ﬂuid and ghost particles. The following describes how to assign the required values of pressure, density and velocity of the ghost particles in order to deﬁne a continuous ﬂow ﬁeld. First, a number of layers, N, of ﬁxed ghost particles are placed on the outer side of a solid boundary to assure that ﬂuid particles always have a full kernel support. The number of layers is chosen to be N = (wh h )/x where wh is the coeﬃcient of the kernel support and x denotes the initial particle spacing. For example, when using a kernel with wh = 2 the smoothing length can be chosen to be h = 1.5x so that N = 3. To obtain a continuous ﬂow ﬁeld, the ghost particles are assigned pressure values obtained by a force balance on either side of the solid wall as pw = f p f Ww f + f ρ f (g + aw ) · (xw − x f ) Ww f f Ww f (12) where the subscripts f and w denote the ﬂuid and the ghost particles, respectively, and aw is the prescribed wall acceleration. The density of the ghost particles is obtained from the equation of state (6) as ρw = ρ0 p w p0 +1 γ1 . (13) This treatment ensures that there is no ﬂuid particle penetration of the solid wall and that the ﬂuid pressure ﬁeld adjacent to the solid walls is smooth. The ghost particles have a constant mass, taken to be equal to the initial mass of a ﬂuid particle and their velocity is calculated as follows. The velocity in the continuity Eq. (7) is set to the speciﬁed wall velocity, vs , (e.g. vs = 0 when the wall is ﬁxed) regardless of the slip condition, while in the momentum Eq. (8) the velocity is treated depending on the wall slip condition. To approximate a free-slip condition of the form v · n = 0, the interaction of a ﬂuid particle with a ghost particle ignores the viscous term and the ghost particle velocities do not have to be calculated. To impose a no-slip condition of the form v = vw , where vw is the velocity assigned to a ghost particle, the velocity ﬁeld is ﬁrst 3 extrapolated onto the ghost particle’s positions by vW v˜ a = b b ab b Wab (14) and the ghost particle’s velocity is vw = 2 vs − v˜ a . (15) 2.4. Dissipation terms Dissipation terms in SPH are introduced to either represent physical dissipation or as artiﬁcial dissipation to overcome spurious numerical effects and stabilise the numerical scheme. The physical dissipation employed here is limited to laminar viscous stress in combination with a no-slip boundary condition. The spurious pressure oscillations, pervasive in weakly compressible SPH simulations, are usually damped through the addition of artiﬁcial, or numerical, dissipation. This dissipation should have no effect on the physical model and, broadly speaking, consist of terms that are added either to the momentum Eq. (8), imitating physical viscosity, or to the continuity Eq. (7). The second group is commonly referred to as δ -SPH schemes. The following sections will describe and assess the various dissipation terms we have investigated. 2.4.1. Physical and numerical viscosity The term a in the momentum Eq. (8) represents dissipation which can either model laminar viscous stresses, or act as an artiﬁcial viscosity term to stabilise the scheme. We use the laminar viscous stress term proposed by Lo and Shao [32] given by a = 4mb b ν (xa − xb ) · ∇aWab ( v − vb ) ( ρa + ρb ) x a − x b 2 ) a (16) where ν is the kinematic viscosity of the ﬂuid. For the artiﬁcial viscosity, we use the term proposed by Monaghan in reference [18] given by ⎧ ⎨ 2m α hc0 μab ∇ W , (v − v ) · (x − x ) < 0 a a b b b (ρa + ρb ) a ab a = (17) ⎩ b 0, ( va − vb ) · ( x a − x b ) ≥ 0 where α is a problem-dependent tuning parameter and μab is deﬁned as μab = ( va − vb ) · ( x a − x b ) . xa − xb 2 (18) This term is utilised here as artiﬁcial dissipation but only used when a free-slip boundary condition is applied. 2.4.2. δ -SPH Antuono et al. [26] proposed the δ -SPH model that adds a dissipative term, Da , to the continuity Eq. (7) of the form Da = 2δ hc0 b Vb ψba ( xb − xa ) ·∇ W xb − xa 2 a ab (19) where δ is a tuning parameter. The analysis of reference [26] shows that values in the range 0 < δ ≤ 0.2 do not have a strong inﬂuence on global ﬂuid evolution. We will use the commonly adopted value of δ = 0.1 [14] in this work which has proven to be effective in reducing the spatial pressure noise. The term ψ ba can take various forms. A second-order form that approximates a Laplacian [33] is given in reference [24] reads ψba = ρb − ρa (20) and proves effective for short term dynamical simulations. In the reminder of this document we take the liberty of referring to this Please cite this article as: M.D. Green, J. Peiró, Long duration SPH simulations of sloshing in tanks with a low ﬁll ratio and high stretching, Computers and Fluids (2018), https://doi.org/10.1016/j.compﬂuid.2018.07.006 JID: CAF 4 ARTICLE IN PRESS [m5G;August 14, 2018;18:50] M.D. Green, J. Peiró / Computers and Fluids 000 (2018) 1–21 Fig. 1. The SPH particles coloured according to the pressure values after one physical second: (top) the origin is placed at bottom left corner; (bottom) the origin is placed at centre of ﬂuid domain. Fig. 2. The SPH particles coloured according to the pressure values after 100 physical seconds: (left) the DCB; (right) the GPBC. Please cite this article as: M.D. Green, J. Peiró, Long duration SPH simulations of sloshing in tanks with a low ﬁll ratio and high stretching, Computers and Fluids (2018), https://doi.org/10.1016/j.compﬂuid.2018.07.006 JID: CAF ARTICLE IN PRESS [m5G;August 14, 2018;18:50] M.D. Green, J. Peiró / Computers and Fluids 000 (2018) 1–21 5 Fig. 3. Pressure of all ﬂuid particles normalised by ρ gH plotted against the particle’s height (y) at a hundred seconds of hydrostatic simulations using different boundary conditions: (a) the DBC; (b) the GPBC. Fig. 4. Total energy normalised by Eref = ρ0 gH 3 over time (s). Showing the DBC and GPBC over the time interval 0 ≤ t ≤ 100. Fig. 5. Energy components normalised by Eref = ρ0 gH 3 over time (s). Showing the DBC and GPBC over the time interval 0 ≤ t ≤ 5: (a) kinetic, (b) potential, (c) internal. term as the Molteni δ -SPH formulation because it can be written in the form (19) although it was proposed before the introduction of the δ -SPH scheme. However, this approximation of the Laplacian is not consistent in locations close to the free surface and it has been shown that the application of such second-order dissipative terms results in a force acting outwards of the ﬂuid [27], which then causes an increase of the ﬂuid volume over long intervals of time. The results presented in Section 4 conﬁrm this behaviour. Antuono et al. [26] proposed a corrected term in the δ -SPH model to address this problem that incorporates a kernel renormalisation for the approximation of the kernel derivative. This model is formulated to guarantee that the total volume of ﬂuid is preserved and thus it is consistent with the hydrostatic pressure ﬁeld and prevents any ﬂuid volume increase. The correction is very similar to the kernel gradient correction introduced in [34] but applied here only to the calculation of the gradient of density in the dissipative term. The term reads ψba = (ρb − ρa ) − 1 ρ · (x − x ) ∇ ρb + ∇ a a b 2 (21) ρ is the corrected gradient of ρ calculated as where ∇ a a ρ = ∇ a (ρb − ρa )La ∇aWab (22) b Please cite this article as: M.D. Green, J. Peiró, Long duration SPH simulations of sloshing in tanks with a low ﬁll ratio and high stretching, Computers and Fluids (2018), https://doi.org/10.1016/j.compﬂuid.2018.07.006 JID: CAF 6 ARTICLE IN PRESS [m5G;August 14, 2018;18:50] M.D. Green, J. Peiró / Computers and Fluids 000 (2018) 1–21 Fig. 6. Pressure ﬁeld at a hundred seconds of hydrostatic simulations using the numerical dissipation schemes: (left) corrected δ -SPH (Antuono), (right) artiﬁcial viscosity (AV). Fig. 7. Pressure ﬁeld at a hundred seconds of hydrostatic simulations using the δ -SPH (Molteni) scheme. Fig. 8. Pressure of all ﬂuid particles normalised by ρ gH plotted against the particle’s height (y) at a hundred seconds of hydrostatic simulations using different numerical dissipation schemes: (a) using artiﬁcial viscosity (AV) with α = 0.01; (b) using the corrected δ -SPH scheme (Antuono) with δ = 0.1. Please cite this article as: M.D. Green, J. Peiró, Long duration SPH simulations of sloshing in tanks with a low ﬁll ratio and high stretching, Computers and Fluids (2018), https://doi.org/10.1016/j.compﬂuid.2018.07.006 ARTICLE IN PRESS JID: CAF [m5G;August 14, 2018;18:50] M.D. Green, J. Peiró / Computers and Fluids 000 (2018) 1–21 7 Fig. 9. Total energy normalised by Eref = ρ0 gH 3 over time (s). Showing the artiﬁcial viscosity (AV) and the corrected δ -SPH (Antuono) schemes over the time interval 0 ≤ t ≤ 100. Fig. 10. Energy components normalised by Eref = ρ0 gH 3 over time (s). Showing the artiﬁcial viscosity (AV) and corrected δ -SPH (Antuono) schemes over the time interval 0 ≤ t ≤ 5: (a) kinetic, (b) potential, (c) internal. with La = 2.5. Time-stepping scheme (xb − xa ) ∇aWab The time-stepping scheme used in this work is a symplectic predictor-corrector scheme [35] given by −1 (23) b 1 ρan+ 2 = ρan + n+ 12 va where the symbol denotes vector cross product. This term dissipates the spatial pressure oscillations very well without causing unphysical effects on the ﬂuid. It also has a stabilising effect that permits the removal of the artiﬁcial viscosity term from the momentum equation in inviscid ﬂows, however they can, and often are, used in tandem. The main drawback of the method is that it is signiﬁcantly more expensive than the traditional SPH dissipation. The implementation requires two additional loops through the ﬂuid particles and the solution of a 2 × 2 or 3 × 3 linear system of equations for each particle. This additional numerical expense is however a small price to pay for the ability to obtain physically realistic pressure ﬁelds. It is important to note that the correction matrix, La , is the identity matrix in the bulk of the ﬂuid when all particles are exactly evenly distributed (such as at the ﬁrst time-step of the simulation) and stays diagonally dominant when the kernel is full. However, close to the solid boundaries and the free surface, particularly in splash regions where there are very few ﬂuid particles in the neighbourhood of the ﬂuid particle of interest, the matrix is no longer diagonally dominant and could even become poorly conditioned. For this reason, we use a LU factorisation with row pivoting to improve the stability of the numerical solution of the linear system. n+ 12 xa t Mn 2 a t n = vna + F 2 a t n = xna + v 2 a ρan+1 = ρan 2 − 2+ ; = −t M n+ 12 a ρa (24) n+ 12 vna +1 = vna + tFa xna +1 = xna + t 2 vna + vna +1 where the pressure is calculated with the equation of state (6) at each half-step. To ensure stability, the time-step is calculated as t = ct min(t f , tcv ) (25) where ct is a coeﬃcient in the range 0.1 ≤ ct ≤ 0.3, and t f = min a tcv = min a h/Fa , h h ( va − vb ) · ( x a − x b ) c0 + maxb xab 2 Here tf is based on the force per unit mass and tcv combines the Courant and viscous time-step restrictions. Please cite this article as: M.D. Green, J. Peiró, Long duration SPH simulations of sloshing in tanks with a low ﬁll ratio and high stretching, Computers and Fluids (2018), https://doi.org/10.1016/j.compﬂuid.2018.07.006 ARTICLE IN PRESS JID: CAF 8 [m5G;August 14, 2018;18:50] M.D. Green, J. Peiró / Computers and Fluids 000 (2018) 1–21 Fig. 11. Linear sloshing wave elevation at the walls ζ (m) over time (s) for the SPH simulation and the analytical solution. Fig. 12. Linear sloshing ﬂuid force in x-axis direction, (N) over time (s) for the SPH simulation and the analytical solution. Fig. 13. Accumulated error of the ﬂuid force in x-axis direction in a linear sloshing simulation using different resolutions. The errors are calculated as the L2 norm of the difference between the simulated and analytical force. The experimental order of convergence (EOC) is also shown. 2.6. Energy calculation in post-processing The Hamiltonian H, a measure of the energy of a system, is given as the sum of kinetic and potential energy Given the particle nature of the SPH method, the energy of the system can be derived from the Hamiltonian of the set of particles [30]. The total energy offers a unique way of assessing the quality of a SPH scheme by measuring its energy dissipation. For instance, in a still water case an ideal numerical scheme should perfectly conserve energy. However due to various choices of formulation, boundary conditions and numerical dissipation schemes, this is rarely the case. Observing the total energy deviation from the ideal static case allows for a qualitative assessment of these choices and provides a rigorous methodology for selecting the best ones. H = Ekin + Epot (26) which for an isolated system in the absence of dissipation due to viscous stresses, is conserved. The kinetic energy for a set of Np discrete particles is given as Ekin = Np a=1 ma v2a 2 (27) where v2a = va · va . The potential energy can be further decomposed into its internal and external components denoted by Eint and Eext , where the internal energy depends on the interactions Please cite this article as: M.D. Green, J. Peiró, Long duration SPH simulations of sloshing in tanks with a low ﬁll ratio and high stretching, Computers and Fluids (2018), https://doi.org/10.1016/j.compﬂuid.2018.07.006 ARTICLE IN PRESS JID: CAF [m5G;August 14, 2018;18:50] M.D. Green, J. Peiró / Computers and Fluids 000 (2018) 1–21 9 The total energy of the system Etot is then the sum of the kig netic energy Ekin (27), the external (potential) energy Eext (28) and the internal energy Eint (32). For the purpose of post-processing, we are only interested in the conservation of energy and not its absolute value, giving the ﬁnal expression for the total energy presented in the following section as g g Etot (t ) = Ekin (t ) + Eext (t ) + Eint (t ) − (Ekin (t0 ) + Eext (t0 ) + Eint (t0 )) (33) where t0 is the initial time of the simulation. Fig. 14. Schematics of the rectangular tank used for sloshing experiments. Taken from [9]. between the discrete particles, and the external energy depends on the interactions with an external system such as gravitational pull. The external potential energy due to gravity is g Eext = Np ma g · x (28) a=1 The internal potential energy can be derived from the ﬁrst law of thermodynamics δ Eint = δ Q + δW (29) where δ Eint is the change in the internal energy of the system, δ Q is the heat transfer into the system, and δ W is the work done on the system. Energy conservation equations and thermal effects are not included in the framework, so for the purposes of this work, δ Q is dropped. The work done on the system is given by the change in the volume of the system, referred to as the pressure– volume work, and is given by δW = −pdV, where p is the pressure of the system and dV the change in volume, so that δ Eint = −pdV (30) From the continuity equation and taking the volume as V = m/ρ , the change in the internal energy can be obtained as δ Eint = mP ρ2 dρ (31) where the pressure is calculated from the equation of state (6). The internal energy is then calculated by integrating δ Eint , which after some algebra can be obtained as Eint (ρ ) = Np a=1 ma p0 ρ0 ( γ − 1 ) γ −1 ρ a ρ0 + ρ0 (γ − 1 ) − γ ρa + Eint (ρ0 ). (32) where γ is the polytropic constant from the equation of state (6). 3. GPU implementation The open-source code DualSPHysics [29] has been modiﬁed to develop the above numerical scheme for GPUs. DualSPHysics is a combined CPU and GPU code based on the SPHysics [35] solver suite. It is written in C++ with the parallel CPU implementation utilising OpenMP and the GPU implementation written in CUDA. All work shown here is restricted to the GPU version. DualSPHysics is speciﬁcally optimised to utilise the maximum eﬃciency of the GPU. However, in the modiﬁcation to the code some sacriﬁces in eﬃciency were made, speciﬁcally by the inclusion of the corrected δ -SPH scheme described in Section 2.4.2. We discuss the need for double precision and assess the eﬃciency of the modiﬁed code in the next sections. 3.1. Precision GPUs are designed mainly for fast rendering of images in computer games and optimized for performance in single precision (32 bit ﬂoating point) calculations. Their performance in double precision (64 bit ﬂoating point) operations is limited as such precision is not required for rendering. The ratio of single to double precision time per ﬂop is very low for most current GPU cards, with typical ratios of 1:24 to 1:32. Only a few dedicated general purpose GPU cards (such as the Nvidia Tesla cards and select Quadro cards) have a more reasonable ratio of 1:3 or 1:2. Therefore the use of single precision is favoured in GPU calculations. However, when a problem exhibits a wide range in scales, single precision is not suﬃcient to achieve the accuracy required [36]. To illustrate the problem a SPH simulation in single precision of static ﬂuid in a tank representing the geometry of the sloshing experiment is studied. The ratio of ﬂuid width to depth in this geometry is very high, with a value of approximately 20, which unearths the problems associated with the use of single precision to describe the locations of particles in the simulation of shallow water sloshing. In this example, the location of the origin is shifted from the centre to the corner of the tank, while all other parameters are left the same. The shift in the location of the origin lead Fig. 15. Rectangular tank normalised displacement for the frequencies f w = 0.32, f w = 0.46 Hz and f w = 0.55 Hz. Please cite this article as: M.D. Green, J. Peiró, Long duration SPH simulations of sloshing in tanks with a low ﬁll ratio and high stretching, Computers and Fluids (2018), https://doi.org/10.1016/j.compﬂuid.2018.07.006 JID: CAF 10 ARTICLE IN PRESS [m5G;August 14, 2018;18:50] M.D. Green, J. Peiró / Computers and Fluids 000 (2018) 1–21 Fig. 16. Water heights normalised by hw over time (s) in the rectangular tank sloshing experiment with excitation frequency f w1 = 0.32 Hz at location H1. Fig. 17. Water heights normalised by hw over time (s) in the rectangular tank sloshing experiment with excitation frequency f w1 = 0.32 Hz at location H2. to numerical errors in different locations of the tank. The effect of the location shift is depicted in Fig. 1. To address this issue, DualSPHysics 4.0 features an option for double precision support for the positions of particles and for updating variables. This approach avoids the scale related numerical problems while keeping the double precision operations to a minimum, and requiring only simple modiﬁcations to the code [37]. 3.2. Performance To assess the performance of the modiﬁed DualSPHysics code, we run a simple test case and compare the results of the original code using the dynamic boundary condition (DBC) with the those obtained with two versions of the modiﬁed code incorporating the ﬁxed ghost particles boundary condition (GPBC): the ﬁrst without and the second with the addition of the corrected δ -SPH scheme. The test consists of an enclosed cubic tank with a side length of 1m and ﬁlled with ﬂuid. Particle spacings x = 10−3 m in 2D and x = 0.01 m in 3D were chosen, which gave rise to 998 001 and 970 299 ﬂuid particles respectively. The DBC simulation has 40 0 0 boundary particles in 2D and 60 002 in 3D, while the GPBC case has 12 024 and 197 326. The computer used for this test has an i7 5960X octa-core CPU (16 threads) over-clocked to 3.9 GHz, 32GB of DDR4 2666 MHz RAM and a Geforce GTX Titan Z GPU. The GTX Titan Z is the very top end gaming card from the previous generation GPUs, code Please cite this article as: M.D. Green, J. Peiró, Long duration SPH simulations of sloshing in tanks with a low ﬁll ratio and high stretching, Computers and Fluids (2018), https://doi.org/10.1016/j.compﬂuid.2018.07.006 JID: CAF ARTICLE IN PRESS [m5G;August 14, 2018;18:50] M.D. Green, J. Peiró / Computers and Fluids 000 (2018) 1–21 11 Fig. 18. Sloshing force (N) over time (s) of the rectangular tank sloshing experiment with excitation frequency f w1 = 0.32 Hz. Fig. 19. Water heights normalised by hw over time (s) in the rectangular tank sloshing experiment with excitation frequency f w2 = 0.46 Hz at location H1. named Kepler, which includes 2 GPUs on one card. Each GPU has 2880 CUDA cores, 6GB of DDR5 memory and a base clock speed of 0.88 GHz. However when only one of the two GPUs is running the clock rate is boosted to 1.1 GHz, as was the case in these tests. This GPU card has a single to double precision cost of ﬂop ratio of 1:3, making it a suitable platform for double precision calculations. Each case was run twice using either single or double precision support in DualSPHysics, and their performance, evaluated in terms of steps per second, is summarised in Table 1. The results show that the GPBC performs as well as the DBC with only a minor differences. The GPBC with δ -SPH shows a decrease in performance where the number of steps per second increases by a factor of 4.5 to 7, depending on the number of dimensions and precision. This large increase in computational cost Table 1 Comparison of performance, measured in steps per second, using single and double precision for different features of the code. Test type DBC GPBC GPBC + δ -SPH 2D 2D 3D 3D 53.1 46.5 11.1 8.4 50.9 40.2 12.3 8.9 10.2 8.9 1.8 1.5 Single Double Single Double is partially due to the two additional sweeps over the ﬂuid domain, but more signiﬁcantly to the indiscriminate use of double precision storage and computations in our GPU implementation of Please cite this article as: M.D. Green, J. Peiró, Long duration SPH simulations of sloshing in tanks with a low ﬁll ratio and high stretching, Computers and Fluids (2018), https://doi.org/10.1016/j.compﬂuid.2018.07.006 JID: CAF 12 ARTICLE IN PRESS [m5G;August 14, 2018;18:50] M.D. Green, J. Peiró / Computers and Fluids 000 (2018) 1–21 Fig. 20. Water heights normalised by hw over time (s) in the rectangular tank sloshing experiment with excitation frequency f w2 = 0.46 Hz at location H2. Fig. 21. Sloshing force (N) over time (s) of the rectangular tank sloshing experiment with excitation frequency f w2 = 0.46 Hz. the corrected δ -SPH scheme. With a view towards improving the performance of the scheme, there is evidence in the literature that a more eﬃcient implementation is feasible. For instance, reference [38] reports an additional computational cost of around 40% when a similar combination of features is used in a two-dimensional CPU implementation. test is a linear sloshing simulation where we compare our SPH results against analytical solutions of potential ﬂow theory and perform a convergence test to verify the formulation using the corrected δ -SPH scheme with no artiﬁcial viscosity. 4. Veriﬁcation We use a long duration hydrostatic test to compare the differences between the DBC and the GPBC formulations. The test consists of a square two-dimensional tank, its sides parallel to the coordinate axes, of length L = 1 m (along the x-axis) ﬁlled with water to a height H = 0.1 m. Simulation parameters were set to dimen3 sional values of water, i.e. ρ0 = 10 0 0 kg/m , γ = 7, c0 = 20vmax where vmax = |gH | and g = (0, g) with g = −9.81 m/s2 . Only arti- To verify the formulation and its implementation in DualSPhysics, we carried out simulations of two test cases. The ﬁrst is a simple still water case where we compare the impact of the boundary condition against the original DualSPHsyics code which only includes the dynamic boundary condition [39]. The second 4.1. Hydrostatic boundary conditions comparison Please cite this article as: M.D. Green, J. Peiró, Long duration SPH simulations of sloshing in tanks with a low ﬁll ratio and high stretching, Computers and Fluids (2018), https://doi.org/10.1016/j.compﬂuid.2018.07.006 JID: CAF ARTICLE IN PRESS [m5G;August 14, 2018;18:50] M.D. Green, J. Peiró / Computers and Fluids 000 (2018) 1–21 ﬁcial viscosity was used, with a numerical viscosity parameter set to α = 0.01. Particle spacing was chosen to be x = 5 × 10−3 m, which gives H/x = 20 and 40 0 0 ﬂuid particles. The simulation was run for 100 s from a hydrostatic initial solution. Fig. 2 shows the pressure ﬁeld at t = 100 s for both cases. Closer inspection shows the difference in the pressure values close to the boundaries, especially close to the bottom of the tank, as well as the distortion of the free surface with the DBC which can be seen to be slightly pushed away from the wall. Fig. 3 shows the pressure value of all the particles normalised by the maximum hydrostatic pressure, p/ρ gH, against their normalised height, y/H, compared to the hydrostatic pressure proﬁle. While both boundary conditions generally agree with the hydrostatic proﬁle, especially towards the free surface (y/H = 1) the DBC shows a much wider spread, meaning that more particles deviate from the correct analytic hydrostatic pressure value. Examining the evolution of the total energy and its components will allow us to elucidate the reasons for such behaviour. Fig. 4 shows the evolution of total energy over time, normalised by a reference energy, Eref = ρ0 H 3 g, as suggested in reference [24]. In the DBC case, energy is far less accurately conserved, showing almost an order of magnitude drop in the total energy in comparison to the GPBC in the ﬁrst second of the simulation, as observed in Fig. 4. To assess the effect of the boundary condition on each of the ﬂow variables, we will analyse the components of the total energy separately. This is possible because the kinetic energy, Ekin (v), the g external potential energy due to gravity Eext (x ) and the internal energy Eint (ρ ) are functions of the velocity, position and density respectively. We will now discuss the evolution of various energy components, shown in Fig. 5, during the ﬁrst 5 s of the simulation. The main cause for the drop in total energy comes from the potential energy, as seen in Fig. 5(b). This is likely caused by the disordered distribution of particles at the free surface, which are noticeably different from their initial positions, as seen in Fig. 2. The initial oscillations are also noticeably larger using the DBC than the GPBC, indicating an overall larger reshuﬄing of the particles positions. Fig. 5(a) shows that the values of the kinetic energy calculated by the DBC scheme are far higher than those obtained by the GPBC version. This indicates that the initial velocity is noticeably higher. In both cases, the kinetic energy gradually drops due to numerical diffusion, and reaches a value close to zero as the particles settle down into a lower energy state. Finally, changes in the internal energy, which depends only on density, can be interpreted as a measure of compressibility. The internal energy oscillations, as seen in Fig. 5(c), are always opposite to the potential energy oscillations and generally cancel each other out. This can be inferred from the far smoother time history of the total energy. Both boundary conditions reach a similar steady value, but two noticeable differences are observed. The ﬁrst observation is that the initial magnitude of the internal energy is far greater using DBC, indicating that the initial pseudo-pressure waves are far stronger. The second observation is that despite the higher magnitude of the oscillations, the DBC method dampens them faster. This damping might indicate that the boundary particles absorb too much internal energy, and hence the values of density and pressure of the adjacent ﬂuid particles are underestimated, also seen in Fig. 2. To conclude, both the local properties of particles in the vicinity of the boundaries and the global energy conservation are evaluated more accurately with GPBC than with DBC. The GPBC method, while not perfect, avoids many of the issues encountered by the use of DBC and it has a comparable computational cost, as shown 13 in Section 3.2. The GPBC scheme will be used for all further studies unless speciﬁcally indicated otherwise. 4.2. Hydrostatic numerical dissipation comparison The hydrostatic tank test is used again to compare the different numerical dissipation schemes discussed in Section 2.4. These schemes are the Molteni et al. δ -SPH scheme [24], the Antuono et al. corrected δ -SPH scheme [26] and the artiﬁcial viscosity scheme [18]. These schemes will be referred to in the following as the Molteni, Antuono and AV schemes, respectively. The main purpose of this test is to assess the volume and energy conservation of the different schemes. We use the same geometry and basic ﬂuid parameters as in the previous hydrostatic case, namely L = 1 m, H = 0.1 m, x = 0.01 m, ρ0 = 10 0 0 kg/m3 , γ = 7 and c0 = 20vmax with vmax = |gH | and g = (0, g) with g = −9.81 m/s2 . Both the Molteni and Antuono schemes use a value of δ = 0.1 and no artiﬁcial viscosity, and the AV case uses a value of α = 0.01. The cases run for a duration of 100 s from a hydrostatic initialisation, allowing the simulation to develop over a long time. The value of the artiﬁcial viscosity coeﬃcient, α , should be as small as possible to provide similar values of the maximum kinetic energy and level of dissipation to the velocity which, in turn, will ensure that the motion of the particles is not overwhelmingly affected by numerical dissipation. Here we have used α = 0.01 which is slightly too dissipative for a fair comparison against the δ -SPH cases with no artiﬁcial viscosity, but lower values of α produced unstable simulations. Figs. 6 and 7 show the pressure ﬁeld for all three cases at the end of the simulation (t = 100 s). The Molteni scheme expands the particle spacing in the top part of the ﬂuid domain, raising the ﬂuid signiﬁcantly above the initial height and higher than the other two cases. In contrast, the Antuono scheme keeps the ﬂuid volume almost constant, with only a slight gap between the ﬂuid and the upper layer of particles. The AV scheme similarly preserves the ﬂuid volume and produces a more regular layer of particles at the top of the ﬂuid domain than the Antuono scheme, while the rest of the particles are slightly less ordered, especially adjacent to the walls. The increase in ﬂuid volume of Molteni δ -SPH scheme renders its use unsuitable for long duration simulations, and hence the reminder of this section will focus on assessing the energy components of the Antuono and AV schemes exclusively. Fig. 8 shows the pressure values plotted against the heights alongside the hydrostatic pressure proﬁle for the AV and Antuono schemes. The artiﬁcial viscosity scheme exhibits a noisier pressure ﬁeld than the Antuono scheme, showing it is less effective in dealing with the spurious spatial density oscillation. To better understand the effects of the numerical dissipation schemes we examine the evolution of the total energy and its components over time. The time history of the normalised total energy is shown in Fig. 9. Both schemes offer a similar ability to conserve the total energy, with the most noticeable difference being the nature of the dissipation. We investigate the behaviour of the kinetic, potential and internal components of the energy, focusing on the initial 5 s of the simulation shown in Fig. 10. The ﬁrst noticeable difference is in the kinetic energy, seen in Fig. 10(a), where the corrected δ -SPH scheme is around twice as large as the artiﬁcial viscosity. While a lower kinetic energy might seem desirable, restricting the motion of the particles artiﬁcially could also affect the dynamics of the ﬂuid if it is modelled as inviscid. On the other hand, the corrected δ -SPH scheme exhibits small oscillations throughout the simulation which are restricted by the artiﬁcial viscosity scheme. These are seen by the small but high frequency oscillations in the total energy in Fig. 9. The changes in total energy however are domi- Please cite this article as: M.D. Green, J. Peiró, Long duration SPH simulations of sloshing in tanks with a low ﬁll ratio and high stretching, Computers and Fluids (2018), https://doi.org/10.1016/j.compﬂuid.2018.07.006 ARTICLE IN PRESS JID: CAF 14 nated by the changes in potential energy which is an order of magnitude larger, as seen in Fig. 10(b). Both schemes show a similar drop in potential energy over time, which corresponds to the positions of the particles after reshuﬄing from their initial ﬁxed positions and into a lower energy state. Both schemes exhibit a similar amplitude in the initial oscillations, however these get damped much faster using the corrected δ -SPH scheme, indicating the overall shift in the particles positions are smaller. The internal energy is a measure of the compressibility of the ﬂuid and its oscillations are opposite and inverse to the potential energy. Fig. 10(c) shows that the corrected δ -SPH scheme is far more effective at damping out these oscillations, it limits the changes in the density ﬁeld, and hence improves the WCSPH approximation to an incompressible ﬂuid. The ability of the corrected δ -SPH scheme to damp the internal energy oscillations makes it effective at producing a smooth density (and hence pressure) ﬁeld and keeping the simulation stable. To conclude, the Molteni scheme causes critical failures by not conserving the ﬂuid volume and potential energy, and cannot be considered a viable numerical dissipation for long duration simulations. Both the AV and Antuono δ -SPH schemes show a similar ability to conserve the total energy. However, the artiﬁcial viscosity scheme is unable to effectively damp the internal energy, and this leaves the density (and hence pressure) ﬁeld noisier than with the corrected δ -SPH scheme. On the other hand, the corrected δ -SPH scheme induces high-frequency low-magnitude oscillations and a redistribution of particles on the free surface is less uniform, however, the bulk of the ﬂuid is left closer to the hydrostatic solution. These ﬁndings suggest that an even higher value of the numerical viscosity parameter, α , is required for the artiﬁcial viscosity scheme to produce a smoother pressure ﬁeld and damp the spatial pressure oscillations, which in turn would make the ﬂuid more viscous than intended and potentially affect the behaviour of the ﬂow. The corrected δ -SPH is capable of maintaining the hydrostatic solution with a smooth and consistent density ﬁeld while not affecting the velocity, making it suitable for both laminar and inviscid ﬂuid models, while also keeping the simulation stable. Hence, the corrected δ -SPH scheme proves to be a better choice for numerical dissipation and is used exclusively to this effect in the reminder of this work. 4.3. Linear sloshing To verify the formulation in a moving frame of reference, a linear sloshing simulation is compared against analytical values based on potential ﬂow theory. The test consists of a square two-dimensional tank with a width of length L = 1 m and ﬁlled with water to a level H = 0.1 m. The tank is subjected to a rectilinear sinusoidal excitation in the xaxis direction which is imposed in a moving frame of reference through a forcing term, namely the external force per unit mass in Eq. (8), of the form fe = (a0 (t ), −g) where g is the acceleration due to gravity and a0 (t ) = A(2π fw )2 sin(2π fw t ) is the lateral acceleration. The displacement amplitude is set to A = 0.5 m and the frequency to fw = 0.3 f0 = 0.1462 Hz where f0 is the ﬁrst natural mode obtained from the linear theory deﬁned by ωm fm = ; ωm = 2π [m5G;August 14, 2018;18:50] M.D. Green, J. Peiró / Computers and Fluids 000 (2018) 1–21 g π (m + 1 ) L tanh π (m + 1 ) L H (34) with m = 0. This large displacement amplitude is speciﬁcally chosen to ensure that the wave height is large enough to be accurately captured with SPH at the lower resolutions tested. The simulation is run for a duration of six sloshing periods, Tw = 1/ fw , which equates to 41.05 s. The excitation parameters are chosen so that a noticeable change in ﬂuid height and sloshing force occur without the need for simulations with millions of particles. However, such values are stretching the linear sloshing assumption. Following reference [40], the analytic solutions for rectilinear sinusoidal excitation in the x-axis direction of the wave elevation ζ (x, t) and sloshing force Fx (t) are given by ζ (x, t ) = ∞ βm (t ) cos mπ m=0 Fx (t ) = ρ HLa0 (t ) + ρ HL2 x L ∞ m=0 β˙ (t ) (35) 1 + (−1 )m+1 m2 (36) where β (t) and β˙ (t ) are obtained by solving the modal equation β̈ (t )m + ωm2 β (t )m = −a0 (t ) 2 π m tanh H (−1 )m − 1 . πm L (37) For purely lateral tank motion, the modal equation only has odd modes and the even modes are zero. Analytic solutions to Eq. (37) exist in some cases, but in general, solutions are found with a numerical method like the Runge–Kutta scheme. In this case, the fourth-order Runge–Kutta scheme [41] is used to solve the equation for the ﬁrst 25 modes, which should provide a suﬃciently accurate solution to Eq. (36). 3 The simulation parameters are ρ0 = 10 0 0 kg/m , γ = 7, and c0 = 10vmax where vmax = |gH | and g = −9.81 m/s2 is the magnitude of the acceleration due to gravity. Particle spacing is set to H/x = 160 generating 256 0 0 0 ﬂuid particles. A free-slip boundary condition was applied using only δ -SPH for numerical dissipation with δ = 0.1. Figs. 11 and 12 show comparisons of the left and right wave heights, and of the sloshing lateral force respectively. Both the wave heights and forces are in good agreement, but there is a noticeable increasing difference between the analytical and SPH results with time. A likely explanation is that the simulation shows a small non-linear behaviour as the assumptions of the linear sloshing are close to their limit for the value of amplitude used. For the convergence study, the simulation was run for one second with values hw /x of 20, 40, 80 and 160, resulting in 40 0 0, 16 0 0 0, 64 0 0 0 and 256 0 0 0 ﬂuid particles respectively. For a given resolution, xn , the accumulated error in the force, ε , over a simulation time, T, is evaluated as ε (xn ) = 1 T T 0 Fa (t ) − Fc (t )2 dt (38) where Fa and Fc denote the analytical and computed forces. To assess the numerical convergence of the scheme with respect to the particle resolution, x, we deﬁne an “experimental” order of convergence (EOC) as EOC = (xn ) log εε( x ) n+1 n log xxn+1 1 hw using a particle spacing xn = 20 for n = 0, 1, 2, 3. Fig. 13 plots 2n the accumulated force error against hw /x, a measure of the numerical resolution in depth, together with the values of the EOC. The method converges, but the EOC varies from case to case due to small oscillations in the SPH forces. Since the EOC is quite low, a further increase in resolution might not improve the accuracy signiﬁcantly enough to justify the additional computational cost. 5. Validation The experiment consists of a rectangular tank with a width of l = 0.59 m, depth of d = 0.335 m and height ht = 0.25 m and Please cite this article as: M.D. Green, J. Peiró, Long duration SPH simulations of sloshing in tanks with a low ﬁll ratio and high stretching, Computers and Fluids (2018), https://doi.org/10.1016/j.compﬂuid.2018.07.006 ARTICLE IN PRESS JID: CAF M.D. Green, J. Peiró / Computers and Fluids 000 (2018) 1–21 Table 2 Speciﬁcations of the rectangular tank. Speciﬁcation Value Tank height ht Tank length l Tank depth d Water height hw First natural mode f0 Stretching ratio l/hw 0.25 m 0.59 m 0.335 m 0.03 m 0.4578 Hz 19.6667 was ﬁlled with water to a level of hw = 0.03m. These dimensions, shown in Fig. 14, give the ﬁrst natural mode as f0 = 0.4578 Hz and the stretching ratio of the tank’s width to the ﬂuid height as l/hw = 19.6667. All these values are summarised in Table 2. The quantities measured for this set-up were water heights and sloshing forces. The water heights were measured using wave gauges placed at two locations, the centre of the left wall of coordinates (0, 0.5d), named H1, and the centre of the tank with coordinates (0.5l, 0.5d) which is referred to as H2. Sloshing forces were measured using tri-axial component force sensor mounted on the frame of the tank. As these sensors measured the forces of both the ﬂuid and the tank, forces from a dry run were subtracted from the raw measurements to give only the sloshing force. 5.1. Long duration 2D simulations The long duration simulations of the rectangular tank were run at three different frequencies: one below the ﬁrst natural mode, one at it, and one above it. These correspond to fw1 = 0.7 f0 = 0.32 Hz, fw2 = f0 = 0.46 Hz and fw3 = 1.202 f0 = 0.55 Hz respectively. The displacement amplitude, A, for each frequency was built up during each run to a maximum value A = 20 mm where it dwelt for a time, as illustrated in Fig. 15. Each simulation was run for 300 s using acceleration data directly from the experiment as input for the prescribed external force in the momentum conservation equation in a moving frame of reference. In the experiment, after dwelling at the maximum amplitude, the amplitude was reduced to a very small value A = 5 × 10−3 mm before the next frequency was run. While this signiﬁcantly reduced the sloshing motion between each frequency and allowed reaching a steady state, the initial conditions did not correspond to rest. This was a signiﬁcant difference with respect to the initial conditions used in the simulations that started from complete rest. This has limited our ability to study the transient phase of the sloshing but, on the other hand, the long duration of these tests has allowed the synchronisation of experiment and simulation before reaching a steady state. Due to the long duration of the experiment at each frequency, the simulations were attempted in two-dimensional space, using a particle spacing x = 5 × 10−4 m. This resolution represents 60 particles in the ﬂuid height, i.e. hw /x = 60, and gives rise to 70 800 ﬂuid particles and 6576 boundary particles. The ﬂuid is wa3 ter, ρ0 = 10 0 0 kg/m , γ = 7 and the numerical speed of sound was chosen as c0 = 20 ghw . The simulation is performed using the corrected δ -SPH scheme with a value of δ = 0.1, a laminar dynamic viscosity ν = 10−6 kg/(ms), and a no-slip boundary condition at the walls of the tank. At this spatial resolution, it is not possible to capture velocity gradients within the boundary layer, but the use of laminar viscosity assists in controlling the behaviour of thin ﬂuid layers travelling up the walls and in reducing the number of stray particles that slide up the walls in the splash region and leave the domain. Despite its apparent simplicity, the rectangular tank test proved to be extremely challenging due to two factors: the low ﬁll ratio [m5G;August 14, 2018;18:50] 15 of the tank and the extended durations of the experiment at each frequency. Regarding the ﬁll ratio, the high stretching ratio of the liquid mass used in the experiment required a small particle spacing x = 5 × 10−4 m to represent the ﬂuid height with suﬃcient resolution. This caused a large variation in scale between the particle spacing and the width of the tank which required the use of double precision to reduce numerical errors, and in turn a signiﬁcant increase of the computational cost. A simulation of 300 s with a smoothing length h = 1.5x, which corresponds to 7.5 × 10−4 m, required approximately 50 million time-steps for each run and took around one week on a Nvidia Tesla K40 GPU card. This number of time-steps is a direct consequence of the time-step restriction of explicit time-stepping schemes. Since SPH simulates advection exactly, any errors introduced during the simulation, e.g. via the boundary conditions, will persist and can accumulate over time. Given the huge number of time-steps required here, this seemly advantageous property of SPH also makes these simulations very sensitive to spurious effects. Despite these diﬃculties, the proposed SPH scheme was capable of reasonably simulating the experiment, as described in detail in the following paragraphs. The ﬁrst case corresponds to a frequency fw1 = 0.32 Hz which is lower than that of the ﬁrst natural mode. Comparison between the experiment and simulation showing the normalised water height at the wall and centre of the tank are given in Fig. 16 and Fig. 17 respectively, and the sloshing force in Fig. 18. This case generated very low ﬂuid motion with no wavefront and was found to be the hardest to simulate as it suffers from many of the known issues found in the hydrostatic case, making it very sensitive to numerical errors. Further, the wave elevation at the centre of the tank is so small that using only 60 particles to represent the height of the ﬂuid means that the difference between the initial and maximum water level is less than the particle spacing. This is evident especially from the comparison during the transient phase. Despite this, the simulation still displays good agreement with the experiment towards the end of the 300 s period for both the water heights and force. However, it should be noted that the force is far smoother than that measured in the experiment, but with the same frequency and similar magnitude. Further, up to approximately 250 s, the wave heights at the centre of the tank (H2) were too low in the simulated case to be properly captured at this resolution. The frequency of the second case is fw2 = 0.46 Hz which corresponds to the ﬁrst natural mode. Comparison between the experiment and simulation showing the normalised water height at the wall and centre of the tank are given in Fig. 19 and Fig. 20 respectively, and the sloshing force in Fig. 21. This case exhibits more complex dynamics: a proper wavefront formed and large splashes were observed when the wave impacted the wall, although there is no clear evidence of wave breaking. This case proved to be the simplest to simulate, with the simulation catching up with to the experiment quite early in the transient stage. After 200 s, the secondary waves are also reasonably captured at both locations. The force shows very similar peaks, with only some discrepancies in the shape, but is in very good overall agreement. The simulated water heights at location H1 show spikes which are far higher than those measured in the experiment. This issue requires further investigation, but it is conjectured that such discrepancy is partially due to discretisation errors in the SPH summations, associated with the ﬂuid splashes where there are only a very few particles involved in the interaction. The third case with a frequency fw3 = 0.55 Hz, which is above the ﬁrst natural mode. Comparison between the experiment and simulation showing the water height at the wall and the centre of the tank are given in Fig. 22 and Fig. 23 respectively, and the slosh- Please cite this article as: M.D. Green, J. Peiró, Long duration SPH simulations of sloshing in tanks with a low ﬁll ratio and high stretching, Computers and Fluids (2018), https://doi.org/10.1016/j.compﬂuid.2018.07.006 JID: CAF 16 ARTICLE IN PRESS [m5G;August 14, 2018;18:50] M.D. Green, J. Peiró / Computers and Fluids 000 (2018) 1–21 Fig. 22. Water heights normalised by hw over time (s) in the rectangular tank sloshing experiment with excitation frequency f w3 = 0.55 Hz at location H1. Fig. 23. Water heights normalised by hw over time (s) in the rectangular tank sloshing experiment with excitation frequency f w3 = 0.55 Hz at location H2. ing force in Fig. 24. This case displays all of the characteristics of the previous case, but with the addition of wave breakage. The results are once again in good agreement for both water heights and force. The biggest difference is that the water heights take much longer to catch up with the experiment. More speciﬁcally, the secondary waves following the wave-front only get resolved towards the end of the 300 s period. Similarly to the previous case, large spikes are observed at location H1 when the ﬂuid splashes. Snapshots representing a full sloshing cycle after reaching the steady state are provided for the cases with frequencies f w2 = 0.46 Hz and fw3 = 0.55 Hz in Figs. 25 and 26, respectively. These images show the free surface from the SPH simulations at similar stages of the sloshing cycle overlaid in red. The snapshots are taken from a video recording of the experiment and the exact times cor- responding to the start of the simulation at the chosen frequencies are unavailable. These snapshots serve to illustrate that the overall shape and features of the ﬂow between the simulations and the experiment show similarities, even after the very long time they take to reach the steady state. The fw1 = 0.32 Hz case is not shown here as the dynamics are so limited that there is no added value in a visual observation of the snapshots. 5.2. Short duration steady-state simulations in 3D The simulation reached a steady state after only a few sloshing cycles, allowing for comparison against the steady-state phase of the ﬂow in the experiment. The period of 10–30 s in the simulation is compared against the period of 280–300 s in the ex- Please cite this article as: M.D. Green, J. Peiró, Long duration SPH simulations of sloshing in tanks with a low ﬁll ratio and high stretching, Computers and Fluids (2018), https://doi.org/10.1016/j.compﬂuid.2018.07.006 JID: CAF ARTICLE IN PRESS [m5G;August 14, 2018;18:50] M.D. Green, J. Peiró / Computers and Fluids 000 (2018) 1–21 17 Fig. 24. Sloshing force (N) over time (s) of the rectangular tank sloshing experiment with excitation frequency f w3 = 0.55 Hz. Fig. 25. Snapshots of the rectangular tank experiment for a whole cycle after reaching a steady state under an excitation frequency of f w2 = 0.46 Hz with the SPH free-surface overlaid in red. (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version of this article.) Fig. 26. Snapshots of the rectangular tank experiment for a whole cycle after reaching a steady state under an excitation frequency of f w3 = 0.55 Hz with the SPH free-surface overlaid in red. (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version of this article.) Please cite this article as: M.D. Green, J. Peiró, Long duration SPH simulations of sloshing in tanks with a low ﬁll ratio and high stretching, Computers and Fluids (2018), https://doi.org/10.1016/j.compﬂuid.2018.07.006 JID: CAF 18 ARTICLE IN PRESS [m5G;August 14, 2018;18:50] M.D. Green, J. Peiró / Computers and Fluids 000 (2018) 1–21 Fig. 27. Water heights normalised by hw over time (s) in the rectangular tank sloshing experiment with an excitation frequency f w2 = 0.46 Hz. Showing the time interval 280 ≤ t ≤ 300 s in the experiment and 10 ≤ t ≤ 30 s in the three-dimensional SPH simulation when both have reached a steady state: (a) water heights at H1; (b) water heights at H2. Fig. 28. Sloshing forces (N) over time (s) in the rectangular tank sloshing experiment with an excitation frequency f w2 = 0.46 Hz. Showing the time interval 280 ≤ t ≤ 300 s in the experiment and 10 ≤ t ≤ 30 s in the three-dimensional SPH simulation when both have reached a steady state. periment, where both are at the steady-state phase. The simulation is run with a fairly low resolution, x = 1.5 × 10−3 m, which due to the stretching ratio resulted in 1 752 780 ﬂuid particles and 904 821 boundary particles. This modest resolution with 20 particles in the ﬂuid height, i.e. hw /x = 20, is as low as it can reasonably be used without degrading the quality of the simulation too much to render it unusable. The ﬂuid parameters are the same as in the two-dimensional case, but using a free-slip condition instead of the no-slip. Only the δ -SPH scheme with δ = 0.1 was used and no artiﬁcial viscosity was added. The free-slip condition can be used here because the no-slip condition has negligible effect on the ﬂuid ﬂow at this resolution. Even though the particle spacing, and hence the smoothing length, is much larger than in the two-dimensional case, the simulation still required approximately 3 million time-steps. Considering the number of particles and the cost of using the corrected δ -SPH scheme in three dimensions, this is still extremely expensive and required approximately 5 weeks to simulate on a Nvidia Tesla K20 GPU card. The water heights in both locations and the sloshing force all show a strong agreement, with most of the secondary waves also captured in a very good manner as seen in Figs. 27 and 28. It is worth mentioning that the splashes in the 3D case are slightly smaller than in the 2D case. This reinforces the assumption for the cause of this discrepancy, as the ﬂuid impacting the wall in 3D has more neighbours to interact with thus reducing the effect of the discretisation error (Fig. 30). To better understand the importance of using the δ -SPH scheme in these type of simulations, an additional simulation was run without δ -SPH and using the artiﬁcial viscosity term (17) with a value of α = 0.02 for numerical stability only. The comparison of the two simulations is shown in Figs. 29 and 31. While the water heights at both locations are very similar and agree reasonably Please cite this article as: M.D. Green, J. Peiró, Long duration SPH simulations of sloshing in tanks with a low ﬁll ratio and high stretching, Computers and Fluids (2018), https://doi.org/10.1016/j.compﬂuid.2018.07.006 JID: CAF ARTICLE IN PRESS M.D. Green, J. Peiró / Computers and Fluids 000 (2018) 1–21 [m5G;August 14, 2018;18:50] 19 Fig. 29. Water heights normalised by hw over time (s) in the rectangular tank sloshing experiment with an excitation frequency f w2 = 0.46 Hz. Showing the time interval 280 ≤ t ≤ 300 s in the experiment and 10 ≤ t ≤ 30 s in the three-dimensional SPH simulations when both have reached a steady state. The ﬁrst simulating is using artiﬁcial viscosity only with α = 0.02 (AV) and the second simulation is using δ -SPH only with δ = 0.1 (δ -SPH): (a) water heights at H1; (b) water heights at H2. Fig. 30. Sloshing forces (N) over time (s) in the rectangular tank sloshing experiment with an excitation frequency f w2 = 0.46 Hz. Showing the time interval 280 ≤ t ≤ 300 s in the experiment and 10 ≤ t ≤ 30 s in the three-dimensional SPH simulation when both have reached a steady state. The ﬁrst simulating is using artiﬁcial viscosity only with α = 0.02 (AV) and the second simulation is using δ -SPH only with δ = 0.1 (δ -SPH). Fig. 31. The SPH particles coloured by density from the 3D simulation of the rectangular tank during the impact with wall at a t = 12.5s: (a) using artiﬁcial viscosity only with α = 0.02: (b) using δ -SPH only with δ = 0.1. Please cite this article as: M.D. Green, J. Peiró, Long duration SPH simulations of sloshing in tanks with a low ﬁll ratio and high stretching, Computers and Fluids (2018), https://doi.org/10.1016/j.compﬂuid.2018.07.006 JID: CAF 20 ARTICLE IN PRESS [m5G;August 14, 2018;18:50] M.D. Green, J. Peiró / Computers and Fluids 000 (2018) 1–21 well with the experiment, notable differences can be observed in the density ﬁeld and the sloshing force. The simulation without δ -SPH shows a much nosier density ﬁeld, as expected, and some variation in the free surface shape, but it is really the overestimation of the forces which truly demonstrates the importance of having proper numerical dissipation. These results clearly validate the ability of the SPH scheme to accurately simulate sloshing under various conditions in both 2D and 3D. They demonstrate the stability of the simulations over very long durations and after millions of time steps, and illustrate the importance of selecting the right formulation for numerical dissipation, justifying the large computational cost of the corrected δ -SPH scheme. It is only with the careful selection of the boundary condition combined with the dissipation methods, and their implementation on GPUs, that has allowed the SPH method to accurately simulate these experiments. 6. Conclusions Our attempt to perform long duration SPH sloshing simulations towards reaching a steady state in a rectangular tank with a low ﬁll ratio, and to validate them against experimental data has highlighted three areas of critical importance for the SPH method to produce stable and accurate results for this type of sloshing problems. These are the appropriate treatment of solid wall boundary conditions, the choice of a suitable numerical dissipation term and the issues associated with the computational implementation of the SPH algorithms on GPUs, namely eﬃciency and sensitivity to round-off error. The treatment of solid wall boundary conditions using the ﬁxed ghost particle method of Adami et al. [15] includes a force balance correction which greatly improves the smoothness of the ﬂuid properties adjacent to the solid walls. It is also simple to implement and numerically eﬃcient. Comparisons with other commonly used SPH wall boundary conditions show that this method significantly reduces the pressure noise close to the solid walls, it does not generate a non-physical gap between the wall and the ﬂuid, and it is better at conserving the total energy. Numerical dissipation is required in weakly compressible formulations of SPH to stabilize the simulations and to reduce the commonly found spurious spatial oscillations of pressure. The corrected δ -SPH scheme proposed in Antuono et al. [26] was identiﬁed as the most suitable form of numerical dissipation for long duration simulations. The hydrostatic and sloshing veriﬁcation tests demonstrated its ability to quickly damp oscillations in potential energy whilst conserving volume and to produce correct sloshing forces. The cost of SPH simulations for realistic problems requires efﬁcient parallelisation. The proposed wall boundary conditions and dissipation terms were implemented in the open-source SPH solver DualSPHysics (version 4) [29] using GPU architectures. Simulations in seemingly simple cases such as low ﬁll-ratio sloshing in tank showed that the use of single precision, as commonly favoured in GPU implementations, resulted in signiﬁcant numerical inaccuracies that rendered long duration simulations invalid. This led us to identify the need for double precision to represent the positions of the particles. The selective use of double precision, together with the additional cost of the corrected δ -SPH scheme, represents a signiﬁcant increase in computational cost. The indiscriminate use of double precision in the solution of the linear system in the corrected δ -SPH scheme is mostly responsible for the ﬁve to seven-fold increase is computational time, depending on the speciﬁc problem and simulation parameters. There is however ample scope for improving the performance of the code by only using double precision when required and employing smarter memory management on the GPU. Further, we should investigate more sta- ble and robust iterative techniques for solving the linear system in single precision which could perform faster in current GPU architectures. Finally, the latest NVIDIA GPU architectures offer up to ﬁve times more double-precision ﬂoating point operations per second, and over three times the memory bandwidth. We estimate that the code could potentially beneﬁt from a speed-up of an order of magnitude or more if the previous enhancements are incorporated. Acknowledgements The help and advice of Drs Stefan Adami and Xiangyu Hu, both of Technische Universität München, with the implementation of the SPH formulation that they proposed is greatly acknowledged. We would also like to thank Jan Eichstädt for developing, as part of his MSc project, the Matlab SPH code of Adami’s formulation that has allowed us to identify and solve the numerical precision issues we experienced with DualSPHysics (version 3) and for collaboration with the development of the post-processing tools. We are very grateful to Drs José Manuel Domínguez Alonso and Alejandro Jacobo Cabrera Crespo, both of Universidad de Vigo, for providing us with a pre-release double precision version of DualSPHysics (version 4). This work was partly sponsored by ESA/ESTEC under statement of work TEC-MPA/2012/In. We thank Dr José Longo, the project monitor, Andrea Passaro and Chiara Lada for their support and advice on sloshing phenomena, for making the experimental data available to us and granting permission for its publication, and for helping us with the interpretation of the experimental and computational data. References [1] Gingold R, Monaghan J. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Mon Not R Astron Soc 1977;181(3):375–89. [2] Lucy L. A numerical approach to testing of the ﬁssion hypothesis. Astron J 1977;82:1013–24. [3] Liu M, Liu G. Smoothed particle hydrodynamics (SPH): an overview and recent developments. Arch Comput Methods Eng 2010;17(1):25–76. [4] Monaghan J. Smoothed particle hydrodynamics and its diverse applications. Annu Rev Fluid Dyn 2012;44:323–46. [5] Monaghan J. Simulating free surface ﬂows with SPH. J Comput Phys 1994;110(2):399–406. [6] Gómez-Gesteira M, Rogers B, Dalrymple R, Crespo J. State-of-the-art of classical SPH for free-surface ﬂows. J Hydraul Res 2010;48(S1):6–27. [7] Reed D, Yu J, Yeh H, Gardarsson S. Investigation of tuned liquid dampers under large amplitude excitation. J Eng Mech 1998;124(4):405–13. [8] Lada C, Such-Taboada M, Ngan I, Grigore L, Appolloni M, Roure S, et al. Experimental sloshing reference test 13th European conf on spacecraft structures, materials & environmental testing, ESTEC, Noordwijk, The Netherlands, ESA SP-727, June 2014. Ouwehand L, editor. Braunschweig, Germany: ESA Communications; 2014. [9] Lada C, Passaro A. Experimental sloshing references test report. Tech. Rep.. ESA-ESTEC TEC-MPA Report; 2016. [10] Colagrossi A, Antuono M, Touzé DL. Theoretical considerations on the free surface role in the smoothed-particle-hydrodynamics model. Phys Rev E 2009;79(5):056701. [11] Bouscasse B, Antuono M, Colagrossi A, Lugni C. Numerical and experimental investigation of nonlinear shallow water sloshing. Int J Nonlinear Sci Numer Simul 2013;14(2):123–38. [12] Gotoh H, Khayyer A, Ikari H, Arikawa T, Shimosako K. On enhancement of incompressible SPH method for simulation of violent sloshing ﬂows. Appl Ocean Res 2014. [13] Price D. Smoothed particle hydrodynamics: things I wish my mother taught me. Advances in computational astrophysics: methods, tools and outcomes, Cefalu, Italy; 2011. arxiv: 1111.1259. [14] Marrone S, Antuono M, Colagrossi A, Colicchio G, Touzé DL, Graziani G. δ -SPH model for simulating violent impact ﬂows. Comput Methods Appl Mech Eng 2011;200(13):1526–42. [15] Adami S, Hu X, Adams N. A generalized wall boundary condition for smoothed particle hydrodynamics. J Comput Phys 2012;231(21):7057–75. [16] Mayrhofer A, Rogers B, Violeau D. Investigation of wall bounded ﬂows using SPH and the uniﬁed semi-analytical wall boundary conditions. Comput Phys Commun 2013;184(11):2515–27. [17] Ferrand M, Laurence D, Rogers B, Violeau D, Kassiotis C. Uniﬁed semi-analytical wall boundary conditions for inviscid, laminar or turbulent ﬂows in the meshless SPH method. Int J Numer Methods Fluids 2012;71(4):446–72. Please cite this article as: M.D. Green, J. Peiró, Long duration SPH simulations of sloshing in tanks with a low ﬁll ratio and high stretching, Computers and Fluids (2018), https://doi.org/10.1016/j.compﬂuid.2018.07.006 JID: CAF ARTICLE IN PRESS M.D. Green, J. Peiró / Computers and Fluids 000 (2018) 1–21 [18] Monaghan J. Smoothed particle hydrodynamics. Annu Rev Astron Astrophys 1992;30:543–74. [19] Espanol P., Revenga M.. Smoothed dissipative particle dynamics. Phys Rev E 67 (2) 026705. [20] Hu X, Adams N. Angular-momentum conservative smoothed particle dynamics for incompressible viscous ﬂows. Phys Fluids (1994–present) 2006;18(10):101702. [21] Bonet J, Kulasegaram S. A simpliﬁed approach to enhance the performance of smooth particle hydrodynamics methods. Appl Math Comput 2002;126(2–3):133–55. [22] Fatehi R, Manzari M. A remedy for numerical oscillations in weakly compressible smoothed particle hydrodynamics. Int J Numer Methods Fluids 2011;67(9):1100–14. [23] Colagrossi A, Landrini M. Numerical simulation of interfacial ﬂows by smoothed particle hydrodynamics. J Comput Phys 2003;191(2):448–75. [24] Molteni D, Colagrossi A. A simple procedure to improve the pressure evaluation in hydrodynamic context using the SPH. Comput Phys Commun 2009;180(6):861–72. [25] Ferrari A, Dumbser M, Toro E, Armanini A. A new 3D parallel SPH scheme for free surface ﬂows. Comput Fluids 2009;38(6):1203–17. [26] Antuono M, Colagrossi A, Marrone S, Molteni D. Free-surface ﬂows solved by means of SPH schemes with numerical diffusive terms. Comput Phys Commun 2010;181(3):532–49. [27] Antuono M, Colagrossi A, Marrone S. Numerical diffusive terms in weakly– compressible SPH schemes. Comput Phys Commun 2012;183(12):2570–80. [28] Koukouvinis PK, Anagnostopoulos JS, Papantonis DE. An improved MUSCL treatment for the SPH–ALE method: comparison with the standard SPH method for the jet impingement case. Int J Numer Methods Fluids 2013;71:1152–77. [29] Crespo A, Domínguez J, Rogers B, Gómez-Gesteira M, Longshaw S, Canelas R, et al. DualSPHysics: open-source parallel CFD solver based on smoothed particle hydrodynamics (SPH). Comput Phys Commun 2015;187:204–16. DualSPHysics v4 is available at: www.dual.sphysics.org. [m5G;August 14, 2018;18:50] 21 [30] Violeau D. Fluid mechanics and the SPH method: theory and applications. Oxford University Press; 2012. [31] Wendland H. Piecewise polynomial, positive deﬁnite and compactly supported radial functions of minimal degree. Adv Comput Math 1995;4(1):389–96. [32] Lo E, Shao S. Simulation of near-shore solitary wave mechanics by an incompressible SPH method. Appl Ocean Res 2002;24(5):275–86. [33] Violeau D. Dissipative forces for Lagrangian models in computational ﬂuid dynamics and application to smoothed-particle hydrodynamics. Phys Rev E 2009;80(3):036705. [34] Bonet J, Lok T-S. Variational and momentum preservation aspects of smooth particle hydrodynamic formulations. Comput Methods Appl Mech Eng 1999;181(1):97–116. [35] Gómez-Gesteira M, Rogers B, Crespo A, Dalrymple R, Narayanaswamy M, Domínguez J. SPHysics - development of a free-surface ﬂuid solver. Part 1: theory and formulations. Comput Geosci 2012;48:289–99. SPHysics is available at: www.sphysics.org. [36] Domínguez J, Crespo A, Gómez-Gesteira M, Rogers B. Simulating more than 1 billion SPH particles using GPU hardware acceleration. In: Proceedings of the 8th international SPHERIC workshop, Trondheim, Norway; 2013. [37] Domínguez J. DualSPHysics: towards high performance computing using SPH technique. Environmental Physics Laboratory, University of Vigo; 2014. Ph.D. thesis. [38] Valizadeh A, Monaghan JJ. A study of solid wall models for weakly compressible SPH. J Comput Phys 2015;300:5–19. [39] Crespo A, Gómez-Gesteira M, Dalrymple R. Boundary conditions generated by dynamic particles in SPH methods. Comput Mater Continua 2007;5(3):173–84. [40] Faltinsen M, Timokha A. Sloshing. Cambridge University Press; 2009. [41] Press W, Flannery B, Teukolsky S, Vetterling W. Runge-Kutta method. In: Numerical recipes in FORTRAN: the art of scientiﬁc computing. Cambridge University Press; 1992. p. 704–16. Please cite this article as: M.D. Green, J. Peiró, Long duration SPH simulations of sloshing in tanks with a low ﬁll ratio and high stretching, Computers and Fluids (2018), https://doi.org/10.1016/j.compﬂuid.2018.07.006

1/--страниц