Impact of complex topology of porous media on phase separation of binary mixtures

See allHide authors and affiliations

Science Advances  22 Dec 2017:
Vol. 3, no. 12, eaap9570
DOI: 10.1126/sciadv.aap9570


Porous materials, which are characterized by the large surface area and percolated nature crucial for transport, play an important role in many technological applications including battery, ion exchange, catalysis, microelectronics, medical diagnosis, and oil recovery. Phase separation of a mixture in such a porous structure should be strongly influenced by both surface wetting and strong geometrical confinement effects. Despite its fundamental and technological importance, however, this problem has remained elusive for a long time because of the difficulty associated with the complex geometry of pore structures. We overcome this by developing a novel phase-field model of two coupled order parameters, the composition field of a binary mixture and the density field of a porous structure. We find that demixing behavior in complex pore structures is severely affected by the topological characteristics of porous materials, contrary to the conventional belief that it can be inferred from the behavior in a simple cylindrical pore. Our finding not only reveals the physical mechanism of demixing in random porous structures but also has an impact on technological applications.


Wetting effects on phase separation in random porous materials of mesoscopic scale (118) are of great technological interest because they are relevant to many important industrial applications, for example, phase separation in porous materials can be used for the coating of the internal surface of a porous material (17), and the stability of phase-separated domain structures is a key to oil recovery processes (16). Thus, there is a high demand for the basic understanding of phase-separation behavior in a complex pore geometry. However, there are fundamental questions that remain to be answered.

The difficulty originates from the complex dynamical interplay between phase separation, surface wetting, and geometrical confinement. For example, phase separation in porous materials cannot proceed indefinitely and is trapped in a metastable state, unlike phase separation in bulk. The key question is, “Which pore characteristics affect the kinetic pathway of phase separation and its arrest to a final metastable state?” The two important characteristics of pores are surface wettability and geometrical structure. The former, which characterizes the friendliness of the pore surface to the components of a mixture, can be grouped into three cases: neutral, perfect, and partial wetting (see Fig. 1). On the other hand, the latter can be characterized by its size and topology. The size effect is rather well understood. It can be characterized by the relationship between two key length scales: the characteristic pore size λ (here, we define it as the pore diameter) and the domain interface thickness, which is comparable to the correlation length of composition fluctuations ξ (strictly, Embedded Image). The correlation length ξ increases when approaching a critical point Tc as ξ = ξ0(|TTc|/Tc)ν, where ξ0 is the molecular size and ν is the critical exponent [≈0.63 for a three-dimensional (3D) system]. Thus, ξ is large near a critical point or for polymers and is usually in the range of 1 to 100 nm. Then, the situation can be classified into two cases in terms of the ratio of λ/ξ: (i) For λ ≤ ξ, the randomness in the pore geometry and/or surface chemistry applies a quenched random field (pinning) to the fluid, resulting in separation into domains larger than the pore size. Thus, the randomness is expected to result in random field Ising model behavior (19, 20), that is, the collective behavior of magnetic spins under random magnetic fields. (ii) On the other hand, for λ ≫ ξ, a simple pore model is proposed to work: The metastability and the slow kinetics of the coarsening process can be explained by the geometric confinement of phase-separated domains in a cylindrical pore (2124). Now we have a consensus that case (ii) is relevant for most practical applications. However, this is not the end of the story, because there is a possibility that not only the size but also the complex topology of a porous material matters. However, the impact of complex topology of a porous material on phase-separation behavior has not been addressed to date because of the following reasons.

Fig. 1 Wetting of a droplet on a flat wall.

(A) The wetting behavior of a droplet on a flat wall for different wetting parameter h1. h1 = 0, −2, and −10 from left to right. (B) The cross section of a droplet and a wall along x = Lx/2 for the case of neutral wetting (h1 = 0).

One reason is the experimental difficulty in studying phase-separation behaviors in a mesoscopic pore, particularly in random porous materials, because not only is the pore radius of Vycor glasses too small for the direct real-space observation but also extracting detailed real-space information from scattering experiments is extremely difficult. There have been detailed experimental studies on phase separation in rather macroscopic pores, which allows real-space observation (24). Only recently, phase separation in more complex geometries has been studied experimentally: Kanamori et al. (17) observed coarsening dynamics of phase-separating mixtures using the scanning electron microscopy technique. Synchrotron x-ray tomography measurements were also applied to phase separation of oil/water mixtures in porous materials (25, 26). However, there have been few experimental studies that can address the role of pore topology in the demixing behavior. Another reason is the difficulty in theoretical and numerical studies because of the complex nature of pore structures. Some numerical works revealed phase separation under surface fields for simple geometries such as a thin film or a simple cylindrical pore by means of the phase-field model (6, 12, 27) and molecular dynamics (MD) simulations (12, 28). However, actual porous materials have much more complicated labyrinth-like pore structures than a simple cylindrical pore, which makes it difficult to reveal how phase-separated domains wet and grow in such complex confining geometries. There were some numerical studies on phase demixing in complicated geometries by the phase-field model (29, 30), the lattice Boltzmann method (31), and MD simulations (32). However, these studies are limited to phase demixing in 2D porous structures and to the rather early stage of demixing. Because there are many metastable states in porous materials under an influence of wetting, it is crucial to understand a nonequilibrium kinetic pathway in which phase-separated domains wet walls and coarsen in random porous materials. It is expected that the spatial dimensionality and the topology of porous media should have significant impacts on the kinetic pathway of phase separation and the resulting structural formation, but their roles have not been understood yet even at a primitive level.

Here, we address this problem by numerically studying phase separation in random porous materials in both 2D and 3D. To this end, we develop a novel phase-field method, which allows us to deal with a variety of porous structures. We reveal an intrinsic difference in pattern evolution between 2D and 3D porous structures and the crucial role of topological characteristics of 3D porous materials in the selection of the final phase-separated structure.


Simulation method

We consider a binary mixture immersed in a porous material, whose surface interacts with the mixture via surface wettability. To study the phase-separation kinetics in porous materials, we use a phase-field model described by two order parameters, which represent the composition field of the binary mixture, ψ(r), and the density field representing the porous material, ρ(r), with their coupling strength, that is, wettability, h1. Here, r is the position vector. Throughout this article, we consider a symmetric (50:50) binary mixture of the average composition Embedded Image. The porous structure is represented by the density field of ρ(r) and its surface expressed by a smooth but sharp interface. The details of this method are explained in the Materials and Methods. This method has several merits: We can produce various porous structures by using either a physical process of demixing or a mathematical function that represents various regular pore structures (33). Furthermore, we can introduce hydrodynamic degrees of freedom (27) straightforwardly in the same manner as our fluid particle dynamics method to simulate colloidal suspensions (34), although here, we consider only diffusion as the transport mechanism. Furthermore, we note that our phase-field model can incorporate charges, additional component, and liquid-crystalline order in a straightforward manner.

We check the validity of our simulation method by studying a well-known problem of wetting of a liquid droplet on a flat solid substrate, as shown below. We prepare complex porous structures by freezing spinodal decomposition concerning the ρ(r) field at a certain timing, which is similar to the preparation method of the so-called Vycor glass, which is formed by spinodal decomposition of a borosilicate glass. The porous structure is characterized by the volume (area) fraction of the solid porous matrix Φ (that is, the mean surface curvature H of the porous material) and the average pore size λ.

Results for a simple wetting geometry

First, we check the validity of our method by studying wetting behaviors of a droplet on a flat wall. As the solid wall, we set the wall position to be y = L1 and y = L2, and the interfacial profile of the wall is given by Embedded Image. As an initial condition, we put a droplet on the bottom wall with the contact angle Embedded Image. Thus, the initial profile isEmbedded Image(1)

Here, the center of droplet rc = (L1, Ly/2), and Rd is the radius of the droplet. In this simulation, we set Lx = Ly = 128 and L1 = 20, L2 = 108, and Rd = 30. Figure 1A shows the stationary wetting behavior of a droplet for the wetting parameters h1 = 0, −2, and −10, respectively. We note that h1 = 0 or −10 should correspond to the neutral wetting, that is, the equilibrium contact angle of 90°, and the complete wetting, that is the equilibrium contact angle of 0°, respectively. This is confirmed in Fig. 1A, indicating the validity of our simulation method. We also show the equilibrium spatial profile of the two fields, ψ and ρ, in the direction perpendicular to the wall in Fig. 1B.

Neutral wetting case

Demixing in 3D porous materials. First, we consider the case of neutral wetting, that is, h1 = 0 (θ = 90°), in a 3D symmetric bicontinuous porous structure (H = 0, that is, Φ = 1/2). Figure 2 shows the time evolution of the concentration field ψ in the porous structure. Here, the black interface shows the surface of the porous structure by the iso-interface of ρ = 0.8. The blue and green interfaces represent the iso-surfaces of the A and B phases, ψ = 0.6 and ψ = −0.6, respectively. Initially, wetting layers are not formed on the wall. Instead, small tortuous domains contact the wall with the contact angle of ≈90° (see the image at t ≈ 40 of Fig. 2A). The temporal change of the structure factor S(k) (where k is the wave number) is shown in Fig. 3A. Because our pore size λ is much larger than the correlation length ξ of the ψ field (λ/ξ ≈ 20), the peak wave number k1 of S(k) is about the same as that in bulk in the initial stage of phase demixing, and the initial growth of composition fluctuations is not affected so seriously by the spatial confinement. Accordingly, the growth exponent seen in Fig. 3B is not so different from the well-known growth exponent in bulk, 1/3 (35). We note that Bailey et al. (36) performed pressure-quench experiments of binary mixtures confined in dilute silica gels and reported the similar domain growth exponent of 1/3.

Fig. 2 Time evolution of the concentration field for a 3D random symmetric porous structure in a neutral wetting case (h1 = 0).

(A) Two demixed phases [the A (green) and B (blue) phases] together with the porous structure surface (black). See also movie S1. (B) Only the A phase but without the wetting layer formed on the surface of the porous structure by drawing the iso-interface of ψ = 0.6 in the pore space where ρ < 0.173. (C) Only the B phase by drawing the iso-interface of ψ = −0.6. Here, the green and blue interfaces represent the iso-concentration surface of ψ = 0.6 for the A phase and that of ψ = −0.6 for the B phase, respectively.

Fig. 3 Structural evolution in k-space for a binary mixture in porous structures.

Time evolution of S(k) in a neutral wetting (h1 = 0) case (A) and the first moment of wave number k1 for various wettability (B) for a 3D symmetric porous structure. Time evolution of S(k) in a neutral wetting (h1 = 0) case (C) and the first moment of wave number k1 for various wettability (D) for a 2D symmetric porous structure. For both 2D and 3D, the characteristic wave number k1 decreases as t−1/3 in the intermediate stage for h1 = 0 but more slowly for h1 ≠ 0.

As domains grow, they start to exceed the pore size and, thus, bridge the pore structure around t ≈ 500. However, even after the formation of these plug structures, domains keep growing, as can be seen in Fig. 2. The motion of each interface of disconnected domains is determined by the difference of the curvature of neighboring domain interfaces. The curvature difference between surfaces of neighboring plugs leads to the difference in the chemical potential between them, which induces the diffusion flux from a domain of high interface curvature to that of a low one. This is basically the same as the evaporation-condensation (or Lifshitz-Slyozov-Wagner) mechanism in bulk phase separation (35). In the bulk case, the domain interface curvature of disconnected droplets is inversely proportional to the droplet radius, and the diffusion continuously takes place from smaller to larger droplets. Unlike this case, in a porous medium, the curvature of the domain interface is determined by the confining pore-wall structure and independent of the domain size (or the plug length). Thus, the coarsening rate, which is significantly slower than that in bulk due to the confinement, further slows down with time, and, eventually, the coarsening almost stops around t ≈ 100,000 when the interface curvatures of all the plugs become almost identical (see Fig. 3B). We may say that the domain structure in a porous media self-organizes to attain a configuration of homogeneous domain interface curvature. This feature can be more clearly seen in the 2D case (see Fig. 4A).

Fig. 4 Time evolution of the concentration field in 2D random symmetric porous material.

(A) A neutral wetting case (h1 = 0). The black region represents the symmetric porous structure, and the red and white domains show the A and B phase, respectively. See movie S2. (B) A complete wetting case (h1 = –10). See movie S4. (C) A partial wetting case (h1 = –2).

Next, we discuss a very interesting phase-separation structure formed under a neutral wetting condition. To see the structure and connectivity of each of the phase-separated two phases separately, we show in Fig. 2 (B and C) only the A (green) and B (blue) phases by the iso-surface of ψ = 0.6 and ψ = −0.6, respectively, for the inside of the pore region where ρ < 0.173. From these images, we can clearly see that a 3D percolated network structure of each phase is formed in the bicontinuous pore structure while the plug state blocks each pore. That is, two percolated network structures are simultaneously formed, thanks to both the symmetry provided by neutral wetting and the symmetry of the pore structure itself (Φ = 0.5, or H = 0). Because the curvature of the domain interface is minimized when it is located in a junction part of the porous network structure, the domain structure is self-organized into this interesting structure by the evaporation-condensation mechanism. We note that this bicontinuous double-network structure in a porous material is formed to minimize the domain interface energy. Reflecting that the characteristic length of this network pattern is about twice of that of the porous structure, we can see a distinct peak in the structure factor at the half of the peak wave number of the porous structure (see Fig. 3A).

This unique double-network phase-separation structure may be useful for many applications, for example, for coating its inner pore surface such that the two percolated surfaces have distinct wetting properties. Furthermore, the percolated nature of each phase should provide unique transport properties.

Demixing in 2D porous materials. To study the role of the dimensionality in pattern formation, we also performed numerical simulations under the same condition as the above 3D case for a 2D symmetric porous structure (Φ = 1/2). We show the time evolution of phase separation in 2D in Fig. 4A, where the red and white domains show phase-separated domains and the black regions represent the porous matrix. The temporal changes in S(k) and k1 are also shown in Fig. 3 (C and D, respectively). The basic behaviors in 2D are similar to that in 3D, but there is a crucial difference in the pattern formed in the late stage between 2D and 3D. In 2D, a network structure is immediately destroyed, a plug state is formed, and then the lengths of plugs increase with time by the evaporation-condensation mechanism (see Fig. 4A). In 3D, on the other hand, bicontinuous double networks of the two phases can be formed, as shown in Fig. 2. This shows that there is an essential effect of the dimensionality of a pore structure on the connectivity of phase-separated domains. The difference simply comes from the fact that in 2D, a bicontinuous porous structure can never be formed even for the symmetric situation unlike in 3D.

Perfect wetting case

Demixing in 3D porous materials. Next, we consider the case of perfect wetting, that is, h1 = −10 or θ = 0°. Figure 5 shows the time evolution of the concentration field in the 3D symmetric bicontinuous pore (H = 0, or Φ = 1/2), which is similar to that in Fig. 2A. The black regions show the internal surface of the solid porous structure. The green and blue surfaces show the more wettable A phase and the less wettable B phase, respectively, by the iso-surfaces of ψ = 0.6 and ψ = −0.6. We also show in Fig. 5C only the blue surface of the less wettable B phase by the iso-surface of ψ = −0.6, whereas in Fig. 5B, we show the more wettable A phase by the iso-interface of ψ = 0.6 for the pore region where ρ < 0.173, as above. Note that the choice of ρ < 0.173 is not to show the surface wetting layer of the A phase to make the A phase in the middle part of the pore visible. In the initial stage just after the quench, the more wettable A phase (ψ = 1) wets the surface of the porous material: In Fig. 5A, we can see that the green surface covers the black wall around t = 3. Then, the layered tube structure composed of two parts (the surface wetting layer of the A phase and the core tube of the B phase) is formed around t = 10 (see Fig. 5A). This process of the layered tube structure formation is a consequence of the surface-directed spinodal decomposition (4, 12). The emergence of a domain structure is much faster than for the neutral wetting case (see Fig. 3B). This is because strong adsorption of the more wettable phase to the wall accelerates domain growth.

Fig. 5 Time evolution of the concentration field for a 3D random symmetric porous structure in a complete wetting case (h1 = –10).

(A) Two demixed phases [the more (green) and less wettable (blue) phases] together with the porous structure surface (black). See also movie S3. (B) Only the more wettable phase but without the wetting layer formed on the surface of the porous structure by drawing the iso-interface of ψ = 0.6 in the pore space where ρ < 0.173. (C) Only the less wettable phase by drawing the iso-interface of ψ = −0.6. Here, the green and blue interfaces represent the iso-concentration surface of ψ = 0.6 for the more wettable phase and that of ψ = −0.6 for the less wettable phase, respectively.

After t ≈ 10, the internal phase-separated structure is formed and then domains grow. During this time regime, the temporal change of k1(t) can be described by the relation t−1/3, as seen in Fig. 3B. This growth exponent is the same as that observed in the growth process of a wetting film on a plain wall. After t ≈ 100, the structure factor at a higher k region does not change. In the late stage, narrow pores are filled by the more wettable A phase. This process can also be seen as a breakup of the tube configuration of the less wettable B phase, but the percolated nature of the B phase is still preserved. In this regime, the temporal change of k1(t) becomes slower and weaker compared to that in the initial stage, and after t = 10,000, the domain coarsening stops (see Fig. 3B). Although the phase-separated domain structure is very different between the neutral and complete wetting cases, the two phases are percolated for both cases: the double-network structure for neutral wetting and the double-tube structure for complete wetting.

Demixing in 2D porous materials. Now, we turn to the corresponding 2D case under the same complete wetting condition (see Fig. 4B). Here, the red and white domains correspond to the more and less wettable phases, respectively. The initial demixing process is similar to that in 3D, as seen in Fig. 4B. After the initial surface wetting stage, in which the three-layer configuration is formed, the inner tube part becomes unstable and capillary bridges are formed. In the entire process, the pore surface is always covered by the wetting layer. Then, the thinning of surface wetting layers and the growth of capillary bridges take place. Then, a crucial difference from the 3D case starts to emerge. Unlike the 3D porous structure, there are appendix parts (terminal parts of tubes) in the 2D porous structure because there is no bicontinuous porous structure for 2D. The appendix parts of the porous structure gradually get filled with the more wettable red (A) domains in the late stage. This is because the total interfacial area between the two phases can be reduced by filling appendix regions, which reduces the total interfacial area of the tube parts. Locally, this is driven by the chemical potential gradient between the interfaces with a more negative curvature of the more wettable A phase filling an appendix and the interface of tube parts with a less negative curvature. This diffusive transport process leads to slow but continuous growth of the domain size, as shown in Fig. 3D. We perform independent simulation runs for h1 = 0 and h1 = −10 to see how reproducible the structure formation is. In the case of partial wetting, the final domain configurations are different from run to run, and thus not reproducible. In the case of complete wetting, on the other hand, the final domain configurations are reproducible. This indicates that the nearly equal domain interface curvatures of plug ends can be realized in many ways for neutral wetting, but appendix parts have deterministic roles in the selection of the domain configuration for complete wetting.

Partial wetting case

Demixing in 3D porous materials. Next, we consider a partial wetting case intermediate between neutral and complete wetting, that is, a case of 0 > h1 > −10. Here, we set h1 = −2 (see Figs. 4C and 6, respectively, for 2D and 3D). As in the case of complete wetting, a thin wetting layer and phase-separated stripe domains are formed in the initial stage, but the thickness and composition difference is smaller compared to the complete wetting case. Then, the less wettable B phase domains are deformed and get contact with the surrounding pore walls. After this process, there is a crucial difference between 3D and 2D. For 3D, a double-network structure can be formed, although it is not so symmetric compared to the neutral wetting case (see Fig. 6). The temporal change of k1(t) is summarized in Fig. 3B for 3D.

Fig. 6 Time evolution of the concentration field for a 3D random symmetric porous structure in a partial wetting case (h1 = –2).

(A) Two demixed phases [the more (green) and less wettable (blue) phases] together with the porous structure surface (black). (B) Only the more wettable phase but without the wetting layer formed on the surface of the porous structure by drawing the iso-interface of ψ = 0.6 in the pore space where ρ < 0.173. (C) Only the less wettable phase by drawing the iso-interface of ψ = −0.6. Here, the green and blue interfaces represent the iso-concentration surface of ψ = 0.6 for the more wettable phase and that of ψ = −0.6 for the less wettable phase, respectively.

As we can see above, the behavior for 0 > h1 > −10 is a complicated combination of the behaviors of complete and neutral wetting. The behavior changes gradually with an increase in the wetting parameter h1 from the behavior of complete wetting to that of neutral wetting.

Demixing in 2D porous materials. For 2D, on the other hand, the elongated domains are disconnected and, thus, the pore structure is divided by disconnected plugs (see Fig. 4C). After this happens, the subsequent growth is similar to the case of neutral wetting. The temporal change of k1(t) is summarized in Fig. 3D for 2D. Similar to the above 3D case, the behavior changes gradually with an increase in the wetting parameter h1 from the behavior of complete wetting to that of neutral wetting.

Demixing in off-symmetric porous structures. So far, we focus on the effects of compositional symmetry on phase demixing in a symmetric bicontinuous porous structure (Φ = 0.5). Next, we turn our attention to the effects of the symmetry of the porous structure. As an example, we study phase separation in an asymmetric porous structure (Φ = 0.65). The results are shown in Fig. 7. We can see that even for neutral wetting, a double-network structure cannot be formed and the B phase is disconnected. For complete wetting, a double-tube structure is formed even in this case. For partial wetting (h1 = −2), the situation is marginal. This indicates that a double-network structure is favored for higher wetting symmetry (that is, a more neutral wetting condition) and higher surface curvature symmetry of the solid matrix.

Fig. 7 Time evolution of the concentration field for an asymmetric bicontinuous porous structure (Φ = 0.65).

(A) Neutral wetting (h1 = 0). (B) Partial wetting (h1 = –2). (C) Complete wetting (h1 = –10).

Comparison with scattering experiments

Now, we compare the temporal change of S(k) between our results and the experimental ones. Lin et al. (37) experimentally studied a water/lutidine mixture of a critical composition confined in a Vycor glass. For small-angle neutron scattering experiments, they made contrast-matching between the mixture and the Vycor glass by appropriately mixing protonated and deuterated water. They interpreted the results as a combination of a lutidine-rich wetting layer coating the internal surfaces of the Vycor glass and random single-phase domains occupying small pores. Our structure factors show the behavior similar to their results. However, it is worth mentioning the following possibility: Contrary to the assumption made in the previous works that the domain size does not exceed the pore size, we reveal that 3D double-percolated network structures can be formed for a neutral wetting condition. This means that the characteristic size of domains along pores can be macroscopic. Here, we note that such a situation was experimentally observed by Iglauer et al. (25) for sandstones containing a residual nonwetting phase. Their synchrotron x-ray tomography showed that the residual phase forms large clusters, spanning over many pores. These facts indicate the limitation of the k-space approach and the importance of real-space observation in revealing the domain structure in a complete geometry.


In summary, we have studied phase-separation behavior in mesoscopic random porous structures by numerical simulations. This allows us to compare the time evolution of the structure factor of binary mixtures confined in a random porous structure between numerical simulations and experiments for the first time. In the initial process of phase separation, the domain size grows as ~t1/3 for both complete and partial wetting conditions. This apparently looks similar to bulk, but the real-space analysis reveals that the domain orientation is strongly affected by the difference in the wettability between the two phases, which cannot be easily inferred from scattering experiments.

We find a significant impact of spatial dimensionality by comparing results of 3D and 2D. The breakdown of wetting symmetry between the two components, that is, a small degree of preferential wetting, leads to the formation of a wetting layer on the surface of pores. This seriously affects the final domain configuration. Only for a neutral wetting condition, the two phases equally wet the pore surface. For 2D, this leads to the formation of a plug configuration. This is the only way to keep the wetting symmetry in 2D. For 3D, on the other hand, there is a novel configuration that allows a system to keep the wetting symmetry, that is, the formation of double-percolated bicontinuous structures. This may provide us with interesting applications. For example, we can put the two distinct surface wetting properties to two channels of percolated bicontinuous pores by using a coating technique (17). We can also fill two types of liquids with very different properties so that each of them is percolated in a bicontinuous pore structure.

We also reveal a strong impact of complex topology of porous media on demixing, which makes phase-separation behavior fundamentally different from that in a simple cylindrical pore. One crucial difference comes from the variation of the pore radius in random porous materials. What is important here is that the assumption that the total composition is fixed, which is valid for a simple pore geometry, can be violated locally for random porous media. In our system, the overall composition averaged over a system is of course fixed during phase separation. However, because the pore radius varies from place to place, the composition is not conserved locally and varies spatially even though it is conserved globally. This feature leads to the formation of a much more heterogeneous complex phase-separation structure in random porous media, unlike a rather periodic structure formed in a cylindrical pore (9, 24). Another important difference is the tortuous structure of a bicontinuous network pore. In a simple cylindrical pore, domain coarsening stops when the plug size becomes comparable to the pore radius. This is because the domain interface curvatures are simply determined by the unique characteristic length scale of the tube, that is, the pore diameter, and thus become homogeneous. On the other hand, the random tortuous porous structure allows domains to further coarsen even after exceeding the local pore radius, because inhomogeneous domain interface curvatures due to variable local pore radius and the presence of junctions keep the evaporation-condensation mechanism active. Our study shows that these fundamental differences lead to drastic differences in the physical mechanism of structural formation between simple and complex pore geometry. Thus, the randomness, topology, and bicontinuity of pore structures are quite important in structural formation and its kinetics, which should be relevant to many industrial applications.

Finally, we also mention a potential application of our method to study fluid dynamics in a complex porous media. For example, the dynamics of imbibition, by which a wetting liquid is drawn into a porous medium, is a very important problem in both technological and geophysical viewpoints (38). Because we can easily incorporate the hydrodynamic degrees of freedom into our method [see the study by Tanaka and Araki (34)], it will be useful to study fluid transport through a porous media.


A new phase-field simulation method using two order parameters

We developed a novel two-order parameter model, which can describe the phase-separation process of binary mixtures in arbitrary pore structures. In this model, we adopted the concentration difference of binary mixtures ψ and the solid-phase composition ρ as the two relevant order parameters. The solid wall was expressed by a smoothly varying field ρ, and thus, the interface between the solid phase and the empty phase was described by the surface with a finite thickness. The usage of a smooth interface to represent a solid object was proposed to describe a solid colloidal particle (34). This is a natural extension of the phase-field model, which is often used for studying crystal growth and melting. However, there was no time evolution of the ρ field, which just acts as a solid wall for the ψ field. In this model, the total free energy of a system is expressed asEmbedded Image(2)whereEmbedded Image(3)

Here, f(ψ, ρ) is the free-energy density of the system, which has three minima: The two of them, (ψ, ρ) = (±1, 0), correspond to the two phase-separated phases, whereas Embedded Image correspond to the solid matrix. The second term in the right-hand side of Eq. 2 describes the surface wetting energy at porous walls. The parameters h1 and g represent the strength of the short-range surface field and the surface enhancement field, respectively. The third term in the right-hand side of Eq. 2 is the gradient term of the concentration field. Here, we consider only a case of short-ranged interactions with the solid wall. In the thin interface limit (ξρ → 0), the total free energy is expressed asEmbedded Image(4)where ψs is the concentration difference at the wall ψs = ψ(0). Then, we described the dynamic evolution of ψ by the following generalized diffusion equationEmbedded Image(5)

This is an extension of the so-called model B equation in the Hohenberg-Halperin notation (39). We set the transport coefficient of the concentration difference as Embedded Image so that it is zero inside the solid matrix.

Using the method described above, we can deal with any type of wall structures. Here, we use a random porous wall structure. We prepared it by numerically solving model B for the bulk spinodal decomposition of a 50:50 binary mixture and stopping the simulation in the late stage at a time when the structure becomes a desired size (29, 33). This model porous structure mimics a Vycor glass, which is made by acid leaching of the boron-rich phase of a spinodally decomposed borosilicate glass. The porosity of this system is 50%. The peak wave number of the structure factor of the porous material is kpk1 ≈ 0.2; thus, the average pore size is about 2π/kp = 31. This is about 30 times of the interfacial thickness of domains in bulk.

For simulations, we used the above equations after scaling the time and space by the interface length, ξ, and the characteristic diffusion time of concentration field, τ = ξ2/D. We solved the scaled kinetic equations by the Euler method with a periodic boundary condition. In experiments for a solid pore immersed in a reservoir liquid mixture, exchange with the reservoir can lead to a very slow change in ψ when we vary the temperature, and capillary phase separation can also happen when the pore size is small. In our simulations, on the other hand, the concentration is strictly fixed due to the periodic boundary condition, and the capillary phase separation does not happen. We set the system size to be 256 and the grid size to be Δx = Δy = Δz = 1. For simplicity, we do not introduce temporal thermal fluctuations. The initial profile of ψ is given by Embedded Image, where ψ0 is the average concentration difference and Δψ is the variance of initial random composition fluctuations. We set ψ0 to be 0 and Δψ = 0.01. We define the interface between two immiscible fluids as a surface where the value of the concentration difference ψ is equal to 0 and at the same time the value of the wall profile parameter ρ is smaller than Embedded Image.

Calculation of the structure factor S(k) and the peak wave number k1

To extract the characteristic lengthscale of phase-separated structures, we calculate the power spectrum of ψ. The structure factor S(k) is obtained by spherically averaging the power spectrum of ψ. The characteristic wave number is calculated asEmbedded Image(6)

Because we set the average concentration difference ψ0 to be 0, the structure factor S(k) obtained in this way should correspond to that experimentally observed by light scattering or neutron scattering under an index-matched condition.


Supplementary material for this article is available at

movie S1. Movie corresponding to Fig. 2A.

movie S2. Movie corresponding to Fig. 4A.

movie S3. Movie corresponding to Fig. 5.

movie S4. Movie corresponding to Fig. 4B.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: Funding: This work was partially supported by Grants-in-Aid for Specially Promoted Research (25000002) from the Japan Society of the Promotion of Science. Author contributions: H.T. conceived and supervised the project, R.S. performed simulations, and R.S. and H.T. discussed and wrote the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article