Modelling free surface flows with smoothed particle hydrodynamics
In this paper the method of Smoothed Particle Hydrodynamics (SPH) is extended to include an adaptive density kernel estimation (ADKE) procedure. It is shown that for a van der Waals (vdW) fluid, this method can be used to deal with free-surface phenomena without difficulties. In particular, arbitr...
Збережено в:
| Дата: | 2006 |
|---|---|
| Автори: | , , |
| Формат: | Стаття |
| Мова: | English |
| Опубліковано: |
Інститут фізики конденсованих систем НАН України
2006
|
| Назва видання: | Condensed Matter Physics |
| Онлайн доступ: | http://dspace.nbuv.gov.ua/handle/123456789/121317 |
| Теги: |
Додати тег
Немає тегів, Будьте першим, хто поставить тег для цього запису!
|
| Назва журналу: | Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| Цитувати: | Modelling free surface flows with smoothed particle hydrodynamics / L.Di G. Sigalotti, J. Daza, A. Donoso // Condensed Matter Physics. — 2006. — Т. 9, № 2(46). — С. 359–366. — Бібліогр.: 28 назв. — англ. |
Репозитарії
Digital Library of Periodicals of National Academy of Sciences of Ukraine| id |
irk-123456789-121317 |
|---|---|
| record_format |
dspace |
| spelling |
irk-123456789-1213172017-06-15T03:03:12Z Modelling free surface flows with smoothed particle hydrodynamics Sigalotti, L.Di G. Daza, J. Donoso, A. In this paper the method of Smoothed Particle Hydrodynamics (SPH) is extended to include an adaptive density kernel estimation (ADKE) procedure. It is shown that for a van der Waals (vdW) fluid, this method can be used to deal with free-surface phenomena without difficulties. In particular, arbitrary moving boundaries can be easily handled because surface tension is effectively simulated by the cohesive pressure forces. Moreover, the ADKE method is seen to increase both the accuracy and stability of SPH since it allows the width of the kernel interpolant to vary locally in a way that only the minimum necessary smoothing is applied at and near free surfaces and sharp fluid-fluid interfaces. The method is robust and easy to implement. Examples of its resolving power are given for both the formation of a circular liquid drop under surface tension and the nonlinear oscillation of excited drops. 2006 Article Modelling free surface flows with smoothed particle hydrodynamics / L.Di G. Sigalotti, J. Daza, A. Donoso // Condensed Matter Physics. — 2006. — Т. 9, № 2(46). — С. 359–366. — Бібліогр.: 28 назв. — англ. 1607-324X PACS: 68.03.-g, 68.03.Cd, 44.10.+i, 47.11.+j DOI:10.5488/CMP.9.2.359 http://dspace.nbuv.gov.ua/handle/123456789/121317 en Condensed Matter Physics Інститут фізики конденсованих систем НАН України |
| institution |
Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| collection |
DSpace DC |
| language |
English |
| description |
In this paper the method of Smoothed Particle Hydrodynamics (SPH) is extended to include an adaptive
density kernel estimation (ADKE) procedure. It is shown that for a van der Waals (vdW) fluid, this method can
be used to deal with free-surface phenomena without difficulties. In particular, arbitrary moving boundaries can
be easily handled because surface tension is effectively simulated by the cohesive pressure forces. Moreover,
the ADKE method is seen to increase both the accuracy and stability of SPH since it allows the width of
the kernel interpolant to vary locally in a way that only the minimum necessary smoothing is applied at and
near free surfaces and sharp fluid-fluid interfaces. The method is robust and easy to implement. Examples
of its resolving power are given for both the formation of a circular liquid drop under surface tension and the
nonlinear oscillation of excited drops. |
| format |
Article |
| author |
Sigalotti, L.Di G. Daza, J. Donoso, A. |
| spellingShingle |
Sigalotti, L.Di G. Daza, J. Donoso, A. Modelling free surface flows with smoothed particle hydrodynamics Condensed Matter Physics |
| author_facet |
Sigalotti, L.Di G. Daza, J. Donoso, A. |
| author_sort |
Sigalotti, L.Di G. |
| title |
Modelling free surface flows with smoothed particle hydrodynamics |
| title_short |
Modelling free surface flows with smoothed particle hydrodynamics |
| title_full |
Modelling free surface flows with smoothed particle hydrodynamics |
| title_fullStr |
Modelling free surface flows with smoothed particle hydrodynamics |
| title_full_unstemmed |
Modelling free surface flows with smoothed particle hydrodynamics |
| title_sort |
modelling free surface flows with smoothed particle hydrodynamics |
| publisher |
Інститут фізики конденсованих систем НАН України |
| publishDate |
2006 |
| url |
http://dspace.nbuv.gov.ua/handle/123456789/121317 |
| citation_txt |
Modelling free surface flows with smoothed particle hydrodynamics / L.Di G. Sigalotti, J. Daza, A. Donoso // Condensed Matter Physics. — 2006. — Т. 9, № 2(46). — С. 359–366. — Бібліогр.: 28 назв. — англ. |
| series |
Condensed Matter Physics |
| work_keys_str_mv |
AT sigalottildig modellingfreesurfaceflowswithsmoothedparticlehydrodynamics AT dazaj modellingfreesurfaceflowswithsmoothedparticlehydrodynamics AT donosoa modellingfreesurfaceflowswithsmoothedparticlehydrodynamics |
| first_indexed |
2025-07-08T19:38:38Z |
| last_indexed |
2025-07-08T19:38:38Z |
| _version_ |
1837108844826722304 |
| fulltext |
Condensed Matter Physics 2006, Vol. 9, No 2(46), pp. 359–366
Modelling free surface flows with smoothed particle
hydrodynamics
L.Di G.Sigalotti∗, J.Daza, A.Donoso
Centro de Fı́sica, Instituto Venezolano de Investigaciones Cientı́ficas, IVIC,
Apartado Postal 21827, Caracas 1020A, Venezuela
Received March 9, 2006, in final form April 17, 2006
In this paper the method of Smoothed Particle Hydrodynamics (SPH) is extended to include an adaptive
density kernel estimation (ADKE) procedure. It is shown that for a van der Waals (vdW) fluid, this method can
be used to deal with free-surface phenomena without difficulties. In particular, arbitrary moving boundaries can
be easily handled because surface tension is effectively simulated by the cohesive pressure forces. Moreover,
the ADKE method is seen to increase both the accuracy and stability of SPH since it allows the width of
the kernel interpolant to vary locally in a way that only the minimum necessary smoothing is applied at and
near free surfaces and sharp fluid-fluid interfaces. The method is robust and easy to implement. Examples
of its resolving power are given for both the formation of a circular liquid drop under surface tension and the
nonlinear oscillation of excited drops.
Key words: gas-liquid and vacuum-liquid interfaces, surface tension and related phenomena, heat
conduction, computational methods in fluid dynamics
PACS: 68.03.-g, 68.03.Cd, 44.10.+i, 47.11.+j
1. Introduction
Gas-liquid flows are examples of typical two-phase fluids encountered in many systems within
both natural and artificial environments, where the liquid is separated from the gas by a free
surface whose shape is determined by the flow field. In particular, liquid surfaces are known to
be in a state of tension because fluid molecules at or near the surface experience uneven forces
of attraction [1]. Hence surface tension is the result of a microscopic, localized force that exerts
itself on the fluid elements at the surface in both the normal and tangential directions. Examples
of free-surface motion induced by tension forces can be found in a variety of flows which are
either natural or artificially constructed. For instance, well-known problems involving liquid-gas
and liquid-vacuum interfaces result from experiments and studies of droplet dynamics, which may
include the formation and evaporation of equilibrium liquid drops [2,3], the nonlinear oscillations
of excited drops [4–7], the breakup of liquid filaments [8,9] and the coalescence of colliding drops
[10–15].
In general, the detailed analysis of free-surface flows involves the use of numerical methods
coupled to appropriate boundary conditions to simulate an arbitrarily moving surface. However,
most of the available methods suffer from difficulties in modelling topologically complex interfaces
having surface tension. Among the earlier numerical schemes, the MAC method [16] was known to
be the most flexible and robust. It relies on finite differences to solve the hydrodynamic equations
and on particles to define the surface. The method was further extended to deal with moving
boundaries [17] and applied to a wide variety of problems, including waterfalls, breaking dams and
two-fluid instabilities [18]. In spite of attempts to simplify it, it remains complicated to program.
Later on, Brackbill et al. [19] developed the continuum-surface-force (CSF) method, which treats
surface tension as a continuous, three-dimensional effect across the interface rather than as a pure
boundary condition. In this model, the interfaces between fluids of different properties, or colours,
∗E-mail: lsigalot@cassini.ivic.ve
c© L.Di G.Sigalotti, J.Daza, A.Donoso 359
L.Di G.Sigalotti, J.Daza, A.Donoso
are represented as transition layers of finite thickness, across which the colour variable varies conti-
nuously. Although this method has been successfully applied to simulate fluid phenomena involving
incompressible fluid flow in low-gravity environments, capillarity and droplet dynamics, it still re-
quires the use of complicated surface-tracking techniques, like the MAC method, to accurately
model the volume force along the surface curvature. Alternative surface-tracking methods have
been proposed that rely on the spine-flux technique [20]. However, a serious drawback of these lat-
ter methods is that they are limited to surface deformations that do not result in double-definition
of the surface along a spine. On the other hand, an application of the CSF model to calculate
interfacial curvatures with Smoothed Particle Hydrodynamics (SPH) was reported by Morris [21].
Moreover, Monaghan [22] showed that SPH can be easily extended to handle free-surface phenom-
ena if the continuity equation, rather than the usual summation interpolant, is used in calculating
the density and if the boundary particles are allowed to move with a corrected velocity.
Here we present an alternative approach to model free-surface flows under tension forces with
standard SPH. The method is modified to account for an adaptive density kernel estimation
(ADKE) procedure similar to that described by Silverman [23]. In this procedure the basic ker-
nel and nearest neighbour SPH approaches are combined in a way that the amount of applied
smoothing is minimized in regions where the density is lower. As a result, near a free surface,
where there is a crippling deficiency of particles, the ADKE method uses much smaller kernels
than in the bulk of the fluid without compromising the SPH stability, leading to a finer resolution
of the free surface compared to other SPH schemes. Coupled to the fully Lagrangian nature of
SPH, the ADKE method maintains very well sharp fluid-fluid interfaces and free surfaces with
no need of employing explicit surface reconstruction. The method was also found to give stable
and accurate results for supersonic compressible flows with strong shocks [24]. In this paper, we
focus on two-dimensional test calculations of the formation of circular drops under surface tension
and its nonlinear oscillations starting from elongated shapes. The paper is organized as follows.
In section 2 we write down the SPH hydrodynamic equations that are solved and provide a brief
description of the method. Then in section 3 we describe the results of a few model calculations
and in section 4 we give the conclusions.
2. SPH equations and numerical method
In contrast to Monaghan’s [22] approach, the density ρa of the particles is computed using the
summation interpolant
ρa =
N
∑
b=1
mbWab , (1)
where mb is the mass of particle b, Wab = W (|xa−xb |, h) is a spherically symmetric kernel function
and h is the so-called smoothing length or width of the kernel. In this work, we adopt the quartic
spline kernel of Lucy [25] so that only nearest particles within a circle of radius h contribute to
the SPH summation in equation (1). The same kernel was used by Nugent and Posch [2] in SPH
simulations of drop condensation. Other smoothing kernels, such as the cubic B-spline [26] or the
quintic spline [21], have been used by other authors. Test simulations of drop condensation with
the quartic and quintic spline kernels have shown that the results are almost independent of the
form of the kernel. For the SPH equations of motion and thermal energy we use a symmetrized
representation that contains the full Newtonian viscous stress tensor Sij and the Fourier’s heat
flux vector qi, namely
dxi
a
dt
= vi
a , (2)
dvi
a
dt
=
N
∑
b=1
mb
(
Sij
a
ρ2
a
+
Sij
b
ρ2
b
)
∂Wab
∂xj
a
, (3)
dUa
dt
=
1
2
N
∑
b=1
mb
(
Sij
a
ρ2
a
+
Sij
b
ρ2
b
)
(
vi
b − vi
a
) ∂Wab
∂xj
a
−
N
∑
b=1
mb
(
qj
a
ρ2
a
+
qj
b
ρ2
b
)
∂Wab
∂xj
a
, (4)
360
Modelling free surface flows
where xi
a denotes the position of the particles, vi
a is their velocity and Ua is their thermal energy
per unit mass. The viscous stress tensor is given by
Sij = −pδij + η
(
∂vi
∂xj
+
∂vj
∂xi
)
+
(
ζ −
2
d
η
)
∂vk
∂xk
δij , (5)
and the heat flux vector is defined according to
qi = −κ
∂T
∂xi
. (6)
In these relations, p is the internal (isotropic) pressure, T is the fluid temperature, δij is the
unit tensor, η and ζ are the coefficients of shear and bulk viscosity, κ is the coefficient of heat
conduction and d is a number equal to 2 or 3 depending on whether the motion is confined to two-
or three-dimensions, respectively.
For the drop model of this paper, we assume a van der Waals (vdW) fluid, where the pressure
and thermal energy obey the constitutive relations
p =
ρk̄BT
1 − ρb̄
− āρ2, (7)
U =
ξ
2
k̄BT − āρ. (8)
In these equations, ξ is the number of degrees of freedom of the molecules and k̄B = kB/m, where
kB is the Boltzmann’s constant and m is the particle mass. Furthermore, ā = a/m2 and b̄ = b/m,
where a controls the strength of the attractive forces between the molecules and b is a constant
due to their finite size.
As was first noted by Nugent and Posch [2], surface tension effects can be simulated with the
aid of equation (7) by separating the cohesive term, −āρ2, from the remainder forces. The same
applies to the energy term, −āρ, in equation (8). The former term contributes with an attractive
central force between the SPH particles, which produces a net acceleration given by
dvi
a
dt
= 2ā
N
∑
b=1
mb
∂WH
ab
∂xi
a
, (9)
while the latter term contributes with an effective heating due to the work done by the cohesive
pressure forces on the liquid within the free surface and takes the SPH form
dUa
dt
= 2ā
N
∑
b=1
mb
(
vi
b − vi
a
) ∂WH
ab
∂xi
a
. (10)
In these equations WH
ab = W (|xa −xb |,H ) has the same functional form of the interpolating kernel
employed for all other terms, except that it now depends on the smoothing range H instead of h.
In the CSF model, surface tension is translated into a force per unit volume, Fs = fsδs , where δs
is a normalized function, which peaks at the interface, and fs is a force per unit area defined as
fs = στn + ∇sσ, (11)
where σ is the surface tension coefficient, n is the unit normal to the interface, τ is the curvature
of the interface and ∇s is the surface gradient. Note that the first term in equation (11) is a
force acting normal to the interface and corresponds to the net surface tension force due to the
local curvature. This force acts to smooth regions of high curvature, in an attempt to reduce the
total surface area and hence the surface energy. In contrast, the second term in equation (11)
acts tangentially to the interface, forcing fluid elements on the surface to move from regions of low
surface tension to regions of higher surface tension. The surface normal in equation (11) is evaluated
361
L.Di G.Sigalotti, J.Daza, A.Donoso
according to n = ∇c/[c], where c is the colour function identifying each fluid in the problem and
[c] is the jump in c across the interface. Moreover, the curvature is calculated using the relation
τ = −∇ · n, whereas δs = |n|. In SPH form, Morris [21] suggested the following representations
na =
N
∑
b=1
mb
ρb
cb∇aWab , (12)
and
(∇ · n)a =
N
∑
b=1
mb
ρb
nb · ∇aWab , (13)
for the normal and the surface curvature, respectively. In particular, the divergence of the surface
normal is proportional to the surface tension force. We note the intrinsic similarity between equa-
tions (9) and (13), implying that the former acts effectively as a true surface tension force due to
the local curvature. It then follows that the cohesive term in equation (7) is all we need to correctly
model surface tension in a vdW fluid. Furthermore, the choice of H in equations (9) and (10) is
determined by the same considerations that led to a substantial improvement of the interface sta-
bility properties with the CSF method employed by Morris [21]. In practice, stability is maintained
when H > 2h. In a fluid bounded by a free surface, the forces represented by equation (9) largely
cancel within the bulk of the fluid, except for a small strip H around the surface where particles
are accelerated in the direction of the inward surface normal. In actual numerical simulations, the
value of H determines the thickness of the interface.
Apart from showing how the vdW cohesive pressure can be used with SPH to model free-surface
phenomena we show why the ADKE procedure is an essential ingredient to maintain both adequate
stability and fine resolution of the interface. The essence of the ADKE method is to construct a
collection of local kernels at the position of the particles in order to allow the smoothing length to
vary from point to point. In order to do so, an initial, or pilot, density estimate, ρ̂a, is calculated
using the summation interpolant in equation (1) with Wab = W (|xa −xb |, h0), where h0 is defined
as some dilation factor of the initial interparticle separation, i.e., h0 = D∆s. Local bandwidth
factors, λa, are next calculated according to
λa = k
(
ρ̂a
ḡ
)
−ε
, (14)
where ḡ is the geometric mean of the pilot estimates given by
log ḡ =
1
N
N
∑
b=1
log ρ̂b , (15)
k is a scaling factor of the order of unity and ε is the sensitivity parameter defined in the interval 0 6
ε 6 1. The adaptive estimator is finally constructed by recalculating the density using equation (1)
with h0 replaced by ha = λah0. As long as ε → 1, the ha’s become more sensitive to the density
distribution, implying a greater difference between the smoothing length in different parts of the
fluid domain. This means that near a sharp interface, where there is a crippling deficiency of
particles, the size of the kernels will be much smaller than in the bulk of the fluid, leading to
a finer resolution of the interface in spite of the fact that H > 2h in order to ensure numerical
stability. Moreover, setting ε = 0, λa → k and the ADKE method reduces to the usual fixed width
kernel approach. Also, note that with the λa’s calculated as in equation (14), the ha’s effectively
control the scale of smoothing applied to the data, thus improving the accuracy and stability near
the interface.
3. Numerical tests
In this section, we consider a series of two-dimensional simulations for the specific case of a
forming vdW liquid drop under surface-tension forces, in which the effects of viscosity and heat
362
Modelling free surface flows
conduction are taken into account. This model represents a sensitive test calculation to quantify
the resolving power of the ADKE method on free-surface flows. Its ability to handle arbitrarily
moving surfaces is also tested against the nonlinear oscillations of an initially deformed drop.
3.1. Formation of equilibrium vdW drops
We first consider the formation of a stable, circular drop using the same vdW parameters
employed by Nugent and Posch [2] and Meleán et al. [3], namely m = 1, ā = 2, b̄ = 0.5 and k̄B = 1.
In these reduced units, the critical point of the vdW fluid occurs for ρcr = 2/3, pcr = 8/27 and
Tcr = 32/27 [27]. For simplicity, all variables are normalized with respect to their critical values.
The coefficients of thermal conductivity, shear and bulk viscosity in reduced units are chosen to
be κ = 5, η = 1 and ζ = 0.1, respectively. In addition, the outer medium is a vacuum so that
there is no heat conduction between the drop and its surroundings. All calculations start with
900 SPH particles of equal mass (ma = m = 1) distributed at the vertices of a square cell array.
Initially the particles are at rest and separated from one another by a distance ∆s = 0.75 along the
x- and y-axes. The choice of these parameters is such that the initial density, as calculated using
equation (1), satisfies the constraint ρ0 < 1/b̄. Each particle is given an initial temperature T0 = 0.2
so that the inequality k̄BT0 > 2āρ0(1 − ρ0b̄)
2 is also satisfied for thermodynamic stability. At this
subcritical temperature, an equilibrium, circular drop with no external atmosphere is expected to
form [2]. Therefore, no periodic boundary conditions are required far away from the initial square
array of particles. All runs were made by setting H = 2h in equations (9) and (10).
�
���������
�
��
��
�
�
��
�
����� �
���
����� �
���
�
�
��
��
�
�
���
���
��
�
����������
���������
Figure 1. Particle positions in the (x,y)-plane
showing the structure of a circular liquid drop, with
D = 3, k = 1 and ε = 0 (top left panel); D = 3,
k = 0.9 and ε = 0.6 (top right panel); D = 4,
k = 0.6 and ε = 0.8 (bottom left panel) and D = 5,
k = 0.5 and ε = 0.8 (bottom right panel).
We consider three distinct sequences of
computations, all starting with identical ini-
tial conditions. Each sequence is defined
by fixing the value of D and varying the
ADKE parameters k and ε. In particular,
we choose D = 3, 4 and 5 so that initially
h0 = D∆s, and values of k and ε in the
intervals [0.5,1.0] and [0.5,0.9], respectively.
A typical evolution proceeds as follows. At
the beginning, the internal pressure is not
enough to balance the surface tension forces,
causing an inward motion of the particles at
and near the corners of the square bound-
ary. As a result, part of the surface ten-
sion energy is converted into internal liquid
movement with a consequent sharp increase
of the kinetic energy. Circularization of the
confined liquid is approximately completed
when about 95 percent of the kinetic energy
is dissipated by viscous forces. During this
stage, the thermal energy also increases due
to the heating produced by internal friction
near the drop surface. The increase of ther-
mal energy soon slows down due to heat con-
duction. The further evolution consists of a
very slow decrease of the kinetic energy accompanied by an equally slow increase of the thermal
energy until the drop reaches thermo-mechanical equilibrium with a circular shape, as shown in
figure 1. With a constant timestep ∆t = 0.005, a typical evolution is completed after about 130,000
computational cycles.
The four panels of figure 1 show the final equilibrium drop as obtained for different choices of
the ADKE parameters. The top left panel corresponds to the case when the ADKE procedure is
dropped (i.e., k = 1, ε = 0). Since for this model the dilation factor is D = 3, the calculation was
made using a fixed kernel approach of width h = 3∆s, resulting in H = 6∆s = 4.5. We may see
363
L.Di G.Sigalotti, J.Daza, A.Donoso
-20
-10
0
10
20
y
t = 0
t = 18
t = 28
t = 35
t = 62
20
10
0
-10
-20
20100-10-20
y
x
t = 75
20100-10-20
x
t = 82
20100-10-20
x
t = 88
20100-10-20
x
t = 100
20100-10-20
x
t = 114
Figure 2. Sequence of times showing the first period of oscillation of a liquid drop when it is
released from an elliptic shape with aspect ratio a/b = 3 (t = 0) at Re ≈ 522. In each panel the
time is shown in reduced units.
that in this case the transition layer separating the drop from the outer vacuum is considerably
thicker compared to the other three models where the ADKE procedure was activated. In fact, the
drop surface is numerically resolved by the outermost four concentric circles of particles of thickness
≈ 4.5. In all other cases where the ADKE was used the outer interface is clearly represented by
only one outer circle of particles. The much better resolution afforded in these latter cases is a
consequence of both the adaptive nature of the method and the fact that a minimum smoothing
is effectively applied to the data at the drop surface.
3.2. Oscillating vdW drops
Next, we test the response of the method on large-amplitude oscillations of vdW drops that are
released from a static elliptic shape. A sequence of equilibrium circular drops was first constructed
with the same parameters as before, except that the shear viscosity was varied in the range 0.012 6
η 6 12, corresponding to Reynolds number in the interval 0.5 < Re < 522. In all cases, the initial
smoothing length was chosen to be h0 = 3∆s, while the ADKE parameters were fixed to k = 0.9
and ε = 0.6. The reference circular drops were deformed into an elliptic shape by homogeneous
flattening under pure shear strain [28], which is achieved with an area-preserving and, hence,
density-conserving transformation of the particle coordinates
x′ =
1
(1 + ε)
x, y′ = (1 + ε)y, (16)
where ε is the elongation. The above transformation converts the circle into an ellipse of semi-major
axis a aligned with the y-axis, as shown in the first panel of figure 2. Here we shall consider only
drop deformations with an initial aspect ratio a/b = 3, where b denotes the semi-minor axis of the
ellipses.
The nonlinear oscillations of a liquid drop represent a more stringent test because they involve
a moving free surface. The time resolved evolution for the case of a drop with Re ≈ 522 (η = 0.012)
is displayed in figure 2 during its first period of oscillation. The elongated drop first contracts along
its major axis as most of its surface energy is being transformed into internal liquid movement.
It then passes through a transient circular shape (t = 28) before reaching a stage of maximum
elongation along the x-axis (t = 62). At this point, most of the undamped internal kinetic energy
is in the form of surface energy. The rim pressure soon exceeds the stagnation pressure inside the
drop, causing it to contract back (along the x-axis) under surface tension. The drop evolves again
to an approximate circular shape before achieving a prolate shape towards the completion of the
first oscillation period at t = 114. At this time, the aspect ratio is ≈ 1.94 and decays further
to ≈ 1.35 by the end of the second period (t = 221). The oscillatory behaviour will continue in
364
Modelling free surface flows
time with progressively lower amplitudes until eventually the drop loses all of its internal kinetic
energy by viscous dissipation. The evolution was followed for about 20 periods by the time the
drop has returned to its equilibrium circular shape. The damping of the oscillations is mostly
due to viscous dissipation and, to some extent, to the finite heat conductivity (κ = 5). A value
of κ this large serves to reduce temperature gradients and density fluctuations within the drop.
However, the intrinsic viscosity which is inherent to particle systems due to the smoothing may
also contribute to dissipation of the drop oscillations even when η = 0. A test simulation with η = 0
and κ = 0 shows that by switching off the ADKE procedure, the rate of dissipation due to the
inherent numerical viscosity is much slower than that due to either finite physical viscosity or heat
conductivity. For comparison, the ADKE method (with k = 0.9 and ε = 0.6) resulted in even slower
dissipation rates, mainly because much less smoothing was actually needed at and near the drop
surface. Evidently, these results show that the method is able to handle moving surfaces without
difficulties. Another important feature is that during surface motion the boundary particles tended
to move from regions of low curvature to regions of higher curvature, responding to the action of
the tangential component of the surface tension force.
4. Conclusions
We have performed test calculations which show that standard SPH can be used to accurately
simulate free-surface flows without difficulty for a van der Waals (vdW) fluid, provided that the
kernel function is modified to include an adaptive density kernel estimation (ADKE) procedure, like
the one proposed by Silverman [23]. When dealing with a vdW-type equation of state, the SPH
representation of the cohesive pressure gradients strongly resembles that employed to calculate
interfacial curvatures with SPH using the continuum-surface-force model [21], implying that surface
tension is correctly simulated by the cohesive part of the vdW pressure. The ADKE procedure
increases the resolving power of SPH and provides a satisfactory representation of free boundaries.
In particular, the ADKE method constructs local bandwidth factors, or kernel estimates, evaluated
at the position of particles such that the width of the interpolating kernel is allowed to vary from
point to point. This results in much smaller sizes of the kernel close to sharp boundaries or in
regions where the density is lower, implying that the amount of smoothing required in those
regions is effectively minimized. This is the key point which allows for a much better description of
free surfaces and sharp fluid-fluid interfaces than the fixed width kernel approach does. Here the
method has been successfully tested for a forming two-dimensional, circular drop under the effects
of surface tension and for the two-dimensional oscillations of an initially elongated drop about a
circular shape. The method is simple to implement and does not require special modifications when
applied to fully three-dimensional problems.
Acknowledgements
This work is financially supported by the Instituto Venezolano de Investigaciones Cient́ıficas
(IVIC). One of us (L. Di G. S.) acknowledges the Fondo Nacional de Ciencia, Tecnoloǵıa e Inno-
vación (Fonacit) of Venezuela for partial support.
References
1. Nelleon M. Mechanics and Properties of Matter. Heineman, London, 1952.
2. Nugent S., Posch H.A., Phys. Rev. E, 2000, 62(4), 4968.
3. Meleán Y., Sigalotti L.Di G., Hasmy A., Comput. Phys. Commun., 2004, 157, 191.
4. Basaran O.A., J. Fluid Mech., 1992, 194, 169.
5. Becker E., Hiller W.J., Kowalewski J., Fluid Mech., 1994, 258, 191.
6. Apfel R.E. et al., Phys. Rev. Lett., 1997, 78(10), 1912.
7. Mashayek F., Ashgriz N., Phys. Fluids, 1998, 10(5), 1071.
8. Schulkes R.M.S.M., J. Fluid Mech., 1996, 309, 277.
9. Notz P.K., Basaran O.A., J. Fluid Mech., 2004, 512, 223.
365
L.Di G.Sigalotti, J.Daza, A.Donoso
10. Jiang Y.J., Umemura A., Law C.K., J. Fluid Mech., 1992, 234, 171.
11. Qian J., Law C.K., J. Fluid Mech., 1997, 331, 59.
12. Menchaca-Rocha A., Huidobro F., Martinez-Davalos A., Michaelian K., Perez A., Rodriguez V.,
Cârjan N., J. Fluid Mech., 1997, 346, 291.
13. Willis K.D., Orme M.E., Exp. Fluids, 2000, 29, 347.
14. Mashayek F., Ashgriz N., Minkowycz W.J., Shotorban B., Int. J. Heat Mass Transfer, 2003, 46, 77.
15. Meleán Y., Sigalotti L.Di G., Int. J. Heat Mass Transfer, 2005, 48, 4041.
16. Harlow F., Welch J.E., Phys. Fluids, 1965, 9, 842.
17. Viecelli J., J. Comput. Phys., 1971, 8, 119.
18. Harlow F., Comput. Phys. Commun., 1988, 48, 1.
19. Brackbill J.U., Kothe D.B., Zemach C., J. Comput. Phys., 1992, 100, 335.
20. Mashayek F., Ashgriz N., J. Comput. Phys., 1995, 122, 367.
21. Morris J.P., Int. J. Numer. Meth. Fluids, 2000, 33, 333.
22. Monaghan J.J., J. Comput. Phys., 1994, 110, 399.
23. Silverman B.W. Density Estimation for Statistics and Data Analysis. Chapman & Hall, London, 1996.
24. Sigalotti L.Di G., López H., Donoso A., Sira E., Klapp J., J. Comput. Phys., 2006, 212, 124.
25. Lucy L.B., Astron J., 1977, 80, 1013.
26. Monaghan J.J., Lattanzio J.C., Astron. Astrophys., 1985, 149, 135.
27. Landau L.D., Lifshitz E.M. Statistical Physics. Addison-Wesley, Massachusetts, 1969.
28. Twiss R.J., Moores E.M. Structural Geology. Freeman, New York, 1992.
����������� �i����� ����������� �����i� i !���"����#�$�������� !i�������%i���
&
.'i (.)i*+,-..i, /.'+0+, 1.'-2-0-
34567 8i9:;:, <454=>4?@=@;:A i5=6:6>6 5B>;CD:E FC=?iFG45@,HCI6. =;7:5@;B 21827, JB7B;B= 1020–K, <454=>4?B
LMNOPQRS 9 TUNUVURW 2006 N., X SYMQMSZRSP[ XO\]^_i – 17 `XiMR^ 2006 N.
a46CF 9b?BFG45Cc dB=6:5;CDCc 4?4;67CF:5Bei;: 7C9I:745C F?f D7BE>DB55f H7Cg4F>7: Cgi5;: BFB-H6:D5CbC fF7B b>=6:5:. hC;B9B5C, iCF?f DB5-F47-<BB?@=CDCc 7iF:5: g4A e46CF eCG4 j49 67>F5CiiD
9B=6C=CD>DB6:=f H7: CH:=i fD:i HCj?:9> Di?@5Cc HCD47E5i. kC;74eB, ?4b;C D7BE>DB6: FCDi?@5C 7>-ECei HCD47E5i, C=;i?@;: HCD47E54D:A 5B6fb 484;6:D5C =:e>?lm6@=f =:?Be: ;Cb49:D5CbC 6:=;>. ni?@-I4 6CbC, 6Cd5i=6@ i =6iA;i=6@ g@CbC e46CF> D:iB, C=;i?@;: Di5 FC9DC?fm ?C;B?@5C 9ei5lDB6: I:7:5>
i5647HC?f56: fF7B 6B;, iC 9B=6C=CD>m6@=f ei5ieB?@54 54CjEiF54 9b?BFG>DB55f HCj?:9> Di?@5:E HC-D47EC5@ i HCD47EC5@ 7C9Fi?> 8?lcF-8?lcF.
a46CF m =6iA;:e i ?4b;:e > 9B=6C=>DB55i. h7:DCFf6@=fH7:;?BF: ACbC D:;C7:=6B55f F?f CH:=> 8C7e>DB55f ;7BH?:5 7iF:5: HiF DH?:DCe HCD47E54DCbC 5B-6fb> 6B 54?i5iA5CbC ;C?:DB55f 9j>FG45:E ;7BH4?@.
opqrsti upstv: wxyz{|}i {x~� i�� ��~-{i��}� i y�����-{i��}�, wxyz{|}zy�� }���� i wxy’�~�}i �y���,�zw�xw{xyi�}i���, x������y���}i �z�x�� y ��}��i�i ����� iy
PACS: 68.03.-g, 68.03.Cd, 44.10.+i, 47.11.+j
366
|