Pattern Formation in Porous Media via the Liesegang Ring

Ind. Eng. Chem. Res. 2004, 43, 3073-3084

3073

Pattern Formation in Porous Media via the Liesegang Ring Mechanism M. I. Lebedeva,† D. G. Vlachos,*,† and M. Tsapatsis‡ Department of Chemical Engineering and Center for Catalytic Science and Technology, University of Delaware, Newark, Delaware 19716, and Department of Chemical Engineering and Materials Science, University of Minnesota, Minneapolis, Minnesota 55455

A model for the spontaneous formation of patterned deposits in porous media, reminiscent of Liesegang rings (LR) formation in gels, is presented. Simulations show that under certain conditions the frontal deposition of a continuous solid film changes to a quasi-periodic pattern formation of bands. Bifurcation analysis of a simpler, skeleton model explains LR formation as instability of a uniformly propagating plane reaction front to a time periodic solution. A theoretical stability criterion is developed that suggests that key parameters in pattern formation are the critical concentration for nucleation and the speed of the autocatalytic growth reaction. While our theory is not aimed at explaining all LR formation experiments in colloidal systems, it can be applied in the case of strong concentration gradients. Using our theory, the uneven spacing law of LR bands is explained as a consequence of the time-varying velocity of the moving reaction front. Some apparent contradictions of previous LR models are also explained, and suggestions for experimentally controlling pattern formation are made. Introduction Patterns during materials deposition have been demonstrated experimentally for more than a century and are known as Liesegang rings (LR) or bands.1,2 The terms periodic, rhythmic, or recurrent precipitation are also often used despite the lack of truly periodic patterns. In most experiments, bands obey an uneven spacing law; i.e., the ratio pn ) xn+1/xn tends to a constant value (typically greater than 1) as n increases. Here, xn denotes the position of the nth band from the interface of the two electrolytes. LR patterns form when a soluble reactant (typically an ion of a salt) diffuses from either the center or the periphery of a medium (often a gel) uniformly filled with a second soluble reactant (typically an ion of another salt) to produce an insoluble substance. There exist a number of models for explaining LR pattern formation,3-15 most of which can be grouped into two classes: prenucleation and postnucleation theories. Prenucleation models essentially employ the mechanism proposed by Ostwald.3 According to these models, when a critical supersaturation is reached, precipitation starts and growth depletes neighboring regions of both ions. So, nucleation occurs as a spatially discontinuous process. Further deposition occurs as a concentration wave propagates, leaving behind bands at space and time intervals essentially determined by diffusion. This theory4-6 assumes a priori existence of sharp bands, and while it neglects reaction kinetics, it is able to predict the spacing law. A priori existence of bands is essential to these models. For example, numerical solution of a prenucleation type of model, which included not only diffusion but also reaction kinetics, did not produce bands.9 * To whom correspondence should be addressed. Tel.: (302) 831-2830. Fax: (302) 831-1048. E-mail: [email protected]. † University of Delaware. ‡ University of Minnesota.

The prenucleation theory has been criticized by supporters of the postnucleation models, e.g., refs 8, 10, and 11. The latter models employ the assumption that a macroscopic pattern develops after an initially spatially homogeneous nucleation of solid particles has occurred because of the instability of particles with respect to their size, a phenomenon driven by free-energy minimization (Oswald ripening). Postnucleation models have led to the conclusion that in prenucleation models spatial discontinuities in nucleation were obtained by unwarranted approximations or ad hoc assumptions. However, it has been noted in ref 12 that “band formation occurs for both continuous and discontinuous nucleation space provided that the nucleation law has a sufficiently steep dependence on the degree of supersaturation past a critical value”. The postnucleation models are able to explain complex pattern formation, which arise in nongradient systems,8 and the evolution of already formed patterns. However, they cannot predict band formation, although there has been an attempt to predict the spacing law by employing additional hypotheses.11 In the absence of concentration gradients of the two electrolytes, a structure may still form. The idea to explain structure formation in this case as a manifestation of Turing instability was proposed first in ref 16. Although that model was not develop afterward, the approach of ref 16 was applied in refs 8 and 10 to model the spontaneous precipitated patterns that arise as a uniform sol ages. More recent models13,14 have taken into account both of the aforementioned processes, i.e., a propagating supersaturation wave and the dissolution of small and growth of large particles. Numerical simulations have confirmed the spacing law for parallel bands’ formation. Furthermore, more complex features, such as dislocations and helices, were explained qualitatively. Recent theoretical work17-20 has been devoted to the derivation of the spacing law that could be compared quantitatively with experiments. Besides the aforementioned theories,

10.1021/ie049957r CCC: $27.50 © 2004 American Chemical Society Published on Web 04/23/2004

3074 Ind. Eng. Chem. Res., Vol. 43, No. 12, 2004

a spinodal decomposition mechanism for LR formation has also been proposed.19,20 The pre- and postnucleation models have been contrasted in a recent review of selected experiments on spatial structure formation.21 The authors find extensive experimentally driven evidence for the postnucleation model. However, they do not analyze band formation in aerogels15,22 and experiments demonstrating propagating Liesegang strata.23-25 The first one is in qualitative agreement with prenucleation theory, and the second is described well by a prenucleation theory model.14 Despite intensive theoretical work on LR formation, there remain a number of experimental observations, which cannot be explained by published models. For example, the spontaneous formation of periodic bands (regularity in band spacing) of nanocrystalline titania in mesoporous Vycor glass using chemical vapor deposition26,27 cannot be explained with current models. While pattern periodicity may generally be a result of the Turing-type instability mentioned above (for a recent extension of the Turing mechanism to materials deposition, see ref 28), time-dependent and more extensive experiments29 clearly show that the number of bands increases with time and, under other conditions, the spacing is consistent with the aforementioned spacing law. Because of these findings in the TiO2 experiments, we focus on the prenucleation LR pattern formation mechanism instead of the Turing mechanism. Because of the large number of parameters of such models, it is difficult to develop strategies for experimental control based solely on simulations. Development of theory can be valuable in that regard. Among George R. Gavalas’ contributions in the area of reaction and transport in porous media are modeling and experimental studies of deposit formation in mesoporous glass by counterdiffusion chemical vapor deposition. The mathematical modeling and theoretical results we present here are motivated by experimental observations first made in his laboratory at Caltech in the early 1990’s26 and attempt to relate spontaneous pattern formation in such experiments with the LR phenomenon. Although such pattern formation is undesirable for membrane formation, because it leads to thick and poorly localized deposits, it may be of practical significance for the fabrication of functionally graded micro- and nanostructures with applications like reflective coatings and photonic crystals. In the present work, motivated by the experiments of ref 27, we suggest a model for the formation of patterns in a porous substrate and present selected simulations demonstrating pattern formation. Then a simpler, skeleton model for LR formation is proposed, on which theoretical bifurcation analysis can be performed. Unlike previous models,8,10,16 we perform stability analysis of a spatially inhomogeneous traveling wave solution rather than a spatially homogeneous solution. Numerical simulations and theoretical analysis demonstrate that Ostwald’s assumption of discontinuous nucleation is necessary (but not sufficient) for our model. Note that our objective is to understand the onset of periodicity of deposition in porous media and, thus, our theory is not aimed at explaining all experiments in traditional colloidal systems discussed in ref 21. Nevertheless, our theory can be applied to the classical LR experiments in the case of “strong gradients”. Our approach allows one to estimate the location where

bands start and explains both the uneven spacing law of band formation and the possibility of nearly periodic patterns. It also suggests a connection with the appearance of helical structures mentioned in refs 13 and 14. Finally, it also provides insight into why the numerical simulations in ref 9 did not result in LR bands and why the predictions of refs 4-6 are right despite neglect of the reaction kinetics. Deposition in Porous Media Prototype Nucleation and Growth Model. This model is the first effort to simulate LR pattern formation in porous media. It is based on the premise that an intermediate species C is formed by transport-chemistry interactions of the main reactants A and B and then nucleates and grows to solid deposit S. The proposed reaction scheme is as follows: k1

A + B 98 C, homogeneous reaction k2

(1)

C 98 S, nucleation

(2)

C 98 S, growth autocatalysis

(3)

S, k3

Reaction (2) is the nucleation step. Furthermore, S catalyzes the heterogeneous growth reaction (3) because of an increase of its surface area with time. As an example of a specific system, in the titania deposition system,27 TiCl4 + 2H2O f TiO2 + 4HCl, our model may have the following analogy: A, H2O; B, TiCl4; C, Tin(OH)n(4-x)Clnx (x ) 1-3); S, TiO2. The reduction in porosity of the substrate due to deposition creates a barrier for diffusion of species in the vicinity of the deposit. This eventually may alter the conditions for the formation of LR especially at nearly pore closure conditions. Previous work on LR formation has not accounted for the effect of pore blockage on the transport of the reactants, which is considered here. We suggest that it may also be important to consider such effects in the modeling of LR pattern formation in gel-type media because of the constantly evolving gel microstructure as the reaction proceeds. Because of the molecular size and dissolution differences, it is common that one of the reactants, say A, can diffuse through the solid deposit as well but with a lower diffusion coefficient compared to that in the open pore. The governing equations describing the transport and kinetics in the deposition within the porous medium are

[

]

∂Ci ∂ui(Ci,) ∂ ) d () + fi(Ci,F), i ) A-C (4) ∂t ∂x i ∂x ∂F ) fS(Ci,F) ∂t

(5)

uA() ) CA + (0 - )K-1CA, dA() ) D0A + (0 - )K-1D1A (6) ui() ) Ci, di() ) Di, i ) B and C

(7)

fi ) -W1, i ) A and B, fC ) W1 - W2 - W3, fS ) W2 + W3 (8) F ) φDF0, ) 0 - φD

(9)

Ind. Eng. Chem. Res., Vol. 43, No. 12, 2004 3075

Here, Ci, i ) A-C, are the molar concentrations of species A-C in pore fluid, F is the deposit density in number of moles of the deposit per unit volume of the porous medium, is the porosity, φD is the deposit volume fraction in the solid phase, F0 is the constant (true) density of the deposit in number of moles of the deposit per unit volume of the deposit, K is the Henry’s constant of component A solubility in the deposit, D0A, Di, i ) B and C, are the diffusion coefficients of components in pore fluid, and D1A is the diffusion coefficient of component A in the deposit. Finally, t denotes time and x is the spatial coordinate. Note that in our notation the diffusivities are reduced by a tortuosity factor. The form of the reaction rate function Wi depends on the type of the corresponding reaction. It is assumed that

W1 ) k1CACB

(10)

W2 ) k2CC2H(Θ2)

(11)

W3 ) k3CCφD2/3H(Θ3)

(12)

Here H(Θi), i ) 2 and 3, are the corresponding step functions of Θ2 ) CC - C0C and Θ3 ) - min, where min is some sufficiently small positive number. The function H(Θ2) implies that nucleation begins when the concentration of C exceeds some critical value C0C, and the function H(Θ3) ensures that is not physically possible to deposit product when the porosity is sufficiently small. Finally, the term φD2/3 in eq 12 describes the dependence of growth rate on the surface area of the solid phase. The effect of the initial and boundary conditions that resemble the experiments of ref 27 along with comparison to the experiments will be discussed elsewhere. Here we focus on conditions that are close to those of the previous models of LR pattern formation. The initial conditions are a uniform concentration of B, C0B, in the domain 0 e x e L, whereas A, C, and S are not initially present. Zero-flux boundary conditions are applied to both species B and C at both ends and to A at x ) L. For all practical purposes, the concentration of A at x ) 0, C0A, is assumed to be constant. The governing equations (4) and (5) are dimensionalized and discretized using a second-order finitedifference scheme in space and an implicit backward Euler scheme in time. Newton’s method and a block tridiagonal algorithm are used for solving the corresponding difference equations. Illustrative Simulation Results. Simulations show two distinct deposition modes. In the first growth mode, termed frontal propagation, a deposit forms as a result of a concentration wave propagating through the entire substrate. This deposit is spatially continuous without any bands, as shown in Figure 1a,b (the dimensionless rate constant and supersaturation are defined as γ′ ) k3/k2C0B) and y′c ) C0c /C0B). In the second growth mode, termed quasi-periodic, bands form from the beginning. The term quasi-periodic, used often to describe oscillations in dynamical systems, reflects some (but not complete) regularity in band spacing. Finally, quasiperiodic deposition often starts after some induction time, as shown in Figure 1c,d; i.e., a mixed mode is

Figure 1. Frontal deposition according to the prototype model (4)-(12) when (a) the autocatalytic step (3) is absent (k3 ) 0) and (b) k3 > 0 but no nucleation barrier exists (C0C ) 0). Examples of transition from frontal to quasi-periodic deposition as time increases are depicted in parts c and d. The induction period for the transition decreases with increasing speed of the autocatalytic step. The parameters are D1A/D0A ) 3.3 × 10-4, DB/D0A ) DC/D0A ) 0.11, C0B/F0 ) 4 × 10-3, C0A/C0B ) 5, K-1 ) 10-7, 0 ) 0.5, min ) 0.01, t0 ) L2/D0A ) 370.4, and k1C0Bt0 ) k2C0Bt0 ) 3.7 × 103 with (a) γ′ ) 0, y′c ) 0.1, (b) γ′ ) 50, y′c ) 0, (c) γ′ ) 50, y′c ) 0.1, and (d) γ′ ) 100, y′c ) 0.1. The dimensionless parameters are y′c ) C0C/C0B, γ′ ) k3/k2C0B, and τ ) t/t0.

observed where, for longer times, the frontal growth changes to quasi-periodic deposition. The observed growth mode depends on the choice of parameters. We notice two of these parameters: the critical concentration C0C and the ratio k3/k2 of autocatalytic growth and nucleation rate constants. Numerous simulations we have performed resulted in LR bands only when both parameters were nonzero and sufficiently large, as shown in Figure 1. The effect of the parameter C0C is expected, and it is in agreement with Ostwald’s nucleation theory. However, it is interesting to note that there appears to be a criticality associated with the ratio k3/k2. Our simulations indicate that the frontal growth changes to quasi-periodic deposition when reaction (3) is sufficiently faster than reaction (2). So, it is impossible to obtain LR patterns when the autocatalytic reaction (3) is neglected. These observations from our simulations are rationalized below with theoretical analysis of a simpler but related model. This criticality is the reason nucleation in the model of ref 9 occurs as a spatially continuous process (reminiscent of our Figure 1a,b) although the sink term (W2 + ηW3) in our model formally differs from the one of ref 9. An increase in the ratio k3/k2, shown in Figure 1c,d, can result from a change in the deposition conditions. For example, if the activation energy of the autocatalytic reaction (3) exceeds the one of the nucleation reaction (2), then the ratio k3/k2 rises with increasing temperature. Such variations may then be the reason for the change from frontal deposition to quasi-periodic mode observed in experiments27 when the deposition temperature was changed. Finally, we should note that for a semibatch system and realistic concentration reagents


and TiO2 density, i.e., realistic C0B/F0, little deposition happens (from pore blocking), a situation that is consistent with experimental observations. Therefore, it is expected that the effect of porosity is small for these simulations. However, pore blocking may happen for systems with a higher C0B/F0 or for continuous-flow, open systems. Skeleton Deposition Model To understand the above results, here we focus on the role of the autocatalytic reaction and nucleation in LR formation by considering a slightly simpler model introduced in ref 30 that, with some assumptions, is theoretically tractable. We simplify the prototype model (4)-(12) while retaining its main characteristics. Thus, we consider two parallel reactions: k1

A + B 98 S

(13)

S, k2

A + B 98 S

(14)

The dimensionless governing equations are as follows:

∂y1 ∂2y1 ) 2 - f(y1,y2,z) ∂t ∂x

(15)

∂2y2 ∂y2 ) D 2 - σf(y1,y2,z) ∂t ∂x

(16)

∂z ) σf(y1,y2,z) ∂t

(17)

and the initial and boundary conditions are the same as those in the prototype model (4)-(12). The following dimensionless quantities are used:

t)

t′D0A CA CB CS x′ , x ) , y1 ) 0 , y2 ) 0 , z ) 0 , 2 L L CA CB CB D)

DB Ks k1C0BL2 , y ) , κ ) , c D0A C0A C0B D0A γ)

C0A k2C0B , σ) 0 k1 C B

Here CS is the concentration of S, Ks is the solubility product, L is the domain length, x′, x and t′, t are the real and dimensionless length and time coordinates, respectively. Other symbols are the same as those in model (4)-(12). In this model, we do not take into account the change of porosity because it does not affect the main conclusions about the role of the autocatalytic reaction. However, realistic modeling requires consideration of this factor, especially near the percolation threshold. Finally, we assume that growth occurs when the product CACB exceeds some critical value Ks, as has often been assumed.4-6,9,12 Thus, the form of the dimensionless reaction function f is

f ) κy1y2[H(θ) + γz] Here H(θ) is a step function

(18)

H ) 1 at θ g 0, H ) 0 at θ < 0, θ ) y1y2 - yc (19) Note that, for simplicity, we have left out details from the surface area dependence of reaction (14) while we retain its autocatalytic nature. However, the power law dependence of autocatalysis is examined near the end of the paper. Typical simulation results using this skeleton model are depicted in Figure 2. Comparison of Figures 1 and 2 indicates that the results of this skeleton model are qualitatively similar to the more complex prototype model described above. For example, both growth modes are observed and the occurrence of bands is often observed as a transition from the frontal deposition mode after some induction time. We focus then on this skeleton model (15)-(19), which is amenable to theoretical analysis. A brief outline of the derivation is given in ref 30. Below, we present a more comprehensive derivation, give further numerical evidence of the theory, and finally extend the model to provide insight into the experimental control of patterns. Bifurcation Analysis of the Skeleton Model Simulations indicate that the solutions of the problems (4)-(12) and (15)-(19) evolve in time and pass through various states, as shown in Figures 1 and 2. For example, the frontal growth mode changes to quasiperiodic for some parameter values. Finally, as t f ∞, the concentrations of species of the skeleton model (15)(19) tend to their stationary state y1S ) 1, y2S ) 0, and zS ) zS(x). As parts a and b of Figure 2 indicate, under some conditions the concentration profiles of reactant B and deposit S have traveling wavelike features; i.e., as the time increases, the profiles shift along the x coordinate parallel to themselves. However, unlike real traveling waves, these waves move with a nonconstant, slowly varying velocity ω(t). As a matter of fact, it is well-known that ω ∼ λ/xt, where λ is some nonlinear function of the initial concentrations C0A and C0B and of the diffusion coefficients DA and DB (see, for example, ref 10). Our simulations show that each such wave, along with its instant velocity ω(t), is asymptotically close (except for translation) to the real traveling wave solution, y2 ) Y2[x - ω(t) τ], z ) Z[x - ω(t) τ], of the system

∂2y2 ∂y2 ) D 2 - σF(y2,z,t) ∂τ ∂x

(20)

∂z ) σF(y2,z,t) ∂τ

(21)

with boundary conditions

y2 ) y- and z ) z- at x f -∞; y2 ) 1 and z ) 0 at x f ∞ (22) where 0 e y- < 1 and z- > 0 are constant, t is a parameter, and ω(t) is the constant front velocity for this solution corresponding to t. F(y2,z,t) ) f{y1[y2(x,t)],y2(x,t),z(x,t)} is determined by solving eqs 15-19 at this fixed t; i.e., because y2(x,t) is monotonic, x is written as a function of the independent variable y2 and y1(x,t) in terms of y2, namely, y1(x,t) ) y1(y2,t). Such a function can be constructed numerically for each value of t. Thus, eqs 20 and 21 are expressed in terms of y2 and z only.


model depends on the character of QSS. Similar behavior is well established for the time evolution of autocatalysis in a chemical reactor.31 In particular, for cubic autocatalysis with catalyst decay

P f A, A + 2B f 3B, A f B, B f C

Figure 2. Numerical solutions using the skeleton model (15)(19) for various values of parameters and times: (a and b) the quasi-stationary-state solutions are the stable monotonic traveling waves; (c) the plane reaction front becomes unstable after an induction period of t ) 0.04; (d) the plane front is unstable from the beginning of its propagation. The parameters are σ ) 5, D ) 0.11, κ ) 370.4, with (a) γ ) 0, yc ) 0, (b) γ ) 0, yc ) 0.1, (c) γ ) 5, yc ) 0.1, and (d) γ ) 5, yc ) 0.14.

the concentrations of A and B tend to an oscillatory attractor around their QSS values upon the onset of instability of the stationary solution. We employ the QSS hypothesis to analyze the growth modes and the transition from frontal propagation to quasi-periodic deposition. If the QSS hypothesis holds in the time interval [t1, t2], then ω(t) changes slowly with time (see Figure 3 after an initial transient) and the solution of (15)-(19) y2(x,t) and z(x,t) at each instant is asymptotically close to the solution Y2 and Z of eqs 2022 with the function F given by y1, y2, and z at this instant (taking into account the translation of Y2 and Z in x). The function F must satisfy some conditions in order for the asymptotic correlations y2(0,t) ≈ y-, z(0,t) ≈ z- and y2(1,t) ≈ 1, z(1,t) ≈ 0 to hold at the ends of the domain 0 e x e 1. Simulations indicate that they hold when the kinetic rate constants are large enough. Thus, below we theoretically focus on the case where κ is large, i.e., diffusion is slow compared to reaction, and introduce a small parameter, µ ) 1/κ, to perform asymptotics. Theoretical Analysis for the Necessity of Autocatalysis and Critical Supersaturation. We start with the case where the autocatalytic step (14) is absent (γ ) 0). Then we can consider eqs 15 and 16 only, and the quasi-stationary solution is the traveling wave solution of eq 20 subject to the boundary conditions (22). It is known32-34 that this solution exists, is unique (except for translation), and is stable under some conditions on F(y2,t); in particular, for each time t ∈ [t1, t2], F must be smooth enough and satisfy

F(0,t) ) 0, F(y2,t) ) 0 at 1 - y0 e y2 e 1 (23) F(y2,t) > 0 at 0 < y2 < 1 - y0, 0 < y0 < 1 (24) In the absence of a critical supersaturation, i.e., when yc ) 0, the step function H(θ) ≡ 1 for y1 g 0, y2 g 0, and

F(y2,t) ≡ φ(y2) ) y1(y2) y2

Figure 3. Front positions (a) and front velocities (b and c) for simulations depicted in Figure 2. The letters a-d indicate the corresponding panels in Figure 2.

This system can be solved with initial functions y02 ) y2(x,t) and z0 ) z(x,t) that are the solution of the problem (15)-(19) at the same t. The time variable is now τ. Simulations reveal that the solutions of eqs 20 and 21 are the same pair of functions y2(x,t) and z(x,t), which move along the x axis with velocity ω(t) as τ increases. Thus, the solutions of eqs 20 and 21, y2(x,t) and z(x,t), may be thought of as quasi-stationary states (QSSs) in the transition from an initial to a final state in the problem (15)-(19). This situation is depicted in Figure 3, which shows the front location and its speed for the conditions of Figure 2. The growth mode predicted using the skeleton

(25)

It follows from comparison principles32 that the solutions y1(x,t) and y2(x,t) of the problem (15)-(19) are nonnegative, monotonic with respect to x, and bounded functions. As mentioned above at t ) 0, y1(x,t) ) 0 and y2(x,t) ) 1 at 0 e x e 1. Because of the continuous dependence on t, the solution of eqs 15-19 satisfies 1 - y0 e y2 e 1, 0 < y0 , 1, and y1 , 1 at 1 - x0 e x e 1 and 0 < x0 , 1 for some time interval [0, t]. Moreover, simulations show that as µ f 0 for 0 e x e xj - k, xj + k e x e 1, and k > 0 the solution of eqs 15-19 tends to the functions yj1(x,t) and yj2(x,t), which satisfy yj1yj2 ) 0. Here, x ) xj defines the coordinate of the internal transition layer where the solution of eqs 15 and 16 moves from one branch to another. Thus, the asymptotics y1(1,t) ≈ 0 and y2(1,t) ≈ 1 take place for a large enough time interval [t1, t2]. Consequently, conditions (23) and (24) hold, and there exists a stable quasistationary solution in the form of a traveling wave with y- ) 0 and z- ) 1. Thus, in the absence of autocatalysis (γ ) 0) and of a critical supersaturation (yc ) 0), frontal deposition takes place, as the example in Figure 2a


Figure 4. Form of function φ of eq 25.

shows. This conclusion for the need of a critical supersaturation is consistent with Oswald’s theory and is an expected result. The above conclusion can be extended to the case where there is a critical supersaturation but no autocatalysis; i.e., when yc > 0, γ ) 0, as the example in Figure 2b shows. In this case, we approximate the step function (19) by some monotonic, smooth function H h (θ) such that 0 < H h (θ) < 1 at -δ1 < θ < 0, H h (θ) ) 1 at θ g0, and H h (θ) ) 0 at θ e -δ1 where 0 < δ1 , 1. It follows from simulations that the function φ has the form shown in Figure 4. One can then see that for some interval [t1, t2], the function φ(y2) from eq 25 has two intersection points with the line φ ) yc - δ1 for y2 ∈ [0, 1], for some values of yc - δ1. Denote these points as y+ and y-. If y+ > y-, then φ′y2(y-) > 0 and φ′y2(y+) < 0. Evidently, y2 ) y- is the stable rest point of eq 20 without the diffusion term. Then the boundary conditions for the traveling wave solution are

y2 ) y- and z ) 1 - y- at x f -∞; y2 ) 1 and z ) 0 at x f ∞ With these boundary conditions, the traveling wave is stable; i.e., frontal propagation occurs. The necessity of criticality in terms of the existence of an autocatalytic step reconciles previously apparently confusing results and is qualitatively consistent with the excellent remark of ref 12 about the sensitivity of patterns on the form of nucleation. In particular, it explains why the numerical simulations in ref 9 did not result in LR bands, whereas the simulations in refs 7 and 12-15 resulted in banded deposits. The model of ref 9 consisted of eqs 15-17, but the sink term f was a function of the nucleation rate and the radius at time t of a particle that was nucleated in the past at time t′ < t. It was assumed that the nucleation and radial growth rates of the precipitate are functions of the supersaturation θ (eq 19) only. So, the simulated frontal regime of precipitation in ref 9 is just a consequence of the choice of the form of the function f and its parameters. The model of ref 12 resulted in LR after modification of the function f of ref 9. Independently, the same modification of the function f was done in ref 7 for the model of type (4) and (5) without taking into account the change of porosity . Also the model of ref 15 has the form (4), (5), and (10)-(12) for a constant porosity but a term φD instead of φD2/3 in W3. With some assumptions, the models of refs 13 and 14 can be written in the form (4) and (5), with the nonlinear terms retaining the characteristics of the functions W2 (eq 11) and W3 (eq 12). All of these models result in LR with the appropriate choice of parameters.

Note that, although the authors of refs 4-6 did not consider kinetic functions, they formulated boundary conditions by taking into account the aftereffects of the fast autocatalytic step. Their approach may be treated as some kind of narrow reaction zone approximation, which we will discuss below. This approximation is valid when both the nucleation and growth steps are fast enough. The work of refs 4-6 did not consider how and when the continuous deposition front loses its stability. It was assumed that a quasi-periodic pattern exists, and the diffusion equations (15)-(16) alone were solved without the reaction term for time intervals between two subsequent bands, which were assumed to form instantaneously. The boundary conditions corresponded to fast autocatalytic step when y2 ) 0 at the position of the precedent band formation. The spacing (and time) laws are obtained using the conditions y1y2 ) yc and ∂φ/ ∂x ) 0 at the position of the second band formation. Development of a Criterion for the Onset of LR Patterns. Next we consider the more general case where both autocatalysis and nucleation occur, i.e., yc > 0 and γ > 0. As mentioned above, simulations such as the ones shown in Figure 2c,d indicate that for some parameter values a wavelike solution becomes unstable, giving rise to spatial patterns. It is interesting to note an analogy with combustion of condensed systems. Experiments have shown that under certain conditions the velocity of flame propagation exhibits periodic pulsation and the burned samples have a layered structure normal to the flame front. The existence of pulsating combustion fronts was analyzed theoretically in ref 35, and the numerical solution of the corresponding boundary-value problem was obtained in ref 36. The model of ref 36 was analyzed in ref 37 using the narrow reaction zone approximation. It was shown that when the function F(y2,z,t) is replaced by a δ function of an appropriate strength, the uniformly propagating plane front loses stability when some parameter exceeds the critical value of φc ) 2 + x5. The analysis in ref 37 shows that a Hopf bifurcation occurs, giving rise to a stable time periodic solution in the ξ domain (see below) that leads to spatial pattern in the x domain. Formally, the model of ref 36 has the form of eqs 20 and 21, with an important difference regarding the function F. In our case, F is not an explicit function of y2 but is the result of the solution of the skeleton model (15)-(19) with corresponding boundary conditions at each time. However, we can still analyze the effect of parameters on the stability of the traveling wave solution of eqs 20 and 21. We apply the result of ref 37. In the narrow reaction zone approximation,33,38 the chemical reaction is taken to occur at the front, say at ξ ) 0, where ξ ) x + Ω(τ,t) is the moving coordinate system and x ) -Ω(τ) is the coordinate of the non-steady moving reaction front. For simplicity in notation, below we leave out the t dependence of y2. Furthermore, σF(y2,z,t) is approximated with ψ(y2) δ(ξ), where δ(ξ) denotes the Dirac δ function. Equations 20 and 21 may then be solved separately in each region (ξ > 0 and ξ < 0), subject to the corresponding jump conditions across the front located at ξ ) 0

∂2y2 ∂z ∂y2 ∂y2 ∂z + Ωτ ) D 2, + Ωτ ) 0, ξ * 0 (26) ∂τ ∂ξ ∂τ ∂ξ ∂ξ


ξ ) 0: [y2] ) 0, D

[ ]

∂y2 - ψ[y2(0,τ)] ) 0, ∂ξ Ωτ[z] ) ψ[y2(0,τ)] (27)

where [a] ≡ a(0+) - a(0-) and Ωτ ) ∂Ω(τ,t)/∂τ. To obtain an expression for ψ, we employ asymptotic analysis with respect to µ ) 1/κ, similarly to ref 39. According to eq 18, F can be written as F(y2,z) ) F1(y2,z) + F2(y2,z), where F1 ) κφ(y2) H h (θ), θ ) φ(y2) - yc, and F2 ) κγφ(y2) z. The narrow reaction zone implies that in front of the reaction front (ξ > 0) nucleation happens (F ) F1), whereas behind it (ξ < 0) only autocatalytic growth occurs, i.e., F ) F2. The front coordinate is defined by y2 ) y-, where y- is the smallest solution of

φ(y2) ) yc - δ1

(28)

We consider the problem (20)-(22) in a moving coordinate system, introduce a stretched coordinate (ξ ) µχ) inside the reaction zone, and obtain the outer and inner expansions of the functions y2(ξ,τ), z(ξ,τ) and ω(t), Ω(τ,t), respectively, with respect to the small parameter µ [Ω(τ,t) ) -ω(t) τ for a steadily propagating reaction front]. When the inner expansions are substituted into eqs 20 and 21, the nonlinear terms F1 and F2 are expanded into powers of µ, and the inner expansions are matched with the outer solutions, the lowest order approximation is (see appendix B)

ψ[y2(0,τ)] ) {2σDycκ[γ[y- - y2(0,τ)] + ∆y]}1/2 ψ(0) ) [2σDycκ(γy- + ∆y)]1/2 ) ω

(29) (30)

where ∆y ) y+ - y- and y+ and y- (y+ > y-) are the solutions of eq 28. Furthermore, we introduce the small perturbation quantity v(ξ,τ) according to y2(ξ,τ) ) y02(ξ) + v(ξ,τ), where y02(ξ) is the steady solution of eqs 26 and 27.

{

-(ω/D)ξ ξ > 0 y02(ξ) ) 1 - e 0 ξ φc ) 2 + x5

(32)

where φ ) -ψ′y2[y02(0)]/ω. Using the last equation for φ and eqs 29 and 30 from the asymptotics of our problem, we easily obtain the form of the stability criterion parameter

φ)

σDγycκ ω2(t)

(33)

When reaction (13) is considered to be fast, the dimensional front velocity, u(t′) ) ω(t) D0A/L, is u(t′) ) λ

x0D0A/t′,10 where0 λ is0 some nonlinear function of D ) DB/ DA and σ ) CA/CB. Then one can write eq 33 in the form

φ)

Ksk2DB 2

u (t′)

)

Ksk2DB t′ λ2D0A

(34)

Equation 34 shows an explicit dependence on t′. With

an increase in time, the width of the reaction zone increases as well, so the applicability of eq 34 is limited by the assumption of a narrow reaction zone. Equation 33, our main theoretical result, in conjunction with eq 32, determines the condition of stability of frontal deposition. For parameter combinations resulting in large values of φ, the wave becomes unstable, giving rise to pattern formation. On the other hand, when φ is sufficiently small, frontal propagation occurs for all times. Fast autocatalytic kinetics and large values of critical supersaturation are prerequisites for band formation for large κ and strong concentration gradients. Several simulations conducted (not shown) for various combinations of dimensionless parameters in ref 30 showed that eq 33 holds nearly quantitatively, providing numerical validity to our theoretical result. The time-varying speed has important ramifications for explaining simulation and certain experimental results. From eqs 15-19 and the corresponding boundary conditions, it becomes apparent that, as time t increases, the quasi-stationary velocity ω(t) decreases because of the diffusion time needed for species A to reach the moving front. Thus, there exists an induction time t* and a corresponding distance from the boundary x* where ω(t*) becomes sufficiently small and eq 32 becomes satisfied. Then the traveling wave becomes unstable, and patterns emerge. The change of the speed of the traveling wave with time explains then the presence of an induction time and length. Estimation of the Location for the Onset of Bands. Because of the nonanalytical form of function F, one can only estimate when the transition from frontal propagation to quasi-periodic deposition occurs. The critical conditions can be estimated from simulation or experiment in the following way using the relations x* ) 2λxt* and ω* ) λ/xt*. Here x*, ω*, and t* denote the location of the front, its speed, and the critical time when the front loses stability. The values of the critical velocities are obtained from eq 33, and the values of λ are estimated using graphs such as Figure 3a obtained from simulation or experiment, i.e., by estimating the front location versus time (λ ≈ 0.8 for Figure 5b,c and λ ≈ 0.9 for Figure 5a). After that, it is possible to estimate xt* ) λ/ω* and finally obtain x*. This idea could be repeated for a few reagent concentrations, using factorial design methods, to develop a low-degree polynomial model that can capture the location where the instability happens versus boundary concentrations. This simple model can then be used to guide experiments for controlling the onset of bands. Alternatively, a full transport-chemistry model, which has been validated against continuous deposition experiments, could be used to guide more challenging experiments aimed at control of bands. Figure 5 along with the so-estimated critical locations, summarized in its caption, is in reasonable agreement with the simulation results about when and where instability occurs. Thus, the proposed approach provides a reasonable way to estimate when bands start and a way to control the location of bands within a substrate by mainly changing the concentrations of the reagents. Comparison to Combustion Processes and Suggestions for Deposition Experiments There exists an interesting difference of gasless autooscillatory combustion and pattern formation discussed


Figure 6. Effect of the front velocity on pattern periodicity using the skeleton model: (a and c) the front velocity is determined from the solution of eqs 15-19 and (b and d) the “convection” term, VA (∂y1/∂x), has been added to the right-hand side of eq 15. “Convection” leads to a slower decrease of the front velocity and more regular patterns. The parameters are σ ) 5, κ ) 370.4, D ) 0.11, and t ) 0.075. For parts a and b, γ ) 10, yc ) 0.15, and ν ) 1, and for parts c and d, γ ) 50, yc ) 0.1, and ν ) 2/3. Figure 5. Dynamics of deposition indicating the induction time and location for the transition from frontal propagation to quasiperiodic growth mode using the skeleton model. The parameters are σ ) 5, κ ) 370.4, and D ) 0.11 with (a) γ ) 10, yc ) 0.075, t* ) 0.022, ω* ) 6.0, x* ) 0.27, (b) γ ) 10, yc ) 0.15, t* ) 0.009, ω* ) 8.5, x* ) 0.15, and (c) γ ) 20, yc ) 0.15, t* ) 0.004, ω* ) 12.0, x* ) 0.11. The vertical arrow in each panel indicates the estimated value of x*. Here the asterisk denotes a critical condition.

here related with the periodicity. The combustion experiments and corresponding models exhibit selfsustained oscillations, whereas deposition obeys usually uneven spacing and time laws. The reason is that the source function F(y2,z,t) and the instant quasi-stationary front velocity ω(t) are functions of time. Moreover, the decrease of ω(t) with t leads to a variable period of oscillations, i.e., the uneven spacing-time laws of LR bands. Note that in combustion models, upon an increase of the bifurcation parameter, period doubling and transition to chaos occur.40 Finally, we note that there exists one more similarity with combustion processes. In combustion of condensed systems, spin regimes have been found to consist of a helical motion of a luminous spot along the side surface of the specimen. The conditions of bifurcation of spin waves from plane waves were formulated in ref 35. A helical deposition structure41 was described in ref 13. A spiral solution was found when the radius of the test tube is large. It may be possible that the appearance of such growth modes may be treated as another bifurcation from a plane reaction front when a bifurcation parameter, connected with the test tube radius, exceeds a critical value. From an experimental standpoint, it is interesting to know whether bands are connected or not to each other and what controls band structure and spacing. We discuss here some possible ideas. It is expected that, for a relatively time-independent source F, regular periodic patterns would develop. This may be achieved in different ways. For example, the introduction of “convec-

tive” terms with constant velocities VA * VB into the governing equations of the prototype and skeleton models for reactants A and B leads to more regular deposition patterns. An example is depicted in Figure 6. A similar effect was described in ref 10, where an electrical field was applied to an interdiffusion precipitation system. Analysis of the effect of convection on a propagating front is a rather complicated subject (e.g., see ref 42) and deserves a separate in-depth study. Simulations show that the space between bands is either empty or filled with less dense deposit. Some examples are depicted in Figure 7. The empty space between bands results from an increase of the parameter φ, as shown in Figure 7a,b, and is associated with an increase in the amplitude of the spatial patterns. In general, the structure of bands depends on the nucleation kinetics (e.g., the specific form of the step function H) and the autocatalytic reaction rate function F2. For example, we have extended the work above to consider the function

F2 ) κγφ(y2) zν

(35)

An approximate estimate revealed that the bifurcation parameter φ is proportional to 1/ν. Simulations show that for ν < 1 the banded deposit is more regular than for νg1, as shown in Figure 7c,d. Conclusions In this paper we have presented a prototype nucleation and growth model for LR pattern formation in porous media. Simulations indicate two growth modes, a frontal propagation where a traveling wave leaves behind a continuous solid film and a quasi-periodic band formation mode where bands of uneven spacing are formed. An induction time for the transition from frontal propagation to band formation has often been observed. A skeleton, simpler model for LR pattern formation was


The concentration z has a jump at ξ ) 0 (where the reaction is concentrated) [z] ) z+ - z- ) -1. It follows from eq 27 that

Ωτ ) -ψ[y2(0,τ)]

(A.3)

and for y02(ξ) (eq 31)

ω ) ψ(0)

(A.4)

We introduce the small perturbation quantity v(ξ,τ) according to

y2(ξ,τ) ) y02(ξ) + v(ξ,τ)

(A.5)

Substituting eq A.5 into eqs 26 and 27 and taking into account eqs A.3 and A.4 and that ψ[y02(0) + v(0,τ)] ) ψ[ y02(0)] + ψ′y2[y02(0)]v(0,τ) + Ο(v2), we obtain 0

Figure 7. Effect of parameters on the structure of patterns using the skeleton model. Panels a and b show that an increase in the solubility product with a concomitant increase in the bifurcation parameter φ makes the space between bands more empty. Panels c and d show that the space between bands becomes empty with a decrease of the stoichiometric coefficient ν of the autocatalytic reaction. The parameters are σ ) 5, κ ) 370.4, D ) 0.11, and t ) 0.15 with (a) γ ) 6, yc ) 0.1, ν ) 1, (b) γ ) 6, yc ) 0.15, ν ) 1, (c) γ ) 7, yc ) 0.1, ν ) 2/3, and (d) γ ) 7, yc ) 0.1, ν ) 1.

subsequently introduced that exhibits qualitatively similar solutions under many conditions. By assuming quasi-steady-state behavior for the traveling deposition waves, bifurcation analysis could be applied to the skeleton model and a theoretical criterion was derived to explain how model parameters control the onset of instability of waves, giving rise to band formation. Our simulations and theory explain some apparent discrepancies of the previous work. Analogies to classic oscillatory behavior of autocatalytic reactions and spatial pattern formation in combustion have been drawn. We have attributed the uneven spacing law to the variation of the speed of the traveling wave with time and suggested that the application of an external field, such as convective flow, can diminish such variations and provide almost periodic structures. Finally, our theory holds only for the case of strong concentration gradients. Future work will focus on bifurcations of other cases and comparison to experiments.

∂y2 ∂v ∂v ∂2v -ω + Rv(0,τ) ) D 2, ξ * 0 ∂τ ∂ξ ∂ξ ∂ξ and

[∂v∂ξ] + Rv(0,τ) ) 0

ξ ) 0: [v] ) 0, D

v ) 0 at ξ f (∞ where

R ) -ψ′y2[y02(0)] ) -ψ′y2(0)

(A.6)

The normal-mode solution has the form eiλ(ω2/D)τθ(ξ), where

{()

θ(ξ) )

for ξ < 0 θ(0) e-(ω/D)lξ R θ(0) -(ω/D)ξ R -(ω/D)qξ for ξ > 0 e + θ(0) 1 + e ω iλ iλω

(

)

with

1 1 l ) (1 - x1 + 4iλ), q ) (1 + x1 + 4iλ) 2 2 λ ) λ(R/ω) is determined from

Acknowledgment This work was supported in part by the National Science Foundation. Appendix A: Instability Criterion We consider the linear stability analysis of the timeindependent solution y02(ξ) of eqs 26 and 27 given by eq 31. The boundary conditions are

y2 ) 0 and z ) 1 at ξ f -∞; y2 ) 1 and z ) 0 at ξ f ∞ (A.1) This solution corresponds to the uniformly propagating plane front with

Ωτ ) -ω

(A.2)

[(ωR) - 4(ωR) - 1] + i(ωR) ) 0

4λ2 - iλ

2

(A.7)

The stationary solution y02(ξ) is unstable if Reiλ > 0. Using eq A.7, this implies that the instability condition occurs when R/ω > 2 + x5. Using eq A.6, one obtains

-

ψ′y2[y02(0)] ω

> 2 + x5

(A.8)

Appendix B: Source Strength ψ(y2) It follows from eq A.3 that the source strength ψ is determined from the nonstationary front velocity Ωτ of the traveling wave solution of eqs 20 and 21. In terms of the transformed coordinate system ξ ) x + Ω(τ), eqs 20 and 21 become


∂y ∂y ∂2y + Ωτ ) D 2 - σF(y,z) ∂τ ∂ξ ∂ξ ∂z ∂z + Ωτ ) σF(y,z) ∂τ ∂ξ

(B.1) (B.2)

where y ≡ y2. First we consider the asymptotic solution for steady propagation when eq A.2 holds, and the governing equations are

-ω

d2y dy ) D 2 - σF(y,z) dξ dξ

(B.3)

dz ) σF(y,z) dξ

(B.4)

-ω

with boundary conditions given by eq A.1. In the limit of large κ

µ)

1 ,1 κ

We then seek expansion of the variables y and z and parameter ω as

(B.5)

the reaction zone is an asymptotically thin layer in the neighborhood of ξ ) 0. As assumed above, one has

{

y ) µy(1) + µ2y(2) + ...

(B.11)

z ) z(0) + µz(1) + ...

(B.12)

ω ) ω(0) + µω(1) + ...

(B.13)

Let the front coordinate be defined by the value χn-, where y ) y- and the nucleation begins. Nucleation occurs at χn- < χ < χn+ with y ) y+ at χ ) χn+. Autocatalytic growth occurs at χ < χn- where 0 < y < y-, and the concentration z takes some value zn at χ ) χn-. For the thin reaction zone

y- ) µv-, y+ ) µv+, v( ) O(1)

(B.14)

We substitute eqs B.11-B.13 into eqs B.3 and B.4 and expand the nonlinear terms in powers of µ to obtain

{

1 [φ0(y(1)) + O(µ)] at v- < y(1) < v+ µ F1 ) 0 at y(1) < v-, y(1) > v+ 1 F2 ) [γφ0(y(1)) z(0) + Ο(µ)] µ

F at ξ > 0 F(y,z) ) F1 2 at ξ < 0 where

with (see eqs 18 and 25)

φ0 ) lim φ(µy(1) + µ2y(2) + ..., µ)

φ(y,µ) H h (θ), θ ) φ(y,µ) - yc F1 ) µ

(B.6)

φ(y,µ) γz F2 ) µ

(B.7)

µf0

The first-order inner problem is

D

Recall that H h (θ) is a smooth approximation of the step function H (eq 19). Thus,

F1 ) 0 at y e y- and y g y+

subject to the boundary conditions (A.1). Thus, to all algebraic orders in the small parameter µ, the outer solution is

{

{

dz(0) ) σγφ0z(0) dχ

(B.17)

at χ < χn- and

D

d2y(1) - σφ0 ) 0 dχ2

-ω(0)

dz(0) ) σφ0 dχ

(B.18)

(B.19)

y(1) f 0 and z(0) f 1 at χ f -∞

(B.20)

y(1) ) v- and z(0) ) zn at χ ) χn-

(B.21)

y(1) ) v+ at χ ) χn+

(B.22)

(B.9)

We seek an asymptotic expansion for the steady front velocity ω by suitable matching of the inner reaction zone solution with the outer one (eq B.9). In the thin reaction zone, the independent variable ξ is stretched according to

ξ ) µχ

-ω(0)

(B.16)

at χn- < χ < χn+. Using eqs B.9-B.12 and B.14, we obtain the boundary conditions

d2y dz dy ) D 2, -ω ) 0, ξ * 0 dξ dξ dξ

-(ω/D)ξ ξ > 0 0 ξ>0 , z) y) 1-e 1 ξ

Pattern Formation in Porous Media via the Liesegang Ring

Recommend Documents