- Email: [email protected]
A Smoothed Particle Hydrodynamics approach for thermo-capillary flows
Accepted Manuscript
A Smoothed Particle Hydrodynamics Approach for Thermo-Capillary Flows Manuel Hopp-Hirschler, Mostafa Safdari Shadloo, Ulrich Nieken PII: DOI: Reference:
S0045-7930(18)30640-6 https://doi.org/10.1016/j.compfluid.2018.09.010 CAF 4005
To appear in:
Computers and Fluids
Received date: Revised date: Accepted date:
24 March 2018 4 September 2018 13 September 2018
Please cite this article as: Manuel Hopp-Hirschler, Mostafa Safdari Shadloo, Ulrich Nieken, A Smoothed Particle Hydrodynamics Approach for Thermo-Capillary Flows, Computers and Fluids (2018), doi: https://doi.org/10.1016/j.compfluid.2018.09.010
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights
flow driven by gradients of the surface tension.
CR IP T
• Smoothed Particle Hydrodynamics (SPH) model for thermo-capillary
• Continuum surface force (CSF) approach including Marangoni forces.
AC
CE
PT
ED
M
AN US
• Convergence study for thermos-capillary flow
1
ACCEPTED MANUSCRIPT
A Smoothed Particle Hydrodynamics Approach for Thermo-Capillary Flows
CR IP T
Manuel Hopp-Hirschler1a , Mostafa Safdari Shadloob , Ulrich Niekena a
AN US
Institute of Chemical Process Engineering, University of Stuttgart, 70199 Stuttgart, Germany b CORIA-UMR 6614, CNRS-University, INSA of Rouen and Normandie University, 76000 Rouen, France
Abstract
Interfacial-driven flows are important phenomena in many processes. In this article, we present a Smoothed Particle Hydrodynamics (SPH) model for thermo-capillary flow driven by gradients of the surface tension. The model
M
is based on the continuum surface force (CSF) approach including Marangoni forces. An incompressible SPH approach using (i) density-invariant divergence-
ED
free (DIDF), (ii) corrected SPH and (iii) particle shifting approaches for multi-phase systems is used for accurate results.
PT
We carefully validate the proposed model using several test cases. First, we demonstrate the effects of corrected SPH and particle shifting approaches
CE
using Taylor-Green vortex. Then, we study single-phase flow problems to validate correct implementation of boundary conditions, momentum and energy
AC
balance using lid-driven cavity, diffusive transport problem, and buoyancydriven cavity test cases. Afterward, we investigate different multi-phase flow problems to validate nor1
corresponding author: [email protected]
Preprint submitted to Computers & Fluids
September 13, 2018
ACCEPTED MANUSCRIPT
mal and tangential component of the surface tension. Finally, we apply the model to thermo-capillary rise of a droplet due to a temperature gradient.
CR IP T
We present a convergence study and compare the results with their counterparts obtained from OpenFoam software as well as Finite Volume method
(FVM) reference from literature. We demonstrate that the proposed model is very accurate for thermo-capillary flow. The simulation results of the current SPH approach will be available online for the community.
AN US
Keywords: Thermal Convection, Marangoni Effect, Multiphase Flow, Smoothed Particle Hydrodynamics (SPH), Meshless Method 1. Introduction
2
Smoothed particle hydrodynamics (SPH) is a meshless Lagrangian particle
3
method to solve partial differential equations. Tough, initially introduced for
4
astrophysical problems, by Lucy [1] and Monaghan and Gingold [2] in 1977,
5
SPH gained popularity in 1990s for fast dynamic flows [3, 4, 5], in 2000s for
6
multi-phase problems [6, 7, 8] and lately in diverse engineering applications
7
[9, 10].
8
The main reason for this popularity, especially for multi-phase flows, is its
9
inherent strengths which some of those may be found in other methods but
10
the combination of all those features can be found only in SPH [9]. Among
11
others one can mention, firstly, the natural distinguishing between different
AC
CE
PT
ED
M
1
12
phases due to holding material properties at each individual particle. This
13
feature helps to extend a code for more than two-phase flows, such as those
14
in liquid-liquid-gas (such as oil-water-gas) and liquid-liquid-solid (oil-water-
15
salt/sediment) [11]. Additionally, SPH inherently captures the interface and 3
ACCEPTED MANUSCRIPT
is not tracking it like those in Level-set method. Therefore, the interface is
17
always sharp, is not smeared, the mass of each phase is rigorously conserved
18
during the simulation and no need for complex coding for large topological
19
change of interface (such as those in folding, merging and break up). This is
20
very important especially when the code is extended for 3D problems. Sec-
21
ondly, SPH naturally incorporate singular forces and coefficient discontinu-
22
ities into the numerical scheme, therefore, the interface boundary condition
23
is inherently implemented. Thirdly, since the field properties’ derivatives
24
will transform to derivatives of some analytical functions, so-called kernel
25
function, they have exact solutions and can be naturally integrated into the
26
scheme. And finally, the convective term in discretisation of the balance equa-
27
tions is absent in the numerical approximation scheme due to its Lagrangian
28
nature. The latter makes SPH a suitable candidate for fast-dynamic flows.
29
For multi-phase flows, the most important interfacial force is the surface ten-
30
sion. It is known that the surface tension may depend on the temperature
31
or concentration. If there is a temperature or concentration gradient tan-
32
gential to the interface between two fluids, then, an additional force known
33
as Marangoni force, arises. However, and in spite of numerus articles on
34
multiphase flow problems in the context of SPH [12, 13, 14, 15], in literature
35
there are only a few publications on Marangoni force (i.e. the gradient of the
36
surface tension) [16, 17]. Marangoni force due to a concentration gradient
CE
PT
ED
M
AN US
CR IP T
16
at the interface was first discussed by Adami et al. [16]. They used the
38
conservative formulation of the continuum surface force (CSF) model, the
39
so-called continuum surface stress (CSS) model, and investigated the motion
40
and surface concentration of a bubble.
AC 37
4
ACCEPTED MANUSCRIPT
In the context of SPH, thermocapillary motion was only recently studied by
42
Tong et al. [17] for 2D problems. Tong et al. presented results for the migra-
43
tion velocity of a single droplet in thermocapillary flow. In this case a droplet
44
is accelerated by gradients of the surface tension coefficients. For a very low
45
resolution, the migration velocity presented by Tong et al. match results
46
of Ma & Bothe [18] but they didn’t provide a convergence study. Tong et
47
al. [17] presented a model based on the non-conservative formulation of the
48
CSF model (the same that we use here). They investigated thermocapillary
49
migration of a droplet where heat transfer was included using an enthalpy
50
balance. Our model is similar to the model of Tong et al. [17] except that
51
we use the energy balance in temperature form and neglect phase transition.
52
Therefore, the current work aims at the proposition and the development of
53
a consistent SPH method that is capable of treating complex flow physics in-
54
cluding Marangoni forces and thermocapillary effects. At the same time this
55
article provides step by step validation test cases that can serve as bench-
56
marks for both SPH and thermo-fluids communities. It is noted that unlike
57
Adami et al. [16] the proposed scheme is based on the non-conservative for-
58
mulation of the CSF model. Additionally, this manuscript extend the results
59
of Tong et al. [17] to 3D problems and it provides more details on the con-
60
vergence of SPH results when dealing with thermocapillary forces.
61
This article is organized as follows: First we introduce the balance equa-
CE
PT
ED
M
AN US
CR IP T
41
tions of the multiphase system. Then we introduce the SPH method and
63
the discrete balance equations, including boundary conditions. Afterwards
64
numerical stability is discussed introducing kernel consistency, particle con-
65
sistency and volume conserving ISPH scheme. Numerical results are divided
AC 62
5
ACCEPTED MANUSCRIPT
into single-phase and multi-phase sections. First, numerical implementation
67
is validated for single-phase flow. Afterwards, surface tension is validated
68
for multi-phase flow and results for thermo-capillary droplet migration are
69
presented and validated against analytical solution and related numerical
70
solutions.
71
2. Balance equations
72
In this section we present the balance equations. A complete thermodynamic
73
consistent form of the component, momentum and energy balance, as found
74
in Thieulot et al. [19], is very complex and has parameters that are not
75
known for the system under consideration. Therefore, we make the following
76
assumptions to only reduce complexity without neglecting important mech-
77
anisms.
78
Firstly, we assume that the fluids are immiscible and incompressible. Sec-
79
ondly, we neglect phase transitions and assume that the surface tension de-
80
pends linearly on temperature. Using these assumptions we introduce the
81
balance equations. Then, we review these derivations in Lagrangian form for
82
our Lagrangian numerical framework.
83
2.1. Continuity equation
84
The continuity equation is [20]
AC
CE
PT
ED
M
AN US
CR IP T
66
dρ + ∇ · (ρ~u) = 0, dt
(1)
85
where ρ and u are the mass density and the mass average velocity of the fluid,
86
respectively. We formulate the continuity equation in a Lagrangian reference
6
ACCEPTED MANUSCRIPT
frame. Therefore, we rewrite Eq. 1 using the material derivative Dρ + ρ∇ · ~u = 0, Dt
88
D Dt
(2)
CR IP T
87
We define an incompressible system with constant density as ∇ · ~u = 0,
(3)
where the density is independent of the pressure, hence, the divergence of
90
the velocity is zero.
91
2.2. Momentum equation
92
The momentum balance of the multi-phase system is [19]
AN US
89
(4)
M
d (ρ~u) + ∇ · (ρ~u~u) = −∇p + ∇ · Πviscous + F~body + ∇ · Πcapillary , dt
with the pressure p, the viscous stress tensor Πviscous , a body force, e.g.
94
gravity F~body = ρg, and the capillary stress tensor Πcapillary . The first three
95
terms on the right hand side are classical terms in the momentum balance
96
as in single-phase flow [20]. The last term in Eq. 4 accounts for the presence
97
of interfaces. With the continuity equation (Eq. 1) the momentum balance
98
for an incompressible, Newtonian fluid reduces to
CE
PT
ED
93
d~u + ρ (~u · ∇) ~u = −∇p + η∇2~u + F~body + ∇ · Πcapillary , dt
(5)
AC
ρ
99
with the dynamic viscosity η.
100
The capillary stress tensor in Eq. 5 arises due to the presence of the interface,
101
which in general is diffuse on molecular scale because density continuously
102
varies between different phases. Here, we consider the sharp interface limit 7
ACCEPTED MANUSCRIPT
of the capillary stress tensor. The sharp interface limit may be seen as
104
projecting the integral over the diffuse interface on the interface area (plane
105
in 3D). A review of a mathematical derivation of the sharp interface limit
106
from a diffuse formulation is shown in [21]. The main result is
CR IP T
103
∇ · Πcapillary = σκ~n ˆ + ∇S σ δ,
(6)
with the surface tension coefficient σ, the curvature κ and the gradient of
108
the surface tension tangential to the surface ∇S σ. Additionally, ~n ˆ is the unit
109
normal to the interface and δ is the Dirac-Delta function
δ=
AN US
107
1 at the interface 0 in the bulk
.
(7)
Finally, the momentum equation in Lagrangian form with a sharp interface
111
is
ρ
ED
M
110
D~u = −∇p + η∇2~u + F~body + σκ~n ˆ + ∇S σ δ. Dt
(8)
2.3. Energy balance
113
The energy balance in temperature form, neglecting viscous dissipation and
114
surface energy contributions, is
AC
CE
PT
112
ρcp
DT = ∇ (λ∇T ) , Dt λ , ρcp
115
with the thermal diffusivity of k =
116
the specific heat capacity at constant pressure, cp .
8
(9)
the thermal conductivity of λ and
ACCEPTED MANUSCRIPT
3. Smoothed Particle Hydrodynamics method
118
As mentioned earlier, SPH method was first introduced by Monaghan and
119
Gingold [2] and Lucy [1] in 1977. Since then, several reviews are published
120
presenting the basics of SPH. For a general introduction to the SPH method
121
and its various applications, we refer to the reviews of [22, 10, 9]. Here, we
122
briefly present the SPH initiatives and the relevant discretizations that is
123
used in the current study.
124
The starting point of SPH is coming from an exact mathematical solution
125
where we can estimate any quantity of A (~x) at a given position ~x by mul-
126
tiplying it with a Dirac-Delta function and taking its integration A (~x0 ) over
127
the whole volume [2]
A(~x0 )δ (|~x − ~x0 |) d~x0 ,
M
A (~x) =
Z
AN US
CR IP T
117
V
(10)
where ~x0 is the position of any point in the volume V . The Dirac-Delta
129
function is defined to be unity for ~x0 = ~x and zero otherwise.
130
In the next step, this analytical equation undergoes two approximations. The
131
first, known as kernel approximation, is done by replacing the Dirac-Delta
PT
function with a weighting function, so-called kernel function W (h, |~x − ~x0 |) Z A (~x) = A(~x0 )W (h, |~x − ~x0 |) d~x0 , (11)
CE
132
ED
128
V
with h being the radius of the support volume, so-called smoothing length.
134
In our case, where an initial regular particle distribution in a Cartesian order
135
with the spacing of L0 is used, we define a constant value of h = 2.05L0 .
136
In this work we use a Wendland C2 kernel function of order 4 with rc = 2h
AC
133
9
ACCEPTED MANUSCRIPT
[23]
q 4 1 − (1 + 2q) for 0 ≤ q ≤ 2 fw 2 W (h, |~x − ~x0 |) = d . h 0 for 2 < q
(12)
CR IP T
137
138
Here, d is the dimension and fw is a normalization constant with fw = 43 ,
139
7 , 21 4π 16π
140
the non-dimensional smoothing length [24]. The first derivative of this kernel
141
function is
for d = 1, d = 2, and d = 3 dimensions, respectively. q =
for 0 ≤ q ≤ 2
AN US
−5q 1 − q 3 2
∂W (h, |~x − ~x0 |) fw = d+1 ∂x h 0
.
|~ x−~ x0 | h
is
(13)
for 2 < q
Next is the particle approximation where the integration in Eq. 11 is replaced
143
by a summation over finite number of interpolation points, known as particles
144
A (~x)a =
M
142
X mb ρb
A(~xb )W (h, |~xa − ~xb |) .
(14)
ED
b
Here, the indexes a and b represent the particle of interest and the particles
146
in the vicinity of a (see Fig. 1) called neighbor particles. mb and ρb are
147
the mass and density of particle b, respectively. In the following the kernel
148
function is abbreviated as W (h, |~xa − ~xb |) = Wab with |~xa − ~xb | = rab being
149
the distance between particle a and b.
150
An example, using Eq. 14 with A = ρ, is the calculation of the density of
151
particle a [22]
AC
CE
PT
145
152
153
ρa =
X
mb Wab .
(15)
b
Note that particle a should also be considered in its vicinity because Waa 6= 0. Due to the particle approximation, the 0th moment (unity condition of the 10
CR IP T
ACCEPTED MANUSCRIPT
AN US
Figure 1: (Left) Scheme of discrete particle distribution around a particle of interest,
particle a, and its particles in the vicinity b within the support of the kernel function. Dashed particles are outside the support domain. (Right) Schematic representation of a kernel function W with a support domain of 2h.
kernel function) is not guaranteed [10]. To guarantee it, we may use a renor-
155
malization function, the shepard kernel, to ensure the unity condition [25, 26].
156
This is very important if the support domain is truncated e.g. at free sur-
157
faces or when surface tension has to be calculated as we will see later. For
158
bulk flow the error is normally in the range of a few percent. We discuss this
159
point in section 5.1.
160
As mentioned before, one strong advantage of the SPH method is that the
161
gradient of a quantity A (~x) can be derived from the gradient of the kernel
162
function [2]. One can formulate the first derivative using [22]
AC
CE
PT
ED
M
154
163
or [6]
X mb
∇A (~x)a =
1 ∇A (~x) ρ
b
a
=
ρb
(Ab − Aa ) ∇W (h, rab ) ,
X
mb
b
Ab + Aa ∇W (h, rab ) , ρa · ρb
11
(16)
(17)
ACCEPTED MANUSCRIPT
164
where the gradient of the kernel function ~rab ∂W , rab ∂x
(18)
CR IP T
∇W (h, rab ) =
is analytically known from Eq. 13. Here, ~rab is the distance vector between
166
particle a and b pointing from particle a to b and rab is the length of this
167
vector. Eq. 16 guarantees that the gradient of a constant field is exactly zero
168
even if numerical errors are present. Linear momentum is only conserved if
169
the 1th moment of the kernel function is exactly zero [10]. Eq. 17 conserves
170
linear momentum exactly as long as the kernel function is symmetric but
171
the gradient of a constant field is only zero if the 1th moment of the kernel
172
function is exactly zero. In literature Eq. 16 is often used for the divergence
173
of a velocity field, where the gradient of a constant field is very important,
174
and Eq. 17 is often used for the pressure gradient because linear momentum
175
is more important in this case.
176
For the second derivative of a quantity, as introduced by Brookshaw [27], a
177
hybrid formulation for the second derivative is used with a first order deriva-
178
tive of the kernel function in combination with a finite difference approxima-
179
tion
(19)
where α is a variable that does not depend explicit on position but e.g.
AC
180
X mb ~rab ~ (~x) = ~a − A ~b , ∇ · α∇A (αa + αb ) 2 ∇a Wab · A ρb rab a b
CE
PT
ED
M
AN US
165
181
may depend on other properties of a fluid. We abbreviate the gradient with
182
respect to particle a as ∇a Wab = ∇a W (h, rab ). In contrast to other formula-
183
tions of the second derivative [28] this formulation proved to be less sensitive
184
to particle disorder. 12
ACCEPTED MANUSCRIPT
4. Discrete balance equations
186
In this section we present the discrete balance equations using the operators
187
from section 3, followed by presenting the projection method used in this
188
work to calculate the pressure.
189
4.1. Continuity equation
190
In SPH, there exist two different ways [22] to formulate the discrete continuity
191
equation. One way is to directly use the discrete continuity equation [29, 30]
192
using Eq. 16
AN US
CR IP T
185
X mb Dρa (~ub − ~ua ) ∇a Wab . = ρa Dt ρ b b
(20)
A disadvantage of this equation is the error of conservation due to integration
194
errors. Another formulation of the discrete continuity equation is [12]
M
193
ED
ρ a = ma
X
Wab ,
(21)
b
where the density is directly calculated using the mass and volume of par-
196
ticle a. Note that the summation of the kernel function represents the re-
197
ciprocal volume 1/Va . Advantages of this formulation are its applicability
198
to multi-phase systems with density differences between the phases and its
199
conservation property. Therefore, we use Eq. 21 throughout this work.
200
4.2. Momentum balance
AC
CE
PT
195
201
In the momentum balance (Eq. 8) we need first-order derivatives for the
202
pressure term, for the color to calculate the normals, for the normals to
203
calculate the curvature and for the surface tension to calculate the Marangoni
204
force. In addition, we need a second-order derivative for the divergence of 13
ACCEPTED MANUSCRIPT
205
the viscous stress tensor.
206
We use the positive formulation of the gradient (Eq. 17) for the pressure
207
force
∇p ρ
a
=
X mb (pb + pa ) ∇a Wab , ρ a ρb b
CR IP T
(22)
208
with the pressure pa and pb of particle a and b, respectively. It is common
209
in literature to use the positive formulation of the gradient especially for
210
the pressure force to ensure linear momentum conservation. For a constant
213
214
AN US
212
pressure field the gradient of the pressure is only zero for a constant pressure P p = 0 or symmetric particle distribution, which results in b ∇a Wab = 0. The discrete formulation of the viscous term using the second-order derivative (Eq. 19) is
η 2 ∇ ~u ρ
=
a
X mb ~rab (ηa + ηb ) 2 ∇a Wab · (~ua − ~ub ) . ρa ρb rab b
(23)
M
211
The body force is straight forward since it only depends on particle a.
216
The last two terms in Eq. 8 account for interface forces. To formulate surface
217
tension at a sharp interface, Brackbill et al. [31] and later, in the context of
218
SPH, Morris [32] introduced the so-called Continuum Surface Force (CSF)
219
approach. In the CSF approach the surface force, that is a force per area, is
220
reformulated as a force per volume. Therefore, the surface tension is numer-
221
ically smoothed out around the sharp interface. We need the normal vector
222
~n at the interface between phases to calculate the magnitude and direction
223
of the interface forces. The normal vector at an interface is calculated from a
224
color field where a unique color is assigned to particles of every phase. We use
225
the negative formulation of the gradient (Eq. 16) to calculate the gradient of
226
the color field because it is important that the gradient vanishes if the color
AC
CE
PT
ED
215
14
ACCEPTED MANUSCRIPT
228
is constant. Otherwise, we get surface tension due to numerical artifacts in the bulk. Hence, interface normal vector ~na is calculated as X mb (Cb − Ca ) ∇C ~na = = ∇a Wab , [C] a ρb [C] b
(24)
CR IP T
227
with the color C and its jump across the interface [C]. Whether the color is a
230
smooth or a sharp function across the interface, the jump of the color function
231
represents the absolute difference between the maximum and minimum of the
232
color function. The normalized normal vector is
AN US
229
~n ~n , ˆ= |~n| 233
(25)
with the magnitude of the normal vector |~n|. Then, the curvature is (26)
M
X mb ~n ˆ b − ~n ˆ a ∇a Wab . κa = − ∇ · ~n ˆ =− ρb a b
We again use the negative formulation because it is important that the cur-
235
vature is zero for equal normal vectors. Note that we truncate the field of
236
normals to avoid numerical artifacts and only consider normals with magni-
237
tude larger than /h with = 0.01.
238
The last term in Eq. 8 accounts for gradients of the surface tension tangential
239
to the interface. We use a projection of the gradient of the surface tension
240
into the tangential direction of the interface
AC
CE
PT
ED
234
241
(∇S σ)a = ∇σ − ∇σ · ~n ˆ ~n ˆ , a
(27)
with the gradient of the surface tension coefficient σ (∇σ)a =
X mb b
ρb
(σb − σa ) ∇a Wab .
15
(28)
ACCEPTED MANUSCRIPT
242
The delta function in Eq. 8 is unity at the interface. In the discrete formu-
243
lation we smooth out the delta function (29)
CR IP T
δa = |~na |.
This leads to a distribution of the surface force in a volume that in the limit is equal to that of a delta function.
Finally, the used discrete momentum balance in the current work is
with f~body =
245
4.3. Energy balance
M
~body F . ρ
244
The discrete energy balance for convective and diffusive heat transfer is X mb ~rab DT = (ka + kb ) 2 ∇a Wab · (Ta − Tb ) , (31) Dt a ρ r b ab b
PT
ED
246
AN US
X mb X mb ~rab D~ua =− (pb + pa ) ∇Wab + (ηa + ηb ) 2 ∇a Wab · (~ua − ~ub ) Dt ρa ρb ρa ρb rab b b |~na | σa κa~n ˆ a + (∇S σ)a , (30) + f~body,a + ρa
using the second-order derivative (Eq. 19) where Ta and λa are the temper-
248
ature and heat conductivity of particle a.
249
4.4. Estimation of pressure
250
The pressure in the momentum balance is still unknown. For incompress-
251
ible fluids there are two ways to calculate the pressure. One way is to use
252
an equation of state for a so-called weakly compressible fluid and is known
253
as WCSPH. It is a stiff equation that relates a small density variation to a
254
pressure variation. One example is the Tait equation [33]. This equation is
AC
CE
247
16
ACCEPTED MANUSCRIPT
commonly used in SPH to simulate the pressure of water [6, 22].
256
An alternative way to determine the pressure for a truly incompressible fluid
257
is to use a projection method. We can apply the Helmholtz-Hodge decom-
258
position [34] on the momentum balance. The momentum balance is divided
259
into a curl-free and divergence-free part. In a predictor step we first calculate
260
an intermediate solution of velocity and density using the curl-free part. The
261
curl-free part of the momentum balance (Eq. 30) consists of all terms except
262
the pressure term. In a second step we solve a Pressure Poisson Equation
263
(PPE) to calculate the pressure and then correct the intermediate velocity
264
using the pressure term to get the final velocity of the time step. In SPH
265
this projection method was introduced by Cummins and Rudman [3] and is
266
usually called ISPH. The details of the projection method are discussed in
267
literature [8, 26, 29, 35] and we only review the discrete equations used in
268
this work.
AN US
M
The curl-free part acceleration ~aa,∗ in the momentum balance is X mb ~rab |~na | ~aa,∗ = (ηa + ηb ) 2 ∇a Wab ·(~ua − ~ub ) +f~body,a + σa κa~n ˆ a + (∇S σ)a . ρa ρb rab ρa b (32)
PT
ED
269
CR IP T
255
Here, the subscript asterisk indicates the intermediate step. We calculate the
271
intermediate velocity using an explicit Euler scheme
CE
270
(33)
with the time step ∆t. Then, we determine the new particle position
AC 272
~ua,∗ = ~uta + ~aa,∗ ∆t
~xa,∗ = ~xta + ~ua,∗ ∆t,
(34)
273
also using an explicit Euler step. The intermediate particle position is only
274
necessary because we calculate the intermediate density using Eq. 21 instead 17
ACCEPTED MANUSCRIPT
of the continuity equation (Eq. 20). When the intermediate density is calcu-
276
lated we put the particles back to the old particle position ~xta and calculate
277
the pressure using the PPE ∇·
279
280
∇p ρ∗
=
∇ · ~u∗ . ∆t
(35)
We use the divergence of the velocity as a condition for incompressibility. The discrete Laplace of the pressure using Eq. 19 is (for α = 1) X mb ~rab ∇p · (pa − pb ) 2 ∇a Wab , = 2 ∇· ρa,∗ a ρb ρa rab b
AN US
278
CR IP T
275
and the right hand side is
1 X mb ∇ · ~ua,∗ = (~ub,∗ − ~ua,∗ ) ∇a Wab . ∆t ∆t b ρb
(36)
(37)
The PPE requires the solution of a linear equation system. We use a sta-
282
bilized version of the biCGStab method from the PETSc library [36] with a
283
boomerAMG pre-conditioner from the HYPRE library [37].
284
Then, with the acceleration of the corrector step ~aa,∗∗ , using the pressure
285
term
ED
M
281
(38)
~ut+1 = ~ua,∗ + ~aa,∗∗ ∆t. a
(39)
The final position of particle a is
AC
287
X mb (pb + pa ) ∇a Wab , ρ a · ρb b
we determine the final velocity with an explicit Euler step
CE
286
PT
~aa,∗∗ = −
= ~xta + ~xt+1 a
~uta + ~ut+1 a ∆t. 2
(40)
288
Here, a semi-implicit step with the velocities at time t and t+1 is used. More
289
details on the projection method in SPH is found in Keller [38]. 18
ACCEPTED MANUSCRIPT
4.5. Boundary conditions
291
Throughout this paper we use impermeable boundary conditions. For veloc-
292
ity we either use slip or no-slip conditions. For temperature and pressure we
293
use Dirichlet or Neumann boundary conditions, while for the gradient of the
294
color C only Neumann boundary conditions are used, because the fictitious
295
boundary condition should not influence the gradient of C. As we will see in
296
Sec. 5.1, we do not need boundary conditions for the normal and the surface
297
tension coefficient.
298
Beside different approaches for boundary treatment in SPH [7, 29, 39, 40,
299
41, 42, 43], we use wall boundary conditions from Cummins and Rudman
300
[3] (mirror particles) at straight walls. Due to the integral character of SPH
301
the wall partially needs to be discretized up to a certain distance from the
302
fluid domain. In the mirror particle approach, the fluid particles are mir-
303
rored in every time step across the boundary. In this way an image of the
304
fluid particles represent the wall particles. This is sketched in Fig. 2. The
305
properties of the wall particles are manipulated to incorporate Dirichlet and
306
Neumann boundary conditions. The boundary condition for any property
307
Ψ ∈ [~u, p, T, C] at the wall is considered in the property of an image particle.
308
The property of an image particle is ~ 0a = 2Ψ ~ w − δw Ψ ~a Ψ
AC
CE
PT
ED
M
AN US
CR IP T
290
309
~ w and the sign with the value of the wall Ψ +1 for Dirichlet boundary conditions δw = −1 for Neumann boundary conditions with Ψ ~w = 0 19
(41)
.
(42)
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 2: Scheme of mirror boundary conditions. The fluid particles are mirrored across
PT
the boundary and the images represent the wall particles (dotted particles). Their properties are manipulated to incorporate Dirichlet and Neumann boundary conditions. Here
AC
CE
no-slip boundary conditions of an impenetrable wall is shown.
20
ACCEPTED MANUSCRIPT
More details can be found in [29, 42]. This technique works well on straight
311
walls but leads to errors on curved boundaries. Nevertheless, in the present
312
test cases we only need straight walls.
313
5. Numerical stability
314
It is known that the unity condition for the SPH kernel function is not guar-
315
anteed for random particle distributions. One way to guarantee the unity
316
condition of the kernel function was proposed by Bonet & Lok [25]. It is
317
called the corrected SPH approach and is presented in section 5.1.
318
In practical applications we observe that a corrected SPH approach is not
319
enough to guarantee stable and accurate simulations. The reason is that
320
in SPH we must satisfy two conditions. The first condition is the already
321
mentioned consistency of the kernel function. The second condition is the
322
consistency of the particle distribution [10]. If the spacing between particles
323
is larger than the interaction radius, the particles become disconnected and
324
even corrected SPH approaches can’t prevent that.
325
In a proper ISPH approach, particles that are initially on a regular Cartesian
326
grid, may move on the same trajectory. In front of a stagnation point, parts
327
of the particles move in different directions and get disconnected when their
328
distance is larger than the cut-off radius. As a result, a large void space is
329
formed between them. Therefore, we need a numerical correction that rear-
AC
CE
PT
ED
M
AN US
CR IP T
310
330
ranges the particles if necessary without influencing the balance equations.
331
In literature there exist two approaches to improve the particle consistency
332
in ISPH. The first approach was introduced by Hu & Adams [35]. It is the
333
so-called combined DIDF. The other approach was introduced by Xu et al. 21
ACCEPTED MANUSCRIPT
[44] and Shadloo et al. [43], used for the first time in the concept of mul-
335
tiphase flows by Shadloo et al. [45, 46] and modified to consider its arts
336
effects for interfacial flows by Lind et al. [47]. It is called the particle shifting
337
technique. Using these approaches improve the particle consistency and also
338
the 0th and 1th moment of the kernel function and the gradient of the kernel
339
function. We summerize both approaches in section 5.2 and 5.3.
340
5.1. Corrected kernel function and corrected gradient of kernel function
341
In this section we present the corrected SPH approach, introduced by Bonet
342
& Lok [25] as mixed kernel and gradient correction. Further details about
343
the implementation and accuracy are found in [38]. The idea of the corrected
344
SPH approach is to renormalize the kernel function and the gradient of the
345
kernel function to guarantee the unity condition for the kernel and the 1th
346
moment
M
AN US
CR IP T
334
X mb
∇a Wab = 0,
(43)
˜ ab = Wab , W Sa
(44)
X mb
(45)
with the Shepard kernel Sa [48]
AC
CE
348
ρb
for the gradient of the kernel. The renormalized kernel function is
PT
347
ED
b
Sa =
b
ρb
Wab ,
349
and the uncorrected kernel function Wab . In the nominator the interaction
350
between all particles including the interaction of particle a with itself need
351
to be considered.
352
In most of the terms in the balance equations we need the gradient of the 22
ACCEPTED MANUSCRIPT
353
kernel function. In the mixed kernel and gradient correction of Bonet & Lok,
354
the corrected gradient of the corrected kernel is
355
(46)
CR IP T
˜ aW ˜ ab = La ∇a W ˜ ab , ∇
with the correction matrix, that ensures symmetric vector direction, !−1 X mb ˜ ab ⊗ ~rab ∇a W , La = ρ b b
(47)
and the corrected gradient of the kernel, that ensures equal magnitude of the
357
vector,
AN US
356
˜ ab = ∇a W 358
with γ=
∇a Wab − γ , Sa
P
mb b ρb ∇a Wab
Sa
.
(48)
(49)
An advantage of the corrected SPH approach is that a constant pressure field
360
has zero gradient because of Eq. 43. In all previous discrete formulations in
361
section 4 we may replace the kernel or gradient of the kernel function by Eqs.
362
44 or 46, respectively, if the corrected SPH approach is used. We explicit
363
exclude Eq. 24 because a corrected SPH kernel leads to wrong normals.
364
In addition we may use the corrected SPH approach if we only consider a
365
subset of the particles in the domain. For the curvature and the gradient of
366
the surface tension we use only particles with normal vectors of magnitude
367
larger than a threshold /h > 0.01 (see section 4.2). In all other cases we use
AC
CE
PT
ED
M
359
368
all particles in the vicinity, including mirrored particle at the boundary.
369
5.2. Density-invariant divergence-free (DIDF) approach
370
Hu & Adams [35] introduced an approach to improve particle consistency of
371
the ISPH method for bulk flow problems. The idea is to split up the time 23
ACCEPTED MANUSCRIPT
integration of velocity and position into two parts in the corrector step of
373
the projection method in ISPH. In a first integration step the left side of the
374
continuity equation 1 Dρ = ∇ · ~u, ρ Dt
375
is used as criterion for incompressibility in the PPE t τ2 ∇p ρ0 − ρ∗,t+1 ∇· = a 0a , 2 ρ∗ a ρa
CR IP T
372
(50)
(51)
where τ is the time step and ρ0a is the initial density of particle a. In this
377
first step, we calculate a pressure to correct the position of the particles on
378
the half time step using an explicit Euler step. The intermediate velocity is
AN US
376
~ua,∗ = ~uta + 0.5 · a~a,∗ ∆t, and the intermediate position is
M
379
380
In a second integration step, the divergence of the velocity is used as criterion for incompressibility in the PPE t ∇p τ∇ · = ∇a · ~ut+1 ∗ . ρ∗ a
(54)
CE
PT
381
(53)
ED
~xa,∗ = ~aa,∗ ∆t.
(52)
In this second step a pressure to correct the velocity of the particles is cal-
383
culated. The final velocity is
AC
382
384
~ut+1 = ~ua,∗ + ~aa,∗∗ ∆t, a
(55)
~xt+1 = ~uta + ~ut+1 aa,∗∗ (∆t)2 . a ∆t + 0.5 · ~ a
(56)
and the position is
24
ACCEPTED MANUSCRIPT
In the DIDF approach we need to solve the PPE two times per time step.
386
Therefore, the numerical effort increases drastically because the solution step
387
of the LES takes approximately 70 % of the hole computation time. Further-
388
more, the approach is not applicable to free surface flows [49]. On the other
389
hand the DIDF approach improves the conservation of volume because both
390
particle position and velocity are considered as criterion for incompressibility.
391
The DIDF approach is used here to improve volume conservation during
392
the simulation. Nevertheless, particles still remain on their trajectory and
393
split up at stagnation points. Therefore, another approach to regularize the
394
particles is necessary.
395
5.3. Particle Shifting
396
Particles in SPH move on their streamline due to their Lagrangian nature.
397
Therefore, particles split near a stagnation point leaving a void space between
398
two streamlines. This leads to large errors in the particle consistency of SPH.
399
The particle shifting technique is introduced to geometrically shift particles
400
slightly away from their streamlines into a void space [44, 43]. As a result
401
the particles are more regularly distributed in the domain. We introduce the
402
shifting algorithm in the following.
403
The idea is to shift particles slightly into a void space after every time step.
404
The total shifting vector ~ra,s is
AC
CE
PT
ED
M
AN US
CR IP T
385
~ a, ~ra,s = CαR
(57)
405
with the arbitrary constant C = 0.04, the shifting magnitude α = |~umax |dt,
406
the magnitude of the maximum particle velocity |~umax |, the time step dt, and 25
ACCEPTED MANUSCRIPT
407
the shifting vector ~a = R
Nb X r¯a2 ˆ ~r , 2 ab rab b
(58)
where ~rˆab , rab and r¯a are the unit distance vector between particle a and b,
409
the distance between both particles and the average particle spacing in the
410
vicinity defined as Nb
1 X r¯a = rab . Nb b
CR IP T
408
(59)
In this approach it is assumed that all particles represent the same volume.
412
Therefore, a geometrically uniform distribution is the best particle distribu-
413
tion.
414
Some properties of the fluid are connected to the particles and their posi-
415
tion. If we shift the particles we need to correct the hydrodynamic variables
416
Ψ to be consistent and not influence the balance equations. Following [44],
417
a Taylor series of the hydrodynamic variable
ED
M
AN US
411
Ψa0 = Ψa + δ~raa0 · (∇Ψ)a + · · ·
(60)
is used to approximate the property of a particle at the new position a0 af-
419
ter shifting. We use a linear interpolation and neglect quadratic and higher
420
terms. Hydrodynamic variables used in this work are the velocity and tem-
421
perature. The density is not a hydrodynamic variable in this context because
422
we calculate the density by summation of masses in the vicinity instead of
AC
CE
PT
418
423
calculating the continuity equation directly. Therefore, we are able to calcu-
424
late the density at any point.
425
The particle shifting technique is independent of the time integration and
426
does not affect the balance equations because we shift the particles only at 26
ACCEPTED MANUSCRIPT
the very end of a time step and correct the particle properties accordingly.
428
The shifting technique in its original form assumes that every particle is an
429
identical interpolation point in the domain. In a multi-phase system this is
430
not valid any more because of the presence of an interface. In SPH the inter-
431
face between two phases is defined by the type (color) of particles. Therefore,
432
the interface is implicitly moved when particles change their position. If we
433
shift the particles, we implicitly shift the interface between phases as well. To
434
overcome this shortcoming we modify Eq. 57 near an interface between two
435
phases and apply shifting only in the tangential direction to the interface.
436
This can be done by only consider the tangential part of Eq. 57. A similar
437
extension was proposed by Lind. et al. [47]. They modified the shifting
438
technique for free surface flows based on a Fickian approach.
439
5.4. Time integration
440
In this work we use a predictor-corrector scheme to integrate the momentum
441
balance equation in time, as introduced earlier in the context of the projection
442
method to calculate the pressure (Sec. 4.4 and 5.2).
443
For the energy balance (Eq. 31) we use an explicit Euler integration scheme
444
during the predictor step.
445
For an explicit integration scheme there exist criteria for the maximum time
446
step because the algorithm is not unconditionally stable. Therefore, the
447
following time step criteria are used throughout the work. The time step
448
criteria due to convective transport (CFL criterion) is
AC
CE
PT
ED
M
AN US
CR IP T
427
∆tCF L = αCF L
27
L0 , umax
(61)
ACCEPTED MANUSCRIPT
450
due to viscous diffusion is ∆tvisc = αdif f
L20 νmax
∆theat = αdif f
L20 , λmax
and due to heat diffusion is
(62)
CR IP T
449
(63)
with αCF L = 0.05 and αdif f = 0.0625. Additionally, umax , νmax and λmax are
452
the magnitude of the maximum velocity, the maximum kinematic viscosity
453
and the maximum heat conductivity, respectively.
454
6. Single-phase flows
455
In this section we validate the momentum balance of the proposed model
456
and its discrete formulation using the SPH method. Since the SPH method
457
itself is well validated in literature we focus on the modifications made in
458
this work. After showing the numerical stability of the proposed method
459
and validation of the diffuse transport of a scalar property, we start with
460
the very common isotherm lid-driven cavity test case [50, 51] to demonstrate
461
accurate implementation of the momentum balance, boundary conditions
462
and numerical stability approaches. It is one of the standard test cases of
463
computational fluid dynamics simulation codes. It enables us to validate no-
464
slip boundary conditions for velocity as well as the wall velocity.
465
Afterwards, we couple the momentum equation with the energy equation.
466
We consider a buoyancy-driven flow in a cavity where the density of the fluid
467
is temperature-dependent. In this test case a circular flow is induced due to
468
a temperature gradient [52, 53, 54]. Note that the validation of diffusive heat
469
transport in shown in section 6.3 for completeness.
AC
CE
PT
ED
M
AN US
451
28
ACCEPTED MANUSCRIPT
6.1. Numerical stability
471
Here, we first compare the three previously introduced approaches (i.e. cor-
472
rected kernel function and corrected gradient of kernel function, abbreviated
473
as corrected SPH in the following, DIDF and particle shifting approach)
474
for numerical stability and combinations of them. One common test case
475
in literature is a two-dimensional Taylor-Green vortex. In this test case all
476
boundaries are periodic and initially a divergence-free velocity profile [55]
AN US
ux = sin(2πx)cos(2πy)
CR IP T
470
477
uy = −cos(2πx)sin(2πy)
(64)
(65)
is imposed. The vortex decays due to viscous dissipation of the fluid. Details
479
about this test case are found in [44, 47]. We consider a case with Reynolds
480
number Re = 1000.
481
The magnitude of velocity and the particle distribution are shown in Fig. 3
482
at time t = 0.24s for different combinations of the previous approaches. Fig.
483
3a shows the results without any numerical stability. We see that particles
484
are aligned on their streamline and large void spaces are formed between
485
the streamlines. The other three cases in the top line (Fig. 3b-d) are cases
486
where the corrected SPH and DIDF approach without particle shifting is
487
used. Except for some small deviations we get the same particle distribution
AC
CE
PT
ED
M
478
488
as in Fig. 3a. We note that the velocity distribution as well as the decay of the
489
kinetic energy are in good agreement with the analytic solution (not shown
490
here), but at later times the solution of the PPE tends to diverge because
491
void spaces between particle get too large, e.g. larger than the cutoff radius 29
ACCEPTED MANUSCRIPT
of the kernel function, and 2D-connectivity of the interpolation is lost. If
493
we apply the particle shifting approach alone or in addition with any of the
494
other stability approaches (Fig. 3e-h) the particle distribution as well as the
495
velocity distribution is very smooth. As shown by many authors [44, 47, 56]
496
the accuracy is largely improved by the particle shifting approach and leads
497
to better particle consistency. The effect of corrected SPH and DIDF is minor
498
in this test case. It is known from literature that the corrected SPH approach
499
improves the calculation of surface tension forces and DIDF improves volume
500
conservation [57, 35]. Since both approaches do not influence the particle
501
consistency, in this work, we use the corrected SPH and particle shifting
502
approach where particle movement or surface tension is involved. The DIDF
503
approach is only used in simulations where volume conservation is crucial,
504
since this is the main advantage of the approach. Its use is highlighted
505
explicitly in the text.
506
6.2. Lid-driven cavity
507
Circular flow in a cavity was theoretically and numerically studied by many
508
researchers e.g. [44, 50, 58, 59] to mention a few of them. An isotherm flow
509
inside a wall bounded cavity, as schematically shown in Fig. 4, driven by a
510
constant velocity on the top of the cavity is considered.
511
Due to the no-slip condition on the top wall, the fluid is accelerated and a
512
circular flow is initiated. In the following we perform simulations for different
513
Reynolds numbers and compare the axial velocity profiles with results from
514
the well validated grid-based solver OpenFOAM [60] version 2.4.0. We in-
515
vestigate different resolutions to show numerical convergence of the solution
516
of SPH.
AC
CE
PT
ED
M
AN US
CR IP T
492
30
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 3: Taylor-Green vortex simulation with different numerical stability approaches, namely (a) without numerical corrections, (b) with corrected SPH approach, (c) with DIDF approach, (d) with DIDF approach and corrected SPH approach, (e) with Particle
M
Shifting approach, (f) with corrected SPH approach and Particle Shifting approach, (g) with DIDF approach and Particle Shifting approach and (h) with corrected SPH approach,
ED
DIDF approach and Particle Shifting approach. Particle distribution shown at time t = 0.24s. The Reynolds number is Re = 1000. Color indicates the magnitude of the velocity of the particles. Red indicates |u| = 1m/s and blue |u| = 0m/s. Colors are distributed by
AC
CE
PT
an unsaturated rainbow color scheme.
Figure 4: Scheme of lid-driven cavity example.
31
ACCEPTED MANUSCRIPT
517
The Reynolds number for the lid-driven cavity is defined as Re =
ρLuw , η
(66)
where L is the length of the cavity, uw is the absolute velocity of the top
519
wall and ρ and η are the density and the dynamic viscosity of the fluid.
520
We perform simulations for Re = 10, Re = 100 and Re = 1000 where we
521
alter the Reynolds number by changing the wall velocity only. The size of the
522
cavity is quadratic and the length is L = 1mm. The density and viscosity are
523
ρ = 1000kg/m3 and η = 0.01P as, respectively. For each Reynolds number
524
we vary the resolution from L0 = 16.67µm to L0 = 4.17µm that corresponds
525
to 60 and 240 particles across L. For this test case we use the DIDF and
526
particle shifting approaches to avoid void spaces in the stagnation points in
527
the upper corners.
528
The boundary conditions for the velocity at the impenetrable walls are no-slip
529
boundary conditions. From physical point of view the boundary conditions
530
for pressure at the wall are Neumann boundary conditions. As a result the
531
reference pressure is not set in the linear equation system (LES). Therefore,
532
infinite solutions exist and an additional bootstrap is necessary. There are
533
two valid options to specify it. One option is to use a constant null space of
534
the krylov subspace solver. If we use this option we cannot use the algebraic
535
multigrid method as preconditioner in the PETSc library. Instead we may use
536
a Jacobi preconditioner. Therefore, the stability of the solution for disturbed
537
particle configuration decreases because the condition of the matrix for the
538
krylov subspace method gets worse. The number of iterations also increases
539
by a factor of 100.
540
Another option is to specify a reference pressure, e.g. p = 0, at a point in the
AC
CE
PT
ED
M
AN US
CR IP T
518
32
ACCEPTED MANUSCRIPT
domain where the pressure gradient is very low or zero. For the lid-driven
542
cavity such a point may be at the bottom corners. A result of that is a
543
spurious pressure in this corner because the specified pressure may not be
544
the best choice. Since the velocity in the corner is nearly zero the effect of
545
this spurious pressure vanishes but the stability of the solution is much better
546
because the algebraic multigrid preconditioner can be used. For that reason
547
we choose the latter option. We show the pressure field later in the discussion.
548
Note that instead of choosing a corner we may choose the bottom wall and
549
apply a constant pressure on the wall. The difference, between specifying the
550
pressure in the corner or the bottom wall, on the pressure is small. Therefore,
551
we choose the bottom left corner here to minimize numerical artifacts with
552
a reference pressure p = 0P a.
553
We use the grid-based solver OpenFOAM [60] version 2.4.0 to prepare a
554
reference solution with two different meshes, a 192×192 regular mesh for
555
Re = 10 and Re = 100 and a 240×240 regular mesh for Re = 1000. The
556
mesh size was chosen after investigation of convergence of the solution. For
557
finer meshes the solution was unchanged for the selected Reynolds number
558
(not shown here). The dimensionless velocity profiles u∗y =
559
at dimensionless position x∗ =
560
For Re = 10 (Fig. 5) the velocity along the vertical and horizontal axis is
561
almost independent of the resolution. We find good agreement between SPH
x L
and y ∗ =
y L
uy uw
and Ux∗ =
ux uw
are shown in Figs. 5, 6 and 7.
CE
PT
ED
M
AN US
CR IP T
541
and OpenFOAM for the horizontal velocity (Fig. 5b), however, small devia-
563
tions at the extreme value for the vertical velocity (Fig. 5a). One indication
564
for convergence is the position of the center of the vortex. The relative error
565
between the position of the vortex for the highest resolution in SPH and the
AC 562
33
CR IP T
ACCEPTED MANUSCRIPT
OpenFOAM (192x192)
0.15
SPH (60x60)
0.1
SPH (90x90)
0.05
SPH (120x120)
0
SPH (240x240)
-0.05
AN US
dimensionless velocity uy* [-]
0.2
-0.1 -0.15 -0.2 0
0.2
0.4
0.6
0.8
1
dimensionless position x* [-]
0.6 0.4 0.2
CE
0 -0.4
M OpenFOAM (192x192)
ED
0.8
PT
dimensionless position y* [-]
1
-0.2
SPH (60x60) SPH (90x90) SPH (120x120) SPH (240x240)
0
0.2
0.4
0.6
0.8
1
dimensionless velocity ux* [-]
AC
Figure 5: Velocity profiles for Re = 10 of a lid–driven cavity: (a) vertical velocity vs. x∗ at y ∗ = 0.5, and (b) horizontal velocity vs. y ∗ at x∗ = 0.5.
34
CR IP T
ACCEPTED MANUSCRIPT
0.2 0.1 0.05 0
OpenFOAM (192x192)
-0.05 -0.1
SPH (60x60)
-0.15
SPH (90x90)
-0.2
SPH (120x120)
-0.25
SPH (240x240)
-0.3 0
0.2
AN US
dimensionless velocity uy* [-]
0.15
0.4
0.6
0.8
1
dimensionless position x* [-]
0.6 0.4 0.2
CE
0 -0.4
M OpenFOAM (192x192)
ED
0.8
PT
dimensionless position y* [-]
1
-0.2
SPH (60x60) SPH (90x90) SPH (120x120) SPH (240x240)
0
0.2
0.4
0.6
0.8
1
dimensionless velocity ux* [-]
AC
Figure 6: Velocity profiles for Re = 100 of a lid–driven cavity: (a) vertical velocity vs. x∗ at y ∗ = 0.5, and (b) horizontal velocity vs. y ∗ at x∗ = 0.5.
35
CR IP T
ACCEPTED MANUSCRIPT
0.4 0.2 0.1 0
OpenFOAM (240x240)
-0.1 -0.2
SPH (60x60)
-0.3
SPH (90x90)
-0.4
SPH (120x120)
-0.5
SPH (240x240)
-0.6 0
0.2
AN US
dimensionless velocity uy* [-]
0.3
0.4
0.6
0.8
1
dimensionless position x* [-]
0.6 0.4 0.2
CE
0 -0.4
M OpenFOAM (240x240)
ED
0.8
PT
dimensionless position y* [-]
1
-0.2
SPH (60x60) SPH (90x90) SPH (120x120) SPH (240x240)
0
0.2
0.4
0.6
0.8
1
dimensionless velocity ux* [-]
AC
Figure 7: Velocity profiles for Re = 1000 of a lid–driven cavity: (a) vertical velocity vs. x∗ at y ∗ = 0.5, and (b) horizontal velocity vs. y ∗ at x∗ = 0.5.
36
ACCEPTED MANUSCRIPT
567
reference is u∗x = 0.60%, u∗y = 0.42% with respect to the reference solution.
568
tween different resolutions in SPH and the reference solution. The rela-
For Re = 100 (Fig. 6) we find good agreement in the velocity profiles be-
CR IP T
566
571
tive error in the position of the vortex for the highest resolution in SPH is u∗x = 0.24%, u∗y = 0.31% with respect to the reference solution.
572
change in the velocity profiles is sharper than before. Therefore, a good
573
agreement between SPH and reference solution is only obtained for the high-
574
est resolution. We see convergence of the vertical and horizontal velocity to
570
For Re = 1000 (Fig. 7) the effect of finer resolution is more obvious. The
AN US
569
577
the reference solution for higher resolution. The relative error in the position of the vortex for the highest resolution is u∗x = 0.60%, u∗y = 0.32% with
578
The 2D lid-driven cavity is investigated in the context of SPH in literature as
579
well [39, 44, 50, 58, 61, 62]. Compared to literature, the present results are
580
within a very good agreement with the reference solution and, therefore, we
581
conclude that the viscous and pressure force as well as the no-slip boundary
582
conditions for the velocity at the wall and numerical stability approaches
583
(DIDF and particle shifting) are well implemented.
584
At the beginning of the section we mentioned that an arbitrary bootstrap for
585
the solution of the Pressure Poisson Equation is needed. We note that the
586
velocity or pressure field is not affected by the artificial bootstrap as shown
ED
M
respect to the reference solution.
PT
576
CE
575
in Fig. 8. The reason is that the gradient of the pressure in the bottom left
588
corner is very small (Fig. 8b). Even if the pressure near the corner is different
589
from the specified bootstrap we obtain a good solution. This highlights that
590
the bootstrap is only necessary for algorithmic reasons during the solution
AC 587
37
AC
CE
PT
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 8: (a) Velocity and (b) pressure fields for Re = 100 with a resolution of L0 = 8.33µm.
38
591
of the linear equation system.
592
6.3. Diffusive transport of a scalar property
593
Diffusive transport of a scalar property φ
CR IP T
ACCEPTED MANUSCRIPT
Dφ = ∇ (α∇φ) , Dt
(67)
with the transport coefficient α appears in the energy equation (Eq. 9). We
595
validate diffusive transport for a general scalar property φ representing the
596
temperature. The discrete form of Eq. 67 was introduced in Eq. 19 and in
597
terms of the scalar variable φ is (∇ · α∇φ)a =
X mb b
ρb
AN US
594
(αa + αb )
~rab ∇a Wab · (φa − φb ) . 2 rab
(68)
In literature this discrete form for diffusive transport of a scalar variable
599
was validated, e.g. in [63, 64, 65]. Here we demonstrate convergence of
600
the present model. A very common test case is diffusion of a solute in one
601
direction. This is similar to heat conduction in a bar. For this test case a
ED
theoretical solution is available [66] 1 x0 − x x0 + x √ √ φ = φ0 erf + erf , 2 2 αt 2 αt
(69)
CE
PT
602
M
598
where x0 defines the region of concentrated solute φ0 at time t = 0. We
604
choose a transport coefficient α = 10−10 m2 /s, ρ = 1000kg/m3 , x0 = 0.5mm
AC
603
605
and
φ0 =
1 x < x0 0 else 39
,
(70)
ACCEPTED MANUSCRIPT
where the length of the domain is xmax = 2.5mm.
607
We show numerical convergence of the present model by increasing the reso-
608
lution for a smoothing length h = 2.05L0 . We consider 3 different resolutions
609
with L0 = 167µm, L0 = 83µm and L0 = 41.7µm. This corresponds to 15,
610
30 and 60 particles along x. We apply Neumann boundary condition at the
611
boundaries.
612
The profile of φ at time t = 1500s is shown in Fig. 9a. For the lowest resolu-
613
tion we see deviations from the theoretical profile especially near x = 0. The
614
profile for the highest resolution matches the theoretical profile very well.
615
The profile of φ for L0 = 41.7µm at different times in shown in Fig. 9b. For
616
all times the agreement is very good. Therefore, we conclude that the im-
617
plementation and the discrete form of the second order derivative are valid,
618
and therefore the heat transfer as well as mass transport in the bulk is valid.
M
AN US
CR IP T
606
ED
619
6.4. Buoyancy-driven cavity
621
Next we investigate non-isotherm single-phase flow to validate the coupling
622
of momentum and energy equation. Buoyancy-driven flow was investigated
623
in the context of SPH in a few articles [67, 68, 69]. Here we investigate a
624
buoyancy-driven cavity. In contrast to a lid-driven cavity, where the fluid is
625
driven by a velocity of the wall, the fluid is accelerated due to a temperature
626
difference between two walls. A schematic sketch is shown in Fig. 10.
627
The left wall is heated, the right wall is cooled and the top and bottom walls
628
are adiabatic. Due to the change in temperature, the fluid density decreases
629
with increasing its temperature. Because of a gravitational force, a circular
630
flow is initiated where denser fluid is accelerated downwards and lighter fluid
AC
CE
PT
620
40
ACCEPTED MANUSCRIPT
CR IP T
0.7
theoretical
φ [-]
0.6
L0 = 167µm
0.5
L0 = 83µm
0.4
L0 = 41.7 µm
0.3
AN US
0.2 0.1 0 0
0.5
1
1.5
2
2.5
position x [mm]
1
theoretical
0.8 0.7
t = 400s t = 1500s
ED
0.6 φ [-]
t = 100s
M
0.9
0.5 0.4 0.3
PT
0.2
t = 3000s
0.1
0
CE
0
0.5
1
1.5
2
2.5
position x [mm]
AC
Figure 9: Profile of φ along the x direction. (a) The profile for different resolutions at time t = 1500s and the theoretical profile are shown. (b) The profile of φ at different times with L0 = 41.7µm and the theoretical profiles are shown.
41
CR IP T
ACCEPTED MANUSCRIPT
Figure 10: Scheme for buoyancy-driven cavity example.
moves upwards. For small density variations the Boussinesq approximation
632
[70] is valid and the system can be regarded as incompressible since the pres-
633
sure is still independent of the density [20, 67]. Therefore, the only coupling
634
term between momentum balance and energy balance is the gravitational
635
acceleration ~g . According to Landau and Lifshitz [20], the difference in the
636
acceleration between two fluid portions with different temperature is
M
AN US
631
(71)
with the thermal expansion coefficient
β=−
1 ∂ρ , ρ ∂T p
(72)
CE
PT
637
ED
F~body = −β~g ∆T, ρ
in the range of 10−4 [1/K] for fluid considered here (i.e. polymers), and the
639
temperature difference
AC
638
∆T = T − T0 ,
(73)
640
with a reference temperature T0 . Note that the hydro-static part of the pres-
641
sure is neglected here because it should be negligible per definition. There42
ACCEPTED MANUSCRIPT
fore, we directly use Eq. 71 as an effective external force in the momentum
643
balance Eq. 8.
644
For natural convection the dimensionless numbers are the Prandtl number
Pr = 645
ν , k
the Rayleigh number
646
and the Gay-Lussac number
AN US
gβL3 ∆T , Ra = νk
CR IP T
642
Ga = β∆T.
(74)
(75)
(76)
The characteristic length L is the size of the cavity. ν and k are the kine-
648
matic viscosity of the fluid and the thermal diffusivity, respectively. The
649
Gay-Lussac number Ga describes the level of density variations caused by
650
temperature [67]. For Ga → 0 the Boussinesq approximation holds [71]. In
651
our case the temperature difference is ∆T = Th − Tc with Th and Tc as the
652
temperature at the hot and cold walls, respectively (see Fig 10).
653
We consider a square cavity of length L under the influence of gravity ~g =
654
[0, −9.81]. In the direction of gravity the cavity is adiabatic. In the other
655
direction a constant temperature is applied at the boundary. On the left side
656
we set the hot temperature to Th = 274.15K and on the right side the cold
657
temperature to Tc = 273.15K. We apply no-slip boundary conditions for the
658
velocity at all boundaries. Initially the temperature in the system is homoge-
659
neous T0 = 273.15K. The dynamic viscosity of the fluid is µ = 0.01P as and
660
kg m the density is ρ = 10 m 3 . The thermal diffusion coefficient is k = 0.00141 s
AC
CE
PT
ED
M
647
2
43
ACCEPTED MANUSCRIPT
and the thermal expansion coefficient is β = 0.071.
662
We investigate the flow in the cavity for two different Rayleigh numbers.
663
The size of the cavity varies with the Rayleigh number and is L = 126.4mm
664
and L = 272.4mm for Ra = 103 and Ra = 104 , respectively. We vary the
665
resolution to show convergence.
666
We compare the profile of the velocity and temperature along the vertical
667
and horizontal axis of the cavity with reference simulations using OpenFOAM
668
[60]. After investigation of the different resolutions in OpenFOAM we choose
669
a 120×120 regular grid for all reference simulations.
670
As already discussed in the previous section we need to set a bootstrap for
671
the pressure to solve the LES. In this case it is more difficult to identify a
672
point with very small gradient in the pressure. Therefore, it was more stable
673
to fix the pressure along the bottom wall at p = 0P a. We used the DIDF
674
and particle shifting approach as for the lid-driven cavity. The dimensionless
675
velocities are u∗y =
676
and y ∗ =
677
and Tmax = Th .
678
First we investigate a buoyancy-driven cavity at Ra = 1000. We vary the
679
resolution from L0 = 4.2mm to L0 = 1.05mm. The vertical and horizontal
680
velocity profiles along the axis at x∗ = 0.5 and y ∗ = 0.5 are shown in Fig. 11.
681
The velocity profiles are almost independent of the resolution and in good
ED
and u∗x =
ux ux,max
at dimensionless positions x∗ =
The dimensionless temperature is T ∗ =
T −T0 Tmax −T0
x L
with T0 = Tc
CE
PT
y . L
uy uy,max
M
AN US
CR IP T
661
agreement with the reference solution. The temperature profile along x∗ is
683
shown in Fig. 13a. Even with higher resolution we see a small deviation from
684
the reference solution. The deviation is in a reasonable range compared to
685
the available literature [67, 68, 69].
AC 682
44
ACCEPTED MANUSCRIPT
Next, we investigate the buoyancy-driven cavity at Ra = 10000. We again
687
vary the resolution. The vertical and horizontal velocity profiles along the
688
axis are shown in Fig. 12. For higher Rayleigh numbers the buoyancy force
689
is stronger. We see more deviations of the velocity in SPH compared to the
690
reference solution for coarse resolutions. With increasing resolution we find
691
convergence and good agreement to the reference solution. The temperature
692
profile shown in Fig. 13b also converges and is in very good agreement with
693
the reference solution.
694
In Fig. 14 we demonstrate for Ra = 10000 that the velocity and tempera-
695
ture field is smooth. We don’t see an effect of the artificial bootstrap for the
696
pressure and, therefore, skip the plot of the pressure here. The temperature
697
distribution is also in good agreement with literature [67]. Therefore, we
698
conclude that the coupling of the momentum and energy equation is valid.
699
Throughout this work we only need buoyancy flow with low and moderate
700
Rayleigh numbers. Therefore, we don’t show results for higher Rayleigh num-
701
bers here. The application of the SPH method to larger Rayleigh numbers
702
is shown in [67, 69].
703
7. Multi-phase systems
704
In this section, we validate the model for multi-phase applications. We dis-
705
cuss the formulation of surface tension including gradients of the surface
AC
CE
PT
ED
M
AN US
CR IP T
686
706
tension tangential to the interface. We assume that the surface tension is
707
temperature-dependent and analyze the surface force in a static and dy-
708
namic case to highlight features of the discretization. The normal part of the
709
surface tension is investigated using a spherical droplet. The tangential part 45
CR IP T
ACCEPTED MANUSCRIPT
OpenFOAM SPH (30x30) SPH (60x60) SPH (90x90) SPH (120x120)
0.5
0
AN US
dimensionless velocity uy* [-]
1
-0.5
-1 0
0.2
0.4
0.6
0.8
1
dimensionless position x* [-]
M
0.8
ED
0.6 0.4 0.2
PT
dimensionless position y* [-]
1
OpenFOAM SPH (30x30) SPH (60x60) SPH (90x90) SPH (120x120)
0
CE
-1
-0.5
0
0.5
1
dimensionless velocity ux* [-]
AC
Figure 11: Velocity profile for Ra = 1000: (a) vertical velocity vs. x∗ at y ∗ = 0.5, and (b) horizontal velocity vs. y ∗ at x∗ = 0.5.
46
CR IP T
ACCEPTED MANUSCRIPT
OpenFOAM SPH (30x30) SPH (60x60) SPH (90x90) SPH (120x120)
0.5
0
AN US
dimensionless velocity uy* [-]
1
-0.5
-1 0
0.2
0.4
0.6
0.8
1
dimensionless position x* [-]
M
0.8
ED
0.6 0.4 0.2
PT
dimensionless position y* [-]
1
OpenFOAM SPH (30x30) SPH (60x60) SPH (90x90) SPH (120x120)
0
CE
-1
-0.5
0
0.5
1
dimensionless velocity ux* [-]
AC
Figure 12: Velocity profile for Ra = 10000: (a) vertical velocity vs. x∗ at y ∗ = 0.5, and
(b) horizontal velocity vs. y ∗ at x∗ = 0.5.
47
CR IP T
ACCEPTED MANUSCRIPT
OpenFOAM SPH (30x30) SPH (60x60) SPH (90x90) SPH (120x120)
0.8 0.6 0.4
AN US
dimensionless temperature T* [-]
1
0.2 0 0
0.2
0.4
0.6
0.8
1
dimensionless position x* [-]
M
OpenFOAM SPH (30x30) SPH (60x60) SPH (90x90) SPH (120x120)
ED
0.8 0.6 0.4
PT
dimensionless temperature T* [-]
1
0.2
CE
0
0
0.2
0.4
0.6
0.8
1
dimensionless position x* [-]
AC
Figure 13: Temperature profile vs. x∗ at y ∗ = 0.5 for (a) Ra = 1000 and (b) Ra = 10000.
48
ACCEPTED MANUSCRIPT
of the surface tension is investigated using a planar interface.
711
Finally, we consider the complete model to investigate thermocapillary flow.
712
We consider the migration of a spherical droplet [72, 18, 17] and investi-
713
gate the stationary migration velocity. We show numerical convergence and
714
compare the SPH method to reference simulations taken from literature.
715
In a multi-phase system we get a non-isotropic pressure contribution, due to
716
the interface, in the momentum balance. This contribution is called capillary
717
stress in general or surface tension in an immiscible system [73]. The form of
718
the capillary stress tensor in the momentum balance accounts for a normal
719
and tangential stresses. Therefore, we may split the capillary stress into a
720
normal and tangential part. The tangential part is known as the Marangoni
721
force and accounts for gradients in the surface energy along the interface.
722
In this work we use the continuum surface force (CSF) model [31, 32] to
723
discretise the capillary stress. This model is well validated in literature for
724
SPH [12, 26, 32, 65, 74, 75]. In this section we focus on the validation of the
725
implementation. We use the corrected SPH approach and particle shifting
726
every time when capillary stress is considered.
727
As introduced in section 4.2 we use a color function to calculate the normal
AN US
M
ED
PT
at the interface between two phases. The color function is sharp 0 if fluid 1 Ca = , 1 if fluid 2
AC
CE
728
CR IP T
710
729
where each particle of a fluid has a constant color Ca .
49
(77)
AC
CE
PT
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 14: (a) Velocity and (b) temperature distributions in a buoyancy–driven cavity for Ra = 10000.
50
ACCEPTED MANUSCRIPT
7.1. Static capillary stress normal to interface
731
We first investigate the normal part of the capillary stress in a static test
732
case. We consider a spherical droplet in a second fluid at rest and compare
733
the pressure jump over the interface with the Young-Laplace equation [76].
734
In 2D the Young-Laplace equation for a pressure jump is ∆pY L =
σ , Rd
CR IP T
730
(78)
with the radius of the droplet Rd . We consider a droplet with a radius
736
Rd = 0.2m inside a quadratic box with size L = 1m filled with a second
737
fluid. Both fluids have the same physical properties. The density, dynamic
738
kg mN viscosity and surface tension are ρ = 1000 m 3 , η = 0.01P as and σ = 72.75 m ,
739
respectively. The resolution of the box is 120×120 particles where 24 particles
740
per droplet radius is used. Initially the particles are in regular Cartesian
741
order. At the boundary we apply no-slip conditions for the velocity. At the
742
top of the box we use Dirichlet boundary conditions and at the other three
743
boundaries Neumann boundary conditions for the pressure. The theoretical
744
pressure jump in the present case is ∆p = 0.36375P a.
745
At the beginning the order of the particles of the droplet isn’t smooth because
746
of the Cartesian order. We wait until an equilibrium configuration of the
747
particles is reached and then investigate the pressure jump. Fig. 15a shows a
748
3-dimensional graph of the pressure distribution in the box at an equilibrium
AC
CE
PT
ED
M
AN US
735
749
configuration.
750
We see that the pressure jump isn’t smooth. At the interface we observe
751
fluctuations where the pressure is higher or lower than the pressure in the
752
bulk phases. This is more obvious in the pressure profile along the radius 51
AN US
CR IP T
ACCEPTED MANUSCRIPT
0.4 0.2
M
Young-Laplace
SPH (120x120)
ED
0.6
PT
pressure p [Pa]
0.8
CE
0
-0.4
-0.2
0
0.2
0.4
radius r [m]
AC
Figure 15: (a) 3D graph and (b) plot of the pressure distribution of a single droplet in rest in a square box with surface tension σ = 72.75 mN m .
52
ACCEPTED MANUSCRIPT
of the droplet (Fig. 15b). These pressure fluctuations are also observed and
754
discussed in literature [65, 77]. The reason for the pressure fluctuations are
755
the spurious currents at the interface caused by unfavorable energetic con-
756
figuration of the particles [78]. The pressure in the bulk phase is constant.
757
Therefore, we quantify the pressure jump across the interface using the con-
758
stant pressure in the bulk phases. The pressure jump in the simulation is
759
∆pSP H = 0.3601P a. The relative error is ∆pY L − ∆pSP H · 100%, ∆pY L
AN US
=
CR IP T
753
(79)
of = 0.99%. This error is in the same order as in literature [32, 65].
761
7.2. Dynamic capillary stress normal to interface
762
Next we demonstrate that the dynamics of the fluid due to capillary stress
763
is valid. A common test case is an oscillating droplet, where an analytic
764
solution of the oscillation period is only available in the limit of a large
765
density ratio [79]. It was shown in literature that the CSF model using a
766
sharp color function accurately reproduces the oscillation period of a droplet
767
for different surface tensions [12, 32, 65, 77]. Numerical reference data for a
768
density ratio of 1 is also available in literature [65].
769
The test case consists of a spherical droplet with radius R = 0.2m in the
770
center of a square box with length L = 1m. At the boundaries we use
771
no-slip conditions for the velocity. The pressure is fixed at the top of the
772
box and Neumann boundary conditions for the pressure are applied at the
773
other boundaries. In both phases, the density, dynamic viscosity and surface
774
tension are ρ = 1kg/m3 , η = 0.05P as and σ = 1N/m, respectively. Initially,
AC
CE
PT
ED
M
760
53
ACCEPTED MANUSCRIPT
t [s]
[65]
≈ 0.4
sharp
0.409
CR IP T
color function
Table 1: Oscillation periods of droplet with density ratio 1. Comparison with estimated period from Adami et al. [65].
776
we impose a divergence-free velocity field [32, 65] x y2 r ux = U0 1− exp − , r0 r0 r r0 and y uy = −U0 r0
AN US
775
x2 r 1− exp − , r0 r r0
(80)
(81)
with u0 = 10m/s and r0 = 0.05m. The amplitude of the oscillation is
778
damped by the viscous force but for small viscosity the oscillation period is
779
not influenced.
780
The results are shown in Tab. 1. We estimate the period as the mean of x and
781
y position of the droplet, as shown in Fig. 7a. The mean is approximately
782
t = 0.4s. We are very close to the period of the reference literature. In
783
addition we compare the contour of the oscillating droplet at different times
784
with a contour from [65]. The particle distribution of the droplet is shown in
785
Fig. 16 where only the particles of the droplet are shown. The surrounding
786
particles are not shown. The particles in the upper half are taken from [65].
787
The particles in the bottom half correspond to the present simulation. The
788
contours match with small deviation for larger times, however, the particle
789
distributions are very similar even for longer periods.
790
Finally, we show the maximum position of the droplet in x and y direction
AC
CE
PT
ED
M
777
54
CR IP T
ACCEPTED MANUSCRIPT
Figure 16: Particle distribution of an oscillating droplet. Top half of every image is taken
from [65]. The bottom half of every image are present results. The surface tension is
AN US
N σ = 1m . From left to rigth, time is t = 0s, t = 08s, t = 16s, and t = 26s, respectively.
over time in Fig. 17, compared to results from [65]. There is a difference
792
between the amplitude of the present results and the literature. The reason
793
for the difference in the amplitude is that in the reference [65] the variation
794
of position is calculated in a different way than in the present simulations.
795
Adami et al. take the average of the mean position in a quarter of the
796
droplet and plot it over time. Here we take the maximum position of a
797
particle in x and y direction. Therefore, the variation is larger in the present
798
simulation but both kind of analysis indicate the frequency. We see that
799
the position oscillates around the average position x = 0.7m and y = 0.7m,
800
respectively. In comparison to the reference, the oscillation profiles in the
801
present simulations are very smooth. A very small drift of the position is
802
seen for t → 1s because the droplet moves a little bit out of the center of
803
the box. This is caused by small fluctuations of the pressure due to spurious
AC
CE
PT
ED
M
791
804
currents.
805
7.3. Static capillary stress tangential to interface
806
In this section we investigate the Marangoni force in a static test case. We
807
consider two layered fluids in an infinitely long channel of width 5.76mm in 55
ACCEPTED MANUSCRIPT
0.74
position [m]
0.72 0.71 0.7 0.69 0.68 0.67 0.66 0.2
0.4
0.6
0.8
1
AN US
0
CR IP T
X position (SPH) Y position (SPH) X position (Adami et al.) Y position (Adami et al.)
0.73
time t [s]
N Figure 17: Position variations over time for a droplet with R = 0.2m and σ = 1 m .
2D. The interface between the fluids is planar. We apply a linear temperature
809
profile along the channel with a slope of ∇T = 200K/m. We assume that
810
the surface tension coefficient depends linear on temperature
M
808
811
ED
σ = σ0 + σT (T − T0 ) ,
(82)
with the surface tension coefficient σ0 at T0 and the coefficient σT = −2mN/(mK). With these parameters the theoretical Marangoni force is ∇S σ = 0.4N/m2 .
813
The only force in the system is the thermocapillary force tangential to the
814
planar interface. We investigate the Marangoni force for different resolutions
815
to show numerical convergence. We investigate the resolutions L0 = 180µm,
816
L0 = 90µm and L0 = 45µm.
817
Because SPH is a smoothed discretization method the Marangoni force is
818
smoothed out around the interface between the fluids. In Fig. 18a we see the
819
profile of the Marangoni force perpendicular to the interface in the channel.
820
We see that the Marangoni force is wider for low resolution and narrower for
AC
CE
PT
812
56
ACCEPTED MANUSCRIPT
CR IP T
3
L0 = 180 µm
L0 = 90 µm L0 = 45 µm
2 1.5 1 0.5 0 2
2.2
2.4
2.6
AN US
force [kN/m3]
2.5
2.8
3
3.2
3.4
3.6
position [mm]
0.401
0.2
0.4
PT
0.3995
analytic
0.1
0.05
0.399
CE
0
50
0.15
error
ED
force [N]
0.4005
error [%]
M
SPH
100
150
0 200
L0 [µm]
AC
Figure 18: (a) Profile of Marangoni force in normal direction to the interface using a resolution of L0 = 180µm, L0 = 90µm and L0 = 45µm. (b) Integrated force and relative
error of the Marangoni force for different resolutions.
57
ACCEPTED MANUSCRIPT
higher resolution. In the limit of infinite resolution the profile tends to a Dirac
822
function with the magnitude of the Marangoni force. The calculated integral
823
of the smoothed Marangoni force of the simulation is shown in Fig. 18b.
824
The Marangoni force in SPH is only slightly different from the theoretical
825
value. The relative error is = 0.014% based on the theoretical value. The
826
calculated Marangoni force is independent of the resolution because we used
827
the corrected SPH approach with regular order of the particles and, therefore,
828
the force is as accurate as possible. In [17] a similar case was presented
829
where the relative error was also constant, but with = 3.65% very large.
830
The reason may be a different application of the corrected SPH approach in
831
[17]. One major difference is that we calculate the shepard kernel using only
832
particles with a normal larger than 0.01/h. We conclude that the Marangoni
833
force is very accurate in the present model.
834
7.4. Thermocapillary droplet migration
835
Thermocapillary motion was studied for several decades theoretically [76, 80],
836
numerically [18, 72] and experimentally [81, 82]. An overview about the work
837
of the 20th century is found in [83] and the references therein. In this section
838
we combine the normal and tangential surface force and investigate thermo-
839
capillary migration of a droplet in 2D and 3D. In thermocapillary flow, a
840
droplet moves because of gradients of the surface tension tangential to the
841
interface, caused by a gradient of the temperature, even in the absence of
842
gravity. For sufficiently low Reynolds numbers, the shape of the droplet is
843
spherical and constant during motion. Here, we first show numerical conver-
844
gence with increasing resolution, investigate the effect of a bounded domain
845
and compare our results to [17] and [18]. Young et al. [76] derived an ana-
AC
CE
PT
ED
M
AN US
CR IP T
821
58
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 19: Schematic setup of thermocapillary migration of a droplet.
lytic solution for the steady-state migration velocity of a droplet in 3D. We
847
compare the present results for 3D migration velocity to the analytic solution.
848
Initially, a spherical droplet (phase 1) with the radius of R is located in
849
the center of a square box with the length of L. In x direction we ap-
850
ply Neumann boundary conditions for the temperature and free slip condi-
851
tions for the velocity. In y direction we apply Dirichlet boundary conditions
852
for temperature and no-slip conditions for the velocity. The density, vis-
853
cosity and thermal diffusion coefficient of the droplet are ρ1 = 250kg/m3 ,
854
η1 = 0.012P as and λ1 = 1.2 · 10−6 W/(mK). The density, viscosity and
ED
PT
CE
thermal diffusion coefficient of the continuous phase are ρ2 = 500kg/m3 ,
AC
855
M
846
856
η2 = 0.024P as and λ2 = 2.4 · 10−6 W/(mK), respectively. The specific heat
857
capacities are cp,1 = 0.5 · 10−4 J/(kgK) and cp,2 = 10−4 J/(kgK). The surface
858
tension is σ = σ0 + σT (T − T0 ) , 59
(83)
ACCEPTED MANUSCRIPT
with σ0 = 10mN/m and σT = 2mN/(mK). The temperature at the bottom
860
of the box is the reference temperature T1 = T0 = 290K. The temperature
861
gradient in y direction is ∇T = 200K/m. The radius of the droplet is
862
R = 1.44mm. The setup of the test case is shown in Fig. 19. Note that
863
we don’t need a special treatment of the interface, e.g. we don’t use an
864
additional pressure contribution at the interface.
865
The dimensionless numbers to characterize thermocapillary migration are the
866
Reynolds number
AN US
ρ2 RUR = 0.72, η2
Re = 867
Marangoni number Ma =
ρ2 cp,2 RUR = 0.72, λ2
and Capillary number
M
868
Ca =
(85)
(86)
PT
σT |∇T |R m = 0.024 . η2 s
(87)
t · UR , R
(88)
U . UR
(89)
CE
The dimensionless time is t∗ =
and the dimensionless velocity is
AC
871
(84)
with the characteristic velocity of UR =
870
η2 UR = 0.0576, σ0
ED
869
CR IP T
859
U∗ =
872
First, we investigate numerical convergence of the migration velocity in 2D.
873
We choose a box length L = 5.76mm which corresponds to a ratio L/R = 4.
874
We vary the resolution from L0 = 180µm to L0 = 45µm. The dimensionless 60
ACCEPTED MANUSCRIPT
0.16 0.14 0.12
L0 = 180 µm L0 = 90 µm
0.08
CR IP T
U* [-]
0.1
L0 = 45 µm
0.06
L0 = 22.5 µm
0.04
Ma et al.
0.02
Tong et al.
0 0
0.5
1
1.5
2
AN US
t* [-]
Figure 20: Droplet migration using different resolutions compared to solution from [18] and [17]. The box length to radius ratio is L/R = 4.
migration velocity U ∗ over dimensionless time t∗ is shown in Fig. 20. Repre-
876
sentative plots of the acceleration due to surface forces, velocity, temperature
877
and pressure are shown in Figs. 21 and 22 for L0 = 45µm. We compare the
878
migration velocity to the results of Tong et al. [17] and Ma et al. [18] where
879
a resolution similar to L0 = 90µm is used.
880
In Fig. 20 we find two major results. Firstly, we see that the migration
881
velocity of the present model is very similar to the migration velocity of Ma
882
& Bothe [18]. Compared to Tong et al. [17] the present results are more
883
accurate. The reason is that the error in the Marangoni force is lower in the
884
present model. Secondly, we see that with increasing or decreasing resolution
885
we over- and underestimate the data of Ma & Bothe. The present model
AC
CE
PT
ED
M
875
886
converges with increasing resolution indicating that the results presented by
887
Tong et al. [17] were not yet converged.
888
If analytic solutions of migration velocity are available, the box size is always
889
assumed infinitely large. This means that the wall effects do not disturb the 61
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 21: Plot of droplet migration using L0 = 45µm particles at t = 0.6s. The box length to radius ratio is L/R = 4. On the left side the acceleration due to surface forces
AC
CE
PT
ED
M
is plotted and on the right side the magnitude of the velocity of shown.
Figure 22: Plot of droplet migration using L0 = 45µm particles at t = 0.6s. The box length to radius ratio is L/R = 4. On the left side the temperature is plotted and on the
right side the pressure is shown.
62
ACCEPTED MANUSCRIPT
0.2 0.18 0.16 0.12 0.1
CR IP T
U* [-]
0.14
L/R = 4
0.08
L/R = 8
0.06 0.04
L/R = 16
0.02
L/R = 32
0 0
0.5
1
1.5
2
AN US
t* [-]
Figure 23: Influence of box size on the migration velocity of a droplet in 2D. The droplet radius is resolved with 16 particles.
flow field. In numerical simulations we always have a bounded domain due to
891
the limited computational resources. Next we investigate the influence of the
892
bounded domain on the migration velocity of the droplet in 2D. We consider
893
the same setup and properties of the fluids as before. Then we increase the
894
box size from L = 5.76mm to L = 46.08mm. The resolution is L0 = 90µm
895
or 16 particles per droplet radius. Initially the droplet is placed in the center
896
of the box.
897
The migration velocity U ∗ over time t∗ is shown in Fig. 23. We see that
898
with increasing box size the steady-state migration velocity increases and
899
converges to a limit of U ∗ ≈ 0.18. For L/R ratios larger than 16 the migra-
900
tion velocity is converged and we may assume that the flow field is approx-
AC
CE
PT
ED
M
890
901
imately unaffected by the boundaries. A similar value but for the 3D case
902
was presented by Ma & Bothe [18] where the steady-state migration velocity
903
is constant for L/R ratios larger than 16.
904
Finally, we consider thermocapillary migration of a spherical droplet in 3D. 63
ACCEPTED MANUSCRIPT
905
For the steady-state migration velocity of a spherical droplet in 3D, Young
906
et al. [76] derived a theoretical value of 2 , η1 λ1 2 + 3 λ2 η2
(90)
CR IP T
U∗ = 2+
with the ratio of the thermal conductivity and dynamic viscosity of the
908
droplet and the continuous phase, where they assumed that the domain is
909
unbounded. But in the simulation we have a bounded domain and therefore
910
we will reach the theoretical value only in the limit of a large L/R ratio.
911
In the 3D model we consider a spherical droplet in a closed box. We use slip
912
boundary condition for velocity in x and z directions and no-slip boundary
913
conditions in y direction. Similarly we use Neumann boundary condition for
914
the temperature in x and z directions and Dirichlet boundary conditions in
915
y direction. The properties of the fluids are the same as in the previous sim-
916
ulations. Initially a spherical droplet is centered in the box. To reduce the
917
numerical effort we use 16 particles per droplet radius with L0 = 90µm. The
918
smoothing length is reduced to h = 1.55L0 because there are more particles
919
in the vicinity in 3D than in 2D and, therefore, the SPH approximation is
920
already good for lower smoothing lengths.
921
A representative snapshot of the 3D simulation is shown in Fig. 24. Only
922
the droplet without the surrounding fluid is shown. The temperature and
923
velocity distributions at the surface of the droplet are indicated by the color
AC
CE
PT
ED
M
AN US
907
924
and the arrows. The velocity at the bottom is directed into the droplet while
925
the velocity at the top is directed out of the droplet. The temperature profile
926
is close to linear.
927
In Fig. 25 we demonstrate that we converge towards the theoretical value 64
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 24: Plot of droplet migration using L0 = 90µm particles at t = 0.6s. The box length to radius ratio is L/R = 4. The color indicates temperature and the arrows indicate the
M
direction of velocity.
by increasing the L/R ratio. We see that the steady-state migration velocity
929
increases with increasing L/R ratio. We underestimate the theoretical value
930
because the domain is still bounded but we are very close for L/R = 8. As
931
we have seen at the beginning of this section the influence of the bounded
932
domain is only negligible for L/R ratio larger than 16.
933
We also compare the results of the present simulations with results of Ma
934
et al. [18]. The present results are slightly closer to the theoretical value
935
than the results of Ma & Bothe. We have not performed simulations for very
936
large L/R ratios because the numerical effort is very high. In the results of
937
Ma & Bothe it seems that there is a constant offset in the migration velocity
938
because the migration velocity is very similar for L/R = 8 and L/R = 16.
939
We expect similar behavior of the present model because we use a very sim-
940
ilar formulation of the CSF model. Therefore, we conclude that the present
AC
CE
PT
ED
928
65
ACCEPTED MANUSCRIPT
0.25 0.2
0.1
CR IP T
U* [-]
0.15
SPH
0.05
Ma et al. Young et al.
0 0
5
10
20
AN US
L/R ratio [-]
15
Figure 25: Influence of box length to droplet diameter ratio for 3D simulation where the droplet diameter is resolved with 32 particles. Comparison with analytic migration velocity [76, 84]
.
model predicts thermocapillary flow very accurate compared to available lit-
942
erature.
943
8. Conclusion
944
We presented an incompressible SPH model for thermo-capillary flow driven
945
by gradients of the surface tension. Rigorous validation of the proposed
946
model shows accurate implementation and convergence of the results.
947
In the first part of the article, we validated the model for single-phase prob-
948
lems. First, we demonstrated the effects of corrected SPH, DIDF and par-
AC
CE
PT
ED
M
941
949
ticle shifting approach on particle distribution using a Taylor-Green vortex.
950
It is shown that only particle shifting enables a homogeneous particle dis-
951
tribution. Then, lid-driven cavity is studied for different Reynolds numbers.
952
The results are compared to reference solutions taken from OpenFOAM. A 66
ACCEPTED MANUSCRIPT
convergence study shows that with increasing resolution the SPH model con-
954
verges to the OpenFOAM results. Next, we validated diffusive mass trans-
955
port, compared it to analytic solution and found very accurate agreement.
956
Afterwards, buoyancy-driven flow is investigated for different Rayleigh num-
957
bers and compared to reference solution of OpenFOAM. Again, a convergence
958
study demonstrates that the results of SPH match the results of OpenFOAM
959
very well. Therefore, we conclude that the model was correctly implemented.
960
In the second part of the article, we investigated multi-phase systems. First,
961
we studied surface tension in detail using static and dynamic test cases.
962
The normal component of the surface tension is validated using the Young-
963
Laplace law at steady-state in 2D. The analytic law is captured very well.
964
Deviations of the pressure at the interface are a result of the CSF model and
965
where previously found in other models as well. The error is approx. 1 %.
966
We studied the oscillation of a 2D droplet and compared the results with
967
reference solution of Adami et al. from literature. We analyzed the time
968
evolution of the oscillation and its deviation to the reference. The oscillation
969
period as well as the contour of the 2D droplet match very well.
970
The Marangoni force is studied separately for a planar interface with a tan-
971
gential, linear temperature gradient. A convergence study is presented and
972
we found that the error is around 0.1 %, and independent of the numerical
973
resolution because of the use of corrected SPH approach for the Marangoni
CE
PT
ED
M
AN US
CR IP T
953
force. This is way better than in Tong et al.
975
Finally, we studied thermo-capillary droplet migration in a temperature field
976
with temperature-dependent surface tension coefficient. First, we present a
977
convergence study for 2D droplet and compare the results to reference of SPH
AC 974
67
ACCEPTED MANUSCRIPT
(Tong et al.) and Finite Volume (Ma et al.). We show that the previously
979
presented results of Tong et al. agree with FVM but have not converged in
980
SPH. Then we study the same case in 3D and compare the results to the an-
981
alytic solution of Young et al. The results of SPH show convergence and are
982
slightly more accurate than the results of FVM taken from Ma et al. There-
983
fore, we conclude that the proposed model for thermo-capillary flow captures
984
the studied cases very well. Compared to the previous work of Tong et al.,
985
this study completes with a convergence study and more accurate results.
986
9. Acknowledgement
987
M. Hopp-Hirschler and U. Nieken acknowledge funding from the German Re-
988
search Foundation within the Collaborative Research Center SFB716. M.S.
989
Shadloo acknowledges funding provided by DAAD programme ”Research
990
Stays for University Academics and Scientists, 2017” through funding ID:
991
57314019 during his visit from Stuttgart University. The access to the French
992
CRIANN (Centre Rgional Informatique et d’Applications Numriques de Nor-
993
mandie) resources, under allocation 2017002 is also acknowledged. This work
994
was granted access to HPC resources of IDRIS under the allocation 2017-
995
100752 made by GENCI (A0022A10103).
996
10. References
AC
CE
PT
ED
M
AN US
CR IP T
978
997
998
999
References [1] L. B. Lucy, A numerical approach to the testing of the fission hypothesis, The Astronomical Journal 82 (1977) 1013–1024.
68
ACCEPTED MANUSCRIPT
[2] R. A. Gingold, J. J. Monaghan, Smoothed Particle Hydrodynamics:
1001
theory and application to non-spherical stars, Mon. Not. R. astr. Soc
1002
181 (1977) 375–389.
1004
1005
1006
1007
1008
[3] S. J. Cummins, M. Rudman, An SPH Projection Method, J. Comp. Phys. 152 (1999) 584–607.
[4] J. J. Monaghan, Simulating free surface flows with sph, Journal of Computational Physics 110 (1994) 399–406.
AN US
1003
CR IP T
1000
[5] J. P. Morris, P. J. Fox, Y. Zhu, Modeling low reynolds number incompressible flows using sph, J. Comp. Phys. 136 (1997) 214–226. [6] A. Colagrossi, M. Landrini, Numerical simulation of interfacial flows by
1010
smoothed particle hydrodynamics, J. Comp. Phys. 191 (2003) 448–475.
1012
[7] J. J. Monaghan,
Smoothed Particle Hydrodynamics,
Reports on
Progress in Physics 68, Issue 8 (2005) 1703–1759.
ED
1011
M
1009
[8] S. Shao, E. Y. M. Lo, Incompressible SPH method for simulating New-
1014
tonian and non-Newtonian flows with a free surface, Advances in Water
1015
Resources 26 (2003) 787–800. [9] M. Shadloo, G. Oger, D. Le Touz´e, Smoothed Particle Hydrodynamics method for fluid flows, towards industrial applications: Motivations,
AC
1017
CE
1016
PT
1013
1018
current state, and challenges, Computers & Fluids 136 (2016) 11 – 34.
1019
[10] M. B. Liu, G. R. Liu, Smoothed Particle Hydrodynamics (SPH): an
1020
Overview and Recent Developments, Arch. Comput. Method. Eng. 17
1021
(2010) 25–76. 69
ACCEPTED MANUSCRIPT
[11] N. Tofighi, M. Yildiz, Numerical simulation of single droplet dynam-
1023
ics in three-phase flows using isph, Computers & Mathematics with
1024
Applications 66 (2013) 525–536.
CR IP T
1022
1025
[12] X. Y. Hu, N. A. Adams, A multi-phase sph method for macroscopic
1026
and mesoscopic flows, Journal of Computational Physics 213 (2006)
1027
844–861.
[13] D. Le Touz´e, G. Oger, B. Alessandrini, Smoothed Particle Hydrodynam-
1029
ics simulation of fast ship flows, in: Proc. 27th Symposium on Naval
1030
Hydrodynamics.
AN US
1028
[14] X. Y. Hu, N. A. Adams, A constant-density approach for incompressible
1032
multi-phase sph, Journal of Computational Physics 228 (2009) 2082–
1033
2091.
M
1031
[15] A. Zainali, N. Tofighi, M. S. Shadloo, M. Yildiz, Numerical investigation
1035
of newtonian and non-newtonian multiphase flows using isph method,
1036
Computer Methods in Applied Mechanics and Engineering 254 (2013)
1037
99–113.
1039
PT
[16] S. Adami, X. Y. Hu, N. A. Adams, A conservative sph method for surfac-
CE
1038
ED
1034
tant dynamics, Journal of Computational Physics 229 (2010) 1909–1926.
[17] M. Tong, D. Browne, An incompressible multi-phase smoothed particle
1041
hydrodynamics (sph) method for modelling thermocapillary flow, Int.
1042
J. Heat Mass Trans. 73 (2014) 284–292.
AC
1040
1043
[18] C. Ma, D. Bothe, Direct numerical simulation of thermocapillary flow 70
ACCEPTED MANUSCRIPT
1044
based on the volume of fluid method, International Journal of Multi-
1045
phase Flow 37 (2011) 1045 – 1058. [19] P. Espanol, C. Thieulot, Microscopic derivation of hydrodynamic equa-
1047
tions for phase-separating fluid mixtures,
1048
Physics 118 (2003).
CR IP T
1046
The Journal of Chemical
[20] L. D. Landau, E. M. Lifshitz, Course of Theoretical Physics: Fluid
1050
Dynamics, volume 6, Elsevier Ltd., Butterwoth Heinemann, 2nd edition,
1051
1987.
1053
1054
1055
[21] M. Hopp-Hirschler, Modeling of Porou Polymer Membrane Formation, Dissertation d93, shaker verlag, Universit¨at Stuttgart, 2017. [22] J. J. Monaghan, Smoothed Particle Hydrodynamics, Annu. Rev. Stron.
M
1052
AN US
1049
Astrphys. 30 (1992) 543–574.
[23] H. Wendland, Piecewise polynominal, positive definite and compactly
1057
supported radial functions of minimal degree, Adv. Comput. Math. 4
1058
(1995) 389–396.
1060
PT
Press, 2012.
[25] J. Bonet, T.-S. L. Lok, Variational and momentum preservation as-
AC
1061
[24] D. Violeau, Fluid Mechanics and the SPH Method, Oxford University
CE
1059
ED
1056
1062
pects of Smooth Particle Hydrodynamic formulations, Comput. Method.
1063
Appl. M. 180 (1999) 97–115.
1064
[26] N. Grenier, M. Antuono, A. Colagrossi, D. Le Touz´e, B. Alessandrini, An
71
ACCEPTED MANUSCRIPT
1065
Hamiltonian interface SPH formulation for multi-fluid and free surface
1066
flows, Journal of Computational Physics 228 (2009) 8380–8393.
1068
[27] L. Brookshaw, A method of calculating radiative heat diffusion in par-
CR IP T
1067
ticle simulations, Proc. ASA 6(2) (1985) 207–210.
[28] R. Fatehi, M. T. Manzari, Error estimation in smoothed particle hy-
1070
drodynamics and a new scheme for second derivatives, Computers and
1071
Mathematics with Applications 61 (2011) 482–498.
AN US
1069
1072
[29] K. Szewc, J. Pozorski, J.-P. Minier, Analysis of the incompressibility
1073
constraint in the smoothed particle hydrodynamics method, Int. J. Nu-
1074
mer. Meth. Eng. 92 (2012) 343–369.
[30] J. Morris, Analysis of smoothed particle hydrodynamics with applica-
1076
tions, Ph.D. thesis, Department of Mathematics, Monash University,
1077
1996.
1080
1081
ED
[32] J. P. Morris, Simulating surface tension with Smoothed Particle Hydrodynamics, Int. J. Numer. Meth. Fluids 33 (2000) 333–353.
[33] P. G. Tait, Report on some of the physical properties of fresh water and
AC
1082
Modeling Surface Tension, J. Comp. Phys. 100 (1992) 335–354.
PT
1079
[31] J. U. Brackbill, D. B. Kothe, C. Zemach, A Continuum Method for
CE
1078
M
1075
1083
sea water, Rept. Sci. Results Voy., H.M.S. Challanger, Phys. Chem. 2
1084
(1888) 1–76.
1085
[34] H. Helmholtz,
u ¨ber integrale der hydrodynamicschen gleichungen,
72
ACCEPTED MANUSCRIPT
1086
welcher der wirbelbewegungen entsprechen, Journal fr die reine und
1087
angewandte Mathematik 55 (1858) 25–55.
1089
[35] X. Y. Hu, N. A. Adams, An incompressible multi-phase sph method, J. Comp. Phys. 227 (2007) 264–278.
CR IP T
1088
[36] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschel-
1091
man, L. Dalcin, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley,
1092
L. C. McInnes, K. Rupp, B. F. Smith, S. Zampini, H. Zhang, H. Zhang,
1093
PETSc Web page, http://www.mcs.anl.gov/petsc, 2016. Version 3.3.
1094
[37] R. Falgout, A. Cleary, J. Jones, E. Chow, V. Henson, C. Baldwin,
1095
P. Brown, P. Vassilevski, U. M. Yang, Hypre Web page, http://acts.
1096
nersc.gov/hypre, 2016.
M
1098
[38] F. Keller, Simulation of the Morphogenesis of Open-porous Materials, Ph.D. thesis, University of Stuttgart, 2015.
ED
1097
AN US
1090
[39] M. Yildiz, R. A. Rook, A. Suleman, Sph with the multiple boundary
1100
tangent method, International Journal for Numerical Methods in Engi-
1101
neering 77 (2009) 1416–1438.
1103
Unified semi-analytical wall boundary conditions for inviscid laminar or turbulent flows in the meshless sph method, Int. J. Numer. Meth. Fluids
AC
1104
[40] M. Ferrand, D. R. Laurence, B. D. Rogers, D. Violeau, C. Kassiotis,
CE
1102
PT
1099
1105
71 (2013) 446–472.
1106
[41] S. Khorasanizade, J. M. M. Sousa, An innovative open boundary treat-
1107
ment for incompressible sph, International Journal for Numerical Meth-
1108
ods in Fluids 80 (2016) 161–180. 73
ACCEPTED MANUSCRIPT
[42] M. Hirschler, P. Kunz, M. Huber, F. Hahn, U. Nieken, Open boundary
1110
conditions for ISPH and their application to micro-flow, Journal of
1111
Computational Physics 307 (2016) 614 – 633.
CR IP T
1109
[43] M. S. Shadloo, A. Zainali, S. H. Sadek, M. Yildiz, Improved incom-
1113
pressible smoothed particle hydrodynamics method for simulating flow
1114
around bluff bodies, Computer methods in applied mechanics and engi-
1115
neering 200 (2011) 1008–1020.
AN US
1112
1116
[44] R. Xu, P. Stansby, D. Laurence, Accuracy and stability in incompressible
1117
SPH (ISPH) based on the projection method and a new approach, J.
1118
Comp. Phys. 228 (2009) 6703–6725.
[45] M. S. Shadloo, A. Zainali, M. Yildiz, Simulation of single mode rayleigh–
1120
taylor instability by sph method, Computational Mechanics 51 (2012)
1121
699–715.
ED
M
1119
[46] M. Shadloo, A. Rahmat, M. Yildiz, A smoothed particle hydrodynamics
1123
study on the electrohydrodynamic deformation of a droplet suspended
1124
in a neutrally buoyant newtonian fluid, Computational Mechanics 52
1125
(2013) 693–707.
CE
PT
1122
[47] S. J. Lind, R. Xu, P. K. Stansby, B. D. Rogers, Incompressible Smoothed
1127
Particle Hydrodynamics for free-surface flows: A generalised diffusion-
AC
1126
1128
based algorithm for stability and validations for impulsive flows and
1129
propagating waves, J. Comp. Phys. 231 (2012) 1499–1523.
1130
1131
[48] D. Shepard, A two dimensional function for irregulary space data, Proceedings of ACM national conference (1968) 517–524. 74
ACCEPTED MANUSCRIPT
[49] X. Hu, personal communication, SPHERIC 2015, Parma, 2015.
1133
[50] U. Ghia, K. N. Ghia, C. T. Shin, High–re solutions for incompressible
1134
flow using thre navier–stokes equations and a multigrid method, J.
1135
Comp. Phys. 48 (1982) 387–411.
1137
1138
1139
[51] R. Schreiber, H. Keller, Driven cavity flows by efficient numerical techniques, J. Comput. Phys. 49 (1983) 310 – 333.
[52] F. Penot, Numerical calculation of 2d natural convection in isothermal
AN US
1136
CR IP T
1132
open cavities, Numer. Heat Transfer 5 (1982) 421–437. [53] O. LeQuere, J. A. C. Humphery, F. S. Sherman, Numerical calculation
1141
of thermally driven 2d unsteady laminar flow in cavities of rectangular
1142
cross section, Numer. Heat Transfer 4 (1981) 249–283.
M
1140
[54] B. Hajra, A review of some recent studies on buoyoncy driven flows in
1144
an urban environment, International Journal of Atmospheric Sciences
1145
2014 (2014) 15, Article ID 362182.
PT
1147
[55] G. I. Taylor, A. E. Green, Mechanism of the production of small eddies from large ones, Proc. R. Soc. Lond. A 158 (1937) 499–521.
CE
1146
ED
1143
[56] M. S. Shadloo, A. Zainali, M. Yildiz, A. Suleman, A robust weakly
1149
compressible sph method and its comparison with an incompressible
AC
1148
1150
sph, International Journal for Numerical Methods in Engineering 89
1151
(2012) 939–956.
1152
[57] M. Huber, F. Keller, W. Sckel, M. Hirschler, P. Kunz, S. Hassanizadeh,
1153
U. Nieken, On the physically based modeling of surface tension and 75
ACCEPTED MANUSCRIPT
1154
moving contact lines with dynamic contact angles on the continuum
1155
scale, Journal of Computational Physics 310 (2016) 459 – 477. [58] E.-S. Lee, C. Moulinec, R. Xu, D. Violeau, D. Laurence, P. Stansby,
1157
Comparisons of weakly compressible and truly incompressible algo-
1158
rithms for the sph mesh free particle method, Journal of Computational
1159
Physics 227 (2008) 8417–8436.
1161
1162
1163
[59] A. Rafiee, Modelling of generalized Newtonian lid-driven Cavity Flow
AN US
1160
CR IP T
1156
using an SPH method, The ANZIAM Journal 49 (2008) 411–422. [60] OpenFOAM Documentation, The OpenFOAM Foundation, 2.4.0 edition, 2015.
[61] M. Basa, N. J. Quinlan, M. Lastiwka, Robustness and accuracy of
1165
sph formulations for viscous flow, International Journal for Numerical
1166
Methods in Fluids 60 (2009) 1127–1148.
ED
M
1164
[62] S. Khorasanizade, J. M. M. Sousa, A detailed study of lid-driven cavity
1168
flow at moderate reynolds numbers using incompressible sph, Interna-
1169
tional Journal for Numerical Methods in Fluids 76 (2014) 653–668.
1170
[63] P. W. Cleary, J. J. Monaghan, Conduction modelling using smoothed
1171
particle hydrodynamics, Journal of Computational Physics 148 (1999)
CE 227–264.
AC
1172
PT
1167
1173
1174
[64] Y. Zhu, P. J. Fox, Smoothed particle hydrodynamics model for diffusion through porous media, Transport in Porous Media 43 (2001) 441–471.
76
ACCEPTED MANUSCRIPT
[65] S. Adami, X. Y. Hu, N. A. Adams, A new surface-tension formulation for
1176
multi-phase sph using a reproducing devergence approximation, Journal
1177
of Computational Physics 229 (2010) 2011–5021.
1178
1179
CR IP T
1175
[66] J. Crank, The mathematics of duiffusion, Oxford University Press, 2nd edition, 1975.
[67] K. Szewc, J. Pozorski, A. Tanire, Modeling of natural convection with
1181
smoothed particle hydrodynamics: Non-boussinesq formulation, Inter-
1182
national Journal of Heat and Mass Transfer 54 (2011) 4807 – 4816.
AN US
1180
[68] A. G. V., B. Firoozabadi, M. Mahdinia, 2d numerical simulation of
1184
density currents using the {SPH} projection method, European Journal
1185
of Mechanics - B/Fluids 38 (2013) 38 – 46.
M
1183
[69] A. Leroy, D. Violeau, M. Ferrand, A. Joly, Buoyancy modelling with in-
1187
compressible sph for laminar and turbulent flows, International Journal
1188
for Numerical Methods in Fluids 78 (2015) 455–474.
ED
1186
[70] J. Boussinesq, Theorie de l’ecoulement tourbillonnant et tumultueux des
1190
liquides dans les lits rectilignes a grande section, volume 1, Gauthier-
1191
Villars et Fils, 1897.
CE
PT
1189
[71] T. Pesso, S. Piva, Laminar natural convection in a square cavity: low
1193
prandtl numbers and large density differences, Int. J. Heat Mass Transfer
AC
1192
1194
1195
1196
52 (2009) 1036–1043.
[72] S. Nas, G. Tryggvason, Thermocapillary interaction of two bubbles or drops, Int. J. Multi. Flow 29 (2003) 1117–1135. 77
ACCEPTED MANUSCRIPT
1198
1199
1200
[73] J. G. Kirkwood, F. P. Buff, The statistical mechanical theory of surface tension, The Journal of Chemical Physics 17 (1949) 338–343. [74] J. J. Monaghan, Smoothed particle hydrodynamics and its diverse ap-
CR IP T
1197
plications, Annu. Rev. Fluid Mech. 44 (2012) 323–346.
[75] M. S. Shadloo, M. Yildiz, Numerical modeling of kelvinhelmholtz in-
1202
stability using smoothed particle hydrodynamics, International Journal
1203
for Numerical Methods in Engineering 87 (2011) 988–1006.
1205
1206
1207
[76] N. Young, J. Goldstein, M. Block, The motion of bubbles in a vertical temperature gradient, J. Fluid Mech. 6 (1959) 350–356. [77] K. Szewc, J. Pozorski, J.-P. Minier, Spurious interface fragmentation in multiphase sph, Int. J. Numer. Meth. Engng 103(9) (2015) 625–649.
M
1204
AN US
1201
[78] S. Heß, V. Springel, Particle hydrodynamics with tessellation tech-
1209
niques, Monthly Notices of the Royal Astronomical Society 406 (2010)
1210
2289–2311.
1213
Royal Society of London 29 (1879) 71–97. [80] R. Balasubramaniam, A.-T. Chai,
Thermocapillary migration of
droplets: An exact solution for small marangoni numbers, Journal of
AC
1214
PT
1212
[79] L. Rayleigh, On the capillary phenomena of jets, Proceedings of the
CE
1211
ED
1208
1215
Colloid and Interface Science 119 (1987) 531 – 538.
1216
[81] G. Wozniak, R. Balasubramaniam, H. P. Hadland, S. R. Subramanian,
1217
Temperature fields in a liquid due to the thermocapillary motion of
1218
bubbles and drops, Experiments in Fluids 31 (2001) 84–89. 78
ACCEPTED MANUSCRIPT
[82] Q. Kang, L. Hu, C. Huang, H. Cui, L. Duan, W. Hu, Experimental
1220
investigations on interaction of two drops by thermocapillary-buoyancy
1221
migration, International Journal of Heat and Mass Transfer 49 (2006)
1222
2636 – 2641.
1224
1225
Drops in Reduced Gravity, Cambridge University Press, 2001. [84] A. Fedosov,
Thermocapillary motion,
arXiv:1303.0243.
AC
CE
PT
ED
M
1226
[83] R. S. Subramanian, R. Balasubramaniam, The Motion of Bubbles and
arXiv preprint (2013)
AN US
1223
CR IP T
1219
79
A Smoothed Particle Hydrodynamics Approach for Thermo-Capillary Flows Manuel Hopp-Hirschler, Mostafa Safdari Shadloo, Ulrich Nieken PII: DOI: Reference:
S0045-7930(18)30640-6 https://doi.org/10.1016/j.compfluid.2018.09.010 CAF 4005
To appear in:
Computers and Fluids
Received date: Revised date: Accepted date:
24 March 2018 4 September 2018 13 September 2018
Please cite this article as: Manuel Hopp-Hirschler, Mostafa Safdari Shadloo, Ulrich Nieken, A Smoothed Particle Hydrodynamics Approach for Thermo-Capillary Flows, Computers and Fluids (2018), doi: https://doi.org/10.1016/j.compfluid.2018.09.010
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights
flow driven by gradients of the surface tension.
CR IP T
• Smoothed Particle Hydrodynamics (SPH) model for thermo-capillary
• Continuum surface force (CSF) approach including Marangoni forces.
AC
CE
PT
ED
M
AN US
• Convergence study for thermos-capillary flow
1
ACCEPTED MANUSCRIPT
A Smoothed Particle Hydrodynamics Approach for Thermo-Capillary Flows
CR IP T
Manuel Hopp-Hirschler1a , Mostafa Safdari Shadloob , Ulrich Niekena a
AN US
Institute of Chemical Process Engineering, University of Stuttgart, 70199 Stuttgart, Germany b CORIA-UMR 6614, CNRS-University, INSA of Rouen and Normandie University, 76000 Rouen, France
Abstract
Interfacial-driven flows are important phenomena in many processes. In this article, we present a Smoothed Particle Hydrodynamics (SPH) model for thermo-capillary flow driven by gradients of the surface tension. The model
M
is based on the continuum surface force (CSF) approach including Marangoni forces. An incompressible SPH approach using (i) density-invariant divergence-
ED
free (DIDF), (ii) corrected SPH and (iii) particle shifting approaches for multi-phase systems is used for accurate results.
PT
We carefully validate the proposed model using several test cases. First, we demonstrate the effects of corrected SPH and particle shifting approaches
CE
using Taylor-Green vortex. Then, we study single-phase flow problems to validate correct implementation of boundary conditions, momentum and energy
AC
balance using lid-driven cavity, diffusive transport problem, and buoyancydriven cavity test cases. Afterward, we investigate different multi-phase flow problems to validate nor1
corresponding author: [email protected]
Preprint submitted to Computers & Fluids
September 13, 2018
ACCEPTED MANUSCRIPT
mal and tangential component of the surface tension. Finally, we apply the model to thermo-capillary rise of a droplet due to a temperature gradient.
CR IP T
We present a convergence study and compare the results with their counterparts obtained from OpenFoam software as well as Finite Volume method
(FVM) reference from literature. We demonstrate that the proposed model is very accurate for thermo-capillary flow. The simulation results of the current SPH approach will be available online for the community.
AN US
Keywords: Thermal Convection, Marangoni Effect, Multiphase Flow, Smoothed Particle Hydrodynamics (SPH), Meshless Method 1. Introduction
2
Smoothed particle hydrodynamics (SPH) is a meshless Lagrangian particle
3
method to solve partial differential equations. Tough, initially introduced for
4
astrophysical problems, by Lucy [1] and Monaghan and Gingold [2] in 1977,
5
SPH gained popularity in 1990s for fast dynamic flows [3, 4, 5], in 2000s for
6
multi-phase problems [6, 7, 8] and lately in diverse engineering applications
7
[9, 10].
8
The main reason for this popularity, especially for multi-phase flows, is its
9
inherent strengths which some of those may be found in other methods but
10
the combination of all those features can be found only in SPH [9]. Among
11
others one can mention, firstly, the natural distinguishing between different
AC
CE
PT
ED
M
1
12
phases due to holding material properties at each individual particle. This
13
feature helps to extend a code for more than two-phase flows, such as those
14
in liquid-liquid-gas (such as oil-water-gas) and liquid-liquid-solid (oil-water-
15
salt/sediment) [11]. Additionally, SPH inherently captures the interface and 3
ACCEPTED MANUSCRIPT
is not tracking it like those in Level-set method. Therefore, the interface is
17
always sharp, is not smeared, the mass of each phase is rigorously conserved
18
during the simulation and no need for complex coding for large topological
19
change of interface (such as those in folding, merging and break up). This is
20
very important especially when the code is extended for 3D problems. Sec-
21
ondly, SPH naturally incorporate singular forces and coefficient discontinu-
22
ities into the numerical scheme, therefore, the interface boundary condition
23
is inherently implemented. Thirdly, since the field properties’ derivatives
24
will transform to derivatives of some analytical functions, so-called kernel
25
function, they have exact solutions and can be naturally integrated into the
26
scheme. And finally, the convective term in discretisation of the balance equa-
27
tions is absent in the numerical approximation scheme due to its Lagrangian
28
nature. The latter makes SPH a suitable candidate for fast-dynamic flows.
29
For multi-phase flows, the most important interfacial force is the surface ten-
30
sion. It is known that the surface tension may depend on the temperature
31
or concentration. If there is a temperature or concentration gradient tan-
32
gential to the interface between two fluids, then, an additional force known
33
as Marangoni force, arises. However, and in spite of numerus articles on
34
multiphase flow problems in the context of SPH [12, 13, 14, 15], in literature
35
there are only a few publications on Marangoni force (i.e. the gradient of the
36
surface tension) [16, 17]. Marangoni force due to a concentration gradient
CE
PT
ED
M
AN US
CR IP T
16
at the interface was first discussed by Adami et al. [16]. They used the
38
conservative formulation of the continuum surface force (CSF) model, the
39
so-called continuum surface stress (CSS) model, and investigated the motion
40
and surface concentration of a bubble.
AC 37
4
ACCEPTED MANUSCRIPT
In the context of SPH, thermocapillary motion was only recently studied by
42
Tong et al. [17] for 2D problems. Tong et al. presented results for the migra-
43
tion velocity of a single droplet in thermocapillary flow. In this case a droplet
44
is accelerated by gradients of the surface tension coefficients. For a very low
45
resolution, the migration velocity presented by Tong et al. match results
46
of Ma & Bothe [18] but they didn’t provide a convergence study. Tong et
47
al. [17] presented a model based on the non-conservative formulation of the
48
CSF model (the same that we use here). They investigated thermocapillary
49
migration of a droplet where heat transfer was included using an enthalpy
50
balance. Our model is similar to the model of Tong et al. [17] except that
51
we use the energy balance in temperature form and neglect phase transition.
52
Therefore, the current work aims at the proposition and the development of
53
a consistent SPH method that is capable of treating complex flow physics in-
54
cluding Marangoni forces and thermocapillary effects. At the same time this
55
article provides step by step validation test cases that can serve as bench-
56
marks for both SPH and thermo-fluids communities. It is noted that unlike
57
Adami et al. [16] the proposed scheme is based on the non-conservative for-
58
mulation of the CSF model. Additionally, this manuscript extend the results
59
of Tong et al. [17] to 3D problems and it provides more details on the con-
60
vergence of SPH results when dealing with thermocapillary forces.
61
This article is organized as follows: First we introduce the balance equa-
CE
PT
ED
M
AN US
CR IP T
41
tions of the multiphase system. Then we introduce the SPH method and
63
the discrete balance equations, including boundary conditions. Afterwards
64
numerical stability is discussed introducing kernel consistency, particle con-
65
sistency and volume conserving ISPH scheme. Numerical results are divided
AC 62
5
ACCEPTED MANUSCRIPT
into single-phase and multi-phase sections. First, numerical implementation
67
is validated for single-phase flow. Afterwards, surface tension is validated
68
for multi-phase flow and results for thermo-capillary droplet migration are
69
presented and validated against analytical solution and related numerical
70
solutions.
71
2. Balance equations
72
In this section we present the balance equations. A complete thermodynamic
73
consistent form of the component, momentum and energy balance, as found
74
in Thieulot et al. [19], is very complex and has parameters that are not
75
known for the system under consideration. Therefore, we make the following
76
assumptions to only reduce complexity without neglecting important mech-
77
anisms.
78
Firstly, we assume that the fluids are immiscible and incompressible. Sec-
79
ondly, we neglect phase transitions and assume that the surface tension de-
80
pends linearly on temperature. Using these assumptions we introduce the
81
balance equations. Then, we review these derivations in Lagrangian form for
82
our Lagrangian numerical framework.
83
2.1. Continuity equation
84
The continuity equation is [20]
AC
CE
PT
ED
M
AN US
CR IP T
66
dρ + ∇ · (ρ~u) = 0, dt
(1)
85
where ρ and u are the mass density and the mass average velocity of the fluid,
86
respectively. We formulate the continuity equation in a Lagrangian reference
6
ACCEPTED MANUSCRIPT
frame. Therefore, we rewrite Eq. 1 using the material derivative Dρ + ρ∇ · ~u = 0, Dt
88
D Dt
(2)
CR IP T
87
We define an incompressible system with constant density as ∇ · ~u = 0,
(3)
where the density is independent of the pressure, hence, the divergence of
90
the velocity is zero.
91
2.2. Momentum equation
92
The momentum balance of the multi-phase system is [19]
AN US
89
(4)
M
d (ρ~u) + ∇ · (ρ~u~u) = −∇p + ∇ · Πviscous + F~body + ∇ · Πcapillary , dt
with the pressure p, the viscous stress tensor Πviscous , a body force, e.g.
94
gravity F~body = ρg, and the capillary stress tensor Πcapillary . The first three
95
terms on the right hand side are classical terms in the momentum balance
96
as in single-phase flow [20]. The last term in Eq. 4 accounts for the presence
97
of interfaces. With the continuity equation (Eq. 1) the momentum balance
98
for an incompressible, Newtonian fluid reduces to
CE
PT
ED
93
d~u + ρ (~u · ∇) ~u = −∇p + η∇2~u + F~body + ∇ · Πcapillary , dt
(5)
AC
ρ
99
with the dynamic viscosity η.
100
The capillary stress tensor in Eq. 5 arises due to the presence of the interface,
101
which in general is diffuse on molecular scale because density continuously
102
varies between different phases. Here, we consider the sharp interface limit 7
ACCEPTED MANUSCRIPT
of the capillary stress tensor. The sharp interface limit may be seen as
104
projecting the integral over the diffuse interface on the interface area (plane
105
in 3D). A review of a mathematical derivation of the sharp interface limit
106
from a diffuse formulation is shown in [21]. The main result is
CR IP T
103
∇ · Πcapillary = σκ~n ˆ + ∇S σ δ,
(6)
with the surface tension coefficient σ, the curvature κ and the gradient of
108
the surface tension tangential to the surface ∇S σ. Additionally, ~n ˆ is the unit
109
normal to the interface and δ is the Dirac-Delta function
δ=
AN US
107
1 at the interface 0 in the bulk
.
(7)
Finally, the momentum equation in Lagrangian form with a sharp interface
111
is
ρ
ED
M
110
D~u = −∇p + η∇2~u + F~body + σκ~n ˆ + ∇S σ δ. Dt
(8)
2.3. Energy balance
113
The energy balance in temperature form, neglecting viscous dissipation and
114
surface energy contributions, is
AC
CE
PT
112
ρcp
DT = ∇ (λ∇T ) , Dt λ , ρcp
115
with the thermal diffusivity of k =
116
the specific heat capacity at constant pressure, cp .
8
(9)
the thermal conductivity of λ and
ACCEPTED MANUSCRIPT
3. Smoothed Particle Hydrodynamics method
118
As mentioned earlier, SPH method was first introduced by Monaghan and
119
Gingold [2] and Lucy [1] in 1977. Since then, several reviews are published
120
presenting the basics of SPH. For a general introduction to the SPH method
121
and its various applications, we refer to the reviews of [22, 10, 9]. Here, we
122
briefly present the SPH initiatives and the relevant discretizations that is
123
used in the current study.
124
The starting point of SPH is coming from an exact mathematical solution
125
where we can estimate any quantity of A (~x) at a given position ~x by mul-
126
tiplying it with a Dirac-Delta function and taking its integration A (~x0 ) over
127
the whole volume [2]
A(~x0 )δ (|~x − ~x0 |) d~x0 ,
M
A (~x) =
Z
AN US
CR IP T
117
V
(10)
where ~x0 is the position of any point in the volume V . The Dirac-Delta
129
function is defined to be unity for ~x0 = ~x and zero otherwise.
130
In the next step, this analytical equation undergoes two approximations. The
131
first, known as kernel approximation, is done by replacing the Dirac-Delta
PT
function with a weighting function, so-called kernel function W (h, |~x − ~x0 |) Z A (~x) = A(~x0 )W (h, |~x − ~x0 |) d~x0 , (11)
CE
132
ED
128
V
with h being the radius of the support volume, so-called smoothing length.
134
In our case, where an initial regular particle distribution in a Cartesian order
135
with the spacing of L0 is used, we define a constant value of h = 2.05L0 .
136
In this work we use a Wendland C2 kernel function of order 4 with rc = 2h
AC
133
9
ACCEPTED MANUSCRIPT
[23]
q 4 1 − (1 + 2q) for 0 ≤ q ≤ 2 fw 2 W (h, |~x − ~x0 |) = d . h 0 for 2 < q
(12)
CR IP T
137
138
Here, d is the dimension and fw is a normalization constant with fw = 43 ,
139
7 , 21 4π 16π
140
the non-dimensional smoothing length [24]. The first derivative of this kernel
141
function is
for d = 1, d = 2, and d = 3 dimensions, respectively. q =
for 0 ≤ q ≤ 2
AN US
−5q 1 − q 3 2
∂W (h, |~x − ~x0 |) fw = d+1 ∂x h 0
.
|~ x−~ x0 | h
is
(13)
for 2 < q
Next is the particle approximation where the integration in Eq. 11 is replaced
143
by a summation over finite number of interpolation points, known as particles
144
A (~x)a =
M
142
X mb ρb
A(~xb )W (h, |~xa − ~xb |) .
(14)
ED
b
Here, the indexes a and b represent the particle of interest and the particles
146
in the vicinity of a (see Fig. 1) called neighbor particles. mb and ρb are
147
the mass and density of particle b, respectively. In the following the kernel
148
function is abbreviated as W (h, |~xa − ~xb |) = Wab with |~xa − ~xb | = rab being
149
the distance between particle a and b.
150
An example, using Eq. 14 with A = ρ, is the calculation of the density of
151
particle a [22]
AC
CE
PT
145
152
153
ρa =
X
mb Wab .
(15)
b
Note that particle a should also be considered in its vicinity because Waa 6= 0. Due to the particle approximation, the 0th moment (unity condition of the 10
CR IP T
ACCEPTED MANUSCRIPT
AN US
Figure 1: (Left) Scheme of discrete particle distribution around a particle of interest,
particle a, and its particles in the vicinity b within the support of the kernel function. Dashed particles are outside the support domain. (Right) Schematic representation of a kernel function W with a support domain of 2h.
kernel function) is not guaranteed [10]. To guarantee it, we may use a renor-
155
malization function, the shepard kernel, to ensure the unity condition [25, 26].
156
This is very important if the support domain is truncated e.g. at free sur-
157
faces or when surface tension has to be calculated as we will see later. For
158
bulk flow the error is normally in the range of a few percent. We discuss this
159
point in section 5.1.
160
As mentioned before, one strong advantage of the SPH method is that the
161
gradient of a quantity A (~x) can be derived from the gradient of the kernel
162
function [2]. One can formulate the first derivative using [22]
AC
CE
PT
ED
M
154
163
or [6]
X mb
∇A (~x)a =
1 ∇A (~x) ρ
b
a
=
ρb
(Ab − Aa ) ∇W (h, rab ) ,
X
mb
b
Ab + Aa ∇W (h, rab ) , ρa · ρb
11
(16)
(17)
ACCEPTED MANUSCRIPT
164
where the gradient of the kernel function ~rab ∂W , rab ∂x
(18)
CR IP T
∇W (h, rab ) =
is analytically known from Eq. 13. Here, ~rab is the distance vector between
166
particle a and b pointing from particle a to b and rab is the length of this
167
vector. Eq. 16 guarantees that the gradient of a constant field is exactly zero
168
even if numerical errors are present. Linear momentum is only conserved if
169
the 1th moment of the kernel function is exactly zero [10]. Eq. 17 conserves
170
linear momentum exactly as long as the kernel function is symmetric but
171
the gradient of a constant field is only zero if the 1th moment of the kernel
172
function is exactly zero. In literature Eq. 16 is often used for the divergence
173
of a velocity field, where the gradient of a constant field is very important,
174
and Eq. 17 is often used for the pressure gradient because linear momentum
175
is more important in this case.
176
For the second derivative of a quantity, as introduced by Brookshaw [27], a
177
hybrid formulation for the second derivative is used with a first order deriva-
178
tive of the kernel function in combination with a finite difference approxima-
179
tion
(19)
where α is a variable that does not depend explicit on position but e.g.
AC
180
X mb ~rab ~ (~x) = ~a − A ~b , ∇ · α∇A (αa + αb ) 2 ∇a Wab · A ρb rab a b
CE
PT
ED
M
AN US
165
181
may depend on other properties of a fluid. We abbreviate the gradient with
182
respect to particle a as ∇a Wab = ∇a W (h, rab ). In contrast to other formula-
183
tions of the second derivative [28] this formulation proved to be less sensitive
184
to particle disorder. 12
ACCEPTED MANUSCRIPT
4. Discrete balance equations
186
In this section we present the discrete balance equations using the operators
187
from section 3, followed by presenting the projection method used in this
188
work to calculate the pressure.
189
4.1. Continuity equation
190
In SPH, there exist two different ways [22] to formulate the discrete continuity
191
equation. One way is to directly use the discrete continuity equation [29, 30]
192
using Eq. 16
AN US
CR IP T
185
X mb Dρa (~ub − ~ua ) ∇a Wab . = ρa Dt ρ b b
(20)
A disadvantage of this equation is the error of conservation due to integration
194
errors. Another formulation of the discrete continuity equation is [12]
M
193
ED
ρ a = ma
X
Wab ,
(21)
b
where the density is directly calculated using the mass and volume of par-
196
ticle a. Note that the summation of the kernel function represents the re-
197
ciprocal volume 1/Va . Advantages of this formulation are its applicability
198
to multi-phase systems with density differences between the phases and its
199
conservation property. Therefore, we use Eq. 21 throughout this work.
200
4.2. Momentum balance
AC
CE
PT
195
201
In the momentum balance (Eq. 8) we need first-order derivatives for the
202
pressure term, for the color to calculate the normals, for the normals to
203
calculate the curvature and for the surface tension to calculate the Marangoni
204
force. In addition, we need a second-order derivative for the divergence of 13
ACCEPTED MANUSCRIPT
205
the viscous stress tensor.
206
We use the positive formulation of the gradient (Eq. 17) for the pressure
207
force
∇p ρ
a
=
X mb (pb + pa ) ∇a Wab , ρ a ρb b
CR IP T
(22)
208
with the pressure pa and pb of particle a and b, respectively. It is common
209
in literature to use the positive formulation of the gradient especially for
210
the pressure force to ensure linear momentum conservation. For a constant
213
214
AN US
212
pressure field the gradient of the pressure is only zero for a constant pressure P p = 0 or symmetric particle distribution, which results in b ∇a Wab = 0. The discrete formulation of the viscous term using the second-order derivative (Eq. 19) is
η 2 ∇ ~u ρ
=
a
X mb ~rab (ηa + ηb ) 2 ∇a Wab · (~ua − ~ub ) . ρa ρb rab b
(23)
M
211
The body force is straight forward since it only depends on particle a.
216
The last two terms in Eq. 8 account for interface forces. To formulate surface
217
tension at a sharp interface, Brackbill et al. [31] and later, in the context of
218
SPH, Morris [32] introduced the so-called Continuum Surface Force (CSF)
219
approach. In the CSF approach the surface force, that is a force per area, is
220
reformulated as a force per volume. Therefore, the surface tension is numer-
221
ically smoothed out around the sharp interface. We need the normal vector
222
~n at the interface between phases to calculate the magnitude and direction
223
of the interface forces. The normal vector at an interface is calculated from a
224
color field where a unique color is assigned to particles of every phase. We use
225
the negative formulation of the gradient (Eq. 16) to calculate the gradient of
226
the color field because it is important that the gradient vanishes if the color
AC
CE
PT
ED
215
14
ACCEPTED MANUSCRIPT
228
is constant. Otherwise, we get surface tension due to numerical artifacts in the bulk. Hence, interface normal vector ~na is calculated as X mb (Cb − Ca ) ∇C ~na = = ∇a Wab , [C] a ρb [C] b
(24)
CR IP T
227
with the color C and its jump across the interface [C]. Whether the color is a
230
smooth or a sharp function across the interface, the jump of the color function
231
represents the absolute difference between the maximum and minimum of the
232
color function. The normalized normal vector is
AN US
229
~n ~n , ˆ= |~n| 233
(25)
with the magnitude of the normal vector |~n|. Then, the curvature is (26)
M
X mb ~n ˆ b − ~n ˆ a ∇a Wab . κa = − ∇ · ~n ˆ =− ρb a b
We again use the negative formulation because it is important that the cur-
235
vature is zero for equal normal vectors. Note that we truncate the field of
236
normals to avoid numerical artifacts and only consider normals with magni-
237
tude larger than /h with = 0.01.
238
The last term in Eq. 8 accounts for gradients of the surface tension tangential
239
to the interface. We use a projection of the gradient of the surface tension
240
into the tangential direction of the interface
AC
CE
PT
ED
234
241
(∇S σ)a = ∇σ − ∇σ · ~n ˆ ~n ˆ , a
(27)
with the gradient of the surface tension coefficient σ (∇σ)a =
X mb b
ρb
(σb − σa ) ∇a Wab .
15
(28)
ACCEPTED MANUSCRIPT
242
The delta function in Eq. 8 is unity at the interface. In the discrete formu-
243
lation we smooth out the delta function (29)
CR IP T
δa = |~na |.
This leads to a distribution of the surface force in a volume that in the limit is equal to that of a delta function.
Finally, the used discrete momentum balance in the current work is
with f~body =
245
4.3. Energy balance
M
~body F . ρ
244
The discrete energy balance for convective and diffusive heat transfer is X mb ~rab DT = (ka + kb ) 2 ∇a Wab · (Ta − Tb ) , (31) Dt a ρ r b ab b
PT
ED
246
AN US
X mb X mb ~rab D~ua =− (pb + pa ) ∇Wab + (ηa + ηb ) 2 ∇a Wab · (~ua − ~ub ) Dt ρa ρb ρa ρb rab b b |~na | σa κa~n ˆ a + (∇S σ)a , (30) + f~body,a + ρa
using the second-order derivative (Eq. 19) where Ta and λa are the temper-
248
ature and heat conductivity of particle a.
249
4.4. Estimation of pressure
250
The pressure in the momentum balance is still unknown. For incompress-
251
ible fluids there are two ways to calculate the pressure. One way is to use
252
an equation of state for a so-called weakly compressible fluid and is known
253
as WCSPH. It is a stiff equation that relates a small density variation to a
254
pressure variation. One example is the Tait equation [33]. This equation is
AC
CE
247
16
ACCEPTED MANUSCRIPT
commonly used in SPH to simulate the pressure of water [6, 22].
256
An alternative way to determine the pressure for a truly incompressible fluid
257
is to use a projection method. We can apply the Helmholtz-Hodge decom-
258
position [34] on the momentum balance. The momentum balance is divided
259
into a curl-free and divergence-free part. In a predictor step we first calculate
260
an intermediate solution of velocity and density using the curl-free part. The
261
curl-free part of the momentum balance (Eq. 30) consists of all terms except
262
the pressure term. In a second step we solve a Pressure Poisson Equation
263
(PPE) to calculate the pressure and then correct the intermediate velocity
264
using the pressure term to get the final velocity of the time step. In SPH
265
this projection method was introduced by Cummins and Rudman [3] and is
266
usually called ISPH. The details of the projection method are discussed in
267
literature [8, 26, 29, 35] and we only review the discrete equations used in
268
this work.
AN US
M
The curl-free part acceleration ~aa,∗ in the momentum balance is X mb ~rab |~na | ~aa,∗ = (ηa + ηb ) 2 ∇a Wab ·(~ua − ~ub ) +f~body,a + σa κa~n ˆ a + (∇S σ)a . ρa ρb rab ρa b (32)
PT
ED
269
CR IP T
255
Here, the subscript asterisk indicates the intermediate step. We calculate the
271
intermediate velocity using an explicit Euler scheme
CE
270
(33)
with the time step ∆t. Then, we determine the new particle position
AC 272
~ua,∗ = ~uta + ~aa,∗ ∆t
~xa,∗ = ~xta + ~ua,∗ ∆t,
(34)
273
also using an explicit Euler step. The intermediate particle position is only
274
necessary because we calculate the intermediate density using Eq. 21 instead 17
ACCEPTED MANUSCRIPT
of the continuity equation (Eq. 20). When the intermediate density is calcu-
276
lated we put the particles back to the old particle position ~xta and calculate
277
the pressure using the PPE ∇·
279
280
∇p ρ∗
=
∇ · ~u∗ . ∆t
(35)
We use the divergence of the velocity as a condition for incompressibility. The discrete Laplace of the pressure using Eq. 19 is (for α = 1) X mb ~rab ∇p · (pa − pb ) 2 ∇a Wab , = 2 ∇· ρa,∗ a ρb ρa rab b
AN US
278
CR IP T
275
and the right hand side is
1 X mb ∇ · ~ua,∗ = (~ub,∗ − ~ua,∗ ) ∇a Wab . ∆t ∆t b ρb
(36)
(37)
The PPE requires the solution of a linear equation system. We use a sta-
282
bilized version of the biCGStab method from the PETSc library [36] with a
283
boomerAMG pre-conditioner from the HYPRE library [37].
284
Then, with the acceleration of the corrector step ~aa,∗∗ , using the pressure
285
term
ED
M
281
(38)
~ut+1 = ~ua,∗ + ~aa,∗∗ ∆t. a
(39)
The final position of particle a is
AC
287
X mb (pb + pa ) ∇a Wab , ρ a · ρb b
we determine the final velocity with an explicit Euler step
CE
286
PT
~aa,∗∗ = −
= ~xta + ~xt+1 a
~uta + ~ut+1 a ∆t. 2
(40)
288
Here, a semi-implicit step with the velocities at time t and t+1 is used. More
289
details on the projection method in SPH is found in Keller [38]. 18
ACCEPTED MANUSCRIPT
4.5. Boundary conditions
291
Throughout this paper we use impermeable boundary conditions. For veloc-
292
ity we either use slip or no-slip conditions. For temperature and pressure we
293
use Dirichlet or Neumann boundary conditions, while for the gradient of the
294
color C only Neumann boundary conditions are used, because the fictitious
295
boundary condition should not influence the gradient of C. As we will see in
296
Sec. 5.1, we do not need boundary conditions for the normal and the surface
297
tension coefficient.
298
Beside different approaches for boundary treatment in SPH [7, 29, 39, 40,
299
41, 42, 43], we use wall boundary conditions from Cummins and Rudman
300
[3] (mirror particles) at straight walls. Due to the integral character of SPH
301
the wall partially needs to be discretized up to a certain distance from the
302
fluid domain. In the mirror particle approach, the fluid particles are mir-
303
rored in every time step across the boundary. In this way an image of the
304
fluid particles represent the wall particles. This is sketched in Fig. 2. The
305
properties of the wall particles are manipulated to incorporate Dirichlet and
306
Neumann boundary conditions. The boundary condition for any property
307
Ψ ∈ [~u, p, T, C] at the wall is considered in the property of an image particle.
308
The property of an image particle is ~ 0a = 2Ψ ~ w − δw Ψ ~a Ψ
AC
CE
PT
ED
M
AN US
CR IP T
290
309
~ w and the sign with the value of the wall Ψ +1 for Dirichlet boundary conditions δw = −1 for Neumann boundary conditions with Ψ ~w = 0 19
(41)
.
(42)
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 2: Scheme of mirror boundary conditions. The fluid particles are mirrored across
PT
the boundary and the images represent the wall particles (dotted particles). Their properties are manipulated to incorporate Dirichlet and Neumann boundary conditions. Here
AC
CE
no-slip boundary conditions of an impenetrable wall is shown.
20
ACCEPTED MANUSCRIPT
More details can be found in [29, 42]. This technique works well on straight
311
walls but leads to errors on curved boundaries. Nevertheless, in the present
312
test cases we only need straight walls.
313
5. Numerical stability
314
It is known that the unity condition for the SPH kernel function is not guar-
315
anteed for random particle distributions. One way to guarantee the unity
316
condition of the kernel function was proposed by Bonet & Lok [25]. It is
317
called the corrected SPH approach and is presented in section 5.1.
318
In practical applications we observe that a corrected SPH approach is not
319
enough to guarantee stable and accurate simulations. The reason is that
320
in SPH we must satisfy two conditions. The first condition is the already
321
mentioned consistency of the kernel function. The second condition is the
322
consistency of the particle distribution [10]. If the spacing between particles
323
is larger than the interaction radius, the particles become disconnected and
324
even corrected SPH approaches can’t prevent that.
325
In a proper ISPH approach, particles that are initially on a regular Cartesian
326
grid, may move on the same trajectory. In front of a stagnation point, parts
327
of the particles move in different directions and get disconnected when their
328
distance is larger than the cut-off radius. As a result, a large void space is
329
formed between them. Therefore, we need a numerical correction that rear-
AC
CE
PT
ED
M
AN US
CR IP T
310
330
ranges the particles if necessary without influencing the balance equations.
331
In literature there exist two approaches to improve the particle consistency
332
in ISPH. The first approach was introduced by Hu & Adams [35]. It is the
333
so-called combined DIDF. The other approach was introduced by Xu et al. 21
ACCEPTED MANUSCRIPT
[44] and Shadloo et al. [43], used for the first time in the concept of mul-
335
tiphase flows by Shadloo et al. [45, 46] and modified to consider its arts
336
effects for interfacial flows by Lind et al. [47]. It is called the particle shifting
337
technique. Using these approaches improve the particle consistency and also
338
the 0th and 1th moment of the kernel function and the gradient of the kernel
339
function. We summerize both approaches in section 5.2 and 5.3.
340
5.1. Corrected kernel function and corrected gradient of kernel function
341
In this section we present the corrected SPH approach, introduced by Bonet
342
& Lok [25] as mixed kernel and gradient correction. Further details about
343
the implementation and accuracy are found in [38]. The idea of the corrected
344
SPH approach is to renormalize the kernel function and the gradient of the
345
kernel function to guarantee the unity condition for the kernel and the 1th
346
moment
M
AN US
CR IP T
334
X mb
∇a Wab = 0,
(43)
˜ ab = Wab , W Sa
(44)
X mb
(45)
with the Shepard kernel Sa [48]
AC
CE
348
ρb
for the gradient of the kernel. The renormalized kernel function is
PT
347
ED
b
Sa =
b
ρb
Wab ,
349
and the uncorrected kernel function Wab . In the nominator the interaction
350
between all particles including the interaction of particle a with itself need
351
to be considered.
352
In most of the terms in the balance equations we need the gradient of the 22
ACCEPTED MANUSCRIPT
353
kernel function. In the mixed kernel and gradient correction of Bonet & Lok,
354
the corrected gradient of the corrected kernel is
355
(46)
CR IP T
˜ aW ˜ ab = La ∇a W ˜ ab , ∇
with the correction matrix, that ensures symmetric vector direction, !−1 X mb ˜ ab ⊗ ~rab ∇a W , La = ρ b b
(47)
and the corrected gradient of the kernel, that ensures equal magnitude of the
357
vector,
AN US
356
˜ ab = ∇a W 358
with γ=
∇a Wab − γ , Sa
P
mb b ρb ∇a Wab
Sa
.
(48)
(49)
An advantage of the corrected SPH approach is that a constant pressure field
360
has zero gradient because of Eq. 43. In all previous discrete formulations in
361
section 4 we may replace the kernel or gradient of the kernel function by Eqs.
362
44 or 46, respectively, if the corrected SPH approach is used. We explicit
363
exclude Eq. 24 because a corrected SPH kernel leads to wrong normals.
364
In addition we may use the corrected SPH approach if we only consider a
365
subset of the particles in the domain. For the curvature and the gradient of
366
the surface tension we use only particles with normal vectors of magnitude
367
larger than a threshold /h > 0.01 (see section 4.2). In all other cases we use
AC
CE
PT
ED
M
359
368
all particles in the vicinity, including mirrored particle at the boundary.
369
5.2. Density-invariant divergence-free (DIDF) approach
370
Hu & Adams [35] introduced an approach to improve particle consistency of
371
the ISPH method for bulk flow problems. The idea is to split up the time 23
ACCEPTED MANUSCRIPT
integration of velocity and position into two parts in the corrector step of
373
the projection method in ISPH. In a first integration step the left side of the
374
continuity equation 1 Dρ = ∇ · ~u, ρ Dt
375
is used as criterion for incompressibility in the PPE t τ2 ∇p ρ0 − ρ∗,t+1 ∇· = a 0a , 2 ρ∗ a ρa
CR IP T
372
(50)
(51)
where τ is the time step and ρ0a is the initial density of particle a. In this
377
first step, we calculate a pressure to correct the position of the particles on
378
the half time step using an explicit Euler step. The intermediate velocity is
AN US
376
~ua,∗ = ~uta + 0.5 · a~a,∗ ∆t, and the intermediate position is
M
379
380
In a second integration step, the divergence of the velocity is used as criterion for incompressibility in the PPE t ∇p τ∇ · = ∇a · ~ut+1 ∗ . ρ∗ a
(54)
CE
PT
381
(53)
ED
~xa,∗ = ~aa,∗ ∆t.
(52)
In this second step a pressure to correct the velocity of the particles is cal-
383
culated. The final velocity is
AC
382
384
~ut+1 = ~ua,∗ + ~aa,∗∗ ∆t, a
(55)
~xt+1 = ~uta + ~ut+1 aa,∗∗ (∆t)2 . a ∆t + 0.5 · ~ a
(56)
and the position is
24
ACCEPTED MANUSCRIPT
In the DIDF approach we need to solve the PPE two times per time step.
386
Therefore, the numerical effort increases drastically because the solution step
387
of the LES takes approximately 70 % of the hole computation time. Further-
388
more, the approach is not applicable to free surface flows [49]. On the other
389
hand the DIDF approach improves the conservation of volume because both
390
particle position and velocity are considered as criterion for incompressibility.
391
The DIDF approach is used here to improve volume conservation during
392
the simulation. Nevertheless, particles still remain on their trajectory and
393
split up at stagnation points. Therefore, another approach to regularize the
394
particles is necessary.
395
5.3. Particle Shifting
396
Particles in SPH move on their streamline due to their Lagrangian nature.
397
Therefore, particles split near a stagnation point leaving a void space between
398
two streamlines. This leads to large errors in the particle consistency of SPH.
399
The particle shifting technique is introduced to geometrically shift particles
400
slightly away from their streamlines into a void space [44, 43]. As a result
401
the particles are more regularly distributed in the domain. We introduce the
402
shifting algorithm in the following.
403
The idea is to shift particles slightly into a void space after every time step.
404
The total shifting vector ~ra,s is
AC
CE
PT
ED
M
AN US
CR IP T
385
~ a, ~ra,s = CαR
(57)
405
with the arbitrary constant C = 0.04, the shifting magnitude α = |~umax |dt,
406
the magnitude of the maximum particle velocity |~umax |, the time step dt, and 25
ACCEPTED MANUSCRIPT
407
the shifting vector ~a = R
Nb X r¯a2 ˆ ~r , 2 ab rab b
(58)
where ~rˆab , rab and r¯a are the unit distance vector between particle a and b,
409
the distance between both particles and the average particle spacing in the
410
vicinity defined as Nb
1 X r¯a = rab . Nb b
CR IP T
408
(59)
In this approach it is assumed that all particles represent the same volume.
412
Therefore, a geometrically uniform distribution is the best particle distribu-
413
tion.
414
Some properties of the fluid are connected to the particles and their posi-
415
tion. If we shift the particles we need to correct the hydrodynamic variables
416
Ψ to be consistent and not influence the balance equations. Following [44],
417
a Taylor series of the hydrodynamic variable
ED
M
AN US
411
Ψa0 = Ψa + δ~raa0 · (∇Ψ)a + · · ·
(60)
is used to approximate the property of a particle at the new position a0 af-
419
ter shifting. We use a linear interpolation and neglect quadratic and higher
420
terms. Hydrodynamic variables used in this work are the velocity and tem-
421
perature. The density is not a hydrodynamic variable in this context because
422
we calculate the density by summation of masses in the vicinity instead of
AC
CE
PT
418
423
calculating the continuity equation directly. Therefore, we are able to calcu-
424
late the density at any point.
425
The particle shifting technique is independent of the time integration and
426
does not affect the balance equations because we shift the particles only at 26
ACCEPTED MANUSCRIPT
the very end of a time step and correct the particle properties accordingly.
428
The shifting technique in its original form assumes that every particle is an
429
identical interpolation point in the domain. In a multi-phase system this is
430
not valid any more because of the presence of an interface. In SPH the inter-
431
face between two phases is defined by the type (color) of particles. Therefore,
432
the interface is implicitly moved when particles change their position. If we
433
shift the particles, we implicitly shift the interface between phases as well. To
434
overcome this shortcoming we modify Eq. 57 near an interface between two
435
phases and apply shifting only in the tangential direction to the interface.
436
This can be done by only consider the tangential part of Eq. 57. A similar
437
extension was proposed by Lind. et al. [47]. They modified the shifting
438
technique for free surface flows based on a Fickian approach.
439
5.4. Time integration
440
In this work we use a predictor-corrector scheme to integrate the momentum
441
balance equation in time, as introduced earlier in the context of the projection
442
method to calculate the pressure (Sec. 4.4 and 5.2).
443
For the energy balance (Eq. 31) we use an explicit Euler integration scheme
444
during the predictor step.
445
For an explicit integration scheme there exist criteria for the maximum time
446
step because the algorithm is not unconditionally stable. Therefore, the
447
following time step criteria are used throughout the work. The time step
448
criteria due to convective transport (CFL criterion) is
AC
CE
PT
ED
M
AN US
CR IP T
427
∆tCF L = αCF L
27
L0 , umax
(61)
ACCEPTED MANUSCRIPT
450
due to viscous diffusion is ∆tvisc = αdif f
L20 νmax
∆theat = αdif f
L20 , λmax
and due to heat diffusion is
(62)
CR IP T
449
(63)
with αCF L = 0.05 and αdif f = 0.0625. Additionally, umax , νmax and λmax are
452
the magnitude of the maximum velocity, the maximum kinematic viscosity
453
and the maximum heat conductivity, respectively.
454
6. Single-phase flows
455
In this section we validate the momentum balance of the proposed model
456
and its discrete formulation using the SPH method. Since the SPH method
457
itself is well validated in literature we focus on the modifications made in
458
this work. After showing the numerical stability of the proposed method
459
and validation of the diffuse transport of a scalar property, we start with
460
the very common isotherm lid-driven cavity test case [50, 51] to demonstrate
461
accurate implementation of the momentum balance, boundary conditions
462
and numerical stability approaches. It is one of the standard test cases of
463
computational fluid dynamics simulation codes. It enables us to validate no-
464
slip boundary conditions for velocity as well as the wall velocity.
465
Afterwards, we couple the momentum equation with the energy equation.
466
We consider a buoyancy-driven flow in a cavity where the density of the fluid
467
is temperature-dependent. In this test case a circular flow is induced due to
468
a temperature gradient [52, 53, 54]. Note that the validation of diffusive heat
469
transport in shown in section 6.3 for completeness.
AC
CE
PT
ED
M
AN US
451
28
ACCEPTED MANUSCRIPT
6.1. Numerical stability
471
Here, we first compare the three previously introduced approaches (i.e. cor-
472
rected kernel function and corrected gradient of kernel function, abbreviated
473
as corrected SPH in the following, DIDF and particle shifting approach)
474
for numerical stability and combinations of them. One common test case
475
in literature is a two-dimensional Taylor-Green vortex. In this test case all
476
boundaries are periodic and initially a divergence-free velocity profile [55]
AN US
ux = sin(2πx)cos(2πy)
CR IP T
470
477
uy = −cos(2πx)sin(2πy)
(64)
(65)
is imposed. The vortex decays due to viscous dissipation of the fluid. Details
479
about this test case are found in [44, 47]. We consider a case with Reynolds
480
number Re = 1000.
481
The magnitude of velocity and the particle distribution are shown in Fig. 3
482
at time t = 0.24s for different combinations of the previous approaches. Fig.
483
3a shows the results without any numerical stability. We see that particles
484
are aligned on their streamline and large void spaces are formed between
485
the streamlines. The other three cases in the top line (Fig. 3b-d) are cases
486
where the corrected SPH and DIDF approach without particle shifting is
487
used. Except for some small deviations we get the same particle distribution
AC
CE
PT
ED
M
478
488
as in Fig. 3a. We note that the velocity distribution as well as the decay of the
489
kinetic energy are in good agreement with the analytic solution (not shown
490
here), but at later times the solution of the PPE tends to diverge because
491
void spaces between particle get too large, e.g. larger than the cutoff radius 29
ACCEPTED MANUSCRIPT
of the kernel function, and 2D-connectivity of the interpolation is lost. If
493
we apply the particle shifting approach alone or in addition with any of the
494
other stability approaches (Fig. 3e-h) the particle distribution as well as the
495
velocity distribution is very smooth. As shown by many authors [44, 47, 56]
496
the accuracy is largely improved by the particle shifting approach and leads
497
to better particle consistency. The effect of corrected SPH and DIDF is minor
498
in this test case. It is known from literature that the corrected SPH approach
499
improves the calculation of surface tension forces and DIDF improves volume
500
conservation [57, 35]. Since both approaches do not influence the particle
501
consistency, in this work, we use the corrected SPH and particle shifting
502
approach where particle movement or surface tension is involved. The DIDF
503
approach is only used in simulations where volume conservation is crucial,
504
since this is the main advantage of the approach. Its use is highlighted
505
explicitly in the text.
506
6.2. Lid-driven cavity
507
Circular flow in a cavity was theoretically and numerically studied by many
508
researchers e.g. [44, 50, 58, 59] to mention a few of them. An isotherm flow
509
inside a wall bounded cavity, as schematically shown in Fig. 4, driven by a
510
constant velocity on the top of the cavity is considered.
511
Due to the no-slip condition on the top wall, the fluid is accelerated and a
512
circular flow is initiated. In the following we perform simulations for different
513
Reynolds numbers and compare the axial velocity profiles with results from
514
the well validated grid-based solver OpenFOAM [60] version 2.4.0. We in-
515
vestigate different resolutions to show numerical convergence of the solution
516
of SPH.
AC
CE
PT
ED
M
AN US
CR IP T
492
30
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 3: Taylor-Green vortex simulation with different numerical stability approaches, namely (a) without numerical corrections, (b) with corrected SPH approach, (c) with DIDF approach, (d) with DIDF approach and corrected SPH approach, (e) with Particle
M
Shifting approach, (f) with corrected SPH approach and Particle Shifting approach, (g) with DIDF approach and Particle Shifting approach and (h) with corrected SPH approach,
ED
DIDF approach and Particle Shifting approach. Particle distribution shown at time t = 0.24s. The Reynolds number is Re = 1000. Color indicates the magnitude of the velocity of the particles. Red indicates |u| = 1m/s and blue |u| = 0m/s. Colors are distributed by
AC
CE
PT
an unsaturated rainbow color scheme.
Figure 4: Scheme of lid-driven cavity example.
31
ACCEPTED MANUSCRIPT
517
The Reynolds number for the lid-driven cavity is defined as Re =
ρLuw , η
(66)
where L is the length of the cavity, uw is the absolute velocity of the top
519
wall and ρ and η are the density and the dynamic viscosity of the fluid.
520
We perform simulations for Re = 10, Re = 100 and Re = 1000 where we
521
alter the Reynolds number by changing the wall velocity only. The size of the
522
cavity is quadratic and the length is L = 1mm. The density and viscosity are
523
ρ = 1000kg/m3 and η = 0.01P as, respectively. For each Reynolds number
524
we vary the resolution from L0 = 16.67µm to L0 = 4.17µm that corresponds
525
to 60 and 240 particles across L. For this test case we use the DIDF and
526
particle shifting approaches to avoid void spaces in the stagnation points in
527
the upper corners.
528
The boundary conditions for the velocity at the impenetrable walls are no-slip
529
boundary conditions. From physical point of view the boundary conditions
530
for pressure at the wall are Neumann boundary conditions. As a result the
531
reference pressure is not set in the linear equation system (LES). Therefore,
532
infinite solutions exist and an additional bootstrap is necessary. There are
533
two valid options to specify it. One option is to use a constant null space of
534
the krylov subspace solver. If we use this option we cannot use the algebraic
535
multigrid method as preconditioner in the PETSc library. Instead we may use
536
a Jacobi preconditioner. Therefore, the stability of the solution for disturbed
537
particle configuration decreases because the condition of the matrix for the
538
krylov subspace method gets worse. The number of iterations also increases
539
by a factor of 100.
540
Another option is to specify a reference pressure, e.g. p = 0, at a point in the
AC
CE
PT
ED
M
AN US
CR IP T
518
32
ACCEPTED MANUSCRIPT
domain where the pressure gradient is very low or zero. For the lid-driven
542
cavity such a point may be at the bottom corners. A result of that is a
543
spurious pressure in this corner because the specified pressure may not be
544
the best choice. Since the velocity in the corner is nearly zero the effect of
545
this spurious pressure vanishes but the stability of the solution is much better
546
because the algebraic multigrid preconditioner can be used. For that reason
547
we choose the latter option. We show the pressure field later in the discussion.
548
Note that instead of choosing a corner we may choose the bottom wall and
549
apply a constant pressure on the wall. The difference, between specifying the
550
pressure in the corner or the bottom wall, on the pressure is small. Therefore,
551
we choose the bottom left corner here to minimize numerical artifacts with
552
a reference pressure p = 0P a.
553
We use the grid-based solver OpenFOAM [60] version 2.4.0 to prepare a
554
reference solution with two different meshes, a 192×192 regular mesh for
555
Re = 10 and Re = 100 and a 240×240 regular mesh for Re = 1000. The
556
mesh size was chosen after investigation of convergence of the solution. For
557
finer meshes the solution was unchanged for the selected Reynolds number
558
(not shown here). The dimensionless velocity profiles u∗y =
559
at dimensionless position x∗ =
560
For Re = 10 (Fig. 5) the velocity along the vertical and horizontal axis is
561
almost independent of the resolution. We find good agreement between SPH
x L
and y ∗ =
y L
uy uw
and Ux∗ =
ux uw
are shown in Figs. 5, 6 and 7.
CE
PT
ED
M
AN US
CR IP T
541
and OpenFOAM for the horizontal velocity (Fig. 5b), however, small devia-
563
tions at the extreme value for the vertical velocity (Fig. 5a). One indication
564
for convergence is the position of the center of the vortex. The relative error
565
between the position of the vortex for the highest resolution in SPH and the
AC 562
33
CR IP T
ACCEPTED MANUSCRIPT
OpenFOAM (192x192)
0.15
SPH (60x60)
0.1
SPH (90x90)
0.05
SPH (120x120)
0
SPH (240x240)
-0.05
AN US
dimensionless velocity uy* [-]
0.2
-0.1 -0.15 -0.2 0
0.2
0.4
0.6
0.8
1
dimensionless position x* [-]
0.6 0.4 0.2
CE
0 -0.4
M OpenFOAM (192x192)
ED
0.8
PT
dimensionless position y* [-]
1
-0.2
SPH (60x60) SPH (90x90) SPH (120x120) SPH (240x240)
0
0.2
0.4
0.6
0.8
1
dimensionless velocity ux* [-]
AC
Figure 5: Velocity profiles for Re = 10 of a lid–driven cavity: (a) vertical velocity vs. x∗ at y ∗ = 0.5, and (b) horizontal velocity vs. y ∗ at x∗ = 0.5.
34
CR IP T
ACCEPTED MANUSCRIPT
0.2 0.1 0.05 0
OpenFOAM (192x192)
-0.05 -0.1
SPH (60x60)
-0.15
SPH (90x90)
-0.2
SPH (120x120)
-0.25
SPH (240x240)
-0.3 0
0.2
AN US
dimensionless velocity uy* [-]
0.15
0.4
0.6
0.8
1
dimensionless position x* [-]
0.6 0.4 0.2
CE
0 -0.4
M OpenFOAM (192x192)
ED
0.8
PT
dimensionless position y* [-]
1
-0.2
SPH (60x60) SPH (90x90) SPH (120x120) SPH (240x240)
0
0.2
0.4
0.6
0.8
1
dimensionless velocity ux* [-]
AC
Figure 6: Velocity profiles for Re = 100 of a lid–driven cavity: (a) vertical velocity vs. x∗ at y ∗ = 0.5, and (b) horizontal velocity vs. y ∗ at x∗ = 0.5.
35
CR IP T
ACCEPTED MANUSCRIPT
0.4 0.2 0.1 0
OpenFOAM (240x240)
-0.1 -0.2
SPH (60x60)
-0.3
SPH (90x90)
-0.4
SPH (120x120)
-0.5
SPH (240x240)
-0.6 0
0.2
AN US
dimensionless velocity uy* [-]
0.3
0.4
0.6
0.8
1
dimensionless position x* [-]
0.6 0.4 0.2
CE
0 -0.4
M OpenFOAM (240x240)
ED
0.8
PT
dimensionless position y* [-]
1
-0.2
SPH (60x60) SPH (90x90) SPH (120x120) SPH (240x240)
0
0.2
0.4
0.6
0.8
1
dimensionless velocity ux* [-]
AC
Figure 7: Velocity profiles for Re = 1000 of a lid–driven cavity: (a) vertical velocity vs. x∗ at y ∗ = 0.5, and (b) horizontal velocity vs. y ∗ at x∗ = 0.5.
36
ACCEPTED MANUSCRIPT
567
reference is u∗x = 0.60%, u∗y = 0.42% with respect to the reference solution.
568
tween different resolutions in SPH and the reference solution. The rela-
For Re = 100 (Fig. 6) we find good agreement in the velocity profiles be-
CR IP T
566
571
tive error in the position of the vortex for the highest resolution in SPH is u∗x = 0.24%, u∗y = 0.31% with respect to the reference solution.
572
change in the velocity profiles is sharper than before. Therefore, a good
573
agreement between SPH and reference solution is only obtained for the high-
574
est resolution. We see convergence of the vertical and horizontal velocity to
570
For Re = 1000 (Fig. 7) the effect of finer resolution is more obvious. The
AN US
569
577
the reference solution for higher resolution. The relative error in the position of the vortex for the highest resolution is u∗x = 0.60%, u∗y = 0.32% with
578
The 2D lid-driven cavity is investigated in the context of SPH in literature as
579
well [39, 44, 50, 58, 61, 62]. Compared to literature, the present results are
580
within a very good agreement with the reference solution and, therefore, we
581
conclude that the viscous and pressure force as well as the no-slip boundary
582
conditions for the velocity at the wall and numerical stability approaches
583
(DIDF and particle shifting) are well implemented.
584
At the beginning of the section we mentioned that an arbitrary bootstrap for
585
the solution of the Pressure Poisson Equation is needed. We note that the
586
velocity or pressure field is not affected by the artificial bootstrap as shown
ED
M
respect to the reference solution.
PT
576
CE
575
in Fig. 8. The reason is that the gradient of the pressure in the bottom left
588
corner is very small (Fig. 8b). Even if the pressure near the corner is different
589
from the specified bootstrap we obtain a good solution. This highlights that
590
the bootstrap is only necessary for algorithmic reasons during the solution
AC 587
37
AC
CE
PT
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 8: (a) Velocity and (b) pressure fields for Re = 100 with a resolution of L0 = 8.33µm.
38
591
of the linear equation system.
592
6.3. Diffusive transport of a scalar property
593
Diffusive transport of a scalar property φ
CR IP T
ACCEPTED MANUSCRIPT
Dφ = ∇ (α∇φ) , Dt
(67)
with the transport coefficient α appears in the energy equation (Eq. 9). We
595
validate diffusive transport for a general scalar property φ representing the
596
temperature. The discrete form of Eq. 67 was introduced in Eq. 19 and in
597
terms of the scalar variable φ is (∇ · α∇φ)a =
X mb b
ρb
AN US
594
(αa + αb )
~rab ∇a Wab · (φa − φb ) . 2 rab
(68)
In literature this discrete form for diffusive transport of a scalar variable
599
was validated, e.g. in [63, 64, 65]. Here we demonstrate convergence of
600
the present model. A very common test case is diffusion of a solute in one
601
direction. This is similar to heat conduction in a bar. For this test case a
ED
theoretical solution is available [66] 1 x0 − x x0 + x √ √ φ = φ0 erf + erf , 2 2 αt 2 αt
(69)
CE
PT
602
M
598
where x0 defines the region of concentrated solute φ0 at time t = 0. We
604
choose a transport coefficient α = 10−10 m2 /s, ρ = 1000kg/m3 , x0 = 0.5mm
AC
603
605
and
φ0 =
1 x < x0 0 else 39
,
(70)
ACCEPTED MANUSCRIPT
where the length of the domain is xmax = 2.5mm.
607
We show numerical convergence of the present model by increasing the reso-
608
lution for a smoothing length h = 2.05L0 . We consider 3 different resolutions
609
with L0 = 167µm, L0 = 83µm and L0 = 41.7µm. This corresponds to 15,
610
30 and 60 particles along x. We apply Neumann boundary condition at the
611
boundaries.
612
The profile of φ at time t = 1500s is shown in Fig. 9a. For the lowest resolu-
613
tion we see deviations from the theoretical profile especially near x = 0. The
614
profile for the highest resolution matches the theoretical profile very well.
615
The profile of φ for L0 = 41.7µm at different times in shown in Fig. 9b. For
616
all times the agreement is very good. Therefore, we conclude that the im-
617
plementation and the discrete form of the second order derivative are valid,
618
and therefore the heat transfer as well as mass transport in the bulk is valid.
M
AN US
CR IP T
606
ED
619
6.4. Buoyancy-driven cavity
621
Next we investigate non-isotherm single-phase flow to validate the coupling
622
of momentum and energy equation. Buoyancy-driven flow was investigated
623
in the context of SPH in a few articles [67, 68, 69]. Here we investigate a
624
buoyancy-driven cavity. In contrast to a lid-driven cavity, where the fluid is
625
driven by a velocity of the wall, the fluid is accelerated due to a temperature
626
difference between two walls. A schematic sketch is shown in Fig. 10.
627
The left wall is heated, the right wall is cooled and the top and bottom walls
628
are adiabatic. Due to the change in temperature, the fluid density decreases
629
with increasing its temperature. Because of a gravitational force, a circular
630
flow is initiated where denser fluid is accelerated downwards and lighter fluid
AC
CE
PT
620
40
ACCEPTED MANUSCRIPT
CR IP T
0.7
theoretical
φ [-]
0.6
L0 = 167µm
0.5
L0 = 83µm
0.4
L0 = 41.7 µm
0.3
AN US
0.2 0.1 0 0
0.5
1
1.5
2
2.5
position x [mm]
1
theoretical
0.8 0.7
t = 400s t = 1500s
ED
0.6 φ [-]
t = 100s
M
0.9
0.5 0.4 0.3
PT
0.2
t = 3000s
0.1
0
CE
0
0.5
1
1.5
2
2.5
position x [mm]
AC
Figure 9: Profile of φ along the x direction. (a) The profile for different resolutions at time t = 1500s and the theoretical profile are shown. (b) The profile of φ at different times with L0 = 41.7µm and the theoretical profiles are shown.
41
CR IP T
ACCEPTED MANUSCRIPT
Figure 10: Scheme for buoyancy-driven cavity example.
moves upwards. For small density variations the Boussinesq approximation
632
[70] is valid and the system can be regarded as incompressible since the pres-
633
sure is still independent of the density [20, 67]. Therefore, the only coupling
634
term between momentum balance and energy balance is the gravitational
635
acceleration ~g . According to Landau and Lifshitz [20], the difference in the
636
acceleration between two fluid portions with different temperature is
M
AN US
631
(71)
with the thermal expansion coefficient
β=−
1 ∂ρ , ρ ∂T p
(72)
CE
PT
637
ED
F~body = −β~g ∆T, ρ
in the range of 10−4 [1/K] for fluid considered here (i.e. polymers), and the
639
temperature difference
AC
638
∆T = T − T0 ,
(73)
640
with a reference temperature T0 . Note that the hydro-static part of the pres-
641
sure is neglected here because it should be negligible per definition. There42
ACCEPTED MANUSCRIPT
fore, we directly use Eq. 71 as an effective external force in the momentum
643
balance Eq. 8.
644
For natural convection the dimensionless numbers are the Prandtl number
Pr = 645
ν , k
the Rayleigh number
646
and the Gay-Lussac number
AN US
gβL3 ∆T , Ra = νk
CR IP T
642
Ga = β∆T.
(74)
(75)
(76)
The characteristic length L is the size of the cavity. ν and k are the kine-
648
matic viscosity of the fluid and the thermal diffusivity, respectively. The
649
Gay-Lussac number Ga describes the level of density variations caused by
650
temperature [67]. For Ga → 0 the Boussinesq approximation holds [71]. In
651
our case the temperature difference is ∆T = Th − Tc with Th and Tc as the
652
temperature at the hot and cold walls, respectively (see Fig 10).
653
We consider a square cavity of length L under the influence of gravity ~g =
654
[0, −9.81]. In the direction of gravity the cavity is adiabatic. In the other
655
direction a constant temperature is applied at the boundary. On the left side
656
we set the hot temperature to Th = 274.15K and on the right side the cold
657
temperature to Tc = 273.15K. We apply no-slip boundary conditions for the
658
velocity at all boundaries. Initially the temperature in the system is homoge-
659
neous T0 = 273.15K. The dynamic viscosity of the fluid is µ = 0.01P as and
660
kg m the density is ρ = 10 m 3 . The thermal diffusion coefficient is k = 0.00141 s
AC
CE
PT
ED
M
647
2
43
ACCEPTED MANUSCRIPT
and the thermal expansion coefficient is β = 0.071.
662
We investigate the flow in the cavity for two different Rayleigh numbers.
663
The size of the cavity varies with the Rayleigh number and is L = 126.4mm
664
and L = 272.4mm for Ra = 103 and Ra = 104 , respectively. We vary the
665
resolution to show convergence.
666
We compare the profile of the velocity and temperature along the vertical
667
and horizontal axis of the cavity with reference simulations using OpenFOAM
668
[60]. After investigation of the different resolutions in OpenFOAM we choose
669
a 120×120 regular grid for all reference simulations.
670
As already discussed in the previous section we need to set a bootstrap for
671
the pressure to solve the LES. In this case it is more difficult to identify a
672
point with very small gradient in the pressure. Therefore, it was more stable
673
to fix the pressure along the bottom wall at p = 0P a. We used the DIDF
674
and particle shifting approach as for the lid-driven cavity. The dimensionless
675
velocities are u∗y =
676
and y ∗ =
677
and Tmax = Th .
678
First we investigate a buoyancy-driven cavity at Ra = 1000. We vary the
679
resolution from L0 = 4.2mm to L0 = 1.05mm. The vertical and horizontal
680
velocity profiles along the axis at x∗ = 0.5 and y ∗ = 0.5 are shown in Fig. 11.
681
The velocity profiles are almost independent of the resolution and in good
ED
and u∗x =
ux ux,max
at dimensionless positions x∗ =
The dimensionless temperature is T ∗ =
T −T0 Tmax −T0
x L
with T0 = Tc
CE
PT
y . L
uy uy,max
M
AN US
CR IP T
661
agreement with the reference solution. The temperature profile along x∗ is
683
shown in Fig. 13a. Even with higher resolution we see a small deviation from
684
the reference solution. The deviation is in a reasonable range compared to
685
the available literature [67, 68, 69].
AC 682
44
ACCEPTED MANUSCRIPT
Next, we investigate the buoyancy-driven cavity at Ra = 10000. We again
687
vary the resolution. The vertical and horizontal velocity profiles along the
688
axis are shown in Fig. 12. For higher Rayleigh numbers the buoyancy force
689
is stronger. We see more deviations of the velocity in SPH compared to the
690
reference solution for coarse resolutions. With increasing resolution we find
691
convergence and good agreement to the reference solution. The temperature
692
profile shown in Fig. 13b also converges and is in very good agreement with
693
the reference solution.
694
In Fig. 14 we demonstrate for Ra = 10000 that the velocity and tempera-
695
ture field is smooth. We don’t see an effect of the artificial bootstrap for the
696
pressure and, therefore, skip the plot of the pressure here. The temperature
697
distribution is also in good agreement with literature [67]. Therefore, we
698
conclude that the coupling of the momentum and energy equation is valid.
699
Throughout this work we only need buoyancy flow with low and moderate
700
Rayleigh numbers. Therefore, we don’t show results for higher Rayleigh num-
701
bers here. The application of the SPH method to larger Rayleigh numbers
702
is shown in [67, 69].
703
7. Multi-phase systems
704
In this section, we validate the model for multi-phase applications. We dis-
705
cuss the formulation of surface tension including gradients of the surface
AC
CE
PT
ED
M
AN US
CR IP T
686
706
tension tangential to the interface. We assume that the surface tension is
707
temperature-dependent and analyze the surface force in a static and dy-
708
namic case to highlight features of the discretization. The normal part of the
709
surface tension is investigated using a spherical droplet. The tangential part 45
CR IP T
ACCEPTED MANUSCRIPT
OpenFOAM SPH (30x30) SPH (60x60) SPH (90x90) SPH (120x120)
0.5
0
AN US
dimensionless velocity uy* [-]
1
-0.5
-1 0
0.2
0.4
0.6
0.8
1
dimensionless position x* [-]
M
0.8
ED
0.6 0.4 0.2
PT
dimensionless position y* [-]
1
OpenFOAM SPH (30x30) SPH (60x60) SPH (90x90) SPH (120x120)
0
CE
-1
-0.5
0
0.5
1
dimensionless velocity ux* [-]
AC
Figure 11: Velocity profile for Ra = 1000: (a) vertical velocity vs. x∗ at y ∗ = 0.5, and (b) horizontal velocity vs. y ∗ at x∗ = 0.5.
46
CR IP T
ACCEPTED MANUSCRIPT
OpenFOAM SPH (30x30) SPH (60x60) SPH (90x90) SPH (120x120)
0.5
0
AN US
dimensionless velocity uy* [-]
1
-0.5
-1 0
0.2
0.4
0.6
0.8
1
dimensionless position x* [-]
M
0.8
ED
0.6 0.4 0.2
PT
dimensionless position y* [-]
1
OpenFOAM SPH (30x30) SPH (60x60) SPH (90x90) SPH (120x120)
0
CE
-1
-0.5
0
0.5
1
dimensionless velocity ux* [-]
AC
Figure 12: Velocity profile for Ra = 10000: (a) vertical velocity vs. x∗ at y ∗ = 0.5, and
(b) horizontal velocity vs. y ∗ at x∗ = 0.5.
47
CR IP T
ACCEPTED MANUSCRIPT
OpenFOAM SPH (30x30) SPH (60x60) SPH (90x90) SPH (120x120)
0.8 0.6 0.4
AN US
dimensionless temperature T* [-]
1
0.2 0 0
0.2
0.4
0.6
0.8
1
dimensionless position x* [-]
M
OpenFOAM SPH (30x30) SPH (60x60) SPH (90x90) SPH (120x120)
ED
0.8 0.6 0.4
PT
dimensionless temperature T* [-]
1
0.2
CE
0
0
0.2
0.4
0.6
0.8
1
dimensionless position x* [-]
AC
Figure 13: Temperature profile vs. x∗ at y ∗ = 0.5 for (a) Ra = 1000 and (b) Ra = 10000.
48
ACCEPTED MANUSCRIPT
of the surface tension is investigated using a planar interface.
711
Finally, we consider the complete model to investigate thermocapillary flow.
712
We consider the migration of a spherical droplet [72, 18, 17] and investi-
713
gate the stationary migration velocity. We show numerical convergence and
714
compare the SPH method to reference simulations taken from literature.
715
In a multi-phase system we get a non-isotropic pressure contribution, due to
716
the interface, in the momentum balance. This contribution is called capillary
717
stress in general or surface tension in an immiscible system [73]. The form of
718
the capillary stress tensor in the momentum balance accounts for a normal
719
and tangential stresses. Therefore, we may split the capillary stress into a
720
normal and tangential part. The tangential part is known as the Marangoni
721
force and accounts for gradients in the surface energy along the interface.
722
In this work we use the continuum surface force (CSF) model [31, 32] to
723
discretise the capillary stress. This model is well validated in literature for
724
SPH [12, 26, 32, 65, 74, 75]. In this section we focus on the validation of the
725
implementation. We use the corrected SPH approach and particle shifting
726
every time when capillary stress is considered.
727
As introduced in section 4.2 we use a color function to calculate the normal
AN US
M
ED
PT
at the interface between two phases. The color function is sharp 0 if fluid 1 Ca = , 1 if fluid 2
AC
CE
728
CR IP T
710
729
where each particle of a fluid has a constant color Ca .
49
(77)
AC
CE
PT
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 14: (a) Velocity and (b) temperature distributions in a buoyancy–driven cavity for Ra = 10000.
50
ACCEPTED MANUSCRIPT
7.1. Static capillary stress normal to interface
731
We first investigate the normal part of the capillary stress in a static test
732
case. We consider a spherical droplet in a second fluid at rest and compare
733
the pressure jump over the interface with the Young-Laplace equation [76].
734
In 2D the Young-Laplace equation for a pressure jump is ∆pY L =
σ , Rd
CR IP T
730
(78)
with the radius of the droplet Rd . We consider a droplet with a radius
736
Rd = 0.2m inside a quadratic box with size L = 1m filled with a second
737
fluid. Both fluids have the same physical properties. The density, dynamic
738
kg mN viscosity and surface tension are ρ = 1000 m 3 , η = 0.01P as and σ = 72.75 m ,
739
respectively. The resolution of the box is 120×120 particles where 24 particles
740
per droplet radius is used. Initially the particles are in regular Cartesian
741
order. At the boundary we apply no-slip conditions for the velocity. At the
742
top of the box we use Dirichlet boundary conditions and at the other three
743
boundaries Neumann boundary conditions for the pressure. The theoretical
744
pressure jump in the present case is ∆p = 0.36375P a.
745
At the beginning the order of the particles of the droplet isn’t smooth because
746
of the Cartesian order. We wait until an equilibrium configuration of the
747
particles is reached and then investigate the pressure jump. Fig. 15a shows a
748
3-dimensional graph of the pressure distribution in the box at an equilibrium
AC
CE
PT
ED
M
AN US
735
749
configuration.
750
We see that the pressure jump isn’t smooth. At the interface we observe
751
fluctuations where the pressure is higher or lower than the pressure in the
752
bulk phases. This is more obvious in the pressure profile along the radius 51
AN US
CR IP T
ACCEPTED MANUSCRIPT
0.4 0.2
M
Young-Laplace
SPH (120x120)
ED
0.6
PT
pressure p [Pa]
0.8
CE
0
-0.4
-0.2
0
0.2
0.4
radius r [m]
AC
Figure 15: (a) 3D graph and (b) plot of the pressure distribution of a single droplet in rest in a square box with surface tension σ = 72.75 mN m .
52
ACCEPTED MANUSCRIPT
of the droplet (Fig. 15b). These pressure fluctuations are also observed and
754
discussed in literature [65, 77]. The reason for the pressure fluctuations are
755
the spurious currents at the interface caused by unfavorable energetic con-
756
figuration of the particles [78]. The pressure in the bulk phase is constant.
757
Therefore, we quantify the pressure jump across the interface using the con-
758
stant pressure in the bulk phases. The pressure jump in the simulation is
759
∆pSP H = 0.3601P a. The relative error is ∆pY L − ∆pSP H · 100%, ∆pY L
AN US
=
CR IP T
753
(79)
of = 0.99%. This error is in the same order as in literature [32, 65].
761
7.2. Dynamic capillary stress normal to interface
762
Next we demonstrate that the dynamics of the fluid due to capillary stress
763
is valid. A common test case is an oscillating droplet, where an analytic
764
solution of the oscillation period is only available in the limit of a large
765
density ratio [79]. It was shown in literature that the CSF model using a
766
sharp color function accurately reproduces the oscillation period of a droplet
767
for different surface tensions [12, 32, 65, 77]. Numerical reference data for a
768
density ratio of 1 is also available in literature [65].
769
The test case consists of a spherical droplet with radius R = 0.2m in the
770
center of a square box with length L = 1m. At the boundaries we use
771
no-slip conditions for the velocity. The pressure is fixed at the top of the
772
box and Neumann boundary conditions for the pressure are applied at the
773
other boundaries. In both phases, the density, dynamic viscosity and surface
774
tension are ρ = 1kg/m3 , η = 0.05P as and σ = 1N/m, respectively. Initially,
AC
CE
PT
ED
M
760
53
ACCEPTED MANUSCRIPT
t [s]
[65]
≈ 0.4
sharp
0.409
CR IP T
color function
Table 1: Oscillation periods of droplet with density ratio 1. Comparison with estimated period from Adami et al. [65].
776
we impose a divergence-free velocity field [32, 65] x y2 r ux = U0 1− exp − , r0 r0 r r0 and y uy = −U0 r0
AN US
775
x2 r 1− exp − , r0 r r0
(80)
(81)
with u0 = 10m/s and r0 = 0.05m. The amplitude of the oscillation is
778
damped by the viscous force but for small viscosity the oscillation period is
779
not influenced.
780
The results are shown in Tab. 1. We estimate the period as the mean of x and
781
y position of the droplet, as shown in Fig. 7a. The mean is approximately
782
t = 0.4s. We are very close to the period of the reference literature. In
783
addition we compare the contour of the oscillating droplet at different times
784
with a contour from [65]. The particle distribution of the droplet is shown in
785
Fig. 16 where only the particles of the droplet are shown. The surrounding
786
particles are not shown. The particles in the upper half are taken from [65].
787
The particles in the bottom half correspond to the present simulation. The
788
contours match with small deviation for larger times, however, the particle
789
distributions are very similar even for longer periods.
790
Finally, we show the maximum position of the droplet in x and y direction
AC
CE
PT
ED
M
777
54
CR IP T
ACCEPTED MANUSCRIPT
Figure 16: Particle distribution of an oscillating droplet. Top half of every image is taken
from [65]. The bottom half of every image are present results. The surface tension is
AN US
N σ = 1m . From left to rigth, time is t = 0s, t = 08s, t = 16s, and t = 26s, respectively.
over time in Fig. 17, compared to results from [65]. There is a difference
792
between the amplitude of the present results and the literature. The reason
793
for the difference in the amplitude is that in the reference [65] the variation
794
of position is calculated in a different way than in the present simulations.
795
Adami et al. take the average of the mean position in a quarter of the
796
droplet and plot it over time. Here we take the maximum position of a
797
particle in x and y direction. Therefore, the variation is larger in the present
798
simulation but both kind of analysis indicate the frequency. We see that
799
the position oscillates around the average position x = 0.7m and y = 0.7m,
800
respectively. In comparison to the reference, the oscillation profiles in the
801
present simulations are very smooth. A very small drift of the position is
802
seen for t → 1s because the droplet moves a little bit out of the center of
803
the box. This is caused by small fluctuations of the pressure due to spurious
AC
CE
PT
ED
M
791
804
currents.
805
7.3. Static capillary stress tangential to interface
806
In this section we investigate the Marangoni force in a static test case. We
807
consider two layered fluids in an infinitely long channel of width 5.76mm in 55
ACCEPTED MANUSCRIPT
0.74
position [m]
0.72 0.71 0.7 0.69 0.68 0.67 0.66 0.2
0.4
0.6
0.8
1
AN US
0
CR IP T
X position (SPH) Y position (SPH) X position (Adami et al.) Y position (Adami et al.)
0.73
time t [s]
N Figure 17: Position variations over time for a droplet with R = 0.2m and σ = 1 m .
2D. The interface between the fluids is planar. We apply a linear temperature
809
profile along the channel with a slope of ∇T = 200K/m. We assume that
810
the surface tension coefficient depends linear on temperature
M
808
811
ED
σ = σ0 + σT (T − T0 ) ,
(82)
with the surface tension coefficient σ0 at T0 and the coefficient σT = −2mN/(mK). With these parameters the theoretical Marangoni force is ∇S σ = 0.4N/m2 .
813
The only force in the system is the thermocapillary force tangential to the
814
planar interface. We investigate the Marangoni force for different resolutions
815
to show numerical convergence. We investigate the resolutions L0 = 180µm,
816
L0 = 90µm and L0 = 45µm.
817
Because SPH is a smoothed discretization method the Marangoni force is
818
smoothed out around the interface between the fluids. In Fig. 18a we see the
819
profile of the Marangoni force perpendicular to the interface in the channel.
820
We see that the Marangoni force is wider for low resolution and narrower for
AC
CE
PT
812
56
ACCEPTED MANUSCRIPT
CR IP T
3
L0 = 180 µm
L0 = 90 µm L0 = 45 µm
2 1.5 1 0.5 0 2
2.2
2.4
2.6
AN US
force [kN/m3]
2.5
2.8
3
3.2
3.4
3.6
position [mm]
0.401
0.2
0.4
PT
0.3995
analytic
0.1
0.05
0.399
CE
0
50
0.15
error
ED
force [N]
0.4005
error [%]
M
SPH
100
150
0 200
L0 [µm]
AC
Figure 18: (a) Profile of Marangoni force in normal direction to the interface using a resolution of L0 = 180µm, L0 = 90µm and L0 = 45µm. (b) Integrated force and relative
error of the Marangoni force for different resolutions.
57
ACCEPTED MANUSCRIPT
higher resolution. In the limit of infinite resolution the profile tends to a Dirac
822
function with the magnitude of the Marangoni force. The calculated integral
823
of the smoothed Marangoni force of the simulation is shown in Fig. 18b.
824
The Marangoni force in SPH is only slightly different from the theoretical
825
value. The relative error is = 0.014% based on the theoretical value. The
826
calculated Marangoni force is independent of the resolution because we used
827
the corrected SPH approach with regular order of the particles and, therefore,
828
the force is as accurate as possible. In [17] a similar case was presented
829
where the relative error was also constant, but with = 3.65% very large.
830
The reason may be a different application of the corrected SPH approach in
831
[17]. One major difference is that we calculate the shepard kernel using only
832
particles with a normal larger than 0.01/h. We conclude that the Marangoni
833
force is very accurate in the present model.
834
7.4. Thermocapillary droplet migration
835
Thermocapillary motion was studied for several decades theoretically [76, 80],
836
numerically [18, 72] and experimentally [81, 82]. An overview about the work
837
of the 20th century is found in [83] and the references therein. In this section
838
we combine the normal and tangential surface force and investigate thermo-
839
capillary migration of a droplet in 2D and 3D. In thermocapillary flow, a
840
droplet moves because of gradients of the surface tension tangential to the
841
interface, caused by a gradient of the temperature, even in the absence of
842
gravity. For sufficiently low Reynolds numbers, the shape of the droplet is
843
spherical and constant during motion. Here, we first show numerical conver-
844
gence with increasing resolution, investigate the effect of a bounded domain
845
and compare our results to [17] and [18]. Young et al. [76] derived an ana-
AC
CE
PT
ED
M
AN US
CR IP T
821
58
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 19: Schematic setup of thermocapillary migration of a droplet.
lytic solution for the steady-state migration velocity of a droplet in 3D. We
847
compare the present results for 3D migration velocity to the analytic solution.
848
Initially, a spherical droplet (phase 1) with the radius of R is located in
849
the center of a square box with the length of L. In x direction we ap-
850
ply Neumann boundary conditions for the temperature and free slip condi-
851
tions for the velocity. In y direction we apply Dirichlet boundary conditions
852
for temperature and no-slip conditions for the velocity. The density, vis-
853
cosity and thermal diffusion coefficient of the droplet are ρ1 = 250kg/m3 ,
854
η1 = 0.012P as and λ1 = 1.2 · 10−6 W/(mK). The density, viscosity and
ED
PT
CE
thermal diffusion coefficient of the continuous phase are ρ2 = 500kg/m3 ,
AC
855
M
846
856
η2 = 0.024P as and λ2 = 2.4 · 10−6 W/(mK), respectively. The specific heat
857
capacities are cp,1 = 0.5 · 10−4 J/(kgK) and cp,2 = 10−4 J/(kgK). The surface
858
tension is σ = σ0 + σT (T − T0 ) , 59
(83)
ACCEPTED MANUSCRIPT
with σ0 = 10mN/m and σT = 2mN/(mK). The temperature at the bottom
860
of the box is the reference temperature T1 = T0 = 290K. The temperature
861
gradient in y direction is ∇T = 200K/m. The radius of the droplet is
862
R = 1.44mm. The setup of the test case is shown in Fig. 19. Note that
863
we don’t need a special treatment of the interface, e.g. we don’t use an
864
additional pressure contribution at the interface.
865
The dimensionless numbers to characterize thermocapillary migration are the
866
Reynolds number
AN US
ρ2 RUR = 0.72, η2
Re = 867
Marangoni number Ma =
ρ2 cp,2 RUR = 0.72, λ2
and Capillary number
M
868
Ca =
(85)
(86)
PT
σT |∇T |R m = 0.024 . η2 s
(87)
t · UR , R
(88)
U . UR
(89)
CE
The dimensionless time is t∗ =
and the dimensionless velocity is
AC
871
(84)
with the characteristic velocity of UR =
870
η2 UR = 0.0576, σ0
ED
869
CR IP T
859
U∗ =
872
First, we investigate numerical convergence of the migration velocity in 2D.
873
We choose a box length L = 5.76mm which corresponds to a ratio L/R = 4.
874
We vary the resolution from L0 = 180µm to L0 = 45µm. The dimensionless 60
ACCEPTED MANUSCRIPT
0.16 0.14 0.12
L0 = 180 µm L0 = 90 µm
0.08
CR IP T
U* [-]
0.1
L0 = 45 µm
0.06
L0 = 22.5 µm
0.04
Ma et al.
0.02
Tong et al.
0 0
0.5
1
1.5
2
AN US
t* [-]
Figure 20: Droplet migration using different resolutions compared to solution from [18] and [17]. The box length to radius ratio is L/R = 4.
migration velocity U ∗ over dimensionless time t∗ is shown in Fig. 20. Repre-
876
sentative plots of the acceleration due to surface forces, velocity, temperature
877
and pressure are shown in Figs. 21 and 22 for L0 = 45µm. We compare the
878
migration velocity to the results of Tong et al. [17] and Ma et al. [18] where
879
a resolution similar to L0 = 90µm is used.
880
In Fig. 20 we find two major results. Firstly, we see that the migration
881
velocity of the present model is very similar to the migration velocity of Ma
882
& Bothe [18]. Compared to Tong et al. [17] the present results are more
883
accurate. The reason is that the error in the Marangoni force is lower in the
884
present model. Secondly, we see that with increasing or decreasing resolution
885
we over- and underestimate the data of Ma & Bothe. The present model
AC
CE
PT
ED
M
875
886
converges with increasing resolution indicating that the results presented by
887
Tong et al. [17] were not yet converged.
888
If analytic solutions of migration velocity are available, the box size is always
889
assumed infinitely large. This means that the wall effects do not disturb the 61
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 21: Plot of droplet migration using L0 = 45µm particles at t = 0.6s. The box length to radius ratio is L/R = 4. On the left side the acceleration due to surface forces
AC
CE
PT
ED
M
is plotted and on the right side the magnitude of the velocity of shown.
Figure 22: Plot of droplet migration using L0 = 45µm particles at t = 0.6s. The box length to radius ratio is L/R = 4. On the left side the temperature is plotted and on the
right side the pressure is shown.
62
ACCEPTED MANUSCRIPT
0.2 0.18 0.16 0.12 0.1
CR IP T
U* [-]
0.14
L/R = 4
0.08
L/R = 8
0.06 0.04
L/R = 16
0.02
L/R = 32
0 0
0.5
1
1.5
2
AN US
t* [-]
Figure 23: Influence of box size on the migration velocity of a droplet in 2D. The droplet radius is resolved with 16 particles.
flow field. In numerical simulations we always have a bounded domain due to
891
the limited computational resources. Next we investigate the influence of the
892
bounded domain on the migration velocity of the droplet in 2D. We consider
893
the same setup and properties of the fluids as before. Then we increase the
894
box size from L = 5.76mm to L = 46.08mm. The resolution is L0 = 90µm
895
or 16 particles per droplet radius. Initially the droplet is placed in the center
896
of the box.
897
The migration velocity U ∗ over time t∗ is shown in Fig. 23. We see that
898
with increasing box size the steady-state migration velocity increases and
899
converges to a limit of U ∗ ≈ 0.18. For L/R ratios larger than 16 the migra-
900
tion velocity is converged and we may assume that the flow field is approx-
AC
CE
PT
ED
M
890
901
imately unaffected by the boundaries. A similar value but for the 3D case
902
was presented by Ma & Bothe [18] where the steady-state migration velocity
903
is constant for L/R ratios larger than 16.
904
Finally, we consider thermocapillary migration of a spherical droplet in 3D. 63
ACCEPTED MANUSCRIPT
905
For the steady-state migration velocity of a spherical droplet in 3D, Young
906
et al. [76] derived a theoretical value of 2 , η1 λ1 2 + 3 λ2 η2
(90)
CR IP T
U∗ = 2+
with the ratio of the thermal conductivity and dynamic viscosity of the
908
droplet and the continuous phase, where they assumed that the domain is
909
unbounded. But in the simulation we have a bounded domain and therefore
910
we will reach the theoretical value only in the limit of a large L/R ratio.
911
In the 3D model we consider a spherical droplet in a closed box. We use slip
912
boundary condition for velocity in x and z directions and no-slip boundary
913
conditions in y direction. Similarly we use Neumann boundary condition for
914
the temperature in x and z directions and Dirichlet boundary conditions in
915
y direction. The properties of the fluids are the same as in the previous sim-
916
ulations. Initially a spherical droplet is centered in the box. To reduce the
917
numerical effort we use 16 particles per droplet radius with L0 = 90µm. The
918
smoothing length is reduced to h = 1.55L0 because there are more particles
919
in the vicinity in 3D than in 2D and, therefore, the SPH approximation is
920
already good for lower smoothing lengths.
921
A representative snapshot of the 3D simulation is shown in Fig. 24. Only
922
the droplet without the surrounding fluid is shown. The temperature and
923
velocity distributions at the surface of the droplet are indicated by the color
AC
CE
PT
ED
M
AN US
907
924
and the arrows. The velocity at the bottom is directed into the droplet while
925
the velocity at the top is directed out of the droplet. The temperature profile
926
is close to linear.
927
In Fig. 25 we demonstrate that we converge towards the theoretical value 64
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 24: Plot of droplet migration using L0 = 90µm particles at t = 0.6s. The box length to radius ratio is L/R = 4. The color indicates temperature and the arrows indicate the
M
direction of velocity.
by increasing the L/R ratio. We see that the steady-state migration velocity
929
increases with increasing L/R ratio. We underestimate the theoretical value
930
because the domain is still bounded but we are very close for L/R = 8. As
931
we have seen at the beginning of this section the influence of the bounded
932
domain is only negligible for L/R ratio larger than 16.
933
We also compare the results of the present simulations with results of Ma
934
et al. [18]. The present results are slightly closer to the theoretical value
935
than the results of Ma & Bothe. We have not performed simulations for very
936
large L/R ratios because the numerical effort is very high. In the results of
937
Ma & Bothe it seems that there is a constant offset in the migration velocity
938
because the migration velocity is very similar for L/R = 8 and L/R = 16.
939
We expect similar behavior of the present model because we use a very sim-
940
ilar formulation of the CSF model. Therefore, we conclude that the present
AC
CE
PT
ED
928
65
ACCEPTED MANUSCRIPT
0.25 0.2
0.1
CR IP T
U* [-]
0.15
SPH
0.05
Ma et al. Young et al.
0 0
5
10
20
AN US
L/R ratio [-]
15
Figure 25: Influence of box length to droplet diameter ratio for 3D simulation where the droplet diameter is resolved with 32 particles. Comparison with analytic migration velocity [76, 84]
.
model predicts thermocapillary flow very accurate compared to available lit-
942
erature.
943
8. Conclusion
944
We presented an incompressible SPH model for thermo-capillary flow driven
945
by gradients of the surface tension. Rigorous validation of the proposed
946
model shows accurate implementation and convergence of the results.
947
In the first part of the article, we validated the model for single-phase prob-
948
lems. First, we demonstrated the effects of corrected SPH, DIDF and par-
AC
CE
PT
ED
M
941
949
ticle shifting approach on particle distribution using a Taylor-Green vortex.
950
It is shown that only particle shifting enables a homogeneous particle dis-
951
tribution. Then, lid-driven cavity is studied for different Reynolds numbers.
952
The results are compared to reference solutions taken from OpenFOAM. A 66
ACCEPTED MANUSCRIPT
convergence study shows that with increasing resolution the SPH model con-
954
verges to the OpenFOAM results. Next, we validated diffusive mass trans-
955
port, compared it to analytic solution and found very accurate agreement.
956
Afterwards, buoyancy-driven flow is investigated for different Rayleigh num-
957
bers and compared to reference solution of OpenFOAM. Again, a convergence
958
study demonstrates that the results of SPH match the results of OpenFOAM
959
very well. Therefore, we conclude that the model was correctly implemented.
960
In the second part of the article, we investigated multi-phase systems. First,
961
we studied surface tension in detail using static and dynamic test cases.
962
The normal component of the surface tension is validated using the Young-
963
Laplace law at steady-state in 2D. The analytic law is captured very well.
964
Deviations of the pressure at the interface are a result of the CSF model and
965
where previously found in other models as well. The error is approx. 1 %.
966
We studied the oscillation of a 2D droplet and compared the results with
967
reference solution of Adami et al. from literature. We analyzed the time
968
evolution of the oscillation and its deviation to the reference. The oscillation
969
period as well as the contour of the 2D droplet match very well.
970
The Marangoni force is studied separately for a planar interface with a tan-
971
gential, linear temperature gradient. A convergence study is presented and
972
we found that the error is around 0.1 %, and independent of the numerical
973
resolution because of the use of corrected SPH approach for the Marangoni
CE
PT
ED
M
AN US
CR IP T
953
force. This is way better than in Tong et al.
975
Finally, we studied thermo-capillary droplet migration in a temperature field
976
with temperature-dependent surface tension coefficient. First, we present a
977
convergence study for 2D droplet and compare the results to reference of SPH
AC 974
67
ACCEPTED MANUSCRIPT
(Tong et al.) and Finite Volume (Ma et al.). We show that the previously
979
presented results of Tong et al. agree with FVM but have not converged in
980
SPH. Then we study the same case in 3D and compare the results to the an-
981
alytic solution of Young et al. The results of SPH show convergence and are
982
slightly more accurate than the results of FVM taken from Ma et al. There-
983
fore, we conclude that the proposed model for thermo-capillary flow captures
984
the studied cases very well. Compared to the previous work of Tong et al.,
985
this study completes with a convergence study and more accurate results.
986
9. Acknowledgement
987
M. Hopp-Hirschler and U. Nieken acknowledge funding from the German Re-
988
search Foundation within the Collaborative Research Center SFB716. M.S.
989
Shadloo acknowledges funding provided by DAAD programme ”Research
990
Stays for University Academics and Scientists, 2017” through funding ID:
991
57314019 during his visit from Stuttgart University. The access to the French
992
CRIANN (Centre Rgional Informatique et d’Applications Numriques de Nor-
993
mandie) resources, under allocation 2017002 is also acknowledged. This work
994
was granted access to HPC resources of IDRIS under the allocation 2017-
995
100752 made by GENCI (A0022A10103).
996
10. References
AC
CE
PT
ED
M
AN US
CR IP T
978
997
998
999
References [1] L. B. Lucy, A numerical approach to the testing of the fission hypothesis, The Astronomical Journal 82 (1977) 1013–1024.
68
ACCEPTED MANUSCRIPT
[2] R. A. Gingold, J. J. Monaghan, Smoothed Particle Hydrodynamics:
1001
theory and application to non-spherical stars, Mon. Not. R. astr. Soc
1002
181 (1977) 375–389.
1004
1005
1006
1007
1008
[3] S. J. Cummins, M. Rudman, An SPH Projection Method, J. Comp. Phys. 152 (1999) 584–607.
[4] J. J. Monaghan, Simulating free surface flows with sph, Journal of Computational Physics 110 (1994) 399–406.
AN US
1003
CR IP T
1000
[5] J. P. Morris, P. J. Fox, Y. Zhu, Modeling low reynolds number incompressible flows using sph, J. Comp. Phys. 136 (1997) 214–226. [6] A. Colagrossi, M. Landrini, Numerical simulation of interfacial flows by
1010
smoothed particle hydrodynamics, J. Comp. Phys. 191 (2003) 448–475.
1012
[7] J. J. Monaghan,
Smoothed Particle Hydrodynamics,
Reports on
Progress in Physics 68, Issue 8 (2005) 1703–1759.
ED
1011
M
1009
[8] S. Shao, E. Y. M. Lo, Incompressible SPH method for simulating New-
1014
tonian and non-Newtonian flows with a free surface, Advances in Water
1015
Resources 26 (2003) 787–800. [9] M. Shadloo, G. Oger, D. Le Touz´e, Smoothed Particle Hydrodynamics method for fluid flows, towards industrial applications: Motivations,
AC
1017
CE
1016
PT
1013
1018
current state, and challenges, Computers & Fluids 136 (2016) 11 – 34.
1019
[10] M. B. Liu, G. R. Liu, Smoothed Particle Hydrodynamics (SPH): an
1020
Overview and Recent Developments, Arch. Comput. Method. Eng. 17
1021
(2010) 25–76. 69
ACCEPTED MANUSCRIPT
[11] N. Tofighi, M. Yildiz, Numerical simulation of single droplet dynam-
1023
ics in three-phase flows using isph, Computers & Mathematics with
1024
Applications 66 (2013) 525–536.
CR IP T
1022
1025
[12] X. Y. Hu, N. A. Adams, A multi-phase sph method for macroscopic
1026
and mesoscopic flows, Journal of Computational Physics 213 (2006)
1027
844–861.
[13] D. Le Touz´e, G. Oger, B. Alessandrini, Smoothed Particle Hydrodynam-
1029
ics simulation of fast ship flows, in: Proc. 27th Symposium on Naval
1030
Hydrodynamics.
AN US
1028
[14] X. Y. Hu, N. A. Adams, A constant-density approach for incompressible
1032
multi-phase sph, Journal of Computational Physics 228 (2009) 2082–
1033
2091.
M
1031
[15] A. Zainali, N. Tofighi, M. S. Shadloo, M. Yildiz, Numerical investigation
1035
of newtonian and non-newtonian multiphase flows using isph method,
1036
Computer Methods in Applied Mechanics and Engineering 254 (2013)
1037
99–113.
1039
PT
[16] S. Adami, X. Y. Hu, N. A. Adams, A conservative sph method for surfac-
CE
1038
ED
1034
tant dynamics, Journal of Computational Physics 229 (2010) 1909–1926.
[17] M. Tong, D. Browne, An incompressible multi-phase smoothed particle
1041
hydrodynamics (sph) method for modelling thermocapillary flow, Int.
1042
J. Heat Mass Trans. 73 (2014) 284–292.
AC
1040
1043
[18] C. Ma, D. Bothe, Direct numerical simulation of thermocapillary flow 70
ACCEPTED MANUSCRIPT
1044
based on the volume of fluid method, International Journal of Multi-
1045
phase Flow 37 (2011) 1045 – 1058. [19] P. Espanol, C. Thieulot, Microscopic derivation of hydrodynamic equa-
1047
tions for phase-separating fluid mixtures,
1048
Physics 118 (2003).
CR IP T
1046
The Journal of Chemical
[20] L. D. Landau, E. M. Lifshitz, Course of Theoretical Physics: Fluid
1050
Dynamics, volume 6, Elsevier Ltd., Butterwoth Heinemann, 2nd edition,
1051
1987.
1053
1054
1055
[21] M. Hopp-Hirschler, Modeling of Porou Polymer Membrane Formation, Dissertation d93, shaker verlag, Universit¨at Stuttgart, 2017. [22] J. J. Monaghan, Smoothed Particle Hydrodynamics, Annu. Rev. Stron.
M
1052
AN US
1049
Astrphys. 30 (1992) 543–574.
[23] H. Wendland, Piecewise polynominal, positive definite and compactly
1057
supported radial functions of minimal degree, Adv. Comput. Math. 4
1058
(1995) 389–396.
1060
PT
Press, 2012.
[25] J. Bonet, T.-S. L. Lok, Variational and momentum preservation as-
AC
1061
[24] D. Violeau, Fluid Mechanics and the SPH Method, Oxford University
CE
1059
ED
1056
1062
pects of Smooth Particle Hydrodynamic formulations, Comput. Method.
1063
Appl. M. 180 (1999) 97–115.
1064
[26] N. Grenier, M. Antuono, A. Colagrossi, D. Le Touz´e, B. Alessandrini, An
71
ACCEPTED MANUSCRIPT
1065
Hamiltonian interface SPH formulation for multi-fluid and free surface
1066
flows, Journal of Computational Physics 228 (2009) 8380–8393.
1068
[27] L. Brookshaw, A method of calculating radiative heat diffusion in par-
CR IP T
1067
ticle simulations, Proc. ASA 6(2) (1985) 207–210.
[28] R. Fatehi, M. T. Manzari, Error estimation in smoothed particle hy-
1070
drodynamics and a new scheme for second derivatives, Computers and
1071
Mathematics with Applications 61 (2011) 482–498.
AN US
1069
1072
[29] K. Szewc, J. Pozorski, J.-P. Minier, Analysis of the incompressibility
1073
constraint in the smoothed particle hydrodynamics method, Int. J. Nu-
1074
mer. Meth. Eng. 92 (2012) 343–369.
[30] J. Morris, Analysis of smoothed particle hydrodynamics with applica-
1076
tions, Ph.D. thesis, Department of Mathematics, Monash University,
1077
1996.
1080
1081
ED
[32] J. P. Morris, Simulating surface tension with Smoothed Particle Hydrodynamics, Int. J. Numer. Meth. Fluids 33 (2000) 333–353.
[33] P. G. Tait, Report on some of the physical properties of fresh water and
AC
1082
Modeling Surface Tension, J. Comp. Phys. 100 (1992) 335–354.
PT
1079
[31] J. U. Brackbill, D. B. Kothe, C. Zemach, A Continuum Method for
CE
1078
M
1075
1083
sea water, Rept. Sci. Results Voy., H.M.S. Challanger, Phys. Chem. 2
1084
(1888) 1–76.
1085
[34] H. Helmholtz,
u ¨ber integrale der hydrodynamicschen gleichungen,
72
ACCEPTED MANUSCRIPT
1086
welcher der wirbelbewegungen entsprechen, Journal fr die reine und
1087
angewandte Mathematik 55 (1858) 25–55.
1089
[35] X. Y. Hu, N. A. Adams, An incompressible multi-phase sph method, J. Comp. Phys. 227 (2007) 264–278.
CR IP T
1088
[36] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschel-
1091
man, L. Dalcin, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley,
1092
L. C. McInnes, K. Rupp, B. F. Smith, S. Zampini, H. Zhang, H. Zhang,
1093
PETSc Web page, http://www.mcs.anl.gov/petsc, 2016. Version 3.3.
1094
[37] R. Falgout, A. Cleary, J. Jones, E. Chow, V. Henson, C. Baldwin,
1095
P. Brown, P. Vassilevski, U. M. Yang, Hypre Web page, http://acts.
1096
nersc.gov/hypre, 2016.
M
1098
[38] F. Keller, Simulation of the Morphogenesis of Open-porous Materials, Ph.D. thesis, University of Stuttgart, 2015.
ED
1097
AN US
1090
[39] M. Yildiz, R. A. Rook, A. Suleman, Sph with the multiple boundary
1100
tangent method, International Journal for Numerical Methods in Engi-
1101
neering 77 (2009) 1416–1438.
1103
Unified semi-analytical wall boundary conditions for inviscid laminar or turbulent flows in the meshless sph method, Int. J. Numer. Meth. Fluids
AC
1104
[40] M. Ferrand, D. R. Laurence, B. D. Rogers, D. Violeau, C. Kassiotis,
CE
1102
PT
1099
1105
71 (2013) 446–472.
1106
[41] S. Khorasanizade, J. M. M. Sousa, An innovative open boundary treat-
1107
ment for incompressible sph, International Journal for Numerical Meth-
1108
ods in Fluids 80 (2016) 161–180. 73
ACCEPTED MANUSCRIPT
[42] M. Hirschler, P. Kunz, M. Huber, F. Hahn, U. Nieken, Open boundary
1110
conditions for ISPH and their application to micro-flow, Journal of
1111
Computational Physics 307 (2016) 614 – 633.
CR IP T
1109
[43] M. S. Shadloo, A. Zainali, S. H. Sadek, M. Yildiz, Improved incom-
1113
pressible smoothed particle hydrodynamics method for simulating flow
1114
around bluff bodies, Computer methods in applied mechanics and engi-
1115
neering 200 (2011) 1008–1020.
AN US
1112
1116
[44] R. Xu, P. Stansby, D. Laurence, Accuracy and stability in incompressible
1117
SPH (ISPH) based on the projection method and a new approach, J.
1118
Comp. Phys. 228 (2009) 6703–6725.
[45] M. S. Shadloo, A. Zainali, M. Yildiz, Simulation of single mode rayleigh–
1120
taylor instability by sph method, Computational Mechanics 51 (2012)
1121
699–715.
ED
M
1119
[46] M. Shadloo, A. Rahmat, M. Yildiz, A smoothed particle hydrodynamics
1123
study on the electrohydrodynamic deformation of a droplet suspended
1124
in a neutrally buoyant newtonian fluid, Computational Mechanics 52
1125
(2013) 693–707.
CE
PT
1122
[47] S. J. Lind, R. Xu, P. K. Stansby, B. D. Rogers, Incompressible Smoothed
1127
Particle Hydrodynamics for free-surface flows: A generalised diffusion-
AC
1126
1128
based algorithm for stability and validations for impulsive flows and
1129
propagating waves, J. Comp. Phys. 231 (2012) 1499–1523.
1130
1131
[48] D. Shepard, A two dimensional function for irregulary space data, Proceedings of ACM national conference (1968) 517–524. 74
ACCEPTED MANUSCRIPT
[49] X. Hu, personal communication, SPHERIC 2015, Parma, 2015.
1133
[50] U. Ghia, K. N. Ghia, C. T. Shin, High–re solutions for incompressible
1134
flow using thre navier–stokes equations and a multigrid method, J.
1135
Comp. Phys. 48 (1982) 387–411.
1137
1138
1139
[51] R. Schreiber, H. Keller, Driven cavity flows by efficient numerical techniques, J. Comput. Phys. 49 (1983) 310 – 333.
[52] F. Penot, Numerical calculation of 2d natural convection in isothermal
AN US
1136
CR IP T
1132
open cavities, Numer. Heat Transfer 5 (1982) 421–437. [53] O. LeQuere, J. A. C. Humphery, F. S. Sherman, Numerical calculation
1141
of thermally driven 2d unsteady laminar flow in cavities of rectangular
1142
cross section, Numer. Heat Transfer 4 (1981) 249–283.
M
1140
[54] B. Hajra, A review of some recent studies on buoyoncy driven flows in
1144
an urban environment, International Journal of Atmospheric Sciences
1145
2014 (2014) 15, Article ID 362182.
PT
1147
[55] G. I. Taylor, A. E. Green, Mechanism of the production of small eddies from large ones, Proc. R. Soc. Lond. A 158 (1937) 499–521.
CE
1146
ED
1143
[56] M. S. Shadloo, A. Zainali, M. Yildiz, A. Suleman, A robust weakly
1149
compressible sph method and its comparison with an incompressible
AC
1148
1150
sph, International Journal for Numerical Methods in Engineering 89
1151
(2012) 939–956.
1152
[57] M. Huber, F. Keller, W. Sckel, M. Hirschler, P. Kunz, S. Hassanizadeh,
1153
U. Nieken, On the physically based modeling of surface tension and 75
ACCEPTED MANUSCRIPT
1154
moving contact lines with dynamic contact angles on the continuum
1155
scale, Journal of Computational Physics 310 (2016) 459 – 477. [58] E.-S. Lee, C. Moulinec, R. Xu, D. Violeau, D. Laurence, P. Stansby,
1157
Comparisons of weakly compressible and truly incompressible algo-
1158
rithms for the sph mesh free particle method, Journal of Computational
1159
Physics 227 (2008) 8417–8436.
1161
1162
1163
[59] A. Rafiee, Modelling of generalized Newtonian lid-driven Cavity Flow
AN US
1160
CR IP T
1156
using an SPH method, The ANZIAM Journal 49 (2008) 411–422. [60] OpenFOAM Documentation, The OpenFOAM Foundation, 2.4.0 edition, 2015.
[61] M. Basa, N. J. Quinlan, M. Lastiwka, Robustness and accuracy of
1165
sph formulations for viscous flow, International Journal for Numerical
1166
Methods in Fluids 60 (2009) 1127–1148.
ED
M
1164
[62] S. Khorasanizade, J. M. M. Sousa, A detailed study of lid-driven cavity
1168
flow at moderate reynolds numbers using incompressible sph, Interna-
1169
tional Journal for Numerical Methods in Fluids 76 (2014) 653–668.
1170
[63] P. W. Cleary, J. J. Monaghan, Conduction modelling using smoothed
1171
particle hydrodynamics, Journal of Computational Physics 148 (1999)
CE 227–264.
AC
1172
PT
1167
1173
1174
[64] Y. Zhu, P. J. Fox, Smoothed particle hydrodynamics model for diffusion through porous media, Transport in Porous Media 43 (2001) 441–471.
76
ACCEPTED MANUSCRIPT
[65] S. Adami, X. Y. Hu, N. A. Adams, A new surface-tension formulation for
1176
multi-phase sph using a reproducing devergence approximation, Journal
1177
of Computational Physics 229 (2010) 2011–5021.
1178
1179
CR IP T
1175
[66] J. Crank, The mathematics of duiffusion, Oxford University Press, 2nd edition, 1975.
[67] K. Szewc, J. Pozorski, A. Tanire, Modeling of natural convection with
1181
smoothed particle hydrodynamics: Non-boussinesq formulation, Inter-
1182
national Journal of Heat and Mass Transfer 54 (2011) 4807 – 4816.
AN US
1180
[68] A. G. V., B. Firoozabadi, M. Mahdinia, 2d numerical simulation of
1184
density currents using the {SPH} projection method, European Journal
1185
of Mechanics - B/Fluids 38 (2013) 38 – 46.
M
1183
[69] A. Leroy, D. Violeau, M. Ferrand, A. Joly, Buoyancy modelling with in-
1187
compressible sph for laminar and turbulent flows, International Journal
1188
for Numerical Methods in Fluids 78 (2015) 455–474.
ED
1186
[70] J. Boussinesq, Theorie de l’ecoulement tourbillonnant et tumultueux des
1190
liquides dans les lits rectilignes a grande section, volume 1, Gauthier-
1191
Villars et Fils, 1897.
CE
PT
1189
[71] T. Pesso, S. Piva, Laminar natural convection in a square cavity: low
1193
prandtl numbers and large density differences, Int. J. Heat Mass Transfer
AC
1192
1194
1195
1196
52 (2009) 1036–1043.
[72] S. Nas, G. Tryggvason, Thermocapillary interaction of two bubbles or drops, Int. J. Multi. Flow 29 (2003) 1117–1135. 77
ACCEPTED MANUSCRIPT
1198
1199
1200
[73] J. G. Kirkwood, F. P. Buff, The statistical mechanical theory of surface tension, The Journal of Chemical Physics 17 (1949) 338–343. [74] J. J. Monaghan, Smoothed particle hydrodynamics and its diverse ap-
CR IP T
1197
plications, Annu. Rev. Fluid Mech. 44 (2012) 323–346.
[75] M. S. Shadloo, M. Yildiz, Numerical modeling of kelvinhelmholtz in-
1202
stability using smoothed particle hydrodynamics, International Journal
1203
for Numerical Methods in Engineering 87 (2011) 988–1006.
1205
1206
1207
[76] N. Young, J. Goldstein, M. Block, The motion of bubbles in a vertical temperature gradient, J. Fluid Mech. 6 (1959) 350–356. [77] K. Szewc, J. Pozorski, J.-P. Minier, Spurious interface fragmentation in multiphase sph, Int. J. Numer. Meth. Engng 103(9) (2015) 625–649.
M
1204
AN US
1201
[78] S. Heß, V. Springel, Particle hydrodynamics with tessellation tech-
1209
niques, Monthly Notices of the Royal Astronomical Society 406 (2010)
1210
2289–2311.
1213
Royal Society of London 29 (1879) 71–97. [80] R. Balasubramaniam, A.-T. Chai,
Thermocapillary migration of
droplets: An exact solution for small marangoni numbers, Journal of
AC
1214
PT
1212
[79] L. Rayleigh, On the capillary phenomena of jets, Proceedings of the
CE
1211
ED
1208
1215
Colloid and Interface Science 119 (1987) 531 – 538.
1216
[81] G. Wozniak, R. Balasubramaniam, H. P. Hadland, S. R. Subramanian,
1217
Temperature fields in a liquid due to the thermocapillary motion of
1218
bubbles and drops, Experiments in Fluids 31 (2001) 84–89. 78
ACCEPTED MANUSCRIPT
[82] Q. Kang, L. Hu, C. Huang, H. Cui, L. Duan, W. Hu, Experimental
1220
investigations on interaction of two drops by thermocapillary-buoyancy
1221
migration, International Journal of Heat and Mass Transfer 49 (2006)
1222
2636 – 2641.
1224
1225
Drops in Reduced Gravity, Cambridge University Press, 2001. [84] A. Fedosov,
Thermocapillary motion,
arXiv:1303.0243.
AC
CE
PT
ED
M
1226
[83] R. S. Subramanian, R. Balasubramaniam, The Motion of Bubbles and
arXiv preprint (2013)
AN US
1223
CR IP T
1219
79