Numerical Simulations of Gas Sorption Experiments in Polymers: Inﬂuence of Aspect Ratio and Pressure Increase Rate on the Determination of Diﬀusion Coeﬃcient

Diﬀusion and sorption phenomena of gases in polymers are of relevance for a number of technological applications, e.g., polymer foaming and membrane technology. The diﬀusion coeﬃcient and the solubility, which depend on temperature and applied gas pressure, can be determined using the buoyancy method. In this work, the inﬂuence of two important parameters are investigated, i.e., the dimensions of the polymer sample and the pressure increase rate, on the accuracy of the determination of the diﬀusion coeﬃcient D using gas sorption measurements. Numerical simulations on eﬀectively two-dimensional diﬀusion phenomena in polymers are performed in order to obtain “ideal” experimental data. This set of “ideal” experimental data for polymer samples with speciﬁed geometrical dimensions is used for determination of the diﬀusion coeﬃcient based on the analytical solution of Fickian diﬀusion in one spatial dimension. Diﬀerent ﬁtting procedures for determination of the diﬀusion coeﬃcient are compared. Because of the ﬁnite time, which is necessary for pressure increase in a sorption experiment, the inﬂuence of pressure increase rate is also studied. This analysis reveals that an aspect ratio larger than 32 and a time for pressure increase smaller than approximately 3 × 10 − 9 m 2 D − 1 is necessary for a reliable determination of diﬀusion coeﬃcient D .


Introduction
Diffusion is a transport phenomenon which is driven by concentration gradients. The diffusion coefficient plays a key role in polymer-gas interactions such as sorption of carbon dioxide (CO 2 ), polymer swelling, polymerization, monomer recycling, stripping, drying, coating, foaming, polymer plasticization, and polymer purification. [1,2] This wide scope of relevance requires a DOI: 10.1002/mats.202100016 profound understanding of diffusion processes in polymers. [3][4][5][6][7] The diffusion coefficient depends on concentration, temperature, and pressure. In general, diffusion is described by Fick's law. Nevertheless, non-Fickian behavior is regularly experienced in polymers below the glass transition temperature due to chain relaxation within the polymer. [8,9] Factors such as the molecular size and physical state of the diffuser, the morphology of the polymer, temperature of media, and gas pressure have an effect on the time scale of diffusion of a gas in a polymer. [10][11][12][13][14] Barrer et al. investigated the edge effect of samples in time-lag measurements. An aspect ratio larger than 10 is ideal for permeability measurements due to a negligible edge effect. [14] Tendulkar et al. determined the composition-dependent diffusivity of carbon dioxide in polyethylene depending on pressure and mass concentration at each temperature. [12] The saturation time as a function of the geometry of the sample was studied. [13] Alentiev et al. investigated the relationship between various transport parameters such as diffusion, permeability, the activation energy of diffusion, and cohesion energy density in glassy polymers. [15] Ushiki et al. measured the solubility and the diffusion coefficient of CO 2 in poly(ethylene-co-acrylic acid) copolymer at high temperature and high pressure. [16] Wessling et al. developed experimental methods by measuring the sorption kinetics at low and high pressures. The authors showed that the determination of sorption kinetics is not possible because of momentous reduction in absolute pressure over time. This change in pressure over time creates a complex time-dependent condition of the differential equation that describes the process of diffusion of penetrating molecules to the depth of a polymer film. In the presence of CO 2 , in particular, extrapolation to high-pressure applications is questionable. [17] Kundra et al. determined the composition-dependent diffusivity and solubility of carbon dioxide in polypropylene at different temperatures and pressures by using the pressure decay method. Pressure has a nonlinear effect on the value of the diffusion coefficient. A rise in pressure increases the frequency of intermolecular collision, which causes gas diffusion; however, increased pressure reduces the intermolecular distances that prevented diffusion.
In addition, the pressure change could lead to some important changes in the structure of the polymer matrix and affect mass transfer. [18] Many researchers have endeavored to explain specific mechanisms by which diffusion occurs in polymeric systems. Nevertheless, there is not a united theory to clarify this phenomenon completely. [19] The value of the diffusion coefficient is a material parameter that can be obtained by sorption measurements that include the time-dependent measurement of the mass of sorbent. This approach provides a direct and sensitive measurement of the gas uptake using a microbalance. It is time consuming and costly to identify the critical processing parameters within the sorption measurement. [1] To determine the diffusion coefficient, different methods have been used such as buoyancy method, [20] time-lag method [21] and kinetics of swelling and motion of diffusion front, [22] and Raman spectroscopy. [23] Recently, the use of simulation methods such as molecular dynamics or finite-elements for design and optimization of gas sorption processes in the field of membrane technology research greatly increased in parallel with the development of such programs and their increased availability. [24][25][26][27] In-house codes, several commercial and open-source computational fluid dynamics (CFD) programs are widely present in the market. Among these programs, COMSOL uses the finite element method (FEM) whereas ANSYS-CFX and OpenFOAM apply the finite volume method (FVM). Each CFD program has its advantages and disadvantages depending on the application field. [28,29] Some researchers have compared different finite element programs. For instance, Jeong et al. specified that FEM solutions are sensitive to the type and quality of the mesh, and require approximately five times more calculation time than FVM codes. [28] Marjani et al. investigated a mechanistic model for the design and optimization of membrane-based gas absorption processes via COMSOL and compared the model predictions with experimental data. [24] For complex problems, a full three-dimensional model may require a heavy computational simulation which makes the analysis time consuming. To this point, three-dimensional modeling usually is avoided if one can reduce the spatial dimensions because of symmetry. Accordingly, 2D-axisymmetric models are used to show the behavior of three-dimensional problems that have high accuracy and have lower costs than complete threedimensional simulations. [30,31] This work is devoted to a numerical analysis of gas sorption experiments. The objective of our study is the analysis of the influence of the three-dimensional geometry on the determination of the diffusion coefficient using the one-dimensional analytical solution based on Fickian diffusion which is generally used for analysis of experimental data. In the present study, numerical simulations in a 2D-axisymmetric geometry are applied to create "ideal" experimental data. Three different fitting methods exist to determine the diffusion coefficient, which are the use of a short-time approximation, the evaluation based on time t 90 until a saturation of 90% is achieved, and the exact solution of the one-dimensional diffusion equation. The sorption curves were evaluated using these three methods based on Fickian diffusion in one dimension. The obtained results reveal quantitatively that the dimensions of the sample influence the numerical determination of the diffusion coefficient. Moreover, the influence of the pressure increase rate is examined because of the finite time for  pressure increase in a sorption experiment. The results of this study are beneficial to improve the knowledge of the optimum measurement parameters in a sorption experiment.

Simulation Details
In this work, COMSOL Multiphysics software (Version 5.6) was used to develop a mathematical model to describe diffusion of a gas through a solid polymer. In this work, a cylindrical sample was chosen which can be prepared, e.g., by compression molding. The method of finite elements was used to solve Fick's second law of diffusion for a disc with cylindrical geometry. The cylindrical sample is shown in Figure 1 where R and h, respectively, denote the radius and thickness of the sample. The simulation was executed using the transport of diluted species in a time-dependent study. It is used to compute the concentration field of a dilute solute in a solvent. The driving force is the gradient of concentration (Fick's law), and the time-dependent local concentration is calculated using the diffusion equation for a species dissolved in a gas, liquid, or solid. The backward differentiation formula method as the time-dependent solver was used. Geometric parameters such as the radius and the aspect ratio of the sample were examined to identify the saturation time of gas diffusion. The influence of the aspect ratio was investigated by varying the radius R in the range from 1 to 20 mm. The thickness h of the sample was constant and equal to 0.5 mm. The geometry of the mesh is shown in Figure 1. The mesh for the whole domain was defined using the parameters in Table 1.
The mapped mesh maps a regular grid defined on a logical unit square onto each domain, which is based on transfinite interpolation. The mapped nodes used by the element ratio 10 000 in the z-axis with symmetric distribution and the element ratio 1000 in the radial direction.
The diffusion equation was solved with MUMPS solver (multifrontal massively parallel sparse direct solver) within a cylindrical coordinate system (r, ϕ, z). In a 2D-axisymmetric simulation the boundary condition in the center is ∂c(r = 0, ϕ, z, t)/∂r = 0 for t ≥ 0. In cylindrical coordinates, the diffusion equation reads where c(r, ϕ, z, t) is the local concentration (%) of the gas at time t and D (m 2 s −1 ) the diffusion coefficient (diffusivity). For the numerical simulations, the time t is rescaled by the diffusion coefficient via t′ = Dt. Furthermore, we have ∂c/∂ϕ = 0 due to symmetry. Thus, Equation (1) simplifies to In this work, r, z are in the range 0 ≤ r ≤ R with R ≤ 20 mm and |z| ≤ 0.25 mm. The total simulation time in the rescaled time unit was 3 × 10 −7 m 2 . The time step was 10 −9 m 2 . In some simulations, the total time varied depending on dimensions and boundary conditions.
The gas solubility c stat in a polymer depends on temperature T and applied pressure p. Henry's law is often used for polymers below the glass transition temperature and implies a linear relation between concentration and pressure p stat : with Henry's constant H. In this work, we assume that the saturation concentration is achieved at all surfaces (top, bottom, and wall of the cylinder) of the sample in the beginning and set as boundary condition a constant surface concentration c 0 = 1% for t ≥ 0 (top and bottom surface and cylinder wall). Because of linearity of the diffusion equation for a constant D value and of the boundary conditions, this can be done without loss of generality. At the beginning of a sorption experiment, the pressure should be increased instantaneously from ambient pressure to a set value, e.g., 50 bar. However, this instantaneous increase cannot be realized in practice. In an experiment, a finite pressure increase rate exists and the device starts to run from zero until it reaches the desired pressure p stat . This pressure increase rate has an effect on the determination of the diffusion coefficient, since the gas diffusion starts from the beginning of pressure increase. The pressure increase rate depends on the required time t p, stat to reach the set pressure and equals p stat /t p, stat . The pressure increase rate in the sample will be taken into account as a timedependent boundary condition for the surface concentration c 0 in the mathematical analysis of the sorption kinetics and in the simulation part.

Result and Discussion
In the first step, we compare our numerical results with the work of Sun and Mark [32] as a reference to illustrate the validity of Fickian diffusion. In their work, CO 2 sorption in poly(ethersulfone) Figure 2. Average concentration c mean as a function of time t as determined by numerical simulation and experiment at a temperature T of 25°C and a pressure p of 57 bar. [30] The sample was a rectangular strip with dimensions of 5 × 10 × 1.5 mm 3 . The value of D was 3 × 10 −12 m 2 s −1 .
as a function of time was measured at a pressure of 57 bar and room temperature (Figure 2). The authors used a rectangular strip with dimensions 5 × 10 × 1.5 mm 3 . The time t 90 which is necessary until the average concentration achieves 90% of its saturation value was 50 h. At saturation, the equilibrium of CO 2 sorption in poly(ethersulfone) reaches 8%. Since diffusion takes place most rapidly toward the thinnest dimension (here parallel to the z-axis), the time-dependent concentration in a pure onedimensional situation is generally used for the further analysis because of simplicity. If the thickness of the rectangular plate is denoted by h, then the analytical solution of the diffusion equation reads [33] c (z, t) = c 0 with |z| ≤ h/2. The average concentration c mean (t) = For a simple rectangular geometry, one finds for the value of t 90 [13] t 90 = 0.21201h 2 ∕D.
In the study of Sun and Mark, the diffusion coefficient is experimentally given by D = 3 × 10 −12 m 2 s −1 . [13] Poly(ethersulfone) is an amorphous polymer. In case of a semicrystalline polymer, the diffusion coefficient also depends on the degree of crystallinity. The mean gas concentration in the rectangular polymer sample, see Figure 2, shows that the numerical simulation with the experimental value of the diffusion coefficient is in a good agreement with the experimental result. At small times, the measured data are smaller than the simulated ones, possibly because of a small gas loss during the weighing procedure.
In order to discuss diffusion in all three spatial dimensions, we consider a cylindrical sample with R = 4 mm and h = 0.5 mm. The role of various boundary conditions on the saturation process is investigated. Figure 3 shows simulation results using four different boundary conditions, namely a constant surface concentration c 0 = 1% at the surfaces with a) r = R and z = ±h/2 , b) only z = ±h/2 , c) r = R and z = +h/2, and d) only r = R. At the other surfaces, a no flux boundary condition holds. Figure 3 reveals the strong effect of different boundary conditions. The gas diffuses along the thickness and the radius. However, for large radii the radial diffusion is much slower compared to the diffusion parallel to the z-axis. Since the most rapid diffusion takes place across the thickness, the diffusion in a cylindrical geometry with a high aspect ratio approximately follows Equation (4). However, the maximum aspect ratio generally is experimentally limited.
For this reason, the influence of the aspect ratio on the saturation time is investigated in this work with a constant surface concentration at all outer boundaries. The evolution of local concentration at various times for z = 0 is presented in Figure 4. The average concentration increases with time to attain its steadystate value of 1%. The concentration gradient in the sample smoothly varies in the sample based on its distance from the sample surface. As it is shown in Figure 5, the counter plot reveals that diffusion takes place starting from the boundaries. The minimum of concentration is achieved in the center of the cylindrical sample.
The calculated average concentration as a function of time for various radii is shown in Figure 6 which reveals the effect of aspect ratio on the diffusion process. A larger sample radius leads to a larger saturation time. For R > 6 mm, the gas concentration almost saturates in the same time range.
In general, there are three methods to determine experimentally the diffusion coefficient (Figure 7), namely (i) the short time approximation, (ii) the determination of time t 90 , and (iii) the fit of the complete analytical solution of Equation (5) to the experimental data. [13] The short time approximation implies fitting the sorption curve using the result for Fickian diffusion Equation (5) in the initial time interval: where c mean is the average concentration of gas adsorbed in the polymer at time t, c 0 is the equilibrium concentration (t → ∞), and h is the thickness of the sample. [2] Equation (7) is an approximation, which is valid at short times (i.e., for values c mean < 0.6c 0 ). To evaluate the short time interval, an average concentration smaller than 30% is considered by fitting the sorption curve.
To demonstrate the dependence of the aspect ratio of the sample on the diffusion coefficient (our results obtained from the simulations), Equation (7) is fitted to the numerical data using the diffusion coefficient as fit parameter. The fitted value of the diffusion coefficient in the rescaled time units is denoted by D fit . The exact value of D fit is unity for the rescaled diffusion Equation (2). The fit curves are shown in Figure 8, and D fit for various radii is presented in Table 2.
The second method to estimate the diffusion coefficient is based on the time t 90 . By evaluating the time t 90 , D fit is gained, and the result is expressed in Table 2.
The third method to determine the diffusion coefficient is based on the exact solution of the diffusion equation in one dimension in the whole time interval. In this case, the diffusion coefficient for various radii R is obtained via fitting Equation (5) (using summands up to the order n max = 13) to the result of the simulations, see Table 2. By fitting an equation to each curve, we found the diffusion coefficient in the range 1.02-1.61. Generally, all three methods yield a fit value larger than unity because of the diffusion in the radial direction which is neglected by the one-dimensional solution. Figure 9 shows the fitted diffusion coefficient D fit as a function of aspect ratio 2R/h for the three methods. The results of three methods are close to each other. All methods follow the same www.advancedsciencenews.com www.mts-journal.de  trend with 2R/h ratios more than 30. However for the same aspect ratio, the result of the t 90 method is closer to the exact value of unity. This method is thus more accurate as well as simpler to estimate for the given parameters. The result of the short time approximation has a higher value than the other two methods be- cause of the initial fast radial gas diffusion. In general, the higher aspect ratio is, the more accurate the result will be. An aspect ratio larger than 32 leads to an experimental error below 10% and is in many cases sufficient for determination of the diffusion     coefficient. Barrer et al. suggest an aspect ratio larger than 10 for time-lag experiments. [14] In these experiments, a cylindrical sample is chosen. The data of the fitted value of the diffusion coefficient can be empirically described by: with the aspect ratio x = 2R/h , see Table 3. In particular, at low aspect ratios the evaluation method strongly influences the numerical determination of the diffusion coefficient. Based on Equation (8), the experimental error of less than 10% is achieved if the aspect ratio fulfills the condition x > 10/b. In this study, the effect of pressure increase rate a is evaluated by taking the time-dependent boundary condition with the rescaled time to reach the set pressure t′ p, stat = 1/a . In Equation (9) we have omitted the unit %. This boundary condition corresponds to a linear increase of pressure. In our simulations, the minimum value of a is 0.034 × 10 8 m -2 which means the steady-state pressure is achieved after 3 × 10 −9 m 2 . The result of our simulation for different pressure increase rates is shown in Figure 10. A low pressure increase rate causes the amount of gas dissolved in the polymer very slowly and almost increases linearly because of the linear pressure increase. By increasing the pressure increase rate, saturation is achieved at earlier times. In summary, the diffusion process at high and low pressure increase rates differs significantly. Obviously, a higher  pressure increase rate yields more accurate results if the common fit Equations (5)-(7) are used. Figure 11 presents the fitted value of diffusion coefficient obtained by analysis of t′ 90 as a function of pressure increase rate a which shows that a sufficiently high pressure increase rate is necessary in order to obtain the diffusion coefficient with a sufficient accuracy. Our simulations show that the pressure increase rate should be larger than 3 × 10 8 m -2 which implies a time for pressure increase being shorter than 3 × 10 −9 m 2 /D. In the Appendix, an analytical consideration yields the value 100/h 2 for a minimum reasonable pressure increase rate. In our simulations we have 100/h 2 = 4 × 10 8 m -2 which agrees well with the results of the numerical simulations in Figure 11.
Because of the finite value of the pressure increase rate, experimentally, the point of time zero is not well defined. This is of high numerical relevance if the short time approximation is used. Figure 12 demonstrates the influence of a "wrong" point of time zero, if, e.g., the time interval of pressure increase is too long. The figure indicates that D fit varies if a small shift of the time scale appears. The other two methods are affected to a lower degree (in particular, the method t 90 ), since they also consider the experimental data at long times and not only at short times.
Our numerical analysis shows the limitations of the analysis of the experimental data. However, in practice experimental restrictions exist. The aspect ratio cannot be extended in an arbitrary manner, and a lower thickness of the sample yields a lower mass which increases the relative error of the weighing procedure. Figure 13 schematically depicts the influence of sample mass on the accuracy of gas sorption measurements. In this Figure, we assume that the sample has a radius R of 0.564 cm and a density Polymer of 1 g cm -3 . If the sample radius is constant, then a higher thickness h of the sample yields a higher mass which increases the measurement time (approximated by Equation 6) which is proportional to h 2 . However, the accuracy of weighing (Δm/m) is inversely proportional to the total mass. The error given by the fit procedure (method time t 90 ) follows from Equation (8) and equals h/(2bR). Because of the finite time for pressure increase, a larger sample is preferred. In summary, in a real experiment a compromise between these considerations has to be made.

Conclusion
In this study, the diffusion of gas molecules in a polymer sample is numerically investigated. The purpose of this work is to show how to reduce the possible experimental error as much as possible by choosing a proper aspect ratio of the sample. The influence of aspect ratio of the sample is one of the main factors that has an effect on the experimentally determined diffusion coefficient D. The diffusion coefficient was evaluated by three different methods. Our analysis reveals that the t 90 method is the most accurate method, since this method is less sensitive to a wrong choice of the point of time zero and is associated with the lowest experimental error. Furthermore, it has the advantage of simple calculation. The pressure increase rate also affects the evaluation of sorption data. In an "ideal" experiment the aspect ratio should be larger than 32 and the time for pressure increase smaller than 3 × 10 −9 m 2 D −1 . For example, for a diffusion coefficient of 3 × 10 −12 m 2 s −1 the time for pressure increase should be less than approximately 15 min.

Appendix A Diffusion Equation with Time-Dependent Boundary Conditions
In the following, a criterion for a minimum pressure increase rate is derived. The influence of pressure increase rate can be studied analytically in one spatial dimension. The exact solution of Equation (1)  Therefore one can define a simple criterion for a minimum reasonable pressure increase rate by a ≥ 100/h 2 ≫ 10/h 2 .