16
1.3.1. Introduction
Droplet evaporation was extensively studied during the last three decades. Deegan et al.
7
studied
the so-called “coffee stain” effect, the well-known phenomena related to the drying droplet of
coffee on a hydrophilic surface, and derived analytically the averaged velocity profile along the z-
axis (normal to the surface). Popov
8,
9
modelled the evaporation of a droplet of colloidal solutions
upon wettable or not wettable surfaces. In particular, Popov
9
deduced a general model taking into
account different values of CAs, humidity levels, interface droplet area and others important
parameters. Hui et Al
10
used the finite element method to figure out the velocity field in a drying
droplet. Afterwards they integrated the Marangoni effect in their model
11
that takes into account the
shear stresses along the droplet-air interface due to the temperature gradient. A programming code
in Matlab was written by me, inspired by these previous results, that uses the finite element method
for the calculations of the flows inside the drying droplets, also on SHSs. The program allows
solving a model in which there are some approximations, i.e. the gravity. The partial differential
equations used for the system are eq. (14) for the evaporation and eq. (18) for the convective flows
inside the droplet. I will explain how to couple the two equations for solving the change of phase
between solid and air. I will also study the dynamics of evaporation of a drying droplet upon a
hydrophilic and a SHS by considering the Popov’s model.
9
1.3.2. Convective flow inside a drying droplet.
A numerical simulation by finite element method was made to obtain the evaporation flux profile
along the free surface of a drying droplet of pure water upon a hydrophilic surface that affects the
inner flow. The model was programmed in Matlab. As written in the previous section, the system
can be considered with an axisymmetric geometry (figure 1.5). Aqueous vapor concentration
changes and becomes steady rapidly respect to the time of evaporation and whenever the droplet
shape changes. For this reason the evaporation was thought to be a quasi-steady-state process
12
and
I solved the diffusion problem given by the gradient of vapor concentration between the free surface
and the atmosphere. The eq. (14) becomes:
7
R. Deegan, O. Bakajin, T. F. Dupont, G. Huber, S. R. Nagel, T. A. Witten. “Capillary flow as the cause of ring stains
from dried liquid drops”. J. F. Institute. Nature(1997), 389, 827-829;
8
Y. O. Popov. “Singularities, universality, and scaling in evaporative deposition patterns”. PhD student’s thesis.
Department of Physics. University of Chicago (US), (2003);
9
Y. O. Popov. “Evaporative deposition patterns: spatial dimensions of the deposits”. Phys. Rev. E (2005), 71, 036313;
10
H. Hu and R. G. Larson. “Analysis of the Microfluid Flow in an Evaporating Sessile Droplet”. Langmuir (2005), 21,
3963-3971;
11
H. Hu and R. G. Larson. “Analysis of the Effects of Marangoni Stresses on the Microflow in an evaporating Sessile
Droplet”. Langmuir (2005), 21, 3972-3980;
0 C
(19)
17
The eq. (19) is also called Laplace equation. The limit of the domain over the free surface of the
droplet is assumed to be 10 times the interface radius (R), (figure 1.6A). The vapor concentration
(Cw) corresponds to that of the surrounding environment. The concentration of vapor at the free
surface is
9
3
2.32 10
S
g
C
mm
,
12
the humidity level (H) at room temperature was set to 40%,
giving a vapor concentration in the air:
9
3
0.92 10
V
g
C H C
mm
. The diffusion coefficient of
water vapor was set to
2
26.1
mm
s
.
12
The domain was divided in triangular elements. The calculus
of the nodes and elements is performed in Ansys. I exported nodes and elements and I solved the
system in Matlab by finite element method. Two meshing were built for the hydrophilic case and
the superhydrophobic one (figure 1.6 and 1.7). For the implementation in Matlab and the solution
of Laplace equation (that gives the evaporation rate) by finite element method, I refer to K.
Petersen.
13
We obtain the weak form of the equation by multiplying both sides by a test function
(i.e. a function which is infinitely differentiable and has compact support), integrating over the
domain , and performing integration by parts.
13
The eq. (19) could be expressed as:
c vdx fvdx
(20)
for all test functions v. The vector f, in eq. (20), was set to 0. To find an approximation to the
solution u, we choose a finite-dimensional space Vh and ask that equation (20) is satisfied only for v
in Vh rather than for all test functions v. Then we look for a function
HH
uV which satisfies:
13
H
for all v V
H
u vdx fvdx
(21)
The assemble matrix is written as
13
:
S
ij i j
TT
T
K dx
(22)
where TH is the set of (mesh) elements and T is an element. In order to assemble the stiffness
matrix, we need to compute integrals of the form:
ij
T
dx
(23)
12
H. Hu and R. G. Larson. “Evaporation of a Sessile Droplet on a Substrate”. J. Phys. Chem. B (2002), 106, 1334-1344;
13
K. Petersen. “Finite Elements for Laplace’s Equation: Implementation”. CFD Reading Group. Courant Institute of
Mathematical Sciences. New York University, (2009);
18
These integrals are computed by making a change of variables to a reference element. For the case
where the mesh elements are triangles, the reference element is the right triangle whose vertices are
at (0,0), (1,0), and (0,1). Every triangle can be transformed into the reference triangle using an
affine transformation, i.e. a transformation of the form:
x A b
(24)
where, α are the coordinates on the reference triangle.
13
For example, if our triangle has its vertices
at (0,0), (h,0), and (h/2, h/2), then the affine transformation is given by:
2
0
2
h
h
A
h
(25)
0
0
B
(26)
When we perform the change of variables we get:
det( ) dx A d
(27)
I implemented all this equations in Matlab. I paid attention about the rules on the derivatives that
are summarized by K. Petersen
13
.
det( ) K A d
(28)
Let’s call
1
BA
. B in my implementation will be given by:
, 2 , 2 , 4 , 2 , 2 ,3 , 4 ,3
,3 , 2 , 4 , 2 ,3 ,3 , 4 ,3
n e i n e i n e i n e i
B
n e i n e i n e i n e i
(29)
where, n and e are two arrays, the first is related to the nodes and the second is related to the
triangular elements I am considering; i is an index that varies from 1 and the total number of
element. The first element of B means the value that is the accessed in the array n, at the position
given by
,4 ,2 ei . Just to be clear, I have to specify that n and e are imported from the Ansys
files that are structured in a precise order that I reduced as:
Node index (ni) nx ny nz Element index ni (:,1) ni (:,2) ni (:,3)
1
...
tot. number
val.
...
...
val.
...
...
val
...
...
...
...
Tot. number
val.
...
...
Val.
...
...
Val.
...
...
Tables. All the nodes in figure 1.6 and 1.7 are indexed. Each node has also the three coordinates related to a Cartesian
system. Since we deal with a 2D geometry all the n z are set to 0 (left). The elements have also a indexed structure. An
element is a triangle; their verteces are given but the indeces of three nodes (right).
19
Let’s call α the matrix containing the vertices coordinates of the reference triangle:
1 0 1
0 1 1
2D Gauss quadrature for triangular domains was used for computing the integral (23) as follows:
11
,,
e i, j 1 ,e i,k 1 K e i, j 1 ,e i,k 1 0.5 det B * B :, j B :,k ;
T
i k j
K
where i is the index related to the local stiffness matrix of each elements (so it varies from one to
the total number of elements). K is a matrix 3x3, so k and j vary from one to three each one.
The solution of the problem is given by:
1
C K f
(30)
The boundary conditions of the problem are assigned as follows:
1) the concentration to the limit of the boundary is C
(whose value was assigned before);
2) the concentration to the free surface is
S
C (whose value was assigned before);
3) no vapor flux was set between the substrate and the environment.
The equation (30) is a linear system. The resolution gives the results shown in figure 1.8A,B,C,D.
Pure water sessile droplet evaporation depends on the wetting properties of the substrate. Here a
low CA provides an evaporation flow greater at the edge and lower at the top surface of the droplet
Figure 1.6.A Case hydrophilic substrate. Triangular mesh of the external ambient. Nodes and elements were created
using the algorithm implemented in Ansys. 10823 nodes and 5270 elements were created for the external
environments. 6855 nodes and 3194 elements were created for the mesh of the droplet (B). The limit of the air domain
is 10 times the radius interface of the droplet. Both figure are edited by Matlab
A
B
20
(figure 1.8B,C,D). On the other hand, when CA is greater than 90° the flow at the edge will result
lower (figure 1.8B).
As Deegan et al. reported
7
, the evaporation affects the inner flows. The evaporation flux
configuration on a hydrophilic surface induces an outward convective flow, pinning of the droplet
and the formation of a coffee-ring residue.
7
Now, I want to derive numerically this outward
convective flow. By using the same approach than before, the linear system by finite element
method for the eq. (23) can be written as:
1
v A p
(31)
Figure 1.8.Vapor concentration upon a droplet on hydrophilic surface (CA =40°) (A) and evaporation flux above the
droplet (B). Arrow line and magnitude of evaporation flux (J) along the free surface
A
C.
D. D.
B. A B
C
D
21
where v is the velocity filed, A is the viscous term and p is the vector of forces acting on the
boundary of the droplet. All terms in eq. (31) are arrays.
As boundary condition, the no-slip holds along the bottom of the droplet at the liquid-glass
interface. On the upper free surface of the droplet, however, the boundary conditions involve a
moving boundary. The shape of the droplet can still be regarded as a spherical cap during
evaporation.
10
Boundary conditions were imposed along the droplet’s spherical free surface: a zero
shear stress was imposed along the tangential direction on the free surface excluding the Marangoni
effect and a kinematic boundary condition was used along the normal direction. As Deegan found
(figure 1.9), I considered the velocity along the r direction proportional to the evaporation flux
magnitude J already computed by solving the diffusion problem (figure 1.8 C,D).
vJ
(32)
As shown in figure figure 1.9, v(r) is perpendicular to the evaporation flux but the evaporation flux
has another direction. I introduced the two components of the evaporation flux (along r and z
direction, see figure 1.5) to determine the velocity profile as boundary condition to the free surface
and solved the stoked flow by FEM.
Figure 1.9 Reprinted from
6
. Vapor leaves at a rate per unit area J. The removed liquid contracts the height h(r)
vertically, vacating the vertically striped region in a short time. The volume of this striped region is equal to the volume
removed by J. But in the shaded annular region the red striped volume is smaller than the volume removed by J there
(red arrows). Thus liquid flows outward to supply the deficit volume: fluid at r sweeps out the blue-striped region in
time. Its volume is the deficit volume; its depth-averaged speed is v(r) is proportional to J.
I finally got the outward flow inside the droplet (figure 1.10A,B). The velocity was determined in
each point of the domain (figure 1.10A). It was considered a droplet of ~ 0.6 l with a CA of 40°.
I defined another geometry and mesh for the case of a droplet on SHS (figure 1.11 A,B) and
solved the Laplace and Navier-Stokes equations for this case. I considered a 5 l droplet with a CA
of ~155°. As boundary conditions, I set the saturated vapor concentration to the free surface
V
C and
the concentration value of the ambient
V
C H C
at a distance equal to 20 times the radius
interface that is considered the limit of the outer domain to the droplet. The substrate is not
considered absorbent and the vapor flux is set to 0.
Ricevi informazioni sui nostri servizi, sulle offerte e non perdere news e consigli su università e lavoro.
