Molecular simulation and mathematical modelling of glass transition temperature depression induced by CO2 plasticization in Polysulfone membranes

A sequence of molecular modelling procedure has been proposed to simulate experimentally validated membrane structure characterizing the effect of CO2 plasticization, whereby it can be subsequently employed to elucidate the depression in glass transition temperature (Tg). Based on the above motivation, unswollen and swollen Polysulfone membrane structures with different CO2 loadings have been constructed, whereby the accuracy has been validated through good compliance with experimentally measured physical properties. It is found that the presence of CO2 constitutes to enhancement in polymeric chain relaxation, which consequently promotes the enlargement of molecular spacing and causes dilation in the membrane matrix. A series of glass transition temperature treatment has been conducted on the verified molecular structure to elucidate the effect of CO2 loadings to the depression in Tg induced by plasticization. Subsequently, a modified Michealis-Menten (M-M) function has been implemented to quantify the effect of CO2 loading attributed to plasticization towards Tg.


Introduction
Polymeric membrane gas separation is an emerging process in industrial application ascribed to its various advantages in comparison to conventional technologies, such as occupying a relatively smaller footprint, chemical free, cost effective, high process flexibility and energy efficiency [1][2][3][4][5]. It has been applied extensively typically in CO2 removal from natural gas [6]. Nonetheless, the drawback of polymeric membrane that reduces its competitive edge is CO2 induced plasticization.
CO2 plasticization is explained as the phenomena whereby the sorption of highly condensable gas molecules into polymer matrix increases the free volume by acting as a "lubricant", contributing to ease of movement for chain molecules to slip over one another and thus causing polymer softening [7,8]. Plasticization increases the polymeric segment mobility and hence permeability of all gas species [9]. However, as free volume increases, the ability of polymeric material to sieve gas molecules for permeation becomes weaker and hence causes reduction in selectivity [10]. Therefore, it is important to elucidate the plasticization phenomena owing to the fact that it can potentially lead to excessive product lost that reduces profitability of the natural gas processing plant. Schematic representation of the plasticization phenomena that increases free volume of the polymeric glassy membrane is provided in Figure 1 [11]. Among the most apparent effect that characterizes the onset of plasticization upon sorption of compressed gases into glassy polymeric membrane is the depression of glass transition temperature, g T . Figure 1 Schematic representation of the plasticization phenomena resulting in facilitated polymer mobility and increased free volume in the polymer matrix, adapted from Kikic et al. [11] The fundamental elucidation of CO2 induced plasticization plays a pivotal role in enhancement of separation proficiency and the development of existing and next-generation polymeric membranes with high performance. In this context, molecular simulation has been proposed as a feasible alternative to provide insights into gas transport characteristic from an atomistic point of view, usually achieved via a coupling of molecular dynamics (MD) and Monte Carlo (MC) technique [12,13], due to the ever-growing improvement in computational technology and simulation tools [14]. Heuchel et al. (2006) pioneered an attempt to characterize molecular models of representative unswollen and swollen states of real polymeric membranes, namely polysulfone (PSF) and poly (ether sulfone) (PES) with and without sorbed carbon dioxide, and validated the methodology through comparison with actual experimental data [15]. Subsequently, detailed atomistic of the two unswollen and swollen reference states of PSF were associated to the common dual sorption model and site distribution model of Kirchheim by Hölck et al. (2006), which permitted the characterization of sorption and thereby swelling induced behavior in glassy polymers [16]. In later work by Hölck et al. (2008), a molecular modeling investigation of dilation effects induced by sorbed gas molecules attributed to plasticization effect in PSF and polyimide 6FDA-TrMPD (P14) has been conducted [17]. Based on a similar motivation, the sorption and swelling of glassy atactic polystyrene (PS) at elevated pressures had been mimicked by Spyriouni et al. (2009) through a new scheme albeit requirement of higher computational procedures and tool to achieve the designated simulation [18]. Zhang et al. (2010) carried out fully atomistic simulations to examine CO2 induced plasticization of a 6FDA-ODA (6FDA = 4, 4-(hexafluoroisopropylidene) diphthalic anhydride, ODA = 4, 4-oxydianiline) polymeric membrane through introduction of different CO2 loadings to achieve a quantitative understanding of the relation among polymer structure, its interaction with CO2 molecules and separation performance [14]. Their work had been extended by Velioğlu et al. (2012), who employed molecular simulation techniques to investigate CO2 induced plasticization in varying fluorinated polyimide membranes (e.g. 6FDA-ODA, 6FPA-DPX and 6FDA-DAM; DPX = 2,5-dimethyl-p-phenylenediamine, DAM = 2,4,6trimethyl-m-phenylenediamine) [19]. Neyertz & Brown (2014) compared the CO2 sorption and plasticization resistant through fully-atomistic molecular model of a meta-linked 6FDA-6FmDA polyimide membrane to its para-linked 6FDA-6FpDA isomer [20,21]. Computational tools were applied to estimate plasticizing effect of CO2 and to calculate CO2-and H2O-loading for three absorbable polyesters: polycaprolactone and two copolymers of poly-D, L-lactid-co-glycolid)-co polyethleneglycol by Knani et al. (2015) [22]. As highlighted earlier, the sorption of gas into polymers contributes to increment in the segmental mobility that results in enhancement of free volume, whereby the effect at the molecular level can be macroscopically observed by the depression of glass transition temperature of the polymer-penetrants mixture. Nonetheless, to date, majority of the simulation model has been confined to construction of polymeric membrane structure ascribed to CO2 plasticization, while molecular simulation of reduction characteristic in glass transition temperature as a direct consequence of plasticization followed by mathematical modeling of the aforementioned phenomena has received less scrutiny. Therefore, the objective of current work is to develop a molecular modeling procedure to simulate an experimentally validated membrane structure characterizing the effect of CO2 plasticization, whence it can be subsequently employed to elucidate the depression in glass transition temperature attributed to different loadings of CO2. A modified Michealis-Menten (M-M) function that has been previously adopted to quantify thickness dependent g T in published work by Kim et al. (2000) [23], has been improvised to characterize the effect of CO2 loading attributed to plasticization towards the g T . To the best of our knowledge, our work is the pioneering effort to address such g T depression in the presence of CO2 plasticization through employment of molecular simulation tool.

Methodology
The methodology is subdivided into two subsections, whereby the first is molecular simulation to construct the molecular structure of the PSF polymeric membrane (with and without CO2 plasticization effect), and the second is procedure for determination of glass transition temperature in the polymeric membranes.

Atomistic Packing Models
The molecular structure has been simulated through adaptation of Materials Studio 8.0 developed by Accelrys Software Inc [24]. Polysulfone (PSF), which is a commercial membrane commonly adapted for gas separation, has been simulated employing a series of molecular dynamics (MD) procedure in Materials Studio 8.0. The repeating unit of the PSF polymeric membrane is provided in Figure 2. The Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies (COMPASS) force field has been adopted consistently. In addition, the Ewald method with an accuracy of 0.001 kcal/ mol has been adopted to describe the electrostatic interactions, while the van der Waals interaction has been characterized via the Lennard-Jones-6-12 function with a cut-offdistance of ~11 Ǻ (spline width of 1 Ǻ and buffer width of 0.5 Ǻ), whereby this value is selected since it is less than half of the expected cell length.

Unswollen Polysulfone (PSF) Polymeric Membrane Structure
A single polymer chain with 20 repeat units is constructed. The initial polymeric chain has been located in the Forcite module of Materials Studio 8.0 and has been subjected to energy minimization and geometry optimization. The COMPASS force field has been adopted alongside the smart algorithm, which is a combination of the steepest descent; adjusted basis set Newton-Raphson (ABNR) and quasi-Newton algorithms in a cascading manner, in order to refine geometry of the initial polymeric chain. Later, the polymeric membrane chain has been folded into Amorphous Cell module adopting Construction task. The polymer chains have been embedded in the hypothetical cell under the periodic condition at initial density corresponding to 70% of the targeted experimental density. Similarly, the COMPASS force field has been adopted to pack the polymeric chains.
Then, the initial constructed atomistic configuration has been subsequently minimized and optimized adopting a series of protocols. Firstly, a 10000 step energy minimization has been conducted in order to remove any undesirable configurations including overlapping and close contact. After this stage, an annealing procedure has been performed on the polymeric structure by adopting the temperature cycle protocol embedded within Forcite module. Throughout the cycle, the system has been heated and cooled back with an interval of 293.15 K between 353.15 K and 653.15 K, which is well above the glass transition temperature of PSF. At each temperature, 100 ps isothermal-isobaric ensemble (NPT) has been conducted. At this step, in order to control the temperature at the designated heating temperature and pressure of 1 atm, the Nose thermostat with Q ratio of 0.01 and Berendsen barostat with decay constant of 0.1 ps have been employed continuously. Subsequently, in order to achieve the ideal structure, which is the lowest energy configuration with the most realistic geometry, a molecular dynamics equilibrium run has been implemented on the amorphous cell structure in the NPT ensemble with a total simulation time of 500 ps. Similarly, the pressure of the system has been maintained at 1 atm while the temperature is fixed at a constant value of 298.15 K with Nose thermostat and Berendsen barostat. Throughout this step, the equation of motions has been integrated by the velocity Verlet algorithm with a time step of 1 fs for all simulation conditions. When approaching the endpoint of the NPT run, an additional 500 ps of Canonical (NVT) ensemble at temperature of 298.15 K has been conducted on the equilibrated polymeric structure. The NPT-NVT molecular treatment has been repeated until changes in the successive density values are within predefined tolerance. The unswollen polysulfone membrane structure is named PSU throughout this study.

Swollen Polysulfone (PSF) Polymeric Membrane Structure
Additional packing model has been prepared for the swollen state of CO2/ polymer system attributed to plasticization effect, which corresponds to plasticization pressure ( pl P ) in Bos et al. (1999) experiment observation of 34 bar [9]. To comply with the PSF polymeric chain of 20 repeat units and solubility coefficient of 1.4 cm 3 (STP)/cm 3 .bar taken from Bos et al. sorption experiment, a number of 16 CO2 molecules are required to be packed into the simulation box in addition to the PSF chain in order to induce the relaxation characteristic attributed to plasticization. A different number of CO2 gas penetrants (i.e. 8, 16, and 24) have been packed to elucidate the effect of CO2 loadings to the glass transition temperature of PSF polymeric membrane structure, whereby they have been named PSU8, PSU16, and PSU24 respectively in the remaining section of the paper. The applied basic techniques of packing and equilibration as described in section 2.1.1 have been employed in this context. After packing the PSF chain with designated CO2 loadings into Amorphous Cell module adopting Construction task, the polymeric membrane structures have been subjected to a similar series of molecular treatment protocol in previous section.

Glass Transition Temperature
In this study, g T s of PSF samples have been determined by mimicking the heating and cooling protocols in laboratory scale adapting a series of thermodynamic treatment in Forcite Module. Firstly, the optimized and equilibrated configuration for each dimension has been subjected to an additional Canonical (NVT) ensemble at 308.15 K with a time step of 1 fs and total simulation time of 10 ps by framing the output every 1000 steps. This procedure is aimed to obtain the trajectory files of PSF polymeric films with 10 frames for each thickness, such that an average g T can be deduced to increase accuracy of the computed value when the series of thermodynamic treatment is iterated while calculating an independent g T for each frame. Then, the individual frame located within the PSF trajectory has been exposed to gentle heating from 300.15 K to 500.15 K, which surpasses that of the bulk glass transition temperature of PSF polymer, with an interval of 1 K. Subsequently, a 100 ps NPT dynamic ensemble has been conducted at the designated heating temperature and operating pressure. Thereafter, the system is cooled down from 500.15 K to 300.15 K with the temperature interval of 1 K while computing density of the structure at each temperature. This protocol is looped over all frames contained in the trajectory file and eventually the values are averaged at the end of the simulations.

Results and Discussion
In this section, the molecular simulation results related to PSF polymeric membranes at unswollen and swollen conditions have been presented, whereby it has been subdivided into two major subcategories, such as 1) molecular structure and density and 2) glass transition temperature.
3.1. Molecular Structure As described in section 2.1, molecular dynamics simulation has been executed for all PSF polymeric films by keeping the operating parameters at the designated laboratory condition while altering the number of CO2 molecules within the membrane structure, while the other configurations are constantly updated in quest of determining the most probable polymeric membrane film with optimized packing and molecular arrangement. The example of molecular structure within the unswollen and swollen PSF membrane with 8, 16 and 24 CO2 loadings are depicted in Figure 3 (a), (b), (c) and (d) respectively. As it is illustrated in Figure 3 (b), (c) and (d), the CO2 gas molecules are randomly dispersed throughout the void space within the PSF polymeric matrix, with the configuration corresponding to the optimized and most probable intermolecular interactions between the diffusants and polymeric chains. Since the system has been initialized from lower density without setting any constrictions throughout the molecular dynamics treatment, the evolution of structure, typically the structure

1234567890
International density, to a stable value provides intuitive reasoning that the polymeric membrane has converged towards the most plausible configuration. Hence, by analyzing physical property of the finalized molecular structure to theoretically attained value and trend with respect to effect of CO2 loadings, which have been measured through either published experimental or molecular simulation work, accuracy of the aforementioned molecular simulation methodology can be validated. Physical properties of the constructed packing models are given in Table 1. In order to verify accuracy of the simulation methodology, the simulated properties have been validated with actual experimental observations. For instance, it is found from Table 1 that the simulated density of PSU is in close agreement with previously reported density of bulk PSF with a value of 1.24 g/ cm 3 by Ahn et al. (2008) [25]. The good compliance between simulated and experimentally observed condition supports the claim that the unswollen PSF polymeric structure has been constructed via high accuracy simulation procedure since the system has been ramped from a low density configuration without confining any constraints and boundaries throughout the molecular dynamics process. Small percentage error of -1.45% between simulated and experimentally observed PSF density can be rationalized through the assumption in molecular simulation, whence cut off distance has been applied throughout the simulation cell that deemed long range molecular interaction to be negligible. As for swollen PSF polymeric membrane structure with presence of CO2, the accuracy has been demonstrated through excellent accordance between the solubility characteristic (e.g. 1.2651 versus 1.38 cm 3 (STP)/cm 3 .bar) in simulated PSU16 and experimentally observed sorption experiment taken from Bos et al. (1999) [9] respectively. In addition, from Table 1, it is depicted that the packing density of molecular structure increases with increment in CO2 loadings. The observation tallies with characteristic of single packing models of PSF, as reported by Heuchel et al. (2006) [15], whereby PSU molecular structure with CO2 exhibits greater density as compared to its pristine unswollen PSU configuration attributed to the presence of CO2 that increases weight of the entire molecular system. On the other hand, density of the membrane matrix experiences a decrement when CO2 loadings are further increased. The reduction in density of membrane matrix has been acclaimed through enhanced interaction between CO2 gas molecule and polymeric matrix, which contributes to augmented swelling and relaxation in polymeric membrane structure [9]. The claim has been supported via incremental in box length that characterizes the degree of swelling, whereby percentage dilation rises with number of CO2, which has been computed and tabulated in Table 1 as well.

Glass Transition Temperature
The variations in specific volume versus temperature for PSF polymeric membranes with different CO2 loadings have been plotted and provided as examples in Figure 4. As it can be seen from Figure 4, all the polymeric membranes experience similar behavior with changes in temperature regardless of CO2 loadings. Initially, the specific volume increases linearly with increment in temperature, and then shows an abrupt alteration in the value before continuing to embark in another linear region. Change in linear relationship is demonstrated through the difference in slope between the two curves, whereby the first at lower temperature is representative of the glassy state region, while the latter describes the rubbery state. The point at which the glassy and rubbery linear correlation meets to form an intercept provides graphical representation of the glass transition temperature, g T . The simulated behavior is consistent to experimental observation reported by Zoller et al. (1978), who investigated the pressuretemperature-volume relationships in bulk PSF over a wide range of operating conditions [26]. Figure 4 Effect of CO2 loadings on specific volume of the polymeric film with respect to alteration in temperature, whereby T1 = 121.6 °C, T2 = 151.2 °C, T3 = 169.9 °C, T4 = 182.5 °C Specific volume of the PSF structure has been found from Figure 4 to be consistently higher for configuration with CO2 gas molecules attributed to plasticization induced dilation as described earlier that causes increment in molecular spacing. The glass transition temperature is found to experience a depression behavior when CO2 loadings are further increased. The reduction in glass transition temperature, g T , can be rationalized through enhanced mobility of polymeric chains in the vicinity of CO2 gas molecules, which constitutes to amplified softening and movement for chain molecules to slip over one another, thus contributing to suppressed g T in comparison to the unswollen PSF polymeric structure.
It has been suggested that the reduction in g T seems to follow the form of growth with saturation, T , is the bulk glass transition temperature,  characterizes a material constant that describes the function growth saturation rate, while CO2 is resemblance of the amount of carbon dioxide gas molecules in simulation box. (1) As depicted in Figure 5, the suggested correlation is able to provide a satisfactory compliance with simulated glass transition temperature, which has been demonstrated through a good R 2 value.

Conclusion
In present work, a molecular modeling procedure has been proposed to construct unswollen and swollen PSF membrane structure with the presence of CO2 in order to capture the effect of plasticization to membrane morphology. Accuracies of both the unswollen and swollen membrane structures have been verified through excellent accordance between simulated and experimental observed physical properties. The good compliance ensures applicability of the simulated membrane structure to be employed to study the impact of CO2 loading on depression of glass transition temperature attributed to plasticization phenomena. A revised form of the Michealis-Menten (M-M) function, which has been previously adopted to quantify thickness dependency in glass transition temperature, has been employed in current work to correlate the effect of CO2 loading to g T . To the best of our knowledge, this work is among the first to tackle the impact of numbers of CO2 towards the depression of g T through molecular simulation technique, which has been employed in proposal of empirical model to address the plasticization induced g T depression. The findings of current work can be sufficiently employed to overcome the barrier, cost and time in preparation and testing of CO2 induced plasticization in membrane at the laboratory scale since it appears to be relatively convenient and straightforward to control operating conditions of the simulation.