Abstract
The development of a methodology for the simulation of structure forming processes is highly desirable. The smoothed particle hydrodynamics (SPH) approach provides a respective framework for modeling the self-structuring of complex geometries. In this paper, we describe a diffusion-controlled phase separation process based on the Cahn-Hilliard approach using the SPH method. As a novelty for SPH method, we derive an approximation for a fourth-order derivative and validate it. Since boundary conditions strongly affect the solution of the phase separation model, we apply boundary conditions at free surfaces and solid walls. The results are in good agreement with the universal power law of coarsening and physical theory. It is possible to combine the presented model with existing SPH models.
1. Introduction
Phase inversion and phase separation play an important role in the preparation of precipitation membranes. Dozens of experiments and simulations have been performed to extend empirical correlations and the understanding of the major aspects of these processes.
Phase separation models based on the Cahn-Hilliard equation [1] describe nucleation and the coarsening dynamics [2] very well. The three consecutive stages of coarsening are diffusion-controlled, viscous-controlled, and inertia-controlled. In the first stage, we use the free energy of a system. Every system tries to minimize its free energy; the strategy is to find a function of the free energy that accurately describes a mixture’s physical behavior. Assuming the local free energy depends on the local concentration and the concentration of the immediate environment, Cahn and Hilliard developed a free energy function based on fundamental thermodynamics. They calculated the chemical potential from the variance of the free energy with respect to density. The chemical potential is the driving force in a linear driving force approach.
After the fundamental work of Cahn and Hilliard, several researchers [2–10] developed phase inversion models, divided into model B and model H in the classification of Hohenberg and Halperin [11]. Weikard [12] studied the dynamics and the stability of the Cahn-Hilliard equation; Saxena and Caneba, Yingrui et al., and others [13, 14] compared these results with experimental data. In the early 1960s, Lifshitz and Slyozov found a universal power law for diffusion-controlled stage [15]. Several numerical methods as, for instance, Lattice Boltzmann method [16], Fast Fourier Transformation [17], or Finite Difference method [18] are used to solve the Cahn-Hilliard equation. The major disadvantage of all these methods is the increasing numerical effort detecting interfaces and geometric changes if diffusion is coupled to hydrodynamic transport. We can reduce the numerical effort of remeshing by using a meshfree simulation method. One promising, meshfree simulation method is the smoothed particle hydrodynamics method (SPH) [19, 20].
For the SPH approach, two models for phase separation are developed. Okuzono [5] first studied phase separation for a binary polymer mixture. His approach is similar to the Cahn-Hilliard approach but differs in the main aspect that diffusion of material is not regarded. Instead of that, he coupled two Navier-Stokes equations. Okuzono distributed different kinds of particles (e.g., particles of types A and B) and described the separation of the particles by viscous and surface tension forces. For the chemical potential, Okuzono used the free energy of the system of Ginzburg-Landau type approximation with an order parameter. Unlike the present approach, Okuzono avoids calculating the second derivative of the density explicitly; instead, he used the first derivative two times in a row. The results of Okuzono are in good agreement with the power law for coarsening dynamics. The model’s major disadvantage is that it only describes phase separation if a phase distribution is known a priori.
Thieulot et al. [21] developed another phase separation model for a binary mixture. It is based on the GENERIC framework. They derived the chemical potential directly from the entropy of the fluid. It takes the same shape originally suggested by Cahn and Hilliard. Their model is thermodynamically consistent, but the effects of diffusive material transport to viscous and inertia forces are neglected.
In this paper, we develop a phase separation model that describes the diffusion-controlled stage using the Cahn-Hilliard equation. It can easily be extended to describe all three stages of coarsening by combination with existing SPH models. As a prerequisite, we introduce an approach to calculate the fourth-order derivative in SPH.
Most models use periodic boundary conditions. But the effects of boundary conditions on the solution of the Cahn-Hilliard equation are not negligible. Therefore, we discuss different kinds of boundary conditions. Since an extension to multicomponent mixtures is straightforward, we consider only a binary fluid mixture in diffusion-controlled stage. Analog to [9], one can extend the model for viscous-controlled and inertia-controlled stages.
In Section 2, we first review the Cahn-Hilliard equation. Subsequently, we introduce the SPH model and an approximation of the fourth-order derivative. After that we discuss different kinds of boundary conditions. In Section 3, we validate the approach considering a diffusion-controlled, binary mixture with periodic boundary conditions. Afterwards, we demonstrate the influence of different kinds of boundary conditions on the solution of the phase separation model and conclude in Section 4.
2. Materials and Methods
The model describes phase separation on a mesoscale. Therefore, we make one fundamental assumption. We assume completed nucleation on microscopic scale for the description with a continuum method. The initial state of the system will contain uniformly distributed nuclei in space.
First, we review a phase separation model based on the Cahn-Hilliard equation. Next, we introduce the Smoothed Particle Hydrodynamics method (SPH) and an approximation for the fourth-order derivative. Then, we apply the SPH method to the phase separation model. For simplicity, we neglect viscous and inertia effects (i.e., the momentum balance). That is valid in case of equal molar volume, density, and negligible excess values. Finally, since boundary conditions have an important influence on the solution, we focus on free surfaces and solid walls.
2.1. Phase Separation Model
Phase decomposition is the result of phase changes as it occurs in, for example, immersion precipitation, condensation, and evaporation. A change of temperature or concentration of a mixture can cause phase decomposition. After exceeding the binodal of the mixture, a fluid mixture switches from stable (miscible) to metastable or unstable (immiscible) state, respectively. In an immiscible system, nucleation occurs spontaneously after exceeding the spinodal. The nuclei grow or disappear depending on their chemical potentials. After infinite time, the system reaches its equilibrium state with homogeneous chemical potential in space.
A phase separation model describes the dynamics of an immiscible fluid mixture based on fundamental thermodynamics. The thermodynamics describes the total free energy of the system. It is the same energy on microscopic, mesoscopic, and macroscopic scale.
In the late 1950s, Cahn and Hilliard [1] developed a model based on statistical thermodynamics and phenomenological observations. The model describes binary phase decomposition very well. They start with a Taylor series expansion of the free energy density around the homogeneous state. Their basic assumption was that the free energy density not only depends on the local concentration but also on the concentration of the immediate environment. That means that they do not neglect higher order gradients of the concentration. As shown later by Fife [8], terms of higher order are necessary for a stable solution in Hilbert space.
The free energy results from Legendre transformation of the internal energy : with and as temperature and entropy. The free energy is also is the volume of the system. We assume an isothermal and isobaric system. Thus, the free energy density only depends on the concentration and its derivatives. As a result of the work of Cahn and Hilliard [1], we use to describe the free energy density. With (3), (2) yields The first term is the free energy density of the homogeneous mixture. It depends only on the local concentration. Cahn and Hilliard supposed using and are Boltzmann constant and temperature. indicates the component.
Two parts represent the energy . One part is that includes the internal energy contribution. The other part is [1] that includes the interface energy contribution. The latter part of (5) is the entropy. One could also use other fundamental equations for based on Flory-Huggins theory [22, 23] or PC-SAFT theory [24] for polymer mixtures.
2.2. Smoothed Particle Hydrodynamics Method
The Smoothed Particle Hydrodynamics method (SPH) is a Lagrangian, particle-based, and meshfree simulation method. Originally developed for astrophysical problems [19, 20], the relevance of SPH increases in engineering science [25].
2.2.1. Kernel Approximation and Second Derivative
In SPH, one calculates a quantity by interpolation of weighted quantities in an appreciable space . The weighting function , called kernel function, approximates a Dirac-Delta function [19]: is the vicinity of . is the smoothing length with, in our case, a constant value of and is the initial particle spacing. The smoothing length represents the width of the Dirac-Delta approximation. Kernel functions are usually spline functions; see [26], for example. In this paper, a quintic spline is used [26].
In discrete formulation of (6), each particle represents one element of fluid. A transition to discrete formulation leads to is the distance between particles and . is the mass of particle . In the following sections, we write instead of .
For example, we could calculate the density with (7) and and obtain But, because of the particle deficiency near boundaries of the computational domain, we introduce a so-called Shepard kernel. The Shepard kernel renormalizes the kernel function that . Therefore, we calculate the density with is the volume of particle .
Brookshaw [27] introduced an approximation of a second derivative: is a property of the fluid, for example, the viscosity. is the distance vector between particles and . Equation (10) is a hybrid formulation for the second derivative using finite differences and the first derivative of the kernel function. In contrast to other formulations of the second derivative, this formulation proved to be less sensitive to particle disorder; we will use this as basis for the fourth-order derivative.
2.2.2. Fourth-Order Derivative
Due to the sensitivity of the second-order and higher order derivatives of the kernel function to particle disorder, it is not possible to do stable simulations using the fourth-order derivative of the kernel function. Therefore, we approximate the fourth-order derivative of a quantity with a biharmonic formulation: This means that the chemical potential as well as the transport equation depends on Laplacian. A similar idea for a second derivative consisting of two first derivatives was published by Watkins et al. [28]. We apply (10) to (11) and get with and are properties of the fluid. Equation (12) is not sensitive to particle disorder and is in good agreement with the power law for the fourth-order derivative. We will discuss this in Section 3.
We determine the error of (12), , analog to Fatehi and Manzari [29]: A comparison with an expanded Taylor series of about leads to and are abbreviated as and . The leading term is proportional to the first derivative. This means that the error of the fourth-order derivative is at least proportional to the error of the first derivative. The leading term in truncation error is at equal particle spacing and at particle disorder is the boundary smoothness, is the particle spacing, and is the displacement of a particle from regular grid.
Even though the approximation, (12), is not exact, it leads to good and reasonable results, as it will be shown in Section 3.
2.2.3. Governing Equations
We consider a binary, isothermal, incompressible, and equimolar fluid mixture. If we neglect viscous and inertia effects, the governing equations are is the density and is the concentration of the mixture. We assume the mobility coefficient to be constant and independent of the concentration. is the chemical potential.
Equation (19) is the general Cahn-Hilliard equation for a multicomponent mixture. The chemical potential of a component is defined by the variance of the free energy to the concentration : If we formulate (19) in SPH using (12), , and , it results in with The model in its presented form is limited to equimolar and incompressible fluid mixtures.
2.3. Boundary Conditions
In a large system, when boundary effects are negligible, it is valid to use periodic boundary conditions to describe phase separation. But there remain systems where boundaries have an effect on the phase separation process.
Free surfaces influence the phase separation of a fluid mixture. There is not a contribution of the boundary to the free energy of the fluid mixture, because of the negligible phase.
A contact angle between two fluids and a solid phase results from solid wall boundary conditions. This was observed in experiments. Now, we consider both kinds of boundaries.
In SPH, we use imaginary particles to model free surfaces by mirroring the particles normal to the boundary. This is important, because of particle deficiency to conserve mass.
Near the solid wall, the influence of the wall free energy on the free energy of the fluid mixture is not negligible. We take this into account with an extension of the free energy: by adding a wall free energy correction term. represents the influence of the free energy density of the wall on the fluids free energy. implies integration over the boundary. In SPH, we reformulate the wall free energy correction term in a volume formulation. We assume that the wall energy is independent of the wall thickness and get In the literature, few approximations of the wall free energy are specified. We use the approximation also used by Yue and Feng [30]; we derive the chemical potential contribution of the wall to the fluid’s chemical potential . It can be written as is the surface tension coefficient and 1, 2, and stand for component 1, component 2, and wall. is a function of the surface tension and the concentration. strongly depends on material properties but only slightly on concentration. is a static contact angle. In the following simulations, we will use independent of concentration, because it is a reasonable choice according to the parameters of the free energy density of the fluid (, , ).
3. Results and Discussion
First, we validate the previously developed model by comparing the time-dependent interfacial area of the system with the universal power law for diffusion-controlled phase separation. Then, we present results of the phase separation with free surface and solid wall boundary conditions and discuss these results.
3.1. Validation
In the early 1960s, Lifshitz and Slyozov found a universal power law for diffusion-controlled coarsening [15]: is the area of the interface and is the time.
We validate our model using a binary mixture (component 1 as and component 2 as ) in two-dimensional space. The computational domain consists of 100 × 100 particles equally spaced with . We use periodic boundary conditions. The model parameters are , , , and . At , the initial concentration is and the uniformly distributed fluctuation is . The time step is .
Figure 1 shows a time series of distribution of the concentration. represents the initial state and represents the equilibrium state of the system. The gray scale represents the composition of the mixture. We represent the concentration of component 1 bright and the concentration of component 2 dark. The time increases from left to right and top down.
(a)
(b)
(c)
(d)
(e)
(f)
The fluid mixture decomposes into two phases. At , the system is in equilibrium state. The equilibrium shape, in two dimensions, is a circle of one phase.
Next, we track the interface size over time and compare it with the power law, (26). We analyzed at different times with Canny’s graphical analysis method [31].
Figure 2 shows the interface size over time . Interface size and time are plotted in double logarithmic reference frame. Crosses in Figure 2(a) represent simulation results for regular particle distribution and the solid line represents the power law.
(a)
(b)
Simulation results and the power law, within resolution error, are in good agreement with regard to the gradient. At times , the interface size is under the resolution of the computational domain. At times , the system is in equilibrium. For quantitative comparison of the surface area calculated with the SPH method, we solved the same equation on a regular grid with grid points discretized with finite differences (FD). The error between SPH and FD is defined as and plotted over time (Figure 2(b)). Within resolution error, the error between both methods is low.
The effect of particle disorder is investigated by repeating the previous simulation with particles randomly displaced in space by of the regular position. In Figure 2(a), points represent the time-dependent surface area versus time. Again, we found good agreement with the power law and quantitative agreement with regular particle distribution. The error to FD (see Figure 2(b)) is within the same range as for regular particle distribution.
3.2. Effects of Boundary Conditions
In the previous simulation, we used periodic boundary conditions. Now, we investigate the effect of free surfaces and solid walls on phase separation.
First, we investigate free surface boundary conditions. The simulation parameters are the previous ones except for the computational domain. In this and the following simulations, we choose the particle number equal to 50 × 50 to reduce the computation time. We only use free surface boundary conditions at all boundaries. As described in Section 2, the gradient of the chemical potential is zero at the boundary: is the normal vector of the boundary. This means that none of the components preferably accumulates at the boundary.
Consider the equilibrium of this system. The system is in equilibrium when the interface size between both phases is minimal. The shape of the equilibrium depends on the initial composition of the fluid mixture. There exist two possible shapes of the equilibrium. One shape is, in two-dimensional case, a circle of one phase. The other shape is a layering of both phases. Now, we analyze the initial composition of the fluid mixture for both shapes.
Consider a circle inside of a square. If layering of both phases is the equilibrium shape, the circumference will be larger than the length of the edge of a square : is the diameter of the circle. The limiting composition from the volume of the circle . If we know the initial composition, the area of the circle at equilibrium is known. We get the perimeter of the circle with simple conversions. If we assume the limiting case , we get two solutions Consequently, the range of the initial composition for layering is .
Figures 3(a)–3(f) show simulation results with free boundary conditions and . From (29), we expect layering of both phases as equilibrium shape. At time , the equilibrium shape actually is a layering (see Figure 3).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
For quantitative comparison of these results, we also calculated the same system using finite differences on a regular grid with 2500 grid points. Figures 3(g)–3(i) show the results at three different times. It is not exactly the same because the random initial conditions are different. Nevertheless, the evolution of phase decomposition is similar and the equilibrium states are identical.
Next, we qualitatively investigate solid wall boundary conditions. As described in Section 2, we use a wall free energy correction term, (25), to describe the effect of a wall on a fluid mixture. For validation of the wall free energy correction term, we consider a droplet of one fluid phase, surrounded by a second fluid phase, on a solid wall. Keep in mind that a contact angle is defined by the intersection of a fluid-fluid interface and a fluid-solid interface. The magnitude of the contact angle depends on the wall free energy correction term. We specify an equilibrium contact angle of the droplet and compare it with the equilibrium contact angle of the droplet in the simulation.
The computational domain and the simulation parameters, except for and , are the previous ones. Figure 4, top left, shows the initial state of the composition. The composition of the smaller square is and of the larger one . is zero. We apply solid wall boundary conditions to all boundaries.
(a)
(b)
(c)
(d)
(e)
(f)
We vary the specified contact angle between and . Figure 4 shows the results. All specified static contact angles are obtained. The interface between the phases is smooth. This is obvious in cases and . In the first case, it seems that the contact angle is not still larger than zero; in the second case, it seems that the droplet moves away from the boundary. Regarding the smoothness of the interface, the droplet is in infinitesimal small contact with the boundary. If we increase the volume of component 2, in case of , the complete wall is wetted.
As a last point, we qualitatively consider phase separation of a binary fluid mixture in a box. The computational domain and parameters are the previous ones. We apply solid wall boundary conditions to all boundaries. The contact angle is . This means that all walls are fully wetted with component 2.
Figure 5 shows the results. At equilibrium, , the solid wall is fully wetted with component 2. Thus, if the circle is centered in the box, the minimum of the free energy of the mixture is obtained. For contact angle (not shown), we obtain the same equilibrium state with changed composition of mixture.
(a)
(b)
(c)
(d)
(e)
(f)
4. Conclusions
We developed a SPH approach of a phase separation model. The model describes diffusion-controlled, equimolar coarsening dynamics based on the Cahn-Hilliard equation. The model is extensible to multicomponent systems and nonequimolar systems, as shown in the literature. Since a diffusion equation can be combined with existing SPH models, see, for example, [32], combination of the present model with existing SPH models is possible.
We calculated the chemical potential of a fluid mixture from free energy and its second derivative. The gradient of the chemical potential is proportional to the diffusion flux; see transport equation (19). In this way, the free energy of immiscible components intrinsically takes care of surface tension effects on coarsening dynamics, because of the second derivative of concentration. Since phases with smaller volume have larger chemical potentials, because of larger curvature, the system is driven from many small volumes to one large volume. Therefore, larger volumes will grow while smaller ones disappear.
The novelty of the present SPH model is a direct Cahn-Hilliard approach using a fourth-order derivative of a quantity . Because of its biharmonic property, we split the fourth-order derivative into two derivatives of second-order.
We compared the time-dependent interface size of the computational domain with the universal power law for diffusion-controlled phase separation, and we found good agreement with the power law.
In the last part, we investigated free surface and solid wall boundary conditions. We validated the boundary condition with the equilibrium state of the system. For solid wall boundary conditions, we introduced a phenomenological approach. This approach is based on a wall free energy correction term. We applied the solid wall boundary conditions to model static contact angles. Finally, we applied solid wall boundary conditions to present model.
The extension to three or more component systems is straightforward, as seen in the literature [9]. Coupling of the transport equation to the momentum balance is possible in straightforward way. For nonequimolar fluid mixtures, variable mass or volume of the SPH particles should be considered.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
The work is supported by the German Science Foundation (DFG) SFB 716 project and within the funding program Open Access Publishing.