Thermodynamic Modeling and Validation of the Temperature Influence in Ternary Phase Polymer Systems

The effect of the temperature, as a process variable in the fabrication of polymeric membranes by the non-solvent induced phase separation (NIPS) technique, has been scarcely studied. In the present work, we studied the influence of temperature, working at 293, 313 and 333 K, on the experimental binodal curves of four ternary systems composed of PVDF and PES as the polymers, DMAc and NMP as the solvents and water as the non-solvent. The increase of the temperature caused an increase on the solubility gap of the ternary system, as expected. The shift of the binodal curve with the temperature was more evident in PVDF systems than in PES systems indicating the influence of the rubbery or glassy state of the polymer on the thermodynamics of phase separation. As a novelty, the present work has introduced the temperature influence on the Flory–Huggins model to fit the experimental cloud points. Binary interaction parameters were calculated as a function of the temperature: (i) non-solvent/solvent (g12) expressions with UNIFAC-Dortmund methodology and (ii) non-solvent/polymer (χ13) and solvent/polymer (χ23) using Hansen solubility parameters. Additionally, the effect of the ternary interaction term was not negligible in the model. Estimated ternary interaction parameters (χ123) presented a linear relation with temperature and negative values, indicating that the solubility of the polymers in mixtures of solvent/non-solvent was higher than expected for single binary interaction. Finally, PES ternary systems exhibited higher influence of the ternary interaction parameter than PVDF systems.


Introduction
Non-solvent induced phase inversion (NIPS) is a synthesis technique commonly used to manufacture asymmetric polymeric membranes for separation processes [1], both in laboratory and industrial scale [2,3]. The mechanism of membrane formation using NIPS consists in the precipitation of a polymer solution when it is introduced in a non-solvent bath by the exchange between solvent and non-solvent or coagulant [4]. The liquidliquid demixing that occurs in the NIPS process, and ultimately the membrane structure, are strongly affected by the thermodynamics and kinetics of the system. Ternary phase diagrams, which are temperature dependent, are employed to depict the thermodynamics and kinetic characteristics of a certain NIPS process. Figure 1 shows a ternary diagram in which the continuous and dashed blue curves represent the thermodynamic binodal and spinodal curves, respectively. The area between both curves is called the metastable region. The point where these curves meet is the critical point. The binodal curve defines the boundary between a homogeneous phase and a heterogeneous phase. When the precipitation pathway (yellow line) crosses the binodal curve (point B) the polymer nucleation takes place and the polymer solution separates into two equilibrium phases, polymer-rich (B ) and polymer-lean (B") phases. These equilibrium compositions are located on the binodal curve and depict a tie-line in the ternary phase diagram. The pathway rium phases, polymer-rich (B′) and polymer-lean (B″) phases. These equilibrium compositions are located on the binodal curve and depict a tie-line in the ternary phase diagram. The pathway between points B-C defines the precipitation of the entire membrane, and finally, from C to D points the rest of the solvent is exchanged with the non-solvent.  [5,6].
The final membrane morphology is described by the distribution and shape of the pores created in the NIPS process. The position of the binodal curve and the rate of solvent/non-solvent exchange modify the pore shape. When the binodal curve is located near the polymer-solvent axis, the precipitation pathway reaches more rapidly the binodal curve so big finger-like pores are expected. Otherwise, when the binodal curve is displaced to the right side of the diagram, smaller sponge-like pores are obtained. The entry point (B) of the precipitation pathway into the heterogeneous phase affects the pore distribution. Depending on the pathway followed during solvent/non-solvent exchange different solid phase morphologies can occur (Figure 1b), i.e., dense membranes (pathway 1) were obtained when the pathway crosses the gel homogeneous region, at sufficiently high initial polymer concentration causing slow exchange due to high viscosity; binodal decomposition occurs if the precipitation pathway enters the heterogeneous phase through the metastable region, to form (pathway 2) open-cellular or closed-cellular membrane morphologies, or (pathway 4) nodular pore distributions; otherwise, spinodal decomposition takes place if the precipitation pathway enters the heterogeneous phase through the heterogeneous phase, obtaining (pathway 3) bicontinuous (lacy) distribution of pores [6]. It is deemed that, due to the significant dependence of membrane morphology on thermodynamics and kinetics of the ternary phase systems, the development of mathematical models able to predict the thermodynamic and kinetic characteristics of different polymer/solvent/non-solvent systems would facilitate the decision on the selection of membrane synthesis variables, such as initial polymer concentration, process temperature or even solvent/non-solvent combination.
Different mathematical models have been developed to describe the thermodynamics of ternary systems, such as Flory-Huggins (FH) and compressible regular solution theory (CRS) [7,8]. FH theory, characterized by its relative mathematical simplicity and good prediction of phase behavior [9], is based on the Gibbs free energy of a mixture and the use of binary interaction parameters. One of the limitations of the FH model is that the binary interaction parameters are traditionally obtained in experimental assays. Many efforts have been made to obtain these parameters from tabulated solubility parameters. The final membrane morphology is described by the distribution and shape of the pores created in the NIPS process. The position of the binodal curve and the rate of solvent/non-solvent exchange modify the pore shape. When the binodal curve is located near the polymer-solvent axis, the precipitation pathway reaches more rapidly the binodal curve so big finger-like pores are expected. Otherwise, when the binodal curve is displaced to the right side of the diagram, smaller sponge-like pores are obtained. The entry point (B) of the precipitation pathway into the heterogeneous phase affects the pore distribution. Depending on the pathway followed during solvent/non-solvent exchange different solid phase morphologies can occur (Figure 1b), i.e., dense membranes (pathway 1) were obtained when the pathway crosses the gel homogeneous region, at sufficiently high initial polymer concentration causing slow exchange due to high viscosity; binodal decomposition occurs if the precipitation pathway enters the heterogeneous phase through the metastable region, to form (pathway 2) open-cellular or closed-cellular membrane morphologies, or (pathway 4) nodular pore distributions; otherwise, spinodal decomposition takes place if the precipitation pathway enters the heterogeneous phase through the heterogeneous phase, obtaining (pathway 3) bicontinuous (lacy) distribution of pores [6]. It is deemed that, due to the significant dependence of membrane morphology on thermodynamics and kinetics of the ternary phase systems, the development of mathematical models able to predict the thermodynamic and kinetic characteristics of different polymer/solvent/nonsolvent systems would facilitate the decision on the selection of membrane synthesis variables, such as initial polymer concentration, process temperature or even solvent/nonsolvent combination.
Different mathematical models have been developed to describe the thermodynamics of ternary systems, such as Flory-Huggins (FH) and compressible regular solution theory (CRS) [7,8]. FH theory, characterized by its relative mathematical simplicity and good prediction of phase behavior [9], is based on the Gibbs free energy of a mixture and the use of binary interaction parameters. One of the limitations of the FH model is that the binary interaction parameters are traditionally obtained in experimental assays. Many efforts have been made to obtain these parameters from tabulated solubility parameters. Hansen's solubility parameters are commonly used because of the huge amount of data available for polymers and solvents. However, solubility parameters for new materials (i.e., green solvents or copolymers) are seldom accessible and need to be calculated with group contribution methods such as Hoftyzer and Van Krevelen [10], Hoy [10], Just [11] or Stefanis [12]. Nevertheless, solubility parameters calculated with this methodology do not give accurate predictions yet.
Polyethersulfone (PES) and polyvinylidene fluoride (PVDF) are polymers classically used to manufacture membranes for different industrial separation applications. Particularly PVDF is a polymer with interest to be used for the development of photocatalytic membranes [4]. Temperature plays an important role in the manufacture of polymeric membranes. From the thermodynamic viewpoint, increasing the temperature of the ternary system will increase the demixing gap [2]. Meanwhile, temperature varies the diffusion parameters of the solvent and non-solvent and therefore, the kinetics of exchange between them. Thus changes in the processing temperature affect the final membrane morphology [1]. Despite the fact that during the industrial synthesis of polymeric membranes the temperature of the polymer solution or the coagulation bath is often set different than room temperature, few works have addressed the experimental analysis of the effect of temperature on the thermodynamics of ternary systems [6,13]. Furthermore, there is a lack of validations of predictive models considering temperature for PVDF and few works address the temperature effect for the synthesis of PES membranes. Kahrs et al., evaluated the influence of temperature on the thermodynamics of PES ternary systems [6,13]. When they compared PES membranes fabricated with N-methyl pyrrolidone (NMP) and 2-pyrrolidone (2Pyr) as the solvents, the SEM images showed an increase in the number of pores when the temperature was increased from 20 to 40 • C, because the reduction of the viscosity favored the nuclei formation in the lean-polymer phase [6].
The aim of this work is on the one hand, to develop an enhanced thermodynamic model based on Flory-Huggins theory that incorporates the effect of temperature to predict the binodal curves of different ternary phase diagrams. On the other hand, the model was validated experimentally, considering four ternary systems, using PVDF and PES as the polymers; N,N-dimethylacetamide (DMAc) and N-methyl pyrrolidone (NMP) as the solvents; and water as the non-solvent/coagulant, at three temperatures in the range 293-313 K. The importance of the binary and ternary interaction parameters on the model predictions was explored and discussed. Additionally, a potential influence of the solid polymer state (rubbery or glassy) on behavior of the thermodynamics of the ternary systems at different temperatures was observed.

Thermodynamic Model
In this study, Flory-Huggins theory was selected to describe the thermodynamic model of a ternary system. FH theory is based on the Gibbs free energy of a mixture (∆G M ), Equation (1). The last term, that represents the ternary interactions of the components, is usually neglected to simplify the calculations [14]. In this work, the effect of neglecting this term will also be addressed.
The binary interaction parameters (g ij or χ ij ) represent the interaction between each pair of components. Solvent-polymer (χ 23 ) and non-solvent-polymer (χ 13 ) are interaction parameters independent of the concentration. The non-solvent-solvent interaction parameter (g 12 ) is a solvent concentration-dependent parameter, function of u 2 , which is defined in Equation (7). The ternary interaction parameter (χ 123 ) was estimated through fitting modeled curves to the experimental cloud points [15].
The equilibrium phases, B and B", connected by the tie-line have the same chemical potential for each component (µ i ), Equation (3), but differ in composition. The chemical potential of each component, Equations (4)-(6), is derived from Equation (1).
In the above equations, s and r are the molar volume relations V 1 /V 2 and V 1 /V 3 , respectively. It has been reported that at sufficiently high molecular weight of the polymer the phase diagram is fairly insensitive to the choice of V 3 (polymer molar volume) provided that V 1 /V 3 is sufficiently small [16]. The last term of Equations (4)-(6) corresponds to the partial derivative of the term for ternary interaction. Most of the literature neglects ternary interaction among the components. In this work, the ternary interaction term will be evaluated under two hypotheses: (1) omitting it or (2) considering it.

Influence of the Temperature in the Model
Temperature is an important variable of the phase inversion process that affects the chemical potential, the solubility parameters (δ), the binary interaction parameters, and the components density (ρ) and thus their molar volume. The next diagram ( Figure 2) schematizes the effect of temperature over the mentioned functions and parameters. parameters independent of the concentration. The non-solvent-solvent interaction parameter (g12) is a solvent concentration-dependent parameter, function of u2, which is defined in Equation (7). The ternary interaction parameter (χ123) was estimated through fitting modeled curves to the experimental cloud points [15]. The equilibrium phases, B′ and B″, connected by the tie-line have the same chemical potential for each component (μi), Equation (3), but differ in composition. The chemical potential of each component, Equations (4)-(6), is derived from Equation (1).

Influence of the Temperature in the Model
Temperature is an important variable of the phase inversion process that affects the chemical potential, the solubility parameters (δ), the binary interaction parameters, and the components density (ρ) and thus their molar volume. The next diagram ( Figure 2) schematizes the effect of temperature over the mentioned functions and parameters. The temperature is directly introduced in the expression of chemical potentials. The molar volume (Vi) changes with the temperature because of the temperature dependence of the density (ρ), Equation (8). Density values of the solvents and non-solvent at different temperatures were retrieved from the database of the software Aspen Plus V9. Table 1 compiles the relevant data and properties of the solvents and non-solvent for the model in this study. The temperature is directly introduced in the expression of chemical potentials. The molar volume (V i ) changes with the temperature because of the temperature dependence of the density (ρ), Equation (8). Density values of the solvents and non-solvent at different temperatures were retrieved from the database of the software Aspen Plus V 9 . Table 1 compiles the relevant data and properties of the solvents and non-solvent for the model in this study.  Over the years, the authors have deepened in the concentration dependence of the binary interaction parameters, which initially were considered constant. Later, Altena and Smolders [17] determined the concentration dependency of g12, and Yilmaz et al. [18] analyzed the influence of concentration on χ23, although both studies agreed to consider χ13 as constant parameter. Currently, it is widely agreed that only g12 is a concentration-dependent parameter [15,[19][20][21][22].
Interaction parameters χ13 and χ23 can be obtained experimentally; light scattering, vapor pressure depression and membrane osmometry are techniques commonly proposed to determine χ23, and equilibrium swelling for χ13. However, reliable estimations can be obtained using Equation (9), defined by Hansen [23] and widely used in the literature [15,[24][25][26][27][28][29]. This expression uses Hansen solubility parameters that are composed of three contributions: dispersive (δd), polar (δp) and hydrogen bonding (δh). Equation (9) also includes the correction factor (αij).  The dependency of solubility parameters with the temperature can be calculated with Equations (10)-(12) [23]. Molar volumes (V and Vref, at the reference temperature of 298 K) for solvents and non-solvent can be found in Table 1. For polymers, the molar volumes in Equations (10)- (12) refer to the molar volume of the monomer (V* and V*ref) calculated using Equations (13) and (14)   Over the years, the authors have deepened in the concentration dependence of the binary interaction parameters, which initially were considered constant. Later, Altena and Smolders [17] determined the concentration dependency of g12, and Yilmaz et al. [18] analyzed the influence of concentration on χ23, although both studies agreed to consider χ13 as constant parameter. Currently, it is widely agreed that only g12 is a concentration-dependent parameter [15,[19][20][21][22].
Interaction parameters χ13 and χ23 can be obtained experimentally; light scattering, vapor pressure depression and membrane osmometry are techniques commonly proposed to determine χ23, and equilibrium swelling for χ13. However, reliable estimations can be obtained using Equation (9), defined by Hansen [23] and widely used in the literature [15,[24][25][26][27][28][29]. This expression uses Hansen solubility parameters that are composed of three contributions: dispersive (δd), polar (δp) and hydrogen bonding (δh). Equation (9) also includes the correction factor (αij).  The dependency of solubility parameters with the temperature can be calculated with Equations (10)-(12) [23]. Molar volumes (V and Vref, at the reference temperature of 298 K) for solvents and non-solvent can be found in Table 1. For polymers, the molar volumes in Equations (10)- (12) refer to the molar volume of the monomer (V* and V*ref) calculated using Equations (13) and (14) [10]. Equation (13)   Over the years, the authors have deepened in the concentration dependence of the binary interaction parameters, which initially were considered constant. Later, Altena and Smolders [17] determined the concentration dependency of g12, and Yilmaz et al. [18] analyzed the influence of concentration on χ23, although both studies agreed to consider χ13 as constant parameter. Currently, it is widely agreed that only g12 is a concentration-dependent parameter [15,[19][20][21][22].
Interaction parameters χ13 and χ23 can be obtained experimentally; light scattering, vapor pressure depression and membrane osmometry are techniques commonly proposed to determine χ23, and equilibrium swelling for χ13. However, reliable estimations can be obtained using Equation (9), defined by Hansen [23] and widely used in the literature [15,[24][25][26][27][28][29]. This expression uses Hansen solubility parameters that are composed of three contributions: dispersive (δd), polar (δp) and hydrogen bonding (δh). Equation (9) also includes the correction factor (αij).  The dependency of solubility parameters with the temperature can be calculated with Equations (10)-(12) [23]. Molar volumes (V and Vref, at the reference temperature of 298 K) for solvents and non-solvent can be found in Table 1. For polymers, the molar volumes in Equations (10)- (12) refer to the molar volume of the monomer (V* and V*ref) calculated using Equations (13) and (14)  Over the years, the authors have deepened in the concentration dependence of the binary interaction parameters, which initially were considered constant. Later, Altena and Smolders [17] determined the concentration dependency of g 12 , and Yilmaz et al. [18] analyzed the influence of concentration on χ 23 , although both studies agreed to consider χ 13 as constant parameter. Currently, it is widely agreed that only g 12 is a concentrationdependent parameter [15,[19][20][21][22].
Interaction parameters χ 13 and χ 23 can be obtained experimentally; light scattering, vapor pressure depression and membrane osmometry are techniques commonly proposed to determine χ 23 , and equilibrium swelling for χ 13 . However, reliable estimations can be obtained using Equation (9), defined by Hansen [23] and widely used in the literature [15,[24][25][26][27][28][29]. This expression uses Hansen solubility parameters that are composed of three contributions: dispersive (δ d ), polar (δ p ) and hydrogen bonding (δ h ). Equation (9) also includes the correction factor (α ij ).  The dependency of solubility parameters with the temperature can be calculated with Equations (10)-(12) [23]. Molar volumes (V and V ref , at the reference temperature of 298 K) for solvents and non-solvent can be found in Table 1. For polymers, the molar volumes in Equations (10)- (12) refer to the molar volume of the monomer (V* and V* ref ) calculated using Equations (13) and (14) [10]. Equation (13) is employed for the rubbery polymer and Equation (14) for glassy polymers. These equations are function of the van der Waals volume (V w ), the temperature and, in the case of glassy polymers, the glass transition temperature (T g ).  Table 3 compiles the chemical structure and parameters relevant for the mathematical model of the monomer unit of the polymers. Table 3. Chemical structure and properties of monomer units of the polymers.

Polymer
Formula Equation (14) for glassy polymers. These equations are function of the van der Waals volume (Vw), the temperature and, in the case of glassy polymers, the glass transition temperature (Tg). Table 3 compiles the chemical structure and parameters relevant for the mathematical model of the monomer unit of the polymers. Equation (18) presents the effect of concentration on the binary interaction parameter between solvent and non-solvent. It is obtained from the relationship, Equation (15), between the excess of Gibbs free energy (G E ), Equation (16), and the Gibbs free energy for a binary mixture (ΔGM), Equation (17) [30].
∆G M RT = x 1 ln ϕ 1 +x 2 ln ϕ 2 +g 12 x 1 ϕ 2 (17) G E is calculated from the activity coefficient (γ) of vapor-liquid equilibrium data, obtained with the modified UNIFAC-Dortmund methodology based on groups contribution ume (Vw), the temperature and, in the case of glassy polymers, the glass transition temperature (Tg). Table 3 compiles the chemical structure and parameters relevant for the mathematical model of the monomer unit of the polymers. Equation (18) presents the effect of concentration on the binary interaction parameter between solvent and non-solvent. It is obtained from the relationship, Equation (15), between the excess of Gibbs free energy (G E ), Equation (16), and the Gibbs free energy for a binary mixture (ΔGM), Equation (17) [30].
∆G M RT = x 1 ln ϕ 1 +x 2 ln ϕ 2 +g 12 x 1 ϕ 2 (17) G E is calculated from the activity coefficient (γ) of vapor-liquid equilibrium data, obtained with the modified UNIFAC-Dortmund methodology based on groups contribution Equation (18) presents the effect of concentration on the binary interaction parameter between solvent and non-solvent. It is obtained from the relationship, Equation (15), between the excess of Gibbs free energy (G E ), Equation (16), and the Gibbs free energy for a binary mixture (∆G M ), Equation (17) [30].
G E is calculated from the activity coefficient (γ) of vapor-liquid equilibrium data, obtained with the modified UNIFAC-Dortmund methodology based on groups contribution of the activity coefficients [31]. This methodology is described in Appendix A. A polynomial expression suggested by Tompa [32] was used to fit the g 12 as a function of the solvent concentration, Equation (19).

Calculation Procedure
In this model there are six variables corresponding to the volume fraction of the three components in the two phases. The set of equations is formed by the three equilibrium relationships for each component (Equation (3)) and the two mass balance equations in each phase, Equation (20).
The system of equations is solved using KNITRO solver in GAMS Development Corporation 27.3.0, USA by minimizing the objective function F, Equation (21), as the equilibrium is found when the chemical potential of each component is the same in both phases. The polymer lean-phase molar volume fraction (φ 3 ) is selected as an independent variable.
As a summary, Table 4 collects the physical parameters of the mathematical model that have been calculated or estimated by fitting to the experimental cloud points.

Cloud Point
The binodal curve is experimentally obtained by cloud point titration. Cloud point experiments were performed for both solvents (DMAc and NMP) using water as the non-solvent and PVDF and PES as the polymers. First, polymers were dried in an oven overnight at 333 K. Then, polymer-solvent binary mixtures were prepared at several polymer concentrations. Polymeric solutions, of an initial volume of 60 mL, were kept under constant magnetic stirring. Drops of water were added using a digital burette (Tittrete ® , Brand GMBH + CO KG, Wertheim, Germany) to a polymer solution until permanent turbidity was obtained for 1 h. Experiments were performed in duplicate at different polymer concentrations in the range between 2 and 15 wt% using DMAc and NMP as solvents. Higher polymer concentrations were not tested as the high viscosity of the solution hindered its adequate stirring. The influence of the temperature was studied at 293, 313 and 333 K. Student's t-test was used to statistically compare the significance of the difference between the experimental cloud points of each ternary system at the three temperatures.

Experimental Cloud Point
Cloud point results are represented in Figure 3 and Table A3 of Appendix A compiles the values. Although the experimental points of the four systems studied appeared very proximate at 293 K the position of the cloud points were significantly different between solvents (i.e., PVDF/DMAc/water vs. PVDF/NMP/water and PES/DMAc/water vs. PES/NMP/water) and between polymers (i.e., PVDF/DMAc/water vs. PES/DMAc/water and PVDF/NMP/water vs. PES/NMP/water). Cloud point experimental data obtained between 293 and 295 K (room temperature) are widely reported in the literature for the four systems studied [7,14,28,[33][34][35][36][37][38][39][40][41][42][43][44][45][46]. Experimental results of this study fall within the range of results previously reported at room temperature. However, it is remarkable the huge dispersion of data reported for PVDF systems, especially in the PVDF/NMP/water system [41,44]. A relationship between the different polymer molecular weight characteristics of the polymers used in these studies with the variation of the cloud points reported was not found. On the other hand, this could be more rationally attributed to the difficulty observed to determine the experimental change in turbidity for PVDF systems compared to the clear change in PES systems, for which literature cloud points are more consistent [35,36].

Experimental Cloud Point
Cloud point results are represented in Figure 3 and Table A3 of Appendix B compiles the values. Although the experimental points of the four systems studied appeared very proximate at 293 K the position of the cloud points were significantly different between solvents (i.e., PVDF/DMAc/water vs. PVDF/NMP/water and PES/DMAc/water vs. PES/NMP/water) and between polymers (i.e., PVDF/DMAc/water vs. PES/DMAc/water and PVDF/NMP/water vs. PES/NMP/water). Cloud point experimental data obtained between 293 and 295 K (room temperature) are widely reported in the literature for the four systems studied [7,14,28,[33][34][35][36][37][38][39][40][41][42][43][44][45][46]. Experimental results of this study fall within the range of results previously reported at room temperature. However, it is remarkable the huge dispersion of data reported for PVDF systems, especially in the PVDF/NMP/water system [41,44]. A relationship between the different polymer molecular weight characteristics of the polymers used in these studies with the variation of the cloud points reported was not found. On the other hand, this could be more rationally attributed to the difficulty observed to determine the experimental change in turbidity for PVDF systems compared to the clear change in PES systems, for which literature cloud points are more consistent [35,36].  with Student's t-test (p < 0.05) was found in all systems comparing points between 293 and 313 K and between 313 and 333 K. Significant difference was also found at 293 K between systems comparing the use of DMAc and NMP with the same polymer and between PVDF and PES with the same solvent.
For PVDF and PES systems, a significant influence of the temperature that displaces the cloud points to the right is observed. However, for PES systems, the displacement of the cloud point curves with the temperature was smaller than for PVDF. PES is classified as a glassy polymer with glass transition temperature (T g = 489-505 K) above the typical temperature used in phase inversion processes (298-333 K), while PVDF is in rubbery solid state (T g = 206-278 K) at this range of processing temperatures. As the cloud points represent a pseudo stable liquid phase point in the proximities to solid precipitation, this difference can be attributed to the different state of the polymers in the solid phase, being PVDF a rubbery-polymer and PES a glassy-polymer. The specific volume of the glassy polymers suffers a thermal expansion when temperature rises but there is a characteristic rigidity of polymer chains. In the rubbery state, there is a higher impact of the temperature on the increase of the specific volume, as chain mobility is more important than for glassy polymers [2]. Hence, the penetration of the solvents and therefore the solubility might be more relevant for polymers in the rubbery state (PVDF) than for glassy polymers (PES). These experimental results indicate that, from the thermodynamic point of view, tuning the NIPS process temperature could be considered an interesting approach to tailor the membrane morphology of rubbery polymers (such as PVDF). However, it might not be a reasonable choice for membranes prepared of glassy polymers such as PES.

Mathematical Modeling
The obtained solubility parameters at different temperatures are compiled in Table A4, Appendix A. It can be seen for all the components that the solubility decreased with the temperature. This reduction of the solubility parameters will affect the values of χ ij used in the simulations, so it is confirmed the necessity of calculating them considering the effect of temperature. Figure 4 presents the g 12 curves as a function of solvent concentration (φ 2 ) for the binary systems (a) DMAc-water and (b) NMP-water calculated using UNIFAC-Dortmund methodology. Table 5 presents the parameters that results from the fitting of the calculated g 12 values to the polynomial Equation (19). It can be seen that the interaction between DMAc and water decreased with the solvent concentration, while the NMP-water behaved as the opposite. Additionally, an increase of the g 12 values with the temperature is observed in both cases.
the cloud points to the right is observed. However, for PES systems, the displacement of the cloud point curves with the temperature was smaller than for PVDF. PES is classified as a glassy polymer with glass transition temperature (Tg = 489-505 K) above the typical temperature used in phase inversion processes (298-333 K), while PVDF is in rubbery solid state (Tg = 206-278 K) at this range of processing temperatures. As the cloud points represent a pseudo stable liquid phase point in the proximities to solid precipitation, this difference can be attributed to the different state of the polymers in the solid phase, being PVDF a rubbery-polymer and PES a glassy-polymer. The specific volume of the glassy polymers suffers a thermal expansion when temperature rises but there is a characteristic rigidity of polymer chains. In the rubbery state, there is a higher impact of the temperature on the increase of the specific volume, as chain mobility is more important than for glassy polymers [2]. Hence, the penetration of the solvents and therefore the solubility might be more relevant for polymers in the rubbery state (PVDF) than for glassy polymers (PES). These experimental results indicate that, from the thermodynamic point of view, tuning the NIPS process temperature could be considered an interesting approach to tailor the membrane morphology of rubbery polymers (such as PVDF). However, it might not be a reasonable choice for membranes prepared of glassy polymers such as PES.

Mathematical Modeling
The obtained solubility parameters at different temperatures are compiled in Table  A4, Appendix B. It can be seen for all the components that the solubility decreased with the temperature. This reduction of the solubility parameters will affect the values of χij used in the simulations, so it is confirmed the necessity of calculating them considering the effect of temperature. Figure 4 presents the g12 curves as a function of solvent concentration (ϕ2) for the binary systems (a) DMAc-water and (b) NMP-water calculated using UNIFAC-Dortmund methodology. Table 5 presents the parameters that results from the fitting of the calculated g12 values to the polynomial Equation (19). It can be seen that the interaction between DMAc and water decreased with the solvent concentration, while the NMP-water behaved as the opposite. Additionally, an increase of the g12 values with the temperature is observed in both cases.   Table 5.
Coefficient of solvent/non-solvent interaction parameters for the equation g 12 (φ 2 ) = a o +a 1 φ 2 +a 2 φ 2 2 +a 3 φ 3 2 +a 4 φ 4 2 +a 5 φ 5 2 +a 6 φ 6 2 . When the value of g 12 increases, the interaction between the solvent and non-solvent is lower, favoring the displacement of the binodal curve to the polymer-non-solvent axis, and therefore extending the region of homogeneous liquid phase [47]. In Figure 4, it can be seen that, overall, NMP-water g 12 presents higher values than for DMAc-water systems (worse interaction). Moreover, for the system DMAc-water, the interaction was favoured at increasing concentrations of DMAc and at lower temperatures, while for the NMP-water system, the high NMP concentrations and high temperatures decreased the interaction between the components. Accordingly, Figure 3 shows the expected displacement to the right of the cloud points with the temperature, and the larger liquid homogeneous region for the NMP systems.

Modeling Binodal Curves Not Considering the Contribution of the Ternary Interaction Term
In this section it is considered a negligible effect of the ternary interaction parameter on the phase diagram (χ 123 = 0). This modeling approach has been the most widely reported in the literature [17,18,27,48]. Firstly, χ 13 and χ 23 parameters were calculated according to Equation (11) considering α 13 = α 23 = 1 ( Table 6). It can be seen that all χ 13 values decreased with temperature. On the other hand, χ 23 values increased with temperature for all binary systems except for PVDF/NMP. Figure 5 shows the solubility parameters in the Hansen space and the radius of interaction between each pair of compounds. Interestingly, PVDF was closer to the solvents than PES, especially PVDF was more soluble in DMAc (Ra = 1.62). Besides, PES was closer to NMP than to DMAc (Ra = 4.06). The solvation could be easier for PVDF since its monomer size was smaller compared with PES, and on the other hand the phenyl groups in PES polymer contribute to the steric hindrance.  When the value of g12 increases, the interaction between the solvent and non-solven is lower, favoring the displacement of the binodal curve to the polymer-non-solvent axis and therefore extending the region of homogeneous liquid phase [47]. In Figure 4, it ca be seen that, overall, NMP-water g12 presents higher values than for DMAc-water system (worse interaction). Moreover, for the system DMAc-water, the interaction was favoure at increasing concentrations of DMAc and at lower temperatures, while for the NMP water system, the high NMP concentrations and high temperatures decreased the inter action between the components. Accordingly, Figure 3 shows the expected displacemen to the right of the cloud points with the temperature, and the larger liquid homogeneou region for the NMP systems.

Modeling Binodal Curves not Considering the Contribution of the Ternary Interac tion Term
In this section it is considered a negligible effect of the ternary interaction paramete on the phase diagram (χ123 = 0). This modeling approach has been the most widely re ported in the literature [17,18,27,48]. Firstly, χ13 and χ23 parameters were calculated ac cording to Equation (11) considering α13 = α23 = 1 ( Table 6). It can be seen that all χ13 value decreased with temperature. On the other hand, χ23 values increased with temperatur for all binary systems except for PVDF/NMP. Figure 5 shows the solubility parameters i the Hansen space and the radius of interaction between each pair of compounds. Interest ingly, PVDF was closer to the solvents than PES, especially PVDF was more soluble i DMAc (Ra = 1.62). Besides, PES was closer to NMP than to DMAc (Ra = 4.06). The solvatio could be easier for PVDF since its monomer size was smaller compared with PES, and o the other hand the phenyl groups in PES polymer contribute to the steric hindrance.   Figure A1 of Appendix A depicts the modeled binodal curves of ternary systems obtained under the above assumptions (i.e., neglecting the ternary interaction parameter). Overall, the model correctly predicts the displacement of the binodal curve to the right at increasing temperature, as it happens with the experimental data. However, the simulated curves did not fit adequately of experimental points. It is observed that all binodal curves were shifted to the left (solvent-polymer axis) in contrast to the experimental cloud points.
Different aspects, such as the polymer having a broad molecular weight distribution, or swelling or plasticizing effects occurring between the polymer and the solvent for binary mixtures have been reported to alter significantly the values of the Flory-Huggins binary interaction parameters [48][49][50]. Previous works have reported the use of a correction factor α ij in the calculation of binary interaction parameters χ 13 and χ 23 , Equation (11), to amend the deviation observed between the modeled and experimental cloud points. Therefore, a fitting procedure to estimate the correction factors was used obtaining the χ 13 and χ 23 values presented in Table 7. Simulated and experimental curves using this approach are depicted in Figure 6. Due to the fitting approach, the simulated curves are now fairly predicting the experimental points, as expected.  Table 6. Binary interaction parameters (χ13 and χ23) at different temperatures from solubility parameters with α13 = α13 = 1 and χ123 = 0.  Figure A1 of Appendix B depicts the modeled binodal curves of ternary systems o tained under the above assumptions (i.e., neglecting the ternary interaction parameter Overall, the model correctly predicts the displacement of the binodal curve to the right increasing temperature, as it happens with the experimental data. However, the simulate curves did not fit adequately of experimental points. It is observed that all binodal curve were shifted to the left (solvent-polymer axis) in contrast to the experimental cloud point Different aspects, such as the polymer having a broad molecular weight distributio or swelling or plasticizing effects occurring between the polymer and the solvent for b nary mixtures have been reported to alter significantly the values of the Flory-Huggin binary interaction parameters [48][49][50]. Previous works have reported the use of a corre tion factor αij in the calculation of binary interaction parameters χ13 and χ23, Equation (11 to amend the deviation observed between the modeled and experimental cloud point Therefore, a fitting procedure to estimate the correction factors was used obtaining the χ and χ23 values presented in Table 7. Simulated and experimental curves using this ap proach are depicted in Figure 6. Due to the fitting approach, the simulated curves are no fairly predicting the experimental points, as expected.   Table 7 shows that for PVDF systems the α 13 correction factors were near to 1, and χ 23 calculation did not require of a correction factor (α 23 = 1). For χ 13 a linear tendency with the temperature is observed with a slope of 0.0033 for both PVDF/DMAc and PVDF/NMP systems. Instead, PES systems needed significant correction in both χ 13 and χ 23 binary interaction parameters. The value of α 13 is the same for both PES systems and it increased with temperature. However, the solvent-polymer interaction, represented by χ 23 , needs a higher correction with α 23 values between 0.1 and 0.5.
Other authors propose different α ij values. Lindving et al. [51] consider a general value of 0.6, independent of the system. Wei et al. [27] estimated fitting values of α 13 = 0.83 and α 23 = 0.08 for both PES/NMP/water and PES/DMAc/water systems at 298 K; these values differ from the results obtained in this study because of the revised expression of g 12 and the solubility parameters for PES employed in the present work. Although the model adequately predicts the experimental points, the use of correction factors entails important limitations and concerns: (1) low polymer concentration points are not well adjusted to the binodal curve; (2) two fitting parameters are needed (α 13 and α 23 ) and (3) no reasonable relation is found in the adjustment of the alpha correction factors, thus limiting its extrapolation to other systems or finding a tendency. In this study, the necessity of a α 23 correction factor for PES systems may be attributed to a higher cosolvency effect between the three components, this phenomenon will be studied in detail in the following section. According to Barth et al. [52], who studied PES/DMF/water systems, the use of α 13 in the calculation of χ 13 points to the occurrence of ternary interactions between polymer, solvent and non-solvent. This approach, that is described below, has not been explored so far for the present ternary systems.

Modeling Binodal Curves Considering the Contribution of the Ternary Interaction Term
To investigate the influence of the ternary interaction term of the systems, the χ 123 ternary interaction parameter was included in the thermodynamic model. Given the lack of any methodology reported in the literature to calculate the ternary interaction parameter (χ 123 ), we decided to estimate it by adjusting the experimental cloud point data to the model, Equations (4)-(6). Figure 7 presents the values of χ 123 and the resulting simulated binodal curves.
It can be seen that the fitted χ 123 values were negative indicating a strong interaction between the components. This ternary interaction was not considered in previous literature approaches to the systems studied, where only binary interactions were deemed. Few authors have taken into account this parameter for polycaprolactone/dimethylformamide/ water system [15] and poly(ethylene-co-vinyl alcohol) (EVAL)/2-propanol (2P)/water system [53]. The negative sign of χ 123 implies a cosolvency between the solvent and water. Cosolvency is a phenomenon that has been broadly reported in ternary systems where polymers were soluble on binary mixtures of two non-solvents EVAL/2P/water and polystyrene/acetone/diethylether [53,54]. In the present study, it is hypothesized that at certain range of solvent/non-solvent compositions, the solubility of the polymer is higher than expected simply considering binary interaction between the components, which can be attributed to a cosolvency effect (solvent/non-solvent mixtures acting as a solvent). Higher absolute values of χ 123 means higher cosolvency influence. PES presents higher degree of cosolvency (higher χ 123 estimated values) than PVDF systems, probably due to its worse solubility on the solvents ( Figure 5). s 2021, 13, x FOR PEER REVIEW 13 of 20 It can be seen that the fitted χ123 values were negative indicating a strong interaction between the components. This ternary interaction was not considered in previous literature approaches to the systems studied, where only binary interactions were deemed. Few authors have taken into account this parameter for polycaprolactone/dimethylformamide/water system [15] and poly(ethylene-co-vinyl alcohol) (EVAL)/2-propanol (2P)/water system [53]. The negative sign of χ123 implies a cosolvency between the solvent and water. Cosolvency is a phenomenon that has been broadly reported in ternary systems where polymers were soluble on binary mixtures of two non-solvents EVAL/2P/water and polystyrene/acetone/diethylether [53,54]. In the present study, it is hypothesized that at certain range of solvent/non-solvent compositions, the solubility of the polymer is higher than expected simply considering binary interaction between the components, which can be attributed to a cosolvency effect (solvent/non-solvent mixtures acting as a solvent).
Higher absolute values of χ123 means higher cosolvency influence. PES presents higher degree of cosolvency (higher χ123 estimated values) than PVDF systems, probably due to its worse solubility on the solvents ( Figure 5). It is observed a linear trend of χ 123 with temperature, in the range 293-333 K, as shown in the fitting equations inserted in Figure 6. The absolute χ 123 value follows the order PVDF/NMP/water > PVDF/DMAc/water > PES/NMP/water > PES/DMAc/water. The increase of χ 123 with temperature indicates that cosolvency remains, although this effect is diminished with temperature as indicated by the reduction in the absolute values of χ 123 . Now the curves adjust better to the experimental cloud points even at low polymer concentrations.
The results obtained in the present study highlight that, despite being the most widely spread methodology, neglecting the ternary interaction among the components of the phase inversion process might be erroneous, particularly for polymers presenting low solubility in traditional solvents. It is also envisaged that the ternary interaction term would have an important contribution in describing the equilibrium of systems that incorporate solvents with low solubility properties, i.e., typical green solvents such as dimethylsulfoxide (DMSO) [55], PolarClean ® [56] or 2Pyr [6].

Conclusions
The influence of temperature in the processing of polymeric membranes by NIPS has not been sufficiently studied. In this work, four systems composed by PVDF and PES as the polymers, DMAc and NMP as the solvents and water as the non-solvent were studied at three temperatures (293, 313 and 333 K). The experimental cloud points of the systems were obtained and used to validate a mathematical model based on the Flory-Huggins theory that included the effect of the temperature.
The experimental cloud points showed that as the temperature increased the solubility region of the ternary systems was enlarged. The displacement of the binodal curve toward the polymer/non-solvent axis was more evident for the PVDF (rubbery polymer) than in the PES (glassy polymer) systems.
For the first time, the thermodynamic model developed in the present work successfully incorporated the effect of temperature to predict the binodal curves of four polymer/solvent/water ternary systems. We found that the model considering only binary interaction parameters critically deviated from experimental data. Therefore, we evaluated the incorporation to the model of a ternary interaction term that includes a ternary interaction parameter χ 123 , which was estimated for each system by fitting the experimental data to the proposed model. Negative values of χ 123 indicated the strong interaction between the three components and the presence of a cosolvency phenomenon that increased the homogeneous region where the polymer is still soluble. In all systems, the magnitude of the ternary interaction parameter decreased linearly at increasing temperatures.
The main disadvantage of the implementation of the ternary interaction term in the model is that χ 123 parameter must be estimated by fitting to experimental data. To expand the universal application of this model, future works should aim at broadening the experimental validation with other polymer systems and the search of methodologies to estimate reliable χ 123 values for different ternary systems.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.
Acknowledgments: Marta Romay is grateful to an FPI contract grant (BES-2017-081112). Arkema Inc. and Sumitomo Chemical Europe Inc. for the supply of polymers.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The UNIFAC modified (Dortmund) methodology is a way to obtain activity coefficient of binary mixtures from group contribution data.
Firstly, contribution groups are assigned to the compounds involved in the mixture. For a multitude of compounds this step can be easily done with the Dortmund data bank tool of assignation http://www.ddbst.com/unifacga.html (accessed on 7 September 2020). Secondly, volume and surface area data of each group and the binary interaction between all the groups involved in a mixture are compiled from http://www.ddbst.com/ PublishedParametersUNIFACDO.html (accessed on 7 September 2020).
Then the modified Dortmund UNIFAC methodology is applied. The activity coefficient γ of each compound is calculated using Equation (A1), this equation is composed by the combinatorial part Equation (A2) and the residual part Equation (A3).
Residual Γ (i) k : residual activity coefficient of group k in a reference solution containing molecules of type i Surface area/mole fraction ratio Volume of the pure component (i) Surface area of the pure component (i) Mole fraction of group m in the mixture Group interaction parameter Ψ mn = exp − a mn +b mn T+c mn T 2 T (A12) Table A1 compiles the information of the group assignment and volume (R k ) and area numbers (Q k ) for the compounds used in this study. Table A2 presents the interaction parameters between the groups listed in Table A1 [31,57].