3444
J. Phys. Chem. 1981, 85,3444-3448
Diffusion with Chemical Reaction and a Moving Boundary J. R. Kirkpatrlck Union Garb& Gorporatlon, Nuclear Dlvlslon, Oak Ridge, Tennessee 37830 (Received: December 3, 198 1; In Final Form: July 17, 198 1)
A solid sample is immersed in an atmosphere containing a corrosive material. The corrosive atoms diffuse through the sample and, at some point, undergo an irreversible combination with the sample material. The resulting compound has considerably different physical properties from the uncorroded sample so that the corroded layers break or otherwise become so porous that the corrosive atmosphere penetrates to the uncorroded portion of the sample. An approximate analytical solution to the equations describing the process is developed under certain assumptions. Expressions for time scales, corrosion front velocity, weight gain, and weight-gain rate are developed. These are compared to published experimental results for the uranium hydride system. Consequences of relaxing some of the assumptions used in the analysis are discussed. Introduction The model described in this paper was originally developed for calculating the hydriding of uranium and was published in a basic form in ref 1. In the years since ref 1was issued, the model has been extended and has been applied to other corrosion reactions besides uranium hydriding.2i3 Through this process of continued use and development, both the utility of the model and its limits have become more apparent. In order for the model to adequately describe a particular corrosion process, the physical system must satisfy a small number of crucial conditions. Some work has been done to show that information may be gained from this model concerning systems which do not meet all of these conditions and to calculate the behavior of such systems. The physical process being studied has the following characteristics. A sample of solid material B is exposed to an atmosphere containing another material A which is soluble in B. The atoms of A diffuse from the surface through the solid B and, at some point, undergo an irreversible chemical reaction nA + m B A,B,. Such a common process has, of course, been examined elsewhere in the literature, for example, ref 4 and 5. The crucial difference between the present paper and most other work is the presence of a moving outer surface. In the present paper, it is assumed that there is a severe enough change in the physical properties of the sample that, when a certain fraction of the solid B has been converted to the AB compound, a layer of the sample breaks or otherwise becomes sufficiently porous so that the atmosphere containing A penetrates to this point. Thus, the diffusion mechanism is bypassed for the damaged layers, and the diffusion problem must be solved with this fracture front as the new outer surface, Tracking the location of the front becomes an important segment of the overall problem of calculating corrosion rates. The linear model described in this paper was developed for hydriding of uranium metal in an atmosphere containing hydrogen. Results from a finite difference solution to this problem were published in ref 6. The linear model
-
(1) Kirkpatrick, J. R. “An Exact Solution of Diffusion and Absorption of Hydrogen in Uranium Followed by Fracture of the Hydride”, Union Carbide Corp., Nuclear Division, Report K/CSD-9, Nov 1977, available from National Technical Information Service, Springfield, VA. (2) Condon, J. B. J . Less-Common Met. 1980, 73, 105. (3) Condon, J. B.; Cristy, S. S. “A Model of Oxidation of Silicon by Oxygen”, Union Carbide Corp., Nuclear Division, Report Y /DW-148, July 1980, submitted to J. Electrochem. Sac. (4) Roberta, G. Met. Sci. 1979, 13, 94. (5) Bird, R. B., e t al. “Transport Phenomena”; Wiley: New York, 1960. 0022-3654/81/2085-3444$01.25/0
is an approximate analytic solution which, so long as the system meets its assumptions, gives the same answers as the finite difference scheme with much less work. The linear model is an approximate solution to the entire corrosion process. It describes the progress of the fracture front, the diffusion of material A ahead of the front, and the inclusion of material A into the AB compound both ahead of and behind the front. The uranium hydride system has a number of characteristics which allow a satisfactory match to experimental data using this model. The more important of these characteristics are the following: (1) The hydride breaks away from the sample (spalls off) by the time that roughly 3% of the uranium has been converted to hydride. (2) The solubility of hydrogen in uranium is low. A mathematical expression will be developed showing what magnitude is to be considered “low”.
Basic Equations The basic mathematical model for the diffusion-plusreaction process is known as “perfusive precipitation”.2 The governing equations for unfractured solids in plane geometry are6
ac = D-a2c + a-au at aX2 at
au/at = - k c u
Ob)
where x is the distance from the surface where the atmosphere containing gaseous A meets the solid B, c is the concentration of atoms of the diffusing corrosive material A, U is the concentration of the immobile species B which is being corroded, D is the diffusivity of A in solid B, a is the stoichiometric ratio of A to B in the AB compound, and k is a rate constant for the corrosion reaction. The initial conditions are
c(x,O) = 0 (IC) U(x,O) = u, (Id) The boundary condition for a semiinfinite slab is c(0,t) = co where co is the solubility concentration of A in solid B. This condition represents the ability of the gaseous A in the atmosphere to saturate the solid lattice at the surface where the two phases meet ( x = 0). The presumption of spa11 requires a modification of the boundary condition. (6) Condon, J. B. J . Phys. Chem. 1975, 79, 392.
0 1981 American Chemical Society
Diffusion with Chemical Reaction
The Journal of Physical ChemWy, Vol. 85,
No. 23, 198 1 3445
TABLE I: Parameters for the U-H System as a Function of Temperature
T, " C
D, mzs-l 1.05 x 1 0 4 5 2.02 x 1.13 X lo-"'
0 150 325
ah,s-l
m 1.4 x 10-9 1.8 X lo-' 1.6 x Xo.01,
10600 1345 447
to.99(Xo,o,),
s
0.0006 0.004 0.013
t1,O s
0.17 0.42 0.68
Assumes 1 atm of hydrogen pressure and UJU0= 0.97.
Initially, the spall front is at the surface. When the corrosion has gone on for a sufficient length of time, the front begins to move into the sample. The boundary condition must account for the assumption that the atmosphere penetrates through the spalled material to the front. The modified boundary condition is c(x,(t),t) = co (le) where x , is the location of the front. Mathematically, x , is defined as the value of x for which U = Uswhere Us is the value of U for which fracture is assumed to occur. In the region x > xs, the U and c curves are found by solving the perfusive precipitation equations l a and l b subject to the boundary condition le. The region x 5 x , is presumed to contain spalled material for which c = cP The corrosion process in this region is calculated by using eq l b by itself. Values of c and U are commonly expressed as mole fractions. However, this is not necessary. The concentrations may be expressed in other units as long as k is in consistent units. If U is a mole fraction, then Uo= 1. Fits for D and k as functions of temperature and for co as a function of temperature and gas pressure are given in ref 6 for the uranium hydride system.
Approximate Solution to Basic Equations There is a solution of eq la-e for a semiinfinite slab under the assumption that akU is a ~ o n s t a n t .It~ is c
co
= -1 e x p [ - x (akU y)l/aI 2
+
(akUt)l/z]
erfc
[
exp[ x( -5-)"] akU erfc
-
[ x:+ 2(Dt)ll2
1
For the uranium hydride system, U varies from Uo to 0.97u0 before the intervention of spall short-circuits the diffusion path. Thus, akU may be assumed to be nearly constant. This solution has an asymptote as t m which is -+
(c/co) = exp [-x (akU / D )1 / 2 ] (3) The time scale for approaching this asymptote is of interest. The second erfc goes to zero rather quickly. The time for the first erfc to approach 99% of its asymptotic value is given by 1.65 + [1.6P + 2 ~ ( a k U / D ) ' / ~ ] ~ / ~ to.99
=
2(ak U)lI2
This is dependent on x . A value of x which is of interest is that for which the asymptotic concentration is 1%of the surface value. This is given by qO1 = 4.6/(akU/D)1/2 Back substitution into eq 4 yields t0.dx0.01)= 6.5/(akU)
(5)
(4')
For the uranium-hydrogen system, these values for ap-
proaching an asymptote are quite small. Table I gives some parameters for representative temperatures. The values of xo.ol which range from micrometers to nanometers imply that the distance that hydrogen can diffuse before being absorbed is very small except at high temperatures. The distance is small because the rate constant k is large while D is small. These small distances suggest that a one-dimensional approximation is quite adequate, even for curved or irregular samples, because the reaction zone is small compared to the radius of curvature. They also suggest that experiments which attempt to measure the reaction rates by exposing powders to a hydrogen atmosphere must use very small particles or risk having the time scale for the chemical reaction masked by the time scale for diffusion. A t the original surface ( x = 0), the concentration is c = co for all t > 0. Thus, the solution to eq l b and Id at x = 0 is
U(0,t)= Uo exp(-kcot) (6) Another time scale of interest is the incubation time, which is that time for which U(0,t)= Us. This is the time after which the spall front can move into the sample. Its value is given by t I = -1n (U,/Uo)/(kco)
(7)
Values of this time for selected temperatures assuming 1 atm of hydrogen gas pressure and Us = 0.97u0 are shown in Table I. When to,%and tI in Table I are compared, the time scale for achieving an asymptotic curve of c vs. x is much smaller than that needed to reach spall concentration. If one takes the ratio of eq 4' to eq 7 with U, = 0.97u0, the time scales have the following relationship:
The value of this ratio is small because the solubility of hydrogen in uranium is small, even at 1-atm hydrogen pressure. Note that the result that the asymptotic c curve is established much quicker than the material can deplete is solely due to the smallness of cos It is independent of either k or D.
Linear Model The spall front moves by depleting the material B, and thus the time scale for front motion is the same as that for depletion to spall concentration. The time scale for establishing the asymptotic c curve has been shown to be much shorter than that for depletion. As a consequence, the c curve ahead of the front can be estimated by using eq 3 with x replaced by x - x,. This result is crucial to the linear model. For the linear model, it is assumed that the asymptotic c curve is reestablished instantaneously whenever the spall front moves, and thus eq 3 with the substitution of x - x , for x gives the c curve for unspalled metal at all times, i.e. c/co = exp[-(x - ~ , ) ( a k U / D ) ~ / ~ l (3') A t time tI,the front begins to move into the unspalled metal pushing the c curve described by eq 3' ahead of it.
3446
Kirkpatrick
The Journal of Physical Chemistry, Vol. 85, No. 23, 1981
The exact solution to eq l b is
per unit area) corresponding to eq 6. These are (9’)
When one remembers that the c curve for t 5 t 1 is given by eq 3 and is constant in time, the U curve at tI found by evaluating the above equation is
U / U o = explln (U,/Uo) exp[-(x - ~ , ) ( a k U / D ) ’ / ~ l )(9) As the spall front moves through the material, it pushes a U curve ahead of it which is given by eq 9. That this is so can be shown by an inductive proof. Assume that eq 9 is correct at some time t and that there is some time interval At which is short enough so that the spall front velocity is nearly constant. Then, the exact solution for U may be integrated to find the change in U from t to t + At. The unknown spall velocity may be eliminated by using the fact that U(x,,t)and U(x,+u,At,t+At) both are equal to U,. The final result is that U(x,t+At) is given by eq 9 with x , replaced by x , t u,At, which is the result that one would expect from eq 9, thus proving the point. For any material whose spall front motion can be described by the advance of a plane for which U = Us, the front velocity after the incubation period is (10) where x, is itself a function of time. This result may be derived from a Taylor series expansion. For Ax and At small enough, the Taylor series for U is
au
U(x,+Ax,t+At) = U(x,,t) + dx AX
+ au At at
But, if Ax is chosen so that x,(t+At) = x , + Ax, then both U terms are equal to U, and can be cancelled. What is left is
au
- Ax ax
au At = O +at
+
e~p[-x(akU/D)’/~l)) dx
wg = a
~
E,
(-l)n+l[-(t/tI> ~ ~ In ~
n=l
-
which in the limit as Ax,At 0 is exactly eq 10. For a material with the U curve from eq 9, the front velocity is a constant given by V, = (1/t1) [ D / ( a ku)l ‘j2
+ c ~ ~ ~ e x p [ - z ( a k U / Ddx) ~ / ~ ]
The solution after considerable manipulation is
which can be rearranged to give Ax At
wg = co(akUD)l12[erf(ahUt)lj2 exp(-ak Ut)/ (aukUt)l12] (12) where wg is weight gain and wg is weight-gain rate (both are per unit area). By t = 4 / ( a k U ) , wg is within 1%of a constant value. Equations 11and 12 apply to unspalled metal during the incubation period (0 < t < tI). The units of wg are length times the units of c. The units of wg are velocity times the units of c. Depending on the units of c, it may be necessary to multiply wg and wg by a conversion factor to get the units measured in the experiment. The weight gain and the weight-gain rate given by eq 11and 12 are valid for any system which can be modeled by the perfusive precipitation equations la-e during the period before depletion of U becomes significant. These equations specifically include the phase during which the asymptotic c curve (eq 3) is being established. Thus, they are more general than the linear model which has as a principal assumption that the asymptotic c curve is established instantaneously. The weight gain and the rate found by applying this assumption are not the same as eq 11and 12. The weight gain is composed of the material A fixed in the AB compound by the reaction plus the material dissolved in the lattice. The gain is found by integrating the terms Uo - U and c over x , where U is found by applying eq 3 to eq 9’ and c is given by eq 3. The integrals are
(10’)
Equation lo’, which says that the front velocity is independent of x and t, is the heart of the linear model-the result that makes it “linear”. Its principal limits are that it requires a system which will assume the asymptotic c curve (eq 3) on a time scale that is small compared to tI and for which U,/ Uo= 1.0. Because Us/Uo= 1.0 for eq 2 and 3 to be valid, it matters little whether U, or lJo is used for U in eq 8-10 and subsequently.
Weight Gain and Weight-Gain Rate The quantity that is measured experimentally is the increase in weight of the sample including the spalled, partially reacted layers. The incubation period is not a time of zero weight gain but, because the linear model requires both that co/aU be small and that U,/Uo = 1.0, the weight gain before t~will be small compared to that afterward. Reference 5 gives expressions for the surface flux of A and time-integrated surface flux (i-e., weight gain
( ~ , t/ ~ ~ ) ~l n /
(nn!)+ coVstI (11’)
The rate is wg = aUoV,(l - U,/U,J
(12’)
The values given by eq 11’and 12’ approach those from eq 11and 12 as time grows. The values are almost identical by t = tI. After tI,the weight gain is composed of the terms from eq 11’ (evaluated a t tI) plus the material A fixed in the already spalled metal which continues to react plus that dissolved in the spalled metal. The last of these is simple as the concentration in spalled metal is assumed to be cW The value of U at any time after the spall front has passed is U(x,t)= U, exp[-kco(t - t,)] where t , is the spall time for that value of x and is given by x = V,(t, - t I ) The weight-gain integral is wg = (eq 11’ a t t I )
+ ~ U ~ & ” -{ I(u,/v,> exp[-kco(t t~ - x/V,)]) dx + coV,(t - t ~ )
The Journal of Physical Chemistry, Vol. 85, No. 23, 198 1 3447
Diffusion with Chemical Reaction
After integration and some manipulation, the gain is
WE
=
n=l
wg = aV8U,[exp(kcoT/Vs) - l]e-kco(t-tr) The rate is wg = aUoVs[l - (U,/Uo)e-kcO(t-tI)] + cOV,
(14)
At large values of time, this rate becomes a constant. Experimental results for uranium hydride show the same behavior. In principle, one may estimate the effect of phenomena which inhibit the early time reaction by extrapolating the linear portion of the weight-gain curve back to zero gain and looking at the difference between this intercept and zero. This graphically determined time interval is sometimes referred to as "induction time". There is a component of this time which comes directly from the dynamics of the system modeled here. This component of induction time will be given the symbol ti,+ The time tid is found by going to the linear portion of the weight-gain curve (t sufficiently large that the exponential in eq 13 and 14 is negligible) and extrapolating back to zero weight gain pretending that the gain was created by a constant rate equal to the asymptote of eq 14 but starting at time tid. Putting the gain from the pretended constant rate on the left-hand side and the actual gain from eq 13 on the right gives (t - tid)(aUoVs + COVs) =
5
+
a ~ ~ ~ , (-1)n+l[-ln ( t ~ (~,/~~)]n/(nn!) n=l
( t - td
-
uJ(kcouo)l
+ coVJ
Solving for tid gives tid = [us/(kcfJuO) + tI(1 -
2 (--l)"+'[-ln ( U s / U o ) ] ~ / ( n n ! ) ~ ] u U o+/ ( co) u U o(15)
n=l
The induction time as given by eq 15 as well as the incubation time from eq 7 are strictly determined by the dynamics of the system described in eq la-e. Experimental results may well show larger times due to the phenomena which may inhibit the early time reaction. Another quantity that is sometimes of interest is the derivative of the logarithm of product with respect to distance a t the spall front. The product concentration is P= Uo-u
(18)
Equation 18 shows the rate going to zero for large time as the remaining B in the fractured powder is converted to the AB compound. There is a discontinuity in both weight gain and rate at t = tf. The discontinuity in rate is easy enough to explain. At time tf, the spall front is no longer moving through the sample and, thus, the term which represents the continued increase in the thickness available for reaction must drop out. The discontinuity in weight gain is an artifact of the approximate nature of the analysis. Evaluating eq 13 and 17 at t = tf, equating the results, and canceling terms, one finds that the first term of eq 13 (the infinite series) and part of the last term are left over. These terms represent the contents of the region ahead of the spall front as given by eq 3 and 9. A t t = t f , this region is no longer there because it has met and been absorbed by the front coming the other way, but eq 13 does not reflect this. Physically, of course, there should be a smooth transition as the s p d fronts approach one another.
Comparison with Experiment Figure 1 shows a comparison between the linear model and data for the uranium hydride system. The data points and the Wicke and Otto' fit on this plot were copied directly from Figure 2 of ref 6. Thus, the figure represents a comparison of my calculations with the plot of experimental data which was published in ref 6. The two calculation lines represent spall front velocity from eq 10' for spall at 3% depletion (U,/Uo = 0.97) and 1.4% depletion (Us/ Uo= 0.986). These lines use fits for k,cotand D from ref 6 and assume a 1-atm hydrogen gas pressure. The quality of the match to the data is rather good up to -250 "C. In discussions with the author of ref 6 concerning reasons for such discrepancies as there are, the point was made that the fits for co and D as functions of temperature are based on data taken at over 300 "C and are not too reliable for lower temperatures. The model diverges from the data for temperatures over -250 " C because the hydride decomposition reaction becomes important for temperatures above this level.
Consequences of Relaxing the Assumptions The perfusive precipitation equations should be applicable to a number of systems. The linear model is restricted to the subset whose characteristics meet the linear assumptions. Some of the intermediate results from the linear model can be used for guidance about the behavior Equation 9 gives the U ahead of the spall front. The of the larger class of systems which can be described by desired derivative is perfusive precipitation. In particular, the space and time scales and the c and U curves for such systems are likely d In (P/Uo)/dx(,,,a = to resemble the linear-model results. Perhaps some ex[ ~ , / (-UUs)l O In ( W U O ) ( ~ ~ U / (16) D ) ~ / ~tensions of the linear model might apply to systems which do not satisfy all of the assumptions. The consequences Given an experimental plot of the product, one can estiof relaxing some of the linear-model approximations merit mate the value of a k U / D a t the spall front by using this study. One of the most important assumptions is that equation. U,/Uo 1.0, which makes akU approximately constant The last phase of the weight gain occurs after all of the for unspalled material. solid B metal has been spalled. The weight gain during Equation 2 is not a solution to eq l a if U is a function this phase can be calculated for a slab of thickness 2T. of x and/or t. Nevertheless, one would hope that there This phase begins at a time r_.
t f = tI
+ T/Vs
The weight gain and the rate for t > tf are
~
~~
(7) Wicke, E.;Otto, K1.Z. Phy. Chern. (Frankfurt am Main) 1962,31, 222.
3440
400
The Journal of Physical Chemistry, Vol. 85,
Klrkpatrick
No. 23, 198 1
the local c/c@ Because U ahead of the spall front had not been depleted as much as is predicted in the theory, the first term of eq 13 was too large so that the overall weight gain was less than that given by eq 13, even with the “correct” value of V,.
t
.
..
e.
2 I
12
I
I
14
16
I 18
I 20
1 I 22 24 llT x lo3 (K-’I
I 26
I 28
I 30
I 32
34
Figure 1. Comparison of model calculations to experiment. The plot of experlmentaldata, together with the line representin the Wicke and Otto7 fit, was copied directly from a plot by Condon.
8
would be an asymptotic c curve resembling eq 3 for systems in which U varies significantly. If this were so, then there should be a U curve ahead of the spall front which might resemble eq 9. If this were so, then the spall velocity would be constant, the weight gain would resemble eq 13, and the rate weight-gain would resemble eq 14. I have a computer program available which solves the set of equations la-e by using finite difference methods. This program is descended from the one used in ref 6 and solves the one-dimensional perfusive precipitation equations using the same numerical algorithm which was discussed in ref 6. Some sample cases were run on this program for uranium hydride at 80 “C with a hydrogen gas pressure of 1torr. The cases studied values of U,/Uo down to 0.25. The calculations showed a constant value for V, and a constant asymptote for weight-gain rate at large time. Referring back to eq 10,the constancy of V, implies that dU/dxlx, is also constant. The calculated value of Vs for U,/Uo = 0.25 was -60% of the value given by using U, for U in eq 10’ but 120% of the value given by using Uofor U. The value of U given by solving eq 16 for U was 0.70U0, which gave a value of V, when substituted into eq lo’, which was within 3% of that calculated. The asymptotic weight-gain rate was almost exactly that given by eq 14 provided that the calculated value of V, was used in that equation. The c curve ahead of the spall front was compared to a curve using eq 3 with the local value of U rather than some constant value. The values of c tended to run low compared to this “ideal” curve. As a consequence, the U curve values were higher than the values given by eq 9 using local U in the term akU. The U curve values were also higher than those which would be given by eq 9 if the e x p [ - x ( ~ k U / D ) ~ /term ~ ] were replaced by
Conclusions The proposed linear model gives expressions for sample weight gain and weight-gain rate through three successive stages. The model predicts an initial time period before the sample begins to spall, followed by a period of constant front velocity during which the weight-gain rate becomes constant, followed by a terminal period for which all of the metal has been spalled and during which the weight-gain rate falls to zero as the remaining B material in the spalled powder converts to the AB compound. For temperatures at which the processes considered in the perfusive precipitation equations are the only ones of importance, the linear model matches experimental results for uranium hydride to a precision which is no worse than that to which the physical constants of the system are known. Acknowledgment. This research was supported by funds from the Oak Ridge Y-12 Plant Development Division. The Y-12 Plant is operated by Union Carbide Corp., Nuclear Division, under Contract No. W-7405eng-26 with the US.Department of Energy. The work was sponsored by Dr. J. B. Condon, with whom I have had many fruitful hours of conversation on this subject. List of Symbols stoichiometric ratio of corrosive material A to material B being corroded in AB compound concentration of corrosive material A C concentration of A at surface or spall front diffuCO sivity of A in solid B D diffusivity of A in solid B k rate constant for A + B reaction P product concentration ( Uo- V) t time time at which all solid B has been spalled tf induction time (eq 15) tid incubation time (eq 7) tI time needed to reach 99% of asymptotic c curve (eq 4 and 4’) half the thickness of a slab sample T concentration of material B being corroded U value of U for which spall occurs initial value of U at zero time velocity of spall front (eq 10 and 10’) weight gain (eq 11, 13, and 17) weight-gain rate (eq 12, 14, and 18) distance from surface location of spall front distance at which asymptotic c curve has dropped to 1%of surface value (eq 5) a
