- Email: [email protected]
Leveling of thixotropic liquids
J. Non-Newtonian Fluid Mech. 166 (2011) 395–403
Contents lists available at ScienceDirect
Journal of Non-Newtonian Fluid Mechanics journal homepage: www.elsevier.com/locate/jnnfm
Leveling of thixotropic liquids S. Livescu ∗ , R.V. Roy, L.W. Schwartz Mechanical Engineering Department, University of Delaware, Newark, DE 19716, USA
a r t i c l e
i n f o
Article history: Received 11 November 2010 Received in revised form 20 January 2011 Accepted 25 January 2011 Available online 23 February 2011 Keywords: Hydrodynamics Leveling Lubrication Thixotropy Thin films
a b s t r a c t In this paper the effect of thixotropy in the hydrodynamic behavior of thin films is studied. The simple problem of leveling on a horizontal substrate is considered. The rheological properties of the material are assumed to evolve over time due purely to changes in its internal structure. These changes are modeled in terms of a single structural variable. Neither elastic nor yielding effects are taken into account. More specifically, two distinct rheological models are considered: the simple model proposed by Moore and the more complex model proposed by Baravanian et al. These models exhibit a large range of variation for the liquid viscosity across the film thickness. After deriving the hydrodynamic equations governing leveling flows with the standard assumptions required by the lubrication approximation and running time-dependent numerical simulations, the nonlinear leveling history of the liquid can be predicted as a function of the initial microstructural state, rheological parameters, and initial disturbance of the liquid free surface. The main effort of this work is devoted to devising approximation schemes which lead to significant simplifications of the governing equations and their numerical computations. By approximating the inverse of viscosity as a monotonic function between its substrate and free-surface values, excellent agreement is found for the film amplitude, irrespective of the values of the rheological parameters of both models. Finally, a linear analysis yields a generalization of the Orchard’s law of leveling for Newtonian liquids to take into account the effect of thixotropy. © 2011 Elsevier B.V. All rights reserved.
1. Introduction Many industrial processes require coating of solid surfaces with thin liquid films, for both decorative and protective reasons. When a liquid film is applied to a flat surface, the initial coating is generally uneven. It is usually desired that the final coating is as uniform as possible. Many of these liquids such as paints, inks, and coatings have thixotropic rheology as their viscosities change with shear. Considering the size of the industries and processes for which the leveling flow of thixotropic liquids is important, it is perhaps surprising that numerical modeling of this type of flow is still in an incipient phase. The hydrodynamic leveling of thin liquid films has been studied extensively since 1920, most particularly for Newtonian liquids (see Dodge [1] and Quach [2], and references therein). Orchard [3] was the first to treat this problem theoretically. He considered the flow under viscous, capillary and body forces. In the limit of small Reynolds number, Orchard used a Fourier series to solve the twodimensional Navier-Stokes equations governing a liquid film with
∗ Corresponding author. Present address: ExxonMobil Upstream Research Company, 3120 Buffalo Speedway, Houston, TX 77056, USA. Tel.: +1 713 431 4578; fax: +1 713 431 6392. E-mail address: [email protected] (S. Livescu). 0377-0257/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.jnnfm.2011.01.010
a sinusoidal free-surface. His linear model predicts the rate of leveling, the time evolution of the free-surface amplitude, and the liquid stress distribution. Orchard obtained good agreement between his linear theory and the numerical solution of the nonlinear evolution equation governing the film thickness. Thin film hydrodynamics studies taking into account nonNewtonian effects such as shear-thinning [4] or visco-plasticity [5] have also received attention. Iyer and Bousfield [4] developed a numerical technique to include a generalized Newtonian rheology in thin-film leveling equations. Their shear-thinning rheology follows Carreau’s model. Their numerical solution of the onedimensional evolution equation for the film thickness is within 2% of the results obtained by numerical simulation of the full twodimensional problem using finite-element techniques. They also performed experiments whose results were in good agreement with their numerical simulations. Keunings and Bousfield [5] studied the effect of viscoelasticity on surface tension driven leveling of thin films. They generalized Orchard’s linear analysis by incorporating viscoelastic rheology with a generalized Maxwell model. Their results show that elasticity retards the leveling process. They also found that, at early times, the nonlinear evolution of the free surface perturbation is similar to that obtained in the Newtonian case by Kheshgi and Scriven [6]. Thixotropy is usually defined as the reversible variation of viscosity with time, whether elastic effects are present or not (see
396
S. Livescu et al. / J. Non-Newtonian Fluid Mech. 166 (2011) 395–403
Barnes [7], Mewis [8] and Mewis and Wagner [9]). According to Barnes [7], thixotropic liquids have gel-like properties which disappear when sheared, but reappear when put to rest. The rheological properties of such liquids depend on the time needed by their microstructure to evolve from one state to another. Structural changes in a flowing thixotropic liquid are due to two competing effects: break-down due to flow stresses and build-up due to collisions of the particles which make up the microstructure. Thixotropy has important consequences on the flow of liquid films. For instance, during brushing of a thixotropic paint, the viscosity should be sufficiently low in order to obtain good spreading. The paint must subsequently possess an optimal rate of structural build-up for good leveling: if the build-up rate is too slow, the viscosity remains low and the paint will have good leveling but will tend to quickly sag; if the build-up rate is too high, the viscosity will become large, sagging will be minimal but leveling will be poor. Another example was presented by Dodge [1] for two latex (highly nonlinear viscosity variation in time) and alkyd (linear viscosity variation in time) exterior paints. The viscosity behavior of any of these two paints cannot accurately be predicted by assuming that the paint is Newtonian. Dodge noted that the degree of viscosity break-down during paint application and the speed the viscosity increase following application will determine the extent of leveling, for a given film thickness and wavelength of brushmarks. There are very few published hydrodynamic studies of thixotropic liquids (see Billingham and Ferguson [10], Cheng and Evans [11], Huynh et al. [12] and Mujumdar et al. [13]) and we know of only one reference which considers the leveling of thixotropic liquids: Cohu and Magnin [14] developed a computational method by stepwise calculations based on Orchard’s analysis in order to incorporate thixotropic rheology with a five-parameter model. Their stepwise method has the limitation that the shear stress is homogeneous within the liquid. However, their numerical calculations of the leveling rate are in good agreement with their experiments. Our goal is to incorporate thixotropic rheology into the hydrodynamic equations governing the flow of a thin film on a horizontal substrate. For this, we consider two different rheological models: Moore’s model [15] and the model of Baravanian et al. (referred to as BQP) [16]. Moore’s model is a natural choice for its simplicity. It assumes that the microstructure of a thixotropic liquid is composed of a large number of “links”. The rheological behavior of such a liquid depends on the number of links formed at a given moment in time. Thus, the microstructure is described in terms of a nondimensional scalar parameter, denoted by . This parameter will reflect the structural build-up and break-down states. The state of build-up is defined as the number of links which are formed at any instant divided by the total number of links and its corresponding state of break-down is defined as the number of links destroyed at that given time divided by the total number of links. It follows that will have the value 0 if the structure is completely brokendown (all links are destroyed) and the value 1 if the structure is completely recovered (all links are restored). Generally, the behavior of a thixotropic material is described by constitutive equations consisting of a rate equation and an equation of state, D = f (, ), ˙ Dt
= (, ˙ ), ˙
(1)
where , ˙ and are the shear stress, shear rate and viscosity, respectively. Here, the operator D/Dt represents the standard material derivative. We assume that elastic effects are not present, or that such effects appear much faster than thixotropic effects. Hence we confine ourselves to modeling purely viscous phenomena and we take the viewpoint that the rheological properties of the mate-
rial evolve over time due purely to changes in its internal structure adequately modeled by a single structural variable . In Moore’s model, the time rate of change of is expressed as the difference of two terms: (i) a break-down rate due to shearing expressed as the product of the structural parameter by some power of the shear rate, and (ii) a build-up rate expressed as some power of (1 − ). For simplicity we take both exponents equal to unity, as in the original equation proposed by Moore, D = a(1 − ) − b|| ˙ Dt
(2)
Thus, it is seen that the shear-induced break-down term results in a decrease of , whereas the build-up term leads to an increase of . The constitutive equations are completely specified by defining the viscosity as a linear function of : in the limit of high shear rate, → 0, the viscosity reaches the value = ∞ ; conversely, the viscosity of a fully recovered structure, → 1, takes the value = 0 , = ∞ + (0 − ∞ )
(3)
where ∞ and 0 are the limiting viscosities at very high and very low shear rates, respectively. Hence Moore’s model is a fourparameter model (∞ , 0 , a and b). BQP’s model uses the shear stress in the rate equation 1 D = Dt ta
1−
e
(4)
The equilibrium value e and the relaxation time ta are determined from the following empirical formulas:
e () = 1 +
p −1 c
,
ta () = ta,∞ −
ta,∞ − ta,0 1+
a,0
q
(5)
where c , p, q, ta,∞ , ta,0 , a,0 are rheological parameters. The viscosity varies between the two limiting values, ∞ and 0 , and it is again only a function of the structural parameter , =
∞ (1 − K)
2
,
K =1−
1/2 ∞
(6)
0
This model requires eight material parameters ∞ , 0 , c , a,0 , p, q, ta,∞ , ta,0 . In order to compare these two models, consider a thin film of thickness h0 and assume that the shear stress is distributed linearly across the film: = 0 (1 − z/h0 ), 0 ≤ z ≤ h0 . In fact such a linear profile can be shown to be a direct consequence of the lubrication approximation of thin films (see the next section for the derivation of the linear shear stress profile in the lubrication approximation limit). We can determine the thickness profile of the equilibrium viscosity e for both models. For Moore’s model we find e = 1 + ıe , ∞
e =
1 1 + AM (1 − z/h0 )/[1 + ıe ]
with ı = 0 /∞ 1 and AM = b 0 /a∞ . For BQP’s model we find 1 e = , 2 ∞ (1 − (1 − ı−1/2 )e )
e =
1 1 + [ABQP (1 − z/h0 )]
p
with ı = 0 /∞ and ABQP = 0 / c . Fig. 1 shows the variation of e /∞ across the thickness for the two models and a range of parameters. Striking differences can be noted. Moore’s model leads to a quasilinear profile between the substrate value and the free-surface value. BQP’s model leads to a highly nonlinear profile: for high values of ABQP (or high values of exponent p) the film displays a thin skin of high viscosity on top of a thick layer of low quasi-uniform viscosity. Thus, for a large range of parameters (both AM and ABQP vary by two orders of magnitude in Fig. 1), the two models have very different viscosity behavior. This conclusion emphasizes our
S. Livescu et al. / J. Non-Newtonian Fluid Mech. 166 (2011) 395–403
film. l is the wavelength of a sinusoidal perturbation and h0 is the mean thickness of the liquid film. Assuming h0 l, the standard lubrication approximation can be adopted to simplify the hydrodynamic equations. It has been used very effectively to model thin layer flows of Newtonian and generalized Newtonian liquids (see for example [17] or [18]). The key assumptions made in the lubrication approximation are: (i) the flow is slow and inertial forces are neglected, (ii) the liquid film is thin, and (iii) free-surface slopes are small. Considerable simplification of the governing equations follows from these assumptions. For a liquid film flowing on a horizontal substrate, the mass conservation equation is
100 90 70
ηe/η∞
A M=1
M BQP
80
A M=10
A M=100
60 ABQP = 0.05
50
ABQP = 0.5
40
ABQP = 5
30
397
20 10 0 0
0.2
0.4
0.6
0.8
1
z/h0 Fig. 1. Variation of the viscosity at equilibrium, e , for the two rheological models: Moore’s model (“M”) and BQP’s model with exponent p = 1. A film of constant thickness h0 is assumed to be under a linear shear–stress profile = 0 (1 − z/h0 ).
∂w ∂u + =0 ∂x ∂z
(7)
and the simplified momentum equations in the x and z directions are ∂xz ∂p ∂p = , = −g , ∂x ∂z ∂z
(8)
where p is the pressure and xz is the shear stress and and g are the liquid density and gravitational acceleration, respectively. Furthermore at the liquid free surface (z = h), a normal stress balance is written as p = − ∂ 2 h/∂ x2 , where is the surface tension and the curvature of the free surface has been linearized. The shear stress and shear rate within the film are, respectively (0 ≤ z ≤ h)
xz =
∂3 h ∂h − g ∂x ∂x3
(h − z) ,
˙ xz =
xz ∂u = ∂z
(9)
The horizontal and vertical components of the liquid velocity are obtained by integration of (9) Fig. 2. Sketch of a thin film of mean thickness h0 perturbed by a sinusoidal ripple of wavelength l.
motivation for using these two rheological models in the present study. While Moore’s model is chosen for its simplicity (the viscosity varies quasi-linearly across the film thickness), BQP’s model is an arbitrary choice displaying a complex viscosity behavior (the viscosity has a highly-nonlinear variation across the film thickness). The lubrication approximation takes into account, among other things, the thinness of the liquid film: after averaging across the thickness of the film, one can derive an evolution equation for the film thickness h(t, x) (in a two-dimensional domain). The difficulty in solving a flow problem for a thixotropic liquid stems from the fact that the flow behavior is coupled to the history of the microstructure which must be determined in both time and space. Calculation of microstructure at every location (x, z) of the liquid domain is prohibitive and thus simplifying schemes are pivotal to efficient hydrodynamic calculations that take advantage of the lubrication approximation. In order to assess the reliability of such schemes, we must nevertheless solve the full set of modeling equations numerically. This paper proceeds as follows. In Section 2, we implement the two thixotropic models into the leveling hydrodynamics. We arrive at a set of evolution equations which we compute numerically in Section 3. In Section 4, we develop two simplifying approximations which present different degrees of accuracy for the two rheological models. Finally in Section 5, we generalize the so-called Orchard law for leveling of Newtonian liquids to account for the thixotropy effect. 2. Hydrodynamic model With the geometry shown in Fig. 2, we choose the x-axis along the substrate and the z-axis pointing toward the free surface of the
u=
∂3 h ∂h 3 − g ∂x ∂x
0
z
h− d ,
z
∂u d ∂x
w=− 0
(10)
With the usual no-slip condition at the substrate, u = w = 0 at z = 0, and no-stress free surface condition, xz = 0 at z = h, integrating (7) in the z direction and using the Leibnitz rule and the kinematic condition at the free surface which states that the liquid does not cross the interface, Dh/Dt = w(h), the mass conservation produces ∂ h/∂t = − ∂ Q/∂ x, where
Q =
h
udz = 0
0
h
∂u h (h − z)dz − [u(h − z)]0 = ∂z
0
h
∂u (h − z)dz ∂z (11)
is the flow rate per unit film width. Finally these equations yield the following evolution equation for the film thickness h(x, t) ∂ ∂h =− ∂t ∂x
∂3 h ∂h 3 − g ∂x ∂x
0
h
2
(h − z) dz
(12)
Eqs. (9), (10) and (12) along with the rheological equations (2) and (3) for Moore’s model (Eqs. (4)–(6) for BQP’s model) fully pose the leveling problem once initial conditions are specified and symmetry conditions are imposed at the endpoints x = 0 and x = l. We assume that at t = 0 the structure is uniform in space and characterized by an initial value i and corresponding viscosity i (i ). Also note that the material derivative D/Dt takes the form ∂/∂t + u∂/∂x + w∂/∂z. For a Newtonian liquid at constant viscosity i , one recovers the expression of the hydrodynamic conductivity
h
2
[(h − z) /]dz = h3 /(3i ). For a thixotropic liquid, determina0 tion of this conductivity requires the evaluation of the integral across the film thickness. Thus, unless further approximations are made, the advantage gained by the lubrication approximation is significantly reduced.
398
S. Livescu et al. / J. Non-Newtonian Fluid Mech. 166 (2011) 395–403
All these equations have been developed for liquids without yield stress. If we assume that the liquid behaves approximately as Bingham plastic, the shear strain relation is nonlinear ∂u = 0, ∂z
if |xz | ≤ y = xz − sign (˙ xz ) y ,
if |xz | ≥ y
(13)
where y is the yield stress. This is the simplest idealization of the shear-thinning pseudo-plastic liquid. Other expressions are available in literature (see for example Quach [2], Barnes [7], Mewis [8], Mewis and Wagner [9], Billingham and Ferguson [10], Cheng and Evans [11], Cohu and Magnin [14], and references therein). According to (9), the shear stress varies linearly across the film thickness from xz = 0, at z = h, to xz = h( ∂ h3 /∂ 3 x − g ∂ h/∂ x), at z = 0. For a fixed value of x, there are two possible cases, depending on the maximum value of the shear stress, |xz |z=0 : (i) |xz |z=0 ≤ y ; in this case, u = w = 0 everywhere in the liquid domain due to the noslip condition at the substrate, i.e., the liquid does not flow; and (ii) |xz |z=0 > y ; in this last case, there is a so-called yield surface defined by hy = h −
y ∂h3 ∂3 x − g ∂h ∂x
(14)
assuming that the denominator is non-zero. Below the yield surface, 0 ≤ z ≤ hy , there is shearing. Above the yield surface, hy ≤ z ≤ h, there is a layer of plug flow moving with u = u(hy ) (independent of z). Thus, the flow rate per unit film width is
hy
Q =
udz + 0
h
hy
udz = 0
hy
∂u hy − z dz + h − hy u(hy ) ∂z
With these notations, the dimensionless evolution equation for the film thickness takes the form
∂ ∂h =− ∂t ∂x
∂3 h ∂h − g ∂x ∂x3
hy
−sign (˙ xz ) y 0
hy
0
H
(H − Z)2 dZ ,
F(H) 0
Z
U = 3F(H) 0
W =− 0
Z
∂U d ∂X
∂3 h ∂h − g ∂x ∂x3
hy
0
h−z dz − sign (˙ xz ) y
hy
0
2 x , l
=
∞
Z=
z , h0
H=
h , h0
T=
t , t0
U=
u , u0
W=
(17)
w , w0
where we use the viscosity ∞ of the broken-down structure for the reference viscosity. The relaxation time associated with the leveling of a Newtonian fluid of viscosity ∞ (also known as Orchard’s time scale [3]) was chosen as time scale t0 =
3∞ l4
(18)
(2 )4 h30
The velocity scales u0 and w0 are then found to be u0 =
l = 2 t0 3∞
2 h 3 0
l
,
w0 =
2 h0 u0 l
(20)
and
H−Z ∂ ∂ ∂ +U +W = A(1 − ) − 3B F(H) ∂T ∂X ∂Z
(22)
with ı=
0 , ∞
A = at0 ,
B=
bl 2 h0
Note that Moore’s model becomes Newtonian in two situations: ı = 1 or A = B = 0. For BQP’s model, the dimensionless rheological equations can be stated as follows: 1 (1 − K)
(23)
2
e =
1−
1 1 + [3˛|F(H)|(H − Z)]p
,
e
(24)
Ta = Ta,1 +
Ta,2 − Ta,1
q
1 + 3ˇ|F(H)|(H − Z) (25)
with
1 dz
The procedure for solving numerically Eqs. (12) and (16) is similar once hy is calculated. This study will be the subject of a future publication. We return now to our model without yield stress (12). In order to obtain dimensionless equations, we use the following nondimensional variables: X=
(19)
(21)
where the horizontal velocity component at the yield surface is u(hy ) =
∂3 H ∂H − Bo ∂X ∂X 3
= 1 + (ı − 1)
(16)
H− d ,
∂ ∂ 1 ∂ +U +W = Ta ∂T ∂X ∂Z
hy − z dz + h − hy u(hy )
F(H) =
They are needed to compute the convective terms in the rate equation for the structural parameter . For Moore’s model, the dimensionless rheological equations requires three dimensionless parameters (ı, A, B) and the initial structural parameter i . The equations are
=
(h − z)(hy − x) dz
where Bo = gl2 / is the Bond number. The dimensionless components of the liquid velocity are
(15) and the evolution equation of the film thickness is
∂ ∂H = −3 ∂T ∂X
ı=
0 , ∞
Ta,1 =
K =1−
ta,∞ , t0
1 ı1/2
Ta,2 =
,
˛=
∞ u0 , h0 c
ˇ=
∞ u0 , h0 a,0
ta,0 t0
This model requires seven dimensionless parameters ı, ˛, ˇ, p, q, Ta,1 , Ta,2 and the initial value i . It becomes Newtonian for ı = 1 or Ta,1 = Ta,2 = 0. 3. Time-dependent numerical simulations The equations derived in the previous section show that timedependent numerical simulation of the leveling of thin films requires the calculation of the microstructure at every location (x, z) of the liquid domain. It would be in the spirit of the lubrication approximation to develop a simplifying scheme which would involve only the x space variable by averaging the rheological properties across the film thickness. In order to develop such a scheme and to assess its reliability, we must nevertheless solve the full set of equations numerically with two simplifying assumptions. First we will assume that Bo 1, i.e., the wavelength l of the free-surface perturbations is much shorter than the capillary length (/g)1/2 , i.e., capillarity dominates gravity in the leveling process. This simplification is supported by observations that little difference exists between paints applied on the floor and on the ceiling [2]. We will
S. Livescu et al. / J. Non-Newtonian Fluid Mech. 166 (2011) 395–403
0.1
0.05
num approx
0.08
399
num approx
0.048
0.06
0.046 0.044
0.04
0.042 0.04
1/Λ
B=2
ε
B=0.2
0.02
0.038 0.036 0.034
0.01
B=20
B=0.2
B=2
B=20
0.032 0.03
0.005
0.028 0
20
40
60
80
100
0
0.2
0.4
0.6
T Fig. 3. Moore’s model: Evolution of the dimensionless amplitude of the free surface, , by time-dependent numerical simulation of the full equations (neglecting convection) for the rheological parameters ı = 100 in (21) and A = 0.04, B = 0.2, 2, 20, and initial structure i = 0 in (22). The dashed curves (“approx”) correspond to approximating the inverse of viscosity as a linear function of Z (see Eq. (27)).
=3 0
2
(H − Z) dZ
(26)
is evaluated at the old time step, T(k) . Thus, at each time step, solving (19) requires the solution of a pentadiagonal system of linear equations with (N + 1) unknowns for the liquid thickness at the mid(k+1) points Hi+1/2 , i = 0 . . . N. For this we use the well-known Thomas algorithm for banded systems. The computational effort is of order O(N) rather than O(N3 ) for the solution of a general system of equations solved by Gaussian elimination. At each time step the rate equations (22) and (24) are treated explicitly and the new values of (k+1) are converted into values of (dimenthe structural parameter i,j (k+1)
1.2
0.15
λi=0 λi=0.2 λi=0.5 λi=1
sionless) viscosity i,j via (21) and (23). At each time step the accuracy of the scheme is verified by checking that the total mass of liquid is conserved. An example of the film amplitude (T) evolution is shown in Fig. 3. The rheology follows Moore’s model with the rheological parameters ı = 100 in (21) and A = 0.04 in (22). The initial profile is H(0, X) = 1 + 0 cos X with initial amplitude 0 = (0) = 0.1. Fig. 3 is a semilogarithmic plot showing the effect of increasing parameter B; as expected, as B is increased, the shear effects on the microstructure are stronger, thus leading to smaller viscosity at any given time and improved leveling. Note that in all cases, the leveling is initially strongly non-Newtonian, and eventually asymptotes to an exponential decay, whose time constant corresponds to the viscosity 0 of the fully recovered structure.
0.06 0.04
ε
also neglect the effect of convection in the evolution of . Simple inspection of (2) shows that convection is negligible for Moore’s model if either u0 /l a or u0 /l bu0 /h is assumed, i.e., if the nondimensional groups A and B satisfy A 1 or B 1. Unless specified otherwise, the initial microstructure is assumed uniformly brokendown: i = (T = 0, X, Z) = 0. Eq. (19) and the rheological equations corresponding to each model can be solved using finite-difference methods. Eq. (19) has the character of a higher-order diffusion equation. To obtain the finite difference form of (19) we use a forward difference in time and centered differences in space. If the liquid domain is discretized as (X, Z) → (Xi , Zj ) = (i x, j z) all derivatives are evaluated at the midpoints (Xi+1/2 , Zj ), where Xi+1/2 is the midpoint between Xi and Xi+1 ; i.e., Xi+1/2 = (i + 1/2) x. ∂ 3 H/∂X3 is evaluated at the new time step, T(k+1) , while the term H
1
Fig. 4. Moore’s model: Variation of the inverse of viscosity 1/ across the film thickness at X = /3 and time T = 10 (same rheological parameters as in Fig. 3). The inverse of viscosity is maximum at the substrate (Z = 0) where the shear is maximum.
0.1
0.8
Z
0.02
0.01
0.004 0
50
100
150
200
T Fig. 5. Moore’s model: Effect of initial structure on the evolution of the film amplitude. The rheological parameters are fixed to ı = 100 (21) and A = 0.04, B = 2, and i = 0, 0.2, 0.5, 1 in (22). The more broken-down the initial structure, the better the leveling.
In Fig. 4 the inverse of viscosity variation across the film thickness is shown for the same values of the rheological parameters. The curves are shown at time T = 10 at a fixed horizontal position, X = /3 (chosen to avoid any obvious symmetry). A few things can be noted. First, the inverse of viscosity is nearly a linear function of Z, a property which we can use to devise a simplifying scheme. Second, all curves reach the same value at the free surface; indeed the stress is zero there, and hence the viscosity always recovers monotonically from ∞ (at t = 0) to 0 (as t → ∞) irrespective of the value of B. Conversely, shear effects are maximum at the substrate (Z = 0), the viscosity is smallest there, and the difference between the free-surface and substrate viscosities increases as B is increased. Understandably, the initial state of the microstructure has a significant effect on leveling. In Fig. 5 the evolution of film amplitude is shown for i = 0, 0.2, 0.5, 1 fixing the rheological parameters to ı = 100, A = 0.04, and B = 2 in (21) and (22), respectively. As i is increased, the initial viscosity is increased, therefore the free surface takes more time to level. Corresponding results for BQP’s model are shown in Figs. 6 and 7 for the parameters ı = 100, ˛ = 1.5, ˇ = 3, p = 0.5, q = 4, Ta,1 = 2.5, Ta,2 = 0.1, and initial structure i = 0 in (24) and (25). The solid curves are obtained by numerical simulation of the full equations. The additional curves correspond to the simplifying schemes introduced in the next section. A notable difference between Moore’s
400
S. Livescu et al. / J. Non-Newtonian Fluid Mech. 166 (2011) 395–403
4.1. Linear inverse of viscosity approximation
0.1 num r=0.5 r=1 r(t)
0.06 0.04
The inverse of viscosity variation across the film thickness in the case of Moore’s model (see Fig. 4) suggests the following scheme: approximate the inverse of viscosity as a linear function of Z in terms of its free-surface value 1/ b (Z = H) and its substrate value 1/ s (Z = 0), i.e.,
ε
0.02 0.01
1 1 Z = + H (T, X, Z) s (T, X)
0.006 0.004
1 1 − b (T, X) s (T, X)
(27)
The evolution equation for H now becomes
0.002
∂H 1 ∂ =− 4 ∂X ∂T
H3
0.0008 0
20
40
60
80
100
T Fig. 6. BQP’s model: Evolution of the film’s amplitude for the rheological values ı = 100, ˛ = 1.5, ˇ = 3, p = 0.5, q = 4, Ta,1 = 2.5, Ta,2 = 0.1, and initial structure i = 0 in (24). “Full” numerical simulations, denoted by “num”, are compared with the simplifying approximation with fixed value of the exponent r = 0.5, 1 and with time-varying exponent denoted by “r(t)”. Best agreement is obtained by fixing r = 0.5 or by using a time-varying exponent.
1 b
+
3 s
(28)
The inverse of viscosity evaluated at the substrate and at the free surface can be related to corresponding values of according to s = 1 + (ı − 1)s ,
b = 1 + (ı − 1)b
(29)
where s = s (T, X, Z = 0) and b = b (T, X, Z = H). To find the evolution equations for s and b from that of (22), let us define (T, X) = (T, X, Z = (T, X))
(30)
Upon using the identity (chain rule of differentiation)
0.11
num r=0.5 r=1 r(t)
0.1 0.09
1/Λ
∂3 H ∂X 3
∂ ∂T
=
∂ ∂ ∂ + , ∂T ∂Z ∂T
0.08
(22) becomes
0.07
∂ ∂T
0.06
=
−U
∂ + ∂X
∂ −W ∂T
∂3 H ∂ + A(1 − ) − 3B 3 (H − Z) ∂Z ∂X
(31) Z=
Thus, on the substrate, = 0, the no-slip condition (U = W = 0) implies that the convection terms drop out. This gives the following rate equation for s
0.05 0.04
∂3 H ∂s s = A(1 − s ) − 3BH 3 ∂T 1 + (ı − 1)s ∂X
0.03 0.02 0.01 0
0.2
0.4
0.6
0.8
1
1.2
(32)
At the free surface we have = H(T, X) and the kinematic condition
Z Fig. 7. BQP’s model: Inverse of viscosity across the film thickness at X = /3 and time T = 10. Comparison between the full numerical simulations (“num”), constant r-scheme (r = 0.5, 1) and time-varying r-scheme (“r(t)”). Parameters are identical to those of Fig. 6.
model and BQP’s model is seen in Fig. 7; the inverse of viscosity does not follow a linear variation across the thickness of the film. This property is also used in the next section in devising a simplifying scheme. 4. Simplifying approximations The numerical results obtained for the inverse of viscosity variation across the film thickness for the two rheological models presented in the previous section suggest two simplifying schemes: linear inverse of viscosity across the film thickness and nonlinear inverse of viscosity with time-dependent dependence across the film thickness. The main advantage of these approximations is that the term (26) is obtained in closed form as a product of H3 and a linear sum of values of the inverse of viscosity evaluated at the substrate and at the free surface. These schemes complement the lubrication approximation depending on the rheological model used. In this section, we will present these approximation schemes and compare their results with the “full” numerical simulation results.
∂H ∂H + Ub = Wb ∂T ∂X must hold. Here Ub = U(T, X, Z = H) and Wb = W(T, X, Z = H) are the components of the free-surface velocity. Upon using the identity ∂ (T, X, H(T, X)) = ∂X
∂ ∂H ∂ + ∂X ∂Z ∂X
Z=H
the following rate equation for b is found ∂b ∂ = −Ub b + A(1 − b ) ∂X ∂T
(33)
The method of characteristics can be used to show that if b (T = 0, X) = i is independent of X, then b (T, X) ≡ b (T) must remain independent of X, i.e., convection plays no role in the evolution of b . The characteristic curves corresponding to (33) are described by {(X(S), T(S)), 0 ≤ S < ∞ }, where dT = 1, dS
dX = Ub , dS
db = A(1 − b ) dS
(34)
Thus, if the structure is initially homogeneous, b varies in time in the same manner for all points of the free surface. Therefore, no convective effects are present in the equations describing the evolution of the substrate and free-surface structural parameters. Convective effects in the bulk of the liquid have not been neglected; yet their effect is made simpler because of the particular properties at the boundaries of the liquid layer.
S. Livescu et al. / J. Non-Newtonian Fluid Mech. 166 (2011) 395–403
Comparison between the results obtained by “full” numerical simulations and the solution of the coupled system of Eqs. (28), (32), and (33) are presented in Figs. 3 and 4 for Moore’s model. The amplitude agreement is very good, even for large values of the break-down parameter B. The agreement for the inverse of viscosity is also good. However, as B is increased, the inverse of viscosity becomes slightly nonlinear. This nonlinearity suggests that the linear inverse of viscosity approximation does not hold for a wide range of parameters or for other rheological models. In fact, the linear inverse of viscosity approximation did not obtain good agreement with the “full” numerical simulations obtained for BQP’s model, for which the inverse of viscosity has a strong nonlinear variation across the film thickness.
fluid domain is set equal to zero
2
H
dX 0
0
∂ −L ∂T
1
dZ = 0
(39)
where we have fixed the weight function to unity. The following equation can then be derived a0 (r)
dr = a1 (r)(r + 1) + [a2 (r) + ac (r) − aL (r)] (r + 1)2 dT
2
a0 =
H
1
0 2
−
s
a1 (r) = 0
1 b
dX
1
(41)
1 − s b
∂H ∂T
+H
∂ ∂T
1 s
r
1 1 − s (T, X) b (T, X)
(35)
Using this nonlinear Z-dependence of the inverse of viscosity leads to the evolution equation of the film thickness ∂H 1 ∂ =− r + 3 ∂X ∂T
H3
∂3 H ∂X 3
r b
+
3 s
r(T )
1 1 − s (T, X) b (T, X)
(37)
1
= −U
∂ ∂ −W + L1 ∂X ∂Z
1
≡L
1
(38)
where the spatial operator L1 is dependent on the rheological model used. The governing equation for r(t) is obtained as in the standard Galerkin method. If the inverse of viscosity (37) is substituted in (38) one obtains a non-zero residual. To make this residual as small as possible, the integral of the weighted residual over the
1
∂ ∂T
b
2
H
dX 0
aL (r) =
0
2
dX 0
This scheme would have the merit of using an exponent r which would vary in time in accordance with the flow and rheology evolution. In a more complex flow problem over a two-dimensional substrate, the exponent r can be made a function of time t and of the substrate variables (x, z). In a leveling flow, it is reasonable to assume that the exponent r is spatially uniform. In problems with significant spatial heterogeneity along the surface of the film, this assumption is unlikely to be acceptable. First note that the equation governing the non-dimensional inverse of viscosity 1/ can always be put in the form ∂ ∂T
0
The question then becomes how to choose exponent r. The numerical results presented in Figs. 6 and 7 obtained for the BPQ model show that a “static” value of r can be found that best fits the “full” numerical simulations (note that equations for the evolution of the boundary values b and s can be derived for the BPQ model in a similar manner). For instance r ≈ 0.5 is found to lead to the best agreement for the particular parameters chosen in Figs. 6 and 7. But this is not particularly appealing; although this approximation obtains more accurate results than the linear inverse of viscosity approximation, a scheme which would determine the optimal value of the parameter r without resorting to “full” numerical simulations must be instead developed. One solution is to treat the exponent r as time-varying and to find the ordinary differential equation that governs its evolution in a systematic way. Thus, the inverse of viscosity can be written as 1 Z 1 = + 1− H (T, X, Z) b (T, X)
H
ac (r) =
(36)
2
a2 (r) =
∂ − ∂T
1 b
dX
(42)
The linear inverse of viscosity approximation (27) could be improved by generalizing it to the nonlinear form 1 1 Z = + 1− H (T, X, Z) b (T, X)
(40)
where
4.2. Time-dependent nonlinear inverse of viscosity approximation
401
(43)
∂ ∂ +W U ∂X ∂Z
H
L1 0
dX
1
dZ
1
dZ
(44)
(45)
Note that a1 (r) = da0 /dT. Coefficients a0 , a1 and a2 are determined numerically once the functions H, s and b are known at any given time. Coefficients aL and ac can be simplified and their expressions for Moore’s model can be found in Appendix A. We apply this scheme by starting with Moore’s model. We choose the parameters ı = 100 in (21) and A = 0.04, B = 20, and initial structure i = 0 in (22) for which a slight discrepancy is seen when comparing the results of the linear inverse of viscosity approximation with the “full” simulations (see Fig. 3). Although the numerical results were quite satisfactory with the static value r = 1 (the static nonlinear approximation reduces to the linear approximation in this case), Moore’s model constitutes a good test case. The initial value r(0) = 1 is chosen, although, in principle any initial value can be taken for an homogeneous initial structure. The errors for the amplitude history along with the corresponding inverse of viscosity 1/ and time-evolution of r are reduced. The maximum error for the inverse of viscosity (at the substrate) is less than 0.2%. The evolution of r shows a rapid initial jump followed by a slow decrease to the asymptotic value r = 1. The chosen initial value for r(t) has little effect on its time evolution, apart for a very short time interval. Also the effect of convection (which we could “turn off” by setting ac = 0) is negligible, unsurprisingly for the chosen value of B. These results confirm the fact that the linear approximation first chosen to approximate 1/ was a reasonable “guess”. We now turn to the BQP model for which the linear approximation is inadequate. The equilibrium curves shown in Fig. 1 also attest to the need for an adjustable exponent r with the rheology. Figs. 6 and 7 show the results of the r(t)-scheme starting with the initial value r(0) = 0.5. Again, r exhibits a sharp increase and then a sharp decrease near the value r = 0.5 followed by a slow decrease towards r = 0.5. Varying the initial value of r gives the same behavior. The agreement with the “full” simulations is quite good. These calculations justify the conclusion that the static value r = 0.5 yields the best results. The r(t)-scheme of course has the advantage of being self-contained: it calculates the value of r which is best-suited for the flow/rheology of the problem at any given time. By turning off again the convection effect (by setting ac = 0), we found that its effect was always negligible for the cases considered here. The
S. Livescu et al. / J. Non-Newtonian Fluid Mech. 166 (2011) 395–403
success of the method is further demonstrated when the rheological parameter p in (25) is varied: the inverse of viscosity exhibits a wide range of variation across the film thickness which is well reproduced by the approximate scheme. Increasing the exponent p leads to higher substrate viscosity at any given time (the freesurface values remain identical), and thus to poorer leveling. Of the other parameters in (25), ˛ has the strongest effect on leveling; increasing ˛ leads to a decrease of the substrate viscosity, and thus to faster leveling. Note that in this study the two simplification schemes have been verified against the “full” numerical results for thixotropic liquids without yield stress. For liquids with yield stress, both approximations can be applied in a similar manner. For example, considering Bingham plastic liquids for which the film thickness evolves according to (16), the viscosity values at the substrate and at the yield surface are needed in the two approximation schemes. Thus, the linear inverse of viscosity approximation becomes 1 1 Z = + Hy (T, X, Z) s (T, X) if 0 ≤ Z ≤ Hy =
1 1 − y (T, X) s (T, X)
1 , if Hy ≤ Z ≤ H y (T, X)
,
1
-1/ε dε/dT
402
B=0.2 B=2 B=20
0.1
0.01 0
20
40
60
1.5 1
(46)
δ=1 δ=10 δ=100
1−
Z Hy
r(T ) ·
1 1 − s (T, X) y (T, X)
1 , if Hy ≤ Z ≤ H y (T, X)
,
(47)
5. Linear analysis
H(T, X) = 1 + (T ) cos(X)
( 1)
(48)
where the amplitude decays linearly according to d = − dT
(49)
Hence the film levels exponentially fast with the time relaxation given by Eq. (18) for a viscosity = ∞ Our goal here is to generalize Eq. (49) for thixotropic liquids obeying Moore’s model. First note that the free-surface viscosity b can be obtained in closed form since convection effects play no role: b = i (ı + (1 − ı)e−AT ) (50) ∞ Second, the substrate viscosity s varies in both T and X. The governing equation (32) suggests the following representation s (T, X) = b (T )(1 + (T )| sin(X)|)
∂ −1 = −Aı − 3B b 2 b ∂T
(52)
b
The rate equation for amplitude is found by substituting H = 1 + cos X and s = b (1 + | sin (X) | ) in the evolution equation (28). To leading order we find
1−
2
0
20
40
60
80
100
T Fig. 9. Moore’s model: Linear decay rate −(1/ )d /dT for varying viscosity ratio ı = 1, 10, 100 in (21) (A = 0.04, B = 2, and i = 0 in (22)).
Eqs. (52) and (53) constitute a generalization of Orchard’s equation (49) for thixotropic liquids. The linear decay rate of − ˙ / is shown in Fig. 8 for Moore’s model, for ı = 100 (21) and A = 0.04 and various values of the parameter B with a broken-down initial structure (0, X, Z) = i = 0 in (22). Increasing the value of B, with fixed parameters A and ı, enhances shear effects, and thus increases leveling rates. In Fig. 9 the decay rate is plotted again for Moore’s model, for A = 0.04 and B = 2 in (22) and varying viscosity ratio, ı = 1, 10, 100 in (21). The leveling rate remains constant for a Newtonian liquid (ı = 1), whereas it decreases with time for a thixotropic liquid due to rebuilding effects and the weakening of shear effects in time. Thus, increasing the viscosity ratio ı worsens leveling. 6. Conclusions
(51)
with | |1. Introducing (51) in (32) we obtain, after retaining only the linear terms,
1 d 1 = dt b
0.1
0.01
Orchard [3] developed a linear theory able to predict the amplitude evolution in time for Newtonian liquids. If the film profile is initially perturbed in a sinusoidal way, its profile takes the form
−
-1/ε dε/dT
if 0 ≤ Z ≤ Hy =
100
Fig. 8. Moore’s model: Linear decay rate −(1/ )d /dT with ı = 100 in (21) and A = 0.04, B = 0.2, 2, 20, and i = 0 in (22).
and the time-dependent nonlinear inverse of viscosity approximation is 1 1 = + (T, X, Z) y (T, X)
80
T
(53)
We have presented a hydrodynamic model for the leveling of thixotropic liquids based on the lubrication approximation. Despite the considerable simplifications which ensue from the lubrication approximation, the microstructure has to be calculated at every location (x, z) of the liquid domain. This would still be quite computationally unwieldy for flow problems of more complexity than leveling flows. In order to make computations more tractable, we developed two approximation schemes. These schemes complement the lubrication approximation depending on the rheological model used. Our first simplification is to approximate the inverse of viscosity as a linear function across the thickness of the film,
S. Livescu et al. / J. Non-Newtonian Fluid Mech. 166 (2011) 395–403
1/ = 1/b + (1/s − 1/b )(1 − z/h), using the boundary values of the viscosity at the substrate, s , and at the free surface, b . This results in considerable simplification since the problem can be reduced to solving for the unknowns (h, s , b ) which are not functions of z. Our second procedure is to use the more general form 1/ = 1/b + (1/s − 1/b )(1 − z/h)r with a “dynamic” exponent r ≡ r(t). This allows a wider range of applicability and more reliable numerical results. We have also developed a linear analysis which yields a generalization of the Orchard’s law of leveling for Newtonian liquids to take into account the effect of thixotropy. According to this analysis, the linear decay rate of − ˙ / depends on the viscosity values at the substrate and the free surface. More efforts are worthwhile. One possible extension of this work would involve rheological models which account for elastic or yield stress effects. This will be the subject of a future publication. Also, efforts for validating the current model against experimental data will be devoted in the future. Acknowledgment This work was supported by The ICI Strategic Research Fund. Appendix A. Coefficients aL and ac for Moore’s Model In order to simplify the expression of the coefficient ac we note that
H
∂ ∂ U +W ∂X ∂Z
0
dZ =
∂ 1 ∂H + b ∂T ∂X
1
H
0
dZ =
W + Z=H
0
H
∂ ∂X
U
U dZ
where we have used the continuity equation ∂ U/∂X = − ∂ W/∂ Z for the first term in the integral. This last expression gives a zero contribution once integrated from X = 0 to X = 2 given the no-flow conditions at the boundaries of the domain. Thus, we obtain the expression
ac = 0
2
1 ∂H dX b ∂T
(A.1)
which is independent on the rheological model used. Note that we can write a2 + ac =
d dt
2
0
(H/ b )dX.
403
The coefficient aL is model-dependent: for Moore’s model the operator L1 is given by L1
1
=A
1 ı − 2
1 ∂3 H 1 1− Z − 2 3 3 H ∂X
− 3BH
H
1
(A.2)
The expression for the integral 0 L1 dZ is then readily obtained after replacing 1/ by 1/ b + (1/ s − 1/ b )(1 − Z/H)r in (A.2). A similar procedure could be followed to obtain aL for BQP’s model. References [1] J.S. Dodge, Quantitative measures of leveling, J. Paint Tech. 44 (1972) 72–78. [2] A. Quach, Polymer Coatings. Physics and mechanics of leveling, Ind. Eng. Chem. Prod. Res. Develop. 12 (1973) 110–116. [3] S.E. Orchard, Surface levelling in viscous liquids and gels, Appl. Sci. Res. A 11 (1962) 451–464. [4] R.R. Iyer, D.W. Bousfield, The leveling of coating defects with shear thinning rheology, Chem. Eng. Sci. 51 (1996) 4611–4617. [5] R. Keunings, D.W. Bousfield, Analysis of surface tension driven leveling in viscoelastic films, J. Non-Newtonian Fluid Mech. 22 (1987) 219–233. [6] H.S. Kheshgi, L.E. Scriven, The evolution of disturbances in horizontal films, Chem. Eng. Sci. 43 (1988) 793–801. [7] H.A. Barnes, Thixotropy – a review, J. Non-Newtonian Fluid Mech. 70 (1997) 1–33. [8] J. Mewis, Thixotropy – a general review, J. Non-Newtonian Fluid Mech. 6 (1979) 1–20. [9] J. Mewis, N.J. Wagner, Thixotropy, Adv. Colloid Interface Sci. 147–148 (2009) 214–227. [10] J. Billingham, J.W.J. Ferguson, Laminar, unidirectional flow of a thixotropic fluid in a circular pipe, J. Non-Newtonian Fluid Mech. 47 (1993) 21–55. [11] D.C.-H. Cheng, F. Evans, Phenomenological characterization of the rheological behaviour of inelastic reversible thixotropic and antithixotropic fluids, Br. J. Appl. Phys. 16 (1965) 1599–1617. [12] H.T. Huynh, N. Roussel, P. Coussot, Aging and free surface flow of a thixotropic fluid, Phys. Fluids 17 (2005) 033101. [13] A. Mujumdar, A.N. Beris, A.B. Metzner, Transient phenomena in thixotropic systems, J. Non-Newtonian Fluid Mech. 102 (2002) 157–178. [14] O. Cohu, A. Magnin, The leveling of thixotropic coatings, Prog. Org. Coatings 28 (1996) 89–96. [15] F. Moore, The rheology of ceramic slips and bodies, Trans. Br. Ceramics Soc. 58 (1959) 470–494. [16] C. Baravanian, D. Quemada, A. Parker, Modelling thixotropy using a novel structural kinetics approach: basis and application to a solution of iota carrageenan, J. Texture Studies 27 (1996) 371–390. [17] M.H. Eres, L.W. Schwartz, R.V. Roy, Fingering phenomena for driven coating films, Phys. Fluids 12 (2000) 1278–1295. [18] D.E. Weidner, L.W. Schwartz, Contact-line motion of shear-thinning liquids, Phys. Fluids 6 (1994) 3535–3539.
Contents lists available at ScienceDirect
Journal of Non-Newtonian Fluid Mechanics journal homepage: www.elsevier.com/locate/jnnfm
Leveling of thixotropic liquids S. Livescu ∗ , R.V. Roy, L.W. Schwartz Mechanical Engineering Department, University of Delaware, Newark, DE 19716, USA
a r t i c l e
i n f o
Article history: Received 11 November 2010 Received in revised form 20 January 2011 Accepted 25 January 2011 Available online 23 February 2011 Keywords: Hydrodynamics Leveling Lubrication Thixotropy Thin films
a b s t r a c t In this paper the effect of thixotropy in the hydrodynamic behavior of thin films is studied. The simple problem of leveling on a horizontal substrate is considered. The rheological properties of the material are assumed to evolve over time due purely to changes in its internal structure. These changes are modeled in terms of a single structural variable. Neither elastic nor yielding effects are taken into account. More specifically, two distinct rheological models are considered: the simple model proposed by Moore and the more complex model proposed by Baravanian et al. These models exhibit a large range of variation for the liquid viscosity across the film thickness. After deriving the hydrodynamic equations governing leveling flows with the standard assumptions required by the lubrication approximation and running time-dependent numerical simulations, the nonlinear leveling history of the liquid can be predicted as a function of the initial microstructural state, rheological parameters, and initial disturbance of the liquid free surface. The main effort of this work is devoted to devising approximation schemes which lead to significant simplifications of the governing equations and their numerical computations. By approximating the inverse of viscosity as a monotonic function between its substrate and free-surface values, excellent agreement is found for the film amplitude, irrespective of the values of the rheological parameters of both models. Finally, a linear analysis yields a generalization of the Orchard’s law of leveling for Newtonian liquids to take into account the effect of thixotropy. © 2011 Elsevier B.V. All rights reserved.
1. Introduction Many industrial processes require coating of solid surfaces with thin liquid films, for both decorative and protective reasons. When a liquid film is applied to a flat surface, the initial coating is generally uneven. It is usually desired that the final coating is as uniform as possible. Many of these liquids such as paints, inks, and coatings have thixotropic rheology as their viscosities change with shear. Considering the size of the industries and processes for which the leveling flow of thixotropic liquids is important, it is perhaps surprising that numerical modeling of this type of flow is still in an incipient phase. The hydrodynamic leveling of thin liquid films has been studied extensively since 1920, most particularly for Newtonian liquids (see Dodge [1] and Quach [2], and references therein). Orchard [3] was the first to treat this problem theoretically. He considered the flow under viscous, capillary and body forces. In the limit of small Reynolds number, Orchard used a Fourier series to solve the twodimensional Navier-Stokes equations governing a liquid film with
∗ Corresponding author. Present address: ExxonMobil Upstream Research Company, 3120 Buffalo Speedway, Houston, TX 77056, USA. Tel.: +1 713 431 4578; fax: +1 713 431 6392. E-mail address: [email protected] (S. Livescu). 0377-0257/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.jnnfm.2011.01.010
a sinusoidal free-surface. His linear model predicts the rate of leveling, the time evolution of the free-surface amplitude, and the liquid stress distribution. Orchard obtained good agreement between his linear theory and the numerical solution of the nonlinear evolution equation governing the film thickness. Thin film hydrodynamics studies taking into account nonNewtonian effects such as shear-thinning [4] or visco-plasticity [5] have also received attention. Iyer and Bousfield [4] developed a numerical technique to include a generalized Newtonian rheology in thin-film leveling equations. Their shear-thinning rheology follows Carreau’s model. Their numerical solution of the onedimensional evolution equation for the film thickness is within 2% of the results obtained by numerical simulation of the full twodimensional problem using finite-element techniques. They also performed experiments whose results were in good agreement with their numerical simulations. Keunings and Bousfield [5] studied the effect of viscoelasticity on surface tension driven leveling of thin films. They generalized Orchard’s linear analysis by incorporating viscoelastic rheology with a generalized Maxwell model. Their results show that elasticity retards the leveling process. They also found that, at early times, the nonlinear evolution of the free surface perturbation is similar to that obtained in the Newtonian case by Kheshgi and Scriven [6]. Thixotropy is usually defined as the reversible variation of viscosity with time, whether elastic effects are present or not (see
396
S. Livescu et al. / J. Non-Newtonian Fluid Mech. 166 (2011) 395–403
Barnes [7], Mewis [8] and Mewis and Wagner [9]). According to Barnes [7], thixotropic liquids have gel-like properties which disappear when sheared, but reappear when put to rest. The rheological properties of such liquids depend on the time needed by their microstructure to evolve from one state to another. Structural changes in a flowing thixotropic liquid are due to two competing effects: break-down due to flow stresses and build-up due to collisions of the particles which make up the microstructure. Thixotropy has important consequences on the flow of liquid films. For instance, during brushing of a thixotropic paint, the viscosity should be sufficiently low in order to obtain good spreading. The paint must subsequently possess an optimal rate of structural build-up for good leveling: if the build-up rate is too slow, the viscosity remains low and the paint will have good leveling but will tend to quickly sag; if the build-up rate is too high, the viscosity will become large, sagging will be minimal but leveling will be poor. Another example was presented by Dodge [1] for two latex (highly nonlinear viscosity variation in time) and alkyd (linear viscosity variation in time) exterior paints. The viscosity behavior of any of these two paints cannot accurately be predicted by assuming that the paint is Newtonian. Dodge noted that the degree of viscosity break-down during paint application and the speed the viscosity increase following application will determine the extent of leveling, for a given film thickness and wavelength of brushmarks. There are very few published hydrodynamic studies of thixotropic liquids (see Billingham and Ferguson [10], Cheng and Evans [11], Huynh et al. [12] and Mujumdar et al. [13]) and we know of only one reference which considers the leveling of thixotropic liquids: Cohu and Magnin [14] developed a computational method by stepwise calculations based on Orchard’s analysis in order to incorporate thixotropic rheology with a five-parameter model. Their stepwise method has the limitation that the shear stress is homogeneous within the liquid. However, their numerical calculations of the leveling rate are in good agreement with their experiments. Our goal is to incorporate thixotropic rheology into the hydrodynamic equations governing the flow of a thin film on a horizontal substrate. For this, we consider two different rheological models: Moore’s model [15] and the model of Baravanian et al. (referred to as BQP) [16]. Moore’s model is a natural choice for its simplicity. It assumes that the microstructure of a thixotropic liquid is composed of a large number of “links”. The rheological behavior of such a liquid depends on the number of links formed at a given moment in time. Thus, the microstructure is described in terms of a nondimensional scalar parameter, denoted by . This parameter will reflect the structural build-up and break-down states. The state of build-up is defined as the number of links which are formed at any instant divided by the total number of links and its corresponding state of break-down is defined as the number of links destroyed at that given time divided by the total number of links. It follows that will have the value 0 if the structure is completely brokendown (all links are destroyed) and the value 1 if the structure is completely recovered (all links are restored). Generally, the behavior of a thixotropic material is described by constitutive equations consisting of a rate equation and an equation of state, D = f (, ), ˙ Dt
= (, ˙ ), ˙
(1)
where , ˙ and are the shear stress, shear rate and viscosity, respectively. Here, the operator D/Dt represents the standard material derivative. We assume that elastic effects are not present, or that such effects appear much faster than thixotropic effects. Hence we confine ourselves to modeling purely viscous phenomena and we take the viewpoint that the rheological properties of the mate-
rial evolve over time due purely to changes in its internal structure adequately modeled by a single structural variable . In Moore’s model, the time rate of change of is expressed as the difference of two terms: (i) a break-down rate due to shearing expressed as the product of the structural parameter by some power of the shear rate, and (ii) a build-up rate expressed as some power of (1 − ). For simplicity we take both exponents equal to unity, as in the original equation proposed by Moore, D = a(1 − ) − b|| ˙ Dt
(2)
Thus, it is seen that the shear-induced break-down term results in a decrease of , whereas the build-up term leads to an increase of . The constitutive equations are completely specified by defining the viscosity as a linear function of : in the limit of high shear rate, → 0, the viscosity reaches the value = ∞ ; conversely, the viscosity of a fully recovered structure, → 1, takes the value = 0 , = ∞ + (0 − ∞ )
(3)
where ∞ and 0 are the limiting viscosities at very high and very low shear rates, respectively. Hence Moore’s model is a fourparameter model (∞ , 0 , a and b). BQP’s model uses the shear stress in the rate equation 1 D = Dt ta
1−
e
(4)
The equilibrium value e and the relaxation time ta are determined from the following empirical formulas:
e () = 1 +
p −1 c
,
ta () = ta,∞ −
ta,∞ − ta,0 1+
a,0
q
(5)
where c , p, q, ta,∞ , ta,0 , a,0 are rheological parameters. The viscosity varies between the two limiting values, ∞ and 0 , and it is again only a function of the structural parameter , =
∞ (1 − K)
2
,
K =1−
1/2 ∞
(6)
0
This model requires eight material parameters ∞ , 0 , c , a,0 , p, q, ta,∞ , ta,0 . In order to compare these two models, consider a thin film of thickness h0 and assume that the shear stress is distributed linearly across the film: = 0 (1 − z/h0 ), 0 ≤ z ≤ h0 . In fact such a linear profile can be shown to be a direct consequence of the lubrication approximation of thin films (see the next section for the derivation of the linear shear stress profile in the lubrication approximation limit). We can determine the thickness profile of the equilibrium viscosity e for both models. For Moore’s model we find e = 1 + ıe , ∞
e =
1 1 + AM (1 − z/h0 )/[1 + ıe ]
with ı = 0 /∞ 1 and AM = b 0 /a∞ . For BQP’s model we find 1 e = , 2 ∞ (1 − (1 − ı−1/2 )e )
e =
1 1 + [ABQP (1 − z/h0 )]
p
with ı = 0 /∞ and ABQP = 0 / c . Fig. 1 shows the variation of e /∞ across the thickness for the two models and a range of parameters. Striking differences can be noted. Moore’s model leads to a quasilinear profile between the substrate value and the free-surface value. BQP’s model leads to a highly nonlinear profile: for high values of ABQP (or high values of exponent p) the film displays a thin skin of high viscosity on top of a thick layer of low quasi-uniform viscosity. Thus, for a large range of parameters (both AM and ABQP vary by two orders of magnitude in Fig. 1), the two models have very different viscosity behavior. This conclusion emphasizes our
S. Livescu et al. / J. Non-Newtonian Fluid Mech. 166 (2011) 395–403
film. l is the wavelength of a sinusoidal perturbation and h0 is the mean thickness of the liquid film. Assuming h0 l, the standard lubrication approximation can be adopted to simplify the hydrodynamic equations. It has been used very effectively to model thin layer flows of Newtonian and generalized Newtonian liquids (see for example [17] or [18]). The key assumptions made in the lubrication approximation are: (i) the flow is slow and inertial forces are neglected, (ii) the liquid film is thin, and (iii) free-surface slopes are small. Considerable simplification of the governing equations follows from these assumptions. For a liquid film flowing on a horizontal substrate, the mass conservation equation is
100 90 70
ηe/η∞
A M=1
M BQP
80
A M=10
A M=100
60 ABQP = 0.05
50
ABQP = 0.5
40
ABQP = 5
30
397
20 10 0 0
0.2
0.4
0.6
0.8
1
z/h0 Fig. 1. Variation of the viscosity at equilibrium, e , for the two rheological models: Moore’s model (“M”) and BQP’s model with exponent p = 1. A film of constant thickness h0 is assumed to be under a linear shear–stress profile = 0 (1 − z/h0 ).
∂w ∂u + =0 ∂x ∂z
(7)
and the simplified momentum equations in the x and z directions are ∂xz ∂p ∂p = , = −g , ∂x ∂z ∂z
(8)
where p is the pressure and xz is the shear stress and and g are the liquid density and gravitational acceleration, respectively. Furthermore at the liquid free surface (z = h), a normal stress balance is written as p = − ∂ 2 h/∂ x2 , where is the surface tension and the curvature of the free surface has been linearized. The shear stress and shear rate within the film are, respectively (0 ≤ z ≤ h)
xz =
∂3 h ∂h − g ∂x ∂x3
(h − z) ,
˙ xz =
xz ∂u = ∂z
(9)
The horizontal and vertical components of the liquid velocity are obtained by integration of (9) Fig. 2. Sketch of a thin film of mean thickness h0 perturbed by a sinusoidal ripple of wavelength l.
motivation for using these two rheological models in the present study. While Moore’s model is chosen for its simplicity (the viscosity varies quasi-linearly across the film thickness), BQP’s model is an arbitrary choice displaying a complex viscosity behavior (the viscosity has a highly-nonlinear variation across the film thickness). The lubrication approximation takes into account, among other things, the thinness of the liquid film: after averaging across the thickness of the film, one can derive an evolution equation for the film thickness h(t, x) (in a two-dimensional domain). The difficulty in solving a flow problem for a thixotropic liquid stems from the fact that the flow behavior is coupled to the history of the microstructure which must be determined in both time and space. Calculation of microstructure at every location (x, z) of the liquid domain is prohibitive and thus simplifying schemes are pivotal to efficient hydrodynamic calculations that take advantage of the lubrication approximation. In order to assess the reliability of such schemes, we must nevertheless solve the full set of modeling equations numerically. This paper proceeds as follows. In Section 2, we implement the two thixotropic models into the leveling hydrodynamics. We arrive at a set of evolution equations which we compute numerically in Section 3. In Section 4, we develop two simplifying approximations which present different degrees of accuracy for the two rheological models. Finally in Section 5, we generalize the so-called Orchard law for leveling of Newtonian liquids to account for the thixotropy effect. 2. Hydrodynamic model With the geometry shown in Fig. 2, we choose the x-axis along the substrate and the z-axis pointing toward the free surface of the
u=
∂3 h ∂h 3 − g ∂x ∂x
0
z
h− d ,
z
∂u d ∂x
w=− 0
(10)
With the usual no-slip condition at the substrate, u = w = 0 at z = 0, and no-stress free surface condition, xz = 0 at z = h, integrating (7) in the z direction and using the Leibnitz rule and the kinematic condition at the free surface which states that the liquid does not cross the interface, Dh/Dt = w(h), the mass conservation produces ∂ h/∂t = − ∂ Q/∂ x, where
Q =
h
udz = 0
0
h
∂u h (h − z)dz − [u(h − z)]0 = ∂z
0
h
∂u (h − z)dz ∂z (11)
is the flow rate per unit film width. Finally these equations yield the following evolution equation for the film thickness h(x, t) ∂ ∂h =− ∂t ∂x
∂3 h ∂h 3 − g ∂x ∂x
0
h
2
(h − z) dz
(12)
Eqs. (9), (10) and (12) along with the rheological equations (2) and (3) for Moore’s model (Eqs. (4)–(6) for BQP’s model) fully pose the leveling problem once initial conditions are specified and symmetry conditions are imposed at the endpoints x = 0 and x = l. We assume that at t = 0 the structure is uniform in space and characterized by an initial value i and corresponding viscosity i (i ). Also note that the material derivative D/Dt takes the form ∂/∂t + u∂/∂x + w∂/∂z. For a Newtonian liquid at constant viscosity i , one recovers the expression of the hydrodynamic conductivity
h
2
[(h − z) /]dz = h3 /(3i ). For a thixotropic liquid, determina0 tion of this conductivity requires the evaluation of the integral across the film thickness. Thus, unless further approximations are made, the advantage gained by the lubrication approximation is significantly reduced.
398
S. Livescu et al. / J. Non-Newtonian Fluid Mech. 166 (2011) 395–403
All these equations have been developed for liquids without yield stress. If we assume that the liquid behaves approximately as Bingham plastic, the shear strain relation is nonlinear ∂u = 0, ∂z
if |xz | ≤ y = xz − sign (˙ xz ) y ,
if |xz | ≥ y
(13)
where y is the yield stress. This is the simplest idealization of the shear-thinning pseudo-plastic liquid. Other expressions are available in literature (see for example Quach [2], Barnes [7], Mewis [8], Mewis and Wagner [9], Billingham and Ferguson [10], Cheng and Evans [11], Cohu and Magnin [14], and references therein). According to (9), the shear stress varies linearly across the film thickness from xz = 0, at z = h, to xz = h( ∂ h3 /∂ 3 x − g ∂ h/∂ x), at z = 0. For a fixed value of x, there are two possible cases, depending on the maximum value of the shear stress, |xz |z=0 : (i) |xz |z=0 ≤ y ; in this case, u = w = 0 everywhere in the liquid domain due to the noslip condition at the substrate, i.e., the liquid does not flow; and (ii) |xz |z=0 > y ; in this last case, there is a so-called yield surface defined by hy = h −
y ∂h3 ∂3 x − g ∂h ∂x
(14)
assuming that the denominator is non-zero. Below the yield surface, 0 ≤ z ≤ hy , there is shearing. Above the yield surface, hy ≤ z ≤ h, there is a layer of plug flow moving with u = u(hy ) (independent of z). Thus, the flow rate per unit film width is
hy
Q =
udz + 0
h
hy
udz = 0
hy
∂u hy − z dz + h − hy u(hy ) ∂z
With these notations, the dimensionless evolution equation for the film thickness takes the form
∂ ∂h =− ∂t ∂x
∂3 h ∂h − g ∂x ∂x3
hy
−sign (˙ xz ) y 0
hy
0
H
(H − Z)2 dZ ,
F(H) 0
Z
U = 3F(H) 0
W =− 0
Z
∂U d ∂X
∂3 h ∂h − g ∂x ∂x3
hy
0
h−z dz − sign (˙ xz ) y
hy
0
2 x , l
=
∞
Z=
z , h0
H=
h , h0
T=
t , t0
U=
u , u0
W=
(17)
w , w0
where we use the viscosity ∞ of the broken-down structure for the reference viscosity. The relaxation time associated with the leveling of a Newtonian fluid of viscosity ∞ (also known as Orchard’s time scale [3]) was chosen as time scale t0 =
3∞ l4
(18)
(2 )4 h30
The velocity scales u0 and w0 are then found to be u0 =
l = 2 t0 3∞
2 h 3 0
l
,
w0 =
2 h0 u0 l
(20)
and
H−Z ∂ ∂ ∂ +U +W = A(1 − ) − 3B F(H) ∂T ∂X ∂Z
(22)
with ı=
0 , ∞
A = at0 ,
B=
bl 2 h0
Note that Moore’s model becomes Newtonian in two situations: ı = 1 or A = B = 0. For BQP’s model, the dimensionless rheological equations can be stated as follows: 1 (1 − K)
(23)
2
e =
1−
1 1 + [3˛|F(H)|(H − Z)]p
,
e
(24)
Ta = Ta,1 +
Ta,2 − Ta,1
q
1 + 3ˇ|F(H)|(H − Z) (25)
with
1 dz
The procedure for solving numerically Eqs. (12) and (16) is similar once hy is calculated. This study will be the subject of a future publication. We return now to our model without yield stress (12). In order to obtain dimensionless equations, we use the following nondimensional variables: X=
(19)
(21)
where the horizontal velocity component at the yield surface is u(hy ) =
∂3 H ∂H − Bo ∂X ∂X 3
= 1 + (ı − 1)
(16)
H− d ,
∂ ∂ 1 ∂ +U +W = Ta ∂T ∂X ∂Z
hy − z dz + h − hy u(hy )
F(H) =
They are needed to compute the convective terms in the rate equation for the structural parameter . For Moore’s model, the dimensionless rheological equations requires three dimensionless parameters (ı, A, B) and the initial structural parameter i . The equations are
=
(h − z)(hy − x) dz
where Bo = gl2 / is the Bond number. The dimensionless components of the liquid velocity are
(15) and the evolution equation of the film thickness is
∂ ∂H = −3 ∂T ∂X
ı=
0 , ∞
Ta,1 =
K =1−
ta,∞ , t0
1 ı1/2
Ta,2 =
,
˛=
∞ u0 , h0 c
ˇ=
∞ u0 , h0 a,0
ta,0 t0
This model requires seven dimensionless parameters ı, ˛, ˇ, p, q, Ta,1 , Ta,2 and the initial value i . It becomes Newtonian for ı = 1 or Ta,1 = Ta,2 = 0. 3. Time-dependent numerical simulations The equations derived in the previous section show that timedependent numerical simulation of the leveling of thin films requires the calculation of the microstructure at every location (x, z) of the liquid domain. It would be in the spirit of the lubrication approximation to develop a simplifying scheme which would involve only the x space variable by averaging the rheological properties across the film thickness. In order to develop such a scheme and to assess its reliability, we must nevertheless solve the full set of equations numerically with two simplifying assumptions. First we will assume that Bo 1, i.e., the wavelength l of the free-surface perturbations is much shorter than the capillary length (/g)1/2 , i.e., capillarity dominates gravity in the leveling process. This simplification is supported by observations that little difference exists between paints applied on the floor and on the ceiling [2]. We will
S. Livescu et al. / J. Non-Newtonian Fluid Mech. 166 (2011) 395–403
0.1
0.05
num approx
0.08
399
num approx
0.048
0.06
0.046 0.044
0.04
0.042 0.04
1/Λ
B=2
ε
B=0.2
0.02
0.038 0.036 0.034
0.01
B=20
B=0.2
B=2
B=20
0.032 0.03
0.005
0.028 0
20
40
60
80
100
0
0.2
0.4
0.6
T Fig. 3. Moore’s model: Evolution of the dimensionless amplitude of the free surface, , by time-dependent numerical simulation of the full equations (neglecting convection) for the rheological parameters ı = 100 in (21) and A = 0.04, B = 0.2, 2, 20, and initial structure i = 0 in (22). The dashed curves (“approx”) correspond to approximating the inverse of viscosity as a linear function of Z (see Eq. (27)).
=3 0
2
(H − Z) dZ
(26)
is evaluated at the old time step, T(k) . Thus, at each time step, solving (19) requires the solution of a pentadiagonal system of linear equations with (N + 1) unknowns for the liquid thickness at the mid(k+1) points Hi+1/2 , i = 0 . . . N. For this we use the well-known Thomas algorithm for banded systems. The computational effort is of order O(N) rather than O(N3 ) for the solution of a general system of equations solved by Gaussian elimination. At each time step the rate equations (22) and (24) are treated explicitly and the new values of (k+1) are converted into values of (dimenthe structural parameter i,j (k+1)
1.2
0.15
λi=0 λi=0.2 λi=0.5 λi=1
sionless) viscosity i,j via (21) and (23). At each time step the accuracy of the scheme is verified by checking that the total mass of liquid is conserved. An example of the film amplitude (T) evolution is shown in Fig. 3. The rheology follows Moore’s model with the rheological parameters ı = 100 in (21) and A = 0.04 in (22). The initial profile is H(0, X) = 1 + 0 cos X with initial amplitude 0 = (0) = 0.1. Fig. 3 is a semilogarithmic plot showing the effect of increasing parameter B; as expected, as B is increased, the shear effects on the microstructure are stronger, thus leading to smaller viscosity at any given time and improved leveling. Note that in all cases, the leveling is initially strongly non-Newtonian, and eventually asymptotes to an exponential decay, whose time constant corresponds to the viscosity 0 of the fully recovered structure.
0.06 0.04
ε
also neglect the effect of convection in the evolution of . Simple inspection of (2) shows that convection is negligible for Moore’s model if either u0 /l a or u0 /l bu0 /h is assumed, i.e., if the nondimensional groups A and B satisfy A 1 or B 1. Unless specified otherwise, the initial microstructure is assumed uniformly brokendown: i = (T = 0, X, Z) = 0. Eq. (19) and the rheological equations corresponding to each model can be solved using finite-difference methods. Eq. (19) has the character of a higher-order diffusion equation. To obtain the finite difference form of (19) we use a forward difference in time and centered differences in space. If the liquid domain is discretized as (X, Z) → (Xi , Zj ) = (i x, j z) all derivatives are evaluated at the midpoints (Xi+1/2 , Zj ), where Xi+1/2 is the midpoint between Xi and Xi+1 ; i.e., Xi+1/2 = (i + 1/2) x. ∂ 3 H/∂X3 is evaluated at the new time step, T(k+1) , while the term H
1
Fig. 4. Moore’s model: Variation of the inverse of viscosity 1/ across the film thickness at X = /3 and time T = 10 (same rheological parameters as in Fig. 3). The inverse of viscosity is maximum at the substrate (Z = 0) where the shear is maximum.
0.1
0.8
Z
0.02
0.01
0.004 0
50
100
150
200
T Fig. 5. Moore’s model: Effect of initial structure on the evolution of the film amplitude. The rheological parameters are fixed to ı = 100 (21) and A = 0.04, B = 2, and i = 0, 0.2, 0.5, 1 in (22). The more broken-down the initial structure, the better the leveling.
In Fig. 4 the inverse of viscosity variation across the film thickness is shown for the same values of the rheological parameters. The curves are shown at time T = 10 at a fixed horizontal position, X = /3 (chosen to avoid any obvious symmetry). A few things can be noted. First, the inverse of viscosity is nearly a linear function of Z, a property which we can use to devise a simplifying scheme. Second, all curves reach the same value at the free surface; indeed the stress is zero there, and hence the viscosity always recovers monotonically from ∞ (at t = 0) to 0 (as t → ∞) irrespective of the value of B. Conversely, shear effects are maximum at the substrate (Z = 0), the viscosity is smallest there, and the difference between the free-surface and substrate viscosities increases as B is increased. Understandably, the initial state of the microstructure has a significant effect on leveling. In Fig. 5 the evolution of film amplitude is shown for i = 0, 0.2, 0.5, 1 fixing the rheological parameters to ı = 100, A = 0.04, and B = 2 in (21) and (22), respectively. As i is increased, the initial viscosity is increased, therefore the free surface takes more time to level. Corresponding results for BQP’s model are shown in Figs. 6 and 7 for the parameters ı = 100, ˛ = 1.5, ˇ = 3, p = 0.5, q = 4, Ta,1 = 2.5, Ta,2 = 0.1, and initial structure i = 0 in (24) and (25). The solid curves are obtained by numerical simulation of the full equations. The additional curves correspond to the simplifying schemes introduced in the next section. A notable difference between Moore’s
400
S. Livescu et al. / J. Non-Newtonian Fluid Mech. 166 (2011) 395–403
4.1. Linear inverse of viscosity approximation
0.1 num r=0.5 r=1 r(t)
0.06 0.04
The inverse of viscosity variation across the film thickness in the case of Moore’s model (see Fig. 4) suggests the following scheme: approximate the inverse of viscosity as a linear function of Z in terms of its free-surface value 1/ b (Z = H) and its substrate value 1/ s (Z = 0), i.e.,
ε
0.02 0.01
1 1 Z = + H (T, X, Z) s (T, X)
0.006 0.004
1 1 − b (T, X) s (T, X)
(27)
The evolution equation for H now becomes
0.002
∂H 1 ∂ =− 4 ∂X ∂T
H3
0.0008 0
20
40
60
80
100
T Fig. 6. BQP’s model: Evolution of the film’s amplitude for the rheological values ı = 100, ˛ = 1.5, ˇ = 3, p = 0.5, q = 4, Ta,1 = 2.5, Ta,2 = 0.1, and initial structure i = 0 in (24). “Full” numerical simulations, denoted by “num”, are compared with the simplifying approximation with fixed value of the exponent r = 0.5, 1 and with time-varying exponent denoted by “r(t)”. Best agreement is obtained by fixing r = 0.5 or by using a time-varying exponent.
1 b
+
3 s
(28)
The inverse of viscosity evaluated at the substrate and at the free surface can be related to corresponding values of according to s = 1 + (ı − 1)s ,
b = 1 + (ı − 1)b
(29)
where s = s (T, X, Z = 0) and b = b (T, X, Z = H). To find the evolution equations for s and b from that of (22), let us define (T, X) = (T, X, Z = (T, X))
(30)
Upon using the identity (chain rule of differentiation)
0.11
num r=0.5 r=1 r(t)
0.1 0.09
1/Λ
∂3 H ∂X 3
∂ ∂T
=
∂ ∂ ∂ + , ∂T ∂Z ∂T
0.08
(22) becomes
0.07
∂ ∂T
0.06
=
−U
∂ + ∂X
∂ −W ∂T
∂3 H ∂ + A(1 − ) − 3B 3 (H − Z) ∂Z ∂X
(31) Z=
Thus, on the substrate, = 0, the no-slip condition (U = W = 0) implies that the convection terms drop out. This gives the following rate equation for s
0.05 0.04
∂3 H ∂s s = A(1 − s ) − 3BH 3 ∂T 1 + (ı − 1)s ∂X
0.03 0.02 0.01 0
0.2
0.4
0.6
0.8
1
1.2
(32)
At the free surface we have = H(T, X) and the kinematic condition
Z Fig. 7. BQP’s model: Inverse of viscosity across the film thickness at X = /3 and time T = 10. Comparison between the full numerical simulations (“num”), constant r-scheme (r = 0.5, 1) and time-varying r-scheme (“r(t)”). Parameters are identical to those of Fig. 6.
model and BQP’s model is seen in Fig. 7; the inverse of viscosity does not follow a linear variation across the thickness of the film. This property is also used in the next section in devising a simplifying scheme. 4. Simplifying approximations The numerical results obtained for the inverse of viscosity variation across the film thickness for the two rheological models presented in the previous section suggest two simplifying schemes: linear inverse of viscosity across the film thickness and nonlinear inverse of viscosity with time-dependent dependence across the film thickness. The main advantage of these approximations is that the term (26) is obtained in closed form as a product of H3 and a linear sum of values of the inverse of viscosity evaluated at the substrate and at the free surface. These schemes complement the lubrication approximation depending on the rheological model used. In this section, we will present these approximation schemes and compare their results with the “full” numerical simulation results.
∂H ∂H + Ub = Wb ∂T ∂X must hold. Here Ub = U(T, X, Z = H) and Wb = W(T, X, Z = H) are the components of the free-surface velocity. Upon using the identity ∂ (T, X, H(T, X)) = ∂X
∂ ∂H ∂ + ∂X ∂Z ∂X
Z=H
the following rate equation for b is found ∂b ∂ = −Ub b + A(1 − b ) ∂X ∂T
(33)
The method of characteristics can be used to show that if b (T = 0, X) = i is independent of X, then b (T, X) ≡ b (T) must remain independent of X, i.e., convection plays no role in the evolution of b . The characteristic curves corresponding to (33) are described by {(X(S), T(S)), 0 ≤ S < ∞ }, where dT = 1, dS
dX = Ub , dS
db = A(1 − b ) dS
(34)
Thus, if the structure is initially homogeneous, b varies in time in the same manner for all points of the free surface. Therefore, no convective effects are present in the equations describing the evolution of the substrate and free-surface structural parameters. Convective effects in the bulk of the liquid have not been neglected; yet their effect is made simpler because of the particular properties at the boundaries of the liquid layer.
S. Livescu et al. / J. Non-Newtonian Fluid Mech. 166 (2011) 395–403
Comparison between the results obtained by “full” numerical simulations and the solution of the coupled system of Eqs. (28), (32), and (33) are presented in Figs. 3 and 4 for Moore’s model. The amplitude agreement is very good, even for large values of the break-down parameter B. The agreement for the inverse of viscosity is also good. However, as B is increased, the inverse of viscosity becomes slightly nonlinear. This nonlinearity suggests that the linear inverse of viscosity approximation does not hold for a wide range of parameters or for other rheological models. In fact, the linear inverse of viscosity approximation did not obtain good agreement with the “full” numerical simulations obtained for BQP’s model, for which the inverse of viscosity has a strong nonlinear variation across the film thickness.
fluid domain is set equal to zero
2
H
dX 0
0
∂ −L ∂T
1
dZ = 0
(39)
where we have fixed the weight function to unity. The following equation can then be derived a0 (r)
dr = a1 (r)(r + 1) + [a2 (r) + ac (r) − aL (r)] (r + 1)2 dT
2
a0 =
H
1
0 2
−
s
a1 (r) = 0
1 b
dX
1
(41)
1 − s b
∂H ∂T
+H
∂ ∂T
1 s
r
1 1 − s (T, X) b (T, X)
(35)
Using this nonlinear Z-dependence of the inverse of viscosity leads to the evolution equation of the film thickness ∂H 1 ∂ =− r + 3 ∂X ∂T
H3
∂3 H ∂X 3
r b
+
3 s
r(T )
1 1 − s (T, X) b (T, X)
(37)
1
= −U
∂ ∂ −W + L1 ∂X ∂Z
1
≡L
1
(38)
where the spatial operator L1 is dependent on the rheological model used. The governing equation for r(t) is obtained as in the standard Galerkin method. If the inverse of viscosity (37) is substituted in (38) one obtains a non-zero residual. To make this residual as small as possible, the integral of the weighted residual over the
1
∂ ∂T
b
2
H
dX 0
aL (r) =
0
2
dX 0
This scheme would have the merit of using an exponent r which would vary in time in accordance with the flow and rheology evolution. In a more complex flow problem over a two-dimensional substrate, the exponent r can be made a function of time t and of the substrate variables (x, z). In a leveling flow, it is reasonable to assume that the exponent r is spatially uniform. In problems with significant spatial heterogeneity along the surface of the film, this assumption is unlikely to be acceptable. First note that the equation governing the non-dimensional inverse of viscosity 1/ can always be put in the form ∂ ∂T
0
The question then becomes how to choose exponent r. The numerical results presented in Figs. 6 and 7 obtained for the BPQ model show that a “static” value of r can be found that best fits the “full” numerical simulations (note that equations for the evolution of the boundary values b and s can be derived for the BPQ model in a similar manner). For instance r ≈ 0.5 is found to lead to the best agreement for the particular parameters chosen in Figs. 6 and 7. But this is not particularly appealing; although this approximation obtains more accurate results than the linear inverse of viscosity approximation, a scheme which would determine the optimal value of the parameter r without resorting to “full” numerical simulations must be instead developed. One solution is to treat the exponent r as time-varying and to find the ordinary differential equation that governs its evolution in a systematic way. Thus, the inverse of viscosity can be written as 1 Z 1 = + 1− H (T, X, Z) b (T, X)
H
ac (r) =
(36)
2
a2 (r) =
∂ − ∂T
1 b
dX
(42)
The linear inverse of viscosity approximation (27) could be improved by generalizing it to the nonlinear form 1 1 Z = + 1− H (T, X, Z) b (T, X)
(40)
where
4.2. Time-dependent nonlinear inverse of viscosity approximation
401
(43)
∂ ∂ +W U ∂X ∂Z
H
L1 0
dX
1
dZ
1
dZ
(44)
(45)
Note that a1 (r) = da0 /dT. Coefficients a0 , a1 and a2 are determined numerically once the functions H, s and b are known at any given time. Coefficients aL and ac can be simplified and their expressions for Moore’s model can be found in Appendix A. We apply this scheme by starting with Moore’s model. We choose the parameters ı = 100 in (21) and A = 0.04, B = 20, and initial structure i = 0 in (22) for which a slight discrepancy is seen when comparing the results of the linear inverse of viscosity approximation with the “full” simulations (see Fig. 3). Although the numerical results were quite satisfactory with the static value r = 1 (the static nonlinear approximation reduces to the linear approximation in this case), Moore’s model constitutes a good test case. The initial value r(0) = 1 is chosen, although, in principle any initial value can be taken for an homogeneous initial structure. The errors for the amplitude history along with the corresponding inverse of viscosity 1/ and time-evolution of r are reduced. The maximum error for the inverse of viscosity (at the substrate) is less than 0.2%. The evolution of r shows a rapid initial jump followed by a slow decrease to the asymptotic value r = 1. The chosen initial value for r(t) has little effect on its time evolution, apart for a very short time interval. Also the effect of convection (which we could “turn off” by setting ac = 0) is negligible, unsurprisingly for the chosen value of B. These results confirm the fact that the linear approximation first chosen to approximate 1/ was a reasonable “guess”. We now turn to the BQP model for which the linear approximation is inadequate. The equilibrium curves shown in Fig. 1 also attest to the need for an adjustable exponent r with the rheology. Figs. 6 and 7 show the results of the r(t)-scheme starting with the initial value r(0) = 0.5. Again, r exhibits a sharp increase and then a sharp decrease near the value r = 0.5 followed by a slow decrease towards r = 0.5. Varying the initial value of r gives the same behavior. The agreement with the “full” simulations is quite good. These calculations justify the conclusion that the static value r = 0.5 yields the best results. The r(t)-scheme of course has the advantage of being self-contained: it calculates the value of r which is best-suited for the flow/rheology of the problem at any given time. By turning off again the convection effect (by setting ac = 0), we found that its effect was always negligible for the cases considered here. The
S. Livescu et al. / J. Non-Newtonian Fluid Mech. 166 (2011) 395–403
success of the method is further demonstrated when the rheological parameter p in (25) is varied: the inverse of viscosity exhibits a wide range of variation across the film thickness which is well reproduced by the approximate scheme. Increasing the exponent p leads to higher substrate viscosity at any given time (the freesurface values remain identical), and thus to poorer leveling. Of the other parameters in (25), ˛ has the strongest effect on leveling; increasing ˛ leads to a decrease of the substrate viscosity, and thus to faster leveling. Note that in this study the two simplification schemes have been verified against the “full” numerical results for thixotropic liquids without yield stress. For liquids with yield stress, both approximations can be applied in a similar manner. For example, considering Bingham plastic liquids for which the film thickness evolves according to (16), the viscosity values at the substrate and at the yield surface are needed in the two approximation schemes. Thus, the linear inverse of viscosity approximation becomes 1 1 Z = + Hy (T, X, Z) s (T, X) if 0 ≤ Z ≤ Hy =
1 1 − y (T, X) s (T, X)
1 , if Hy ≤ Z ≤ H y (T, X)
,
1
-1/ε dε/dT
402
B=0.2 B=2 B=20
0.1
0.01 0
20
40
60
1.5 1
(46)
δ=1 δ=10 δ=100
1−
Z Hy
r(T ) ·
1 1 − s (T, X) y (T, X)
1 , if Hy ≤ Z ≤ H y (T, X)
,
(47)
5. Linear analysis
H(T, X) = 1 + (T ) cos(X)
( 1)
(48)
where the amplitude decays linearly according to d = − dT
(49)
Hence the film levels exponentially fast with the time relaxation given by Eq. (18) for a viscosity = ∞ Our goal here is to generalize Eq. (49) for thixotropic liquids obeying Moore’s model. First note that the free-surface viscosity b can be obtained in closed form since convection effects play no role: b = i (ı + (1 − ı)e−AT ) (50) ∞ Second, the substrate viscosity s varies in both T and X. The governing equation (32) suggests the following representation s (T, X) = b (T )(1 + (T )| sin(X)|)
∂ −1 = −Aı − 3B b 2 b ∂T
(52)
b
The rate equation for amplitude is found by substituting H = 1 + cos X and s = b (1 + | sin (X) | ) in the evolution equation (28). To leading order we find
1−
2
0
20
40
60
80
100
T Fig. 9. Moore’s model: Linear decay rate −(1/ )d /dT for varying viscosity ratio ı = 1, 10, 100 in (21) (A = 0.04, B = 2, and i = 0 in (22)).
Eqs. (52) and (53) constitute a generalization of Orchard’s equation (49) for thixotropic liquids. The linear decay rate of − ˙ / is shown in Fig. 8 for Moore’s model, for ı = 100 (21) and A = 0.04 and various values of the parameter B with a broken-down initial structure (0, X, Z) = i = 0 in (22). Increasing the value of B, with fixed parameters A and ı, enhances shear effects, and thus increases leveling rates. In Fig. 9 the decay rate is plotted again for Moore’s model, for A = 0.04 and B = 2 in (22) and varying viscosity ratio, ı = 1, 10, 100 in (21). The leveling rate remains constant for a Newtonian liquid (ı = 1), whereas it decreases with time for a thixotropic liquid due to rebuilding effects and the weakening of shear effects in time. Thus, increasing the viscosity ratio ı worsens leveling. 6. Conclusions
(51)
with | |1. Introducing (51) in (32) we obtain, after retaining only the linear terms,
1 d 1 = dt b
0.1
0.01
Orchard [3] developed a linear theory able to predict the amplitude evolution in time for Newtonian liquids. If the film profile is initially perturbed in a sinusoidal way, its profile takes the form
−
-1/ε dε/dT
if 0 ≤ Z ≤ Hy =
100
Fig. 8. Moore’s model: Linear decay rate −(1/ )d /dT with ı = 100 in (21) and A = 0.04, B = 0.2, 2, 20, and i = 0 in (22).
and the time-dependent nonlinear inverse of viscosity approximation is 1 1 = + (T, X, Z) y (T, X)
80
T
(53)
We have presented a hydrodynamic model for the leveling of thixotropic liquids based on the lubrication approximation. Despite the considerable simplifications which ensue from the lubrication approximation, the microstructure has to be calculated at every location (x, z) of the liquid domain. This would still be quite computationally unwieldy for flow problems of more complexity than leveling flows. In order to make computations more tractable, we developed two approximation schemes. These schemes complement the lubrication approximation depending on the rheological model used. Our first simplification is to approximate the inverse of viscosity as a linear function across the thickness of the film,
S. Livescu et al. / J. Non-Newtonian Fluid Mech. 166 (2011) 395–403
1/ = 1/b + (1/s − 1/b )(1 − z/h), using the boundary values of the viscosity at the substrate, s , and at the free surface, b . This results in considerable simplification since the problem can be reduced to solving for the unknowns (h, s , b ) which are not functions of z. Our second procedure is to use the more general form 1/ = 1/b + (1/s − 1/b )(1 − z/h)r with a “dynamic” exponent r ≡ r(t). This allows a wider range of applicability and more reliable numerical results. We have also developed a linear analysis which yields a generalization of the Orchard’s law of leveling for Newtonian liquids to take into account the effect of thixotropy. According to this analysis, the linear decay rate of − ˙ / depends on the viscosity values at the substrate and the free surface. More efforts are worthwhile. One possible extension of this work would involve rheological models which account for elastic or yield stress effects. This will be the subject of a future publication. Also, efforts for validating the current model against experimental data will be devoted in the future. Acknowledgment This work was supported by The ICI Strategic Research Fund. Appendix A. Coefficients aL and ac for Moore’s Model In order to simplify the expression of the coefficient ac we note that
H
∂ ∂ U +W ∂X ∂Z
0
dZ =
∂ 1 ∂H + b ∂T ∂X
1
H
0
dZ =
W + Z=H
0
H
∂ ∂X
U
U dZ
where we have used the continuity equation ∂ U/∂X = − ∂ W/∂ Z for the first term in the integral. This last expression gives a zero contribution once integrated from X = 0 to X = 2 given the no-flow conditions at the boundaries of the domain. Thus, we obtain the expression
ac = 0
2
1 ∂H dX b ∂T
(A.1)
which is independent on the rheological model used. Note that we can write a2 + ac =
d dt
2
0
(H/ b )dX.
403
The coefficient aL is model-dependent: for Moore’s model the operator L1 is given by L1
1
=A
1 ı − 2
1 ∂3 H 1 1− Z − 2 3 3 H ∂X
− 3BH
H
1
(A.2)
The expression for the integral 0 L1 dZ is then readily obtained after replacing 1/ by 1/ b + (1/ s − 1/ b )(1 − Z/H)r in (A.2). A similar procedure could be followed to obtain aL for BQP’s model. References [1] J.S. Dodge, Quantitative measures of leveling, J. Paint Tech. 44 (1972) 72–78. [2] A. Quach, Polymer Coatings. Physics and mechanics of leveling, Ind. Eng. Chem. Prod. Res. Develop. 12 (1973) 110–116. [3] S.E. Orchard, Surface levelling in viscous liquids and gels, Appl. Sci. Res. A 11 (1962) 451–464. [4] R.R. Iyer, D.W. Bousfield, The leveling of coating defects with shear thinning rheology, Chem. Eng. Sci. 51 (1996) 4611–4617. [5] R. Keunings, D.W. Bousfield, Analysis of surface tension driven leveling in viscoelastic films, J. Non-Newtonian Fluid Mech. 22 (1987) 219–233. [6] H.S. Kheshgi, L.E. Scriven, The evolution of disturbances in horizontal films, Chem. Eng. Sci. 43 (1988) 793–801. [7] H.A. Barnes, Thixotropy – a review, J. Non-Newtonian Fluid Mech. 70 (1997) 1–33. [8] J. Mewis, Thixotropy – a general review, J. Non-Newtonian Fluid Mech. 6 (1979) 1–20. [9] J. Mewis, N.J. Wagner, Thixotropy, Adv. Colloid Interface Sci. 147–148 (2009) 214–227. [10] J. Billingham, J.W.J. Ferguson, Laminar, unidirectional flow of a thixotropic fluid in a circular pipe, J. Non-Newtonian Fluid Mech. 47 (1993) 21–55. [11] D.C.-H. Cheng, F. Evans, Phenomenological characterization of the rheological behaviour of inelastic reversible thixotropic and antithixotropic fluids, Br. J. Appl. Phys. 16 (1965) 1599–1617. [12] H.T. Huynh, N. Roussel, P. Coussot, Aging and free surface flow of a thixotropic fluid, Phys. Fluids 17 (2005) 033101. [13] A. Mujumdar, A.N. Beris, A.B. Metzner, Transient phenomena in thixotropic systems, J. Non-Newtonian Fluid Mech. 102 (2002) 157–178. [14] O. Cohu, A. Magnin, The leveling of thixotropic coatings, Prog. Org. Coatings 28 (1996) 89–96. [15] F. Moore, The rheology of ceramic slips and bodies, Trans. Br. Ceramics Soc. 58 (1959) 470–494. [16] C. Baravanian, D. Quemada, A. Parker, Modelling thixotropy using a novel structural kinetics approach: basis and application to a solution of iota carrageenan, J. Texture Studies 27 (1996) 371–390. [17] M.H. Eres, L.W. Schwartz, R.V. Roy, Fingering phenomena for driven coating films, Phys. Fluids 12 (2000) 1278–1295. [18] D.E. Weidner, L.W. Schwartz, Contact-line motion of shear-thinning liquids, Phys. Fluids 6 (1994) 3535–3539.