1. Introduction

It is well known that the phenomenon of heat conduction with the phase change due to melting-freezing occurs in many areas of current practical interest, e.g. the solidification of castings, fluid flow in porous media, shock waves in gas dynamics, cracks in solid mechanics, the preservation of foodstuffs, the penetration of frost into the earth, the ablation of space missiles due aerodynamic heating, statistical decision theory, cryosurgery, meteorology, astrophysics and plasma physics. As the space domain where the heat conduction equation is to be solved changes with time, these are known as moving boundary or Stefan problems (^{Carslaw and Jaeger, 1959}; ^{Mori and Akari, 1976}; ^{Ozisik, 1994}).

On the other hand, the phenomena of heat conduction and solid state diffusion are physically and mathematically similar as long as both phenomena occur due the existence of temperature and concentration gradients, respectively. In this way, these phenomena present a certain correspondence among their parameters, variables, mathematical equations and boundary conditions normally assumed. This likeness was first verified by Fick in 1855, taking into account the heat transfer analysis developed by Fourier in 1822 (^{Crank, 1975}). ^{Santos and Garcia (1998)} have presented an analytical heat transfer model which permits the characterization of the thermal behavior during cylindrical or spherical solidification. The analytical method, whose performance was evaluated by comparing theoretical predictions with experimental results, is based on geometric correction functions and variable changes that are introduced into the differential equation which represents the Fourier's law of heat conduction. The modified differential equation is solved using a closed form analytical solution, written in terms of error function, when equations are obtained that are able to describe the kinetics of solidification as well as the temperature distribution at any moment.

In this way, regarding Santos and Garcia's work, the main purpose of this paper is to propose an analytical method able to estimate times and positions of the moving boundary, as long as it has a solute concentration profile during the solid state diffusion phenomenon in a one-dimensional two-phase Stefan problem in a cylindrical radial system.

2. Analytical method based on change of variables

The problem of inward solid state diffusion in an infinitely long circular cylinder with moving boundary may be expressed by:

The boundary conditions are assumed according to the reference system adopted in Figure 1, that is:

Where *C (r,t)* denotes the concentration of one of the present phases at a radial distance *r* from the center at a time *t , C _{0}* , and

*C*

_{ξ}are the initial and moving interface solute concentrations respectively;

*r*

_{ξ}is the position of the moving interface;

*r*is the cylinder radius and

_{0}*D*is the diffusion coefficient, which is assumed to be constant.

In the development of the analytical method for the analysis of the solid state diffusion in a cylindrical radial system with phase change, med the following assumptions will be assure:

The solute atomic flux is directional cylindrical radial.

The diffusion front is macroscopically cylindrical radial.

The atomic diffusion coefficient in both phases does not depend on the solute concentration.

The solute concentration remains constant at the material surface.

The solute concentration remains constant at the moving boundary.

No contact resistance is assumed at the interface between the solute atomic flux and material surface.

In order to obtain an analytical solution for cylindrical radial geometries, as commented before, a geometric correction function must be determined and changes of variables must be performed taking into account the curvature of this radial system in such manner that when introduced in the differential equation that represents Fick's second law, make possible its application in the cylindrical radial system. The application in this modified equation of a well-known closed form solution, written in terms of error function, will allow to obtain a set of equations that will describe the displacement of the diffusion interface and the evolution of the solute concentration profile in binary systems, as shown in Figure 1.

In the work developed by ^{Santos and Garcia (1998)}, the authors have shown that times in both systems investigated must be plotted as a function of a one-dimensional relationship, which takes into account in each position, the reduction of heat transfer area in the moving interface, when it is compared with the heat exchange area of a maximum radius, which remains constant during the process. The most suitable relationship was that provided by the quotient between the volume "*V*" of material in each position and the initial area (of maximum radius) "*A _{0}*". Then we may write:

for cylinders:

for slabs:

According to the experimental results obtained, Santos and Garcia have observed that the solidification time in the cylindrical configuration was, for the same value of ( *V*_{ξ}*/A _{0}* ), approximately, two times greater than that verified in the slab bodies at the end of the process. On the other hand, when they are compared to a cylinder of radius

*r*with a plate of length

_{0}*r*using the Equation 7 for each case, for the maximum values of (

_{0}*V*

_{ξ}

*/A*) at the end of the solidification, the follow relationship is derived:

_{0}

Thus, if a geometric correction function has to cover the entire range of displacement of the solid/liquid interface toward the center of a cylinder, other differences of geometric behavior during solidification have to be analyzed, taking as reference the standard solidification of the equivalent slab (^{Santos and Garcia, 1998}). The ratio between volume of liquid and total volume of the body is certainly a variable that differs strongly, and is given by:

for slabs:

for cylinders:

The geometric correction function must provide values ranging from one at the beginning of the solidification, up to 2 at the end of the process. Then, an accurate analysis suggests the following function:

The time modified by geometric function is:

On the other hand, the concentration variation is described by the transient mass diffusion equation (Fick's second law) in the form:

Changing the variables *x* e *t* e by the variables *R* e τ, defined respectively by the Equations 7, 14 and 15, a general equation is obtained to describe solid state diffusion in the considered coordinate system:

The general solution of Equation 16 is given by:

in which *B _{1}* and

*B*are constants.

_{2}Substituting the initial and boundary conditions in the Equation 17 the follow equation is obtained:

Thus, the equation which describes the solute concentration profiles in the cylindrical radial system is given by:

The diffusion times are then described by the following expression:

Finally, the Equation 21 is defined by Equations 6 and 18:

that is the transcendental equation necessary to provide the value of λ.

3. Numerical method

The equations that describe the solid-state directional diffusion phenomenon in slabs and cylinders were numerically solved by the application of the modified variable time-step method, also called MVTS method (^{Gupta and Kumar, 1980}; ^{Gupta and Kumar, 1983}). This algorithm uses a finite difference method for the equation discretization and the Stefan condition, represented by Equation 6, for the time interpolation. The MVTS method is an extension of the method of Douglas and Gallie (^{Douglas and Gallie, 1955}), called EDG method for one-dimensional problems. The concentration variation for the case for a slab is described by Fick's second law (Equation 15) as well as by equation bellow:

Applying in the Equation 15 and 22 a whole implicit central discretization scheme, yields:

where,

For cylinders, in order to avoid instability, Gupta and Kumar (^{Gupta and Kumar, 1983}) have made the following change of variables:

The initial and boundary conditions (2), (3) and (4), become:

Substituting the new variables found in the Equations 1 and 6, the result is:

and similarly for the Equation 33:

4. Results and discussion

Theoretical positions of diffusion front from the cylinder surface until its center as a function of time, obtained by the application of Equation 20, in a cylindrical radial system under different initial conditions are shown in Figure 2 compared with the simulated results furnished by the numerical model MVTS. In both cases, a good agreement has been observed between the analytical results and numerical values.

In Figure 3 are shown the final concentration profiles as a function of position in the same cylindrical radial systems through the application of both the proposed analytical solution, Equation 19, and the numerical method MVTS.

In the same way, a very good agreement between the analytical and numerical results can be verified. In spite of the agreement observed between the values above, enough to validate them, a set of experimental results would be very important. As a practical application example of the analytical method developed, we may regard the thermochemical treatment of the carburizing of a cylindrical part, occurring at 950º C, where the phases I and II are ferrite and austenite, respectively, in the presence of carburizing gas, that is, *CO + CO _{2}* or

*CH*.

_{4}+ H_{2}5. Conclusions

The proposed analytical method, based on geometric correction functions and variables changes that are introduced in the differential equation which represents Fick's second law of solid state diffusion, was capable of satisfactorily predicting the inward solid state diffusion with phase change in cylindrical parts. It had its performance evaluated by a comparison with simulated predictions furnished by a numerical model, and a good agreement has been observed between analytical results and numerical values. Thus, the proposed theoretical equations were capable of quite approximately representing both the position of diffusion interface as a function of time and the final concentration profiles as a function of position in a one-dimensional two phases Stefan problem.