#### Abstract

A model and its numerical implementation are presented which describe the formation of mixed Zn/ZnO droplets by rapid cooling of a mixture of zinc vapor and oxygen. The model incorporates a nucleation step producing pure metal droplets of a fixed size and three condensation-like processes determining the further growth of the droplets by adding metal atoms and oxygen. Properties to characterize the resulting droplet population are obtained from the model. Examples such as the number density of particles, their average chemical composition and surface area, or the particle mass distribution are presented. The model is verified on a qualitative level by comparing with experimental data the influence of the quench rate and the initial partial pressures of the reactants on the average chemical composition of the product. The model predicts, conforming qualitatively with experimental findings, that increasing the quench rate has little influence on the chemical composition of the products but that decreasing the initial partial pressures of zinc results in less oxidized products, thus, a higher average zinc content.

#### 1. Introduction

PSI’s Laboratory of Solar Technology became interested in modeling the reaction of zinc vapor with oxygen to yield a mixture of zinc and zinc oxide during its quest to convert concentrated solar radiation into storable and transportable chemicals, the so-called solar fuels. To achieve this conversion thermo-chemical cycles [1–3] are applied. At present, the Zn/ZnO-based cycle [4–7] is considered a promising candidate. In this cycle, ZnO is thermally decomposed into zinc vapor and oxygen at temperatures above 2000 K in a solar thermochemical reactor [8, 9]. After separation from oxygen, zinc can be either converted directly into electricity in a zinc/air battery or it can be processed with steam yielding hydrogen [10–12]. In both cases the initial ZnO is reformed and can be used again in the high temperature solar step. Presently, the separation of the products of the high temperature step, that is, the stoichiometric mixture of zinc vapor and oxygen, is attempted by subjecting the mixture to a rapid cooling by injecting cold inert gas. The product of the fast quench is a mixture of partially oxidized zinc particles. The average chemical composition of the product varies with the quench conditions such as the initial partial pressures of the reactants or the quench rate. Zinc yields up to 94% have been achieved [13, 14] in highly diluted zinc/oxygen mixtures.

To advance our understanding of the quench process and to identify process conditions that result in a maximum energy conversion efficiency a kinetic model to describe the time evolution of a mixture of zinc vapor and oxygen subjected to rapid cooling was developed and implemented as computer program. We chose not to apply a model based on chemical nucleation such as has been applied for example, to the plasma processing of silicon films [15] because a reaction scheme that relies on the formation of small Zn_{x}O_{y} molecules [16–18] seems questionable. Although theoretical [19, 20] and experimental evidence [21, 22] exists that small Zn_{x}O_{y} clusters are stable molecules, ZnO or ZnO_{2} has been identified until now only under rather exotic conditions. It was reported [23] that Zn atoms do not react with O_{2} at ambient conditions due to the highly endothermic heat of formation of molecular ZnO of about 340 kJmol^{−1}.

#### 2. Theoretical

In the following, a model that describes the reactions that occur in a mixture of zinc vapor and oxygen subjected to rapid cooling is developed. Here, the focus is on the general traits of the theory. For a thorough treatment including details of the mathematical and numerical aspects the reader is referred to [24].

##### 2.1. Reaction Scheme

The model is based on the assumption that a mixture of zinc vapor and oxygen is metastable below the decomposition temperature of zinc oxide ( K) and thus does not react in absence of surfaces. This assumption has been validated several times [5, 25]. The surface necessary for reactions is provided in the model by small zinc droplets formed by homogeneous nucleation of zinc vapor (reaction 0, Figure 1). Once zinc droplets of a fixed size have formed, surface is available and condensation of zinc becomes possible on pure zinc droplets (reaction 1, Figure 1) as well as on partially oxidized droplets (reaction 1′). If zinc is exposed at the surface, oxygen molecules start to react with the free metal surface reactions 3 and 3′, Figure 1) forming a monolayer of zinc oxide on the droplet. It is obvious that reactions 3 and 3′ represent a gross simplification of the reality. Formation of an oxide monolayer involves several distinct steps, such as the adsorption of oxygen molecules, their dissociation, and, finally, the formation of Zn–O bonds. However, such a simplification is necessary to keep the model simple enough to be numerically treatable and will require that the rate constant of reactions 3 or 3′ corresponds to the slowest of its elementary steps. The last surface reaction (reaction 2, Figure 1) describes the addition of zinc atoms to the oxide layer. During this reaction surface oxygen atoms are transformed into embedded oxygen atoms and free metal surface area is recreated.

Droplets that occur in the reaction scheme are characterized by three numbers: the total number of zinc atoms (), the total number of oxygen atoms (), and the number of embedded oxygen atoms (). The droplets are assumed to be spherical with volume where and denote the atomic volume of liquid zinc and the volume occupied by a ZnO unit in zinc oxide, respectively. Furthermore, the surface area occupied by an oxygen atom as obtained from the crystal structure of ZnO is required to calculate the oxidized and the unoxidized surface areas, respectively,

Reactions 2 and 3′ together describe the direct formation of ZnO from the gas phase with a growth in direction of its crystallographic axis that is characterized by alternating zinc and oxygen layers. Note that the bulk oxidation of liquid zinc is not part of the model. It involves the diffusion of either zinc or oxygen through the oxide layer being formed and is a much slower process than the direct formation from the gas phase. Contradictory data has been reported as whether the diffusion of zinc [26, 27] or of oxygen [28, 29] in the oxide layer is rate limiting. The oxidation of solid zinc, however, seems to be limited by zinc diffusion [30].

##### 2.2. Rate Laws and Rate Constants

Rate laws and rate constants are required to describe the kinetics of each of the reactions 0–3 of Figure 1. For reaction 0 the stationary nucleation rate of the classical nucleation theory [31] was chosen. A numerical scaling factor () for the nucleation rate is introduced which assumes the role of the rate constant. Other expressions for [32, 33] could be substituted. For reactions 1–3, we use as rates the product of the active surface area , the sticking probability , and a collision frequency. The latter is obtained from the statistical description of an ideal gas and depends on the partial pressure , the temperature , and the monomer mass : with denoting net rates. Thus, (4) represents the net condensation rate of zinc as it contains , the difference between the partial pressure of zinc, and its saturation vapor pressure. Furthermore, contains the step function to prevent reevaporation of zinc. Therefore, no sink for particles is present and the total number of particles in the system will grow monotonically. For and the saturation pressures of zinc and oxygen respectively, above a ZnO layer are assumed to be negligible. As seen from the Ellingham diagram [34] the oxygen partial pressure above ZnO at the boiling point of zinc is of the order of Pa. We then assume that the equilibrium pressure of zinc above ZnO is likewise low in agreement with measurements of the equilibrium pressures of Zn and Se above ZnSe [35]. Temperature and partial pressures depend on time, and the rates can be separated into a time-dependent part and into a part that depends on the chemical composition of the droplets. However, this can be achieved only if the Kelvin effect is neglected. The Kelvin effect would only affect the evaporation of zinc. It can be included by multiplying (4) with . Neglecting the Kelvin effect seems justifiable because the zinc droplets become surface oxidized to a large extent soon after their creation, and the enhanced evaporation is only in effect on the metallic surface.

Although reaction 3 implicitly contains the dissociation of an O–O bond, does not seem to reflect this as it does not contain an activation term such as that present in, for example, an Arrhenius-type rate law. Nevertheless, it is functionally equivalent to an Arrhenius rate. One can identify the collision frequency in (6) with the preexponential factor and the sticking probability with the exponential term as the latter represents the fraction of successful collision in a thermal population. A notable difference between and an Arrhenius-type rate law lies in the temperature dependence. In the collision frequency is temperature dependent while is a constant. In contrast, the temperature dependence in the Arrhenius rate law is solely in its exponential term.

Microscopically reactions 1, 2, and 3 follow the scheme as exemplified for reaction 2: where the reaction current is the product of the rate and the number density of the droplets . The resulting time evolution is a coupled system of first order linear differential equations, one equation for each type of droplet occurring. A closer analysis [24] reveals a hierarchic structure of five different classes of droplets: (class 1), (class 2), (class 3), (class 4), and (class 5). Class 1 represents the just nucleated droplets containing zinc atoms while, for example, class 4 represents surface oxidized droplets containing no embedded oxygen atoms (). Only droplets within the same class are coupled through differential equations with the number densities of lower classes appearing as source terms. The hierarchic structure is an important guideline in the solution strategy sketched in the following.

Coagulation, that is, the fusion of two droplets, cannot be treated within the model presented as no unique reaction can be derived for this process:

While the number of zinc and oxygen atoms in the fused droplet corresponds to the sum of the respective numbers in the parent droplets, no unique relationship for the distribution of the oxygen atoms into surface and embedded atoms exists. Their distribution depends on the collision geometry and is arbitrary. If coagulation must be included, it is suggested to use the result of the model presented as input for a separate model treating only the coagulation process [36]. This is justified as coagulation occurs on a much slower time scale than condensation; that is, [37]. As soon as condensation ceases, the chemical composition of the individual droplets remains constant and the distinction of surface and bulk oxygen atoms, that is, index , becomes unimportant.

It is important to realize that the quantity of interest is not the number density that would deliver an unmanageable detailed amount of information but rather more global properties such as the total number density of droplets in a certain -space or the average chemical composition of the product. Generally speaking, such properties are sums of the form

##### 2.3. Solution Strategy

For the reasons stated previously, and because of the lack of a truncation strategy for the infinite system of differential equations, solution strategies heading directly for the density are not promising. Therefore, we chose a different approach. First, the problem is reformulated for continuous variables . This leads to partial differential equations for , specific for each class, except for the single equation (see (3)) for class 1 that remains unchanged. The density functions of the lower classes become boundary values of the densities of the higher classes according to the class hierarchy mentioned. Properties (see (9)) are now expressed by integrals that are again specific for each class. In a second step, the solutions of the partial differential equations for each of the classes 2–5 are expressed by their characteristic curves. The latter are solutions of a system of ordinary differential equations. The class hierarchy induces child parent relations among the characteristic curves as shown in Figure 2. Finally, the property integrals, until now formulated in the -space, must be transformed according to the classwise different parametric representations of the solutions of the partial differential equations. The transformed integrals obtained are based on a small number of time-like parameters with the simple integration domains for classes 2 and 3, for class 4, and for class 5.

##### 2.4. Calculation of Time Resolved Properties of the Droplet Ensemble

The numerical calculations of the characteristic curves and of the property integrals are based on a finite number of time steps:

A new set of characteristic curves starts at each time . The ensemble of these characteristic curves forms a tree with a rapidly growing number of branches proportional to as reported in Figure 2. Although may become quite big, the structure remains manageable because all integrals contain factors exponentially decreasing with time. This allows for an efficient “pruning” strategy strongly reducing the initial growth of the number of characteristic curves, that is, the number of ordinary differential equations that must be simultaneously integrated.

In Figure 3 the influence of pruning on the execution time and memory requirements as well as the resulting absolute and relative errors of properties ( is used as indicator) is reported. The functions denote the exponentially decreasing factors in the integrands. Pruning implies that curves with were discarded. Applying discards curves when their contribution drops below 13.5% of their initial value and introduces a maximum error of 685 Pa (14.1%). Using as threshold, that is, , decreases the error to 9.6 Pa (0.141%) and reduces the computation time by a factor of almost 1800 from about 10.5 h if no pruning is applied to 20 s. The memory requirements that scale directly with the number of curves required drops by a factor of about 2500.

**(a)**

**(b)**

Since only a finite set of characteristic curves is at disposition, the property integrals must be approximated by weighted sums over the available functional values of the integrands. It is gratifying for this task that the integration domains in the characteristic representation are simple. The choice of the weights requires some considerations because the time increments for are determined by the accuracy requirements of the computing algorithm and will not be equal. The integration domains (classes 2, 3), (class 4), and (class 5) are, therefore, subdivided into elementary domains of different sizes and shapes by the partition 10 of the time interval . To attribute equal weights to the corners of a given elementary domain is the best choice and assures that locally linear functions are exactly integrated. As for the total integration several elementary domains participate at a given point of intersection. To determine the weights of these points, therefore, requires some book keeping [24].

The chosen strategy transforming a discrete problem into a continuous one that in turn must be discretized for its numerical evaluation may appear indirect at first sight and requires a considerable amount of analytical work. However the analysis provides the proper truncation rules and numerically efficient means to calculate the desired physical properties.

#### 3. Examples

In this section the power of the model presented is illustrated by analyzing the time evolution of a stoichiometric mixture of zinc vapor (11000 Pa) and oxygen (5500 Pa) subjected to a cooling rate of −10000 Ks^{−1}. Reactants consumed are thought to be replaced by inert gas or, alternatively, the total pressure of the system decreases. Because numerical values for the rate constants of the four processes involved in the model are not known, we arbitrarily set in (3) and in (4)–(6). Additional input parameters required were and m^{2}, the latter corresponding to the surface area covered by a single oxygen atom. Pruning was applied on a very conservative level by setting the thresholds .

In Figure 4 the time evolution of the saturation , the nucleation rate , and the total number of particles formed per volume unit are reported. It can be seen from this figure that nucleation only becomes significant once the saturation becomes sufficiently large, here larger than about 4.2. Once particles have been formed, zinc is consumed, begins to decrease, and nucleation ceases very rapidly. Two different values for are reported in Figure 4. One has been calculated as a property of the droplet ensemble according to (9) with while the other value has been calculated from the time integral of the nucleation rate. It can be shown [24] that both values should be identical. In reality, however, there are numerical errors associated with both values. Pruning results in underestimating the total number of particles formed. On the other hand, is represented by an extremely steep function that covers almost 18 orders of magnitude. Its numerical integration is certainly demanding. We thus consider the deviation of about a factor of 2.4 to be acceptable, especially in view of the small errors introduced in as reported in Figure 3. The deviation between the two values can be reduced by augmenting the pruning thresholds and reducing the integration steps, however at the expense of a considerable increase of computation time.

In Figure 5 the time evolution of the partial pressures of the reactants zinc and oxygen is reported. These values are calculated from the initial partial pressures and the number of zinc atoms and oxygen molecules consumed ( and , resp.). Comparing with Figure 4, it is evident that both partial pressures begin to decrease significantly only after has obtained its peak value. This again illustrates the extremely strong dependence of the nucleation rate on the supersaturation . While zinc is consumed completely, not all available oxygen is consumed. At the end of the reaction, at s, remains at 1115 Pa; that is, oxygen is only consumed by about 80%. This is a consequence of the reaction scheme that allows oxygen to be consumed only by reaction with the metallic surface of the droplets. In other words, the oxide layer grows only by forming ZnO from the gas phase. This reaction requires both reactants, zinc and oxygen. The growth of the oxide layer by bulk oxidation is not included in the model as explained earlier. It requires diffusion of either zinc [26, 27] or oxygen [28, 29] through the oxide layer, a process that occurs on a much longer time scale. The fact that oxygen is not completely consumed can also be seen from the chemical composition () reported in the same figure. The degree of oxidation only reaches values close to 80% towards the end of the reaction. During the nucleation phase it is lower, around 65%, due to the ongoing nucleation of pure zinc droplets.

Alternatively, the incomplete consumption of oxygen can be understood by analyzing the various contributions to the total surface area of the droplet ensemble as reported in Figure 6. Initially both, , the total metallic surface per volume unit (), and , the total oxidized surface per volume unit () grow in parallel. As zinc is consumed, ceases to grow and starts to decrease until it finally becomes negligible at s. Around this time has becomes negligible and the metallic surface of the droplets cannot be reformed anymore. As soon as , the oxygen consumption rate (6) becomes zero; that is, oxygen consumption stops as soon as the surface of the droplet ensemble is completely oxidized.

Up to here, only properties that are represented by a single number (scalar properties) were presented. More complicated properties, such as distribution functions, can also be obtained. This is illustrated by the distribution of particle masses in Figure 7 where four distributions are reported corresponding to the times indicated by small vertical arrows in Figures 4–6. For Figure 7 the number of particles per volume unit (; ) as well as the particle mass (; ) needs to be calculated. The mass axis is divided into an arbitrary number of small intervals and the tuples are binned accordingly. The distribution is scaled to ensure that its integral matches the total number of particles per volume unit.

As the nucleation rate peaks (curve 1, Figure 7), the distribution is strongly asymmetric dominated by mass kg, the mass of the nucleating droplet containing 86 zinc atoms. Already at 10.4 ms (2) when both and start to decrease, the distribution has become bimodal, one maximum at the -cluster and one around kg containing grown particles. Later, as nucleation has ceased (3), the distribution becomes symmetric finally approaching a log-normal distribution (4). Note that unlike in a single component system the pmd cannot be converted into the corresponding particle size distribution (psd): the psd has to be calculated from (9) as no unique relation between the mass of a particle and its diameter exists due to the variable chemical composition of particles.

#### 4. Comparison with Selected Experimental Results

Rate constants for reactions 0–3 are required for a quantitative verification of the model presented. Only then quantitative predictions can be obtained and compared with experimental data. Presently, values for these rate constants are not available. Therefore, we apply and as in the previous section and limit ourselves to a qualitative comparison. Any property of the partly oxidized zinc product can be obtained from the model in principle, but only little suitable experimental data is currently available. In the following the average chemical composition of the product is considered. Its dependence on the quench rate and on the initial partial pressures of the reactants is calculated for quench rates in the range 100000 Ks^{−1} to 5000 Ks^{−1} at an initial dilution of and, for a quench rate Ks^{−1}, for dilutions . These dilutions correspond to initial partial pressures in the ranges Pa; to Pa; . The results are summarized in Figure 8.

**(a)**

**(b)**

From Figure 8 we see that within the parameter range chosen the zinc yield can only be influenced little. A definitive trend towards higher zinc yields with increasing dilution exists. Dilution of a stoichiometric mixture of zinc and oxygen with inert gas by a factor 5 increases the zinc yield from 17.7% to 19.8%, a relative increase of only 12%. No clear trend of the zinc yield is observed by varying quench rate. The latter finding was explained earlier [38] by analyzing the path followed by in the phase diagram of zinc during the quench as follows. First, the net rate of condensation of zinc is proportional only to the difference of and the corresponding equilibrium pressure (see (4)). The reactions that comprise the growth of the oxide layer, however, are favored as the equilibrium pressures of both zinc and oxygen are essentially zero above ZnO. Second, the growth of the oxide layer, once initiated, continues also at for the same reasons. Both effects result in comparably small zinc yields almost independent of the quench rate because on a relative scale the time spent in the region of the liquid phase where zinc condensation is in effect remains the same for all quench rates.

Similarly, the increasing zinc yield observed when working at lower partial pressures can be understood from Figure 9. We note that , the slope of the liquid-gas coexistence line, increases H faster than exponential in the phase diagram (more than linear in Figure 9). Increasing the dilution shifts the quench curves to lower temperatures and lower pressures. The curves recross the phase boundary where its slope is smaller (compare the points for and indicated by circles in Figure 9). Therefore, we expect the time the system spends in the region of liquid zinc to increase with increasing dilution which results in higher zinc yields since condensation of zinc can only occur during this time. As the slope of the solid-gas phase boundary is much smaller than the one of the liquid-gas boundary, we expect much higher zinc yield if nucleation is delayed until resublimation occurs. In this case, we can even hope that the recrossing into the region of gaseous zinc can be avoided completely.

Experimentally, a trend towards higher zinc yields (as determined by X-ray powder diffraction) with increasing dilution was initially observed with the H solar thermogravimetric analyzer (Solar TG) setup [14]. If the amount of hydrogen gas evolved upon treating the product with HCl is used as measure for the zinc yield [39], the zinc yield is generally about 30% lower than by the XRD method. Earlier data [14] obtained from the ZIRRUS thermochemical reactor [9] are included in the evaluation here. Together, they provide a more consistent data set that is summarized in Figure 10. For both experiments the temperature drop is induced by turbulent mixing with cold inert gas. The temperature gradient experienced by the zinc/oxygen mixture cannot be directly measured in both cases. Therefore, the quench rate was determined from averaging the temperature history along selected streamlines obtained from CFD calculations. A clear trend towards higher zinc yield is observed in Figure 10 as predicted by the model. While a weak trend towards higher zinc yield with higher quench rate was observed in [14], no such trend is reported in H [39].

#### 5. Summary and Conclusions

We have developed a theoretical model to describe the time evolution of a mixture of zinc vapor in presence of oxygen subjected to conditions where homogeneous nucleation of the metal is induced. Once zinc particles have nucleated, they grow by either condensation of the metal or by the direct growth of an oxide layer from the gas phase. In traditional nucleation and condensation models, the individual particles are characterized by a single number, the number of monomers they contain. Thus, these models only describe the formation and growth of particles with a fixed chemical composition mirroring the one of the monomer. In our case, the average chemical composition of the particles can vary and particles are characterized by three numbers (). This vastly increases the complexity of the problem. Calculating the properties of the product not from but by summing the contribution of the individual characteristic curves allows to formulate a rational criterion to select characteristic curves that do not contribute in contrast to, for example, sectional models [36, 40] and leads to a constant number of characteristic curves shortly after the nucleation phase. Examples of properties that can be obtained from the model are presented. The theory has been used to predict the dependence of the average chemical composition of the product on quench rate and on the initial partial pressures of the reactants with experimental results. Both, theory and experiment, find that the chemical composition of the product is largely independent of the quench rate but that a higher zinc content is found when the reactants are diluted with inert gas.

#### Acknowledgment

This work was supported by the Swiss Federal Office of Energy (SFOE).