7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
The impact of inkjet droplets on a paperlike structure
M. DoQuang* and A. Carlson* and G. Amberg*
Linn6 Flow Centre, Department of Mechanical, Royal Institute of Technology, Stockholm, 100 44 SWEDEN
minh@mech.kth.se
Keywords: Multiphase flow, wetting surface, droplet impact, CahnHilliard equation
Abstract
Inkjet technology has been recognized as one of the most successful and promising microsystem technologies. The
wide application areas of printer heads and the increasing demand of high quality prints are making ink consumption
and print seethrough important topics in the inkjet technology. In the present study we investigate numerically
the impact of ink droplets onto a porous material that mimics the paper structure. The mathematical framework is
based on a free energy formulation, coupling the CahnHilliard and Navier Stokes equations, for the modelling of the
twophase flow. The case studied here consists of a multiphase flow of airliquid along with the interaction between a
solid structure and an interface. In order to characterize the multiphase flow characteristics, we investigate the effects
of surface tension and surface wettability on the penetration depth and spreading into the paperlike structure.
Introduction
Inkjet technology has been recognized as one of the
most successful and promising microsystem technolo
gies. The wide application areas of print heads and the
increasing demand of high quality prints are making ink
consumption and print seethrough important topics in
the inkjet technology. In order to advance to a cheaper
production of high quality prints, fundamental issues
about the physics concerning the impact of ink droplets
on paper structures needs to be resolved.
Many phenomena with a complex physics takes
place from a droplet are ejected from the print head,
until it has obtained its equilibrium form in the paper
structure. In the first stage the droplet shoots out of the
print head with a velocity typically ranging between
1 20 m s1 and it will often, as it travels towards
the paper, form a tail that consist of small satellite
droplets (DoQuang et al. (2010)). It impacts the
paper, which has a heterogeneous surface of cellulose
fibers. These fibers are chemically treated, often with
different degrees of wettability. After the droplet
has impacted onto the surface it can either penetrate
or recoil, depending on the papers roughness and
wettability. Experiments by Kannangara et al. (2006)
showed that the forced spreading and droplet recoil after
impact on commercial paper surfaces depended on their
wettability. The surfaces were found to inherit a dual
nature, and behaved as hydrophobic upon first contact
with the impacting droplet allowing no ink to penetrate.
Thereafter, during the droplet recoil process it behaved
like a hydrophilic surface. If the droplet fully recoils
or if it is sufficiently large it can make a splash so that
the ink spreads over a larger part of the paper, thus
reducing the prints quality. How to trigger the splash
of a droplet has been a matter of intense investigation
in the literature (Yarin (2006)), where the wettability,
surface structure (Xu et al. (2007)) and environmental
pressure (Xu et al. (2005)) have all been identified as
key parameters to trigger or suppress splashing.
After the droplet has impacted onto the paper, wetting
will dominate the droplet infiltration into the paper
structure. Since the surface often has a heterogeneous
structure, consisting of fibers with different wetting
properties, that also adsorbs the liquid poses additional
challenges to the already complex wetting physics that
one observed on relatively smooth solid surfaces see
Blake (2006). Modaressi and Gamier (2002) found
in experiments that the droplet evolve into the paper
as a function of two sequential phenomena. First, the
droplet spreads into the material, forming its footprint
in the paper, until it reaches its pseudoequilibrium
contact angle. After the droplet has reached its pseudo
equilibrium state it starts to adsorb into the bulk of the
paper material.
Numerical simulations of the impact of a droplet
on a porous surface were performed by Alam et al.
(2007) with the Volume of Fluid method, Hyvaluoma
et al. (2006) with the Lattice Bolzmann method and by
Reis et al. (2004) with a markercell method. These
have both in common that the droplet size was much
larger than the characteristic roughness of the porous
media. Alam et al. (2007) examined the effect of surface
structures and found that a sustained pressure outside
the porous media increased the adsorption depth as a
function of time.
Here, we adopt the Phase Field method to numeri
cally investigate the impact of an inkdroplet onto a
paperlike structure. The case studied here consists of a
multiphase flow of airliquid along with the interaction
between a solid structure and an interface. We focus on
the initial regime, before adsorption of the liquid into
the bulk surface material, and seek to characterize the
pseudoequilibrium regime as reported by Modaressi
and Gamier (2002). A small droplet with the same size
as the characteristic surface roughness is considered
as it impacts a web of cellulose fibers, mimicking the
paper structure. By only changing the fibers wettability
we show that the droplet can either penetrate or bounce
as it impacts the paperlike structure.
Governing equations
Several authors have previously demonstrated the appli
cability for the diffuseinterface method to describe two
phase flows (Anderson et al. (1998); Jacqmin (1999)).
DoQuang and Amberg (2009) and DoQuang et al.
(2010) has demonstrated the capability of this method
for the simulations of liquidgas systems. Here, we will
briefly describe the main ideas and list the governing
equations, for a mixture of two Newtonian fluids.
In the phasefield model, the order parameter or
phasefield 9, is has a distinct equilibrium value rep
resenting the two phases, but it changes rapidly but in
a smooth fashion between the two equilibrium states
across the interface. Here takes the value = 1 in liquid
phase and 4 1 in gas phase. The free energy of the
system is described by a GinzburgLandau expansion of
the free energy of the system (Cahn and Hilliard (1958)),
= / (aq) 'V ) d + gl(j (a)dF, (1)
where a and B are constants that are related to the sur
face tension and interface thickness. y(9) represents
here the bulk energy and takes the form of a double
well potential function, with two minima = 1 corre
sponding to the two stable phases. 7 is here represented
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
by,
S( + 1)( 1)2. (2)
4
The second term in equation (1) describes the inter
face energy. This term associates with variations of the
phase field 4 and contributes the freeenergy of the in
terfacial region, which defines the surface tension coef
ficient,
a f /+ ( do 2dx
OO dx
220
3
3V2 
The free energy at the solid surface dF is formulated
by the surface energy contribution from the three inter
faces appearing; solid (s), gas (g) and liquid (1), and
g(9, a) is a function varying smoothly between the sur
face energies as8 and a,, (Jacqmin (1999)).
By taking the variational derivative of the free energy,
F, with respect to the order parameter 4 and perform
some algebra transformations, we obtain the Cahn and
Hilliard (1958) equation,
at
+ (u V)< = V (nV(, '"(4) aV2 )) (4)
In this equation, the interface is not captured by a
sharp interface. It uses, 4, a finite thickness, smooth
transition region to distinguee the different phases. Here,
K is the constant mobility and r is the chemical potential,
defined as
11 P Cv2". (5)
Once the phase field is calculated, the physical prop
erties such as the density and the viscosity are calculated
as follows,
1+9 1 
P =P + P Ps
P P 2 2
1+ 1 
2 + 3 2 '
where pi, p, and p/, p, are the densities and viscosities
of the liquid and gas phase, respectively.
The fluid flow is described by the NavierStokes equa
tions for an incompressible flow.
du
p(9 + (Vu))
at
Vp + V (p(Vu + VTu))
7VO, ((8)
V u 0,
where p denotes the density, u the velocity vector, f the
viscosity, and p the pressure. The last term in equation
(8) is expressing the potential form of the surface tension
force, proposed by Jacqmin (1999).
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Boundary conditions
In phase field theory, the wetting boundary condition for
the interface is set via the balance of the free energy dis
tribution between the different phases. By making the
assumption that the interface is at local equilibrium as it
wets the surface, the boundary condition becomes Vil
lanueva and Amberg (2006),
aVp. n + a cos(0S)g'(p) 0, (10)
where 8, is the static equilibrium contact angle. Here the
contact angle is related to the surface tension coefficients
a through the Young's equation: a cos(0,) = acr sl.
The assumption of local equilibrium at the solid sur
face has been a widespread assumption in phase field
wetting simulations, which has proven to be successful
in describing numerous physical phenomena involving
moving contact lines (Jacqmin (1999); Villanueva and
Amberg (2006); DoQuang and Amberg (2009); Do
Quang et al. (2010)). Recently, Carlson et al. (2009)
has included the dissipative mechanism into the bound
ary condition of the Phase Field framework. It allows
for the nonequilibrium wetting contact angle in rapid
wetting. Such dissipative effects are assumed to be neg
ligible here, thus applying the assumption of local equi
librium (Villanueva and Amberg (2006); DoQuang and
Amberg (2009)) as the liquid wets the solid surface us
ing the boundary condition given in eq.(10).
Nondimensionalization
The governing equations are made dimensionless based
on the characteristic parameters of the flow, giving the
dimensionless variables,
I X U tu, P
S= = u, Lt p = (11)
L, L, P U 21
Lc Uc Lc Pc c
where Lc is the characteristic length taken to be the
droplet radius, Uc is the characteristic velocity taken to
be the initial velocity of the ink droplet. pc is the char
acteristic density defined as the water density. Dropping
the primes, the dimensionless equations are
Du
DL
V u
Do
DL
Vp + 1 V. (p(y)(Vu+ VTU))
Re
1+V^ (1
+ Ca. Cn.RcnV, (L,
1V
1V (Vl),
Pe
SC' V2 p. (15)
Note that incompressibility does not imply that the
density is constant, only that the density is indepen
dent of pressure, which is a good approximation when
ever flow speeds are small compared to speeds of sound.
Also, note that the Peclet number in eq.(14) is large, due
to the small value of the diffusion coefficient. Eq.(14)
then essentially states that y, and thus density, is con
stant along a streamline, which is consistent with the as
sumption of incompressibility in eq.(13).
The dimensionless parameters are the Capillary num
ber Ca, Reynolds number Re and Peclet number Pe and
Cahn number Cn,
UL, pcUcLc cUc
Pe = OCn = Re = ,Ca =
D Lc R e
(16)
where pc is the characteristic viscosity taken to be the
liquid ink viscosity, D is the difusivity of liquid vapour
in air, = /a/ is the interface thickness. The Peclet
(Pe) number expresses the ratio between advection and
diffusion. The Cahn (Cn) number expresses the ratio
between the interface width and the characteristic length
scale. The Reynolds (Re) number expresses the ratio
between the inertia and the viscous force. The Capil
lary (Ca) number expresses the ratio between the vis
cous and the surface tension force.
Numerical treatment
The numerical simulations were carried out using fem
Lego (Amberg et al. 1999), a symbolic tool to solve par
tial differential equations with adaptive finite element
methods. The partial differential equations, boundary
condition, initial conditions, and the method of solv
ing each equation are all specified in a Maple work
sheet. The CahnHilliard equation is treated as a cou
pled system for the chemical potential n and the com
position p. Both the chemical potential and the com
position equations are discretised in space using piece
wise linear functions and discretised in time using an
implicit scheme. The coupled nonlinear algebraic sys
tem of r and p is solved by an exact Newton's method.
Within each Newton iteration, the sparse linear system
is solved by unsymmetric multifrontal method (UMF
PACK), Davis (2004).
To ensure mesh resolution along the vicinity of the
interface, an adaptively refined and derefined mesh is
used with an adhoc error criterion function,
e V2 < tol.
J^,
The implementation of the mesh adaptivity can be de
scribed as follows. At each mesh refinement step, an
element fIk is marked for refinement if the element
size is still larger than the minimum mesh size allowed,
h > hmi,, and it does not meet the error criterion (17). c
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
is an ad hoc parameter. In the case that an element meets
the error criterion, it is marked for derefinement unless
it is an original element. At the next refinement step, el
ements containing hanging nodes are marked for refine
ment. The refinement/derefinement stops if and only
if no element is marked for refirement/derefinement.
More details about this scheme can be found in Vil
lanueva and Amberg (2006); DoQuang et al. (2007);
DoQuang and Amberg (2009).
The NavierStokes equations are solved using a pro
jection method for variable density that was introduced
by Guermond and Quartapelle (2000). The Navier
Stokes equations are also discretized in space using
piecewise linear functions with the convective term
treated as a semiimplicit term which allows a longer
time step in the computations. The linear system
is solved by the generalized minimal residual method
(GMRES).
Numerical results and discussion
The performance and convergence of the method was
tested on different problems where the motion was
driven by surface tension. The gravitational forces are
supposed to be small in our future applications. The
gravity was therefore set to zero in all the computations
presented here.
The Laplace law
We have measured the pressure jump for different mesh
spacing and Cn numbers. The Cn number gives the ra
tio between the width of the diffuse interface and the
characteristic length scale in the flow, here being the
droplet diameter d. The results are summarized in ta
ble(l), where we have kept the Ca 1 Pe 3 103
and Re = 1 fixed. These dimensionless numbers
gives an analytical pressure difference (AP)analytical =
8/2/3. The numerical domain has an extension of
[2d x 2d x 2d] and an equidistant mesh has been applied.
Table(l) is summarizing the relative error between the
computed and analytical pressure prediction for different
Cn numbers and mesh spacings, after eight time steps.
It is noted that the correct pressure is immediately ob
tained with good agreement between the numerical and
analytical solution. One trend in table(l) is that the error
in pressure depends on the numerical resolution of the
interface. Another observation is that the correct pres
sure jump is obtained even with wide interfaces.
Droplet oscillations
The dynamic behavior of the surface tension model has
been verified by validating the numerical simulations
Cn
0.015
0.04
0.04
0.06
0.08
Ax 0.003 0.013 0.02 0.02 0.027
Perror 0.06% 0.6% 2.0% 0.6% 0.9%
Table 1: Deviation between the numerical and analyt
ical pressure for different Cn numbers and
mesh resolutions. Ax is the mesh spacing and
Perror is defined as the relative error between
the analytically and the numerically predicted
pressure jump, P.. = 100. (1 (AP)an a ).
1 015
1 01 
1 005
0 995
099
0985
0 10 20 30 40 50
Figure 1: Evolution of the droplet radius.
against an analytical expression for droplet oscillations
in the absent of gravity. The droplet has a density pl and
viscosity pi submerged in an external fluid with a den
sity P2 and viscosity P2. In cylindrical coordinates the
droplet radius is given by
r = Ro (1 /4 + (P,(cos0)) (18)
where Ro is the initial droplet radius, P, is the Legendre
polynomial of order n, and < 1. Fyfe et al. (1988) ex
tended the linear Rayleigh's theory for small amplitude
oscillations on cylindrical jets and introduce an analyti
cal expression for infinitesimal amplitude oscillations of
an incompressible, inviscid droplet. The frequency c for
the droplet oscillation is given by
n n
2 .n (19)
P1 + p2 Ro
Several simulation were performed with different den
sity ratios, where the nondimension parameters have
been kept constant as; Re = 200, Ca = 0.01, Ro = 1,
n 2, = 0.01. The evolution of the radius in time
is shown in fig. 1. Table 2 shows the analytical result
for the oscillation frequency (eq.(19)), numerical solu
tion and their relative error for different density ratios.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Figure 2: Impact droplet on a paperlike structure with the wetting contact angle 0 = 300. (A) at time reference
t = 0 s, (B) t = 17 s, (C) t =34 ps and (D) t =85 ps
P1lP2
1
0.1
0.01
0.001
Analytical w
1.1892
1.6035
1.6734
1.6810
Numerical w
1.1758
1.5921
1.6644
1.6790
Error %
1.13
0.71
0.54
0.12
Table 2: The analytical and numerical oscillation fre
quency (w) and their relative error for different
density ratios pl/p2.
It should be noted that the conservation of mass m
Jf pdV is recorded in time, with only a 6.339 x 10 6 %
mass variation from the initial to the final time.
The impact of the inkjet droplets
Fig.(2, 3) shows two events as the ink droplet impacts
the paperlike structure. The cellulose fibers are mim
icked in the simulation by a set of solid circular rods,
see fig.(2, 3). A nonevaporating droplet with a diameter
of R 27 pm impacts onto the two surfaces and it has
an initial speed, o = 6 m s 1. The droplet has a
density p 880 kg m3, viscosity p 0.01 Ns m 2
and surface tension coefficient a 0.032 N m1
Fig.(2) shows the temporal evolution of the droplet
shape as it impacts a surface of cellulose fibers with
an equilibrium contact angle 0 30. As the droplet
interacts with the fibers a thin air layer separates the
solid surface and the droplet interface. The air layer
diffuses into the droplet, and the interface wets the
solid, see Fig.(2b). A liquid jet propagates into the
paper structure, which finally touches the second layer
of fibers (Fig.(2c)). Immediately, the interface forms a
contact line with a large apparent contact angle, and the
liquid spreads across the fiber. The liquid that penetrates
through the two fibers continues to spread onto the
papers top fiber layer. This along with the spreading
on the second layer, as the droplet relaxes toward its
equilibrium shape, squeezes the liquid through the film
formed in between the two layers (fig.(2c)) resulting
in a film breakup Fig.(2d). A small secondary drop is
deposited on the fiber in the second layer and the rest of
the droplet forms a footprint on the top layer covering
the gap formed between the two horizontally adjacent
fibers, see Fig.(2d).
By changing the surface chemistry of the fibers
they might become less wettable, resulting in a very
different droplet characteristics. Fig.(3) shows the
droplet shape at four snapshots in time, as it impacts a
paper consisting of fibers with an equilibrium contact
angle 0 180. First the droplet takes a similar shape
in the structure, as observed on a more wettable surface
Fig.(2b). However, the air layer separating the solid the
interface is retained as it is favored by the surface. As
a consequence, the liquid will not wet the second layer
in the structure and surface tension contracts the liquid
film into a energy minimizing shape, see Fig.(3c). On
the top layer the droplet spreads in a similar fashion as
observed on a flat surface, where a liquid rim is formed.
A second phase in the impact dynamics takes then place
as the droplet has decelerated on the surface, as its
initial inertial energy is converted into surface energy.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Figure 3: Impact droplet on a paperlike structure with the wetting contact angle 0 = 1800. At time reference t = 0 Is
(A), 17 ps (B), 34 ps (C) and t = 85 ps (D).
The capillary force then contracts the droplet into its
energy minimizing shape, so that the droplet bounces
back leaving no ink on the paper surface, see Fig(3d).
Conclusions
In the present paper we have presented twodimensional
simulations of ink droplets impacting on a paperlike
surface. By only changing the wettability of the
cellulose fibers droplets were either deposited on the
paper or could bounce off. This identifies the wettability
to be key parameter in order to obtain high quality prints.
The droplet spreading on a highly wettable surface
(0 30) shows that secondary droplets can be
deposited on the second fiber layer. The droplet was
found to significantly decelerate as it penetrated the web
of fibers, in compare with the results on a hydrophobic
surface (0 = 180). This points out that energy is
stored at the contact line, indicating that the presence
of contact lines could even influence macroscopic
parameters like the friction coefficient in porous media.
Acknowledgements
Computer time provided by SNIC (Swedish National In
frastructure for computing) is gratefully acknowledged.
References
P. Alam, M. Toivakka, K. Backfolk, and P. Sirvio. Im
pact spreading and absorption of newtonian droplets on
topographically irregular porous materials. CHEMICAL
ENGINEERING SCIENCE, 62:31423158, 2007. ISSN
00092509.
Gustav Amberg, Robert TOnhardt, and Christian Win
kler. Finite element simulations using symbolic com
puting. Mathematics and Computers in Simulation, 49
(45):257274, 1999.
DM Anderson, GB McFadden, and AA Wheeler.
Diffuseinterface methods in fluid mechanics. Annual
Review of Fluid Mechanics, 30:139165, 1998.
Terence D. Blake. The physics of moving wetting lines.
Journal of Colloid and Interface Science, 299(1): 113,
2006. doi: DOI 10.1016/j.jcis.2006.03.051.
John W. Cahn and John E. Hilliard. Free energy of a
nonuniform system. i. interfacial free energy. The Jour
nal of Chemical Physics, 28(2):258267, 1958.
Andreas Carlson, Minh DoQuang, and Gustav Am
berg. Modeling of dynamic wetting far from equilib
rium. Physics of Fluids, 21(12):1217014, 2009.
T. A. Davis. Algorithm 832: Umfpack v4.3 
an unsymmetricpattern multifrontal method. ACM
TRANSACTIONS ON MATHEMATICAL SOFTWARE,
30(2):196199, 2004. ISSN 00983500.
M. DoQuang, W. Villanueva, I. SingerLoginova, and
G. Amberg. Parallel adaptive computation of some time
dependent materialsrelated microstructural problems.
Bulletin of the Polish Academy of SciencesTechnical
Sciences, 55(2):229237, 2007.
Minh DoQuang and Gustav Amberg. The splash of a
solid sphere impacting on a liquid surface: Numerical
simulation of the influence of wetting. Physics of Fluids,
21(2):022102, 2009.
Minh DoQuang, Laurent Geyl, GOran Stemme, Wouter
van der Wijngaart, and Gustav Amberg. Fluid dynamic
behavior of dispensing small droplets through a thin liq
uid film. Microfluidics and Nanofluidics, 2010. doi:
10.1007/s104040090547x.
D. E Fyfe, E. S Oran, and M. J Fritts. Surface tension and
viscosity with lagrangian hydrodynamics on a triangular
mesh. Journal of Computational Physics, 76(2):349
384, 1988.
JL Guermond and L Quartapelle. A projection fem for
variable density incompressible flows. Journal of Com
putational Physics, 165(1):167188, 2000.
J. Hyvaluoma, P. Raiskinmaki, A. Jasberg, A. Koponen,
M. Kataja, and J. Timonen. Simulation of liquid pene
tration in paper. PHYSICAL REVIEW E, 73(3):036705,
Mar 2006.
D Jacqmin. Calculation of twophase navierstokes
flows using phasefield modeling. Journal of Compu
tational Physics, 155(1):96127, 1999.
D Kannangara, HL Zhang, and W Shen. Liquidpaper
interactions during liquid drop impact and recoil on pa
per surfaces. Colloids and Surfaces APhysicochemical
and Engineering Aspects, 280(13):203215, 2006.
H Modaressi and G Gamier. Mechanism of wetting and
absorption or water droplets on sized paper: Effects of
chemical and physical heterogeneity. Langmuir, 18(3):
642649, 2002.
NC Reis, RF Griffiths, and JM Santos. Numerical sim
ulation of the impact of liquid droplets on porous sur
faces. Journal of Computational Physics, 198(2):747
770, 2004.
W. Villanueva and G. Amberg. Some generic capillary
driven flows. International Journal of Multiphase Flow,
32(9):10721086, 2006.
L Xu, WW Zhang, and SR Nagel. Drop splashing on
a dry smooth surface. Physical Review Letters, 94(18),
2005.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Lei Xu, Loreto Barcos, and Sidney R. Nagel. Splashing
of liquids: Interplay of surface roughness with surround
ing gas. Physical Review E, 76(6), 2007.
AL Yarin. Drop impact dynamics: Splashing, spreading,
receding, bouncing... Annual Review of Fluid Mechan
ics, 38:159192, 2006.
