Moldi(Tm): a Fluid Permeation Model to Calculate the Annulus Composition in Flexible Pipes

Résumé — MOLDI™ : un modèle de perméation pour calculer la composition de fluides dans l’annulaire de conduites flexibles — Les pipelines flexibles employés dans la production d’hydrocar-bures offshore utilisent des gaines de polymères pour assurer l’étanchéité interne et externe. Ces couches de polymères laissent cependant diffuser certains gaz qui peuvent fragiliser, par corrosion, les couches métalliques situées entre les deux gaines et affecter la durée de vie d’une conduite flexible. Cette fragili-sation ne s’opère qu’en présence d’eau liquide, qui doit donc aussi être prédite. Cette publication décrit la modélisation par éléments finis des phénomènes de diffusion de fluides à travers les conduites flexibles. Le gradient de température dans la structure, élément important dans les vitesses de diffusion, est pris en compte. Cette Abstract — MOLDI™: A Fluid Permeation Model to Calculate the Annulus Composition in Flexible Pipes — In the offshore petroleum industry, flexible pipes, used as risers or flowlines, are made with internal and external polymer sheaths that ensure internal fluid and outer seawater integrity. These polymers exhibit permeability towards some gases which can create potential damage mechanisms (CO 2 corrosion and hydrogen attack) and reduce the lifetime of the steel layers located between the internal and the external plastic sheath. These damage mechanisms are associated with water condensation, which must therefore be predicted. The paper describes a finite element model named used to predict diffusion of gases through layers of flexible pipe versus time. The temperature gradient in the pipe is taken into account, as permeation rates depend on temperature. It is coupled with a thermodynamic flash calculation that gives vapor and liquid phases, compositions and volumes in the annulus. Associated mathematical models and resolution methods are examined, as well as the coupling of diffusion and phase prediction. MOLDI™ was programmed under Matlab, with a friendly graphical user interface that helps input and output of data and results. This paper illustrates the possibilities of this software.


INTRODUCTION
In the offshore petroleum industry, flexible pipes, used as risers or flowlines, are made with internal and external polymer sheaths that ensure internal fluid and outer seawater integrity (Figs. 1 and 2).These polymers exhibit permeability towards some gases, especially towards methane (CH 4 ), hydrogen sulfide (H 2 S), carbon dioxide (CO 2 ), and water (H 2 O) that can be present in the bore.They can create potential damage mechanisms (CO 2 corrosion and hydrogen attack) and reduce the lifetime of the steel layers present in the annulus, defined as the space located between the internal and the external plastic sheath.These damage mechanisms are associated with water condensation, which must therefore be predicted.
This paper describes a finite element model MOLDI™ used to predict diffusion of gases through layers of flexible pipe versus time.The temperature gradient in the pipe is taken into account, as permeation rates depend on temperature.All types of pipe layers can be modeled, and the shielding effect of metallic layers on polymer sheaths is considered.It is coupled with a thermodynamic flash calculation that gives vapor and liquid phases, compositions and volumes in the annulus.This flash calculation is appropriate for water condensation.Many of the simplified assumptions of a previous paper (Frost and Buchner, 1996) on this subject are here removed.
However, to reduce the calculation time, a simplified 2D, axisymmetrical geometry is considered in the finite element model.This geometry, with related assumptions and equations of heat and mass diffusion, is first explained.The prediction of gas and liquid phases in the annulus is then presented.Associated mathematical models and resolution methods are examined, as well as the coupling of diffusion and phase prediction.
Around those models, MOLDI™ uses five databases of material characteristics, fluid properties, diffusion data, field conditions and solutions.It is programmed under Matlab, with a friendly graphical user interface that helps input and output of data and results.Some details on databases and results illustrate the possibilities of this software.
Results can then be used as an input for a corrosion model, to estimate the life of the flexible pipe.

MODELING
Fundamental parameters and equations used to obtain the annulus composition are defined in the following.
The flexible pipe geometry is first described.It is shown how the three-dimensional structure of the flexible pipe can be transformed into a two-dimensional one.
Hypotheses of the models, mathematical formulations of heat and mass transport, and phase prediction are then presented, with associated boundary and initial conditions.

Geometrical Model
As shown in Figure 1, a flexible pipe used as a riser is a complex three-dimensional structure typically made of the following layers (from inside to outside): -a carcass, which is a metallic construction used to prevent, totally or partially, collapse of the pipe internal sheath due for example to external pressure;  -an internal pressure sheath, which is an extruded polymer layer that ensures internal fluid containment; -a pressure vault, which is a structural layer with a laying angle close to 90°, to increase the resistance of the pipe to internal and external pressure; -tensile armor layers, which are structural layers with a laying angle typically between 25°and 55°, and consist of helically wound wires.They are used to sustain tensile loads and internal pressure.Tensile armor layers are counter-wound in pairs; -anti-wear layers, which are non-metallic extruded thermoplastic sheaths or tape wrappings, used to minimize wear between structural layers; -an outer sheath, which is an extruded polymer layer that protects the pipe against penetration of seawater and external damages.
The annulus is defined as the space located between the internal and the external plastic sheath.
To reduce the calculation time, a simplified 2D, axisymmetrical (r, z) geometry (r: radius, z: distance along the pipe axis) is considered in the model for mass diffusion.This is obvious for thermoplastic cylindrical sheaths.Helices (pressure vault, tensile armors layers, anti-wear layers) are transformed into rings as shown in Figure 3.To keep the shielding effect of helices on sheaths in the 2D model, the width b' of a ring is obtained by making equal the surfaces of the ring and the helices on their lower substrate.Then, if b is the real width of a helix, and α its laying angle: (1) The gap j' between two rings is also obtained from the real gap j between two helices with: (2) Thicknesses and radii of rings and initial helices are the same.
For heat diffusion, and for a given layer made of free space and rings, it is considered that the temperature in the free space is equal to the temperature of the rings of the layer.Heat diffusion is then solved in a 1D model (no variations versus z), where temperature depends only on the radius r (Fig. 4).
Those models allow us to transform in a simple way the real geometry into an equivalent axisymmetrical one.In Figure 5, an example of plane equivalent structure with three polymer layers and two armor layers is given.
A triangular mesh of the structure (Georges, 1991) is automatically obtained from the Partial Differential Equation (PDE) Toolbox of Matlab (PDE Toolbox, 1996;Matlab Users' Guide, 1996).
Steel components are considered as impermeable to fluids that diffuse through the pipe.Then they are not meshed.Only the free volume between steel helices (or rings) is meshed.
A layer of free volume is introduced in the 2D model between two layers of rings to avoid closed volumes (and fluid accumulation) that superposition of rings could create.

Physical Parameter Hypotheses
For a given pressure p and temperature T, the fluid quantity absorbed by the surrounding at saturation is measured by using the solubility coefficient S (Comyn, 1986;Flaconnèche, 1995).The fluid flow from the surrounding is measured by using the mass diffusion coefficient D. S and D depend on the temperature following Arrhenius' law: (3) and: (4) S 0 , E s , D 0 , E D are characteristic constants for each pair of polymer/gas.They are fitted by a least squares optimization method (Gauder, 1981), on experimental data.R is the constant of perfect gas.

Heat Transport
In general, the temperature T in a material is governed by Fourier's Equation (Crank, 1968): (5) where ρ, c and k are respectively the density, the specific heat and the thermal conductivity of the material.
In the present model, it is considered that the steady state of heat diffusion is obtained faster than that for fluid diffusion: the heat transient state is neglected.Hence, the temperature depends only on the radial parameter, and the heat equation can be written as: (6) This equation is true for a homogeneous material, without discontinuities.Wherever discontinuities (with normal n) occur (between the layers), the general conditions are the following: -continuity of heat flux density: -no perfect contact: where, throughout, exponents + and -denote the right and left limits of a discontinuous function at the discontinuity, respectively.
The constant ξ ech is the heat transfer coefficient on the considered interface.It is not used in the present version of MOLDI™: a modified value of heat conductivity of steel components was adopted instead of this coefficient.Then: Boundary conditions are directly given by working conditions.If the internal and external temperatures of the flexible pipe are known and assumed constant, and if r int and r ext are respectively the interior radius and the exterior radius of the pipe, the boundary conditions are given by: The initial conditions are:

Mass Transport
If C is the fluid mass concentration, diffusion is governed by Fick's law (Crank, 1968): (13) In the present version of MOLDI™, only linear mass diffusion is considered, assuming: As the diffusion coefficient D depends only on temperature, Fick's law can be written in cylindrical co-ordinates for the 2D model of the pipe: (15) On the discontinuity surfaces (of normal n), i.e. between the materials, the following conditions must be respected: -continuity of the flux density of concentration:  -continuity of the pressure: (17) Equality ( 17) uses Henry's law (for more details see Benali et al., 2001): The boundary conditions are the concentrations on the internal and external walls.They are determined by assuming that the skins of the sheaths are saturated by the fluid at the environment pressure (equal in this case to the partial pressure of the considered fluid), on the external and internal radii of the pipe, so that: Everywhere else, on the boundary ∂Ω of the domain, the flux is null: Initially, the concentration is null in the domain except on the external and internal walls of the pipe: The previous mass transport equations were established for a single component.In the case of mixing with n f components, the mass transport problem is treated independently for each component, giving

Prediction of Gas and Liquid Phases in the Annulus
The vapor-liquid equilibrium of n f components is determined in the annulus at a given temperature T and pressure p (or volume V annulus ), with the mole number n i (1 ≤ i ≤ n f ) of each component i obtained from the diffusion calculation.
The Søreide and Whitson (1992) modification of Peng-Robinson Equation Of State (EOS) is used to predict aqueous and vapor compositions.The following set of non-linear equations has to be solved: -vapor and liquid EOS, with state = V or L: -equilibrium condition, expressed as identity of vapor and liquid fugacities: -mass conservation: -stoichiometric equations: (26) -volume conservation if the volume is known: where: The fugacity f i state is calculated using: where: (29) The equation parameters a state and b state are determined from pure component data (critical pressure, critical temperature and acentric factor, all available in the component database) and mixing rules.Binary interaction parameters for the a parameter are defined for each couple of components that can be present in the annulus (Carroll and Mather, 1995;Valderrama and Obeid-Ur-Rehman, 1988).Søreide and Whitson (1992) propose a modification of the pure component a(T) parameter for water and suggest using different binary interaction parameters in the two possible physical states.
The molar volume of water is corrected using a volume translation (Péneloux et al., 1982), to improve the liquid (aqueous) phase volume prediction.

Analytical Expression of the Temperature
Heat diffusion is governed by Equation ( 6).The solution of this equation is given by: where a i and b i are constant on each layer i (1 ≤ i ≤ n layer , with n layer the number of layers) located between the radii r = r i and r = r i + 1 .The values of the coefficients a i and b i are determined by using the boundary conditions on the flexible pipe and by using the interface conditions: -Equation ( 30) and perfect contact conditions at the interfaces (Eq.( 9)) give: -heat flux density continuity (Eq.( 7)) gives: where χ is a constant.Then: and for i = 2 to n layer .

Gas and Liquid Phases Prediction in the Annulus (Flash Calculation)
Pressure-temperature calculations are performed when the annulus pressure is given or when venting occurs.In other cases, volume-temperature calculations are considered, and the annulus pressure is a result of the equilibrium calculation.
The lowest temperature in the annulus T annulus , determined from the analytical expression of the temperature, is used as the temperature for phase prediction.
Successive substitutions (Gramajo, 1986), followed by Newton's method with quadratic convergence, are used to obtain the solution.To illustrate the successive substitutions, a pressure-temperature calculation is considered.Unknowns are the mole numbers n i state (1 ≤ i ≤ n f ) of each component i in each state (with state = V or L).After a given initialization of these unknowns, they are successively determined: , a state and b state from the mixing rules, v L and v V from EOS (23), f i V and f i L from Equation (28).New values x i new of x i are then determined from: Then, new values of n i state are deduced from: The method is repeated until |f i L -f i V | < ε, where ε is the desired precision.It can be accelerated close to the solution by switching to Newton's method.

Finite Element Method for Numerical Solution of the Concentration
The finite element method (Bank, 1990) available in the PDE Toolbox (1996) of Matlab, imposes a P 1 discretisation.This means that for solving the system given by Equations ( 15)-( 17), the unknowns must be continuous over the full domain.
However, the concentrations depend on the solubility in each medium, and therefore, are not continuous.In order to overcome this difficulty, a variable change is performed according to Equation (18).Then, the resolution of the variational problem (Johnson, 1987) defined by Fick's law is equivalent to solving it by using Henry's law (Eq.( 18)), where the pressure p becomes the unknown that is continuous over the full domain.

Variational Formulation
Let Ω = Ω 1 ∪ Ω 2 ∪ … ∪ Ω N be the domain of the flexible pipe of length H (Fig. 5) that is investigated, ∂Ω its boundary, and its partition p eq (M, t) is equal to the pressure in the free volume between rings, and is equal to the equivalent pressure in the plastic sheaths.
Let P D : Γ D × R + → R and P 0 such that: One can define also a function p by: ∀ M ∈ Ω, ∀ t, p(M, t) = P eq (M, t) -P 0 (M) (33) A variational problem in pressure can be formulated by: Find p regular functions such that: (34) Let E t the set of functions enough regular such that: To establish the variational formulation, the differential Equation ( 34) is multiplied by a function p t ∈ E t and the equality obtained on each sub-domain Ω i is integrated by parts.

Discretisation
The discretisation chosen in the PDE Toolbox (1996) of Matlab is a P 1 type element.Let: -Ω h = {K l } 1 ≤ l ≤ n , a triangulation of Ω; -E tH , a discrete sub-space of E t , for the functions p tH linear on each K l and null on Γ D ; -{a i } 1 ≤ i ≤ n , the set of triangulation nodes; -and {φ i } 1 ≤ i ≤ n , functions connecting the nodes, such that ∀ i, φ i is linear by element and φ i (a j ) = δ ij where δ ij is the Kronecker tensor.Choosing p tH (M, t) = φ i (M), and decomposing p H on the φ i 's base, the ordinary differential equations are then obtained: (41) with: The resolution of this system is then done by a classical method (Johnson and Erikson, 1991).

Iterative Time Coupled Problems
At a given time step t, the diffusion problem is first solved for each component considered in the problem.This gives each component concentration C i t (i = 1 to n f ) in the annulus free volume, from which the component mole number n i t is deduced as follows: (44) where S t i -∆t and p t i -∆t are respectively the component solubility and partial vapor pressure at the previous step time.
The thermodynamic equilibrium is then solved and gives the component fugacity f t i and partial vapor pressures p t i , which are used as annulus initial condition of the next step (t + ∆t ) diffusion problem.The solubility S t i is determined as: (45)

Venting
The closed annulus can be vented if the annulus pressure p is greater than a venting pressure P vent .Component mole numbers n i t are modified to obtain an equilibrium at the venting closure pressure P vent closure and the annulus free volume V annulus .P vent and P vent closure are defined by the user as inputs of the problem.
The following method is used to model the venting: -the solution of the flash calculation with temperature T annulus , pressure P vent closure and mole numbers n i t is determined.It gives an artificial volume V of the mixing, made of a liquid volume V L and/or a gas volume V V ; -the new mole numbers are corrected to obtain a volume V annulus (or a variation of the volume ∆V = V -V annulus ) from one of the following three cases: • if V L = 0 or V V = 0 (only one phase), then: ).
, and V V > ∆V, then: where n i t, V is the mole number of the ith component in the vapor phase; • if V L ≠ 0 or V V ≠ 0 (two phases), and V V < ∆V, then: where n i t, L is the mole number of the ith component in the liquid phase.
One can note that, with this method, the vapor phase is first vented.

Fugacity and Pressure
The pressure p can easily be replaced in MOLDI™ (especially in Henry's law and the mass transport equation) with the fugacity f of the component, which is equal in the vapor and liquid phases.Then, the boundary conditions ( 19) and (20) become: 19b) and: In this case, of course, all solubility coefficients S that are used in MOLDI™ must have been determined as functions of the fugacity of the components.

Programming
The general algorithm of the heat-mass diffusion and thermodynamic equilibrium, that constitutes the kernel of MOLDI™, is given in Figure 7.There are two main loops: one loop for each time step, and one loop on mass transport calculations for each component, followed by a gas-liquid (flash) calculation.
The kernel is programmed with Matlab language Version 5.0, using the PDE Toolbox (1996) of Matlab for finite element calculations.Note that some modules of PDE were modified for this application, especially because the problem was solved by using the pressure as unknown.

VALIDATION
MOLDI™ was validated on several theoretical cases: -for heat transport, the analytical solution was compared to MOLDI™ solution on multi-layers geometry; -for mass transport, MOLDI™ has been used for reproducing permeation measurements performed at IFP Figure 8 Flux Q (in moles per cm 2 ) of CO 2 through PVDF, as measured and calculated with MOLDI™.(Flaconnèche, 1995).In this case, a very large pipe radius was artificially used to model the diffusion across a plane membrane, in isothermal conditions.The comparison between MOLDI™ results and experimental data (for example, CH 4 diffusion in a plastified PA11 membrane, CO 2 gas permeation through PVDF ) was encouraging.An example is shown in Figure 8.It presents the flux Q (in moles per cm 2 ) of CO 2 through a plane circular membrane of PVDF (thickness: 0.001150 m, radius: 0.02 m) as a function of time.The constant temperature is 100°C and the pressure of CO 2 is 2.5 MPa. Figure 8 shows good correlation between experimental and numerical results; Component database.-a numerical validation of the phase equilibrium calculation was performed using an IFP existing software.
Coflexip Stena Offshore and IFP are also performing validation tests on simplified and real flexible pipes.Simplified pipes consist of an internal plastic sheath with a pressure vault, both inserted in a steel cylinder.These validations will be presented in a future paper.

DATABASES DESCRIPTION
In this section, the five MOLDI™ databases supplying the heart of the geometrical, diffusion and thermal-dynamical model are presented.

The Component Database
It contains the characteristics of the components considered in the diffusion problem.Each component is described by its: -chemical formula; -molecular weight (kg/mol); -critical temperature (T c ) (K and °C); -critical pressure (P c ) (MPa); -acentric factor.
Figure 9 shows a window from this database.

The Material Database
It contains the characteristics of the materials used in flexible pipes.
Each material is described by its type (polymer or steel) and thermal conductivity (W/m•K).

The Component/Material Database
It contains the data needed for permeation calculations, i.e. the diffusion (D) and solubility (S) coefficients for each couple component/material.Figure 10 shows a window from this database.

The Flexible Pipe Database
The flexible pipe database allows to create, copy or modify a flexible structure for permeation calculation.
Only mono-annular structures are allowed.The internal carcass is not considered as a layer in MOLDI™.The following data are used to create a structure: -free volume of the annulus per meter of pipe (cm 3 /m); -pipe internal diameter (that is the internal diameter of the internal plastic sheath); -a description of each layer.Two types of layers are considered: • continuous layers (plastic sheaths) for which the material and thickness of the layer must be given; • cross-layers (steel wires, tapes and anti-wear tapes) for which the material, the thickness, and the description of the layer must be input.The width of the wire or tape, the number of wires or tapes by pitch and the angle of helices must also be given.An example of this database is given in Figure 11.

The Field Conditions Database
The field conditions database allows to define a case study made out of the combination of the following items: -a flexible pipe; -a given number of components that will permeate through the pipe; -their partial pressure in the pipe; -internal and external pressures and temperatures; -a simulation time; -a venting pressure.
The flexible pipe and the components must have been previously defined in the flexible pipe and components databases.
Note that the calculation performed with a flexible pipe can be quite long, depending on the number and the type of the layers in this pipe.To reduce the calculation time, MOLDI™ offers the possibility to generate automatically a simplified structure from the original flexible pipe.This simplified structure is defined with an internal plastic sheath, a steel vault and an external plastic sheath as follows: -from the original pipe, the internal and external plastic sheaths, the free volume of the annulus, and the shielding effect of the steel vault on the internal plastic sheath are kept; -the steel vault thickness becomes the thickness of the original annulus; -the steel vault thermal conductivity is modified to obtain similar temperature distributions in the internal and external plastic sheaths in the original and simplified pipes.The model is reduced to half of a ring and gap of the pressure vault, as shown in Figure 12.Calculations performed with a simplified structure are slightly less precise (mainly for the transient state) but significantly faster when compared to a full pipe model.

Graphical User Interface of MOLDI™
The graphical user interface of MOLDI™ uses all the functionalities of Matlab graphical user interface.All inputs in databases and outputs of solutions are introduced or displayed with user-friendly windows.
The results concern all the calculated variables, and a wide range of visualization possibilities are offered.Some examples of visualization are given in Figures 13-17.They concern a flexible pipe where four components (methane, carbon dioxide, water and hydrogen sulfide) are diffusing from the internal sheath.Examples of visualization show: -the mesh used for calculation (length versus radius).One can distinguish the internal plastic sheath, the annulus (made of two armor layers and two polymer anti-wear tapes) and the external plastic sheath (Fig. 13); -the temperature profile in the structure (2D) versus radius, in °C (Fig. 14); -the annulus pressure, a 2D graph, shows the variation of total pressure (bar) in the annulus versus time and, partial pressure of each component (Fig. 15).The partial pressure of each component can be displayed separately.For example, Figure 16 shows the water partial pressure.Drops of pressure are due to venting that opens at 10 bar (1 MPa) and closes at 7 bar (0.7 MPa).The partial pressure of water increases rapidly when there is only a gas phase, until the liquid phase appears.From that point on, its partial pressure remains stationary; -the volume variation of the condensed phase in the annulus versus time (Fig. 17).This volume is given in cm 3 per meter of flexible pipe.In the example, the condensed volume increases with time, as the steady state is not yet reached.
Several 3D representations are also possible: pressures, concentrations and fluxes versus length and radius in the whole structure or in a given layer, for a given time.Figure 18 is an example of this 3D display: the total pressure (in bar) is shown at the end of calculation.

General Chart of the Software MOLDI™
The general chart of the software (Fig. 19) shows the interactions between databases.

CONCLUSION AND PERSPECTIVES
A 2D finite element model, MOLDI™, is developed in order to calculate the annulus composition of a flexible pipe versus time, including the following aspects: -permeation of bore components through the thermoplastic layers, namely CH 4 , CO 2 , H 2 S and H 2 O; -calculation of the temperature gradient through the structure and variation of permeation coefficients with temperature; -shielding effect of steel vault on the internal plastic sheath; -thermodynamic phase equilibrium in the annulus allowing water condensation calculation; -gas release out of the annulus through gas-release valves.
Programmed under Matlab, with a user-friendly graphical user interface, the software MOLDI™ offers a flexible environment for engineering and research, and rigorous numerical methods describing mass, heat diffusion and annular composition in flexible pipes.MOLDI™ is based on interactive, rich databases for offshore specialists, as shown in the paper.
In future, comparison between pressure and fugacity alternatives versus experimental results will be performed.
MOLDI™ represents a significant improvement in flexible pipe design.Figure 19 A general chart of the software.

Figure 1 Flexible
Figure 1 Flexible pipe structure.
Figure 2 Schematic permeation in the pipe.
Figure 5Example of the axisymmetrical geometrical model.

Figure 6
Figure 6Boundary conditions for temperature.
the part of the boundary where concentration is imposed (Dirichlet condition), Γ N the part of the boundary where flux concentration is imposed (Newman condition).{Γ D , Γ N } is a partition of ∂Ω.S(M) and D(M) are the solubility and the diffusivity coefficient at point M at each time.The component concentration at M is defined by:

Figure 7
Figure 7Chart organization of the kernel MOLDI™ working.
Figure 11 Flexible pipe database.

Figure 12
Figure 12Simplified structure of the flexible pipe.

Figure 14
Figure 14Temperature profile in the structure.

Figure 17 Condensed
Figure 17Condensed water volume in the annulus versus time.
Figure183D plot of pressure in the pipe, for a given time (1 bar = 0.1 MPa).