Version of Record: https://www.sciencedirect.

com/science/article/pii/S0376738821007547
Manuscript_50082ca1cbdf849dbd910d17c42933aa

1 A review on models and simulations of membrane formation via phase

2 inversion processes
3
4 Yuanhui Tang1,2,3, Yakai Lin2*, David M. Ford3*, Xianghong Qian4, M. Rosario Cervellere5, Paul C.
5 Millett5, Xiaolin Wang2*
1
6 School of Chemical and Environmental Engineering, China University of Mining and Technology, Beijing
7 100083, P.R.China
2
8 Beijing Key Laboratory of Membrane Materials and Engineering, Department of Chemical Engineering,
9 Tsinghua University, Beijing, China
3
10 Ralph E. Martin Department of Chemical Engineering, University of Arkansas, USA
4
11 Department of Biomedical Engineering, University of Arkansas, USA
5
12 Department of Mechanical Engineering, University of Arkansas, USA
13
*
14 Corresponding author: Xiaolin Wang: xl-wang@tsinghua.edu.cn, David M. Ford: daveford@uark.edu and
15 Yakai Lin: yk_lin@tsinghua.edu.cn

16 Abstract

17 Phase inversion processes, in which a polymer is transformed from a thin solution film into a solid matrix
18 upon exposure to another gas/liquid phase or a temperature change, represent the most widely used
19 manufacturing technique for polymeric membranes. However, the level of understanding of the connection
20 between process variables and final membrane structure remains surprisingly empirical. Modeling efforts have
21 largely been based on heuristics that correlate experimentally observed pore structures with process conditions,
22 using thermodynamic phase diagrams and macroscopic transport models. Recent developments have added
23 mesoscopic phase field approaches and molecular simulations to the array of modeling techniques applied to
24 this problem. These latest developments hold the promise of delivering accurate, directly visualized
25 representations of membrane pore structures formed by a given process. This article firstly provides an overview
26 of the capabilities and limitations of the different modeling techniques and then focuses on how they have been
27 developed and applied over the last 40 years to aid our understanding of membrane formation via phase inversion
28 processes, especially nonsolvent induced phase separation (wet-casting), vapor/evaporation induced phase
29 separation (dry-casting), and thermally-driven. Meanwhile, current challenges and future prospects, such as
30 linking techniques across length and time scales as well as capturing the detailed effects of polymer
31 crystallization, are also discussed. This review aims to offer inspiration for further progress in the field of
32 modeling development for membrane preparation by phase inversion processes.
33
34 Keywords: Modeling; Membrane formation; Simulation; Phase inversion; Polymer solution;
35

36 Contents
37 Abstract .................................................................................................................................................. 1
38 1. Introduction ....................................................................................................................................... 2
39 2. Model and simulation method ............................................................................................................ 4
40 2.1. Macroscopic transport models .................................................................................................. 6
41 2.2.1 NIPS (wet-casting process): isothermal precipitation ....................................................... 7
42 2.2.2 VIPS (dry-casting): coupled heat and mass transfer ......................................................... 8
43 2.2.3 TIPS ............................................................................................................................. 10
44 2.2. Mesoscopic phase-field (PF) approach .................................................................................... 11
45 2.3. Molecule/particle-based simulations ....................................................................................... 12
46 2.3.1 Molecular dynamics (MD) ............................................................................................ 12
47 2.3.2 Dissipative particle dynamics (DPD) ............................................................................. 13
48 2.3.3 Monte Carlo (MC) ........................................................................................................ 14
49 3. Developments and applications of the dynamic modeling tools ......................................................... 15
50 3.1. Macroscopic transport models ................................................................................................ 15

© 2021 published by Elsevier. This manuscript is made available under the Elsevier user license
https://www.elsevier.com/open-access/userlicense/1.0/
 1 3.1.1 NIPS (wet-casting) ........................................................................................................ 16
2 3.1.2 VIPS (dry-casting) ........................................................................................................ 17
3 3.1.3 TIPS ............................................................................................................................. 19
4 3.2. Mesoscopic PF approach ........................................................................................................ 20
5 3.3. Molecule/particle-based simulations ....................................................................................... 21
6 3.3.1 MD ............................................................................................................................... 21
7 3.3.2 DPD ............................................................................................................................. 23
8 3.3.3 MC ............................................................................................................................... 24
9 3.4. Challenges and outlook........................................................................................................... 25
10 4. Conclusions ..................................................................................................................................... 26
11 References ............................................................................................................................................ 27
12

13 1. Introduction

14 Phase inversion (PI) processes, in which the thermodynamic state of a homogeneous polymer solution is
15 changed by contacting with another phase (liquid or vapor), and possibly changing temperature as well, to
16 promote the formation of a solid phase, are the most commonly used techniques for manufacturing microporous
17 polymeric membranes [1-8]. Depending on the state and property of the contacting phase, PI processes can be
18 divided into nonsolvent induced phase separation (NIPS) [2] (more commonly wet casting, or immersion
19 precipitation), vapor/evaporation induced phase separation (VIPS) (more commonly dry casting) [4] and
20 thermally induced phase separation (TIPS) [6]. As described in recent review papers [8, 9], for commonly
21 employed polymers such as cellulose acetate (CA), polyethersulfone (PES), polysulfone (PSf), polyvinylidene
22 fluoride (PVDF), the resulting membrane structure can range from nearly symmetric, as in microfiltration
23 membranes, to highly asymmetric, as in the integrally skinned nanofiltration and reverse osmosis membranes.
24 A long-term goal in this field is to be able to engineer the system and process to achieve the desired structure
25 and pore size distribution in the resulting membrane. Key variables in the PI processes, including solvent and
26 nonsolvent species, compositions of the phases, and temperatures of the polymer solution and the contacting
27 phase, must be carefully considered and controlled to achieve this goal [10-14].
28 Membrane formation is a complex multiscale process governed by thermodynamic potential differences,
29 heat, and mass transfer rates, and the kinetics of phase formation. For simplicity, one might consider three key
30 scales. At the molecular scale, solvent molecules and polymer segments interact and move on the nanometer
31 (nm) and nanosecond (ns) length and time scales, respectively. At the mesoscale, phase domains and other
32 ordered structures with characteristic length scales of 10~100 nm form over time scales of 10~1,000 ns. At the
33 production scale, the membrane is normally several micrometers (µm) thick with pore sizes on the order of
34 10~1000 nm, and the formation operation may take several seconds (s) or even hours [15-20]. In the recent two
35 decades, various techniques have been adopted to observe and characterize the phase separation, the domain
36 evolution, the membrane thickness, and the chemical constitution evolution via the PI processes, including
37 optical microscopy equipped with a temperature controller [16, 20-25], water contact angle measurements [22],
38 Raman confocal spectrometry [22], Fourier transform infrared spectroscopy microscopy [22], laser displacement
39 sensor [26, 27], and light-scattering measurements [21, 28]. However, considering the large number of variables,
40 process speed, and complexity of the interactions, PI processes are difficult to observe and characterize via
41 experiment methods alone. Therefore, modeling is a particularly useful tool to significantly increase our
42 understanding of membrane formation.
43 Firstly, one can easily understand that the membrane formation induced by energy variation and nonsolvent
44 intrusion should involve three closely related transport phenomena: fluid dynamics, heat transfer, and mass
45 transfer. In physics and engineering science, the three transport phenomena normally occur simultaneously and
46 the mathematical tools for describing them are similar [29]. The bases for the mathematical models are the
47 fundamental physical and chemical laws, such as the laws of conservation of mass, energy, and momentum,
48 which are normally expressed by different kinds of continuity equations that can be regarded as a stronger, local
49 form of the conservation laws [30]. The continuity equations underlie more specific transport equations such as
50 the convection-diffusion equation[31], Boltzmann transport equation[32, 33], and Navier–Stokes equation [34].
51 Necessary assumptions should be made to realize an engineering compromise between a rigorous description

2
 1 and getting an answer that is good enough. Combining the specific transport equations, essential assumptions,
2 and the specific characteristics of the membrane formation, appropriate mathematic models can be established
3 to obtain the system concentration profiles [35]. The most useful results of these transport models are so-called
4 ‘‘composition’’ or ‘‘mass-transfer’’ paths that are time-dependent concentration profiles plotted directly on
5 phase diagrams of the polymeric/solvent/nonsolvent systems, and are used to predict (i) the delay time between
6 immersion and phase separation and (ii) the polymer concentration profile in the film during the initial demixing.
7 Until recently, models that were much closer to real temporal-spatial scale and conditions have been gradually
8 developed to provide more essential results and predictions [36], although the results of transport models are
9 limited in directly reflecting the phase separation morphologies and membrane microstructure.
10 Secondly, because the membrane formation via PI involves phase separation and polymer solidification,
11 some researchers believed that the phase field (PF) approach that was initially proposed to study the phase
12 separation and phase inversion, would be quite powerful in membrane formation modeling. The original idea of
13 the PF approach can be treated back to Gibbs who implicitly entertained the idea of a diffuse interface, which
14 pointed that any density of an extensive quantity (e.g., the mass density) between two coexisting phases varied
15 smoothly from its value in one phase to its value in the other [37]. Then focusing on the transition phenomena,
16 Landau employed this notion and firstly introduced an additional parameter to label the different phases in their
17 theory on the absorption of liquid helium [38, 39]. Subsequently based on the variational principles of
18 irreversible thermodynamics as founded by Onsager [40, 41] and the thermodynamic theory of capillarity
19 derived by van der Waals [42], Cahn and collaborators rigorously developed a nonlinear diffusion theory,
20 considering the phase separation process of a "nonuniform" system and then obtained the well-known Cahn-
21 Hilliard and Allen–Cahn equations [43-48]. Within their contribution to the field, the PF model is extremely
22 powerful in obtaining the kinetics of the separation and the morphology of the separated phases in real fluids,
23 which included applications with emphasis on special processes that lead to commercially interesting materials
24 such as micro-dispersions and membranes. Compared to the macroscopic transport models, PF methods are more
25 intuitive to produce simulations of the membrane formation processes at the temporal-spatial scale of several
26 micrometers and seconds and predict a variety of structural features that are qualitatively similar to experimental
27 observations.
28 Thirdly, since the component selection and molecular structure play significant roles in the membrane
29 formation via PI, and the aforementioned two methods (transfer and PF models) treat atoms and molecules
30 implicitly, some people tended to adopt molecular-scale or particle-based simulations to study the membrane
31 formation. The first computer simulation of liquid was carried out at the Los Alamos national laboratories in
32 1953 [49]. Nowadays with the rapid development of computer hardware, the field of computer simulation has
33 enjoyed rapid advances. The molecular Molecular dynamics (MD), dissipative particle dynamics (DPD), and
34 Monte Carlo (MC) are three molecular/particle-based approaches that have commonly been employed to
35 simulate the membrane formation via the PI processes. Considering the spatiotemporal scale of an MD
36 simulation is only limited on the nanoscale, now rather than being used to study the earliest stages of the phase
37 separation during the membrane formation, MD is more successful to be an analytical tool to study the affinity
38 and diffusion between components. DPD can be regarded as a coarse-grained MD simulation method while it
39 employs particles that substantially represent a set of several adjacent repeat units along a polymer chain or a
40 group of several solvent molecules [50]. DPD can access length and time scales at or beyond the μm and μs
41 levels, which makes it suitable to simulate the phase separation process and study the membrane pore formation.
42 Unlike MD and DPD, MC simulation is built on probabilistic concepts developed from the ensemble theory of
43 statistical mechanics, and no velocity variables are incorporated. The spatiotemporal scale of an MC simulation
44 is comparable to that in DPD since it is ultimately limited by the ability to identify the set of particles (or
45 configuration space variables).
46 Since the 1970s, macroscopic transport models, mesoscopic PF models, and molecular/particle-scale
47 simulations have been proposed to describe the transfer phenomena and membrane formation in the evaporation
48 and precipitation steps of the PI processes. The macroscopic transport models and the PF models represent “top-
49 down” approaches, while the molecular-scale simulations, which include molecular dynamics (MD), dissipative
50 particle dynamics (DPD), and Monte Carlo (MC), are “bottom-up” in the sense that the system behavior is
51 emergent from a set of particle interaction potentials [51-55]. This is why DPD and MC are classified to the
52 category of molecular-scale simulation, although they have much larger or even mesoscale spatiotemporal
53 scales. This article aims to present a full and comprehensive review of the literature on the macro-, meso- and

3
 1 molecular-scale modeling approaches of membrane formation via the phase inversion processes. Based on an
2 introduction to classical thermodynamic diagrams of the polymer solution system, Section 2 provides an
3 overview of the dynamic models at the three different scales. Detailed equations and results of several key
4 models are presented to give readers a basic understanding of their predictive capabilities and limitations. In
5 section 3 we focus more on how the models have been developed over time and impacted membrane formation
6 research. Finally, outstanding challenges and conclusions are presented.
7 2. Model and simulation method
8 Thermodynamic phase diagrams have traditionally been the primary source of guidance for understanding
9 the outcome of the phase inversion processes. Such diagrams can be determined directly from experimental
10 techniques like cloud-point [56-61] and crystallization temperature measurements [62-64] or computed from
11 free energy models like Flory-Huggins (FH) theory [65-70].

12
13 Fig.1 A temperature-composition phase diagram for a binary polymer-solvent system, relevant to the TIPS
14 processes. Thermal quenching through different regions of the phase diagram (Paths A, B, and C) yields different
15 porous structures. Note that the ends of the tie lines show the compositions of the two phases that exist in
16 equilibrium with each other at this temperature.
17 The simplest example is a temperature-composition phase diagram for a binary polymer-solvent system, as
18 shown schematically in Fig. 1. Of the different PI processes discussed in Section 1, only TIPS can be
19 comprehensively described by this type of diagram, but it serves as a good starting point for understanding more
20 complex situations. The idea here is that thermal quenches along different vertical cooling paths will produce
21 different membrane morphologies due to the influence of thermodynamics. Three different possibilities can be
22 described, corresponding to the three arrows in Fig. 1 [6]. Passing through the liquid-liquid equilibrium (LLE)
23 spinodal curve by quenching or a rapid cooling before entering the solid-liquid equilibrium (SLE) region, as
24 shown in Path (A), leads to bi-continuous, open-pore morphologies due to the process of spinodal decomposition
25 in which rapid local phase separation results in the formation of liquid-phase microstructures that coarsen over
26 time and govern the ultimate pore structure. Passing through the LLE binodal curve by a moderate cooling, as
27 shown in Path (B), will lead to cellular morphologies due to the process of nucleation and growth of a polymer-
28 lean liquid phase; the rate at which these solvent-rich nuclei grow and coalesce, as determined by their inherent
29 kinetics and the rate at which the polymer solidifies, will determine the size and openness of the cellular structure.
30 Finally, passing through the SLE curve, as shown in Path (C), leads to semi-crystalline structures due to the
31 process of nucleation and growth of a solid phase out of the liquid; the pore structure is governed by the size and
32 connectedness of the solid polymer particles, with the most common particle geometry being spherulitic.
33 To address NIPS (including VIPS), a phase diagram that includes at least three components - polymer,
34 solvent, and nonsolvent - is necessary. For systems that involve a fourth additive component, such as a
35 nonsolvent or an additive polymer, a tetrahedral phase diagram should be considered [71-73]. A schematic of a

4
 1 typical ternary phase diagram for a ternary system including a semicrystalline polymer, solvent, and nonsolvent,
2 adapted from Ref. [1], is shown in Fig. 2 (A). This type of phase diagram can help to understand and quantify
3 the phases that are formed when a liquid polymer-solvent solution (somewhere between the solvent vertex and
4 Point P) is mixed with a nonsolvent in a certain ratio, as in a NIPS process. Except for the homogeneous region
5 (I) where all components are miscible, and a region where liquid-liquid demixing occurs (II), other phase
6 separations (III, IV, and V) can also be observed. Usually, the focus is on the LLE region (II) and the
7 crystallization curve (PQ). A composition somewhere in the region of P-Q-polymer will contain a crystalline
8 pure polymer that is in equilibrium with a composition somewhere on the crystallization line PQ. However, a
9 real NIPS process is somewhat more complex since solid-liquid demixing occurs in addition to liquid-liquid
10 demixing. One can imagine studying a hybrid TIPS/NIPS process by adding temperature as an out-of-plane axis,
11 thus creating a three-dimensional phase diagram by superposition of a series of isothermal ternary systems, as
12 sketched in Fig. 2(B) from Ref. [74].

13
14 Fig. 2 (A) A schematic phase diagram for a ternary polymer-solvent-nonsolvent system at a fixed temperature,
15 relevant to the NIPS process. This is a general example of a semicrystalline polymer that includes the possibility
16 of five different regions: homogeneous liquid (I), liquid-liquid coexistence (II), solid-liquid (S-L) coexistence
17 (III), solid-liquid-liquid coexistence (IV), and another S-L coexistence (V). Other possible features include gels
18 or glasses, although these are not true thermodynamic phases. Adapted from Ref. [1] with permission. (B) A
19 schematic phase diagram for a ternary polymer-solvent-nonsolvent system at different temperatures. Reproduced
20 from Ref.[74] with permission.
21 As mentioned above, free energy models like the FH theory can be used to predict the LLE portions of
22 phase diagrams like the ones shown above by applying the thermodynamic principles associated with the
23 spinodal (stability limit) and binodal (equal chemical potential) curves. There are many resources describing the
24 procedures for computing the liquid-liquid coexistence part of a binary phase diagram, like Fig. 1 [65, 66], or a
25 ternary phase diagram, like Fig. 2(A), from the FH model [75-78]. The solid-liquid (S-L) curves in these phase
26 diagrams mainly depend on the experimental detection due to the complexity of the polymer crystallization and
27 solidification. Although recently, some computationally complex models like lattice cluster theory have been
28 employed to compute the phase diagram [79].
29 For the systems consisting of more than two components, the idea here is that system composition changes
30 along different composition profiles (transfer paths) will produce different membrane morphologies. Schematic
31 examples are shown in Fig. 3. In the case of Fig. 3(A), the composition profile crosses into the binodal region,
32 and probably inside the spinodal region, which would lead to an instantaneous demixing. This is likely to produce
33 a membrane with a highly porous substructure (with finger-like macrovoids) and a thin skin top layer. In the
34 case of Fig. 3(B), the composition profile doesn’t cross the binodal line, so the mechanism of demixing is
35 different and membrane formation would take much longer. This is likely to produce a membrane with a
36 comparatively dense top layer and sponge-like substructure.

5
 1
2 Fig. 3 Sketches of two different concentration profiles from the top (T) to bottom (B) of a casting film
3 immediately after immersion into a coagulation bath and undergoing NIPS processes, generated from a transport
4 model and superimposed on the ternary phase diagram of the polymer-solvent-nonsolvent system. Adapted from
5 Ref.[2] with permission
6 Although thermodynamic phase diagrams have been regarded as a key ingredient for understanding the
7 phase inversion process, one must be aware that such diagrams represent only the equilibrium properties of a
8 system, while the membrane formation is always quite far from an equilibrium process. The resultant membrane
9 structure is highly dependent on dynamic factors like cooling rates, demixing rates, aggregation kinetics, gelation,
10 and solidification. In the next subsection, a brief overview will be presented to explain how the aforementioned
11 thermodynamic phase diagrams are integrated with different dynamic modeling tools to predict and explain
12 observed phenomena in membrane formation.
13 2.1. Macroscopic transport models
14 Since the 1970s, a series of macroscopic transport models have been constructed to explore the mechanism
15 of microporous structure formation via the PI processes. Figure 4 displays schematic representations of the
16 casting polymer solution and the surrounding “medium” for the different PI processes. According to different
17 phase states of the bath and transport phenomena involved, the models can be divided into three categories,
18 including the NIPS (wet-casting), VIPS (dry-casting), and TIPS process as reviewed below. Here a short
19 introduction is given to describe the characteristics of the three processes. In the wet-casting process, there is a
20 polymer dissolved in a solvent (sometimes with a small amount of nonsolvent) to form a solution facing a liquid
21 nonsolvent bath. Because normally the nonsolvent and the solvent are miscible, mass transfer of these two
22 components would occur between the polymer solution and bath. Accordingly, in the transport model, normally
23 only solvent and nonsolvent fluxes need to be considered because the polymer’s fraction was simply determined
24 by conservation, and for the casting solution the polymer concentration is normally none due to its relatively low
25 diffusion. This process is usually modeled as isothermal, as the solution and the nonsolvent are assumed to be
26 initially at the same temperature, and the mass transfer between the two liquid phases is not assumed to have an
27 obvious thermal effect. The dry-casting process is a little complex and actually, it refers to two processes: vapor
28 induced phase separation and evaporation induced phase separation. The former is a process that a casting
29 solution film consisting of a polymer and a solvent, is placed in a vapor nonsolvent phase which may initially
30 contain solvent to avoid too fast solvent evaporation from the solution. Membrane formation occurs because of
31 the penetration (diffusion) of nonsolvent into the cast film; While the latter is a process that the polymer is
32 dissolved in a mixture of solvent and nonsolvent where the solvent is more volatile than the nonsolvent. The
33 composition shifts during evaporation to a higher nonsolvent and polymer content, leading to polymer
34 precipitation and the formation of a denser membrane. Whereas considering the two processes from the
35 perspective of transport modeling or simulation, it would be reasonable to find that the calculation system is
36 similar: The polymer/solvent solution (sometimes with a small amount of nonsolvent) is facing a vapor phase,
37 so there is solvent (and possibly nonsolvent) evaporating into the gas phase or nonsolvent penetrating the casting
38 solution. The evaporation process has a considerable thermal effect, and hence coupled heat and mass transfer

6
 1 should be considered in the transport model. We note that this situation also applies during solvent evaporation
2 steps of wet-casting and TIPS. In the TIPS process, a polymer solution at a high temperature is contacted with a
3 cooler liquid that is immiscible with both the solvent and the polymer. There is no mass transfer due to
4 immiscibility, so only heat transfer is considered in the transport model. In the following part of this section,
5 three specific models will be briefly introduced as representative examples for modeling the three membrane
6 formation processes described above.

7
8 Fig. 4. Schematic diagrams of the casting polymer solution and the surrounding “medium” for the different PI
9 processes: A, B, and C represent the NIPS (wet-casting), VIPS (dry-casting), and TIPS processes, respectively.
10 2.2.1 NIPS (wet-casting process): isothermal precipitation
11 For models of the isothermal wet casting process as shown in Fig. 4 (A), herein a general and rigorous
12 approach developed by Tsay and McHugh in 1990 [80] (here referred to as the TM model) is shown as a classic
13 example. The casting film and bath were treated as finite and infinite domains, respectively, and correspond to
14 the blue and yellow regions in Fig. 4. The variable l(t) denoted the (moving) position of the interface between
15 them. With assumptions of constant partial specific volume and no polymer in the bath domain, continuity
16 equations of the system could be written as Eq. (1) and Eq. (2)
17
18 = + (1),
19 = + (2),
20 for the casting film and Eq. (3)
21 = (3)
22 for the coagulation bath. Here, subscripts 1, 2, p, and b refer to the nonsolvent (1), the solvent (2), the casting
23 film side and coagulation bath side, respectively, and refers to the mutual diffusion coefficient of the solvent
24 and nonsolvent. The polymer didn’t appear explicitly in Eqs. (1) to (3), as its volume fraction was simply
25 determined by conservation. The symbols and represent the partial specific volume and volume fraction of
26 component i, respectively, and the are the appropriate phenomenological diffusion coefficients for the
27 ternary system.
28 In the TM model, the boundary conditions at the film/bath interface were handled as follows by the
29 introduction of a differential equation for the position of the interface, l(t),
30
! !
"
= = at % = &(()
! !
31 (4).
32
33 The variable l is then used as the basis for a coordinate transformation on z and thus becomes directly
34 incorporated in the final formulation of Eqs. (1) to (3). Therefore, continuity of chemical potential of the solvent
35 and nonsolvent (local equilibrium) across this interface completed the model, along with no-flux conditions at
36 the film/substrate interface and the limiting condition of bulk bath composition as % → ∞. Overall conservation
37 of mass for the solvent and nonsolvent, along with the restriction of the polymer to the film side, led to the
38 following equations
39
40 , (%, ( = 0) = ,/ for 0 ≤ % ≤ & (5),
41 , (%, ( = 0) = ,/ for 0 ≤ % ≤ & (6),

7
 1 (%, ( = 0) = / for 0 ≤ % ≤ & (7),
2 =0 at %=0 (8),
3 =0 at %=0 (9),
4 = / as %→∞ (10),
5
6 where 0 denotes the initial composition. The solutions of the model comprised the volume fractions of the three
7 species in position and time. Polymer concentration evolutions (an example is shown in Fig.5(A)) were
8 superposed on a ternary phase diagram, like that in Fig. 3, to give composition paths in the cast film Fig.5 (B).
9 Herein, the results based on calculation (Fig. 5(C)) agree with the sketches of Fig. 3 well, pointing out that the
10 CA/water/acetone system and CA/water/dioxane system had different mass paths, which meant the former
11 tended to favor sponge-type formation, while the latter favored finger-type membranes. However, for both the
12 water/DMSO and water/DMF systems, realistic solutions to the diffusion equations could not be obtained due
13 to violation of the numerical stability criteria in which D11D22 < D12D21 results. Therefore, based on the change
14 of polymer concentration profiles obtained by the model results [80] and the mass transfer paths ((an example
15 can be found in Fig.5(B))), effects of different factors including initial casting solution compositions, molecular
16 weight and bath compositions on the membrane formation and the structures were interpreted and predicted.
17 Besides, since there are no many temperature differences between the polymer solution and the coagulation, the
18 heat transfer model is not typically needed in wet casting.

19
20 Fig. 5 Evolution of CA concentration profiles in the casting film along with different immersing times (A) and
21 the corresponding mass transfer paths (B) for 10 vol % CA-acetone solution immersed in water with an initial
22 thickness of 200 μm. (C) Mass transfer paths and phase diagrams for different solvent: (a) acetone, (b) dioxane,
23 (c) DMF, (d) DMSO. Reproduced from Ref. [80] with permission.
24 2.2.2 VIPS (dry-casting): coupled heat and mass transfer
25 For models of the VIPS (dry-casting) (or the evaporation steps in wet-casting or TIPS), the ‘bath’ in contact
26 with the casting film (Fig. 3B) is normally humid air. The latent heat effects associated with evaporation of the
27 solvent and nonsolvent cause surface cooling and affect the mass transfer. Herein, the model established by
28 Krantz and Greenberg et al. in 1994 (here referred to as the KG model) is presented as the prototypical example
29 of a coupled heat and mass transfer approach [81, 82]. Mass transfer within the solution was described by
30 diffusion equations for the nonsolvent and solvent species, as in the TM model of wet casting, but the KG model
31 incorporated a simplified form of Bearman’s friction-based theory to describe the diffusion of the ternary system.
32 The governing equations were
33
5 5 5
34 =6 8 +9 (11),
7
5 5 5
35 =6 8 +9 (12),
7
36
37 where the fi and gi are phenomenological parameters that can be obtained from FH thermodynamic theory and
38 the diffusion coefficients of the ternary system, : denotes the mass fraction of component i, and ρ is the solution
39 mass density. Note that i=1, 2 denote the nonsolvent and the solvent; as in the TM model, the polymer mass

8
 1 fraction didn’t appear directly in these equations as it was simply determined by conservation. The
2 phenomenological parameters depended on temperature, which was obtained by a coupled heat transfer model
3 described below. At the solution/gas interface, located at z=l(t), the Neumann boundary conditions were
4 expressed as
5
5 5
6 8 +9 6 (1 6 : )<= + : <= = 0 (13),
5 5
7 8 +9 + : < 6 (1 6 :
= )<= =0 (14),
8
9 where <= denotes the i-th component mass flux into the gas phase. In the KG model, the <= was obtained by
10 multicomponent film theory and standard mass transfer correlations. The differential equation governing l was
11 simply based on the mass being lost by the solution, as
12
>? ">?
13 = 6 7
(15)
14
15 As in the case of the liquid/liquid interface for the TM model, the variable l was used as the basis for a coordinate
16 transformation on z and thus became directly incorporated in the final formulation of Eqs. (10) and (11). The
17 continuity of the chemical potential of the solvent and nonsolvent across this interface completed the KG model,
18 along with no-flux conditions at the interface with the substrate. For the heat transfer, the KG model employed
19 Fourier’s Law within the casting solution as
20
@ @
21 =A (16),
22
23 where α is the thermal diffusivity of the polymer solution (assumed independent of composition) and T is the
24 temperature. The boundary condition at the solution/gas interface that incorporated the evaporation of solvents,
25 convective and radiative heat transfer was expressed as
26
@ D∎ JK >? MN O >? MNO
27 =C (GH 6 G)I + (GHL 6 G) 6 6 at z = &(t) (17).
F F F F
28
29 where ℎ∎ is the heat-transfer coefficient, GH is the ambient temperature, Q is the Stefan-Boltzmann constant, R
30 is the emissivity of the polymer solution, and ΔT U is the latent heat of vaporization of the i-th component. The
31 heat transfer coefficient was obtained by a standard empirical correlation. The boundary condition at the
32 substrate/solution interface for the heat transfer model was not a simple no-flux condition, as it was for the mass
33 transfer model because the solid substrate generally conducted heat. In the KG model, a separate Fourier’s Law
34 model was used for the substrate, with the continuity of temperature and heat flux employed at the
35 solution/substrate interface and a no-flux condition employed at the lower boundary of the substrate itself.
36 In the original KG papers, the model was solved numerically using a finite element method, and numerical
37 solutions were obtained for polymer solutions of water/acetone/ CA cast on a glass support. Figure 6(A) shows
38 a comparison between the predicted and measured surface temperature profiles of the casting solution, which
39 proved that the KG model accurately predicted the time when the minimum temperature occurred, although it
40 underpredicted the minimum by 5°C. Fig.6(B) presents predicted mass transfer paths of different initial
41 compositions, which demonstrated that the mass transfer paths could predict the possibility of the phase
42 separation, the inception time and duration of the phase separation, and the type of morphology which should
43 result from a particular casting process. These model results showed the power of the model in predicting the
44 concentration and temperature of the system.

9
 1
2 Fig. 6 (A) predicted (-) and measured (●) surface temperature as a function of time for a dense polymer film
3 cast from a water/acetone/CA solution. The initial water and acetone mass fractions, film thickness, and
4 temperature are :VH/
WX = 0.02, :H[W \>W = 0.83, _/ = 246 μm and G/ = 24 C, respectively. (B) Predicted
/ o

5 mass fractions for the top (○) and bottom (●) surfaces of the cast polymer solutions whose initial conditions
6 were: (1) :VH
/
WX = 0.02, :H[W \>W = 0.83, _/ = 246 μm and G/ = 24 C; (2) :VH WX = 0.1, :H[W \>W =
/ o / /

7 0.75, _/ = 194 μm and G/ = 23.5 C; (3) :VH WX = 0.20, :H[W \>W = 0.70, _/ = 266 μm and G/ = 24 oC.
o / /

8 Consecutive points are three seconds apart and the initial solution composition is denoted by the symbol ★.
9 Reproduced from [82] with permission.
10 2.2.3 TIPS
11 The TIPS process is a little simpler than the previous wet-casting and dry-casting processes. Normally the
12 polymer/solvent (or called dilute) solution at a high temperature is immersed into a cooling bath after short
13 evaporation of the solvent to induce a polymer concentration gradient in the casting solution; however, the
14 concentration gradient normally only exists near the surface of the polymer solution. For the evaporation part,
15 the models are very similar to the dry casting models (KG model), while the main difference is the temperature
16 of TIPS is always relatively higher, which makes it a little difficult to estimate the corresponding parameters. If
17 the casting solution doesn’t include a nonsolvent, the binary diagram should be used instead of the previous
18 triangle one. For the quenching part, since only heat transfer exists in the process due to the poor interaction
19 between the bath and the polymer solution, the model would be similar to the cooling models of the dry casting
20 process, and the main difference is the boundary condition at the casting solution/bath interface because the bath
21 is normally a liquid. Compared to the wet- and dry-casting process, obvious temperature gradients can be induced
22 and should be regarded as the main reason for the resultant asymmetric cross-section structure. Fig. 7 shows
23 temperature profiles of the casting solution in the quenching stage via the TIPS process, considering a cylindrical
24 configuration to describe a hollow fiber membrane [83]. The figure illustrates that the cooling rates of different
25 parts of the solution are different, which would significantly affect the final structure.

10
 1
2 Fig. 7 The temperature profiles of the polymer solution at various quenching time: (a) 0.1 s; (b) 0.3 s; (c) 0.5 s.
3 In the figure, T and Tq represent the local temperature and the temperature of the coagulation bath. The symbol
4 ζ denotes a dimensionless position coordinate. Adopted from [83] with permission.
5 2.2. Mesoscopic phase-field (PF) approach
6 The PF approach is a powerful tool for modeling membrane formation processes due to its ability to predict
7 morphological and structural evolution during the phase separation, on time and length scales commensurate
8 with the phase domains and their evolution [84-89]. PF approach evolves continuous field variables that
9 represent volume-averaged molecular concentrations. Such field variables provide information such as the local
10 composition within a particular phase at any location, as well as the positions of interfaces between unique phase
11 domains, where concentration transiting from one value to another. This is an outstanding advantage of the PF
12 approach compared to the above continuum methods since the phase domains can be simulated and visualized.
13 In the PF approach, the transition of a field variable across an interface occurs continuously within a finite width;
14 hence, PF methods are often categorized as diffuse-interface methods. There are several broad review articles
15 and textbooks that summarize the numerous applications of the PF approach in many scientific and engineering
16 fields [90-92].
17 At present, two basic equations represent the pillars of the PF technique: the Cahn-Hilliard (CH) equation
18 to evolve conserved field variables, and the Allen-Cahn (AC) equation to evolve non-conserved field variables.
19 These two equations can be used independently or in a combined fashion, depending on the problem. Generally,
20 efforts to simulate polymeric membrane formation have utilized the CH equation alone, therefore we will limit
21 our discussion to it. The CH equation is essentially a diffusion equation for multi-component mixtures that are
22 informed by an assumed thermodynamic model for the free energy of mixing [43] as
23
lmno ( )
24 = g ∙ ijk 6 2pg q. (18).
25
26 Here, is a conserved field variable which, for polymer membrane formation research, is typically used to
27 represent the local volume fraction of polymer in a polymer-solvent solution. The mobility of the field variable
28 is represented by j, and the interfacial energy of the hetero-phase interfaces is scaled by the parameter p. The
29 thermodynamic energy of mixing 8r s can be chosen from a quantitative theoretical model (such as the FH
30 model, in which 8r s ( ) = tur ), or based on a simpler, qualitative polynomial equation.
31 While the thermodynamic components of the model are accounted for with 8r s , the kinetic component is
32 associated with the mobility parameter, j. Depending on the problem, the mobility may be assigned a single
33 uniform value, or as a function of parameters like temperature or concentration. As outlined by Zhu et al.[18]
34 diffusion in the polymer solution is highly complicated and different mechanisms have been adopted to explain
35 it, each with its advantages and limitations. The model of Duda and Vrentas, a widely used model that is based
36 on free volume theory (FVT), is accurate over a wide range of systems and concentrations but difficult to
37 implement due to the numerous model parameters [55-59][93-95]. In 2017, Bouyer et al. conducted a systematic
38 study to investigate the influence of different mobility models, including a constant, a slow (Rouse) mode [96],
39 a fast mode [97], and a model based on the FVT of Vrentas and Duda [94, 98-104], on the simulated patterns
40 and the growth law during the early and intermediate stages of phase inversion in the 2D model [105]. The
41 exponents of growth should be affected by the initial composition and the quenching conditions. However, the

11
 1 results showed that constant mobility could induce a faster pattern growth whatever the quenching conditions,
2 which demonstrated the importance to choose a realistic mobility model for the phase inversion. Recently
3 another mobility model based on the de Gennes diffusion mechanism was adopted to conduct the TIPS phase
4 separation simulations for PVDF/DPC (diphenyl carbonate) system [106].
5 By evolving the field variables in space and time (generally using a discrete numerical method such as a
6 finite-difference, a finite-element, or a spectral method), the PF approach provides a mesoscopic description of
7 multi-phase microstructures, progressing from the nucleation stage to the growth stage and into the coarsening
8 stage. Droplets or bicontinuous morphologies naturally emerge during simulations depending on the average
9 (integrated) value of the field variable, and complex interface processes such as coalescence or pinch-off events
10 occur automatically [107]. With the development of computational power and software, the numerical resolution
11 of the complex PDEs that must be involved during PF modeling work has become simple and efficient.
12 It should be noted that the PF approach is not based on a discrete molecular representation of the system.
13 Molecular interactions are considered only in a mean-field approximation, therefore making it difficult to include
14 effects such as individual chain conformations and micro-scale variations of the composition. The PF approach
15 is limited at predicting very early nucleation kinetics [104]; however, once a secondary phase has been
16 established, they are very good at evolving it during growth and/or coalescence stages.
17 2.3. Molecule/particle-based simulations
18 Since the 1970s, several molecule-based simulation methods have been adopted to investigate the
19 microporous membrane formation process, and these models are intrinsically different from the macroscopic
20 continuum model and the PF model. Typically, assumptions are made about the forces or energetic interactions
21 between segments of matter, and the trajectories of these segments over time and space are directly computed.
22 Outputs are direct representations of structures, which can be analyzed in the context of phenomena of interest,
23 such as phase domain formation and growth. Molecular dynamics (MD), dissipative particle dynamics (DPD),
24 and Monte Carlo (MC) are three molecular/particle-based approaches that have commonly been used to simulate
25 the membrane formation process via the phase inversion method. These are all methods in which a particle
26 represents a finite bit of matter at the molecular scale, but they differ in the precise scale and the method of
27 sampling particle configurations. This section presents an overview of the principles of the three types of models.
28 2.3.1 Molecular dynamics (MD)
29 In MD, the trajectories of individual atoms through time and space are computed by a stepwise numerical
30 integration of Newton’s equations of motion. Equation (23) shows the Velocity-Verlet method, one of the most
31 commonly used MD integrators.
32
33 vw (( + ∆() = vw (() + yw (()∆( + zw (()(∆()
34 yw (( + ∆() = yw (() + {zw (() + zw (( + ∆()|∆( (19)
35
36 The acceleration of atom i, represented by zw , is computed from the net force on the atom divided by its mass,
37 as zw = 8w ~} . In MD, the forces are all conservative and may arise from both bonded (e.g. bond stretching,
38 angle bending) and nonbonded (e.g. van der Waals, Coulomb) interaction. The form of the pairwise Lennard-
39 Jones model, commonly employed for van der Waals forces in MD simulations, is given in Eq. (20),
40
ƒ
Jn‚ Jn‚
41 • = 4€ •C I 6C I „ (20),
X n‚ X n‚

42
43 where • is the interaction potential between particle i and particle j, v is the magnitude of the vector from the
44 center of particle i to the center of particle j, Q and € are the size parameter and energy parameter associated
45 with the interaction between the two molecules. The parameters Q and € are the source of chemical specificity
46 in the model and they have values that depend on the chemical identity of a given pair of atoms.
47 In Eq. (20), the acceleration at time t is used to compute the velocity yw and position vw at the next point in
48 time ( + ∆(, as shown in the equations. The time step ∆( generally must be on the order of one femtosecond (fs)
49 to properly conserve energy and momentum. In practice, MD simulation runs of ~107-time steps using ~106

12
 1 atoms may be completed in a reasonable amount of time with a modern desktop computer or small cluster, thus
2 putting a practical limit of ~10 ns on the time scale and ~10 nm on the length scale of the phenomena that may
3 be studied. The MD technique is now readily available to the research community through open-source codes
4 like LAMMPS [108], NAMD [109, 110], and commercial products like Materials Studio [111] and J-OCTA
5 [112]. Fig. 8 shows an MD representative of a polymer solution/extruded solvent interface using Materials Studio.
6 The mass transfer phenomenon at the interface during the extrusion stage of the membrane formation can be
7 studied based on the MD simulations [113].

8
9 Fig. 8 An MD simulation model of interfacial behavior between PVDF, DPC, and propylene carbonate (PC).
10 Reproduced from [113] with permission
11 2.3.2 Dissipative particle dynamics (DPD)
12 DPD was originally proposed by Hoogerbrugge and Koelman [114] and later refined by Groot and Warren
13 [50]. Like Eq. (19) of MD, DPD is also based on computing the trajectories of individual particles in space via
14 stepwise numerical integration in time. Unlike MD, DPD employs particles that are substantially larger than an
15 atom but still on the molecular scale. For example, a particle might represent a set of several adjacent repeat
16 units along a polymer chain or a group of several solvent molecules. DPD can thus be considered as a kind of
17 coarse-grained MD method that provides access to larger time and length scales, which is the way we present it
18 here. However, we also note that DPD can be considered as a “quick and dirty” minimalist approach to rigorously
19 incorporating thermal fluctuations and hydrodynamics into mesoscale fluid simulations, and we refer the
20 interested reader to recent review papers [115, 116] for that viewpoint.
21 The forces in DPD are generally assumed to be pairwise, so that 8w = ∑ 8w where 8w is the force on particle
22 i due to particle j. This pairwise interactive force has three contributions: the conservative force (†w ‡ ), the
23 dissipative force (†w ), and the random force (†w ˆ ). All of these forces act within a certain cutoff radius, rc, which
24 is chosen as the length unit in the simulations. The conservative force is given by
25
A {1 6 v ⁄v[ |‹w (v < v[ )
†w ‡ = ‰
0 (v ≥ v[ )
26 (21),
27
28 where vw is the vector from the center of particle i to the center of particle j, v = Ž vw Ž, and ‹w = vw ~v . The
29 parameter A is the conservative force parameter characterizing the interaction between particle i and particle j.

13
 1 This parameter A is the source of chemical specificity in the DPD force model, just as the parameters Q and €
2 are the sources of chemical specificity in the Lennard-Jones MD force model. According to Groot and Warren,
3 the interaction parameter between the same type of bead A equals to 25•• G, while the parameter A between
4 different types of Particle i and j are chosen according to the linear relation with the interaction parameter in the
5 FH theory by the following equation [50],
6
7 A = A + 3.497‘ (’ = 3), (22),
8
9 where ‘ represents the FH interaction parameter between different species, which is a very important parameter
10 in polymer solution thermodynamics, and can be calculated based on solubility parameters [117] or experiments
11 [118]. Then, integrating the force to r will produce potential energy of interaction, analogous to that shown for
12 Lennard-Jones in Eq. (20). A comparison between the two potential models is shown in Fig. 9. This comparison
13 shows that the potential of the DPD force field is purely repulsive and much softer than the repulsive component
14 of the Lennard-Jones potential. It is the softness of this potential that permits the use of a relatively long time
15 step in DPD. The figure also shows that the particle size in DPD is much larger than that of MD. Although there
16 is no attractive component explicit in the DPD potential, the effective attraction between particles can be
17 achieved by reducing the magnitude of A .
18 The remaining two forces between the interacting particles are dissipative †w and random forces †w ˆ , and
19 they couple together to form a momentum-conserving thermostat. Due to the coarse-grained particles and the
20 soft potential, DPD can access length and time scales at or beyond the μm and μs levels, which makes it one of
21 the most promising simulation methods to investigate the phase separation and membrane pore formation. Fig.10
22 shows a morphology evolution obtained by DPD simulation on a PES membrane formation process via NIPS
23 with a zwitterionic copolymer PES-block-polycarboxybetaine methacrylate (PES-b-PCBMA) as the blending
24 additive [119]. As shown in the figure, DPD can produce the molecular spatial distribution to reveal the
25 membrane formation mechanism and membrane structure information such as porosity and hydrophilicity.
26

27
28 Fig. 9 Potential vs. distance curves for typical particle-particle force models used in MD and DPD simulations.
29 See Eqs. (20) and (21). In the figure, rc is set to 2 nm, which is a common value for DPD, and the parameters of
30 Eq. (20) are set to σij = 0.34 nm and ε/kB=120 (K), which are commonly ascribed to the argon atom. The potential
31 unit is the figure is kBT and T is chosen to be 273K.

14
 1
2 Fig. 10 Morphology evolution obtained by the DPD simulation on a PES membrane formation process via NIPS
3 with a zwitterionic copolymer PES-b-PCBMA as the blending additive. The time step are at (a) step 0; (b) step
4 20,000; (c) step 100,000; (d) step 200,000; (e) step 400,000. Tan beads represent PES; orange beads represent
5 the PES segments in PES-b-PCBMA copolymer; cyan for methyl methacrylate segments (the backbone of
6 PCBMA). Solvent beads, water beads, and carboxybetaine (the zwitterionic segment of PCBMA) segment beads
7 are omitted for the sake of clarity. Adopted from Ref. [119] with permission
8 2.3.3 Monte Carlo (MC)
9 In contrast to MD and DPD, MC simulation is built on probabilistic concepts developed from the ensemble
10 theory of statistical mechanics and there is no velocity variable used. Instead of moving through time, the system
11 samples configurations in its position variables according to the Metropolis criterion
12
13 “ = ‹”“{6∆• ⁄•G| (23)
14
15 where pij, k, and T represent the probability of making a transition from configuration i to configuration j, the
16 Boltzmann’s constant and the system temperature, and ∆• is the change in potential energy associated with
17 this transition, which is computed from particle positions using an integrated form of the conservative force
18 functions like those in Eq. (20). The fundamental unit in an MC simulation is generally an atom, molecule, or
19 coarse-grained segment, as in MD or DPD. MC simulation can be done in continuous space or on a lattice. The
20 length scale accessible to an MC simulation is ultimately limited by the ability to identify the set of particles (or
21 configuration space variables) that obey the transition probabilities prescribed by Eq. (23), so it will be similar
22 to the length scales accessible in MD and DPD. Fig. 11 presents a time-lapse graphical rendition of the polymer-
23 lean droplet growth and resulting polymer-rich solidified matrix during non-isothermal quenching of an L–L
24 TIPS system that has gone through droplet coarsening, thus pore size and pore size distribution can be calculated
25 according to the final morphologies to compare to experimental results [120].
26 One advantage of MC is that it can often explore configuration space much more efficiently than MD and
27 DPD since MC moves are not tied to small steps in time, which also brings about a disadvantage of MC that its
28 dynamic properties cannot be directly predicted. However, there are “kinetic” versions of MC in which rates are
29 assumed for each of the types of allowed fundamental move in configuration space; a relevant example is the
30 bond fluctuation model for simulating polymer chain dynamics on a lattice [53, 121].

15
 1
2 Fig.11 A time-lapse graphical rendition of the polymer-lean droplet growth and the polymer-rich matrix
3 solidification lead from isolated droplets to a co-continuous structure. (a–d) show the simulated growth of the
4 droplets with time in the periodic box, while (e) shows the polymer matrix phase at the end of the simulation.
5 Adopted from Ref. [120] with permission

6 3. Developments and applications of the dynamic modeling tools

7 From the 1970s to the present, hundreds of modeling research have been conducted to deepen the
8 understanding of membrane formation and predict the membrane properties. This section presents a detailed and
9 comprehensive review of the development history and application of the dynamic modeling tools mentioned in
10 Section 2 on the membrane research.
11 3.1. Macroscopic transport models
12 The state-of-the-art of continuum modeling involves the formulation of highly accurate yet tractable
13 transport models to predict concentration and temperature profiles in the casting solution as well as useful
14 composition paths on thermodynamic phase diagrams. The calculation results agreed fairly with the experimental
15 data, especially in the membrane and skin layer thickness as well as the precipitation period. These model results
16 have facilitated related researchers to obtain a full and comprehensive understanding of membrane preparation
17 and performance regulation. Also, due to the model construction, people have gradually realized that it is the
18 driving force of the mass transfer and heat transfer phenomena that distinguish the aforementioned various
19 membrane formation processes and also determine the resultant membrane morphologies and performance. This
20 has encouraged researchers to consider the membrane formation process based on a new perspective of a
21 tridimensional phase diagram like Fig. 2(B) [74] and then helped to explore new fabrication methods that
22 combine the two driving forces to regulate the membrane structure.
23 3.1.1 NIPS (wet-casting)
24 For models of the isothermal wet casting process that doesn’t involve heat transfer, the earliest research
25 should be traced back to the one-dimensional steady-state diffusion model developed by Cohen et al. in 1979
26 [122]. They firstly introduced a position coordinate measured in terms of the volume m of polymer per unit area
27 of the membrane between the moving interface (as shown in Fig.3(A)) and the point of observation, as shown
28 in the formula below,
29
30 }(%) = –r\U >š > WXlH[W — (%
˜
)™% ˜ (24),
31

16
 1 where — denotes the polymer volume concentration. This position coordinate was used to avoid the
2 complications arising from the motion of the membrane-bath interface and it was widely adopted in the following
3 models. Many assumptions were also employed to simplify the model construction and numerical solutions were
4 gained to depict a composition path on the ternary phase diagram. Although the first model had some
5 inaccuracies due to too many assumptions employed, which was pointed out by Wijmans and McHugh later [75,
6 123], it built a foundation for subsequent studies. Then from 1985 to 1988 in view of a full and rigorous analysis
7 on the thermodynamics of the ternary phase behavior, McHugh et al. adopted a pseudo-binary model to enable
8 decoupling of a specific mixing rule. Thus, the analytical solution was obtained and the effects of different
9 boundary conditions on the polymer volume fraction profiles were studied [75, 124, 125]. At the same period,
10 Reuvers et al. made a considerable advance to the aforementioned Cohen’s model to cover the mass transfer of
11 the casting solution and the coagulation bath [126, 127]. They firstly proposed a diffusive layer existing at the
12 moving interface and then provided a comprehensive analysis that two types of membranes can be formed
13 according to the composition path change in an initial short time. This model was considered to be a standard
14 for many successive studies due to the detailed and appropriate description of the polymeric system. However,
15 the model’s computation was only valid for very short transfer times due to the use of an infinite thick casting
16 solution. In 1992, Radovanovic and coworkers modified one equation of the coagulation mass transfer of
17 Reuvers’s model [127] and developed a new way to obtain the model parameters based on experimental flux
18 data. The predictability of the computational results was qualitatively verified by the membrane structure
19 evolution of the PSf-N,N-dimethylacetamide-2-propanol system [128, 129].
20 In 1990, Tsay and McHugh developed the rigorous TM model as mentioned in Section 2.2.1. They
21 abandoned Eq. (24) and the model laid another foundation for the following research. Based on the model results,
22 the effects of different factors including initial casting solution compositions, molecular weight, and bath
23 compositions on the membrane structures were interpreted and predicted. Subsequently in 1994, by employing
24 solvent-polymer binary diffusion coefficients and a moving interface velocity, Cheng et al. improved Reuvers’s
25 model and the TM model to extend the computational time until the system reached the spinodal line [130, 131].
26 The model results could help to predict various membrane structures including skin formation, uniform, and
27 graded cellular porous, which agreed with experimental results, although it failed in predicting the skin rupture
28 and the finger morphology. In 2001, considering the difficulty of obtaining the model parameters in the previous
29 models, Fernandes et al. proposed a simple model based on the Fick diffusion equations, which could be regarded
30 as a rough and simple way to follow for further practical and industrial applications [132]. Also, in 2001, Karode
31 et al. developed an improved algorithm to numerically solve the governing equations of Reuvers’s model to
32 avoid any intervention in terms of initial guesses for the interfacial composition [133]. Meanwhile, Kim and
33 coworkers considered the equilibrium boundary condition such as Eq.(4) could not describe the mass transfer
34 process of the systems that possessed a very rapid demixing rate, and then developed a new boundary condition
35 for polyurethane (PU)/dimethyl formamide (DMF)/water system [134]. The previous many models restricted
36 their calculation to the spinodal curve and the phase separation onset, whereas Kim’s model extended the
37 simulation time range until after the phase separation [134]. Almost 10 years later in 2010, based on the TM
38 model [80], Krantz and Hwang et al. developed a new model to incorporate “densification-induced convection”,
39 which means densification contribution to the mass-transfer flux, by adopting a new rigorous method for
40 obtaining an explicit equation for the mass-average velocity that arises owing to densification [135]. If the casting
41 solution was relatively thick, the calculation could provide a reasonable explanation for the macro-void
42 formation of the sub-layer. In 2017, Khansary and coworkers tried to build a rigorous thermo-kinetic model for
43 the isothermal mass transfer process by incorporating a compressible regular solution free energy model [36].
44 This model employed a new insight (as shown in Fig. 12) that the casting solution should be regarded as a
45 polymer-rich and a polymer-lean phase, for which different governing equations and the thermodynamic
46 equation would be solved after the phase separation. Besides, the enforced equilibrium criteria at the moving
47 interface were improved by a kinetic approach so that the composition of the interface was no longer fixed to
48 the vitrification point. This model should be regarded as a big step that the homogenous assumption for the
49 casting solution was replaced by multiple phases consideration.
50

17
 1
2 Fig. 12 Handling and calculation scheme for inner phase-separated points within the casting film wherever and
3 whenever it emerges. Reproduced from Ref. [36] with permission
4 3.1.2 VIPS (dry-casting)
5 For models of the dry-casting as well as the evaporation steps of TIPS, the earliest model should be traced
6 back to the one-dimensional unsteady-state diffusion model upon evaporation of volatile solvent acetone from
7 CA/acetone solution proposed by Anderson and Ullman in 1973 [136]. This model calculated solvent profiles
8 and polymer contraction within the thin interface layer based on many assumptions including a semi-infinite
9 casting film thickness, a specified constant surface concentration, negligible film shrinkage, and isothermal mass
10 transfer. Subsequently, in 1981, Castellari and Ottani relaxed the assumptions of semi-infinite film thickness and
11 specified constant surface concentration, and then proposed a modified diffusion model allowing for uniform
12 film shrinkage due to solvent loss [137]. In 1986 Krantz and coworkers developed a fully predictive model that
13 incorporated a semi-empirical correlation for the binary-diffusion coefficient and convective contribution to the
14 mass transfer flux because of local film thinning arising both from direct solvent loss and the excess volume of
15 mixing contraction effects [138]. Using dimensionless analysis, they claimed that the gas phase mass transfer
16 coefficient should be regarded as the most influential process parameter for skinning formation. In 1991,
17 following Krantz’s model, Tsay and McHugh adopted a more predictive correlation of the binary mutual
18 diffusion coefficient based on FH theory and the free-volume theory [139]. They incorporated the surface solvent
19 evaporation cooling effect that significantly influenced the results whereas was excluded by the aforementioned
20 models.
21 In 1994, Krantz and Greenberg et al. established the KG model mentioned in Section 2.2.2 [81, 82].
22 Focusing on CA/acetone/water system, this model was capable of predicting the concentration and temperature
23 profiles, instantaneous casting solution thickness, densified “skin” thickness, and composition path, although
24 ignoring a convection term. As a subsequent work, in 1995 Krantz et al. improved the KG model by adding a
25 convective term to the governing equations to make the calculation results agree more with the experimental
26 data [140]. In 1997, Matsuyama et al. studied the membrane formation and structure development aiming at
27 CA/acetone/non-volatile nonsolvent systems both experimentally and theoretically [141]. The evaporation
28 model was similar to the model proposed by Tsay and McHugh [139]. Based on a comparison between the
29 experiments and the model calculations, they concluded that the asymmetric structures obtained from these
30 systems were found to be attributable to the difference in the increasing rates of the polymer fractions at the
31 gas/liquid surface. However, the mass paths in the research were nearly straight due to the neglect of the
32 nonsolvent diffusion. Subsequently in 1999, focusing on the PVDF/DMF/water system, Matsuyama et al.
33 modified the previous evaporation model by adding a nonsolvent diffusion flux and analyzed the phase
34 separation induced by the penetration of water numerically and experimentally [142, 143]. The membrane
35 morphologies evolution could be explained based on the calculation results although with ignoring the
36 evaporation cooling effect.
37 In 2004 and 2005, Altinkaya and Ozbas investigated the consistency between the model computation and
38 experimental results [144, 145]. Different from the KG model, herein the temperature through the casting

18
 1 solution and the substrate was uniform and the heat transfer was approximated by a lumped parameter approach,
2 which was reasonable because the resistance to heat transfer in the gas phase was much greater than that in the
3 polymer solution or the substrate layer, which was widely applied in other modeling work [146, 147]. In 2006,
4 Krantz and Hwang et al. improved the KG model (Eqs. (11 to 17)) by incorporating a convective transport term
5 owing to density changes [148]. The model was similar to the TM model [135] since they came from the same
6 research group, in which a modified equation-of-state for the mass density was employed. In 2006, Yip and
7 McHugh incorporated a diffusion model that was derived from the friction-based diffusion model proposed by
8 Alsoy and Duda [146]. This model allowed the diffusion of nonsolvent into the film from the atmosphere, which
9 was the model’s strength. In 2010, Krantz and Greenberg et al. conducted a full sensitivity analysis to assess the
10 relative importance of the various casting properties and parameters, including the FH interaction parameters,
11 the friction coefficients, and thermal properties on the basis of their previous model [81, 82], hence a
12 phenomenological model was proposed to provide a set of design heuristics for the dry-casting process [81, 149].
13 In 2010, Bouyer et al. conducted modeling and experimental research to analyze membrane morphologies
14 of poly(ether imide) (PEI)/1-methyl-2-pyrrolidone (NMP)/water system [150]. This work didn't only test the
15 effectiveness of several various diffusion formalisms of the polymer solution mentioned in Section 2.2 but also
16 linked the self- and mutual-diffusion coefficients. It should be noted that in 2013 Bouyer and Pochat-Bohatier
17 presented a comparison between experimental kinetic data and the corresponding simulation results during the
18 VIPS process for the ternary PEI/NMP/water system [151]. This is the first work to directly validate modeling
19 composition profiles by in situ near-infrared spectroscopy measurement (NIRS) data. The schematic of the
20 process performed in the deep cell dedicated to in situ Near-infrared spectroscopy (NIRS) measurement can be
21 found in Fig. 13 (A) while the comparison between the model predictions and experimental data is shown in Fig.
22 13(B), which demonstrates that the the agreement between the numerical results and the experimental data was
23 fairly acceptable. Besides, the numerical model was also firstly shown to overestimate the nonsolvent diffusion
24 rate for higher polymer concentration, however this could be improved by the addition of a modified expression
25 of the average hole-free volume.

26
27 Fig. 13 (A): Schematic of the VIPS process performed in the deep cell dedicated to in situ Near-infrared
28 spectroscopy (NIRS) measurement; (B) Comparison between numerical predictions and experimental data for
29 the initial PEI concentration was 12 wt.%, the temperature was 40oC, and the relative humidity was 75%.
30 Reproduced from Ref. [151] with permission
31
32 3.1.3 TIPS
33 From 1998 to 2000, Lloyd et al. conducted a series of theoretical, modeling, and experimental work to
34 investigate the effects of different membrane formation conditions on the membrane structure [152-155]. The
35 evaporation model was similar to the TM model but only incorporating an extra convective term, and the
36 quenching model was similar to the model shown from Eqs. (16) and (17). The component volume concentration
37 and temperature profile evolution were calculated and proven to agree with the membrane microstructure and
38 membrane thickness obtained from the corresponding experimental observation. At the same time, Barton and

19
 1 McHugh calculated the evolution of pore size gradients in the casting solution via a combination of the simple
2 heat transport (Eq. (16)) and kinetics of the droplet growth rate, which was the first TIPS modeling work to
3 predict the membrane pore morphologies [156]. In 2006, Li and Krantz et al. incorporated Lloyd’s model [152]
4 with the polymer crystallization Avrami theory to predict spherulite sizes prepared by TIPS solid-liquid phase
5 separation process as shown in Path (C) of Fig.1 [157, 158]. In 2016, Tang and Wang et al. established a modified
6 Maxwell-Stefan model to describe the mass and heat transfer in the air gap and coagulation bath stages of the
7 TIPS process, taking into account a cylindrical configuration [83]. The model form was similar to Lloyd’s
8 models [152] and an example of the results can be found in Fig.7. In 2017, based on the FH theory, a specific
9 diffusion formalism for dilute systems, and external mass transfer in free convection, a numerical model were
10 developed by Bouyer et al. to investigate optimized operating conditions for the PVA/water system, because of
11 an interest to use water as a solvent instead of traditional organic solvents [36].
12 During the development of the mass transport models, much effort has been spent to obtain mutual diffusion
13 coefficients ( ) or friction parameters (› ) as well as the thermal diffusivity A of the ternary systems to solve
14 the diffusion models. In many models reviewed in Section 3.1 [71], the mass transfer between the components
15 was always expressed by relations derived from the Bearman statistical mechanical theory [159], which showed
16 that the chemical potential gradients could be written as: œ• ⁄œ% = ∑—Ÿ ž › (y 6 y ) where œ• ⁄œ% was a
17 one-dimensional chemical potential gradient of component i, while ž , › , and y were the molar concentration
18 of component j, the friction coefficient between molecule i and j, and the velocity of component i. Nevertheless,
19 it was hard to get › for a specific system. Therefore, different assumptions were applied to estimate the friction
20 parameters. For miscible fluids, the friction factors could be derived from the mutual-diffusion coefficients in
21 the limit of quasi-binary diffusion [81, 160]. Considering the solvent/polymer systems, Vrentas and Duda
22 proposed to estimate self-diffusion coefficients based on the free-volume theory, which included thirteen
23 independent parameters and molecular transport description by jumps of small molecules into voids, also called
24 hole-free volume [161]. Several formalisms were proposed to express the mutual-diffusion coefficients in terms
25 of self-diffusion coefficients and gradients of chemical potential [146, 162, 163], which aimed to simplify the
26 Bearman equation that involved individual frictional forces between each molecule of the multicomponent
27 system. These mechanisms were also applied to the Fickian diffusion [141, 144]. Recently Bouyer et al. tested
28 the effectiveness of several various formalism theories including Alsoy and Duda, Dabral, Zielinski and Hanley,
29 as well as Price and Romdhane’s, based on a comparison between the modeling results and experimental data
30 [150]. The results showed that the Dabral formalism gave the best model predictions with regards to the
31 membrane morphologies and the cross-diffusion term Dij could not be neglected.
32 In summary, during the period from the 1970s to the present, the continuum models have been gradually
33 improved by the elimination of erroneous assumptions and the application of more physical considerations of
34 the transfer problems. The latest model can allow film shrinkage, evaporative cooling, coupled transfer and
35 incorporates reliable diffusion theory as well as complex boundary conditions especially at the polymer
36 solution/air interface. It has also been acknowledged that the mass paths could be regarded as an accurate tool
37 to reason how mass transfer affecting membrane morphology. However, due to the difficulty in the measurement
38 of composition and temperature change in the casting solution during the membrane formation, few experimental
39 methods were obtained to directly verify the theoretical calculations. Also, as reviewed in this section, the
40 transport model is weak in dealing with spinodal decomposition, and membrane pore morphology is still hard to
41 predict. This issue is unlikely to be resolved by just reformulating continuum-level models, and then people tend
42 to find new methods such as PF approaches or molecule-based simulations to study the membrane formation
43 process further, which will be discussed in the next sections.
44 3.2. Mesoscopic PF approach
45 Below, we review the literature associated with PF models applied specifically to the membrane formation
46 process. We acknowledge there are a variety of PF studies on the somewhat related topic of phase separation in
47 binary polymer blends [164-167], however, we wish to limit our review to polymer/solvent or
48 polymer/solvent/nonsolvent models.
49 The earliest application of the CH equation to specifically simulate the polymer membrane formation
50 process via TIPS can be attributed to Caneba and Soong in 1985 [168, 169]. They conducted one-dimensional
51 simulations of the TIPS process in a polymer-solvent system at various locations relative to a cooling surface,

20
 1 using the FH and FVT models for the thermodynamic and kinetic descriptions, respectively. Their results
2 estimated pore sizes as a function of membrane depth away from the cooling surface and demonstrated the
3 versatility of this approach. Soon after, two-dimensional simulations followed that and investigated the
4 temperature-induced spinodal decomposition of polymer-solvent systems quenched into the unstable region of
5 the phase diagram, including works by Chan and Rey and Barton et al. [170, 171]. These studies assumed
6 isotropic quenches and focused on the growth and coarsening rates of the polymer-rich and polymer-poor phases.
7 Lee et al. then applied PF techniques to understand the effect of temperature gradients (i.e. anisotropic quenches)
8 on the phase morphology, firstly with one-dimensional [89] and then two-dimensional simulations [107].
9 Anisotropic domain morphologies were analyzed, and were found to depend on the rate of quenching in both
10 the lower and higher temperature regions, as shown in Fig.14 (a, b). Kukadiya et al. [172] extended this study
11 with additional two-dimensional simulations of non-uniform temperature fields, and Chan conducted a similar
12 investigation corresponding to concentration gradients [173]. Recently, Mino et al. conducted three-dimensional
13 simulations of the TIPS process, including the effects of a polymer concentration gradient that leads to an
14 anisotropic structure, as shown in Fig.14 (c to f) [174]. In 2019, Millett et al. conducted large-scale three-
15 dimensional computer simulations to investigate the TIPS process for the PVDF/DPC system and predict the
16 complex networks of porosity. Although the domain size is an order of magnitude smaller than typical polymer
17 sheet membranes, this work demonstrated the simulation tool is fast becoming a valuable tool to predict pore
18 structure throughout an entire membrane cross-section[106].
19

20
21 Fig. 14 In (a-b), late-stage progress of a 2D in phase field simulation of TIPS for a uniform quench case (a) and
22 temperature gradient case (b) during the phase separation. Reproduced with permission from Ref. [107]. Images
23 (c-f) are a 3D anisotropic morphology evolution during the phase separation under a linear concentration profile
24 from = 0.25 at the bottom and = 0.35 at the top. Reproduced with permission from Ref. [174]
25 The earliest attempt to simulate the NIPS process was made by Barton and McHugh, who chose two
26 common membrane-forming systems: CA/acetone/water and PES/dimethyl sulfoxide/water [52]. To simulate
27 the NIPS process, a ternary PF model was required. Barton and McHugh varied the relative quantities of non-
28 solvent and polymer to understand their influence on the structure growth rates [52]. Saxena and Caneba carried
29 out a similar study, using poly(methacrylic acid)/methacrylic acid/water [175]. Zhou and Powell extended this
30 concept to two-dimensional and three-dimensional simulations, in which the systems were initialized to have
31 two separate domains associated with the non-solvent and the polymer solution, with simulations capturing the
32 exchange of solvent and non-solvent elements across the interface, which lead to anisotropic phase separation
33 [176]. Recently based on a hypothetical fluid mixture model with symmetric miscibility gap, Hopp-Hirschler
34 predicted different resultant pore morphologies, including the sponge pores, finger pores, and lamella structure,
35 depending on the initially reduced mass fraction and the velocity of the precipitation front [177]. Immersion
36 precipitation has also been investigated using other somewhat similar mesoscale methods, including the lattice-
37 Boltzmann method applied by Akthakul et al. [178] and a multi-fluid model by the Fredrickson group[179].

21
 1 Both approaches capture the convective fluid flow in addition to diffusion, as opposed to the purely diffusive
2 kinetics of the standard CH equation.
3 The PF approach has also been used to investigate the liquid-to-solid phase transformations that accompany
4 the membrane formation process, which is related to the path (C) of Fig.1. Because liquid-to-solid phase
5 transformations are non-conserved in the fraction of liquid (or solid) phase fraction, these models have utilized
6 the AC equation. Xu et al. investigated polymer crystallization and the variety of dendritic growth patterns that
7 emerge at different crystallization temperatures (the model was capable of simulating a single crystal orientation)
8 [180]. Later, Asle Zaeem et al. expanded the model to include multiple crystal nuclei with multiple crystal
9 orientations that led to polycrystalline packings of spherulites [181]. Recently, Wang et al. used a PF method to
10 study flow-induced crystallization of polymers with multiple crystal structures forming simultaneously [182,
11 183].
12 With the development of the computational power and resources, the numerical resolution of the complex
13 PDEs that has to be involved during the macroscopic transport and PF modeling work has become simple and
14 efficient due to the help of some powerful finite element software such as COMSOL [105] or Matlab [172].
15 Compared to the continuum models on the macroscale, PF methods are particularly suited to produce simulations
16 of membrane formation processes at the length scale of several micrometers, namely the length scale of the pore
17 structures. Published studies have demonstrated the ability to simulate isotropic and anisotropic temperature
18 quenches associated with the TIPS process, as well as immersion precipitation processes associated with the
19 NIPS process. Although the molecular-scale details of the polymer distributions are implicit, the mean-field
20 approximations of polymer distributions seem to be adequate to predict a variety of structural features that are
21 qualitatively similar to experimental observations. PF models are dependent on the accuracy of the mean-field
22 thermodynamic and kinetic theories they employ, and it can be a challenge to obtain specific material-dependent
23 parameters associated with a specific polymer/solvent/non-solvent system, as discussed in Section 2.2.
24 Therefore, more effort can be contributed to make connections between the PF mobility models and the
25 microscale molecular-based simulation methods, which also helps to improve the accuracy of the simulation
26 results.
27 3.3. Molecule/particle-based simulations
28 Since the 1970s, several molecule/particle-based simulation methods like MD, DPD, and MC have been
29 adopted to investigate the microporous membrane formation process. Compared with the continuum and PF
30 approach, molecular-scale modeling has an inherent advantage that the only assumption is at the level of the
31 particle force field. Its main disadvantage is that so far due to the limitation of the computational power, the
32 simulated process still can’t be followed on appropriate length and time scales, especially time. However, with
33 the rapid development of computational hardware and advanced algorithms, the particle-based approaches will
34 provide a solid validation to theoretical characterization and a natural complement to experimental observation.
35 3.3.1 MD
36 Considering that a complete phase inversion membrane formation process normally involves time scales
37 on the order of seconds and length scales on the order of microns, it’s clear that an MD simulation will not be
38 able to follow a complete process. However, MD can be used to study the earliest stages of phase separation
39 during membrane formation. The earliest work in this area can be traced back to the research accomplished by
40 Mruzik and coworkers in 1978 [55]. The MD simulations were conducted to understand the spinodal
41 decomposition mechanism for a single Lennard-Jones fluid system undergoing a quenching process into the
42 vapor-liquid coexistence region. Structure factor, iso-surface morphology, and average density were used to
43 characterize the phase separation process; the simulations showed an initial exponential growth of density
44 fluctuations that agreed with a generalized diffusion theory of the structure factor dynamics. Since then, in 1995
45 Leptoukh and coworkers [184] used MD to study the phase behavior of two-dimensional fluid mixtures after an
46 asymmetric quench into the spinodal region of the phase diagram. For the binary fluid system, different relations
47 between the growth exponents and the reduced temperature were observed for different volume fractions, and
48 the reshaping of droplets played an important role in inducing long-range droplet flow. In 1999 Liu and
49 coworkers [185] carried out Brownian dynamics simulations (similar to MD, but with an implicit solvent) to
50 study the spinodal decomposition in a quenched polymer solution, also with an LJ potential, in two and three
51 dimensions. The simulation results demonstrated the existence of network pattern formation at early to

22
 1 intermediate time regimes for sufficiently dense polymer solutions, and especially for longer polymer chains, as
2 shown in the case of t = 400 in Fig. 15 (A). This network structure changes to droplets at later times to minimize
3 the interfacial energy, as seen in later times of Fig. 15 (A). Instead of simulating the whole phase separation,
4 recently people tend to adopt MD simulation to study the affinity between the polymer and some additives. In
5 2018, Liu et al. applied MD to study the compatibility change of PU/polyvinyl chloride (PVC) system with a
6 variation of PU/PVC mass ratio to optimize the PVC amount for the PU membrane preparation via NIPS [186].
7 The simulation results showed that the compatibility between PU and PVC is the best with a mass ratio of 90
8 /10. Then, Yonezu’s group carried out MD simulations to clarify the fouling mechanism [112]. As shown in Fig.
9 15 (B), the simulation results illustrated that hydration layers with a higher water density was formed by the
10 polyvinyl pyrrolidone (PVP) addition and the repulsion force for benzene from the water molecules was larger
11 for the PSf + PVP membrane as compared with the PSf membrane. This year, Matsuyama’s group employed a
12 combination of experimental results and MD simulations to investigate the mass transfer between the extruded
13 solvent propylene carbonate (PC) and PVDF/DPC solution by using a triple-orifice spinneret in the TIPS process
14 [113], as shown in Fig. 8. The MD calculation results elucidated the mutual penetration of PC and the PVDF
15 solution from the contacting interface occurred and increased as increasing the contacting time, which
16 subsequently affected the phase separation mechanism of the PVDF solution near the contact interface and the
17 resultant membrane structure.
18 With the development of the computational power and application of the commercial and open-source MD
19 software, now people no longer use the MD method to simulate the phase separation due to its limitation in
20 Spatio-temporal scales but tend to use it to reveal the interaction between the system components and calculate
21 some key thermodynamic parameters such as solubility parameters, to assist and verify the experimental results
22 to optimize the membrane preparation parameters. In the future due to the advantage in precision, MD will have
23 an increasingly wide utilization in the fields of membrane preparation optimization and modeling parameter
24 exploration, especially for DPD and PF.
25

26
27 Fig. 15 In (A), snapshots of evolutions in 2D MD simulation of TIPS for a 30% polymer solution with chain
28 length N=100 for times t = 400, 1200, 2400, and 4000, respectively. Reproduced with permission from Ref.
29 [185]; In (B), MD simulation models of water molecules on (a) PSf membrane and (b) PSf with PVP additive.
30 Reproduced with permission from Ref. [112]
31 3.3.2 DPD
32 To our knowledge, Lv and coworkers in 2008 [54] were the first to employ DPD to simulate the immersion
33 precipitation (wet casting) process, using a two-dimensional model. From 2009 to 2020, Tang et al. have
34 established a DPD simulation methodology and published a series of papers on the investigation of the membrane
35 formation process via TIPS [187-192]. They focused on specific polymer-solvent systems to provide some
36 detailed information that could be connected to experimental observation. The simulation systems have been
37 extended into millions of particles. The box length and the simulation time can reach 0.2 μm and 2 μs [191, 192].
38 Different membrane preparation conditions, including quenching temperature, initial polymer concentration,
39 additives, have all been covered to get a full understanding of the membrane formation process via the PI process.

23
 1 Fig. 16 (A) and (B) show the morphology evolutions of the phase separation via the NIPS and TIPS processes,
2 respectively. The latest research has found that pore size and phase separation mechanism can be directly
3 disclosed based on the simulation results, and multiscale simulation that is closer to a real process, should be
4 paid more attention to in the future. In 2016 and 2019, Li et al. have applied DPD to simulate the membrane
5 morphologies and gas permeability of poly(4-methyl-1-pentene) (PMP) during TIPS, which to our knowledge
6 is the first such work to simulate the fabrication and design of gas-diffusive membranes [193, 194]. Fig. 16 (C)
7 presents the 3-D morphology of the diffusion of fluid (hypothetical) particles through the constructed membrane
8 structures with 30 wt% PMP. The results demonstrated that the more connective open-pore structure in Fig. 16
9 (C, a) provided less transport resistance, while the weakened pore connectivity in Fig. 16 (C, b) and Fig. 16 (C,c)
10 provided greater transport resistance because of the additional blockage from the polymer domain along the
11 pathway. In 2019, Zhou and coworkers adopted DPD simulation to study the effect of the additive (a zwitterionic
12 copolymer) PES-B-PCBMA on the PES membrane formation via NIPS, as shown in Fig. 10 [119]. The
13 simulation results proved that the hydrophilic PCBMA segments enriched on the membrane surface by surface
14 segregation and exhibited pH-responsive behavior, which was attributed to the deprotonation of the carboxylic
15 acid group.
16 Although DPD has been acknowledged by researchers to be a promising technique to understand the
17 membrane formation mechanism, there are still some challenges: firstly, the spatiotemporal scales have been
18 enormously expanded into several μms and μss in 3 dimensions, nevertheless, it is still very hard to simulate a
19 real membrane pore since the pore size of the microfiltration membranes is around 0.1 μm and the time for the
20 formation of membrane pores may exceed several seconds for a typical TIPS process and an instantaneous
21 demixing. Secondly, the membrane formation is a consequence of the phase separation and polymer
22 solidification, so how to determine the termination of the phase separation and get the resultant polymer matrix
23 should be figured out in the future; furthermore, the coarse-grain algorithm needs further elaboration to construct
24 a more clear and accurate mapping between the DPD particles and the microscale atoms or molecules.
25

26
27 Fig. 16 (A) DPD simulation morphology evolution of the PES/NMP/water systems with the initial polymer
28 volume fractions (12 v/v) at three different time steps via the NIPS process. Reproduced from Ref. [192] with
29 permission; (B) DPD simulation morphology evolution for the polymer solution at different quenching times.
30 Reproduced from Ref. [191] with permission; (C) Morphology of the diffusion of hypothetical liquid beads
31 through a PMP-diluent system after a 2500τ phase separation evolution ((a) PMP-dioctyl phthalate (DOP)

24
 1 system; (b) PMP- diphenyl ether (DPE) system; (c) PMP- dibutyl phthalate (DBP) system). Reproduced from
2 Ref. [193] with permission.
3 3.3.3 MC
4 The application of MC to the problem of membrane formation is a little different from MD and DPD. From
5 1995 to 1997, the original two MC simulation studies on the membrane formation via the phase inversion method
6 including NIPS and TIPS were conducted by Termonia et al. in two and three dimensions [195, 196]. Based on
7 a lattice algorithm, the systems employed a cubic box of (90nm)3 and a molecular weight of 140,000, so
8 membrane formation mechanisms at the mesoscale level were captured. In 2012, Zhang and his coworker
9 employed a new MC diffusion model to simulate the phase separation during the wet-spinning process. The
10 results could predict composition paths that were able to impose on the phase diagram so that different fiber
11 structures ranging from dust-like, finger-like to sponge-like morphologies are obtained. Nevertheless, the whole
12 separation process was not captured and hypothetical parameters were only included in the work [197]. From
13 2011 to 2013, Jiang et al. employed a polymer bond fluctuation lattice model cooperating with MC to simulate
14 polymeric ultrafiltration membrane formation via the NIPS process in two dimensions [198, 199]. Based on the
15 new modified MC simulation model, the addition of amphiphilic block copolymers to the anti-fouling polymeric
16 ultrafiltration membrane formation was studied. These results firstly displayed a clear two-dimensional phase
17 separation process of the NIPS method (as shown in Fig. 17), then proved that the asymmetrical structure induced
18 by the fast exchange between the solvent and the nonsolvent could be adjusted by controlling the mechanism of
19 the phase separation, such as changing the interaction between the solvent and the nonsolvent to adjust the mutual
20 diffusion or adding some amphiphilic block copolymers to change the thermodynamics of the systems. Moreover,
21 in a different type of approach, an off-lattice MC simulation based on mesoscopic particles was employed to
22 simulate the growth of the phase domains and predict membrane pore sizes for some real polymer-diluent
23 systems via the TIPS process and the results were verified by experiments [120, 200].
24 These research works have proved that MC is an effective and efficient way to get an understanding of
25 portions of the membrane formation process via the phase inversion method. However, quantitative predictions
26 of the process dynamics and the mapping between the lattice MC model and the real polymer-solution systems
27 still represent challenges in the future.

28
29 Fig. 17 Morphologies of phase structure at MC step = 2.3×105 using the different initial concentrations at the
30 temperature ω= 0.624: (a) φ = 0.2; (b) φ = 0.3; (c) φ = 0.4; (d) φ = 0.5. Here the temperature ω = 1/kT and T is
31 the temperature while k is the Boltzmann constant. Note that in the figure, different colors represent various
32 species: polymer (red), solvent(cyan), and nonsolvent(yellow). Reproduced with permission from Ref. [198].
33
34 As mention at the beginning of Section 3.3, molecule/particle-based simulation methods including MD,
35 DPD, and MC have an inherent advantage that the only assumption is at the level of the particle force field,

25
 1 which is inevitably and significantly influences the accuracy of the simulation results. For MD and MC, there
2 are kinds of molecular or atomistic force fields that can be applied depending on the computational complexity
3 and characteristics of the target systems [201]. Abundant knowledge and theory have been accumulated in the
4 field of liquid simulations by the MD and MC methods, and interested readers may refer to review articles and
5 textbooks that summarize the principles and theory of the two approaches [202]. As introduced in Section 2.3.2,
6 the repulsion parameter A of the conservative force that can be obtained based on the FH interaction parameter
7 ‘ is the only source of chemical specificity in the DPD force model, and also the only connection between the
8 simulated and real system. Therefore, for the classic DPD method proposed by Warren [50], the interaction
9 parameters of different species ‘ related to the temperature should be the only important factor to strongly
10 affect the accuracy of the DPD simulation results. Whereas, additional force fields have also been proposed to
11 enhance the connection between the DPD simulated system and the reality, with a sacrifice of computational
12 speed [203-205]. In general, if one wants to promote the predictability of the molecule/particle-based simulation
13 results, a more precise force field and corresponding parameters should be firstly considered.
14 The simulated phase separation in Section 3.3 is mainly related to path (A) of Fig. 1. So far it is still hard
15 to simulate the process of Path B by the simulation tools described here, whereas there are some polymer
16 crystallization simulation studies in connection with the path (C). MD simulation has been employed to simulate
17 a single polymer chain crystallization from a dilute solution during a cooling process [206, 207]. The chain
18 folding and conformation transformation were observed and characterized during the quenching and cooling
19 processes. Meanwhile, due to the merit of computational efficiency, MC has also been used to perform
20 simulations on the crystallization of entangled polymer chains in dilute solution [208, 209]. The results showed
21 that temperature and the stiffness of polymer chains had significant effects on the crystallization rate and the
22 simulated morphology. Although these simulation studies didn’t aim to describe a membrane formation process,
23 the results can still be considered as relevant for the membrane preparation related researchers, since the path
24 (C) is just polymer crystallization from a semi-dilute polymer solution, and the polymer crystallization process
25 must be put into consideration for the path (A)and the path (B).
26 3.4. Challenges and outlook
27 Membrane technology is now well-developed and widely used in numerous industrial processes, and
28 notably, it can greatly benefit from improved intrinsic separation properties of membranes such as selectivity
29 and permeability. As described above, with a half century’s efforts each modeling and simulation method has
30 gradually matured and can provide some important information about the membrane formation mechanism. It
31 can be known from the summary above, much experimental and modeling research has been conducted for the
32 classic membrane systems such as CA, PES, PVDF, and the common solvents, etc. So, the modeling work may
33 not provide much improvement to those classic polymer systems. In the future, in response to the international
34 call of “sustainable development and environmentally friendly”, there is a real drive towards making membranes
35 more sustainable by looking at alternatives to the current unsustainable solvents such as NMP and DMAc, which
36 means new green solvents and novel membrane formation methods to prepare new membrane materials. In
37 addition, these new membrane materials and enhanced manufacturing methods should have the potential to push
38 the trade-off curve to higher permeability (or permeance) and selectivity values and possibly provide close to
39 defect-free membranes [9], which is also significant for many membrane-manufacturers at industrial scale along
40 the world using the PI processes. Modeling and simulation cannot only help people to get a deeper understanding
41 of the mechanism of membrane formation but also allow people to study the effects of the solvent molecular
42 structure and intermolecular interaction on phase separation and membrane formation. Therefore, it can be
43 inferred that more advanced and precious modeling and simulations could guide the experimental work into a
44 fast lane, including the solvent and additive selection, process control and prediction, et al., especially with the
45 help of the much faster computers that have become common recently.
46
47 Table 1 Key issues on different scales for the membrane formation via PI
Key issues in the membrane formation Transfer model Phase field model Molecular simulation
Molecular information Not involved Not involved Involved
Phase domain evolution Not involved Involved Involved
Concentration/temperature field Quantitatively Quantitatively Semi-quantitatively

26
 Not directly
Characteristics of membrane structure Directly involved Directly involved
involved
Semi-quantitatively
Membrane performance Not involved Not involved
involved
CA, PSf, PES, PMMA, PVDF,
PMMA CA, PES, PSf, PVDF, PP, PMP
Real polymer system involved or
(polymethyl polystyrene (PS), (poly (4-methyl-1-
addressed
methacrylate), polypropylene pentene))
PU, PEI (PP)
1
2 Whereas, at present, we think that there are still two critical issues that need to be fixed in the next few
3 years. The first would be developing more powerful and rigorous models. Table 1 lists key issues of greatest
4 interest in the membrane formation process, including molecular interaction, phase domain evolution,
5 concentration or temperature field, Characteristics of membrane structure, and membrane performance as well
6 as the polymers that have been addressed. As listed in Table 1, the three kinds of dynamic models have their
7 superiorities and limitations in dealing with these issues. Although it seems that molecular simulation is
8 preferable due to almost involving all the issues, it should be noted that at present, the spatiotemporal scale of
9 state-of-the-art molecular simulation still cannot meet the real membrane formation scale, which makes the
10 results are not accurate enough and only can be regarded as semi-quantitative. Each of the aforementioned
11 modeling methods has its unique superiority in dealing with the problems on different scales, whereas they also
12 possess different limitations in the membrane formation mechanism exploration respectively: the membrane
13 morphologies and performance cannot be directly predicted only by the latest continuum models; Also, variation
14 in the molecular structure is hardly added into the continuum and PF models; Finally, the temporospatial scale
15 of the molecule-based simulation is still far from the real membrane manufacture procedure. Considering
16 membrane formation is a complicated multiscale process, currently, people are starting to combine the different
17 models mentioned above together to obtain more precise results and improve predictability. For example, Tang
18 et al. adopted cooling rates obtained by the macroscopic mathematical models to determine the simulation steps
19 of different parts of the system to study the TIPS process, which can overcome the thermostatic limitation of the
20 DPD method [191]. Also, Li et al. firstly conducted DPD simulations to obtain the porous PMP membrane
21 structure and then applied MD simulations to investigate the gas permeability of the resultant membranes [194].
22 Additionally, Hopp-Hirschler et al. employed the nonsolvent composition profile to identify the boundary of the
23 phase field model and study the mechanism of membrane pore formation via the NIPS process [177].
24 The second thing to be fixed is until now for sake of reducing derivation and calculation cost, most of the
25 previous models just focused on a hypothetical polymer system instead of a real membrane formation system.
26 Only several classic polymer/solvent/nonsolvent systems like PVDF/DPC [106], PU/DMF/water [134], and
27 CA/acetone/non-volatile [81, 82] systems have been detailly analyzed, as listed by Table 1. Whereas numerous
28 polymer systems have been applied to fabrication membranes in labs and industries, and the resultant membranes
29 present various structures and performance. Differences between the membrane formation processes of the
30 various complicated systems and that of the classic systems remain ambiguous, which obviously block people
31 to realize controllable preparation of the membranes with higher performance. Next, people should develop the
32 existing models to map these real and novel polymer systems to produce results that can give direct instructions
33 to membrane production. Besides, some hydrophilic polymers such as PVP, PEG (polyethylene glycol), and
34 PMMA (polymethyl methacrylate) have been added into some classic polymer systems to improve the
35 hydrophilicity and water permeability of the membranes [210-212]. Whereas up to now, only a few modeling
36 studies are available to study the two-polymer system [73, 199], and this should be attributed to the difficulty in
37 characterizing the interaction and getting proper parameters of two polymers. People basically rely on many
38 experiments to get the parameters, while the same parameter that comes from different experiments sometimes
39 always cannot agree with each other. So, in the future, the problem may be relaxed by the wide application of
40 powerful and large-scale simulations. People will obtain reliable and consistent results based on force fields
41 derived from atomistic interaction.

27
 1 4. Conclusions

2 This review presented a full and integrated summary of the literature in the macroscale, meso- and
3 molecular-scale modeling approaches of membrane formation by phase inversion over nearly fifty years. It
4 began with a basic description of different modeling approaches on membrane formation via the phase inversion
5 technique. Secondly, based on an introduction to the principal conception of the thermodynamics of the polymer
6 solution, fundamental information of the dynamic modeling tools including macroscale continuum transfer
7 models, mesoscopic PF theory as well as molecular-scale simulations on the membrane formation have been
8 summarized. Subsequently, the model development and the applications of the modeling work on the three scales
9 have been presented entirely. Finally, current challenges and outlooks are proposed and discussed.
10 The continuum transfer model can produce concentration and temperature profile evolution of the whole
11 system over the full time of a process. These profiles are then superimposed as composition paths on
12 thermodynamic phase diagrams to give heuristic concepts about the membrane pore structure, which agree fairly
13 with the experimental data. The PF model predicts morphological and structural evolution during the phase
14 separation on broader time and length scales, and it may play a key role in making a connection between the
15 molecular-based methods and the macroscale transfer models in the future. The molecular simulations can
16 involve the microscale molecular structure and outputs are direct representations of structures, which can be
17 analyzed in desired properties. The results of these different models have provided a comprehensive
18 understanding of the membrane formation mechanism via phase inversion and helped the researchers to design
19 more advanced membrane fabrication techniques by combining the driving forces of different phase inversion
20 processes.
21 However, current limitations and challenges also exist. The continuum models have received the most
22 attention, therefore it is much more developed and powerful in this area. However, more rigorous thermo-kinetic
23 models that don’t only deal with mass transfer within a homogenous region, but also can predict the composition
24 change within the phase separation regions, should be developed in the future. The PF tends to be strong on the
25 growth and coarsening, but not so much on the molecular level and the fine structure. Besides, molecular-scale
26 simulations are limited by time and length scales, and this limitation is disappearing gradually with the rapid
27 development of modern computational power. Based on the present modeling research, we can predict that in
28 the next few years, firstly appropriate multiscale modeling with efficient coarse-grained potentials obtained from
29 the molecular tools and covering advantages of large scale modeling methods may be developed to reproduce
30 the pore structures on certain systems of industrial interest; secondly, more model works considering real systems
31 will emerge to provide practical instructions for the industrial field.
32
33 Acknowledgments
34
35 Funding for this work was provided by the National Science Foundation Industry / University Cooperative
36 Research Center for Membrane Science, Engineering and Technology (MAST) (NSF IIP 1361809) and the
37 University of Arkansas. This research is also supported by the Open Project of the State Key Laboratory of
38 Chemical Engineering (No. SKL-ChE-19A02) and the Fundamental Research Funds for the Central Universities
39 (No. 2020YQHH05 and No. 2021YJSHH25).
40

41 References

42 [1] M. Mulder, Phase inversion membranes, in: Encyclopedia of Separation Science,

43 Academic Press, Oxford, 2000, pp. 3331-3346.
44 [2] G.R. Guillen, Y.J. Pan, M.H. Li, E.M.V. Hoek, Preparation and characterization of
45 membranes formed by nonsolvent induced phase separation: A review, Industrial & Engineering
46 Chemistry Research, 50 (2011) 3798-3817.
47 [3] B.S. Lalia, V. Kochkodan, R. Hashaikeh, N. Hilal, A review on membrane fabrication:
48 Structure, properties and performance relationship, Desalination, 326 (2013) 77-95.

28
 1 [4] A. Venault, Y. Chang, D.-M. Wang, D. Bouyer, A review on polymeric membranes and
2 hydrogels prepared by vapor-induced phase separation process, Polymer Reviews, 53 (2013) 568-
3 626.
4 [5] D.M. Wang, J.Y. Lai, Recent advances in preparation and morphology control of
5 polymeric membranes formed by nonsolvent induced phase separation, Current opinion in
6 chemical engineering, 2 (2013) 229-237.
7 [6] J.F. Kim, J.H. Kim, Y.M. Lee, E. Drioli, Thermally induced phase separation and
8 electrospinning methods for emerging membrane applications: A review, AIChE Journal, 62
9 (2016) 461-490.
10 [7] M. Liu, S.H. Liu, Z.L. Xu, Y.M. Wei, H. Yang, Formation of microporous polymeric
11 membranes via thermally induced phase separation: A review, Frontiers of Chemical Science and
12 Engineering, 10 (2016) 57-75.
13 [8] A. Kausar, Phase inversion technique-based polyamide films and their applications: A
14 comprehensive review, Polymer-plastics technology and engineering, 56 (2017) 1421-1437.
15 [9] S.P. Nunes, P.Z. Culfaz-Emecen, G.Z. Ramon, T. Visser, G.H. Koops, W. Jin, M.
16 Ulbricht, Thinking the future of membranes: Perspectives for advanced and new membrane
17 materials and manufacturing processes, Journal of Membrane Science, 598 (2020).
18 [10] B. Wang, Z. Lai, Finger-like voids induced by viscous fingering during phase inversion
19 of alumina/PES/NMP suspensions, Journal of Membrane Science, 405 (2012) 275-283.
20 [11] M.M. Pendergast, E.M.V. Hoek, A review of water treatment membrane
21 nanotechnologies, Energy & Environmental Science, 4 (2011) 1946-1971.
22 [12] F. Liu, N.A. Hashim, Y. Liu, M.R.M. Abed, K. Li, Progress in the production and
23 modification of PVDF membranes, Journal of Membrane Science, 375 (2011) 1-27.
24 [13] Y. Lin, Y. Tang, H. Ma, J. Yang, Y. Tian, W. Ma, X. Wang, Formation of a bicontinuous
25 structure membrane of polyvinylidene fluoride in diphenyl carbonate diluent via thermally induced
26 phase separation, Journal of Applied Polymer Science, 114 (2009) 1523-1528.
27 [14] P.L. Hanks. Post coarsening effects on membrane microstructure. The University of
28 Texas at Austin, 2008.
29 [15] Q. Li, Z.L. Xu, L.Y. Yu, Effects of mixed solvents and PVDF types on performances of
30 PVDF microporous membranes, Journal of Applied Polymer Science, 115 (2010) 2277-2287.
31 [16] G.R. Guillen, G.Z. Ramon, H.P. Kavehpour, R.B. Kaner, E.M. Hoek, Direct microscopic
32 observation of membrane formation by nonsolvent induced phase separation, Journal of
33 Membrane Science, 431 (2013) 212-220.
34 [17] P. Graham, A. McHugh, Kinetics of thermally induced phase separation in a
35 crystallizable polymer solution, Macromolecules, 31 (1998) 2565-2568.
36 [18] L. Masaro, X. Zhu, Physical models of diffusion for polymer solutions, gels and solids,
37 Progress in Polymer Science, 24 (1999) 731-775.
38 [19] G.D. Phillies, Universal scaling equation for self-diffusion by macromolecules in
39 solution, Macromolecules, 19 (1986) 2367-2376.
40 [20] X. Li, C. Chen, J. Li, Formation kinetics of polyethersulfone with cardo membrane via
41 phase inversion, Journal of Membrane Science, 314 (2008) 206-211.
42 [21] H. Matsuyama, S. Kudari, H. Kiyofuji, Y. Kitamura, Kinetic studies of thermally induced
43 phase separation in polymer–diluent system, Journal of Applied Polymer Science, 76 (2000) 1028-
44 1036.
45 [22] P. Menut, Y.S. Su, W. Chinpa, C. Pochat-Bohatier, A. Deratani, D.M. Wang, P. Huguet,
46 C.Y. Kuo, J.Y. Lai, C. Dupuy, A top surface liquid layer during membrane formation using vapor-

29
 1 induced phase separation (VIPS) - Evidence and mechanism of formation, Journal of Membrane
2 Science, 310 (2008) 278-288.
3 [23] S.H. Yoo, C.K. Kim, Effects of the diluent mixing ratio and conditions of the thermally
4 induced phase-separation process on the pore size of microporous polyethylene membranes,
5 Journal of Applied Polymer Science, 108 (2008) 3154-3162.
6 [24] J.S. Taurozzi, C.A. Crock, V.V. Tarabara, C-60-polysulfone nanocomposite membranes:
7 Entropic and enthalpic determinants of C-60 aggregation and its effects on membrane properties,
8 Desalination, 269 (2011) 111-119.
9 [25] H.C. Yang, Q.Y. Wu, H.Q. Liang, L.S. Wan, Z.K. Xu, Thermally induced phase
10 separation of poly (vinylidene fluoride)/diluent systems: Optical microscope and infrared
11 spectroscopy studies, Journal of Polymer Science Part B: Polymer Physics, 51 (2013) 1438-1447.
12 [26] T. Ishigami, K. Nakatsuka, Y. Ohmukai, E. Kamio, T. Maruyama, H. Matsuyama,
13 Solidification characteristics of polymer solution during polyvinylidene fluoride membrane
14 preparation by nonsolvent-induced phase separation, Journal of Membrane Science, 438 (2013)
15 77-82.
16 [27] T. Ishigami, Y. Nii, Y. Ohmukai, S. Rajabzadeh, H. Matsuyama, Solidification Behavior
17 of Polymer Solution during Membrane Preparation by Thermally Induced Phase Separation,
18 Membranes, 4 (2014) 113-122.
19 [28] M. Shang, H. Matsuyama, M. Teramoto, D.R. Lloyd, N. Kubota, Preparation and
20 membrane performance of poly(ethylene-co-vinyl alcohol) hollow fiber membrane via thermally
21 induced phase separation, Polymer, 44 (2003) 7441-7447.
22 [29] W.E.S. Robert Byron Bird, Edwin N. Lightfoot, Transport phenomena, Rev. 2nd ed ed.,
23 New York : John Wiley & Sons, 2007.
24 [30] L. Euler, Principes généraux du mouvement des fluides, Mémoires de l'académie des
25 sciences de Berlin, 11 (1757) 274-315.
26 [31] S. Chandrasekhar, Stochastic problems in physics and astronomy, Reviews of modern
27 physics, 15 (1943) 1.
28 [32] L.P. Pierre Degond, Giovanni Russo, Modeling and Simulation in Science, Engineering
29 and Technology, 1 ed., Birkhäuser, Boston, MA, Toulouse cedex, 2004.
30 [33] W.T. E.G.D. Cohen, The Boltzmann equation theory and applications, Springer-Verlag
31 Wien, New York, 1973.
32 [34] P.K. Grzegorz Łukaszewicz, Navier–Stokes Equations An Introduction with
33 Applications, Springer International Publishing Switzerland, 2016.
34 [35] W.L. Luyben, process modeling, simulation, and control for chemical engineers, 2nd ed.
35 ed., McGraw-Hill publisbing company, 1989.
36 [36] M.A. Khansary, A. Marjani, S. Shirazian, On the search of rigorous thermo-kinetic
37 model for wet phase inversion technique, Journal of Membrane Science, 538 (2017) 18-33.
38 [37] J.W. Gibbs, Sciences, A method of geometrical representation of the thermodynamic
39 properties by means of surfaces, Transactions of Connecticut Academy of Arts, 2 (1873) 382-404.
40 [38] L.D.L.D. Ter-Haar, Collected Papers of L.D. Landau, Pergamon, 1965.
41 [39] H. Emmerich, The diffuse interface approach in materials science : thermodynamic
42 concepts and applications of phase-field models, Springer-Verlag, New York 2003.
43 [40] L. Onsager, Reciprocal relations in irreversible processes. I, Physical Review, 37 (1931)
44 405-426.
45 [41] L. Onsager, Reciprocal relations in irreversible processes. II, Physical Review, 38 (1931)
46 2265-2279.

30
 1 [42] J.S. Rowlinson, Translation of J.D. van der Waals' "the thermodynamic theory of
2 capillarity under the hypothesis of a continuous variation of density", Journal of Statistical Physics,
3 20 (1979) 197-244.
4 [43] J.W. Cahn, J.E. Hilliard, Free energy of a nonuniform system. I. Interfacial free energy,
5 The Journal of chemical physics, 28 (1958) 258-267.
6 [44] J.W.C. J.E. Hilliard, On the nature of the interface between a solid metal and its melt,
7 Acta Metallurgica, 6 (1958) 772-774.
8 [45] J.W. Cahn, Free energy of a nonuniform system. 2. thermodynamic basis, Journal of
9 Chemical Physics, 30 (1959) 1121-1124.
10 [46] J.W. Cahn, J.E. Hilliard, Free energy of a nonuniform system. 3. nucleation in a 2-
11 component incompressible fluid, Journal of Chemical Physics, 31 (1959) 688-699.
12 [47] S.M. Allen, J.W. Cahn, Ground state structures in ordered binary-alloys with second
13 neighbor interactions, Acta Metallurgica, 20 (1972) 423-&.
14 [48] S.M. Allen, J.W. Cahn, Correction to ground-state of FCC binary ordered alloys with
15 first and second neighbor pairwise interactions, Scripta Metallurgica, 7 (1973) 1261-1264.
16 [49] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, E. Teller, Equation of
17 state calculations by fast computing machines, Journal of Chemical Physics, 21 (1953) 1087-1092.
18 [50] R.D. Groot, P.B. Warren, Dissipative particle dynamics: bridging the gap between
19 atomistic and mesoscopic simulation, Journal of Chemical Physics, 107 (1997) 4423.
20 [51] G.T. Caneba, Analysis of polymer membrane formation through spinodal
21 decomposition, Polymer Engineering and Science, 31 (1991) 879-885.
22 [52] B. Barton, A. McHugh, Kinetics of thermally induced phase separation in ternary
23 polymer solutions. I. Modeling of phase separation dynamics, Journal of Polymer Science Part B
24 Polymer Physics, 37 (1999) 1449-1460.
25 [53] R. Larson, L. Scriven, H. Davis, Monte Carlo simulation of model amphiphile‐oil–
26 water systems, The Journal of chemical physics, 83 (1985) 2411-2420.
27 [54] X.-L. Wang, H.-J. Qian, L.-J. Chen, Z.-Y. Lu, Z.-S. Li, Dissipative particle dynamics
28 simulation on the polymer membrane formation by immersion precipitation, Journal of Membrane
29 Science, 311 (2008) 251-258.
30 [55] M.R. Mruzik, F.F. Abraham, G.M. Pound, Phase separation in fluid systems by spinodal
31 decomposition. II. A molecular dynamics computer simulation, The Journal of chemical physics,
32 69 (1978) 3462-3467.
33 [56] D.R. Lloyd, K.E. Kinzer, H. Tseng, Microporous membrane formation via thermally
34 induced phase separation. I. Solid-liquid phase separation, Journal of Membrane Science, 52
35 (1990) 239-261.
36 [57] D.R. Lloyd, S.S. Kim, K.E. Kinzer, Microporous membrane formation via thermally-
37 induced phase separation. II. Liquid—liquid phase separation, Journal of Membrane Science, 64
38 (1991) 1-11.
39 [58] S.S. Kim, D.R. Lloyd, Microporous membrane formation via thermally-induced phase
40 separation. III. Effect of thermodynamic interactions on the structure of isotactic polypropylene
41 membranes, Journal of Membrane Science, 64 (1991) 13-29.
42 [59] S.S. Kim, G. Lim, A.A. Alwattari, Y.F. Wang, D.R. Lloyd, Microporous membrane
43 formation via thermally-induced phase separation. V. Effect of diluent mobility and crystallization
44 on the structure of isotactic polypropylene membranes, Journal of Membrane Science, 64 (1991)
45 41-53.

31
 1 [60] G. Lim, S.S. Kim, Q. Ye, Y.F. Wang, D.R. Lloyd, Microporous membrane formation
2 via thermally-induced phase separation. IV. Effect of isotactic polypropylene crystallization
3 kinetics on membrane structure, Journal of Membrane Science, 64 (1991) 31-40.
4 [61] A.A. Alwattari, D.R. Lloyd, Microporous membrane formation via thermally-induced
5 phase separation. VI. Effect of diluent morphology and relative crystallization kinetics on
6 polypropylene membrane structure, Journal of Membrane Science, 64 (1991) 55-67.
7 [62] S.S. Kim, D.R. Lloyd, Thermodynamics of polymer/diluent systems for thermally
8 induced phase separation: 1. Determination of equation of state parameters, Polymer, 33 (1992)
9 1026-1035.
10 [63] S.S. Kim, D.R. Lloyd, Thermodynamics of polymer/diluent systems for thermally
11 induced phase separation: 2. Solid-liquid phase separation systems, Polymer, 33 (1992) 1036-
12 1046.
13 [64] S.S. Kim, D.R. Lloyd, Thermodynamics of polymer/diluent systems for thermally
14 induced phase separation: 3. Liquid-liquid phase separation systems, Polymer, 33 (1992) 1047-
15 1057.
16 [65] C. Qian, S.J. Mumby, B.E. Eichinger, Phase Diagrams of Binary Polymer Solutions and
17 Blends, Macromolecules, 24 (1991) 1655-1661.
18 [66] W.R. Burghardt, Phase Diagrams for Binary Polymer Systems Exhibiting both
19 Crystallization and Limited Liquid-Liquid Miscibility, Macromolecules, 22 (1989) 2482-2486.
20 [67] K.S. McGuire, A. Laxminarayan, D.R. Lloyd, A simple method of extrapolating the
21 coexistence curve and predicting the melting point depression curve from cloud point data for
22 polymer-diluent systems, Polymer, 35 (1994) 4404-4407.
23 [68] P.J. Flory, Principles of Polymer Chemistry, George Banta Publishing Company,
24 Menasha Wisconsin, 1953.
25 [69] D. Bouyer, O. M'Barki, C. Pochat-Bohatier, C. Faur, E. Petit, P. Guenoun, Modeling the
26 membrane formation of novel PVA membranes for predicting the composition path and their final
27 morphology, AIChE Journal, 63 (2017) 3035-3047.
28 [70] P.T.P. Aryanti, D. Ariono, A.N. Hakim, I.G. Wenten, Flory-Huggins Based Model to
29 Determine Thermodynamic Property of Polymeric Membrane Solution, Journal of Physics:
30 Conference Series, 1090 (2018) 012074 (012013 pp.)-012074 (012013 pp.).
31 [71] R. Boom, T. Van den Boomgaard, C. Smolders, Equilibrium thermodynamics of a
32 quaternary membrane-forming system with two polymers. 1. Calculations, Macromolecules, 27
33 (1994) 2034-2040.
34 [72] R. Horst, B. Wolf, Phase diagrams calculated for quaternary polymer blends, The Journal
35 of chemical physics, 103 (1995) 3782-3787.
36 [73] V.P. Khare, A.R. Greenberg, W.B. Krantz, Vapor-induced phase separation - effect of
37 the humid air exposure step on membrane morphology Part I. Insights from mathematical
38 modeling, Journal of Membrane Science, 258 (2005) 140-156.
39 [74] J.T. Jung, J.F. Kim, H.H. Wang, E. Di Nicolo, E. Drioli, Y.M. Lee, Understanding the
40 non-solvent induced phase separation (NIPS) effect during the fabrication of microporous PVDF
41 membranes via thermally induced phase separation (TIPS), Journal of Membrane Science, 514
42 (2016) 250-263.
43 [75] A. McHugh, L. Yilmaz, The diffusion equations for polymer membrane formation in
44 ternary systems, Journal of Polymer Science Part B: Polymer Physics, 23 (1985) 1271-1274.

32
 1 [76] F.W. Altena, C.A. Smoulders, Calculation of Liquid-Liquid Phase Separation in a
2 Ternary System of a Polymer in a Mixture of a Solvent and a Nonsolvent, Macromolecules, 15
3 (1982) 1491-1497.
4 [77] J.-Y. Lai, S.-F. Lin, F.-C. Lin, D.-M. Wang, Construction of Ternary Phase Diagrams in
5 Nonsolvent/Solvent/PMMA Systems, Journal of Polymer Science Part B: Polymer Physics, 36
6 (1997) 607-615.
7 [78] P.v.d. Witte, P.J. Dijkstra, J.W.A.v.d. Berg, J. Feijen, Metastable Liquid-Liquid and
8 Solid-Liquid Phase Boundaries in Polymer-Solvent-Nonsolvent Systems, Journal of Polymer
9 Science Part B: Polymer Physics, 35 (1997) 763-770.
10 [79] K. Langenbach, M. Fischlschweiger, S. Enders, Prediction of the solid–liquid–liquid
11 equilibria of linear and branched semi-crystalline poly-ethylene in solutions of diphenyl ether by
12 Lattice Cluster Theory, Molecular Physics, 114 (2016) 2717-2723.
13 [80] C. Tsay, A. McHugh, Mass transfer modeling of asymmetric membrane formation by
14 phase inversion, Journal of Polymer Science Part B: Polymer Physics, 28 (1990) 1327-1365.
15 [81] S.S. Shojaie, W.B. Krantz, A.R. Greenberg, Dense polymer film and membrane
16 formation via the dry-cast process Part I. Model development, Journal of Membrane Science, 94
17 (1994) 255-280.
18 [82] S.S. Shojaie, W.B. Krantz, A.R. Greenberg, Dense Polymer Film and Membrane
19 Formation via the dry-cast process Part II. Model validation and morphological studies, Journal of
20 Membrane Science, 94 (1994) 281-298.
21 [83] H.h. Lin, Y.h. Tang, T.y. Liu, H. Matsuyama, X.-l. Wang, Understanding the thermally
22 induced phase separation process via a Maxwell–Stefan model, Journal of Membrane Science, 507
23 (2016) 143-153.
24 [84] L.Q. Chen, Phase-field models for microstructure evolution, Annual review of materials
25 research, 32 (2002) 113-140.
26 [85] W.J. Boettinger, J.A. Warren, C. Beckermann, A. Karma, Phase-field simulation of
27 solidification, Annual review of materials research, 32 (2002) 163-194.
28 [86] I. Steinbach, Phase-field models in materials science, Modeling and Simulation in
29 Materials Science and Engineering, 17 (2009) 1-31.
30 [87] S. Asai, S. Majumdar, A. Gupta, K. Kargupta, S. Ganguly, Dynamics and pattern
31 formation in thermally induced phase separation of polymer–solvent system, Computational
32 Materials Science, 47 (2009) 193-205.
33 [88] L.T. Yan, X.M. Xie, Numerical simulation of substrate effects on spinodal
34 decomposition in polymer binary mixture: Effects of the surface potential, Polymer, 47 (2006)
35 6472-6480.
36 [89] K.W.D. Lee, P.K. Chan, X. Feng, A computational study into thermally induced phase
37 separation in polymer solutions under a temperature gradient, Macromolecular Theory and
38 Simulations, 11 (2002) 996-1005.
39 [90] N. Provatas, K. Elder, Phase-field methods in materials science and engineering, John
40 Wiley & Sons, 2011.
41 [91] H. Emmerich, The diffuse interface approach in materials science: thermodynamic
42 concepts and applications of phase-field models, Springer Science & Business Media, 2003.
43 [92] E.B. Nauman, D.Q. He, Nonlinear diffusion and phase separation, Chemical Engineering
44 Science, 56 (2001) 1999-2018.
45 [93] S.-U. Hong, Prediction of polymer/solvent diffusion behavior using free-volume theory,
46 Industrial & Engineering Chemistry Research, 34 (1995) 2536-2544.

33
 1 [94] J.M. Zielinski, J. Duda, Predicting polymer/solvent diffusion coefficients using free‐
2 volume theory, AIChE Journal, 38 (1992) 405-415.
3 [95] J. Vrentas, J. Duda, H.C. Ling, Self‐diffusion in polymer‐solvent‐solvent systems,
4 Journal of Polymer Science Part B: Polymer Physics, 22 (1984) 459-469.
5 [96] P. De Gennes, Dynamics of entangled polymer solutions. I. The Rouse model,
6 Macromolecules, 9 (1976) 587-593.
7 [97] E.J. Kramer, P. Green, C.J.J.p. Palmstrøm, Interdiffusion and marker movements in
8 concentrated polymer-polymer diffusion couples, 25 (1984) 473-480.
9 [98] J. Vrentas, J. Duda, Diffusion in infinitely dilute polystyrene solutions, Journal of
10 Applied Polymer Science, 20 (1976) 1125-1131.
11 [99] J. Vrentas, J. Duda, Diffusion in polymer–solvent systems. II. A predictive theory for
12 the dependence of diffusion coefficients on temperature, concentration, and molecular weight,
13 Journal of Polymer Science Part B: Polymer Physics, 15 (1977) 417-439.
14 [100] J. Vrentas, J. Duda, Diffusion in polymer—solvent systems. I. Reexamination of the
15 free‐volume theory, Journal of Polymer Science Part B: Polymer Physics, 15 (1977) 403-416.
16 [101] J. Vrentas, C. Vrentas, A new equation relating self-diffusion and mutual diffusion
17 coefficients in polymer-solvent systems, Macromolecules, 26 (1993) 6129-6131.
18 [102] J. Vrentas, J. Duda, H.C. Ling, Free‐volume theories for self‐diffusion in polymer–
19 solvent systems. I. Conceptual differences in theories, Journal of Polymer Science Part B: Polymer
20 Physics, 23 (1985) 275-288.
21 [103] J. Vrentas, J. Duda, H.C. Ling, A.C. Hou, Free‐volume theories for self‐diffusion in
22 polymer–solvent systems. II. Predictive capabilities, Journal of Polymer Science Part B: Polymer
23 Physics, 23 (1985) 289-304.
24 [104] L. Gránásy, T. Börzsönyi, T. Pusztai, Nucleation and bulk crystallization in binary
25 phase field theory, Physical Review Letters, 88 (2002) 206105.
26 [105] H. Manzanarez, J. Mericq, P. Guenoun, J. Chikina, D. Bouyer, Modeling phase
27 inversion using Cahn-Hilliard equations–Influence of the mobility on the pattern formation
28 dynamics, J Chemical Engineering Science, 173 (2017) 411-427.
29 [106] M.R. Cervellere, Y.-h. Tang, X. Qian, D.M. Ford, P.C.J.J.o.M.S. Millett, Mesoscopic
30 Simulations of Thermally-Induced Phase Separation in PVDF/DPC Solutions, (2019).
31 [107] K.-W.D. Lee, P.K. Chan, X. Feng, Morphology development and characterization of
32 the phase-separated structure resulting from the thermal-induced phase separation phenomenon in
33 polymer solutions under a temperature gradient, Chemical Engineering Science, 59 (2004) 1491-
34 1504.
35 [108] M.D. Vo, D.V. Papavassiliou, Physical adsorption of polyvinyl pyrrolidone on carbon
36 nanotubes under shear studied with dissipative particle dynamics simulations, Carbon, 100 (2016)
37 291-301.
38 [109] J.C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R.D.
39 Skeel, L. Kale, K. Schulten, Scalable molecular dynamics with NAMD, Journal of Computational
40 Chemistry, 26 (2005) 1781-1802.
41 [110] Y. Wang, C.B. Harrison, K. Schulten, J.A. McCammon, Implementation of accelerated
42 molecular dynamics in NAMD, Computational science & discovery, 4 (2011) 015002.
43 [111] B. Prathab, V. Subramanian, T. Aminabhavi, Molecular dynamics simulations to
44 investigate polymer–polymer and polymer–metal oxide interactions, Polymer, 48 (2007) 409-416.

34
 1 [112] K. Miyamoto, D. Ikeshima, T. Furutani, H. Xiao, A. Yonezu, X. Chen, On the surface
2 hydrophilization of a blended polysulfone membrane: atomic force microscopy measurement and
3 molecular dynamics simulation, Surface Topography-Metrology and Properties, 7 (2019).
4 [113] C. Fang, S. Rajabzadeh, H.-C. Wu, X. Zhang, N. Kato, M. Kunimatsu, T. Yoshioka, H.
5 Matsuyama, Effect of mass transfer at the interface of the polymer solution and extruded solvent
6 during the air gap on membrane structures and performances in TIPS process using triple-orifice
7 spinneret, Journal of Membrane Science, 595 (2020).
8 [114] P.J. Hoogerbrugge, J.M.V.A. Koelman, Simulating Microscopic Hydrodynamic
9 Phenomena with Dissipative Particle Dynamics, Europhysics Letters, 19 (1992) 155-160.
10 [115] P. Espanol, P.B. Warren, Perspective: Dissipative particle dynamics, Journal of
11 Chemical Physics, 146 (2017).
12 [116] Z.Y. Lu, Y.L. Wang, An introduction to dissipative particle dynamics, Methods in
13 molecular biology (Clifton, N.J.), 924 (2013) 617-633.
14 [117] C.M. Hansen, Hansen solubility parameter a user’s handbook, Second Edition, CRC
15 Press, Boca Raton, 2007.
16 [118] Y. Tang, M. Li, Y. Lin, L. Wang, F. Wu, X. Wang, A novel green diluent for the
17 preparation of poly (4-methyl-1-pentene) membranes via a thermally-induced phase separation
18 method, Membranes, 11 (2021) 622.
19 [119] J. Huo, Z. Chen, J. Zhou, Zwitterionic Membrane via Nonsolvent Induced Phase
20 Separation: A Computer Simulation Study, Langmuir, 35 (2019) 1973-1983.
21 [120] P.L. Hanks, D.R. Lloyd, Deterministic model for matrix solidification in liquid–liquid
22 thermally induced phase separation, Journal of Membrane Science, 306 (2007) 125-133.
23 [121] I. Carmesin, K. Kremer, The bond fluctuation method: a new effective algorithm for
24 the dynamics of polymers in all spatial dimensions, Macromolecules, 21 (1988) 2819-2823.
25 [122] C. Cohen, G. Tanny, S. Prager, Diffusion‐controlled formation of porous structures
26 in ternary polymer systems, Journal of Polymer Science Part B: Polymer Physics, 17 (1979) 477-
27 489.
28 [123] J. Wijmans, F. Altena, C. Smolders, Diffusion during the immersion precipitation
29 process, Journal of Polymer Science: Polymer Physics Edition, 22 (1984) 519-524.
30 [124] L. Yilmaz, A. McHugh, Analysis of nonsolvent–solvent–polymer phase diagrams and
31 their relevance to membrane formation modeling, Journal of Applied Polymer Science, 31 (1986)
32 997-1018.
33 [125] L. Yilmaz, A. McHugh, Modelling of asymmetric membrane formation. I. Critique of
34 evaporation models and development of a diffusion equation formalism for the quench period,
35 Journal of Membrane Science, 28 (1986) 287-310.
36 [126] A. Reuvers, C. Smolders, Formation of membranes by means of immersion
37 precipitation: Part II. the mechanism of formation of membranes prepared from the system
38 cellulose acetate-acetone-water, Journal of Membrane Science, 34 (1987) 67-86.
39 [127] A. Reuvers, J. Van den Berg, C. Smolders, Formation of membranes by means of
40 immersion precipitation: Part I. A model to describe mass transfer during immersion precipitation,
41 Journal of Membrane Science, 34 (1987) 45-65.
42 [128] P. Radovanovic, S.W. Thiel, S.-T. Hwang, Formation of asymmetric polysulfone
43 membranes by immersion precipitation. Part I. Modelling mass transport during gelation, Journal
44 of Membrane Science, 65 (1992) 213-229.
45 [129] P. Radovanovic, S.W. Thiel, S.-T. Hwang, Formation of asymmetric polysulfone
46 membranes by immersion precipitation. Part II. The effects of casting solution and gelation bath

35
 1 compositions on membrane structure and skin formation, Journal of Membrane Science, 65 (1992)
2 231-246.
3 [130] L.P. Cheng, Y.S. Soh, A.H. Dwan, C.C. Gryte, An improved model for mass transfer
4 during the formation of polymeric membranes by the immersion‐precipitation process, Journal
5 of Polymer Science Part B: Polymer Physics, 32 (1994) 1413-1425.
6 [131] A. McHugh, C. Tsay, B. Barton, J. Reeve, Comments on a model for mass transfer
7 during phase inversion, Journal of Polymer Science Part B: Polymer Physics, 33 (1995) 2175-
8 2179.
9 [132] G. Fernandes, J. Pinto, R. Nobrega, Modeling and simulation of the phase‐inversion
10 process during membrane preparation, Journal of Applied Polymer Science, 82 (2001) 3036-3051.
11 [133] S.K. Karode, A. Kumar, Formation of polymeric membranes by immersion
12 precipitation: an improved algorithm for mass transfer calculations, Journal of Membrane Science,
13 187 (2001) 287-296.
14 [134] Y.D. Kim, J.Y. Kim, H.K. Lee, S.C. Kim, A new modeling of asymmetric membrane
15 formation in rapid mass transfer system, Journal of Membrane Science, 190 (2001) 69-77.
16 [135] H. Lee, W.B. Krantz, S.-T. Hwang, A model for wet-casting polymeric membranes
17 incorporating nonequilibrium interfacial dynamics, vitrification and convection, Journal of
18 Membrane Science, 354 (2010) 74-85.
19 [136] J. Anderson, R. Ullman, Mathematical analysis of factors influencing the skin thickness
20 of asymmetric reverse osmosis membranes, Journal of Applied Physics, 44 (1973) 4303-4311.
21 [137] C. Castellari, S. Ottani, Preparation of reverse osmosis membranes. A numerical
22 analysis of asymmetric membrane formation by solvent evaporation from cellulose acetate casting
23 solutions, Journal of Membrane Science, 9 (1981) 29-41.
24 [138] W.B. Krantz, R.J. Ray, R.L. Sani, Theoretical study of the transport processes occurring
25 during the evaporation step in asymmetric membrane casting, Journal of Membrane Science, 29
26 (1986) 11-36.
27 [139] C. Tsay, A. McHugh, Mass transfer dynamics of the evaporation step in membrane
28 formation by phase inversion, Journal of Membrane Science, 64 (1991) 81-92.
29 [140] L. Tan, W.B. Krantz, A.R. Greenberg, R.L. Sani, Studies of convective transport in
30 evaporative casting of dense polymer films, Journal of Membrane Science, 108 (1995) 245-255.
31 [141] H. Matsuyama, M. Teramoto, T. Uesaka, Membrane formation and structure
32 development by dry-cast process, Journal of Membrane Science, 135 (1997) 271-288.
33 [142] H. Matsuyama, M. Teramoto, R. Nakatani, T. Maki, Membrane formation via phase
34 separation induced by penetration of nonsolvent from vapor phase. I. Phase diagram and mass
35 transfer process, Journal of Applied Polymer Science, 74 (1999) 159-170.
36 [143] H. Matsuyama, M. Teramoto, R. Nakatani, T. Maki, Membrane formation via phase
37 separation induced by penetration of nonsolvent from vapor phase. II. Membrane morphology,
38 Journal of Applied Polymer Science, 74 (1999) 171-178.
39 [144] S.A. Altinkaya, B. Ozbas, Modeling of asymmetric membrane formation by dry-casting
40 method, Journal of Membrane Science, 230 (2004) 71-89.
41 [145] S.A. Altinkaya, H. Yenal, B. Ozbas, Membrane formation by dry-cast process: Model
42 validation through morphological studies, Journal of Membrane Science, 249 (2005) 163-172.
43 [146] S. Alsoy, J.L. Duda, Modeling of multicomponent drying of polymer films, AIChE
44 Journal, 45 (1999) 896-905.
45 [147] Y. Yip, A.J. McHugh, Modeling and simulation of nonsolvent vapor-induced phase
46 separation, Journal of Membrane Science, 271 (2006) 163-176.

36
 1 [148] H. Lee, S.R. Chaudhuri, W.B. Krantz, S.-T. Hwang, A model for evaporative casting
2 of polymeric membranes incorporating convection due to density changes, Journal of Membrane
3 Science, 284 (2006) 161-172.
4 [149] W.B. Krantz, A.R. Greenberg, D.J. Hellman, Dry-casting: Computer simulation,
5 sensitivity analysis, experimental and phenomenological model studies, Journal of Membrane
6 Science, 354 (2010) 178-188.
7 [150] D. Bouyer, W. Werapun, C. Pochat-Bohatier, A. Deratani, Morphological properties of
8 membranes fabricated by VIPS process using PEI/NMP/water system: SEM analysis and mass
9 transfer modelling, Journal of Membrane Science, 349 (2010) 97-112.
10 [151] D. Bouyer, C. Pochat‐Bohatier, Validation of mass‐transfer model for VIPS process
11 using in situ measurements performed by near‐infrared spectroscopy, AIChE Journal, 59 (2013)
12 671-686.
13 [152] H. Matsuyama, S. Berghmans, D.R. Lloyd, Formation of anisotropic membranes via
14 thermally induced phase separation, Polymer, 40 (1999) 2289-2301.
15 [153] H. Matsuyama, S. Berghmans, M.T. Batarseh, D.R. Lloyd, Effects of thermal history
16 on anisotropic and asymmetric membranes formed by thermally induced phase separation, Journal
17 of Membrane Science, 142 (1998) 27-42.
18 [154] P.M. Atkinson, D.R. Lloyd, Anisotropic flat sheet membrane formation via TIPS:
19 thermal effects, Journal of Membrane Science, 171 (2000) 1-18.
20 [155] P.M. Atkinson, D.R. Lloyd, Anisotropic flat sheet membrane formation via TIPS:
21 atmospheric convection and polymer molecular weight effects, Journal of Membrane Science, 175
22 (2000) 225-238.
23 [156] B. Barton, A. McHugh, Modeling the dynamics of membrane structure formation in
24 quenched polymer solutions, Journal of Membrane Science, 166 (2000) 119-125.
25 [157] D. Li. Modeling of the solid-liquid thermally induced-phase-separation (TIPS)
26 membrane formation process. University of Colorado at Boulder., 2003.
27 [158] D. Li, W.B. Krantz, A.R. Greenberg, R.L. Sani, Membrane formation via thermally
28 induced phase separation (TIPS): Model development and validation, Journal of Membrane
29 Science, 279 (2006) 50-60.
30 [159] R.J. Bearman, On the molecular basis of some current theories of diffusion1, The
31 Journal of Physical Chemistry, 65 (1961) 1961-1968.
32 [160] B.F. Barton, A.J. McHugh, Kinetics of thermally induced phase separation in ternary
33 polymer solutions. I. Modeling of phase separation dynamics, Journal of Polymer Science Part B-
34 Polymer Physics, 37 (1999) 1449-1460.
35 [161] J.S. Vrentas, J.L. Duda, Diffusion in Polymer-Solvent Systems. I. Reexamination of
36 the Free-Volume Theory, Journal of Polymer Science: Polymer Physics Edition, 15 (1977) 1125-
37 1131.
38 [162] P.E. Price Jr, I.H. Romdhane, Multicomponent diffusion theory and its applications to
39 polymer‐solvent systems, AIChE Journal, 49 (2003) 309-322.
40 [163] J.M. Zielinski, B.F. Hanley, Practical friction ‐ based approach to modeling
41 multicomponent diffusion, AIChE Journal, 45 (1999) 1-12.
42 [164] A. Chakrabarti, R. Toral, J.D. Gunton, M. Muthukumar, Spinodal decomposition in
43 polymer mixtures, Physical Review Letters, 63 (1989) 2072.
44 [165] G. Brown, A. Chakrabarti, Phase separation dynamics in off‐critical polymer blends,
45 The Journal of chemical physics, 98 (1993) 2451-2458.

37
 1 [166] T.L. Tran, P.K. Chan, D. Rousseau, Morphology control in symmetric polymer blends
2 using two-step phase separation, Computational Materials Science, 37 (2006) 328-335.
3 [167] Y. Li, R. Shi, C. Wang, X. Liu, Y. Wang, Phase-field simulation of thermally induced
4 spinodal decomposition in polymer blends, Modelling and Simulation in Materials Science and
5 Engineering, 20 (2012) 075002.
6 [168] G.T. Caneba, D.S. Soong, Polymer membrane formation through the thermal-inversion
7 process. 1. Experimental study of membrane structure formation, Macromolecules, 18 (1985)
8 2538-2545.
9 [169] G.T. Caneba, D.S. Soong, Polymer membrane formation through the thermal-inversion
10 process. 2. Mathematical modeling of membrane structure formation, Macromolecules, 18 (1985)
11 2545-2555.
12 [170] P.K. Chan, A.D. Rey, Computational analysis of spinodal decomposition dynamics in
13 polymer solutions, Macromolecular Theory and Simulations, 4 (1995) 873-899.
14 [171] B. Barton, P. Graham, A. McHugh, Dynamics of spinodal decomposition in polymer
15 solutions near a glass transition, Macromolecules, 31 (1998) 1672-1679.
16 [172] S.B. Kukadiya, P.K. Chan, M. Mehrvar, The Ludwig‐Soret Effect on the Thermally
17 Induced Phase Separation Process in Polymer Solutions: A Computational Study, Macromolecular
18 Theory and Simulations, 18 (2009) 97-107.
19 [173] P.K. Chan, Effect of concentration gradient on the thermal-induced phase separation
20 phenomenon in polymer solutions, Modelling and Simulation in Materials Science and
21 Engineering, 14 (2006) 41.
22 [174] Y. Mino, T. Ishigami, Y. Kagawa, H. Matsuyama, Three-dimensional phase-field
23 simulations of membrane porous structure formation by thermally induced phase separation in
24 polymer solutions, Journal of Membrane Science, 483 (2015) 104-111.
25 [175] R. Saxena, G.T. Caneba, Studies of spinodal decomposition in a ternary polymer‐
26 solvent‐nonsolvent system, Polymer Engineering & Science, 42 (2002) 1019-1031.
27 [176] B. Zhou, A.C. Powell, Phase field simulations of early stage structure formation during
28 immersion precipitation of polymeric membranes in 2D and 3D, Journal of Membrane Science,
29 268 (2006) 150-164.
30 [177] M. Hopp-Hirschler, U. Nieken, Modeling of pore formation in phase inversion
31 processes: Model and numerical results, Journal of Membrane Science, 564 (2018) 820-831.
32 [178] A. Akthakul, C.E. Scott, A.M. Mayes, A.J. Wagner, Lattice Boltzmann simulation of
33 asymmetric membrane formation by immersion precipitation, Journal of Membrane Science, 249
34 (2005) 213-226.
35 [179] D.R. Tree, K.T. Delaney, H.D. Ceniceros, T. Iwama, G.H. Fredrickson, A multi-fluid
36 model for microstructure formation in polymer membranes, Soft Matter, 13 (2017) 3013-3030.
37 [180] H. Xu, R. Matkar, T. Kyu, Phase-field modeling on morphological landscape of
38 isotactic polystyrene single crystals, Physical Review E, 72 (2005) 011804.
39 [181] M.A. Zaeem, S. Nouranian, M.F. Horstemeyer, Simulation of Polymer Crystal Growth
40 with Various Morphologies Using a Phase-Field Model, in: AIChE Annual Meeting, American
41 Institute of Chemical Engineers, Pittsburgh, PA, 2012.
42 [182] X.D. Wang, O. Jie, S. Jin, Z. Wen, A phase-field model for simulating various
43 spherulite morphologies of semi-crystalline polymers, Chinese Physics B, 22 (2013) 106103.
44 [183] X. Wang, J. Ouyang, W. Zhou, Z. Liu, A phase field technique for modeling and
45 predicting flow induced crystallization morphology of semi-crystalline polymers, Polymers, 8
46 (2016) 230.

38
 1 [184] G. Leptoukh, B. Strickland, C. Roland, Phase separation in two-dimensional fluid
2 mixtures, Physical Review Letters, 74 (1995) 3636-3639.
3 [185] H. Liu, A. Bhattacharya, A. Chakrabarti, Network domain structure in phase-separating
4 polymer solutions, Journal of Chemical Physics, 111 (1999) 11183-11191.
5 [186] M. Liu, H. Shen, X. Li, S. Wang, Molecular Simulation on Compatibility of
6 Polyurethane and Polyvinyl Chloride and Properties of PU/PVC Blend Membranes, Polymer
7 Materials Science & Engineering, 34 (2018) 114-120.
8 [187] Y.D. He, Y.H. Tang, X.L. Wang, Dissipative particle dynamics simulation on the
9 membrane formation of polymer–diluent system via thermally induced phase separation, Journal
10 of Membrane Science, 368 (2011) 78-85.
11 [188] Y.H. Tang, Y.D. He, X.L. Wang, Effect of adding a second diluent on the membrane
12 formation of polymer/diluent system via thermally induced phase separation: Dissipative particle
13 dynamics simulation and its experimental verification, Journal of Membrane Science, 409-410
14 (2012) 164-172.
15 [189] Y. Tang, Y. He, X. Wang, Three-dimensional analysis of membrane formation via
16 thermally induced phase separation by dissipative particle dynamics simulation, Journal of
17 Membrane Science, 437 (2013) 40-48.
18 [190] Y.H. Tang, Y.D. He, X.L. Wang, Investigation on the membrane formation process of
19 polymer-diluent system via thermally induced phase separation accompanied with mass transfer
20 across the interface: Dissipative particle dynamics simulation and its experimental verification,
21 Journal of Membrane Science, 474 (2015) 196-206.
22 [191] Y.H. Tang, H.H. Lin, T.Y. Liu, H. Matsuyama, X.L. Wang, Multiscale simulation on
23 the membrane formation process via thermally induced phase separation accompanied with heat
24 transfer, Journal of Membrane Science, 515 (2016) 258-267.
25 [192] Y.-h. Tang, E. Ledieu, M.R. Cervellere, P.C. Millett, D.M. Ford, X. Qian, Formation
26 of polyethersulfone membranes via nonsolvent induced phase separation process from dissipative
27 particle dynamics simulations, Journal of Membrane Science, 599 (2020).
28 [193] X. Huang, W.P. Wang, Z. Zheng, X.L. Wang, J.L. Shi, W.L. Fan, L. Li, Z.B. Zhang,
29 Dissipative particle dynamics study and experimental verification on the pore morphologies and
30 diffusivity of the poly (4-methyl-1-pentene)-diluent system via thermally induced phase
31 separation: The effect of diluent and polymer concentration, Journal of Membrane Science, 514
32 (2016) 487-500.
33 [194] T. Tang, M. Xu, T. Ling, X. Huang, S. Huang, W. Fan, L. Li, Computer simulation and
34 experimental verification of morphology and gas permeability of poly (4-methyl-1-pentene)
35 membranes: effects of polymer chain and diluent extractant, Journal of Materials Science, (2019)
36 1-14.
37 [195] Y. Termonia, Molecular modeling of phase-inversion membranes: effect of additives
38 in the coagulant, Journal of Membrane Science, 104 (1995) 173-180.
39 [196] Y. Termonia, Molecular modeling of structure development upon quenching of a
40 polymer solution, Macromolecules, 30 (1997) 5367-5371.
41 [197] Y. Han, J. Wang, H. Zhang, Effects of kinetics coefficients on ternary phase separation
42 during the wet spinning process, Journal of Applied Polymer Science, 125 (2012) 3630-3637.
43 [198] X.H. He, C.J. Chen, Z.Y. Jiang, Y.L. Su, Computer simulation of formation of
44 polymeric ultrafiltration membrane via immersion precipitation, Journal of Membrane Science,
45 371 (2011) 108-116.

39
 1 [199] S.Y. Liu, Z.Y. Jiang, X.H. He, Computer simulation of the formation of anti-fouling
2 polymeric ultrafiltration membranes with the addition of amphiphilic block copolymers, Journal
3 of Membrane Science, 442 (2013) 97-106.
4 [200] K. Kamide, H. Iijima, H. Shirataki, Thermodynamics of formation of porous polymeric
5 membrane by phase separation method II. particle simulation approach by Monte Carlo method
6 and experimental observations for the process of growth of primary particles to secondary
7 particles, Polymer Journal, 26 (1994) 21-31.
8 [201] U. Dinur, A.T. Hagler, "New approaches to empirical force fields" Reviews in
9 computational chemistry, New York, NY : VCH Publishers, New York, NY, 1990.
10 [202] M.P. Allen，, D.J. Tildesley, Computer Simulation of Liquids, Oxford University
11 Press, New York, 2017.
12 [203] M. Lisal, J.K. Brennan, W.R. Smith, Mesoscale simulation of polymer reaction
13 equilibrium: Combining dissipative particle dynamics with reaction ensemble Monte Carlo. I.
14 Polydispersed polymer systems, Journal of Chemical Physics, 125 (2006).
15 [204] A.V. Berezkin, Y.V. Kudryavtsev, Hybrid approach combining dissipative particle
16 dynamics and finite-difference diffusion model: Simulation of reactive polymer coupling and
17 interfacial polymerization, Journal of Chemical Physics, 139 (2013).
18 [205] M.B. Liu, G.R. Liu, L.W. Zhou, J.Z. Chang, Dissipative particle dynamics (DPD): An
19 overview and recent developments, Archives of Computational Methods in Engineering, (2014).
20 [206] J. Ramos, J. Martinez-Salazar, Computer modeling of the crystallization process of
21 single-chain ethylene/1-hexene copolymers from dilute solutions, Journal of Polymer Science Part
22 B-Polymer Physics, 49 (2011) 421-430.
23 [207] L. Larini, D. Leporini, A manifestation of the Ostwald step rule: Molecular-dynamics
24 simulations and free-energy landscape of the primary nucleation and melting of single-molecule
25 polyethylene in dilute solution, Journal of Chemical Physics, 123 (2005) 11.
26 [208] L. Xu, Z.Y. Fan, H.D. Zhang, H.S. Bu, Simulation of crystallization of entangled
27 polymer chains, Journal of Chemical Physics, 117 (2002) 6331-6335.
28 [209] C.M. Chen, P.G. Higgs, Monte-Carlo simulations of polymer crystallization in dilute
29 solution, Journal of Chemical Physics, 108 (1998) 4305-4314.
30 [210] G.D. Kang, Y.M. Cao, Application and modification of poly(vinylidene fluoride)
31 (PVDF) membranes - A review, Journal of Membrane Science, 463 (2014) 145-165.
32 [211] S. Zhang, L. Wu, Research progress of hydrophilic modification of polyvinylidene
33 fluoride membranes, Chemical Industry and Engineering Progress, 35 (2016) 2480-2487.
34 [212] X. Tan, D. Rodrigue, A review on porous polymeric membrane preparation. Part I:
35 Production techniques with polysulfone and poly (vinylidene fluoride), Polymers, 11 (2019) 1160.
36

40