Spinodal decomposition of binary mixtures with composition-dependent heat conductivities

Chemical Engineering Science 63 (2008) 2402 – 2407 www.elsevier.com/locate/ces

Spinodal decomposition of binary mixtures with composition-dependent heat conductivities Dafne Molin, Roberto Mauri ∗ Dipartimento di Ingegneria Chimica, Chimica Industriale e Scienza dei Materiali, Università di Pisa, via Diotisalvi 2, I-56127 Pisa, Italy Received 29 May 2007; received in revised form 23 January 2008; accepted 26 January 2008 Available online 6 February 2008

Abstract Spinodal decomposition of a very viscous binary mixture is simulated, assuming that the fluid is bounded by two walls that are instantaneously quenched below the critical temperature. As usual, the mixture phase separates by forming dendritic structures whose size grows with time. While in conventional spinodal decomposition the morphology of the mixture is isotropic, here we show that, when the two components of the mixture have different heat conductivities, the configuration of the mixture may become anisotropic, depending on whether heat propagates slower or faster than mass. Specifically, for large heat conductivity differences, when the Lewis number is small, dendrites tend to align parallel to the isothermal lines; when it is large, dendrites progressively align perpendicular to the isothermal lines. This behavior can be explained observing that the morphology of the mixture is the result of two competing effects: on one hand, the mixture phase separates as soon as its temperature reaches its critical value and therefore it tends to follow the isothermal lines; on the other hand, the system tend to maximize heat transport, and therefore dendrites tend to span between the two quenched walls, i.e. perpendicular to the isothermal lines. This explanation is confirmed showing that when the imposed temperature difference between the side walls (and therefore the heat flux at equilibrium as well) is increased, dendrites tend to align across the walls more and more. 䉷 2008 Elsevier Ltd. All rights reserved. Keywords: Phase change; Separation; Transport processes; Simulation

1. Introduction The objective of this work is to study numerically how the morphology of a very viscous binary system evolves after its walls are instantaneous quenched below the critical temperature. In conventional spinodal decomposition, the mixture is quenched instantaneously to a temperature that is uniform everywhere and, therefore, it phase separates by forming isotropic bicontinuous structures (for a review on spinodal decomposition, see Gunton et al., 1983). A different result is obtained when we simulate spinodal decomposition in a slab whose walls are quenched instantaneously to the same temperature. In that case, as shown by Vladimirova et al. (1998), the morphology of the system during phase separation consists of dendrites aligned along the longitudinal direction (i.e. parallel to the quenched ∗ Corresponding author. Tel.: +39 050 511248; fax: +39 050 511266.

E-mail addresses: [email protected] (D. Molin), [email protected] (R. Mauri). 0009-2509/$ - see front matter 䉷 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2008.01.028

walls), forming typical striped pattern as in the experiments of Sagui and Desai (1994). In the present work, instead, we simulate the phase separation process of a mixture that is confined in a slab whose walls are quenched to two different temperatures, assuming, in addition, that the heat conductivities of the two phases are different from each other. In that case, since the system will tend to a steady state where entropy production, and therefore heat transfer, will be maximized, the morphology of the mixture should get as close as possible to that corresponding to straight tubular domains spanning the whole width of the channel. This effect should become more relevant as the temperature difference between the walls and the heat conductivity ratio are increased. To make this point clearer, we chose to consider systems where convection is negligible, so that heat and mass transport are exclusively diffusive, as it happens in polymer solutions, like a polybutadiene/polyisoprene mixture and alloys, like a basaltic melt. We believe, however, that the general features that are explored here are also valid for low viscosity mixtures.

D. Molin, R. Mauri / Chemical Engineering Science 63 (2008) 2402 – 2407

The binary mixture is described using the diffuse interface model, that is based on the pioneering work by Van der Waals (1979) at the end of the 19th century, then extended by Landau and Lifshitz (1980) by developing the mean field theory and finally was applied to binary mixtures by Cahn and Hilliard (1958) at the end of the 1950s. The model shows that during the early stages of the phase separation, initial instabilities grow exponentially, forming single-phase micro-domains whose size corresponds to the fastest-growing mode of the linear regime. During the late stage of the process, the system consists of well-defined patches in which the average concentration is not too far from its equilibrium value. At this point, both analytical calculations (Lifshitz and Pitaevsky, 1984) and dimensional analysis (Siggia, 1979) predict a growth law of the typical domain size with time, R(t) ∼ t 1/3 , due to the Brownian coagulation of droplets. Naturally, while, when the mixture is unbounded, the process of phase separation is isotropic, when the mixture is confined between two walls, the single phase domains will form anisotropic patterns. In fact, this is the subject of the present work. The paper is divided as follows: in Section 2 we derive the governing equations using the diffuse interface model, while in Section 3 we present and discuss the numerical results. 2. The governing equations 2.1. The diffuse interface approach Consider a homogeneous mixture of two species A and B with molar fractions xA =  and xB = 1 − , respectively, which are kept at temperature T and pressure P. For simplicity, we assume that the molecular weights, specific volumes and viscosities are constant and that they are the same for the two species, namely, MA = MB = Mw , VA = VB = V , A = B = ; in particular, that means that we assume the mixture to be incompressible. The thermodynamic molar Gibbs free energy, gth is a coarse-grained free energy functional given by gth () = gid () + gex (),

(1)

where gid is the Gibbs free energy of an ideal mixture, gid () = RT[ ln  + (1 − ) ln(1 − )],

(2)

where R is the gas constant, while gex is the excess part of the free energy, gex () = RT(1 − ),

(3)

where  is a function of T (it cannot depend on P, since the mixture is assumed to be incompressible). This expression, generally referred to as the one-parameter Margules correlation, can been derived considering the intermolecular potentials between identical and different neighboring particles (Sandler, 1999; Mauri et al., 1996). The ideal part (1) of the free energy is proportional to the temperature and represents an entropic contribution, while the excess part is an enthalpic term and depends on inter-particle potentials. In fact, assuming that our mixture is regular, its entropy equals that of an ideal mixture

2403

(Sandler, 1999), i.e. gex does not depend on T, so that  ∝ 1/T . Accordingly, considering that the free energy has an inflection point at the critical point, i.e. when T = TC and  = C = 1/2, then we see that C = 2, and therefore the Margules parameter can be expressed as 2TC . (4) T In order to take into account spatial inhomogeneities, Cahn and Hilliard (1958) introduced the generalized specific molar free energy g, given by the expression

=

g(, ∇) = gth () + gNL (∇),

(5)

where gNL (∇) = 21 a 2 (∇)2 ,

(6)

is a nonlocal molar free energy due to changes in composition, where a represents the typical length of spatial inhomogeneities. As shown in Molin and Mauri (2007), a depends on the temperature as  √ 2TC a = aˆ (7) = aˆ , T where aˆ is a characteristic length, independent of the temperature. As shown by Van der Waals (1979), as the surface tension  is the energy stored in the unit area of the interface separating two phases at local equilibrium, we obtain  1 RT 2 RTa = (∇)2 dl ∼ a , (8) 2 Mw Mw where we have considered that far from the critical point, the molar fraction difference between the two phases at equilibrium, ()eq , and the nondimensional thickness of the interface √ region at local equilibrium, /a ≈  − 2, are both of O(1). This shows that the characteristic length a can be determined once surface tension is known (see discussion in Vladimirova et al., 2000). From the molar free energy, we may define a generalized chemical potential difference ˜ as ˜ =

(g/RT) j(g/RT) j(g/RT) = −∇ ·  j j(∇)  = 0 + ln + (1 − 2) − a 2 ∇ 2 , 1−

(9)

where 0 = (gB − gA )/RT. 2.2. The equations of motion Imposing that the number of particles of each species is conserved and that the mixture is incompressible, we obtain the continuity equations (Vladimirova et al., 2000), d j = + v · ∇ = −∇ · j , dt jt

(10)

where d/dt = j/jt + v · ∇ is the material derivative, v is the average velocity of the mixture, while j is the diffusive mass

2404

D. Molin, R. Mauri / Chemical Engineering Science 63 (2008) 2402 – 2407

flux. The latter is proportional to the gradient of the generalized chemical potential gradient as j = −(1 − )D∇, ˜

(11)

where D is a composition-independent diffusion coefficient. Finally, substituting Eq. (9) into Eq. (11), we obtain j = − D∇ + D(1 − )[a 2 ∇∇ 2  + 2∇ + (2 − 1)∇].

(12)

Here, the term D∇ on the RHS represents the regular diffusion flux, while the last term vanishes for small concentrations of either solvents, i.e. when ( → 0 or 1), and for ideal mixtures, i.e. when a =  = 0. Note that the a 2 term is always stabilizing, with a being a function of temperature through Eq. (7). Eq. (10) must be coupled to the traditional heat equation, c

dT ˙ + ∇ · jq = q, dt

(13)

where  is the mass density and c the specific heat (which are both assumed to be uniform, here), jq is the diffusive heat flux, while q˙ is the energy dissipation term. Assuming a simple Fourier constitutive relation for the heat flux, we have jq = −k()∇T ,

(14)

where k() is the composition-dependent heat conductivity. Assuming a linear dependence between heat conductivity and composition, with kA = k and kB = k denoting the heat conductivities of the two species, we obtain k() = k[ + (1 − )].

(15)

Since the problem is invariant upon changing  and  into (1−) and −1 , here we will assume, without loss of generality, that 0 < 1. In general, Eqs. (10) and (13) should be coupled to the equation of momentum conservation. This equation here does not play any role, though, as we assume that heat and molar convection fluxes are negligible compared to diffusion fluxes, that is both molar and thermal Peclet numbers are small, i.e. NP eM = V a/D>1 ˆ and NP eT = V a/ >1, ˆ where = k/c is the heat diffusivity of species A and V is a characteristic ve−1 locity. Note that NP eT = NP eM NLe , where NLe = /D is the Lewis number, which for most polymer melts ranges from 104 to 106 . Accordingly, NP eM >1 is the leading assumption here, which allows us to state that v = 0, i.e. d/dt = j/jt, in Eqs. (10) and (13). As shown in previous works (see for example Vladimirova et al., 1999), V ≈ (aRT)/(M ˆ w ), so that this assumption can also be written as (aˆ 2 RT)/(Mw D)>1, which is satisfied for very viscous mixtures. The energy dissipation term in Eq. (13) is the sum of viscous dissipation and heat of mixing (Molin and Mauri, 2007),    ˙ q˙ = q˙v + h˙ mix = |∇v|2 + hex , (16) Mw where hex is the energy that is generated when one mole of mixture is mixed together isothermally. Now, since

from classical thermodynamics (Sandler, 1999) hex = −RT2 [j(gex /RT)/jT ]P , , we see that, for regular mixtures, when gex does not depend on T, we have hex = gex = RT(1 − ).

(17)

From these expressions, we see that q˙v = O(V 2 /aˆ 2 ), while h˙ mix =O(RTD/Mw aˆ 2 ). Therefore, defining the ratio between viscous dissipation and heat of mixing as a modified Eckert number, i.e. NEc ≡ q˙v /h˙ mix ≈ (V 2 MW )/(RTD), here we assume that NEc >1, which is satisfied in polymer melts, so that q˙ ≈ h˙ mix . Now, comparing in Eq. (13) the diffusive term with the heat generation term, we see that O(∇ · jq /q) ˙ = NLe c, ˆ where cˆ = Mw c/R is the nondimensional specific heat. Therefore, since for polymer melts this ratio is ?1, here we will neglect the heat generation term altogether, so that the heat equation becomes dT = (1 − )∇ · ∇T + [ + (1 − )]∇ 2 T . dt

(18)

At this point, scaling Eqs. (10) and (18) in terms of rˆ =

1 r aˆ

and

t˜ =

D t aˆ 2

(19)

we obtain the following equations governing the concentration field, , and the Margules parameter  (which is inversely dependent on the temperature T through Eq. (4)): j ˜ · {∇ ˜ − (1 − )[(2 + ∇˜2 )∇ ˜ =∇ jt˜ ˜ + (2 − 1)∇]} and −1 j NLe ˜



2 ˜ 2 = [ + (1 − )] ∇˜ 2  − (∇)  jt ˜ · ∇. ˜ + (1 − )∇

(20) 

(21)

In these equations, we remind, all convective effects have been neglected, which means that we have assumed: (a) both molar and thermal Peclet numbers are small, i.e. NP eM >1 and NP eT >1; (b) the modified Eckert number is small, >1. In addition, heat generation has been neglected as NEc ˆ For example, in a typical polywell, assuming that NLe c?1. mer melt, like a polybutadiene/polyisoprene mixture, we have aˆ ≈ 10−6 cm, MW ≈ 102 g mol−1 ,  ≈ 104 g cm−1 s−1 and D ≈ 10−6 cm2 s−1 , so that we find NP eM ≈ 10−2 ?NP eT , ≈ 10−2 , while the last condition, N c?1, NEc is obviously Le ˆ amply satisfied. 3. Numerical results The governing equations (20) and (21) were solved on a uniform two-dimensional square grid with constant width ((x (i) , y (j ) ) = (ix, j y), i = 1, . . . , Nx , j = 1, . . . , Ny ) and time discretization (t (n) = nt, n = 0, 1, 2, . . .). The physical dimensions of the grid and the time step were chosen such that x = y = 2a, while t/(a 2 /D) ≈ 0.01.0.001. A convergence test with respect to grid refinement was carried out as in

D. Molin, R. Mauri / Chemical Engineering Science 63 (2008) 2402 – 2407

Lamorgese and Mauri (2005), showing that such grid size is appropriate. The time step was determined semi-empirically, in order to maintain the stability of the numerical scheme. Note that the nonlinearity of the equations prevents a rigorous derivation of the stability constraints on t, but one can roughly estimate that the size of t will scale as O(x 4 , y 4 ), which is the order of the highest operator in the discretized system (see discussions in Lamorgese and Mauri, 2005). The space discretization was based on a cell-centered approximation of both the concentration and the inverse temperature variables,  and . The spatial derivatives on the RHS of Eqs. (20) and (21) were discretized using a straightforward second-order-accurate approximation. The time integration from t (n) to t (n+1) was achieved in two steps. First, we solved Eq. (21) to determine the Margules parameter  at time t (n+1) using (t (n) ) and (t (n) ); then Eq. (20) was solved, determining the concentration field (t (n+1) ) using (t (n) ) and (t (n+1) ). This makes the entire scheme stable both in time and in space. Eqs. (20) and (21) were solved with periodic boundary conditions along the y-axis, while at the other boundaries, i.e. the walls along the x-axis, we assumed that the material flux is zero, while the temperature is imposed. Initially, at time t < 0, the concentration field is uniform, with =0.5, while the Margules parameter is  = 1.9, corresponding to a temperature above its critical value. A random noise is superimposed to the composition field, with   =0 and ( )2 1/2 =10−2 , assuming that it is uncorrelated both in space and in time (Lamorgese and Mauri, 2006). Our results appear to be identical to those of another set of simulations, where the random noise was introduced initially and then removed. In fact, as shown in Vladimirova et al. (1998), the strength of the random noise determines a delay, after which spinodal decomposition starts to occur, independent

2405

of the random noise strength, provided that the Ginzburg criterion is satisfied (see Lamorgese and Mauri, 2006). That means assuming that the mean square concentration fluctuations are much smaller than the average concentration variation; in our case, since we are far from the critical point, this criterion is satisfied. The simulations were carried out using different square grids, ranging from 128 × 128 to 1024 × 1024, to ensure that the results are independent of the domain size. Time was measured using a timescale 105 aˆ 2 /D, which is typically of the order of 1–10 s. The results depend on three parameters: (a) the Lewis number, NLe = /D, expressing the ratio between thermal and mass diffusivities; (b) the heat conductivity ratio,  = kA /kB ; (c) the quenching depth. First, we studied the effects of the first two parameters, i.e.  and NLe , while the quenching depth was kept constant. Accordingly, we assumed that at time t = 0 the two walls are cooled down to two different temperatures, Tw1 on the left side wall and Tw2 < Tw1 on the right side wall, corresponding to w1 =2.1 and w2 =3.0, both well within the two phase region of the phase diagram. This ensures the presence of a temperature gradient within the system, so that a heat flux is always present. In Fig. 1, we show some snapshots of our simulations for different values of the conductive ratio and the Lewis number, at different time steps. The cases with small Lewis number, although perhaps not very common in practice, have been included because, as will be shown, they help to understand the phenomenon. Comparing the results for  = 0.001 with those for =0.1 we see that, as expected, the smaller is , the smaller is the average heat conductivity of the mixture, and therefore the shorter it takes to phase separate. In addition we see that, keeping constant the heat conductivity ratio (here  = 0.001),

Fig. 1. Evolution of the concentration field for different values of the conductive ratio, , and the Lewis number, NLe . The snapshots are taken at different times, where t is expressed in 105 aˆ 2 /D units.

2406

D. Molin, R. Mauri / Chemical Engineering Science 63 (2008) 2402 – 2407

when NLe = 0.01 the dendritic structures that form during the phase separation tend to align along the y-axis, while when the Lewis number is increased to NLe =100 they tend to align along the x direction. In fact, the main result of our simulations is that the values of the parameters  and NLe strongly influence the anisotropy of the system morphology. The isotropy of any given discrete concentration field i,j , with i, j = 1, 2, . . . , N, can be quantified by the isotropy coefficient as follows: =

Nx − N y Nx + N y

(22)

with

N

N 1  Nx = (N − ni,j ni+1,j ) ; 2 j =1 i=1 ⎡ ⎤ N N   1 ⎣ (N − ni,j , ni,j +1 )⎦ , Ny = 2

(23)

Fig. 2. The isotropy coefficient at steady state for different values of the conductive ratio, , and the Lewis number, NLe .

where the index ni,j is defined as equal to +1 when i,j > 0.5 and equal to −1 when i,j < 0.5. Clearly, Nx and Ny indicate how many times, moving along the x- and the y-axis, respectively, we cross an interface. Accordingly, −1   + 1, with = 0 when the morphology is isotropic, = +1 when it is composed of straight lines along the y-axis, and = −1 when it is composed of straight lines along the x-axis. Note that the isotropy coefficient has such a simple definition because we know that the anisotropy can only be along the x- or the y-axis. At this point, it should be pointed out that the value of the isotropy coefficient fluctuates considerably, since in any typical run Nx and Ny are of O(10). In order to evaluate the extent of such uncertainty, we considered a large number of spinodal decomposition simulations with periodic boundary conditions that we have carried out in the past. As expected, we found that the corresponding values of fluctuate around = 0, due to the system isotropy, with mean square displacement = ±0.05. In Fig. 2, the value of the isotropy coefficient corresponding to the final, steady state configuration of the mixture is plotted as a function of  and NLe . The height of the error bars was determined based on the above-mentioned fluctuations. Again, without loss of generality, we have considered only the cases with  < 1: the cases with  > 1 can be derived from these = N /. results with  = 1/ and NLe Le First of all, we see that, as expected, for =1 the morphology is isotropic for any value of the Lewis number, as the two phases have the same heat conductivity. On the other hand, when the difference between the conductivity of the two phases increases, the morphology of the system strongly depends on the value of the Lewis number. Specifically, we see that when the Lewis number is very small, the dendrites tend to align along the y-axis, while, as the Lewis number increases, they turn their orientation toward the x-axis. In fact, when NLe >1, heat propagates much slower than mass; accordingly, as soon as the temperature of the slab reaches its critical value, the mixture phase separates immediately. Since the isothermal lines (not shown here) are in all cases almost parallel to the y-axis, it is

therefore not surprising that, in this case, the dendrites will be also aligned along the same direction. On the other hand, when NLe ?1, heat propagates much faster than mass; our results suggest that, in this case, the system tends to maximize the heat flux, forming lines of the more conducting species spanning between the two quenched walls, along the x-axis. Flux maximization can be seen as a consequence of the fact that the steady state reached by a system with given, constant constraints must maximize its entropy production. Studying a similar problem, Martiouchev and Seleznev (2000) showed that the morphologies selected during crystallization tend to follow the same rule, maximizing the entropy production. The law of maximum entropy production was recently demonstrated in great generality by Dewar (2005) using nonequilibrium statistical mechanics, showing that it applies to all cases where macroscopic fluxes are free to vary under imposed constraints. Despite its name, this law generalizes Onsager’s (1931) principle of least energy dissipation to cases where systems are far from equilibrium. This explanation was confirmed by the results obtained by varying the quenching depth. First, we studied the phase separation occurring when the two walls are quenched to the same temperature, i.e. when the mixture tends to an equilibrium state, with no entropy production. In this case, the system corresponds to the specular double of that obtained when the left wall is quenched, while the right wall is insulated. Accordingly, as shown in Fig. 3, the mixture tends to phase separate by forming y-oriented dendrites near the insulated wall (i.e. near the center of the channel), in agreement with Vladimirova et al. (1998). Therefore, we may conclude that it is only in the presence of a steady state heat flux, i.e. when the two walls are quenched to two different temperatures, that anisotropy might develop perpendicular to the isothermal lines, that is, in our case, along the x-axis. This happens when NLe ?1, so that the tendency of the system to maximize the heat flux can prevail over the competing tendency to phase change as soon as the temperature falls

i=1

j =1

D. Molin, R. Mauri / Chemical Engineering Science 63 (2008) 2402 – 2407

2407

Fig. 3. Evolution of the concentration field for NLe = 0.1 and  = 0.001 when the left wall is quenched, while the right wall is insulated.

below its critical value (therefore forming y-oriented dendrites). In fact, when, with =0.001 and NLe =100, we keep w1 =2.1 constant on the left side wall, while w2 on the right side wall is changed from 3.0 to 5.0 (corresponding to an unrealistically large temperature quench), then we see that at steady state the mean isotropy coefficient changes from −0.22 to −0.53. Acknowledgments This work was supported by Ministero Italiano dell’Università e della Ricerca (MIUR). References Cahn, J.W., Hilliard, J.E., 1958. Free energy of a nonuniform system. I. Interfacial free energy. Journal of Chemical Physics 28, 258–267. Dewar, R.C., 2005. Maximum entropy production and the fluctuation theorem. Journal of Phisics A 38, 371–381. Gunton, J.D., San Miguel, M., Sahni, P.S., 1983. In: Domb, C., Lebowitz, J.L. (Eds.), Phase Transition and Critical Phenomena, vol. 8. Academic Press, London. Lamorgese, A.G., Mauri, R., 2005. Nucleation and spinodal decomposition of liquid mixtures. Physics of Fluids 17, 034107/1-10. Lamorgese, A.G., Mauri, R., 2006. Mixing of macroscopically quiescent liquid mixtures. Physics of Fluids 18, 044107/1-11. Landau, L.D., Lifshitz, E.M., 1980. Statistical Physics, Part I. Pergamon Press, New York. (Chapters 142–153).

Lifshitz, E.M., Pitaevsky, L.P., 1984. Kinetic Physics. Pergamon Press, New York. (Chapter 100). Martiouchev, L.M., Seleznev, V.D., 2000. Maximum-entropy production principle as a criterion for the morphological-phase selection in the crystallization process. Doklady Physics 45, 129–131. Mauri, R., Shinnar, R., Triantafyllou, G., 1996. Spinodal decomposition in binary mixtures. Physical Review E 53, 2613–2623. Molin, D., Mauri, R., 2007. Enhanced heat transport during phase separation of liquid binary mixtures. Physics of Fluids 19, 074102/1-10. Onsager, L., 1931. Reciprocal relations in irreversible processes. II. Physical Review 38, 2265–2279 (Note that Eq. (5.9) already expresses the principle of maximum entropy production). Sagui, C., Desai, R.C, 1994. Kinetics of phase separation in two-dimensional systems with competing interactions. Physical Review E 49, 2225–2244. Sandler, I.S., 1999. Chemical and Engineering Thermodynamics. third ed. Wiley, New York. (Chapter 7). Siggia, E.D., 1979. Late stages of spinodal decomposition in binary mixtures. Physical Review A 20, 595–605. Van der Waals, J.D., 1979. The thermodynamic theory of capillarity under the hypothesis of a continuous variation of density. Reprinted from: Dutch in Journal of Statistical Physics 20, 167–244. Vladimirova, N., Malagoli, A., Mauri, R., 1998. Diffusion-driven phase separation of deeply quenched mixtures. Physical Review E 58, 7691–7699. Vladimirova, N., Malagoli, A., Mauri, R., 1999. Two-dimensional model of phase segregation in liquid binary mixtures. Physical Review E 60, 6968–6977. Vladimirova, N., Malagoli, A., Mauri, R., 2000. Two-dimensional model of phase segregation in liquid binary mixtures with an initial concentration gradient. Chemical Engineering Science 55, 6109–6118.