close

Вход

Забыли?

вход по аккаунту

?

j.compfluid.2018.07.006

код для вставкиСкачать
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 fill
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 fill 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 difficulties and common pitfalls of such simulations
and propose a SPH method to deal with them effectively. The formulation combines an efficient 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 fluid flow 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 fluid flow, e.g. free-surface and multiphase
flows, and for structural modelling. Two recent general reviews of
the state-of-the-art SPH are [3,4]. The first application of SPH to
free-surface flows 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 first 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 filled 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 fill
ratios, several difficulties are encountered which arise specifically 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
efficient computer implementation to produce accurate simulations
at a reasonable computing cost.
A review of the literature indicated that wall boundary treatments such as the fixed 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.compfluid.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 fill ratio and high stretching,
Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.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 significant consideration when dealing with long duration
simulations.
SPH schemes generally require some form of artificial 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 artificial viscosity [18], based on
the Neumann–Richtmyer artificial viscosity, which is added to the
momentum conservation equation. It is effective in reducing oscillations appearing from pseudo-sound waves generated by the fluid
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 fluctuations in the density (and hence also
pressure) field 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 artificial dissipation schemes applied to
the density field 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 artificial 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 fluid 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 significant changes in free-surface position and fluid
properties along its height, and also a regular distribution of particles to ensure SPH accuracy. This results in a very large number
of fluid 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 modified to include the desired features
and the simulations run on GPUs.
The proposed combination of the fixed ghost particle method
[15] and the corrected δ -SPH scheme [26], and their implementation in the efficient open-source GPU code DualSPHysics is verified 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 fill-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 fill-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 fill ratio and high stretching,
Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.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 coefficient
α d are 47π and 1621π in two and three dimensions, respectively. The
variables Da and a represent diffusion terms. The first term is
referred to as δ -SPH which is used to numerically dissipate fluctuations in the density field. The second term is a viscous term,
acting as either a laminar viscous stress term for laminar flows or
a numerical dissipation to stabilise the scheme for inviscid flows.
The forms adopted for these terms are described in more detail
and their use justified 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 fluid 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
fixed ghost particle method proposed by Adami et al. [15] is employed. In this method ghost particles are placed outside the fluid
domain and their properties obtained by interpolating a force balancing term between the fluid and ghost particles. The following
describes how to assign the required values of pressure, density
and velocity of the ghost particles in order to define a continuous
flow field.
First, a number of layers, N, of fixed ghost particles are placed
on the outer side of a solid boundary to assure that fluid particles
always have a full kernel support. The number of layers is chosen to
be N = (wh h )/x where wh is the coefficient 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 flow field, 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 fluid 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 fluid particle penetration of
the solid wall and that the fluid pressure field adjacent to the solid
walls is smooth. The ghost particles have a constant mass, taken to
be equal to the initial mass of a fluid particle and their velocity is
calculated as follows. The velocity in the continuity Eq. (7) is set to
the specified wall velocity, vs , (e.g. vs = 0 when the wall is fixed)
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 fluid 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 field is first
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 artificial 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 artificial, 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 artificial 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 fluid.
For the artificial 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 defined as
μab =
( va − vb ) · ( x a − x b )
.
xa − xb 2
(18)
This term is utilised here as artificial 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
influence on global fluid 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 fill ratio and high stretching,
Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.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 fluid 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 fill ratio and high stretching,
Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.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 fluid 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 fluid [27], which then causes an increase of the fluid volume over long intervals of time. The results
presented in Section 4 confirm 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 fluid
is preserved and thus it is consistent with the hydrostatic pressure
field and prevents any fluid 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 fill ratio and high stretching,
Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.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 field at a hundred seconds of hydrostatic simulations using the numerical dissipation schemes: (left) corrected δ -SPH (Antuono), (right) artificial viscosity
(AV).
Fig. 7. Pressure field at a hundred seconds of hydrostatic simulations using the δ -SPH (Molteni) scheme.
Fig. 8. Pressure of all fluid 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 artificial 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 fill ratio and high stretching,
Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.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 artificial 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 artificial 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 fluid. It also has a stabilising effect that permits the removal of the artificial viscosity
term from the momentum equation in inviscid flows, however they
can, and often are, used in tandem. The main drawback of the
method is that it is significantly more expensive than the traditional SPH dissipation. The implementation requires two additional
loops through the fluid 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 fields. It is important to note that
the correction matrix, La , is the identity matrix in the bulk of the
fluid when all particles are exactly evenly distributed (such as at
the first 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 fluid particles in the neighbourhood of the fluid 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 coefficient 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 fill ratio and high stretching,
Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.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 fluid force in x-axis direction, (N) over time (s) for the SPH simulation and the analytical solution.
Fig. 13. Accumulated error of the fluid 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 fill ratio and high stretching,
Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.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 final 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 first 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 modified 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 specifically optimised to utilise the maximum efficiency of the
GPU. However, in the modification to the code some sacrifices in
efficiency were made, specifically by the inclusion of the corrected
δ -SPH scheme described in Section 2.4.2. We discuss the need for
double precision and assess the efficiency of the modified 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 floating point) calculations. Their performance in double precision (64 bit floating point) operations is limited as such precision
is not required for rendering. The ratio of single to double precision
time per flop 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 sufficient to achieve the accuracy required [36].
To illustrate the problem a SPH simulation in single precision of
static fluid in a tank representing the geometry of the sloshing experiment is studied. The ratio of fluid 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 fill ratio and high stretching,
Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.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 modifications to the code [37].
3.2. Performance
To assess the performance of the modified 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 modified code incorporating the
fixed ghost particles boundary condition (GPBC): the first 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 filled with fluid. 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 fluid 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 fill ratio and high stretching,
Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.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 flop 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 fluid domain, but more significantly 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 fill ratio and high stretching,
Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.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 efficient 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 flow theory and perform a convergence test to verify the formulation using the corrected δ -SPH scheme with no artificial viscosity.
4. Verification
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) filled 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 first
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 fill ratio and high stretching,
Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.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
ficial 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 fluid particles. The simulation
was run for 100 s from a hydrostatic initial solution.
Fig. 2 shows the pressure field 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 profile. While both boundary
conditions generally agree with the hydrostatic profile, 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 first second of the simulation, as observed in Fig. 4.
To assess the effect of the boundary condition on each of the
flow 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 first 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 reshuffling 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 first
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 fluid 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 specifically 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 artificial 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 fluid 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 artificial 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 artificial viscosity coefficient, α , 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 artificial viscosity, but lower values of α produced
unstable simulations.
Figs. 6 and 7 show the pressure field 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 fluid domain, raising
the fluid significantly above the initial height and higher than the
other two cases. In contrast, the Antuono scheme keeps the fluid
volume almost constant, with only a slight gap between the fluid
and the upper layer of particles. The AV scheme similarly preserves
the fluid volume and produces a more regular layer of particles
at the top of the fluid domain than the Antuono scheme, while
the rest of the particles are slightly less ordered, especially adjacent to the walls. The increase in fluid 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 profile for the AV and Antuono
schemes. The artificial viscosity scheme exhibits a noisier pressure
field 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 first noticeable difference is in
the kinetic energy, seen in Fig. 10(a), where the corrected δ -SPH
scheme is around twice as large as the artificial viscosity. While a
lower kinetic energy might seem desirable, restricting the motion
of the particles artificially could also affect the dynamics of the
fluid 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 artificial 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 fill ratio and high stretching,
Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.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 reshuffling from their initial fixed 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 fluid 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 field, and
hence improves the WCSPH approximation to an incompressible
fluid. The ability of the corrected δ -SPH scheme to damp the internal energy oscillations makes it effective at producing a smooth
density (and hence pressure) field and keeping the simulation stable.
To conclude, the Molteni scheme causes critical failures by not
conserving the fluid 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 artificial viscosity
scheme is unable to effectively damp the internal energy, and this
leaves the density (and hence pressure) field 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 fluid is left closer to the hydrostatic solution.
These findings suggest that an even higher value of the numerical viscosity parameter, α , is required for the artificial viscosity
scheme to produce a smoother pressure field and damp the spatial
pressure oscillations, which in turn would make the fluid more viscous than intended and potentially affect the behaviour of the flow.
The corrected δ -SPH is capable of maintaining the hydrostatic solution with a smooth and consistent density field while not affecting
the velocity, making it suitable for both laminar and inviscid fluid
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 flow theory.
The test consists of a square two-dimensional tank with a width
of length L = 1 m and filled 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 first natural
mode obtained from the linear theory defined 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 specifically 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 fluid 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 first 25 modes, which should provide a sufficiently 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 fluid 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 fluid 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 define 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 significantly 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 fill ratio and high stretching,
Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.07.006
ARTICLE IN PRESS
JID: CAF
M.D. Green, J. Peiró / Computers and Fluids 000 (2018) 1–21
Table 2
Specifications of the rectangular tank.
Specification
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 filled with water to a level of hw = 0.03m. These dimensions,
shown in Fig. 14, give the first natural mode as f0 = 0.4578 Hz
and the stretching ratio of the tank’s width to the fluid 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 fluid 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 first 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 significantly reduced the sloshing motion between each frequency and allowed reaching a steady state,
the initial conditions did not correspond to rest. This was a significant 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 fluid height, i.e. hw /x = 60, and gives rise to
70 800 fluid particles and 6576 boundary particles. The fluid 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 fluid 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 fill ratio
[m5G;August 14, 2018;18:50]
15
of the tank and the extended durations of the experiment at each
frequency.
Regarding the fill 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 fluid height with sufficient 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 significant
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 difficulties, the proposed SPH scheme was capable of reasonably simulating the experiment, as described in detail in the following paragraphs.
The first case corresponds to a frequency fw1 = 0.32 Hz which
is lower than that of the first 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 fluid 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 fluid 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 first 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 fluid 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 first 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 fill ratio and high stretching,
Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.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 specifically, 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 fluid 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 flow 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 flow 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 fill ratio and high stretching,
Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.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 figure 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 figure 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 fill ratio and high stretching,
Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.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 fluid particles and
904 821 boundary particles. This modest resolution with 20 particles in the fluid 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 fluid 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 artificial viscosity was added. The free-slip condition can
be used here because the no-slip condition has negligible effect
on the fluid flow 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 fluid 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 artificial 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 fill ratio and high stretching,
Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.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 first simulating is using artificial
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 first simulating is using artificial 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 artificial 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 fill ratio and high stretching,
Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.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 field and the sloshing force. The simulation without
δ -SPH shows a much nosier density field, 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
fill 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 efficiency and sensitivity to
round-off error.
The treatment of solid wall boundary conditions using the fixed
ghost particle method of Adami et al. [15] includes a force balance correction which greatly improves the smoothness of the fluid
properties adjacent to the solid walls. It is also simple to implement and numerically efficient. 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 fluid,
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 identified as the most suitable form of numerical dissipation for long duration simulations. The hydrostatic and sloshing verification 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 efficient 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 fill-ratio sloshing in tank
showed that the use of single precision, as commonly favoured in
GPU implementations, resulted in significant 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 significant 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
five to seven-fold increase is computational time, depending on the
specific 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 five times more double-precision floating point operations per
second, and over three times the memory bandwidth. We estimate
that the code could potentially benefit 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 fission 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 flows 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 flows. 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 flows. 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 flows. 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 flows using
SPH and the unified semi-analytical wall boundary conditions. Comput Phys
Commun 2013;184(11):2515–27.
[17] Ferrand M, Laurence D, Rogers B, Violeau D, Kassiotis C. Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows 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 fill ratio and high stretching,
Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.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 flows. Phys Fluids (1994–present)
2006;18(10):101702.
[21] Bonet J, Kulasegaram S. A simplified 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 flows 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 flows. Comput Fluids 2009;38(6):1203–17.
[26] Antuono M, Colagrossi A, Marrone S, Molteni D. Free-surface flows 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 definite 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 fluid
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 fluid 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 scientific 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 fill ratio and high stretching,
Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.07.006
Документ
Категория
Без категории
Просмотров
0
Размер файла
6 691 Кб
Теги
2018, 006, compfluid
1/--страниц
Пожаловаться на содержимое документа