Mathematical Modelling of the Permeation of Gases in Polymers

Résumé — Modélisation mathématique de la perméation des gaz dans les polymères — Ce travail concerne la caractérisation des coefficients de transport de gaz dans les polymères à partir d’expériences de perméation. Dans un premier temps, la démarche de modélisation physique et mathématique est présentée. Deux modèles de dépendance en concentration sont utilisés pour décrire la diffusion des gaz dans le polymère. La démarche de modélisation mathématique aboutit à la construction d’un problème d’optimisation non classique. Un travail complémentaire à cette démarche concerne l’étude de sensibilité des paramètres qui caractérisent le coefficient de diffusion et la concentration maximale de gaz, inconnue dans les essais de perméation. Cette méthode d’optimisation ne dépend pas de l’initialisation des paramètres et permet de les déterminer de façon unique. Abstract — Mathematical Modelling of the Permeation of Gases in Polymers — In this work, we focus on the experiments of permeation in order to characterise the transport coefficients and how they depend on the gas concentration. First, a physical and mathematical modelling approach is presented. Two models are used in order to describe the diffusion of gases in polymers. The mathematical modelling leads thereafter to the construction of an unusual optimisation problem. A work complementary to this approach relates to the study of parameter sensitivity, which characterises the coefficient of diffusion and the maximum gas concentration, unknown in the permeation tests. This method of optimisation does not depend on the initialisation of the parameters and makes it possible to determine a unique set of parameters.


INTRODUCTION
The understanding of the mechanisms of damage in polymer sheaths during gas decompression passes by the knowledge of the gas transport phenomena in polymers, by the study of the influence of gas absorption on material properties, and by modelling the material behaviour during a decompression. The mathematical theory of diffusion (Crank, 1968) in an isotropic system is based on the assumption of proportionality between the diffusing flow of molecules (that is the quantity of species that cross a membrane by unit of time and surface) and the concentration gradient between both faces of the membrane. It is Fick's first law: (1) D is the diffusion coefficient (cm 2 /s).
In the case of unidirectional diffusion, by expressing the fact that the quantity of gas held by unity of volume of polymer is proportional to the increase of the concentration with time, it is possible to write the equation: (2) When the interactions between the polymer chains and the gas molecules are strong, for example, during the diffusion of carbon dioxide (CO 2 ) in fluoride polymers, the diffusion coefficient depends on the gas concentration inside the material. In the case of a polymer-gas system, Relation (3) links the coefficient of diffusion to thermodynamical variables of the system: The objective of this study is to describe, in the case of permeation experiments, the variation of D according to the concentration, to the temperature T and to the pressure p.
If D 0 (T, p) represents the gas diffusion coefficient in the upstream side of the polymer, one can write, in a general way, D as follows: where φ(C, T, p) is a function of coupling of the concentration C with the temperature T and the pressure p.
As indicated in Klopffer and Flaconnèche (2001), the function φ(C, T, p) can take two forms: -an exponential form: -a linear form: Diffusion of gas in polymer membrane of ω section and l thickness.
The experimental device ) which describes gas permeation through a polymer membrane allows to perform a permeation test on a membrane of thickness l and of constant section ω (Fig. 1). Only one face of the sample is in contact with the studied gas.
The analysis will be limited to the unidirectional case. This choice is justified by the fact that the thickness l of the membrane is very small with regard to its diameter. At a temperature T and a pressure p fixed, the gas penetrates and crosses this membrane. A gradient of concentration is set up inside according to the direction (0x), according to the following unidirectional equations: The first boundary condition of System (7) represents the concentration of gas saturation in the polymer. It is expressed in cm 3 of gas, measured in the STP conditions (standard of temperature and pressure), by cm 3 of polymer. For multiple gas-polymer systems, this concentration is unknown. In practice (Flaconnèche, 1995), the quantity of diffused molecules crossing the membrane is measured at x = l. This quantity is defined by: φ Gradient of concentration Gas quantity crossing the membrane in function of time.
When the time t is large enough, the experimental curve of Q(l, t) is a straight line, the slope a(T, p) of which depends on the temperature T and on the pressure p (Fig. 2). The intersection of this line with the time axis defines the time lag τ(T, p).
To validate the general model (4), the objective consists in showing that it exists only one triplet X = (D 0 (T, p), β(T, p), C ∞ (T, p)) ∈ℜ 3 which minimises the function: where: and C, the gas concentration in the polymer, is a solution of Fick's equation: The function φ(C, T, p) is defined according to the choice of the model, exponential or linear: (12)

STUDY OF SENSITIVITY OF THE PARAMETERS
By integrating the first equation of System (7) from the variable space x to l, one obtains: By integrating again from 0 to l, it comes: As: (15) according to the boundaries conditions of System (7), it occurs: So, Equation (14) becomes: By integrating this equation from 0 to the variable time t, one can obtain: The first term of the left-hand member can be written as: (19) Consequently: When the time t is large enough, the steady state is reached and concentration becomes independent of time. Then, the curve of Q(l, t) in function of time is a line of equation: and: The line intersection with the t-axis is given by Equation (23). It defines in general terms the time lag by: (23)

Application in the Exponential Model
When the diffusion D(C, T, p) follows the exponential law, then: In that case, the steady state is reached if and only if: The concentration of gas in steady state is then given by the expression: T p D T p yC y t y

Application in the Linear Model
By a similar procedure to that used for the exponential model, System (27) can be proven, in the linear model case (6) (see the bottom of the page).
In that case, the functions studied are: This leads to: -properties of the parameter D 0 (T, p) in the linear case: -properties of the parameter C ∞ (T, p) in the linear case:

NUMERICAL APPROXIMATION
After the step of modelling developed in the previous section, the parameters, D 0 (T, p), β(T, p) and C ∞ (T, p) are calculated. The resolution of the optimisation problem (9) is presented. At first, a finite difference approximation, using an order 2 centred scheme, is used (Smith, 1969), for space discretisation of Fick's equation. This allows to transform the search for the gas concentration in the polymer into the resolution of a system of differential equations in time. An estimate of the solution of this system allows afterward to obtain analytically the expression of the quantity of gas crossing the membrane, according to the global parameter β(T, p) C ∞ (T, p) and to the experimental data which characterise the steady state. Finally, one supplies a Gauss-Newton iterative algorithm for the optimisation of the difference between the observed values and the calculated values.

Numerical Calculation of the Quantity of Gas Crossing the Polymer
The object of this section is to present the numerical method to calculate the quantity of gas crossing the polymer. According to Equation (20): β β β 299 Oil & Gas Science and Technology -Rev. IFP, Vol. 56 (2001), No. 3 A numerical approximation of this quantity first requires to estimate the integral:

D T p l T p T p C T p T p C T p T p C T p T p l b T p T p C T p T p C T p
(30) By using the change of variables: (31) System (7) is equivalent to: By using a second-order finite-difference method, with constant time steps, the resolution of System (32) is equivalent to the resolution of the following common differential system: such as: where: C h i (t) represents an estimate of (x i , t) C 0 is a vector of ℜ N+1 verifying C 0 (1) = 1 and C 0 (i) = 0, i = 1, ..., N + 1 ϕ is the function defined by: • Exponential case: • Linear case: where: For the time discretisation of System (33), the Simulink libraries (User's Guide Matlab 5.3, 1999) allows us to use the Euler implicit method with fixed time step. By applying a trapezoidal rule, one obtains: where C k i is an estimation of . The numerical quantity of gas crossing the membrane can be written: • Linear case: Here the optimisation method is presented, it allows us to determine the optimal values of the parameters minimising the objective function: The resolution of Problem (38) represents an essential difficulty related to the fact that no information on the parameters D 0 (T, p), β(T, p) and C ∞ (T, p) is given. To go beyond this difficulty, it can be noticed that these three unknown parameters are related to the product β(T, p)C ∞ (T, p). Indeed, according to Section 1: Then, Equation (38) turns to a problem of optimisation of a single parameter. It consists in looking for X * = β(T, p)C ∝ (T, p) ∈ ℜ * minimising the objective function: (40) where: • Exponential case: • Linear case: In both cases, is an approached solution at time t k of the function C h (t) solution of the differential system: with the function ϕ defined by: • Exponential case: • Linear case: To study Problem (40), the least-square non-linear method is used. It consists in: Find X * ∈ ℜ * such as: (42) with: The resolution of Equation (42) is performed by using the Gauss-Newton iterative method (Dennis et al., 1983;Powel, 1970).

RESULTS OF OPTIMISATION
In this section, the optimal values, which minimise the difference between the experimental values of the quantity of gas crossing the polymer and the corresponding numerical results, are presented. The permeation tests are relative to the (CO 2 , polyethylene (PE)) couple. The experimental data associated to these tests are shown in Table 1.
The deviation of the results from the experimental values is represented by a percentage of error, noted E r , in the following way: with: E r percentage of error for a transport coefficient, D 0 (T, p), β(T, p) or C ∞ (T, p) σ standard deviation Coeff mean value of the coefficient at given temperature and pressure. Tables 2 and 3 allow to recapitulate the optimal values D 0 (T, p), β(T, p) and C ∞ (T, p) obtained from the optimisation method developed, as well as the mean values of the percentage of error.
The coefficients of solubility (cm 3 (STP)/cm 3 ·MPa), of diffusion (cm 2 /s) and of permeability (cm 3 (STP)·cm/ cm 2 ·s·MPa) are then obtained by the following equations (Crank, 1968): According to the model used, linear or exponential, the values of these coefficients are summarised in Tables 4 and 5.  By using the values of the percentage of error, in Tables 2 to 5, one can conclude that the exponential model represents better the permeation tests in this case than the linear model.
It is interesting to clarify that the results of the calculated optimal values are independent of the initialisation of the parameters. This shows the efficiency of the optimisation method developed. In Table 6, the comparison between the values of the permeability calculated numerically and measured is given. Mean values are very close, with similar percentages of error.

CONCLUSION
In this work, a mathematical and numerical method is built, to characterise in a unique way the transport coefficients of a gas in a polymer. This method supposes that the diffusion follows Fick's unidirectional law. It is essentially based on a study of sensitivity, which allows to obtain results in agreement with the measures, in the steady state, as well as in the transient state.
The numerical method is based on a not classical optimisation technique, which allows to converge in a unique way to the optimal values, independently on the initialisation of the parameters.
It is possible to exploit this method of optimisation to determine general models which take into account not only the influence of the concentration but also that of the temperature and pressure (Benali et al., 2001).