Phase-field simulation for the formation of porous microstructures due to phase separation in polymer solutions on substrates with different wettabilities

The porous microstructure has been widely observed in a variety of polymer solutions that have been broadly applied in many industry fields. Phase separation is one of the common mechanisms for the formation of the porous microstructure in binary polymeric mixtures. Previous studies for the formation of porous microstructures mostly focus on the separation of the bulk phase. However, there is a paucity of investigation for the phase separation of polymer mixtures contacting the solid substrate. When the polymeric liquid mixtures interact with the solid substrate, the wetting boundary condition has to be taken into account. In this work, we present a phase-field model which is coupled with the wetting boundary condition to study the phase separation in binary polymer solutions. Our consideration is based on the polymerization-induced phase separation, and thermally induced phase separation by using the Flory–Huggins model. By taking the wetting effect into account, we find that polymer droplets spontaneously occur in the microstructure, even though the bulk composition is outside the spinodal region. This phenomenon is caused by the surface composition resulting from the wetting effect that was often overlooked in literature. For the phase separation in the binary polymer mixture, we also study the impact of the temperature gradient on the microstructural evolution. The porosity, the number of droplets, and the mean radius of the droplets are rationalized with the temperature gradient.


Introduction
Phase separation is known as a common mechanism for the formation of the porous microstructure of polymeric membranes [1]. The porous microstructures resulting from phase separation have found wide applications in various fields [2,3]. In the formation process of porous polymeric microstructures, a phase separation process takes place in which the polymer solution separates into two immiscible phases, known as spinodal decomposition (SD) [2,4]. The SD mechanism is based on a double well free energy density as a function of the composition, where there are two inflection points. Subjected to thermal fluctuations, the polymer solution is energetically unstable. The polymer solution separates into a polymer-rich phase and a polymer-lean phase when the initial composition is between the inflection points, as demonstrated by the linear stability analysis [5]. From the aspect of thermodynamics, the double well free energy density and the resulting SD appear when the enthalpy contribution to the system free energy density is comparable to or greater than that of the entropy. This is commonly observed in the polymer solution since the increase in the degree of polymerization reduces the entropy according to the lattice mean-field theory, such as the Flory-Huggins model. Specifically speaking, the enthalpy is related to the difference between the solventpolymer binding energy and the summation of solvent-solvent and polymer-polymer bonding energies. When the effect of the entropy becomes less pronounced due to polymerization, the enthalpy dominates the free energy minimization, and the polymer-polymer and solvent-solvent bonds are energetically favorable. The evolution of SD involves the minimization of both bulk and interfacial energies, and it is typically governed by the so-called Cahn-Hilliard equation, which is a kind of Euler-Lagrange equation within the scope of the minimization problem [6].
Previous studies on the SD of polymer solutions have primarily focused on the simulation of the bulk phase separation, for example, using the phase-field method [7] and the lattice Boltzmann method [8]. Although the previous methods have advantages, such as the influence of temperature and initial composition on the SD droplet sizes and the droplet nucleation rate, yet the influence of the substrate is not taken into account. In particular, the wetting condition between the immiscible liquids resulting from phase separation and the solid substrate has rarely been considered. In addition, the lattice-Boltzmann simulations [8] study the phase separation of a symmetric fluid mixture and the amplitude of compositional variations near the substrate. This study investigates the influence of temperature and initial composition on an asymmetric fluid mixture other than the wettability of the bicontinuous structure on solid substrates.
SD improves the properties of polymeric materials, such as mechanical, thermal, chemical, and/or electrical properties [9]. The phase separation can be initiated by three major methods: (a) thermally induced phase separation (TIPS). In this process, SD occurs either by quenching the system from a high to a low temperature for upper critical point phase diagram [10,11], or by heating the system from a low to a high temperature for lower critical point phase diagram [12]. (b) Non-solvent induced phase separation (NIPS). In this approach, the diffusion of the non-solvent species between a non-solvent bath and the polymer solution leads to the passthrough of the spinodal line of the polymer solution [3,13]. (c) Polymerization-induced phase separation (PIPS). In this way, the polymerization of one or more components increases the molecular weight of the components with time. The polymer with a relatively high degree of polymerization becomes immiscible with the solvent, and the phase separation is observed [2,14,15]. In this study, PIPS and TIPS are focused on. We refer to [16,17] for NIPS.
In the current work, the Cahn-Hilliard phase-field method is used to simulate the microstructure evolution and wettability of the polymeric liquid mixture on solid surfaces. In the Cahn-Hilliard approach, a double-well potential is introduced for immiscible fluid phases, and the simulation of phase separation is coupled with the wetting process. The model characterizes the phase state with the composition φ. The time evolution of the parameter φ is such as to reduce the free energy functional of the system, while the volume of the droplet is naturally conserved [18]. This free energy minimization drives the time evolution of the microstructures [19,20].
We model the phase separation of a polymer with an average degree of polymerization N and a solvent based on the Flory-Huggins model. Due to the paucity of the thermodynamic databases for polymer solutions, the Flory-Huggins model describes the free energy density of the system directly deduced from the statistic mechanics and has been quite successfully adopted to predict the behavior of polymer-solvent systems [4,21]. In this study, we focus on the TIPS and PIPS methods of the binary polymer solution. The NIPS method requires thermodynamic database of ternary system [22,23], which is not in the scope of this paper.
When a liquid and a solid are in contact, the liquid droplet spreads/contracts on the solid substrate to reach an equilibrium contact angle known as Young's contact angle. Using Young's law, the contact angle measures the wettability of the droplet on the substrate [24]. Changing the droplet wettability, as well as the contact angle can influence the phase separation process of the polymer solution, which gives rise to distinct material properties, such as the porosity and the volume fraction of the polymer, which will be discussed in the following sections.
In this work, several phenomenological parameters are taken into account, such as wall free energy [25], which represents the attractive and repulsive interaction between the wall and the liquid. In order to determine the wettability of the liquid on the solid wall, the surface free energy is minimized to obtain the boundary condition and to calculate the equilibrium contact angle [25,26]. The calculated contact angle θ, resulting from the boundary condition is consistent with Young's law, which is applied to immiscible fluids. Young's equation accounts for the interfacial tension between the liquids and the solid surface (γ 1 and γ 2 ) and the interfacial tension between the immiscible liquids (σ). The equation is written as [27]:

Modeling
The Flory-Huggins model [4,28] is adopted to depict the free energy density of a polymer mixture with an average degree of polymerization DP = N: The space-and time-dependent variable φ(x, t) depicts the volume composition of the polymer species. The parameter T is the temperature, R g is the universal gas constant, and v m is the molar volume. The first two terms of equation (1) represent the entropy change of the solution, due to the mixing of the polymer with an average degree of polymerization N in the solvent. The last term expresses the latent heat of mixing. The Flory parameter χ characterizes the interaction between the two components in the liquid solution. If the solution of the polymer and the solvent is energetically favored, then χ < 0; if χ > 0, the polymer and solvent interaction is repulsive, then the polymer-polymer and solvent-solvent combinations are energetically favored [29]. To model immiscible liquids in the polymer solution, namely, polymer-rich and polymer-poor phases, the Flory parameter is set to be 1.5 for an average degree of polymerization N = 5. Hence, the free energy density possesses a double-well shape, and the two liquids are separated by an energy barrier, as shown in figure 1.
The present work focuses on a polymer mixture deposited on a solid substrate, where the interaction between the liquid mixture and the substrate has to be taken into account. By adding the interaction between the liquid and the substrate, i.e. the wall free energy, the free energy functional of the system is expressed as: (2) where Ω is the domain occupied by the system and Γ denotes the contact area of the liquid phases with the substrate. The first integration in equation (2) represents the bulk free energy, where f(φ, T) determines the equilibrium composition of the polymer species in the polymer-poor phase φ L1 and the one in the polymer-rich phase φ L2 (see figure 1). The parameter κ scales the gradient energy density, which is determined by the interfacial tension σ of the polymer-rich and polymer-poor phases.
The second integration in equation (2) represents the wall free energy contribution to the system, whereby the liquid composition on the substrate φ s deviates from the equilibrium compositions of the bulk, φ L1 , φ L2 . The wall free energy depicts the attractive and repelling interaction between the liquid and the substrate. This interaction is short-range, which includes the bonding energy and results in the interfacial tension [29,30]. We refer to [31] for the consideration of the longrange interaction.
The Cahn-Hilliard model is adopted for the time evolution of the composition φ within the liquid region Ω of the bulk, which reads: where M defines the mobility, and ξ refers to the applied noise distribution. At times t, . . . , t ′ ∈ [0, t t ] (t t : total time), the Gaussian noise distribution in the system with a stochastic amplitude a is described in a very general stochastic process, in terms of the position vector ⃗ x [32]: At equilibrium, the chemical potential µ equates its equilibrium value µ e (see figure 1) as follows in the entire bulk domain Ω: Differing from the bulk region, the free energy functional on the substrate Γ is assigned as the following formulation: To minimize the free energy functional on the substrate, the following equation is obtained by using the divergence theorem: The vector n denotes the normal vector of the substrate surface. For the time evolution of the composition φ in the substrate region Γ, the following governing equation is adopted: where τ is the kinetic parameter, which controls the dynamic wetting behavior of the system. In the equilibrium state, the free energy functional F s reaches its minimum and the substrate composition φ s fulfills the equilibrium condition as: 2κ ∇φ s · n − ∂γ/∂φ s = 0.
Here, the wall free density γ is written as a second-order polynomial function of the substrate composition φ s [24,29]: where ω 0 , ω 1 , and ω 2 are material parameters. Considering the equilibrium state between the bulk and the substrate regions with the combination of equations (4) Figure 1. (a) The figure represents the free energy density with respect to the composition φ, at the temperature T < Tc; φ L1 , φ L2 are the equilibrium compositions in the two liquid phases, and φsp 1 , φsp 2 denote the spinodal points. (b) Shows the Φ diagram, at a temperature higher than the critical temperature, (c) shows the diagram below Tc, and (d) shows relatively lower temperatures than Tc. φs 1 and φs 2 represent the surface compositions that remain in equilibrium with the L 1 and L 2 phases, respectively. and (6), we determine the surface composition at the flat substrate (n = [0, 1, 0]) as: Multiplying ∂φ/∂y on both sides of equation (7) and integrating it along the y direction, from y to ∞, we obtain: is called the excess free energy density. Substituting equation (2) into equation (6), we obtain the surface composition φ s , which fulfills For the determination of the surface composition, the intersection points between the two curves Φ and Π are evaluated, where Φ := 2 √ κ∆f and Π := ω 1 + ω 2 φ s . Figure 1 illustrates the intersection points between the curves Φ and Π at different temperatures, which are higher and lower than the critical temperature. In figures 1(b)-(d), the two points Φ = 0 depict the equilibrium bulk composition. When the liquid phase is in contact with the substrate (φ s1 and φ s2 ), the intersection points between the two curves of Φ and Π denote the surface compositions. When the substrate is in contact with a φ L1 -rich liquid mixture, the point φ s1 denotes the surface composition. When the substrate is in contact with a φ L2 -rich liquid phase, φ s2 represents the surface compositions. The intersection points enable the evaluation of the interfacial tensions, which determine the wettability of the fluids [29]. Figure 1(b) represents the Φ diagram with respect to the composition and its intersection points at a temperature higher than the critical temperature, which is the maximum temperature at which two immiscible liquid phases transform into a uniform phase (T > T c ). Figure 1(c) shows a temperature below the critical temperature (T < T c ), and figure 1(d) illustrates the intersection points at a temperature that is relatively lower than the critical temperature (T ≪ T c ). At higher temperatures, the system is a homogeneous liquid phase and there is no miscibility gap. At the temperatures below the critical temperature, where the liquid decomposes into two immiscible liquid phases, two surface compositions of the liquid (φ s1 ) and polymer (φ s2 ) phases are determined in contact with the substrate [29].

Surface tension and contact angle
First, we calculate the theoretical liquid-liquid interfacial tension between L 1 − L 2 as: Substituting equation (2) into equation (9) and changing the integral variable from y to φ, by using the chain rule, the value of σ can be written as: Second, the liquid-substrate interfacial tension between L 1 (L 2 ) and the substrate is reckoned as the integration from the substrate position with y = 0 to the L 1 bulk position y = y b1 (L 2 bulk position y = y b2 ), Here, φ s 1 and φ s 2 represent the surface compositions that remain in equilibrium with L 1 and L 2 , respectively. From equation (12), we notice that the liquid-substrate interfacial tension γ i , i = 1, 2, is contributed by the excess free energy, altogether with the wall energy density. The excess free energy represents the free energy between the substrate with the non-uniform composition φ si and the bulk of the liquids φ Li . The excess free energy is similar to the interfacial tension σ, given by equation (10), but with different integral boundaries. Another second contributor, the wall free energy γ(φ si ), is derived by the interaction between the liquid and the solid substrate.
Finally, Young's law is one of the most important methods in this study to investigate the wetting behavior of the fluids on the solid substrate, using the interfacial tension between them. Young's law is written as [30]: The Young's law is achieved via the boundary condition, equation (6).
Evaluating the wetting behavior of the fluids requires analyzing the free energy density of the system and the stability of the phases with respect to thermal fluctuations at each temperature. In figure 1(a), the diagram of the free energy versus the composition is shown at a temperature below the critical temperature. The red and blue points represent the binodal and spinodal points, respectively. According to the spinodal concept, defined by Cahn and Hilliard [6], the solution within the spinodal region is unstable ∂ 2 f/∂φ 2 < 0, whereas the solution outside the spinodal region is stable and ∂ 2 f/∂φ 2 >0. According to equation (1), the entropy contribution of polymer chains to the whole system is inversely proportional to the degree of the polymerization, and increasing N decreases the entropy of the system. When the entropy becomes smaller, enthalpy becomes more pronounced. Enthalpy represents the repulsive interaction between the polymer and the solvent. The repulsive interaction energy between the polymer and the solvent propels the phase separation. Therefore, the larger degree of polymerization makes SD more pronounced. The second derivative of the free energy density f, with respect to the composition at the spinodal points, is zero (∂ 2 f/∂φ 2 = 0). The binodal points are joined with a common tangent, corresponding to the components with homogenized chemical potential at these compositions. Figures 1(b)-(d), the intersection points between Φ and Π, from the modeling section, shows the variation of the free energy density per phase composition at temperatures higher and lower than the critical temperature.
The equilibrium between the polymer-poor phase (L 1 ) and polymer-rich phase (L 2 ) determines the respective equilibrium compositions, namely φ L 1 and φ L 2 , which can be calculated by The spinodal compositions φ sp1 and φ sp2 on the free energy density curve in figure 1 are determined by

Porous microstructure and performances of the polymers
The volume fraction and the porosity of the polymer solution are examined in this study. According to the following equation, the porosity ψ of each microstructure is proportional to the volume fraction of the polymer-rich phase in the microstructure: where φ M is the material volume fraction. In the following chapters, the detailed properties of the polymer solution are described.

Nondimensionalization
In this work, three reference parameters, namely the characteristic length x * = 1 × 10 −6 m, the characteristic surface tension σ * = 1 × 10 −2 N m −1 , and the characteristic diffusivity D * = 1 × 10 −9 m s −2 are chosen to nondimensionalize the simulation parameters (table 1). The governing (equation (3)) for the time evolution of the composition in the bulk region is written as: The nondimensionalized variables are employed to calculate the free energy of the system (equation (1)). The other physical parameters and their scaling factors are defined in the following (table 1): The symbol * tags the scaling factor of the corresponding quantities with the physical unit. The tilde~labels the nondimensionalized physical quantities. In terms of the nondimensionalized parameters, the mass conservation (equation (3)) of the system with a nondimensonalized free energy functional is rewritten as: where the mobility M is assigned with Onsager's relationship as:M The evolution of φ on the substrate is subjected to the following nondimensionalized boundary condition on the substrate (equation (6)):

Contact angle
According to equation (13) and the interfacial energies between the polymer droplet and the substrate, Young's law is applied to the polymer solution by using its equilibrium composition at the appropriate temperature to measure the contact angle. The simulation of  to equation (13). The color bar indicates the polymer composition for all simulation snapshots in this work. The red squares depict the contact angle as a function of the temperature, determined from the resulting microstructures of the simulation, as schematically illustrated in figure 2(f). The calculation method to measure the contact angle of the simulation results is applied according to article [29]. The height h and the base radius r of the droplet are calculated, and the contact angle is expressed as θ = arctan(h/r) × 360/π. Figure 2(e) demonstrates a relatively small deviation (up to 5 • ) for the contact angle between the theoretical values and the simulation results.

Microstructure with different initial composition and temperature
In this section, the computational simulation of the binary systems during SD is performed on 2D domains with a size of 300 × 300. The boundary condition is the same as that mentioned in section 4. Two systems are considered, one containing a polymer with an initial composition of φ 0 = 0.1 and one with φ 0 = 0.2 in a solvent. We consider isothermal conditions for each composition and perform simulations in a temperature range from T = 0.9 to T = 1.2. DP is set to be the constant N = 5 during the simulation. The free energy functional of the system is minimized with respect to the initial composition of the polymer (equation (3)). Figure 3 shows the phase diagram of the polymer-solvent, represented by the dimensionless temperature T versus the polymer composition φ. The red binodal line in the diagram encloses a two-phase region (equation (15) The noise is added to trigger the phase separation in the matrix at the beginning of the simulation. In later stages, when the polymer-rich and polymer-lean phases reach the equilibrium composition, interfacial energy minimization dominates the microstructure evolution in lieu of the phase separation due to the thermal noise [33]. In addition, as demonstrated in [34], the noise effect does not affect the scaling law of the domain growth in the process of SD.
For an initial composition of φ 0 = 0.1, we simulate five cases for two scenarios: (a) outside the spinodal region, with the respective temperatures of T = 1.2, T = 1.05, and T = 1.0, and (b) inside the spinodal region, with lower temperatures of T = 0.95 and T = 0.9 ( figure 3). In addition, the effect of the substrate on SD is also examined at each temperature, once with and once without a substrate, with the same input parameters. In the simulations, the influence of substrate surface roughness is ignored. The surface roughness leads to advancing and receding contact angles and requires large-scale simulations [35], and it is not in the scope of this study.
To validate the phase diagram in figure 3(a 1 )-(e 1 ), we first perform the polymeric phase separation without the presence of the substrate. As illustrated in the rightmost panel of figure 3, with the homogeneous initial polymer composition φ 0 = 0.1, SD occurs only when the temperature T = 0.9 falls below the dotted blue spinodal line in the phase diagram. At other higher temperatures, the thermal fluctuation is dissipated due to the absence of phase separation driving force and therefore does not lead to the formation of polymer droplets.
Next, in the presence of a solid substrate (gray), we perform the same simulations as in the previous cases. As shown in figure 3(a 2 )-(e 2 ), no polymer droplets are produced at T > 0.9. However, we observe the formation of polymer layers that come into direct contact with the substrate. With descending temperature, the thickness of the polymer film increases, and its polymer composition rises up, which is manifested by its gradually reddish color. The occurrence of this layer is attributed to the local equilibrium on the substrate, as described by equation (8). In contrast to the mass equilibrium, this surface equilibrium composition on the substrate is determined by the bulk free energy and the wall free energy, expressed in figures 1(b)-(d) by the right intersection points φ s 2 . At higher temperatures, the existence of this intersection point gives rise to the aggregation of polymer species on the substrate, causing the surface composition. With decreasing temperature, φ s 2 becomes larger, and the polymer layer on the substrate turns to be thicker. Besides, in the snapshots of figure 3, the blue gradient above the substrate indicates the origin of the polymer layers. As the wall energy reduces the free energy of the system, polymer species remain preferential and accumulate on the substrate. Moreover, the surface tension effect comes into play at lower temperatures, and the polymer layer breaks into separate droplets to minimize the surface energy.
In figure 4, we present the details of the temporal evolution of the domains in the presence of the substrate and the selected temperature range. The color bar indicates the polymer composition in the whole domain, where the droplets represent the polymer-rich phases. At T = 1.2, the time evolution illustrates the formation of a layer of the polymerrich phase at t = 0.4 s reaches the surface composition. At T = 1.05, the first two time stages show a composition gradient between the substrate and the bulk until the substrate reaches the equilibrium surface composition. This process towards equilibrium leads to the formation of a layer with a high polymer composition on the surface. The subsequent time evolution t = 0.4 s displays a thick polymer layer on the substrate, in contrast to the thin polymer layer at T = 1.2. The thickness of the polymer layer depends on the surface composition upon the temperature. At the same time, the bulk composition is not saturated enough to promote SD in the bulk. As the temperature decreases to T = 1.0, the surface composition moves inside the spinodal region, and SD occurs. At t = 0.04 s, SD occurs spontaneously on the substrate, and the droplets grow larger over time. According to the phase diagram (figure 3), the bulk composition at T = 1.0 is outside the spinodal region; hence there is no SD in the bulk.
At T = 0.95, the bulk composition does not show evidence of SD, yet the phase diagram is considered for the system without substrate. The system decreases the free energy by forming polymer droplets on the substrate. At T = 0.9 and t = 0.02 s, SD spontaneously occurs on the substrate as a result of the surface-directed SD, and at t = 0.04 s, the number of droplets increases [36]. At T = 0.9, the local composition gradient causes a delay in the formation of SD in bulk, which As the temperature decreases, the number of droplets evidently increases in the domains. Especially at later times, more droplets tend to form when reducing the temperature. Figure 7(a) indicates that the volume fraction of the lowest temperature maintains a more considerable increment than higher temperatures. The phase separation occurs mainly at low temperatures, which is consistent with the literature (equation (16)) [37,38].
The two phenomena, Ostwald ripening and coalescence are observed in figure 4. The merging between the two small droplets and a neck, known as droplet coalescence, results in the formation of an enlarged polymer droplet. In the process of Ostwald ripening, the big droplet grows at the expense of the small droplet. This kind of growth is caused by the dependence of the solubility on the mean curvature of the droplet, known as the Gibbs-Thomson effect. The high composition of the small droplet and the low composition of the big droplet lead to a diffusion flux between the droplets, causing the growth of the big droplet [39]. At T = 0.9 and t = 0.04 s, a number of small droplets form on the substrate. At t = 0.2 s, the number of droplets decreases, and the radii of the droplets on the substrate increase due to coalescence and Ostwald ripening. Fig, the microstructures at t = 0.2 s show that with decreasing temperature, the effect of coalescence on the substrate is more pronounced [40].

With different temperatures in
The simulation results of φ 0 = 0.2 (figure 5) show that SD forms earlier and faster compared to φ 0 = 0.1, as temperature decreases. This is because of the higher driving force, compared to that of φ 0 = 0.1. In figure 5, it is also noticeable that the SD droplets at higher temperatures look blurry, and the polymer-rich phases have a lower polymer composition than at low temperatures, where the droplets look brighter. Figure 6(a) shows at the beginning of the simulation that droplet formation occurs more frequently at low temperatures than at higher temperatures. Subsequently, the droplet sizes decrease with time. Moreover, figure 6(b) shows that the mean radius of the formed droplets varies at the considered temperatures in the first 0.05 seconds of the simulation. In the first 0.05 s, the mentioned figure illustrates a peak of droplet growth in isothermal simulations. The mean radius of the formed droplets increases with the temperature until T = 1.2, which reaches the highest mean radius 13.8 µm of the five considered temperatures. Subsequently, until the end of the simulation, the droplets in all temperatures grow monotonically and approach a similar size. Figure 5 shows that the polymer droplets are smaller, denser, and more compact as the temperature decreases. At the initial stages t < 0.04 s, the polymer droplets occur earlier at lower temperatures than at higher temperatures. This fast formation of polymer droplets is due to the large driving force which is proportional to |∂ 2 f/∂φ 2 | and the high mobility M which is expressed as (v m /R g )D(φ 0 (1 − φ 0 )/T (equation (18))) at low temperatures. The high mobility at low temperature results from the assumption of a constant diffusion coefficient D. A temperature dependent D may lead to a different relation between mobility and temperature. As the time evolves, the composition reaches the equilibrium values. At this stage, mobility decreases with a decrease in the temperature, attributing to the dependence of the equilibrium composition φ L1 on the temperature. This behavior outside the spinodal region is contrary to the mobility-temperature relation inside the spinodal region and consistent with literature [41]. The small mobility at low temperature leads to slow Ostwald ripening and sluggish coalescence, resulting in a compact microstructure with small droplets (t = 0.2 s and 0.4 s). Figure 6 shows that at T = 1.2, droplet formation starts later than at the lower temperatures, which corresponds with figure 5, which shows that no droplets form in the first two time steps at this temperature.
In the polymer solution with the initial composition of φ 0 = 0.2, the time evolution of the SD formation shows some phenomena such as Ostwald ripening and coalescence. Both phenomena are observed in figure 5 (highlighted with green squares), especially at the temperature T = 0.95. At T = 1.05, the first figure at t = 0.02 s, the worm-like structure, shows the formation of SD. Especially near the substrate, the droplets form a regular orientation to the substrate, which is more pronounced than the rest of the droplets. At the next time evolution, t = 0.04 s, a merger of the two small droplets occurs with a neck, known as droplet coalescence, forming an enlarged polymer droplet. The droplet coalescence in polymer solutions depends on the composition of polymer droplets in a solvent. Figures 7(a) and (b) show the time evolution of the volume fraction of the polymer droplet at a different initial composition of φ 0 = 0.1 and 0.2, respectively. It can be clearly seen that for φ 0 = 0.1, only the case with T = 0.9 enters the spinodal region. Therefore, the volume fraction of the  polymer droplet continues to increase with phase separation. When T > 0.9, the small increment of the volume fraction is attributed to the formation of the polymer film on the substrate, as discussed at the beginning of this section. For the higher initial polymer composition φ 0 = 0.2, the phase diagram shows that all temperature setups enter the spinodal region and that the final polymer volume fraction converges to the same value 23%. In addition, a special delay in the polymer droplet production is observed in T = 1.2. This can be explained by a weaker driving force for phase separation at high temperatures.

Temperature gradient in the domain
In this section, we consider the influence of the temperature gradient on the microstructure evolution. The simulations of the binary systems are performed on a 2D domain with a size of 300 × 900. A linear temperature gradient is applied in the y direction, perpendicular to the substrate with an initial composition of φ 0 = 0.1. The temperature on the substrate (y = 0) is set as T sub = 0.8 and increases in the y dimension until it reaches the preset value T top at the top of the domain.
As a comparison for an isothermal condition T = 0.8 without temperature gradient, figure 8(a) shows that the SD droplets with different sizes are dispersed over the entire domain. In figure 8(b), the temperature increases linearly in the y direction, until it reaches T top = 0.9. Nevertheless, SD still occurs throughout the region, where the droplets have a larger average size due to the temperature rise. A similar tendency is observed in figure 8(c), where the droplets grow bigger as the temperature increases. However, since the peak temperature T top = 1.0 is above the spinodal region, no polymer droplets are generated in the upper region. As the temperature gradient increases, the droplet size distribution also becomes more irregular. The droplet radius is much larger in the high temperature region than in the lower cold area. Close to the substrate, this causes the volume fraction of the polymer to change with distance from the substrate, which will be discussed later. Moreover, figures 8(d) and (e) show that the influence of the temperature gradient on the droplet morphology is more evident.
To compare the microstructure evolution between the two scenarios with and without substrate, the same linear temperature gradient as in figure 8 is applied to the region without substrate.    shows that the temperature gradient has a considerable effect on the phase separation. The effects are related to the distribution of polymer-rich droplets over the domains and their sizes. It is clearly seen that the temperature gradient affects the volume fraction, porosity, and polymer performance. In the following, the mentioned properties are studied individually.  To understand the role of the temperature gradient in the development of the microstructure, the domains are evenly divided into six subdomains. In each subdomain, different properties such as the number of droplets, the mean droplet radius, and the porosity are calculated. Figure 10 shows the domain sectioning for two representative simulations with substrate (figures 8(a) and (c)). Domain A refers to an isothermal simulation at T = 0.8 and domain B visualizes the droplet distribution for a temperature gradient of T grad = 0.2. In the domain with the constant temperature, the number of droplets grows so fast in the first subdomain near the substrate that it takes a few seconds for the droplets in the upper subdomains to nucleate. As mentioned in the last section, the high mobility in the spinodal region causes the substrate to absorb polymers and spontaneously reach the surface composition. The composition gradient near the substrate causes the approach of the bulk composition towards the spinodal composition, leading to the formation of SD. The number of droplets in the first subdomain indicates multiple phenomena, such as a declination around t = 0.05 s. Then the number of droplets increases until t = 0.15 s. Figure 11 display the microstructure evolution in subdomain 1 from figure 10. In figure 11, the droplets in subdomain 1 only form at a certain distance from the substrate. The large droplets straight at the substrate cause a composition gradient, thereby the droplets form in bulk at a distance from the substrate. Figure 11 explains the fluctuation of the number of droplets in the first subdomain, where the green arrows at t = 0.08 s refer to the coalescence of the droplets on the substrate at t = 0.12 s. The green circles at t = 0.2 s represent the shrinkage of the droplets due to Ostwald ripening, which causes the droplets around them to grow larger. These two phenomena lead to a sudden increase of the droplet sizes in subdomain 1 (figure 10). The mean radius of the droplets in the same subdomain becomes larger as the number of droplets increases until the droplets reach a maximum mean radius in the time considered. Except for subdomain 1, the number of droplets in each subdomain increases monotonically until a certain value is reached after a certain time. Over time, the porosity of the domain decreases simultaneously for all subdomains. There is no significant deviation between the porosity curves.
The set of diagrams in B is related to the domain with a linear temperature gradient of T top = 1.0. In this domain, the number of droplets decreases monotonically from the lower to the upper subdomains. No droplets form in the last two subdomains (5 and 6). This observation shows that the number of droplets decreases with increasing temperature. In the later stage, the mean radius of the droplets in the subdomains with high temperature is larger than in the subdomains with low temperature. The effect of temperature on the SD droplets is that the high temperature suppresses the droplet formation, thus, the droplets are more likely to form in the later time steps.
For all subdomains, the porosity decreases with time. The high temperature delays the porosity initiation and decreases the porosity rate with time. Figure 12 shows two domains without substrate, corresponding to figures 9(a) and (c). One domain represents an isothermal simulation at T = 0.8 (A), while the other one represents a simulation with a temperature gradient of 0.2 (B). Between the subdomains of the isothermal domain, there is no significant variation in the number of droplets. The time evolution of the mean droplet radius and porosity also shows the same behavior. Two domains illustrated in figures 10(A) and 12(A) show quite similar quantities for the number and mean radius of the droplets, and the porosity, except for subdomain 1. In subdomain 1, the interaction between the liquid and substrate takes place and accelerates the formation of SD droplets near the substrate.
In the domain B, droplets are formed asynchronously in each subdomain. In the fourth subdomain, it takes longer for droplets to form, and the number of droplets is much lower than in the first three subdomains. According to the applied temperature gradient, the fourth subdomain has a relatively high temperature, which leads to larger droplet size. In the domain without substrate, the highest number of droplets ( figure 12(B)) is lower than in the domain with the substrate ( figure 10(B)), and the substrate leads to a higher number of droplets with small sizes. The highest number of droplets belongs to the first subdomain with a relatively low temperature. With time evolution, increasing the temperature leads to the formation of fewer droplets with larger sizes. The temperature influences the porosity, and the increasing temperature causes a declination rate of the porosity with time compared to relatively low temperatures. Also, the microstructure in the first subdomain for the two configurations with ( figure 10(B)) and without substrate ( figure 12(B)) show that the substrate forces the droplets to form earlier.

Conclusion
In summary, this study presents the time evolution of the porous microstructures in the presence of substrate in comparable conditions without substrate. In the presence of the substrate, the substrate accelerates the formations of SD in such a way that the substrate reaches the surface composition earlier than the same condition without substrate. This study also shows that lowering the temperature leads to a shift in the surface composition towards the spinodal region. This process results in the formation of droplets on the substrate, even if the bulk composition is outside the spinodal region. In isothermal simulations of a binary polymer solution, at low temperatures inside the spinodal region, the volume fraction of the polymer-rich phase increases faster with time. The Cahn model equation (4) [6] has been used to simulate the wetting phenomenon, where SD occurs in immiscible liquids until the liquid reaches its minimum energy at equilibrium. The equilibrium contact angles of the formed droplets and their wetting behavior on the solid substrate are evaluated. Using Young's law, the interfacial tension between the substrate and the droplet determines that the resulting contact angle is consistent with the theory.
In a further study, we investigate the effect of linear temperature gradients on the formation of SD in the direction perpendicular to the substrate. By fixing the temperature on the substrate, the temperature gradient in the y direction significantly influences the number of polymer-rich droplets, their size, and the porosity of the domains. It is also shown that in the domains near the substrate, the substrate accelerates the formation of droplets. Therefore, spinodal regions at low temperatures are more favorable for SD droplet formation than at high temperatures, where SD is suppressed.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.