Abstract
Analytical solutions for the transient single-phase and two-phase inward solid-state diffusion and solidification applied to the radial cylindrical and spherical geometries are proposed. Both solutions are developed from the differential equation that treats these phenomena in Cartesian coordinates, which are modified by geometric correlations and suitable changes of variables to achieve closed-form solutions. The modified differential equations are solved by using two well-known closed-form solutions based on the error function, and then equations are obtained to analyze the diffusion interface position as a function of time and position in cylinders and spheres. These exact correlations are inserted into the closed-form solutions for the slab and are used to update the roots for each radial position of the moving boundary interface. The predictions provided by the proposed analytical solutions are validated against the results of a numerical approach. Finally, a comparative study of diffusion in slabs, cylinders, and spheres is also presented for single-phase and two-phase solid-state diffusion and solidification, which shows the importance of the effects imposed by the radial cylindrical and spherical curvatures with respect to the Cartesian coordinate system in the process kinetics. The analytical model is proved to be general, as far as, a semi-infinite solution for diffusion problems with phase change exists in the form of the error function, which enables exact closed-form analytical solutions to be derived, by updating the root at each radial position the moving boundary interface.
1. Introduction
Diffusion is a process which leads to an equalization of concentrations within a single phase. The laws of diffusion connect the rate of flow of the diffusing substance with the concentration gradient responsible for this flow [1]. The need to study this physical phenomenon may be justified by the great influence it has on the solute enrichment rate and concentration profiles of this element in the solid, therefore, being important in determining the structure and physical, mechanical, and metallurgical properties of obtained products. The diffusion rate, generally in a transient regime, becomes the mathematical analysis of this phenomenon difficult since it leads to differential equations to present nonlinear boundary conditions making it difficult to obtain exact analytical solutions [2] requiring the establishment of physical and/or mathematics simplifying hypotheses from the real conditions so that such solutions may be made viable. A limited number of exact, closed-form solutions to the multiphase, diffusion-controlled, and moving-interface problem exists in the literature. The exact, closed-form solutions generally consider the outermost phases to be infinite in extent as well as assume an interfacial diffusion rate function. It is through use of such boundary conditions that closed-form solutions are readily obtained.
The analytical methods proposed for the study of solid-state diffusion are generally based on infinite series which may offer some difficulty in the face of practical applications since their mathematical resolutions normally occur through nontrivial processes. These methods are mostly limited to studying diffusion in binary systems with planar diffusion geometry, so despite the importance of the subject, exact analytical solutions for cylindrical and spherical geometries have not yet been obtained since they present greater mathematical difficulties arising from the complexity of the mathematical equations as well as from the assumed boundary conditions. Although the numerical methods require a new reprocessing each time a parameter is introduced or modified and need the use of computational resources with the programming involved the more complex the greater the degree of precision desired, these methods lead in the most of the time to a greater proximity of cases observed in practice allowing, for example, more real boundary conditions to be admitted. The numerical methods are the ones that have been used most extensively for the study of solid-state diffusion in systems with cylindrical and spherical geometries. Numerical solutions demonstrate sometimes a certain complexity; however, they present as one of their fundamental advantages the fact that they represent problems of practical interest for which it would not be possible to obtain analytical solutions. It is well known that the conductive heat transfer process [2–5] and solid-state diffusion are physically similar since both occur due to the existence of a temperature gradient and a concentration gradient, respectively, and therefore show a certain correspondence between their physical parameters, variables, mathematical equations, and boundary conditions normally assumed. Thus, many of analytical and numerical solutions proposed for the mathematically analogous solidification problem are directly applicable to diffusion problems.
The concept behind the term “moving boundary problems,” according to Crank [6], is associated with problems whose solutions of differential equations must satisfy certain boundary conditions of a corresponding domain, referred as boundary value problems. Despite this, in many cases of interest, the domain boundary is not known beforehand as the process develops but needs to be determined as a part of the solution. The term “moving boundary problems” is normally applied when the domain boundary is stationary, and a steady-state problem exists. On the contrary, the terminology “free boundary problems” is related to the time-dependent problems and the boundary needs to be determined as a function of time and space. The terms “free” and “moving” are carried on in different ways. Some authors prefer to encompass both kinds of problems in a singular term, i.e., “free boundary problems,” and few ones used to apply “moving boundary” as a more general case. The moving boundary problems require a solution of a parabolic-type differential equation. The practical applications are fluid flow in porous medium, diffusion, and heat conduction with chemical reactions or phase change and, additionally, to help to understand many other free and moving boundary problems, such as collapse of dams, shock waves in gas dynamics, and cracks in solid mechanics.
Jost [1] and Crank [7] provided a good literature survey concerning free and moving boundary problems as well as several studies dealing with phase-change problems. In this sense, an important change of variable technique was presented by Tien and Churchill [8] for the solidification of cylinders, which can be applied to minimize the nonlinearity arisen from the decreasing of the cylindrical radius as the phase change gets close to the center of the cylinder. Tanzilli and Heckel [9] developed a numerical method for treating the two-phase, diffusion-controlled, and moving-interface homogenization problem involving finite boundaries. Solutions were obtained for planar, cylindrical, and spherical diffusion geometries. The effects of parameters such as the interdiffusion coefficients of both phases, the concentrations of each phase at the interface, mean composition, initial size of the dissolving phase, and initial composition of the two phases were considered. The results are displayed to show both the position of the interface and the degree of homogenization as a function of diffusion time. Typical examples to which the results can be applied are also discussed. The numerical solutions devised by Tenney and Unnam [10] for two-phase concentration-dependent diffusion coefficients were calculated for planar, cylindrical, and spherical geometries to compare the effect of interface geometries with those caused by concentration-dependent diffusion coefficients, and two methods of averaging D were considered to determine the best averaging method for different types of D-variations. The effects of interface-location criteria on mass conservation and convergence of interface location, diffusion coefficient variation in the alpha and beta phases of a two-phase binary alloy system, effect of D variation in a cylindrical couple on beta-phase thickness, and geometry and D-variation effects on the degree of homogenization were also determined. The mathematical model proposed by Bongartz et al. [11] permits the prediction of carbon concentration profiles in carburized high-temperature alloys. It is assumed that a proportion of the carbon which diffuses into the alloy reacts with elements such as chromium to form carbide precipitates. The amount of carbon remaining in solution is determined from the solubility product of the carbide and only this carbon in solution is able to diffuse through the alloy matrix, and thus the carbide precipitation process reduces the rate of carburization. The model can be used to adapt carbon concentration profiles from one geometrical configuration (planar, cylindrical, and spherical and diffusion geometries) to another. Unnam [12] developed an analytical method which shows that the variation in solute concentrations during diffusion into parts with both cylindrical and spherical geometries, with radii from a certain dimension, can be represented with relative precision by analytical solutions developed for the planar geometry. The author assumes that the diffusion coefficient does not depend on the solute concentration. The report presented by Agren [13] describes a numerical method to treat diffusional reactions in alloys where a thermodynamic model, complete enough to represent phase equilibria, as well as diffusion coefficients are available. The method combines a rigorous treatment of the thermodynamic equilibrium locally at a moving phase interface with an ambitious treatment of the diffusion inside the one-phase regions and contains no formal restrictions on the number of components. Different kinds of phases and substitutional and interstitial solutions and phases, with stoichiometry constraints, e.g., carbides and nitrides, can be handled. Castleman and Wo [14] devised a numerical solution for the set of differential equations characterizing interdiffusion in a ternary system in which the partial diffusion coefficients are linear functions of the concentrations. The concentration-penetration curves obtained accurately reflect the behavior of an interdiffusing ternary ideal solid solution provided that the tracer diffusion coefficients are constant and the ratios of the latter are not very different from unity. It is shown that, in an ideal interdiffusing ternary solid solution, zero-flux plane positions cannot correspond exactly to the points of intersection of the diffusion path and iso-activity lines drawn through the terminal alloys of the diffusion couple, but approach the latter in the limiting case. The numerical method proposed by Yang et al. [15] for studying the diffusion of hydrogen in metals takes into account the retention effects caused by imperfections existing in their crystalline structure. The diffusion coefficient depends on the hydrogen concentration and considers both the hydrogen that diffuses in the atomic form and retains in the molecular form by various types of defects such as dislocations, microcracks, grain boundaries, and microporosities. The analytical solutions developed by Moreira [16] to represent the kinetics and the solute concentration during transient single-phase solid-state diffusion in slabs, cylinders, and spheres are based on the geometric correlations and variables changes in order to represent Fick’s second law. The resulting modified equation is solved by utilizing a well-known analytical solution for semi-infinite slabs, based on the error function. The solution permits to analyze diffusion times and solute concentration profiles in slabs and during inward solid-state diffusion in cylinders and spheres. Swaminathan and Voller [17] presented two fixed-grid enthalpy methods applied to phase change problems regarding the apparent heat capacity and source-based methods in order to derive a general enthalpy method that encompasses both techniques. Zerroukat and Chatwin [18] conducted comparisons among several numerical solution schemes for moving boundary problems based on the finite different method, such as modified variable time step method and enthalpy formulation techniques, and discretization strategies and explicit and implicit discretization for computation of single and multiple moving boundaries, and Voller [19] derived a similarity solution for solidification of an alloy, regarding a complete coupling of solute and temperature fields in the mushy region, eutectic reaction, macrosegregation, microsegregation, and fluid flow in Cartesian coordinates. The analytical heat transfer model devised by Santos and Garcia [20] 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 the geometric correlations and variables changes that are introduced in the differential equation which represents Fourier’s law of heat conduction. The modified differential equation is solved using a closed-form analytical solution, written in terms of error function, where obtained equations are able to describe the kinetics of solidification as well as the temperature distribution during the solidification process. Ferreira [21] proposed an exact solution for the inward two-phase solid-state diffusion in cylinders based on an exact solution for slabs which was corrected by a geometric correlation in order to provide the updates of roots for each radial position. The obtained root equation by the substitution in the solution based on the error function and the similarity variable into the moving boundary equation behave like a transcendental equation that must be solved for each radial position of the moving interface. The study conducted by Ferreira et al. [22] aimed to investigate the effect of alloy solute content, melt superheat, and metal-mold heat transfer on inverse segregation during upward solidification of Al-Cu alloys. The experimental segregation profiles of analyzed alloys were compared with theoretical predictions furnished by analytical and numerical models. The analytical model was based on an analytical heat-transfer model coupled with the classical local solute redistribution equation proposed by Flemings and Nereo [23], while the numerical model was proposed by Voller [24] with some changes introduced to take into account different thermophysical properties for the liquid and solid phases, time variable metal-mold interface heat-transfer coefficient, and a variable space grid to assure the accuracy of results without raising the number of nodes. It was observed that both numerical and analytical predictions conform with the experimental segregation measurements. Under the assumptions of zero surface energy, no crystal anisotropy, and infinite mobility, Voller [25] developed an analytical solution for the solidification of a solid seed in an undercooled binary melt in Cartesian, radial cylindrical and spherical coordinates. Li et al. [26] adopted a fixed-grid source-based method, originally presented to simulate the temperature fields for melting-solidification phase change processes, to simulate diffusion-controlled dissolution and solidification. The method solves a unique diffusion equation for the different phases and the moving interfaces using implicit time integration. Compared with previously developed models, it is not only simpler in numerical formulation and procedure but also more convenient to extend to many phases and high-dimensional problems. The authors validate the model using experimental data and previous modeling predictions for several systems available from the existing literature. A fully implicit two-dimensional moving-mesh finite element simulation model was devised by Ghoneim and Ojo [27] to investigate the influence of grain boundaries in polycrystalline solids on diffusion-controlled liquid-solid transition during transient liquid phase (TLP) bonding. The model was found to conserve solute, and its calculated solutions were unconditionally stable and in good agreement with experimental results. Numerical calculations and experimental verification showed that enhanced intergranular diffusivity had a minimal effect on the time required to achieve complete diffusion-induced solidification in cast superalloys. Diffusional mass transport, in which Fick’s laws of diffusion can be used [28], is almost always included in the control of drug release out of a dosage form [29–31]. In several cases, drug diffusion is the predominant step [32–34]; however, it also plays an important role in combination with polymer degradation/matrix erosion [35–37] or polymer swelling [38–40]. The boundary conditions concern the conditions for diffusion at the boundaries of the drug delivery system. If the device dimensions are constant with time (no significant swelling or dissolution/erosion), the boundaries are called stationary. In contrast, in case of time-dependent device dimensions, the boundary conditions are called “moving.” If the device swells significantly, the boundaries move outwards and if the system dissolves/erodes significantly, they move inwards. In this sense, Siepmann and Siepmann [41] gave an overview on the current state of the art of mathematical modeling drug release from delivery systems, which are predominantly controlled by diffusional mass transport. The respective equations for a broad range of devices were indicated, including reservoir and matrix systems, exhibiting or not an initial excess of drug and the geometry of slabs, spheres, and cylinders. The assumptions the models were based on as well as their limitations were pointed out. Practical examples illustrate the usefulness of mathematical modeling of diffusion-controlled drug delivery. More recently, Castillo and Morin [42] have performed a mathematical model for drug dissolution-diffusion in nonerodible nor swellable devices in order to analyze existence, uniqueness, and regularity properties of the system. The authors have obtained a coupled nonlinear system which contains a parabolic equation for the dissolved drug and an ordinary differential equation for the solid drug, which is assumed to be distributed in the whole domain into microspheres which can differ in size. Simulations illustrating some features of the solutions and a good agreement with laboratory experiments are presented. The work conducted by Mishra et al. [43] has analyzed the hyperbolic heat conduction in a 1D planar, cylindrical, and spherical geometry using the lattice Boltzamnn method (LBM). Temporal temperature distributions were analyzed for thermal perturbation of a boundary by suddenly raising its temperature as well as also by imposing a constant heat flux to it. Wave-like temperature distributions in the medium were obtained when constant temperature boundary condition was used. To check the accuracy of the LBM results, the problems were also solved using the finite difference method (FDM). LBM and FDM results were compared exceedingly well. The main objective of the numerical model proposed by Olaye and Ojo [44] was to investigate diffusion-controlled interphase-interface migration kinetics. The model is a fully explicit numerical method, incorporates constant and variable diffusion coefficients, and is more accurate, direct, and faster than implicit models because it does not involve nontrivial assumptions that decrease the accuracy of implicit models. The authors use their model to modify the classical Heckel’s criteria for second-phase growth, to predict factors that affect the growth rate, homogenization time, and extent of second-phase growth, in a two-phase alloy. The study performed by Lin and Tsay [45] has constructed a more generalized dissolution-diffusion model for drug release from a spherical matrix system loaded with dispersed undissolved drug particles of different shapes. The effect of the geometrical shape of the dispersed drug particles on the drug release profile was examined for systems with different ratios of dissolution rate/diffusion rate and ratios of solubility/total drug loading value. Finally, Ghanbar and Ojo [46] have presented a new finite difference numerical model in order to investigate the importance of cylindrical and spherical solid-liquid interface geometry on the kinetics of isothermal solidification during diffusion brazing or transient liquid phase repair process. The obtained results showed that nonplanarity of the solid-liquid interface can have significant effects on the kinetics of solid-liquid phase transformation during the repair process which is crucial to the understanding and optimization of this process in nonplanar systems.
In the preceding paragraphs, several special cases were treated by mathematical methods; however, these methods are not sufficient for all cases of practical importance. Therefore, based on the geometric correlation proposed by Moreira [16] as well as on the exact solution developed by Wagner (unpublished work) reported by Furzeland [47] for a semi-infinite slab, in this study is presented an exact solution which is able to estimate diffusion times and position of the moving boundary, as well as the solute concentration profiles during inward transient solid-state diffusion in radial cylindrical and spherical geometries.
The current approach represents a method for a set of solutions for the transient diffusion (heat/mass) with phase change (reaction) for cylindrical and spherical coordinate systems, independently of the number of moving boundaries. There is no other analytical closed-form solution in the literature to deal with such problems. The proposed method remained two decades unpublished, as the numerical solution effort to solve properly the initial solid-state diffusion problem with such a high degree of accuracy imposed by the radial cylindrical coordinate system was considerably difficult to achieve due to nonlinearities arisen in the moving boundary. The numerical resolution at that time allowed us to conclude that the proposed analytical solution for the cylinder was approximate. And only recently, by chance, the authors decided to numerically solve the diffusion times and concentration profile, and the comparison with the analytical results have proved to be a close-form analytical solution for cylinder, which in this paper is extended to deal with the sphere for single- and two-phase problems. During the time of the Master of Science Thesis [21], only an approximate numerical solution was obtained by applying the Modified Variable Time Step Method (MVTS) in a cylindrical coordinate system. In this paper, a variation of the enthalpy method proposed by Swaminathan et al. [17] was applied, and a more accurate numerical result was obtained allowing to demonstrate that both analytical solutions consist of closed-form solutions. Today, only approximate solutions for cylindrical and spherical phase-change problems are available. The technological appeal of energy-storage and pharmaceutical applications are the current motivations. The simplicity of both proposed solutions and their main feature as low-computational demanding resources, allow engineers and scientists to implement them directly in the microcontrollers or low-profile single-board computers. That is not the case for the numerical schemes found in the literature to deal with moving boundary problems.
2. Mathematical Formulation
The mathematical formulation here introduced is found in Jost [1] for the solid-state diffusion and Crank [7] solidification of pure metals. The solution to be derived aims to analyze the changes promoted by the transient solid-state diffusion in cylindrical and spherical radial geometries and imposed by the effects of the evolving curvature, as compared to the atomic flux of Cartesian coordinates. In developing the analytical solutions, the following assumptions are taken:(a)The atomic flux is unidirectional Cartesian, radial cylindrical, or radial spherical(b)The diffusion front is macroscopically plane, cylindrical, or spherical(c)The diffusion coefficient or thermal diffusivity is concentration independent(d)The solute concentration remains constant at the surface of the semi-infinite slab, cylinder, or sphere(e)The solute concentration remains constant at the moving boundary(f)The considered systems are isotropic(g)The contact resistance in the interface between the atomic flux of solute and the material surface is neglected
2.1. Single-Phase Solid-State Diffusion
According to the assumptions assumed, the transient solid-state diffusion for unidirectional Cartesian [1] and radial cylindrical and radial spherical geometries can be represented by Fick’s second law, for single-phase solid-state diffusion with phase change, as shown in Figure 1, and the governing equation and boundary conditions can be written as follows:for ,for ,where is the initial concentration, is the surface concentration, is the solute concentration at the moving interface , is the atomic diffusion coefficient, is the spatial coordinate, is the time variable, for slab, for cylinders, and for spheres.
2.2. Single-Phase Solidification
In the case of the thermal diffusion equivalent problem [7], it is called single-phase solidification. The reference coordinate for this problem is shown in Figure 2. For this problem at and , the temperature at the moving interface is given by :for ,for ,where is the solid phase density, is the latent heat of fusion per volume, is the melting temperature, is the surface temperature, is the temperature at the moving interface , is the thermal diffusivity, and, finally, is the thermal conductivity of the solid phase.
2.3. Two-Phase Solid-State Diffusion
Considering now two-phase solid-state diffusion with phase change [1], the reference coordinate for this problem is shown in Figure 3:for ,for ,where is the solute concentration at the moving boundary interface , is the solute concentration at the interface , and and are diffusion coefficients of the phases I and II, respectively.
2.4. Two-Phase Solidification
In the case of two-phase solidification [7], the reference coordinate for this problem is shown in Figure 4. For this problem, at and , the temperature at the moving interface is given by . The boundary conditions will be left in the form of establishing an equivalence between the solid-state diffusion and heat conduction with phase change phenomena:for ,for ,where is the liquid phase thermal conductivity and is the initial temperature of the liquid phase.
During the solid-state diffusion or the heat conduction with phase change in a slab, the surface at the diffusion front remains constant until the end of the diffusion species/thermal is achieved. Otherwise, in the radial cylindrical and spherical systems, due the effects imposed by the geometric curvature, the surface varies as far as the process advances. Such feature influences directly the flux of species/heat. Therefore, the development of geometrical correlations and the use of a suitable change of variables will enable the use of the partial differential diffusion equation in Cartesian coordinates to represent the radial cylindrical and spherical coordinates.
2.5. One-Dimensional Relation
Santos and Garcia proposed, for the analysis of transient solidification of metals under radial cylindrical [20, 48] and spherical [20] heat flux conditions, geometrical correlations to deal with the geometric curvature. These authors stated that the solidification time must be analyzed as a function of the one-dimensional correlation that considers, for each radial position, the decrease in the surface area of heat exchange at the solid/liquid interface as compared to that of the slab, which remains constant. The correlation that has been shown to be more appropriate was the relation between the volume of metal affected by transformation in each radial position and the surface area of the metal/mold interface :
For slabs,for cylinders:and for spheres:where the variable is the radius of a cylinder and a sphere, is the position of moving interface, is the one-dimensional relation of transformed volume in relation to the initial surface area, for Cartesian, for radial cylindrical, and for radial spherical coordinate systems.
2.6. Geometric Correlation
In 1991, Moreira [16] proposed an analytical model representing the solid-state diffusion in binary single-phase systems for radial cylindrical and spherical coordinates. The representative equations were derived by the application of geometric correlations and a corresponding change in variables able to encompass the curvatures of the radial cylindrical and spherical coordinates. These geometric correlations are introduced into the solution of the partial differential for the semi-infinite slab, based on the well-known error function, in order to permit the diffusion interface position and the evolution of the solute concentration profiles during inward solid-state diffusion in cylinders and spheres to be determined. These correlations obtained by Moreira for cylinders and spheres are the following. For cylinders,and for spheres,
2.7. Change of Variables
In the radial cylindrical and spherical systems of coordinates, the quantity in the transient solid-state diffusion for Cartesian coordinates will be exchanged by , so to be applied to radial coordinates. On the contrary, the introduction of the geometric correlation in the partial differential equation, will be made through the change of variable. From to a modified variable T, which encompasses the geometric correlation:
2.8. Deriving the Analytical Solution for a Single-Phase Diffusion
Figure 5 shows a schematic representation of solid-state inward diffusion in a cylinder and a sphere.
(a)
(b)
(c)
Applying the change of variables in terms of , , and in equations (15)–(21), the mathematical formulation of the problem of solid-state diffusion can be rewritten in the following form:
For ,and for ,where
A general solution in terms of the error function, for solid-state diffusion of a semi-infinite slab, was firstly proposed by C. Wagner [47], regarding the performed change of variables, which gives
Considering equation (29) and the boundary condition given by equation (24),substituting the boundary condition of equation (25) into equation (29),
Then, rearranging equation (31) asintroducing the geometric correlation into the similarity variable and rearranging equation (33) in order to obtain an expression in terms of the diffusion time, i.e., by substituting equation (21) into equation (33),
Consequently, the diffusion time can be expressed asby inserting equation (34) into equation (31):
Finally, by combining equations (30) and (36) into equation (29), an expression for the solute concentration profile can be obtained as follows:
Equation (37) presents the transient solute concentration profile for cartesian, cylindrical, and spherical coordinate systems. The derivative of equation (37) with respect to the spatial variable “,” at the moving interface position, yields
Introducing equation (33) in equation (38),
On the contrary, using equation (33) to derive with respect to T, provides
Finally, by combining equations (39) and (40) with the equation of the moving interface, equation (27), an expression is obtained to calculate and update the roots for each radial position, i.e., each position in the one-dimensional relation :
In the case of the single-phase solidification,where is the specific heat of the solid phase.
2.9. Diffusion Times
The diffusion times for radial cylinder and radial sphere coordinates can be obtained by applying into equation (35) the one-dimensional relation of volume affected by the transformation and initial surface area, equations (17) and (18), and the geometric correlation for cylinder and sphere, equations (19) and (20), respectively:
The derived equation for the diffusion times of inward diffusion in a cylinder can be expressed as
In the case of single-phase solidification,and for the atomic diffusion in a sphere,and, finally, for single-phase solidification,
2.10. Concentration Profiles
The diffusional concentration profiles for inward diffusion in a cylinder and a sphere can be described by equation (37), and multiplying the erf function argument by the term provides
Substituting equation (15) for and , the following expression can be determined as a function of one-dimensional relations and , respectively:
By inserting equation (17) for the cylinder and equation (18) for the sphere in equation (46), the concentration profiles for cylinders and spheres can be derived asfor single-phase solidification inward a cylinder,for diffusion inward a sphere, the solute profile is given byand, for the single-phase solidification,
Now, equations (43) and (44) for diffusion times and equations (47) and (48) for the solute concentration profiles for cylinders and spheres, respectively, can be carried out by solving equation (41) for each radial position. In other words, for positions in the solution domain, which represents one-dimensional correlation positions , it implies now there will be roots to be calculated from equation (41).
2.11. Deriving the Analytical Solution for Two-Phase Diffusion
Applying the change of variables in terms of , , and in equations (7a)–(14a), the mathematical formulation of the problem of solid-state diffusion can be rewritten in the following form:andFor ,and for ,
A general solution in terms of the error function, for solid-state diffusion of a semi-infinite slab, was firstly proposed by C. Wagner [1], regarding the performed change of variables, which givesand
Substituting the boundary condition equations (53) and (54) into equations (57) and (58) asandintroducing the geometric correlation into the similarity variable ,and rearranging equation (61) in order to obtain an expression in terms of the diffusion time,and, by substituting equations (21) into (62) provides
Rearrange equation (59),and equation (60),where
Finally, combining equations (64) and (57) provides an expression for the solute concentration profile for 0 and which can be obtained as follows:and, similarly, considering equations (65) and (66) and equation (60),
Equations (67) and (68) represent the transient solute concentration profile for Cartesian, cylindrical, and spherical coordinate systems. The derivative of equations (67) and (68) with respect to the spatial variable “,” at the moving interface position, , yields
Introducing equation (61) into equation (69),
Introducing equations (61) and (66) into equation (71),
On the contrary, using equation (61) to derive with respect to T provides
Finally, by combining equation (70), equations (72) and (73) with the equation of the moving interface, and equation (56), an expression is obtained to calculate and update the roots for each radial position , i.e., each position in the one-dimensional relation :
Considering the two-phase solidification,
2.12. Diffusion Times for Two-Phase Diffusion
The diffusion times for radial cylinder and radial sphere coordinates can be obtained by applying into equation (63) the one-dimensional relation of volume affected by the transformation and initial surface area, equations (17) and (18), and the geometric correlation for cylinder and sphere, equations (19) and (20), respectively:
The derived equation for the diffusion times of inward diffusion in a cylinder can be expressed asand, for two-phase solidification,and for inward diffusion in a sphere,and, in similar manner for two-phase solidification,
2.13. Concentration Profiles for Two-Phase Diffusion
The diffusional concentration profiles for inward diffusion in a cylinder and a sphere can be described by equations (67) and (68), and multiplying the erf function argument by the term provides
Substituting equation (15) for and , the following expression can be determined as a function of one-dimensional relations and , respectively:and, similarly for equation (79),
By inserting equation (18) for the cylinder and equation (19) for the sphere in equations (80) and (81), the concentration profiles for cylinders can be derived as follows.
For ,and, for one phase-solidification for the region and, for ,and, for one phase-solidification for the region
The case of concentration profiles for spheres provides for ,and, for two-phase solidification,and, ,and, finally for two-phase solidification,
Now, the set of equations for the concentration profile, equations (82) and (83), for cylinder and, equations (84) and (85), for sphere, and diffusion times and equations (76) and (77) for the solute concentration profiles for cylinders and spheres, respectively, can be carried out by solving equation (74) for each radial position, that is, for positions in the solution domain, which represents one-dimensional correlation positions . By consequence, regarding the correlation, now there will be roots to be calculated from equation (74).
3. The Numerical Method
Aiming to verify the validity of the analytical solutions for diffusion times and concentration profiles for radial cylindrical and spherical coordinate systems, a numerical method based on the well-known enthalpy method [49] were applied to a slab, a cylinder, and a sphere. The complexity of solving the moving interface and its high nonlinearity, due to the rise of a discontinuity in the concentration field between the regions and , was circumvented by rewriting the solid-state diffusion partial differential equation in terms of a concentration density field [22], e.g.,for ,for ,where , for the slab, , for the cylinder, and , for the sphere. The diffusion coefficient, in this case, is now taken as in .
The equation to deal with the concentration density field is similar to those applied for the heat conduction with phase change [6, 18] as follows:
The concentration “” can be extracted from the concentration density field via
4. Results and Discussion
In this section, we will be presenting the solid-state diffusion with phase change and the equivalent solidification for single-phase and two-phase processes. The present model is general in nature for cylinders and spheres, regarding the diffusional radial flux of atoms/heat, as in the case of equations (1)–(14) which were considered for atomic and thermal diffusion with phase change. The solution regarding two moving boundaries can also be derived. Two-phase diffusion was chosen as an additional application example. If an analytic solution for diffusion problem exists for semi-infinite slab in terms of error function, for instance, one that considers convective boundary conditions [6], this method can be applied with no-restriction.
4.1. Single Phase
Figure 6 presents the geometric correlations derived by Moreira [16] for solid-state diffusion of slabs, cylinders, and spheres as a function of the dimensionless interface position undergone to diffusion. Ferreira [21] applied the geometric correlation “” to derive an exact analytical solution for inward diffusion in a cylinder. As can be noticed, at the end of the diffusional process, these correlations assume for slabs, for cylinders, and finally for spheres.
Aiming to verify the validity of the analytical solutions of diffusion times and concentration profiles for radial cylindrical and spherical coordinates, a numerical method based on the well-known enthalpy method [49] were applied to a slab, a cylinder, and a sphere, as shown in Figure 7 for single-phase atomic diffusion with phase change and Figure 8 for single-phase solidification. The space grid used were divided in 100 representative elementary volumes (REVs) and . A good agreement is observed for all cases between analytical results and numerical scatters. As the values for diffusion time equation (94), solute concentration, and temperature profiles, and equation (95) are the same in dimensionless form, there is no need to provide a further set of comparisons with the numerical method for the case one-phase solidification. That is the reason why the behavior of Figures 7 and 8 are the same, only the scale of the problems differs from one another:where and are the dimensionless time and concentration or temperature, respectively.
Figures 9 and 10 present analytical simulations for a slab, a cylinder, and a sphere at dimensionless interface positions of 30% of for atomic diffusion and single-phase solidification. In similar manner, Figure 11 present for atomic flux and Figure 12 present for solidification at a position of 70% of . As expected, considering Figures 8 and 11, the lowest diffusion rate is observed for the slab and the highest for the sphere, as the amount of material for an atom to be diffused decrease as the diffusion thickness of the cylinder and sphere increase. This fact is directly related to the geometric correlation factor , which according to equation (20) changes the diffusion kinetics to conform the geometric curvature to the radial cylindrical and spherical coordinate systems.
Figure 13 shows the diffusion times as a function of the one-dimensional relation representing the diffusion interface position for the three geometries. Both cylinder and sphere provide the higher times as compared to that of the slab. This observation can be easily explained in terms of the thickness of material subjected to the diffusion process. For the same value of the one-dimensional relation , the slab has the greatest amount of material remaining, followed by the cylinder and sphere, respectively.
Figure 14 provides the similarity variable equations (24)–(33) for slab, cylinder, and sphere. In the case of the slab, the root is unique for all the domains, as it does not present any geometric curvature to exert influence on the diffusion kinetics. Nevertheless, for cylinders and spheres, this is not the case. For inward diffusional process, considering both radial geometries, the curvature imposes an acceleration of the process as far as the amount of remaining material not undergone to the atomic diffusion decreases as the radius decreases. Then, for each one-dimensional relation for the radial coordinate systems, the root must be updated at each radial position.
Figure 15 shows similarity variables for the slab, cylinder, and sphere as a function of the dimensionless interface position . As can be observed, Figure 7 is very similar to Figure 6. The reason is that, in equation (41), where the roots are updated at each radial position, the only difference among the three geometries is the geometric correlation . This observation permits to state that equation (33) changes the kinetics of the diffusion process to encompass the geometric curvatures of the radial geometries. This is an important remark of the solution, as far as, once provided a one-dimensional relation and a geometric correlation and an analytical solution for a constant cross-section geometry, an exact analytical solution for a generalized coordinate geometry problem, based on these assumptions can be derived.
4.2. Two Phases
Aiming to validate the analytical method developed to obtain transient radial cylindrical and spherical diffusional flux with phase change solutions as a more general solution, once error function solution exists for a semi-infinite slab, by applying a one-dimensional correlation to encompass the kinetic of radial geometries by the geometric update of the roots provided by equation (74) for each radial position of the moving interface. By considering the proposed solution, a solid-state diffusion in the two-phase system and a two-phase solidification problem was chosen, as shown by the governing equations (7) to (14) and the solution derived from equations (49) to (85). The diffusion times for the atomic diffusion and for the solidification problems are presented in Figures 16 and 17, where it can be noticed that atomic diffusion times and solidification times have the same features, but the solid-state diffusion is a thousand times faster than solidification. While performing upward solidification of pure Al in a 100 mm height mold, Ferreira et al. [22], even considering there was a high transient global heat transfer coefficient at the bottom, not assumed here, but instead a Dirichlet (prescribed temperature) boundary condition at the heat cooling surface (chill), times around 70 seconds were observed for pure aluminum, what agrees for the times provided by the Neumann solution (semi-infinite slab) represented by the black line in Figure 17.
Figures 18 and 19 present analytical simulations for a slab, a cylinder, and a sphere at dimensionless interface positions of 30% of for two-phase atomic diffusion and two-phase solidification. On the contrary, Figures 20 and 21 represent the phenomena at a position of 70% of the moving interface. Considering that, in the solid-state diffusion case, atoms are being transported into the sample from the surface, while in the solidification problem, heat is extracted from the bulk to the surface. In the first case, the slab concentration profile is the highest, while in the second case it is the lowest.
In the similar way, Figure 22 provides the similarity variable equation (61) for the slab, cylinder, and sphere for the case of two-phase diffusion and solidification phenomena. The same behavior is noticed for the slab, as shown in Figure 14. The root is unique for all the domain, as it does not present any geometric curvature to influence the diffusion kinetics. Nevertheless, for cylinders and spheres, this is not the case. The curvatures at positions close to 0.01 concentration, for the cylindrical and spherical solutions are due to the existence of a profile in the second phase. For inward diffusional process considering both radial geometries, both the curvatures impose an acceleration of the process, as far as the amount of remaining material not reached by the moving boundary interface decreases with the radius. This implies that, for each one-dimensional relation for the radial coordinate systems, the root must be updated at each radial position of the moving interface, for equations similar to equation (74), acting as a transcendental equation.
Figure 23 shows similarity variables for the slab, cylinder, and sphere as a function of the dimensionless interface position for the two-phase diffusion. The behavior is similar to that of Figure 15 for one-phase diffusion. The kinetic implies for the radial geometries that the similarity variable is two folds that of the slab for a cylinder and three folds for a sphere at the end of the process. That is why it is normally observed that the Neumann solution is enough to be applied providing accurate predictions at positions very close the surface of cylinders and spheres, as far as the influence of geometry on the kinetic can be neglected.
5. Conclusions
Exact analytical solutions for cylinders and spheres were derived for transient one-phase and two-phase solid-state atomic diffusion and solidification regarding inward radial cylindrical and spherical geometries. The exact analytical solutions for both radial coordinate systems were derived from two solutions for a semi-infinite slab based on the error function, performing a suitable change in variables and the application of a geometric correlation factor and the corresponding one-dimensional volume-surface area relations in order to achieve a closed-form solution, in this case, for single-phase and two-phase phase-change problems. The geometric correlation changes the basic partial differential equation for space and time variables, in order modify the diffusion coefficient for each one-dimensional relation position to encompass the effect of the curvature of the radial geometries on the kinetics of the semi-infinite slab solution. The model proved to be general, as far as, a semi-infinite solution for diffusion problems with phase-change exists for a slab, the kinetic of the radial geometries enables exact analytical solutions of inward diffusion in radial geometries to be derived, by updating the root for each radial position of the moving boundary interface by equations similar to equations (41) and (74), acting exactly as a transcendental equation, providing an update for the roots, as observed in the behavior of similarity variable as a function of the solute concentration for single-phase and two-phase diffusion.
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors acknowledge the financial support provided by FAPERJ (The Scientific Research Foundation of the State of Rio de Janeiro-Brazil), CAPES (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior-Brazil-Finance Code 001), and CNPq-Brazil (National Council for Scientific and Technological Development).