Viscoelastic Phase Separation Model for Ternary Polymer Solutions
aa r X i v : . [ c ond - m a t . s o f t ] F e b Viscoelastic Phase Separation Model for Ternary Polymer Solutions
Viscoelastic Phase Separation Model for Ternary Polymer Solutions
Kenji Yoshimoto and Takashi Taniguchi a) Department of Chemical Engineering, Kyoto University, Kyoto 615-8510, Japan (Dated: 22 February 2021)
When a polymer solution undergoes viscoelastic phase separation, the polymer-rich phase forms a network-like structure even if it is a minor phase. This unique feature is induced by polymer dynamics, which areconstrained by the temporal entanglement of polymer chains. The fundamental mechanisms of viscoelasticphase separation have already been elucidated by theory and experiments over the past few decades; however,it is not yet well understood how viscoelastic phase separation occurs in multicomponent polymer solutions.Here, we construct a new viscoelastic phase separation model for ternary polymer solutions that consist ofa polymer, solvent, and nonsolvent. Our simulation results reveal that a network-like structure is formed inthe ternary bulk system through a phase separation mechanism similar to that observed in binary polymersolutions. A difference in dynamics is also found in that the solvent, whose affinity to the polymer is similarto that to the nonsolvent, moves freely between the polymer-rich and water-rich phases during phase separa-tion. These findings are considered important for understanding the phase separation mechanism of ternarymixtures often used in the manufacturing of polymeric separation membranes.PACS numbers: 82.60. Lf, 83.80. Rs, 83.60. BcKeywords: Polymer, Ternary Solution, Phase Separation, Viscoelastic, Model, Simulations
I. INTRODUCTION
When a binary liquid mixture undergoes phase separa-tion, the minor component usually forms droplets in thematrix of the major component. However, if the mixtureis composed of solvent and polymer and if the polymerchains are long enough to be entangled, the polymer so-lution may exhibit unique phase separation, termed vis-coelastic phase separation.
For example, in the earlystage of viscoelastic phase separation, the major com-ponent, i.e. , the solvent, forms small droplets.
Inthe late stage, the polymer-rich phase forms a network-like continuous structure even if it is the minor phase. These characteristics are due to a large difference in dy-namic behavior between the solvent and the polymer. In the case of a binary liquid, each component diffusesdown the concentration gradient, with speeds nearly thesame between the two components. In viscoelastic phaseseparation, however, the entangled polymer chains movemuch more slowly than the solvent. In addition, the dy-namic mechanical response differs between the polymerand the solvent. The entangled polymer deforms instan-taneously as an elastic body does but also flows like a liq-uid over a period much longer than the relaxation time.These differences, which are referred to as dynamic asym-metries, may affect the mechanism and morphology ofphase separation through tight coupling with concentra-tion diffusion.
In addition to the binary polymer solution, phase sep-aration of the ternary polymer solution, which consists ofpolymer, solvent, and nonsolvent ( e.g. , water), is partic-ularly important when preparing a polymer membrane.For example, in the nonsolvent-induced phase separation a) Electronic mail: [email protected] (NIPS) process, which is widely used for the manufac-ture of polymer membranes, the polymer is dissolved ina solvent that is soluble both in the polymer and in wa-ter. By immersing the polymer solution in water, thesolvent is partially replaced by the water. The resulting ternary polymer solution is thermodynamically unstabledue to a strong repulsion between the water and the poly-mer and immediately separates into water-rich (= ma-jor) and polymer-rich (= minor) phases. The water-richphase forms small droplets that become pores, while thepolymer-rich phase forms a network-like structure as amatrix. Both phases are essentially similar to those ob-served in the viscoelastic phase separation of a binarypolymer solution.To provide insight into the phase separation mecha-nism in the NIPS process, some phase separation modelshave been proposed for the ternary polymer solution.
Note, however, that the viscoelastic nature of the entan-gled polymer is not explicitly considered in any models.For example, the free energy of the ternary polymer so-lution has been described as a sum of the Flory-Hugginsmixing energy and the interfacial energy, neglecting theelastic energy generated by deformation of the entangledpolymer. Additionally, the relaxation of the entangledpolymer is not considered in the existing models. Indeed,it has been demonstrated with binary polymer solutionsthat both elastic and relaxation behavior of the entangledpolymer plays a crucial role in the formation of the frozenpolymer-rich phase at the early stage and of the network-like polymer-rich structure at the late stage.
In this paper, we develop a basic model for the phaseseparation of a ternary polymer solution that explicitlyincludes the elastic and relaxation characteristics of theentangled polymer. Our approach is based on the the-ory of viscoelastic phase separation developed for binarypolymer solutions.
First, we describe the free en-iscoelastic Phase Separation Model for Ternary Polymer Solutions 2ergy of the ternary polymer solution as a functional of thevolume fractions of each component and the deformationtensor of the entangled polymer chains (Sec. I A and I B).Then, we prepare a set of equations for evolving the vol-ume fractions and the deformation tensor over time (Sec.I C). The coupling between the volume fractions and thedeformation tensor occurs through the polymer velocity,which is determined from a balance among the thermo-dynamic, hydrodynamic, and elastic forces. Three dif-ferent expressions for the polymer velocity are examined:a thermodynamic term only (diffusion case), thermody-namic + hydrodynamic terms (viscous case), and ther-modynamic + hydrodynamic + elastic terms (viscoelas-tic case) (Sec. II). The first two cases are representativesof the previous models, and the last case is the focus ofdevelopment in this study. Numerical calculations areperformed to simulate the spinodal decomposition of abulk ternary polymer solution and the change in mor-phology of a square-shaped polymer-rich phase immersedinto a water-rich bath. By comparing the results ob-tained from the three cases, we highlight some importantaspects of the viscoelasticity of an entangled polymer onthe phase separation of a ternary polymer solution.
A. Field variables: φ and W A ternary mixture of polymer, solvent, and water isdescribed by three volume fractions: φ p , φ s , and φ w .The subscripts p, s, and w denote polymer, solvent, andwater, respectively. Each volume fraction varies with po-sition r and time t , whereas the sum of the three volumefractions is always unity due to incompressibility; X α =p , s , w φ α = 1 . (1)In this paper, subscript α (or β ) is used to denote one ofthe three components: polymer, solvent, and water.Polymer chains are assumed to be long enough to havetemporal entanglements, as schematically illustrated inFig. 1. The entangled polymer chains form a tempo-ral network structure that behaves like an elastic objectfor instantaneous, infinitesimal deformation. In contrast,the network structure flows like liquid under a constantload as a result of the stretch and orientation relaxationsof entangled polymer chains by reptation motion anda constraint release mechanism. Such viscoelastic be-havior of polymer entanglements can be included in themodel by introducing another field variable, the so-calledconformation tensor W . A microscopic definition of W is given as follows. Let us consider an infinitesimal vol-ume element that contains n p ( ≫
1) polymer chains. Asschematically shown on the right-hand side of Fig. 1, allpolymer chains are assumed to be Gaussian composed of N statistical segments with a size of a . The ℓ -th poly-mer chain ( ℓ = 1 , · · · , n p ) has Z ℓ ( ≫
1) entangled seg-ments whose positions are denoted by r ℓ , ..., r ℓZ ℓ . Thevectors connecting two successive entangled segments, !" ! " !" FIG. 1. Schematic representation of the temporal networkstructure of entangled polymer chains. R ℓm (= r ℓm +1 − r ℓm , where m = 1 , ..., Z ℓ − R ℓm is a √ N e , where N e isthe average number of segments existing between the twoadjacent entangled segments along a chain. The normal-ized dyadic tensor R ℓm R ℓm /a N e represents the changein shape of the network structure from equilibrium. Theaverage of R ℓm R ℓm /a N e over all the entangled segmentsand chains existing in the same infinitesimal volume ele-ment is defined as the conformation tensor W ; W ≡ n p " n p X ℓ =1 Z ℓ − Z ℓ − X m =1 R ℓm R ℓm a N e . (2)The right-hand side of Eq. (2) becomes an identity ten-sor I at equilibrium. Eq. (2) provides microscopic in-sight into the origin of polymer viscoelasticity. However,calculation of W from Eq. (2) is computationally de-manding since it requires an explicit description of allpolymer chains and entanglements. Alternatively, onecan treat W as a continuous field variable and evolveit through a constitutive equation. In this approach,microscopic details of the polymer chains are discarded,whereas any local deformation of the network structurecan be described by the same macroscopic constitutiverelation arising from numerous entanglements of polymerchains. As described below, we take the latter approach,treating W as a continuous, frame-invariant deformationtensor and evolving W by a simple viscoelastic constitu-tive equation. B. Thermodynamic properties
Hereafter, we make all variables and equations dimen-sionless to generalize the model. For a ternary systemwith constant temperature T and volume V , the free en-ergy F can be described as a functional of the volumefractions and conformation tensor; F [ φ p , φ s , φ w , W ] = Z V d r ( f mix + f int + f ela ) , (3)where φ w = 1 − φ p − φ s . The first integrand in theright-hand side of Eq. (3), f mix , represents a mixing freeiscoelastic Phase Separation Model for Ternary Polymer Solutions 3energy as f mix = φ p N ln φ p + φ s ln φ s + φ w ln φ w + χ ps φ p φ s + χ sw φ s φ w + χ wp φ w φ p , (4)where χ αβ is the Flory-Huggins parameter for the affinitybetween two different components, α and β ( = α ). Thesecond integrand f int is the interfacial energy given by f int = κ ps | ∇ ( φ p − φ s ) | + κ sw | ∇ ( φ s − φ w ) | + κ wp | ∇ ( φ w − φ p ) | , (5)where κ αβ is an interfacial parameter for a pair of compo-nents α and β ( = α ). The last integrand f ela representselastic energy arising from deformation of the temporalnetwork composed of entangled polymer chains. Amongseveral different formulae proposed for f ela , we select theone that contains only the shear modulus of the polymersolution; f ela = G W − I ) : ( W − I ) , (6)where I is the identity tensor. The shear modulus G isassumed to vary with the polymer concentration basedon the scaling theory; G = G φ . (7)The scaling of Eq. (7) stems from the characteristiclength between the adjacent entangled segments, ξ . Ac-cording to the book of de Gennes , the shear modu-lus of entangled polymer solution can be expressed as G ∼ k B T /ξ . Since ξ scales with 1 /φ p , G is inverselyproportional to φ as shown in Eq. (7). The coefficient G will be used as a parameter to characterize the elas-ticity of the polymer solution.The chemical potential µ α ( α = p, s) and the stresstensor σ are derived from the functional derivative ofthe free energy F in Eq. (3); µ α = δFδφ α = ∂f∂φ α − ∇ · ∂f∂ ( ∇ φ α ) , (8) σ = 2 W · δFδ W = 2 G W · ( W − I ) , (9)where f ≡ f mix + f int + f ela . The detailed expressions of µ α are summarized in Appendix A. C. Time evolution of the field variables
The volume fractions are evolved with the equationsof continuity; ∂φ p ∂t = − ∇ · ( φ p v p ) , ∂φ s ∂t = − ∇ · ( φ s v s ) . (10)Note that φ w is determined by the compressibility con-dition, Eq. (1). In this study, the random noise term arising from the thermal fluctuations is not includedin Eq. (10) so as to clarify the viscoelastic effects onthe domain dynamics at the late stage of phase sepa-ration. In general, the thermal fluctuations have littleinfluence on the domain dynamics after the formation ofsharp interfaces . The velocities of the polymer and sol-vent, v p and v s , respectively, can be obtained by solvingthe equations of motion derived based on the concept of stress division . The resulting expressions are writtenas v p = v − φ p (cid:26) L pp (cid:18) ∇ µ p − ∇ · σ φ p (cid:19) + L ps ∇ µ s (cid:27) , (11) v s = v − φ s (cid:26) L sp (cid:18) ∇ µ s − ∇ · σ φ s (cid:19) + L ss ∇ µ s (cid:27) , (12)where L αβ is the component of the transport coefficientmatrix L and where v is the mean velocity defined as v = v p φ p + v s φ s + v w φ w . (13)The mean velocity v satisfies the Stokes equation;0 = − ∇ p − X i =p , s φ i ∇ µ i + η ∆ v + ∇ · σ , (14)where p and η denote the pressure and viscosity of theternary mixture, respectively. In this study, we simplifyEqs. (11) and (12) by assuming that the diagonal compo-nents of L are constants and the off-diagonal componentsare zero; v p = v − L pp φ p (cid:18) ∇ µ p − ∇ · σ φ p (cid:19) , (15) v s = v − L ss φ s ∇ µ s . (16)Equation (15) has exactly the same expression as thatderived for a polymer-solvent binary system. A new ad-dition to the ternary system is Eq. (16), which allows thesolvent to move separately from the polymer.To evolve the conformation tensor simultaneously withthe volume fractions described above, we employ a consti-tutive equation, the so-called upper-convected Maxwellequation, D W Dt = ( ∇ v p ) · W + W · ( ∇ v p ) T − τ ( W − I ) , (17)where D/Dt represents the substantial derivative(= ∂/∂t + v p · ∇ ) and ( ∇ v p ) T denotes the transpose ofthe velocity gradient tensor. The relaxation time τ is as-sumed to have the same dependency on φ p as the shearmodulus defined in Eq. (7); τ = τ φ , (18)with a constant τ .iscoelastic Phase Separation Model for Ternary Polymer Solutions 4 TABLE I. Simulation parameters.Parameter ValueNumber of statistical segments per chain: N χ ps , χ sw , χ wp φ p0 , φ s0 , φ w0 κ ps , κ sw , κ wp d D t L pp , L ss τ , G , η
10, 100, 0.1 φ s φ w φ p Polymer WaterSolvent
FIG. 2. Phase diagram for the ternary system. The thinand thick lines represent the spinodal and binodal curves,respectively, and they merge on the critical point ( × ). Thesimulation starts from a homogeneously mixed state ( (cid:4) ) andeventually reaches equilibrium composed of polymer-rich ( • )and water-rich ( ◦ ) phases. The final compositions of thesephases can be connected with a straight line (dashed line)passing through the initial composition. x y (a) (b) x φ α FIG. 3. (color online) Polymer-rich domain at equilibrium.(a) Two-dimensional plot of the polymer volume fraction φ p .Color is assigned to each square grid based on the value of φ p . For example, red and blue colors represent polymer-richand water-rich phases, respectively. (b) Profile of the volumefraction along the dashed line in (a). All three volume frac-tions, φ p , φ s , and φ w , are shown here with the symbols ofcircle, square, and diamond, respectively. II. SIMULATIONSA. Parameters
The parameters used in the simulations are summa-rized in Table I. The Flory-Huggins parameters χ αβ are chosen from the literature to represent a ternarymixture composed of polyvinylidene difluoride (PVDF),dimethylformamide (DMF), and water. Note that χ ps and χ sw depend on the concentration. However, forbrevity, they are assumed to be constant. The numberof statistical segments per chain, N , is set at a relativelysmall number, 10, to increase numerical stability and toperform large-scale simulations.Figure 2 illustrates a phase diagram of the ternary sys-tem with N and χ αβ given in Table I. The binodal andspinodal curves are drawn using the mixing free energydensity given in Eq. (4). The initial polymer volumefraction, φ p0 is set to be 0.15, which a typical value usedin the NIPS process . The initial volume fraction ofthe water, φ w0 , is arbitrary chosen to be 0.25, assumingthat ∼
30% of the solvent in the bulk polymer solutionis instantaneously substituted with the water. With thegiven initial composition, the ternary bulk system is al-ready located at the inside of spinodal region so that it isspontaneously decomposed into polymer-rich and water-rich phases. The final compositions of the two phases, φ ′ α, eq and φ ′′ α, eq , are obtained from the fully equilibratedmorphology shown in Fig. 3.The interfacial parameters κ αβ are set arbitrarily to1.0. Then, the grid size d is adjusted to 2.0, with whichthe interfacial regions are expressed by 5 or 6 grid points[see Fig. 3(b)]. For the time-related parameters, thetransport coefficients, L pp and L ss , are set to unity.Then, the time step ∆ t is increased to the maximum,0.025, with which the iterative calculations can be per-formed stably over a long course of simulation time.Finally, the two viscoelastic parameters, τ and G ,are set to 10 and 100, respectively. The viscosity of theternary system, η , is set at 0.1, which is estimated from η = Gτ with φ p = φ p0 . B. Three typical cases
To grasp some fundamental behaviors of phase sepa-ration in the ternary system, we perform simulations forthe following three cases. The same equations and pa-rameters are used in all cases, except that expressions forthe polymer and solvent velocities are modified depend-ing on the cases.(I)
Diffusion case. The velocity of component α (=p,s)is induced only by the gradient of the chemical potentialof the same component; v α = − ( L α /φ α ) ∇ µ α . (19)Neither the mean velocity v nor the stress tensor σ iscalculated here.iscoelastic Phase Separation Model for Ternary Polymer Solutions 5(II) Viscous case. In addition to case (I), the hydrody-namic contribution is taken into account; v α = v − ( L αα /φ α ) ∇ µ α . (20)The Stokes equation, Eq. (14) is used to solve v , but theelastic contribution, ∇ · σ , is excluded.(III) Viscoelastic case. The polymer velocity is expressedby Eq. (15), including all the diffusion, hydrodynamic,and elastic contributions. The solvent velocity is repre-sented by Eq. (16) which is identical to Eq. (20).
III. RESULTS AND DISCUSSIONA. Phase separation in ternary mixtures
Figure 4 illustrates the distributions of the polymervolume fraction φ p obtained at different times. All sim-ulations are started from a homogeneous mixture, asshown on the leftmost image. Uniform random noiseranging between − t = 100 and becomes clearer at t = 250. Thisstage is referred to as the linear (or early) stage of spin-odal decomposition, where the polymer-rich and water-rich phases become denser while maintaining their spa-tial frequencies. Once the polymer-rich phase is bro-ken into pieces through further densification and localshrinkage ( t = 500), smaller pieces of the polymer-richphase are gradually absorbed into larger domains by anevaporation-condensation mechanism ( t > t = 500, where the polymer-rich domains become largerand more rounded than those in case (I). This is becausethe deformation of the polymer-rich domains is acceler-ated by hydrodynamic flow (see Sec. III B for the detailsof the hydrodynamic effects). Once the discretizationand circularization of the polymer-rich domains is settled( t = 1000), the hydrodynamic flow is diminished, andtherefore, the morphology change is significantly sloweddown. Later, the polymer-rich domains are graduallycoarsened through the evaporation-condensation mecha-nism, similar to case (I).In case (III), a bicontinuous structure is clearly seenat t = 2500, which is slower than the other two cases bya factor of 10. It is also noteworthy that the time scaleof this delay is much longer than the relaxation time τ , i.e. , 10. In the viscoelastic phase separation, wheneverthe polymer-rich domain is deformed by the thermody-namic force, the elastic force is instantaneously generatedas a counterforce. Then the residual and new thermo-dynamic force acts again on the polymer-rich domain,generating the new elastic force in the polymer-rich do-main (see Sec. III-B for the detailed mechanisms). Since this cycle occurs iteratively and continuously throughoutthe entire phase separation process, it takes a consider-ably long time for the polymer-rich domain to be fullyrelaxed. Meanwhile, the solvent and water are gradu-ally squeezed out from the polymer-rich phase. As aresult, a network-like structure is formed with irregu-larly shaped polymer-rich domains ( t = 5000). Someframeworks of the network-like structure are maintainedover time, whereas others are merged into larger ones( t > χ wp = 2 . χ sw = 0 .
3) is almost the same as that between thesolvent and polymer ( χ ps = 0 . t = 250 in cases (I) and (II) (resultsare not shown here due to space restriction). Althoughthe solvent’s free movement does not change any basicfeatures of the phase separation, it becomes particularlyimportant for the NIPS process where the solvent needsto move across the interface between the polymer solu-tion and the water .In the rest of the paper, we focus on the analysis ofpolymer-rich phases to clarify some important roles ofpolymer viscoelasticity in phase separation in ternarysystems.The morphology change can be quantified with thecharacteristic wavenumber q defined as q = R d q | q | S ( q ) R d q S ( q ) , (21)where q is the wave vector. The structure factor S ( q ) isobtained from the square of the Fourier transform of thepolymer’s volume fraction, i.e., h| ˆ φ p ( q ) | i , where h ( · · · ) i stands for the statistical average of ( · · · ). Figure 6 shows q of the polymer morphology sampled over a time pe-riod of 2 . × (including the images in Fig. 4). Incase (I), q is almost constant at t < q ∝ t − / . The former cor-responds to the initial stage of spinodal decomposition,and the latter indicates the domain growth driven bythe evaporation-condensation mechanism. In case (II),a sharp decrease in q is observed from t = 300 to t = 800,with a rate of q ∝ t − / . The exponent of − / andfluid mixture where the numerical simulations wereperformed in two dimensions including the hydrodynamiceffects. Note, however, that the hydrodynamic-drivendomain growth generally scales with t − . This dis-iscoelastic Phase Separation Model for Ternary Polymer Solutions 6
FIG. 4. (color online) Phase separation in the bulk ternary system. Each image represents a two-dimensional distribution ofpolymer volume fraction φ p within the simulation box of size L box × L box = 512 . The color is assigned based on the value of φ p , as shown in the scale bar. x y φ p j p x y x y φ s j s x y x y φ w j w x y FIG. 5. (color online) Morphology and flow of each com-ponent at an early stage of phase separation in case (III)( t = 2500): (top) two-dimensional distribution of the volumefraction φ α ( α = p,s,w), and (bottom) volumetric flux j α within the area enclosed by dashed lines in the top image. Acolor is assigned to each grid point based on the value of φ α .Arrows on the bottom images represent a vector field of j α ,whose lengths are magnified by 3 × for clarity. crepancy is possibly due to the difference in dimension,the limitation of system size and time scale, or the com-bination of the hydrodynamic and diffusion effects .In the later time, q remains almost unchanged with thedisappearance of the hydrodynamic flow. At t > , q follows the same trend as case (I), indicating thatthe polymer-rich domains regrow under the evaporation-condensation mechanism. In case (III), q shows a slightdecrease over time, followed by a sharp drop as q ∝ t − / at t > − / t ¯ q −1/3−2/3 Diffusion Case (I)Viscous Case (II)Viscoelastic Case (III)
FIG. 6. (color online) Change in the characteristic wavenum-ber q over time t . unstable parts in the network-like structure.Figure 7 illustrates a time evolution of the circularlyaveraged structure factor S ( q ) ( q : wavenumber) at anearly stage of spinodal decomposition, calculated from S ( q ) = 12 πqdq Z | q | = q + dq | q | = q d q h (cid:12)(cid:12) ˆ φ p ( q ) (cid:12)(cid:12) i . (22)In cases (I) and (II), S ( q ) shows a peak at q max ≈ S max increases exponentially over time,which is one of the typical characteristics in the earlystage of spinodal decomposition. Note that q max canbe estimated numerically by solving the eigenequationsderived in Appendix B. In case (III), no sharp peakis formed, indicating that the phase separation is sup-pressed by the elastic effects. Another feature of case(III) is that S ( q ) decreases over time at q > q max . Thisis due to smoothing of the random noise added to the ini-tial volume fractions. An opposite trend is seen in cases(I) and (II) since the random noise is already diminishedand the interface is sharpened over time.In addition to the morphological change, the phasecompositions also change differently in the three cases.iscoelastic Phase Separation Model for Ternary Polymer Solutions 7 −10 −5 S Diffusion Case (I)t=100t=150t=200t=250 0 100 200 30010 S m a x −10 −5 S Viscous Case (II)t=100t=150t=200t=250 0 100 200 30010 S m a x −10 −5 S q Viscoelastic Case (III)t=100t=150t=200t=250 0 100 200 30010 S m a x t FIG. 7. (color online) Structure factor S at the early stageof spinodal decomposition. (Left) A log-log plot of S as afunction of the wavenumber q . The arrow points q = 20(= q max ), where S reaches the maximum, S max , in cases (I)and (II). (Right) A semilog graph of S max as a function oftime t . Figure 8 illustrates the polymer volume fraction aver-aged over the polymer-rich domains ( φ p > φ p0 ) and thepolymer volume fraction averaged over the water-rich do-mains ( φ p < φ p0 ), denoted by φ ′ p and φ ′′ p , respectively.The interfacial regions are excluded from the averaging.In the early stage of spinodal decomposition, φ ′ p and φ ′′ p of case (I) overlap with those of case (II). In the laterstage, φ ′ p and φ ′′ p in case (I) asymptotically approach theequilibrium values, φ ′ p , eq and φ ′′ p , eq , respectively (see Fig.3 for the values of φ ′ p , eq and φ ′′ p , eq ). On the other hand,in case (II), the phase compositions reach equilibrium al-most ten times faster than in case (I). In case (III), φ ′ p and φ ′′ p are constrained to the initial value φ p0 (=0.15)over a relatively long period, and then they exhibit asharp change similar to those in the early stage of cases(I) and (II). All these trends are consistent with thoseobserved in Figs. 4-7, except that in case (III), φ ′ p settlesinto a lower value (0.35) than φ ′ p , eq (0.48). This can alsobe seen as a color difference in the polymer-rich domainsat t > φ ′ p still keeps increasingat a considerably slow rate; it will take an extremely longtime for φ ′ p to reach to φ ′ p , eq . t φ ′ p , φ ′′ p FIG. 8. (color online) Averaged polymer volume fractions inthe polymer-rich phase ( φ ′ p ) and in the water-rich phase ( φ ′′ p ):circle for diffusion case (I), diamond for viscous case (II), andsquare for viscoelastic case (III). The solid and dashed linescorrespond to φ ′ p and φ ′′ p , respectively. t = 2 . × t = 7 . × t = 12 . × x y t = 0 φ p Viscoelastic Case (III)ViscousCase (II)Diffusion Case (I)
FIG. 9. (color online) Change in shape of the polymer-richphase. For all cases, the polymer-rich phase is set to be asquare at t = 0, as shown on the leftmost image. Each imagerepresents a two-dimensional distribution of φ p within thesimulation box of L box × L box = 64 . Colors are assignedbased on the scale bar. B. Driving forces for morphology change
To clarify a relationship between the morphologychange of the polymer-rich domains and the polymer ve-locity v p , we set up another initial condition, as illus-trated in the leftmost image in Fig. 9. The compositionsof the polymer-rich and water-rich phases, φ ′ α and φ ′′ α ( α = p, s, w), are estimated from the bicontinuous struc-ture formed at t = 2500 in case (III): ( φ ′ p , φ ′ s , φ ′ w ) =(0.194, 0.584, 0.222), and ( φ ′′ p , φ ′′ s , φ ′′ w ) = (0.102, 0.617,0.281). The size of the square-shaped polymer-rich phaseis determined in a way that the averaged volume fractionsover the system become equal to the initial volume frac-tions used in the previous section, φ α . With this setting, φ ′ p and φ ′′ p eventually reach φ ′ p , eq and φ ′′ p , eq , as seen in Fig.8. Figure 9 shows the polymer morphology at t = 0, 250,iscoelastic Phase Separation Model for Ternary Polymer Solutions 8 FIG. 10. (color online) Polymer velocity v p in the polymer-rich phase at t = 2500 in Fig. 9. The velocity is representedby arrows on top of the distribution of φ p . The velocity fieldon the leftmost image corresponds to v p , and the ones onthe second, third, and fourth left columns represent the meanvelocity v (not calculated in case (I)), the polymer velocitydriven by the gradient of chemical potential, v p ,µ , and thatdriven by the divergence of the stress tensor, v p ,σ (calculatedonly in case (III)), respectively. For clarity, the number of ar-rows is reduced from 32 ×
32 to 16 ×
16, and only those locatedinside of the polymer-rich phase ( φ p > .
15) are shown here.All arrows are magnified uniformly by 200; however, some ofthem are still too small to be apparent. t = 250). The polymermovement is plotted in Fig. 10, where the majority ofthe polymer velocity v p points to the polymer-rich cor-ner regions. Note that in case (I), v p is equal to thevelocity driven by the gradient of the polymer poten-tial, v p ,µ (= − L pp /φ p ∇ µ p ), as defined in Eq. (19). Thepolymer-rich phase constantly exudes the solvent and wa-ter until it reaches equilibrium. As a result, the polymer-rich phase is continuously shrunk and transformed intoa circular shape ( t = 1250) through a rounded squareshape ( t = 750). In case (II), the morphology changeoccurs similarly to that in case (I); however, the cornersof the polymer-rich square become rounded earlier thanin case (I). This is mainly due to the hydrodynamic flow,which is not considered in case (I). As defined in Eq.(20), the polymer velocity v p in case (II) is composedof the mean velocity v and the diffusion-driven veloc-ity v p ,µ . For instance, each velocity at t = 250 can beseen in Fig. 10. Similar to case (I), v p ,µ points to thepolymer-rich corners. On the other hand, v exhibits alarge flow that enhances the rounding of the polymer-rich phase. Once the polymer-rich phase is transformedinto a circular shape ( t = 750), the hydrodynamic ef-fect is diminished, and the change in the composition of W xx − x y −0.01 0 0.01 W xy x y −0.01 0 0.01 W yy − x y −0.01 0 0.01 FIG. 11. (color online) Conformation tensor W at t = 2500:(left) W xx −
1, (center) W xy , and (right) W yy −
1. A color isassigned using the scale bar on the bottom side. The negative(or positive) regions in W ii − i = x or y ) indicate wherethe polymer is locally compressed (or stretched) along the i direction. the polymer-rich phase is governed by the same diffusionprocess as in case (I). In case (III), the square shape ofthe polymer-rich phase is maintained even at t = 12500.As expressed in Eq. (15), v p in case (III) is the sum of v , v p ,µ , and the velocity driven by the divergence of thestress tensor, v p ,σ (= L pp /φ ∇ · σ ). Notably, the flowof v p ,σ in Fig. 10 is almost perfectly opposite to that of v p ,µ . This can be interpreted as a balance between theexternally acting thermodynamic force and the internallygenerated elastic force. Indeed, the polymer-rich phaseis under internal compression, as shown in Fig. 11; theelastic force is generated to counteract the compressivestrain. Note that according to the Stokes equation asshown in Eq. (14), v is induced only by solvent move-ment when v p ,µ is completely offset by v p ,σ . Therefore, v in case (III) becomes much smaller than in case (II),resulting in a large reduction in any polymer movements.The magnitude of the polymer velocity, | v p | , as shownin Fig. 12, is calculated by taking an average of the mag-nitudes of the polymer velocity within the polymer-richphase ( φ p > . | v p | in the early stageshows a gradual decrease over time (100 ≤ t ≤ ≤ t ≤ | v p | ( t ≥ t = 100 to t = 200 and anotherfrom t = 400 to t = 600. The first plateau is associatedwith the rounding of the polymer-rich phase driven byhydrodynamic flow. The second plateau overlaps with | v p | in case (I), indicating that the morphology changeof the polymer-rich phase is driven by the diffusion pro-cess, similar to case (I). In case (III), the initial value of | v p | is approximately 5 times smaller than that of cases(I) and (II) due to the elastic effects mentioned above.Since the elastic force is constantly generated to balanceout the thermodynamic force, | v p | remains the same or-iscoelastic Phase Separation Model for Ternary Polymer Solutions 9 −8 −7 −6 −5 −4 −3 −2 −1 t | v p | Diffusion Case (I)Viscous Case (II)Viscoelastic Case (III)
FIG. 12. (color online) Change in the magnitude of the poly-mer velocity, | v p | , over time t . The arrows on the top x -axisindicate where the data of Fig. 10 are sampled. der of magnitude over the entire simulation time. Thisexplains why the morphology of the polymer-rich phaseremains square, keeping the polymer volume fraction rel-atively low inside. IV. CONCLUSION
We constructed a new model for the viscoelastic phaseseparation of ternary polymer solutions consisting of en-tangled polymers, solvents, and nonsolvents. The effectsof both elasticity and relaxation of the entangled poly-mer were incorporated into the model by adding elasticenergy to the free energy of the bulk ternary system andby introducing relaxation time into the constitutive equa-tion. Our numerical simulations of the bulk ternary poly-mer solution demonstrated that the polymer-rich phasewas frozen at the early stage of viscoelastic phase sep-aration, whereas it was transformed into a network-likestructure and gradually relaxed over a long time at thelate stage. Both phenomena were essentially the same asthose observed in the binary polymer solution. A newfinding was that the solvent, which was soluble to thepolymer and the nonsolvent, moved across the polymer-rich and water-rich phases during the viscoelastic phaseseparation of the polymer ternary solution. As a finalnote, we simplified our model by setting some modelparameters ( e.g. , mobility coefficients L αβ and the vis-cosity η ) to be constant to clarify the viscoelastic ef-fects of the entangled polymer on ternary phase sepa-ration. In reality, these parameters vary with the con-centration and affect the phase separation to some de-gree. Indeed, some previous studies showed in simula-tions that the concentration dependency of L αβ and η could be crucial for the formation of the glassy polymer-rich domains and asymmetric polymeric membranesformed in the NIPS process . It should be em-phasized here that both viscous and elastic features ofthe polymer solution are fundamentally essential for theviscoelastic phase separation; a simple addition of the concentration-dependent parameters into the diffusion orviscous model may not be sufficient to reproduce thefrozen and recovering features of the viscoelastic phaseseparation simultaneously. We are currently investigat-ing the combinatory effect of the viscoelasticity of theentangled polymer, the concentration-dependent modelparameters, and the thermal fluctuations on the phaseseparation of ternary polymer solutions. ACKNOWLEDGMENTS
This work was partially supported by JSPS KAKENHIGrant Number 19H01862. We also acknowledge the su-percomputer systems of the Institute for Chemical Re-search and Academic Center for Computing and MediaStudies, Kyoto University.
Data Availability
The data that supports the findings of this study areavailable within the article. M. Doi and S. F. Edwards,
The Theory of Polymer Dynamics (Oxford University Press, 1986). H. Tanaka, J. Phys. Condens. Matter , 207 (2000). T. Taniguchi, J. Soc. Rheol. Jpn. , 27 (2004). H. Tanaka, Phys. Rev. Lett. , 787 (1996). A. Onuki and T. Taniguchi, J. Chem. Phys. , 5761 (1997). T. Taniguchi and A. Onuki, Phys. Rev. Lett. , 4910 (1996). B. Zhou and A. C. Powell, J. Membr. Sci. , 150 (2006). Y. Mino, T. Ishigami, Y. Kagawa, and H. Matsuyama, J. Membr.Sci , 104 (2015). D. R. Tree, K. T. Delaney, H. D. Ceniceros, T. Iwama, and G. H.Fredrickson, Soft Matter , 3013 (2017). H. Tanaka, Phys. Rev. E , 4451 (1997). T. Araki and H. Tanaka, Macromolecules , 1953 (2001). R. G. Larson,
Constitutive Equations for Polymer Melts and So-lutions (Butterworth-Heinemann, 2013). P. J. Flory, J. Chem. Phys , 51 (1942). A. Onuki, J. Phys. Soc. Jpn. , 3423 (1990). P.-G. de Gennes,
Scaling Concepts in Polymer Physics (CornellUniv. Press, 1979). S. T. Milner, Phys. Rev. E , 3674 (1993). K. Fukawatase, K. Yoshimoto, and M. Ohshima, Jpn. J. Appl.Phys. , 06FE01 (2015). G. Brown and A. Chakrabarti, J. Chem. Phys. , 2451 (1993). A. Onuki, J. Non-Cryst. Solids , 1151 (1994). H. Matsuyama, M. Teramoto, R. Nakatani, and T. Maki, J.Appl. Polym. Sci. , 159 (1999). D. R. Tree, T. Iwama, K. T. Delaney, J. Lee, and G. H. Fredrick-son, ACS Macro Lett. , 582 (2018). J. U. Garcia, T. Iwama, E. Y. Chan, D. R. Tree, K. T. Delaney,and G. H. Fredrickson, ACS Macro Lett. , 1617 (2020). H. Tanaka and T. Araki, Chem. Eng. Sci. , 2108 (2006). H. Furukawa, Phys. Rev. E , 1150 (1997). E. D. Siggia, Phys. Rev. A , 595 (1979). T. Koga, K. Kawasaki, M. Takenaka, and T. Hashimoto, PhysicaA , 473 (1993). E. Scholten, L. M. C. Sagis, and E. van der Linden, Macro-molecules , 3515 (2005). B. F. Barton, P. D. Graham, and A. J. McHugh, Macromolecules , 1672 (1998). Y. Tanaka, I. Noda, and M. Nagawasa, Macromolecules , 2220(1985). Y. Tanaka, M. Umeda, and I. Noda, Macromolecules , 2257(1988). iscoelastic Phase Separation Model for Ternary Polymer Solutions 10 Appendix A: Chemical potential
The chemical potentials µ p and µ s are expressed as µ p = µ ∗ p − µ ∗ w , µ s = µ ∗ s − µ ∗ w , (A1)where µ ∗ p = 1 N (ln φ p + 1) + χ ps φ s + χ wp φ w − κ ps (∆ φ p − ∆ φ s ) − κ wp (∆ φ p − ∆ φ w )+ 34 G φ p2 ( W − I ) : ( W − I ) , (A2) µ ∗ s = ln φ s + 1 + χ ps φ p + χ sw φ w − κ ps (∆ φ s − ∆ φ p ) − κ sw (∆ φ s − ∆ φ w ) , (A3) µ ∗ w = ln φ w + 1 + χ sw φ s + χ wp φ p − κ sw (∆ φ w − ∆ φ s ) − κ wp (∆ φ w − ∆ φ p ) . (A4) Appendix B: Analysis of the characteristic wavelength internary spinodal decomposition
Here, we assume that ternary phase separation occursby a simple diffusion mechanism, as in case (I). The con-tinuity equation for the polymer is written as ∂φ p ∂t = L pp ∇ · (cid:20) ∇ δFδφ p (cid:21) = L pp ∆ ( h p + K p1 ∆ φ p + K p2 ∆ φ s ) , (B1)where h p ≈ ∂f mix ∂φ p (cid:12)(cid:12)(cid:12)(cid:12) φ α + ∂ f mix ∂φ (cid:12)(cid:12)(cid:12)(cid:12) φ α ( φ p − φ p0 ) , (B2) K pp = − κ ps − κ sw + 4 κ wp , (B3) K ps = κ − κ sw + 2 κ wp . (B4)At the early stage, each volume fraction deviates slightlyfrom the initial values; φ α ( r , t ) = φ α + u α ( r , t ) , (B5)where u α ( r , t ) ( α =p, s) denotes a small perturbation.By substituting Eq. (B5) into Eq. (B1), the polymercontinuity equation is expressed with respect to u p ; ∂u p ∂t = L pp ∆ ∂ f mix ∂φ (cid:12)(cid:12)(cid:12)(cid:12) φ α u p + X α =p , s K p α ∆ u α ! . (B6) Similarly, the continuity equation for the solvent can bedescribed as ∂u s ∂t = L ss ∆ ∂ f mix ∂φ (cid:12)(cid:12)(cid:12)(cid:12) φ α u s + X α =p , s K s α ∆ u α ! , (B7)where K sp = κ − κ + 2 κ , (B8) K ss = − κ − κ + κ . (B9) −0.03−0.02−0.010.000.010.02 k − λ k (a) −1.2−1.0−0.8−0.6−0.4−0.20.00.2 k − λ k (b) FIG. 13. Frequency-dependent coefficient − λk : (a) − λ k and (b) − λ k . The circular mark on (a) denotes the maxi-mum point, (0.248, 0.0158). After taking the Fourier transform, Eqs. (B6) and (B7)can be written in matrix form as ddt (cid:20) ˆ u p ˆ u s (cid:21) = − k Ω (cid:20) ˆ u p ˆ u s (cid:21) (B10)where k is the angular frequency (=2 πq/L box ) and ˆ u α isthe Fourier transform of u α . The matrix Ω is defined as Ω = L pp (cid:16) ∂ f mix ∂φ (cid:12)(cid:12)(cid:12) φ α − K pp k (cid:17) − L pp K ps k L ss (cid:16) ∂ f mix ∂φ (cid:12)(cid:12)(cid:12) φ α − K sp k (cid:17) − L ss K ss k . (B11)The eigenvalues of Ω , λ , can be obtained by solvingdet | Ω − λ I | = 0 , (B12)which can be negative or positive depending on the valueof k . In the eigenspace, Ω on the right-hand side of Eq.(B10) can be replaced with a diagonal matrix composedof λ . Therefore, the sign of − λk is essential to determinewhether fluctuation of the volume fractions at a givenfrequency k is amplified or diminished over time. Thetwo eigenvalues calculated from Eq. (B12) are plotted inFig. 13. Here, − λ k forms a positive peak at k = 0 . ≡ k max ), whereas − λ k always becomes negative at k >
0. This indicates that the fluctuation of the volumefraction at k = k max may be amplified at the fastest rate.Note that the wavenumber corresponding to k max is 20.2;the result here agrees with q maxmax