- Email: [email protected]
Solidification of a liquid film on a vertical wall
Chemical
Engineering Science, 1975, Vol. 30, pp. 327-334. PergamonPress.
PrintedinGreat Britain
SOLIDIFICATION OF A LIQUID FILM ON A VERTICAL WALL J. ISENBERG and C. GUTF’INGER Laboratory for Coating Technology, Faculty of Mechanical Engineering, Technion-Israel Institute of Technology, Haifa, Israel (Received 5 March 1974;accepted 10September 1974) Abstract-Solidification during heat transfer to a draining film has been investigated. The mathematical formulation of the problem resulted in a coupled set of three no&&r partial differential equations in temperature, velocity, and film thickness as functions of position, time, and problem parameters. The set was solved by an implicit finite difference method for two sets of boundary conditions.A correlationfor the solidifiedfilmthicknessas a function of
position and physicalpropertiesof the film was found. The problem of drying of paint due to evaporation of volatile solvent is left for future effort.
INTRODUCTION The process of dip coating of solid objects with liquid
melts is interesting from the standpoint of heat transfer. Here, the liquid film adhering to the solid object drains while cooling and solidifying, thus resulting in a problem of unsteady motion and heat transfer with a large change in physical properties, especially viscosity. In the present paper we analyze this process for glass-like materials, that is to say, materials that have a gradual transition from liquid to solid that occurs within a range of temperatures, rather than at a single point. With materials such as these, the transition from liquid to solid takes place through the gradual increase in viscosity from some low value to a value large enough, so as to allow one to consider the material a solid, for heat transfer analysis. Indeed, window panes that have been removed from old European churches were thicker at the bottom than at the top, indicating that glass flows even at room temperatures. In the present problem, however, these super-cooled liquids will not differ in behavior from a true solid from the viewpoint of heat transfer. In addition, many polymer melts have solidification characteristics similar to those of glass. The present work is the second in a series of papers undertaken to investigate the process of dip-coating of objects with paint. The first stage of this investigation[l] was confined to the study of heat (or mass transfer) in a draining film without allowing for solidification or mass loss due to evaporation. In the present study, we account for the change in viscosity of the fihn resulting in solidification. In order to isolate the effect of variable viscosities, the other physical properties, such as thermal diffusivity, were intentionally assumed constant. It should be noted, however, that since the method of solution is numerical, variable properties can be incorporated directly into the solution as it is marched out.
PROBLEMSTATEMENT
We address ourselves to the problem of a draining fihn of an incompressible Newtonian fluid with a temperaturedependent viscosity. Figure 1 depicts the coordinate system used. The problem is defined by the following system of equations:
au,av=, ax
ay
Momentum
Energy aT aT aT at+uz+vdy=a($s$).
(3)
Film mass balance
as -=-at
a ax
s [I I. o
udy
As seen from Eq. (2) the term au/at was omitted from the momentum balance. This is standard procedure in treating creeping flows such as the present one. The film thickness dependence on time enters the problem through Eq. (4). Although there are several isothermal unsteady state solutions in the literature [7,8], it was shown by Gutfinger and Tallmadge that for all practical purposes the
327
.I. ISENBERGand C.
328
GUTFINGER
quasi-steady state solution is accurate. However, one cannot drop the unsteady state term from the energy equation; thus Eq. (3) contains three independent variables. We now assume a viscosity-temperature relation of the form known in the literature[2] as the Andrade equation [3]: v(T) = A exp (BIT).
(5)
This equation is known to be reasonably accurate in correlating data. T, in Eq. (5) is the absolute temperature. For this reason the temperature, T, in Eq. (3) will also be the absolute temperature. Since the fihn is thin in comparison to the wall height, L, one may show, using arguments similar to those of Ostrach[9], that the ratio of vertical to transverse conduction is of the order (S/L)2, thus rendering vertical conduction negligible. Eq. (3), therefore, simplifies to: $T+ a$+
2 ,aT= aa. 8~ ay
(6) I--I!I\ x
The above system of equations will now be recast into dimensionless forms by the following transformations:
Fig. 1. Schematicof the coordinatesystem.
The complete derivation of Eqs. (8)-(12) is given in (7)
Appendix A*(7)-( 12)will now be transformed into a more Equatrons tractable system by introducing a stream function, satisfying continuity, defined by the following relations:
I u&Y A an
Here the subscripts ‘0’ and ‘CO’ refer toinitial and intinite time, respectively. The dimensionless system of equations is now: --aU aX $+
n &AaiJ+ .l av_ o Aaxaq Aaq
(8)
~=-ag,3&!?%! aX AaXan’
(s’ UAdq)];
&;= A’(1- q)/p
(14)
The final form of the system of governing equations is:
ae v ae 1 a28 u-+--=-raX A a17 A a7 +z[ U$-&
(13)
ae+l ag ae+Aa(7jqs-+)ae 1 a% ~=ii5y ar AaqaX A aX
(19
q= A’(1- q)@ 817
(16)
(9) (10) where /3 is a function of 0 as defined by Eq. (II), and
(11)
a7
where C = B/To and $=-&(I,‘UAdq).
aA -=--
(12)
a& ax
(17)
where & is the value of $I at q = 1, and physically represents the flow rate past a given X. The boundary
Solidificationof a liquidfilmon a verticalwall
an asymptotic solution to be employed to solve for 8 in the vicinity of this point. A detailed description of the method used is given by Fox 1101.To apply this method it was assumed that at the beginning of the process the temperature and velocity profiles are identical to those for the constant viscosity modelll], the case of /3 = 1. A similarity variable 5 =7*/X was employed and an asymptotic solution for small values of 8 was obtained for Eqs. (15)-(U):
conditions for the system are: +=Oand--0
a+ 877
atq=O
0(x, II, 0) = 1
atX>O,OSn41
O(O,977) = 0
at 7>O,OSqSl
e(x, 0,T) = 0
at 7>O,X>O
$J,r)=O
at 7>0, X>O.
329
8=erf($)=erf(q) Here, Eq. (19) states that the film is initially at a uniform high temperature 6 = 1, while Eqs. (21) and (22) represent the solid wall kept at 0 = 0, and the insulated free surface, respectively. The insulated boundary condition at the free surface was chosen in preference to a convective boundary condition.because it results in a more simple analysis, while still retaining all the interesting facets of the problem. Moreover, the resulting solution provides an upper bound to the time of solidification. The system of Eqs. (15)-(17) with the boundary conditions (18)-(22) were solved numerically as described below. The set of Eqs. (15)-(17) is a coupled system in 13, $, and A. Equation (15) is initial-valued in X and T, and boundary-valued in TJand hence may be referred to as a parabolic-hyperbolic equation. Equation (16) is, in effect, an ordinary differential equation in 7, with the dependent variables $/A’ and 8, the latter entering through p. Equation (17) when rewritten in the form: (23) is recognized as a first-order hyperbolic equation in A with a nonlinear coefficient. METHOD OF SOLUTION The system of Eqs. (15)-(17) is nonlinear and coupled, and, therefore, it was solved iteratively. A finite difference scheme that was backward in both X and 7 was formulated for Eq. (15). This scheme was chosen in preference to the explicit formula[41 because of its unconditional stability. Equation (16) was solved by conventional finite difference methods for ordinary differential equations. Equation (17) was solved by an implicit scheme for first-order hyperbolic equations known as Wendrotf’s formula[5]. The solution to the system was found to be insensitive to variations in the grid size. The finite difference equations and the method of solution are outlined in Appendix B. The computation began by using an asymptotic solution. A discontinuity in temperature exists at the corner of the boundary as X+0 and T-0. This requires
(24) (25)
(26) In the computations Eqs. (24)-(26) were used until 5 = ?/X reached a value of 0.005 at which point the computation reverted to the original Eqs. (15)-(17). Test runs with other transitional values of 5 showed no significant changes in the results. RESULTS From the solution of Eqs. (15)-(17) one obtains temperature, velocity, and film thickness profiles as functions of time, for various values of the parameters C and T-/T0 which appear in the viscosity relationship, Eq. (11). For the sake of convenience, the parameter C was replaced by a new parameter /3- defined as B.=E=exp[C(g-l)].
(27)
Here v- is the viscosity corresponding to that attained at T+ m, (i.e. room temperature, T, or 0 = 0). For an input value of /3- one may calculate C by using Eq. (27). The pertinent information that can be obtained from the numerical solution described above may be summarized as follows: U= U(x, q, 7; TJTo, /L)
(28)
0= 0(x, 7,~; T-/To, P-)
(29)
A = A(X, T; T-/T,, &).
(30)
Some typical results of the form of Eqs. (28) and (29) are presented in Figs. 2 and 3, respectively. These figures were obtained for the values p- = ld’ and T-/T,, = O-143, which were calculated from viscosity-temperature data given for commercial soda-lime-silica glass (hereafter referred to as ‘glass’) in the literature[6]. The physical properties for this glass are given in Table 1. Figure 2 shows velocity profiles for solidifying glass at a
330
J.
ISENBERG
and C. GUTPINGER
value of X = 1000, which approximately corresponds to x = l/2 ft. The velocities predicted by the constant viscosity model are obtained by applying Eq. (13) to Eq. (26) which yields the parabolic relation: U=A*q(l-42).
(31)
As seen in the figure, the actual velocity profiles are parabolic near the free surface but decay rapidly toward zero while moving toward the wall.
Figure 3 shows representative temperature profiles f solidifying glass for three values of X: 100, 1000, al 10,000. In comparing these results with those for tl constant viscosity model[l], it is seen that the profil have a similar shape. However, for the constant viscosi case the profiles fall more rapidly to zero. For exampl the temperature at X = 1000and 7 = 1 at time 7 = 55 fro Fig. 3b is 0 = 0.675. For the constant viscosity model tl temperature for the same values of the independe variables is 0 = 0408. This lower rate of cooling is due
Table 1. Viscosity temperature data for soda-limesilica glass
x= 1000
U
2,0-
I.O-
0
0.2
0.4
0.8
9
100
“C
lo2 3.16 x lo’ 10” 10’ 10’ lo6 10’ lo8 lo’ 1O’O 10” lOI2 10L3 10” 10”
1451 1295 1178 1013 903 823 764 716 674 639 609 583 559 539 523
Note: Viscosity data below 523°C were extrapolated using the Andrade equation, Eq. (5).
I.0
Fig. 2. Typical velocity profiles for solidifying glass. a:X-
Viscosity (P)
b:X=lOOO
c:x= 10000
Fig. 3. Representative temperature profiles for solidifying glass for different values of X.
Solidificationof a liquidfdmon a verticalwall
IO".
I
I
10-j
I
o”
I
331
I
I
IO2
0’
IO3
IO4
X Fig. 4. Filmthicknessprofilesfor solidifyingglassfor variousvaluesof 7. the lower velocities characterizing the variable viscosity model. Figure 4 gives dimensionless film thickness for glass as a function of dimensionless distance from the leading edge X for various values of time T. It is seen that when the film behaves as a liquid as represented by T = 0.01 and the upper portion of the curve for T = 0.10, the slope on the log-log scale is l/2, in accordance with Eq. (29, the solution for the constant viscosity model [7]. By contrast, the final solidification profile given by r greater than 216 has a slope of l/4. This result was rather unexpected and was, therefore, further investigated for other combinations of /3- and T-/To. The results shown in Fig. 5 all give the same slope of l/4. This indicates a correlation for final
I
I
I
*,
b$3,_10s5
I I
6’
lo”
I
I
IO'
I
o2
A, = FX”4
(32)
or, in dimensional form, as: 114
(7 1 S-=F where
F = F(&, Tm/T,,).
(34)
The correlation obtained in Eq. (33) for the final film thickness is very interesting as it indicates a simple dependence of this final thickness on the problem parameters. Due to the nonlinearity of the problem, it is difficult to assess the generality of this relationship. It should, however, be emphasized that the correlation holds for the two extremes in boundary conditions treated in this paper. In Fig. 6 a semi-log plot of F against T-/T0 is shown for two values of /L. Equation (33) together with Fig. 6 define the range of variation of the film thickness within the limits of practical interest. The method of solution given in this study could obviously be used for different types of boundary conditions generating results equivalent to those above. To demonstrate this concept, the boundary at the wall and free surface were interchanged. The new boundary conditions replacing Eqs. (21) and (22) are
I
IO3
IO'
X Fig. 5. Solidiliedprofilesfor variouscombinationsof &and TAT,.
CES Vol. M, No. 3-E
film thickness vs X, of the form:
$X,0,7)=0 0(X,1,7)=0
at T>O,X>O at
7>O,X>O.
(35) (36)
332
J. ISENBERG and C.
GUTFINGER
Fig. 6. Correlation variable F as a function of T-/T, for two values of /L.
X
Fig. 7. Film thickness protiles for glass solidifying on an insulated wall for various values of 7.
Figure 7 is the modified boundary condition equivalent to Fig. 4 presented earlier for film thickness profiles for various values of r. Here under identical conditions, the thicknesses obtained are equal to or lower than in the previous case. However, the curve for final film thickness to which all the curves in Fig. 7 converge is again a straight line of l/4 slope. Figure 8 provides information on solidification time for glass for the two cases discussed above. The criterion of solidification was taken arbitrarily as: /3 = o-9/3_.
(37)
Under this condition, the material can be considered a solid for all practical purposes. The resulting curves are essentially straight lines with the same slope of about l/2. SUMMARY
(1) The problem of a solidifying film draining down a vertical wall has been formulated for a Newtonian fluid. (2) The resulting formulation consisted of a coupled set of three nonlinear partial differential equations. (3) The set of equations was solved by an implicit finite difference method for two sets of boundary conditions. (4) The results obtained include a l/4 power correlation
333
Solidificationof a liquid 8lm on a vertical wall
Fig. 8. Solidificationtime for glass as a function of X, for two types of boundaryconditions.
that gives the final solidified thickness as a function of position and material properties. (5) The present paper provides a step forward toward the mathematical analysis of coating processes. NOTATION
A,3 coefficients in viscosity-temperature relation, Eq. (5) c coefficient, BIT0 F correlating function for A, vs X gravitational constant ! wall height t time T absolute temperature u, v velocity components in x and y directions, respectively u, v dimensionless velocity components in X and q directions, respectively, defined by Eq. (7) distance from leading edge ; dimensionless distance from leading edge, I/3 kl(~~o)l x
Y transverse coordinate
Greek symbols thermal diiusivity dimensionless viscosity, v/v0 film thickness dimensionless film thickness dimensionless transverse coordinate, y/6 dimensionless temperature, (T - T-)1(TO- T-) kinematic viscosity similarity variable, 7*/X dimensionless time, (g’a/A)“‘t stream function defined by Eqs. (13) and (14) stream function at free surface
Subscripts 0 initial, 7 =0 m final, 7-m
[II Isenberg J. and Gutfinger C., Znr..Z.Heat Mass Tmnsfer, 1973 16 505.
PI Andrade E. N. da& Endeavour, 195413 117.
and r31 Reid R. C. and Sherwood T. K., The Properties ofGases Liquids, p. 431, 2nd Ed., McGraw-Hill, New York, 1966. I41 Hellums J. D. and Churchill, S. W., Znr. Dev. in Heat Transfer, l%l 5 985. VI Mitchell A. R., ComputationalMethods in Partial Differential Equations. p. 165,John Wiley, London, 1%9. 161Cole H. (supervisor), Viscosiry Temperature Relations in Glass. p. 24, International Commission on Glass, S. C. MaisonDedition,Belgium,1970. Guttinger C. and TaBmadge J. A., A.Z.Ch.E.Jl, 196410774. Denson C. D., Trans. Sot. Rheology, 197216 697. Ostrach S., Role of analysis in the solution of complex physical problems. Third International Heat Tmnsfer Conference, Chicago, Blinois, 1966. ml Fox L. (editor),Numerical Solution of Ordinaryand Partial DifFcrentialEauations. DD. 252-254,Addison-Wesley. _ Readin;, Massachusetts, l%y t111 WendrofI B., 3. Sot. Zndust. appl. Math., 19688 549. . APPENOIX A
Derivation of governingdimensionless equations In deriving the dimensionless equations we transform.from x, y, r, coordinates to X, 1). 7, coordinates defined by Eq. (7). The transformations of x and y are:
(A.11 (A.3 (A.31
J. ISENEIERG and C. GUTFINGER
334
convenience we shah rewrite the system (15)-(17):
The time derivative is transformed as follows: ;=
(g)-$+
(~)-&=$;-;+
(A.4)
The term &/at is now replaced by the expression in Eq. (4). The resulting expression, put in dimensionless form, is: &= ($)“‘[$+;-&(1.‘UAdg);].
The equation of continuity [Eq. (I)] is obtained with the help of (A.l) and (A.2) resulting in: aV_llahaU,laU=O ax h ax a7 A a7
$tF$=O
(A.5)
(8)
.
where A, B, C, E, and F are nonlinear variable coefficients, The solution region is divided into a grid of mesh size specl by L, hx, and h,. For each point the finite differ< approximations to (B.l), (B.2), and (B.3), are, respectively: n+,
n
et+,.,
0 ,+,.,h,
Integrating Eq. (2) with the boundary condition &lay = 0 at y = 8 gives: (A.6) Here v(T) is given by Eq. (5), which can be put in dimensionless form as:
. 11
(A.7)
Substituting (A.7) into (A.6) and applying Eq. (7) yields the dimensionless momentum equation:
and A~:,‘-A;-~+A;+l-Ai.tF~(A;:,‘-Al”+itA;+,-AI”)=O. x
erJ= ce*‘e’” E;= A*(1- n)&?.
(
(10)
The film mass balance, Eq. (4), is put in dimensionless form using Eqs. (A.3) and (7): $=-&(I,‘IJAdq).
(
To determine the stability of Eq. (B.4) we must investigate propagation of errors in the r and X directions. To do this we the von Neumann stability criterion. We assume that the roundoff errors can be represented 1 Fourier series. For the r propagation of errors, a typical terr the series is represented by:
where L=d/-1 and { is the amplification factor. To en: stability and convergence, the absolute value of f must not exe unity. Substituting c for 0 in (B.4) and dropping the indices from coefficients, we obtain:
(12)
One should note that the second term on the left hand side of Eqs. (A.l) and (A.5)does not enter into (12)as 6 is only a function of X and r, and not of n. The energy equation, Eq. (6), is put in dimensionless form using Eqs. (A.l), (A.2), (A.3), (A.5), and (A.7), resulting in:
(9) APPENDIX B Finite difference approximations of the governing equations We denote the coordinates X and q by the subscripts i and j,
respectively, and r by the superscript n. For the sake of
= 2C$e’b(cos y - 1) .I
which can be readily shown to give a maximum absolute vahr unity for t, when p = y = 0. A similar conclusion can be obta in studying the propagation of errors in the X-direction. ‘methodis, therefore, unconditionally stable for all combinat of h, hx, and b,. The finite difference equation was solved I conventional method known as the Thomas algorithm. Equation (B.5) is solved by applying the boundary condit + = ag/aq = 0 at q = 0. At q = 0 this yields Jr0= 0, +, = h,,’ The ensuing points are solved by recurrence via Eq. (B.5). As noted earher, the stability proof for Eq. (B.6) is giver Wendroff [ll].
Engineering Science, 1975, Vol. 30, pp. 327-334. PergamonPress.
PrintedinGreat Britain
SOLIDIFICATION OF A LIQUID FILM ON A VERTICAL WALL J. ISENBERG and C. GUTF’INGER Laboratory for Coating Technology, Faculty of Mechanical Engineering, Technion-Israel Institute of Technology, Haifa, Israel (Received 5 March 1974;accepted 10September 1974) Abstract-Solidification during heat transfer to a draining film has been investigated. The mathematical formulation of the problem resulted in a coupled set of three no&&r partial differential equations in temperature, velocity, and film thickness as functions of position, time, and problem parameters. The set was solved by an implicit finite difference method for two sets of boundary conditions.A correlationfor the solidifiedfilmthicknessas a function of
position and physicalpropertiesof the film was found. The problem of drying of paint due to evaporation of volatile solvent is left for future effort.
INTRODUCTION The process of dip coating of solid objects with liquid
melts is interesting from the standpoint of heat transfer. Here, the liquid film adhering to the solid object drains while cooling and solidifying, thus resulting in a problem of unsteady motion and heat transfer with a large change in physical properties, especially viscosity. In the present paper we analyze this process for glass-like materials, that is to say, materials that have a gradual transition from liquid to solid that occurs within a range of temperatures, rather than at a single point. With materials such as these, the transition from liquid to solid takes place through the gradual increase in viscosity from some low value to a value large enough, so as to allow one to consider the material a solid, for heat transfer analysis. Indeed, window panes that have been removed from old European churches were thicker at the bottom than at the top, indicating that glass flows even at room temperatures. In the present problem, however, these super-cooled liquids will not differ in behavior from a true solid from the viewpoint of heat transfer. In addition, many polymer melts have solidification characteristics similar to those of glass. The present work is the second in a series of papers undertaken to investigate the process of dip-coating of objects with paint. The first stage of this investigation[l] was confined to the study of heat (or mass transfer) in a draining film without allowing for solidification or mass loss due to evaporation. In the present study, we account for the change in viscosity of the fihn resulting in solidification. In order to isolate the effect of variable viscosities, the other physical properties, such as thermal diffusivity, were intentionally assumed constant. It should be noted, however, that since the method of solution is numerical, variable properties can be incorporated directly into the solution as it is marched out.
PROBLEMSTATEMENT
We address ourselves to the problem of a draining fihn of an incompressible Newtonian fluid with a temperaturedependent viscosity. Figure 1 depicts the coordinate system used. The problem is defined by the following system of equations:
au,av=, ax
ay
Momentum
Energy aT aT aT at+uz+vdy=a($s$).
(3)
Film mass balance
as -=-at
a ax
s [I I. o
udy
As seen from Eq. (2) the term au/at was omitted from the momentum balance. This is standard procedure in treating creeping flows such as the present one. The film thickness dependence on time enters the problem through Eq. (4). Although there are several isothermal unsteady state solutions in the literature [7,8], it was shown by Gutfinger and Tallmadge that for all practical purposes the
327
.I. ISENBERGand C.
328
GUTFINGER
quasi-steady state solution is accurate. However, one cannot drop the unsteady state term from the energy equation; thus Eq. (3) contains three independent variables. We now assume a viscosity-temperature relation of the form known in the literature[2] as the Andrade equation [3]: v(T) = A exp (BIT).
(5)
This equation is known to be reasonably accurate in correlating data. T, in Eq. (5) is the absolute temperature. For this reason the temperature, T, in Eq. (3) will also be the absolute temperature. Since the fihn is thin in comparison to the wall height, L, one may show, using arguments similar to those of Ostrach[9], that the ratio of vertical to transverse conduction is of the order (S/L)2, thus rendering vertical conduction negligible. Eq. (3), therefore, simplifies to: $T+ a$+
2 ,aT= aa. 8~ ay
(6) I--I!I\ x
The above system of equations will now be recast into dimensionless forms by the following transformations:
Fig. 1. Schematicof the coordinatesystem.
The complete derivation of Eqs. (8)-(12) is given in (7)
Appendix A*(7)-( 12)will now be transformed into a more Equatrons tractable system by introducing a stream function, satisfying continuity, defined by the following relations:
I u&Y A an
Here the subscripts ‘0’ and ‘CO’ refer toinitial and intinite time, respectively. The dimensionless system of equations is now: --aU aX $+
n &AaiJ+ .l av_ o Aaxaq Aaq
(8)
~=-ag,3&!?%! aX AaXan’
(s’ UAdq)];
&;= A’(1- q)/p
(14)
The final form of the system of governing equations is:
ae v ae 1 a28 u-+--=-raX A a17 A a7 +z[ U$-&
(13)
ae+l ag ae+Aa(7jqs-+)ae 1 a% ~=ii5y ar AaqaX A aX
(19
q= A’(1- q)@ 817
(16)
(9) (10) where /3 is a function of 0 as defined by Eq. (II), and
(11)
a7
where C = B/To and $=-&(I,‘UAdq).
aA -=--
(12)
a& ax
(17)
where & is the value of $I at q = 1, and physically represents the flow rate past a given X. The boundary
Solidificationof a liquidfilmon a verticalwall
an asymptotic solution to be employed to solve for 8 in the vicinity of this point. A detailed description of the method used is given by Fox 1101.To apply this method it was assumed that at the beginning of the process the temperature and velocity profiles are identical to those for the constant viscosity modelll], the case of /3 = 1. A similarity variable 5 =7*/X was employed and an asymptotic solution for small values of 8 was obtained for Eqs. (15)-(U):
conditions for the system are: +=Oand--0
a+ 877
atq=O
0(x, II, 0) = 1
atX>O,OSn41
O(O,977) = 0
at 7>O,OSqSl
e(x, 0,T) = 0
at 7>O,X>O
$J,r)=O
at 7>0, X>O.
329
8=erf($)=erf(q) Here, Eq. (19) states that the film is initially at a uniform high temperature 6 = 1, while Eqs. (21) and (22) represent the solid wall kept at 0 = 0, and the insulated free surface, respectively. The insulated boundary condition at the free surface was chosen in preference to a convective boundary condition.because it results in a more simple analysis, while still retaining all the interesting facets of the problem. Moreover, the resulting solution provides an upper bound to the time of solidification. The system of Eqs. (15)-(17) with the boundary conditions (18)-(22) were solved numerically as described below. The set of Eqs. (15)-(17) is a coupled system in 13, $, and A. Equation (15) is initial-valued in X and T, and boundary-valued in TJand hence may be referred to as a parabolic-hyperbolic equation. Equation (16) is, in effect, an ordinary differential equation in 7, with the dependent variables $/A’ and 8, the latter entering through p. Equation (17) when rewritten in the form: (23) is recognized as a first-order hyperbolic equation in A with a nonlinear coefficient. METHOD OF SOLUTION The system of Eqs. (15)-(17) is nonlinear and coupled, and, therefore, it was solved iteratively. A finite difference scheme that was backward in both X and 7 was formulated for Eq. (15). This scheme was chosen in preference to the explicit formula[41 because of its unconditional stability. Equation (16) was solved by conventional finite difference methods for ordinary differential equations. Equation (17) was solved by an implicit scheme for first-order hyperbolic equations known as Wendrotf’s formula[5]. The solution to the system was found to be insensitive to variations in the grid size. The finite difference equations and the method of solution are outlined in Appendix B. The computation began by using an asymptotic solution. A discontinuity in temperature exists at the corner of the boundary as X+0 and T-0. This requires
(24) (25)
(26) In the computations Eqs. (24)-(26) were used until 5 = ?/X reached a value of 0.005 at which point the computation reverted to the original Eqs. (15)-(17). Test runs with other transitional values of 5 showed no significant changes in the results. RESULTS From the solution of Eqs. (15)-(17) one obtains temperature, velocity, and film thickness profiles as functions of time, for various values of the parameters C and T-/T0 which appear in the viscosity relationship, Eq. (11). For the sake of convenience, the parameter C was replaced by a new parameter /3- defined as B.=E=exp[C(g-l)].
(27)
Here v- is the viscosity corresponding to that attained at T+ m, (i.e. room temperature, T, or 0 = 0). For an input value of /3- one may calculate C by using Eq. (27). The pertinent information that can be obtained from the numerical solution described above may be summarized as follows: U= U(x, q, 7; TJTo, /L)
(28)
0= 0(x, 7,~; T-/To, P-)
(29)
A = A(X, T; T-/T,, &).
(30)
Some typical results of the form of Eqs. (28) and (29) are presented in Figs. 2 and 3, respectively. These figures were obtained for the values p- = ld’ and T-/T,, = O-143, which were calculated from viscosity-temperature data given for commercial soda-lime-silica glass (hereafter referred to as ‘glass’) in the literature[6]. The physical properties for this glass are given in Table 1. Figure 2 shows velocity profiles for solidifying glass at a
330
J.
ISENBERG
and C. GUTPINGER
value of X = 1000, which approximately corresponds to x = l/2 ft. The velocities predicted by the constant viscosity model are obtained by applying Eq. (13) to Eq. (26) which yields the parabolic relation: U=A*q(l-42).
(31)
As seen in the figure, the actual velocity profiles are parabolic near the free surface but decay rapidly toward zero while moving toward the wall.
Figure 3 shows representative temperature profiles f solidifying glass for three values of X: 100, 1000, al 10,000. In comparing these results with those for tl constant viscosity model[l], it is seen that the profil have a similar shape. However, for the constant viscosi case the profiles fall more rapidly to zero. For exampl the temperature at X = 1000and 7 = 1 at time 7 = 55 fro Fig. 3b is 0 = 0.675. For the constant viscosity model tl temperature for the same values of the independe variables is 0 = 0408. This lower rate of cooling is due
Table 1. Viscosity temperature data for soda-limesilica glass
x= 1000
U
2,0-
I.O-
0
0.2
0.4
0.8
9
100
“C
lo2 3.16 x lo’ 10” 10’ 10’ lo6 10’ lo8 lo’ 1O’O 10” lOI2 10L3 10” 10”
1451 1295 1178 1013 903 823 764 716 674 639 609 583 559 539 523
Note: Viscosity data below 523°C were extrapolated using the Andrade equation, Eq. (5).
I.0
Fig. 2. Typical velocity profiles for solidifying glass. a:X-
Viscosity (P)
b:X=lOOO
c:x= 10000
Fig. 3. Representative temperature profiles for solidifying glass for different values of X.
Solidificationof a liquidfdmon a verticalwall
IO".
I
I
10-j
I
o”
I
331
I
I
IO2
0’
IO3
IO4
X Fig. 4. Filmthicknessprofilesfor solidifyingglassfor variousvaluesof 7. the lower velocities characterizing the variable viscosity model. Figure 4 gives dimensionless film thickness for glass as a function of dimensionless distance from the leading edge X for various values of time T. It is seen that when the film behaves as a liquid as represented by T = 0.01 and the upper portion of the curve for T = 0.10, the slope on the log-log scale is l/2, in accordance with Eq. (29, the solution for the constant viscosity model [7]. By contrast, the final solidification profile given by r greater than 216 has a slope of l/4. This result was rather unexpected and was, therefore, further investigated for other combinations of /3- and T-/To. The results shown in Fig. 5 all give the same slope of l/4. This indicates a correlation for final
I
I
I
*,
b$3,_10s5
I I
6’
lo”
I
I
IO'
I
o2
A, = FX”4
(32)
or, in dimensional form, as: 114
(7 1 S-=F where
F = F(&, Tm/T,,).
(34)
The correlation obtained in Eq. (33) for the final film thickness is very interesting as it indicates a simple dependence of this final thickness on the problem parameters. Due to the nonlinearity of the problem, it is difficult to assess the generality of this relationship. It should, however, be emphasized that the correlation holds for the two extremes in boundary conditions treated in this paper. In Fig. 6 a semi-log plot of F against T-/T0 is shown for two values of /L. Equation (33) together with Fig. 6 define the range of variation of the film thickness within the limits of practical interest. The method of solution given in this study could obviously be used for different types of boundary conditions generating results equivalent to those above. To demonstrate this concept, the boundary at the wall and free surface were interchanged. The new boundary conditions replacing Eqs. (21) and (22) are
I
IO3
IO'
X Fig. 5. Solidiliedprofilesfor variouscombinationsof &and TAT,.
CES Vol. M, No. 3-E
film thickness vs X, of the form:
$X,0,7)=0 0(X,1,7)=0
at T>O,X>O at
7>O,X>O.
(35) (36)
332
J. ISENBERG and C.
GUTFINGER
Fig. 6. Correlation variable F as a function of T-/T, for two values of /L.
X
Fig. 7. Film thickness protiles for glass solidifying on an insulated wall for various values of 7.
Figure 7 is the modified boundary condition equivalent to Fig. 4 presented earlier for film thickness profiles for various values of r. Here under identical conditions, the thicknesses obtained are equal to or lower than in the previous case. However, the curve for final film thickness to which all the curves in Fig. 7 converge is again a straight line of l/4 slope. Figure 8 provides information on solidification time for glass for the two cases discussed above. The criterion of solidification was taken arbitrarily as: /3 = o-9/3_.
(37)
Under this condition, the material can be considered a solid for all practical purposes. The resulting curves are essentially straight lines with the same slope of about l/2. SUMMARY
(1) The problem of a solidifying film draining down a vertical wall has been formulated for a Newtonian fluid. (2) The resulting formulation consisted of a coupled set of three nonlinear partial differential equations. (3) The set of equations was solved by an implicit finite difference method for two sets of boundary conditions. (4) The results obtained include a l/4 power correlation
333
Solidificationof a liquid 8lm on a vertical wall
Fig. 8. Solidificationtime for glass as a function of X, for two types of boundaryconditions.
that gives the final solidified thickness as a function of position and material properties. (5) The present paper provides a step forward toward the mathematical analysis of coating processes. NOTATION
A,3 coefficients in viscosity-temperature relation, Eq. (5) c coefficient, BIT0 F correlating function for A, vs X gravitational constant ! wall height t time T absolute temperature u, v velocity components in x and y directions, respectively u, v dimensionless velocity components in X and q directions, respectively, defined by Eq. (7) distance from leading edge ; dimensionless distance from leading edge, I/3 kl(~~o)l x
Y transverse coordinate
Greek symbols thermal diiusivity dimensionless viscosity, v/v0 film thickness dimensionless film thickness dimensionless transverse coordinate, y/6 dimensionless temperature, (T - T-)1(TO- T-) kinematic viscosity similarity variable, 7*/X dimensionless time, (g’a/A)“‘t stream function defined by Eqs. (13) and (14) stream function at free surface
Subscripts 0 initial, 7 =0 m final, 7-m
[II Isenberg J. and Gutfinger C., Znr..Z.Heat Mass Tmnsfer, 1973 16 505.
PI Andrade E. N. da& Endeavour, 195413 117.
and r31 Reid R. C. and Sherwood T. K., The Properties ofGases Liquids, p. 431, 2nd Ed., McGraw-Hill, New York, 1966. I41 Hellums J. D. and Churchill, S. W., Znr. Dev. in Heat Transfer, l%l 5 985. VI Mitchell A. R., ComputationalMethods in Partial Differential Equations. p. 165,John Wiley, London, 1%9. 161Cole H. (supervisor), Viscosiry Temperature Relations in Glass. p. 24, International Commission on Glass, S. C. MaisonDedition,Belgium,1970. Guttinger C. and TaBmadge J. A., A.Z.Ch.E.Jl, 196410774. Denson C. D., Trans. Sot. Rheology, 197216 697. Ostrach S., Role of analysis in the solution of complex physical problems. Third International Heat Tmnsfer Conference, Chicago, Blinois, 1966. ml Fox L. (editor),Numerical Solution of Ordinaryand Partial DifFcrentialEauations. DD. 252-254,Addison-Wesley. _ Readin;, Massachusetts, l%y t111 WendrofI B., 3. Sot. Zndust. appl. Math., 19688 549. . APPENOIX A
Derivation of governingdimensionless equations In deriving the dimensionless equations we transform.from x, y, r, coordinates to X, 1). 7, coordinates defined by Eq. (7). The transformations of x and y are:
(A.11 (A.3 (A.31
J. ISENEIERG and C. GUTFINGER
334
convenience we shah rewrite the system (15)-(17):
The time derivative is transformed as follows: ;=
(g)-$+
(~)-&=$;-;+
(A.4)
The term &/at is now replaced by the expression in Eq. (4). The resulting expression, put in dimensionless form, is: &= ($)“‘[$+;-&(1.‘UAdg);].
The equation of continuity [Eq. (I)] is obtained with the help of (A.l) and (A.2) resulting in: aV_llahaU,laU=O ax h ax a7 A a7
$tF$=O
(A.5)
(8)
.
where A, B, C, E, and F are nonlinear variable coefficients, The solution region is divided into a grid of mesh size specl by L, hx, and h,. For each point the finite differ< approximations to (B.l), (B.2), and (B.3), are, respectively: n+,
n
et+,.,
0 ,+,.,h,
Integrating Eq. (2) with the boundary condition &lay = 0 at y = 8 gives: (A.6) Here v(T) is given by Eq. (5), which can be put in dimensionless form as:
. 11
(A.7)
Substituting (A.7) into (A.6) and applying Eq. (7) yields the dimensionless momentum equation:
and A~:,‘-A;-~+A;+l-Ai.tF~(A;:,‘-Al”+itA;+,-AI”)=O. x
erJ= ce*‘e’” E;= A*(1- n)&?.
(
(10)
The film mass balance, Eq. (4), is put in dimensionless form using Eqs. (A.3) and (7): $=-&(I,‘IJAdq).
(
To determine the stability of Eq. (B.4) we must investigate propagation of errors in the r and X directions. To do this we the von Neumann stability criterion. We assume that the roundoff errors can be represented 1 Fourier series. For the r propagation of errors, a typical terr the series is represented by:
where L=d/-1 and { is the amplification factor. To en: stability and convergence, the absolute value of f must not exe unity. Substituting c for 0 in (B.4) and dropping the indices from coefficients, we obtain:
(12)
One should note that the second term on the left hand side of Eqs. (A.l) and (A.5)does not enter into (12)as 6 is only a function of X and r, and not of n. The energy equation, Eq. (6), is put in dimensionless form using Eqs. (A.l), (A.2), (A.3), (A.5), and (A.7), resulting in:
(9) APPENDIX B Finite difference approximations of the governing equations We denote the coordinates X and q by the subscripts i and j,
respectively, and r by the superscript n. For the sake of
= 2C$e’b(cos y - 1) .I
which can be readily shown to give a maximum absolute vahr unity for t, when p = y = 0. A similar conclusion can be obta in studying the propagation of errors in the X-direction. ‘methodis, therefore, unconditionally stable for all combinat of h, hx, and b,. The finite difference equation was solved I conventional method known as the Thomas algorithm. Equation (B.5) is solved by applying the boundary condit + = ag/aq = 0 at q = 0. At q = 0 this yields Jr0= 0, +, = h,,’ The ensuing points are solved by recurrence via Eq. (B.5). As noted earher, the stability proof for Eq. (B.6) is giver Wendroff [ll].