Frédéric Jacquemin

Alain Vautrin

Modelling of the moisture concentration field due to cyclical hygrothermal conditions in thick laminated pipes

Introduction

The article deals with the calculation of the moisture concentration in thick laminated pipes submitted to temperature and relative humidity cycles of the same time period τ . To approach the durability of composite structures submitted to this type of hygrothermal loading, the starting point is to determine the moisture concentration field.

It is first assumed that the kinetic of moisture sorption phenomenon in every ply of the laminated cylinder follows a Fick's law whose diffusion coefficient depends on the temperature through a usual Arrhénius' law. Making the assumption that the thermal equilibrium is instantaneously reached, the temperature is uniform over the thickness of the pipe wall and it is equal to the boundary temperature at any time t . Therefore, the moisture diffusion coefficient does not depend on the spatial coordinates and is only a time dependent periodic function of period τ .

We propose a general approach to solve this problem based on the existence of a transient solution within the pipe and a fluctuating solution in cylindrical regions close to the inner and outer lateral surfaces on which the periodical conditions are applied. After some period of time, it is shown that the transient solution reaches a permanent state and the fluctuating solution becomes periodic with a period τ . The extent of the fluctuating solution depends on the diffusion coefficient, the temperature and the period of the cycles. The transient solution is obtained by using an analytical method and the fluctuating solution is determined through a finite difference scheme.

We reveal the particular effect of cyclic hygrothermal loading on the moisture concentration through examples for thick laminated pipes whose plies have identical (homogeneous case) or different (heterogeneous case) moisture diffusion coefficients.

Proof of the existence of two domains within a cylinder

Let us first address the very simple case of an infinite solid cylinder of radius b to describe the evolution of the moisture concentration due to cyclical hygrothermal conditions of period τ . The target here is to show that the oscillations induced by cyclical conditions become negligible at a certain distance e 0 from the lateral surface.

The diffusion coefficient D(t) depends on the temperature through a usual Arrhénius' law. We make the reasonable assumption of instantaneous thermal equilibrium; finally D(t) is only a uniform spatial function of time period τ .

The moisture concentration field c(r, t) is the solution of the equation set below, with Fick's equation ( 1) and boundary and initial conditions (2):

∂c ∂t = D(t) ∂ 2 c ∂r 2 + 1 r ∂c ∂r , 0 r < b, t > 0, (1) c(b, t) = c 0 + c 1 sin(ωt), c(r, 0) = c 0 , ( 2
)
where c 0 is the initial uniform moisture concentration over the cylinder and c 1 the amplitude of the cyclical moisture concentration on the surface. c 0 and c 1 are constants.

The theoretical problem has now to be changed into a new one where the diffusion coefficient is constant in order to obtain a close form solution. Therefore, we shall consider the average of the diffusion coefficient over a period and in this purpose let us introduce the following change of variable (3):

u(t) = t 0 D(q) dq τ 0 D(q) dq -1 = D(t) D(τ ) ⇒ du = D(t) D(τ ) dt. ( 3
)
The variable u is dimensionless and stands for the number of cycles. Therefore, the system to be solved becomes:

∂c ∂u = D(τ ) ∂ 2 c ∂r 2 + 1 r ∂c ∂r , r b, u > 0, (4) c(b, u) = c 0 + c 1 sin(ωu), c(r, 0) = c 0 . ( 5
)
The solution satisfying the system (4) and ( 5) [START_REF] Carslaw | Conduction of Heat in Solids[END_REF] is:

c(r, u) = c 0 + f (r, u) + g(r, u) (6) with f (r, u) = c 1 M 0 (ω r) M 0 (ω b) sin ωu + θ 0 (ω r) -θ 0 (ω b) , ω = ω D(τ ) 1/2 , g(r, u) = 2D(τ )c 1 b ∞ n=1 exp -D(τ )α 2 n u α n ωJ 0 (rα n ) D(τ ) 2 α 4 n + ω 2 J 1 (bα n ) . ( 7
)
The final solution is composed of two terms. First, a permanent and uniform solution c 0 and second of fluctuating solution f (r, u) + g(r, u), which will tend towards a periodic solution f (r, u) when time increases. Therefore, after some time we shall obtain:

c(r, u) = c 0 + f (r, u). ( 8
)
The functions M 0 (r) and θ 0 (r) are defined with Bessel's series I 0 (r) [START_REF] Watson | A Treatise on the Theory of Bessel Functions[END_REF]: r) .

I 0 ri 1/2 = ber(r) + i bei(r) = M 0 (r)e iθ 0 (
(9)

The amplitude of the periodic solution versus ω and for b = 10 mm is depicted in Fig. 1. It decreases from the surface and becomes practically null at a distance e 0 . Thus it is easy to define a practical criterion such as at a distance e 0 from the surface the amplitude of the periodic solution can be considered as negligible: This criterion (10), taking ε = 0.01, leads to the distance e 0 which has been previously estimated by [START_REF] Verchery | Moisture Diffusion in Polymer Matrix Composites with Cyclic Environmental Conditions[END_REF] in the case of a semi-infinite plate:

M 0 (ω (b -e 0 )) M 0 (ω b) ε. (10)
e 0 = 2 π τ 0 D(t) dt. ( 11
)
Obviously, e 0 depends on ω (Fig. 1), that is on the diffusion coefficient, the temperature, through the Arrhénius' law, and the period τ of the cycles. Finally, we have proved that after some period of time the solution comprises a permanent term, operating over the interior of the cylinder [0, be 0 ], and a periodic one defined within [be 0 , b].

In the following, taking into account this result, the solution for a thick laminated pipe will be sought assuming that the two following regions exist:

-a region within the pipe wall where are operating first the transient solution and then the permanent solution; in the previous case of the solid cylinder c 0 stands for both the transient and the permanent solutions; -a region very close to the lateral surfaces where prevails the fluctuating solution which will tend towards the periodic solution when time increases.

We will prove in the following that e 0 also determines the boundary between the two latter regions, where the fluctuating solution and the transient solution respectively are prevailing.

Pipe submitted to cyclic hygrothermal conditions

Position of the problem

Let us now consider a laminated pipe whose outer and inner radii are a and b, respectively, and composed of n orthotropic plies (Fig. 2). Any ply i of it is a cylinder whose inner and outer radii are r i and r i+1 , respectively. The present approach will be conducted assuming the plies are homogeneous and characterised by homogenised effective properties, whatever the material microstructure is.

The moisture concentration in ply i is solution of the following set of Eqs. ( 12) and ( 13): The diffusion coefficient D i (t) of ply i and the boundary concentrations are periodic functions of period τ . The first equality:

∂c i ∂t = D i (t) ∂ 2 c i ∂r 2 + 1 r ∂c i ∂r , a < r < b, t > 0, i = 1 to n, ( 12
) c i (r i , t) = α i+1 c i+1 (r i , t), D i (t) ∂c i (r i , t) ∂r = D i+1 (t) ∂c i+1 (r i , t) ∂r , c(a, t) = c a (t) and c(b, t) = c b (t), c(r, 0) = 0. (13)
c i (r i , t) = α i+1 c i+1 (r i , t) (14)
results from the continuity of the chemical potential at interface i between materials i and i + 1 [START_REF] Weitsman | Moisture in composites: sorption and damage[END_REF]. The form of the potential is not generally known and it is necessary to make a reasonable assumption to perform the simulation. It is classically assumed for glassy polymers that the concentrations c i (r i , t) and c i+1 (r i , t) are simply proportional, α i+1 being a constant depending on the two materials [START_REF] Stastna | Transport Properties in Polymers[END_REF]. This assumption results from the fact the levels of the concentrations c i (r i , t) and c i+1 (r i , t) remain sufficiently small during the absorption process. Since the relation is valid at any time, it is valid at saturation as well, therefore we obtain:

c i (r i ) sat = α i+1 c i+1 (r i ) sat . ( 15
)
It means that at saturation of the laminate, every ply will be saturated and no more moisture transport will take place inside of the medium, that is in agreement with the physical sense.

Several difficulties arise when solving this problem:

-the diffusion coefficients depend on the temperature and therefore are functions of time; -the diffusion coefficients vary from ply to ply (through thickness heterogeneity of the laminate); -the boundary conditions are cyclical functions of time.

To solve this problem, we use two different methods:

-an explicit finite difference method to obtain the fluctuating solution, then the periodic solution, within the cylindrical regions of thickness e 0 , near the lateral surfaces; -an analytical method to calculate the transient solution then the permanent solution within the interior of the pipe.

The transient solution is valid on the whole cylinder except close to the surfaces where the oscillations are imposed. The fluctuating solution will ensure the continuity with the transient solution at a distance e 0 from the surface.

Determination of the transient solution

In this part, we propose a close form solution of the transient solution. The initial problem is transformed into a classical problem with time independent boundary conditions and diffusion coefficients. For this reason, the transient solution which is finally obtained does not respect the periodic boundary conditions and is valid only at a distance e 0 from the surface.

First, let us consider the same change of variable as the one previously utilised for the infinite solid cylinder (3) for every ply i:

u i (t) = t 0 D i (q) dq τ 0 D i (q) dq -1 = D i (t) D i (τ ) ⇒ du i = D i (t) D i (τ ) dt. ( 16
)
The diffusion coefficients satisfy an Arrhénius' law for each ply:

D i (t) = A i exp B i T (t) , D i+1 (t) = A i+1 exp B i+1 T (t) . ( 17
)
4

When the activation energies B i and B i+1 of the diffusion process in plies i and i + 1 are equal, the reduced time u i (t) and u i+1 (t) are equal (18):

B i = B i+1 ⇔ D i (t) D i (τ ) = D i+1 (t) D i+1 (τ ) ⇔ u i (t) = u i+1 (t) ∀i. (18)
This assumption prevails for instance when the plies of the laminate are made up of the same polymer. Homogeneous cylinders obviously satisfy this requirement. In that case one single reduced time u applies to the laminate as a whole so we can use the previous change of variable (16). The problem is now rewritten using the variable u:

u i (t) = u(t) ∀i. ( 19
)
The problem ( 12) and ( 13) becomes:

∂c i ∂u = D i (τ ) ∂ 2 c i ∂r 2 + 1 r ∂c i ∂r , a < r < b, u > 0, i = 1 to n, ( 20
)
c i (r i , u) = α i+1 c i+1 (r i , u), D i (τ ) ∂c i (r i , u) ∂r = D i+1 (τ ) ∂c i+1 (r i , u) ∂r , c(a, u) = c a (u) and c(b, u) = c b (u), c(r, 0) = 0, ( 21
)
c a (u), c b (u) are cyclic with period 1. The next step consists in transforming the boundary conditions ( 21), which still depend on the reduced time u, into time independent conditions. For this purpose, we introduce the so-called mobile average concentration:

ci (r, u) = u u-1 c i (r, q) dq. ( 22
)
The mobile average concentration is solution of the following system ( 23) and ( 24):

∂ ci ∂u = D i (τ ) ∂ 2 ci ∂r 2 + 1 r ∂ ci ∂r , a < r < b, u > 0, i = 1 to n, ( 23
) ci (r i , u) = α i+1 ci+1 (r i , u), D i (τ ) ∂ ci (r i , u) ∂r = D i+1 (τ ) ∂ ci+1 (r i , u) ∂r , c(a, u) = ca with ca = 1 0 c a (q) dq, c(b, u) = cb with cb = 1 0 c b (q) dq. ( 24
)
Applying the Laplace transform to system (23) and ( 24) and by using the residue theory to get the solution in the time space, we finally obtain the mobile average concentration:

ci (r, u) = 1 ∆ s A i + B i ln r + ∞ m=1 2 exp -ω 2 m u ω m ∆ u (ω m ) A i J 0 (β i ω m r) + B i Y 0 (β i ω m r) , ( 25
)
where J 0 and Y 0 are Bessel's functions of order zero, ∆ u and ∆ s are determinants of 2n × 2n matrices

a 11 = J 0 (β 1 ω ā), a 12 = Y 0 (β 1 ω ā), a 2n2n-1 = J 0 (β n ω), a 2n2n = Y 0 (β n ω), a 2i2i-1 = J 0 (β i ω ri ), a 2i2i = Y 0 (β i ω ri ), a 2i2i+1 = -α i+1 J 0 (β i+1 ω ri ), a 2i2i+2 = -α i+1 Y 0 (β i+1 ω ri ), a 2i+12i-1 = D i (τ )β i ωJ 1 (β i ω ri ), a 2i+12i = D i (τ )β i ωY 1 (β i ω ri ), a 2i+12i+1 = -α i+1 D i+1 (τ )β i+1 ωJ 1 (β i+1 ω ri ), a 2i+12i+2 = -α i+1 D i+1 (τ )β i+1 ωY 1 (β i+1 ω ri ), for i = 1 to n -1; b 11 = 1, b 12 = ln ā, b 2n2n-1 = 1, b 2n2n = 0, b 2i2i-1 = 1, b 2i2i = ln ri , b 2i2i+1 = -α i+1 , b 2i2i+2 = -α i+1 ln ri , b 2i+12i = D i (τ )/r i , b 2i+12i+2 = -α i+1 D i+1 (τ )/r i , for i = 1 to n -1; g 1 = ca , g 2n = cb .
Finally, a recursive scheme is used to express the transient concentration as a function of the mobile average concentration. Deriving ( 22) with respect to u, we obtain:

c i (r, u) = ∂ ci ∂u (r, u) + c i (r, u -1). ( 26
)
If N is the number of cycles necessary to reach the permanent state in ply i:

c i (r, N + 1) = c i (r, N), ( 27
)
considering the expression ( 26), the average concentration in ply i for N + 1 cycles is expressed as follows ( 28):

∂ ci ∂u (r, N + 1) = 0 ⇔ ci (r, N + 1) = 1 ∆ s A i + B i ln r , ( 28
)
where ci (r, N + 1) is defined by:

ci (r, N + 1) = N+1 N c i (r, u) du = c i (r, N). ( 29
)
Therefore, the permanent solution in ply i is expressed by:

c i (r, N) = 1 ∆ s A i + B i ln r . ( 30
)
To start the recurrence, we consider the expression (26) for N cycles:

c i (r, N) = ∂ ci ∂u (r, N) + c i (r, N -1); (31)
for (Nk) cycles (31) becomes:

c i (r, N -k) = ∂ ci ∂u (r, N -k) + c i (r, N -k -1); ( 32
)
so it results from ( 31) and ( 32) that the concentration for (Nk) cycles is:

c i (r, N -k) = c i (r, N) - ∂ ci ∂u (r, N) -• • • - ∂ ci ∂u (r, N -k + 1). ( 33
)
We obtain the general form of the transient solution substituting in (33) the expressions of the mobile average concentration (25) and of the permanent concentration (30):

c i (r, N -k) = 1 ∆ s A i + B i ln r + ∞ m=1 ω 2 m 2 exp -ω 2 m N ω m ∆ u (ω m ) A i J 0 (β i ω m r) + B i Y 0 (β i ω m r) + • • • + ∞ m=1 ω 2 m 2 exp -ω 2 m (N -k + 1) ω m ∆ u (ω m ) A i J 0 (β i ω m r) + B i Y 0 (β i ω m r) . ( 34
)
Finally, we obtain the expression of the transient moisture concentration for (Nk) cycles:

c i (r, N -k) = 1 ∆ s A i + B i ln r + ∞ m=1 2ω m ∆ u (ω m ) k-1 n=0 exp -ω 2 m (N -n) A i J 0 (β i ω m r) + B i Y 0 (β i ω m r) . ( 35
)

Examples

We consider a thick walled pipe made up of 5 plies of equal thickness, the inner and outer radii being respectively 10 mm and 30 mm. Three lamination sequences are studied:

-[55, -55, 55, -55, 55], i.e. ±55 degrees cross-ply unidirectional carbon fibre reinforced polymer T300/5208, the fibre angle being defined versus the longitudinal axis; -[resin, -55, 55, -55, 55], the internal ply being composed of pure resin referenced N 5208 (liner); -[55, -55, 55, -55, resin], the pure resin ply being the outer ply.

The moisture concentration does not depend on the fibre orientations within the plies, so far it is only dependent on the transverse ply diffusion properties. So from the diffusion point of view, the first case can be regarded as homogeneous, since the same diffusion coefficient holds for all the plies, and the two others cases being heterogeneous. It should be pointed out that to compute the internal state of stress the fibre orientation must be considered [START_REF] Jacquemin | Transient thermoelastic stresses in an anisotropic laminated hollow cylinder[END_REF] therefore, the three pipes should be considered as heterogeneous from the mechanical point of view.

We consider temperature and relative humidity cycles of 4 week period (Figs. 3 and4) and select particular points to compute the moisture concentration along the analysis (Figs. 3 and4, 5). Those points correspond to deep changes either in relative humidity or temperature.

Case of diffusion -homogeneous pipes

The hygroscopic properties of the composite T300/5208 are reported in Table 1. In the simulation the activation energy of the composite is assumed to be equal to the activation energy of the neat resin N 5208 [START_REF] Loos | Moisture absorption of graphite -epoxy composition immersed in liquids and in humid air[END_REF], since the diffusion process relies on the resin properties only.

The transient moisture concentration through the thickness the [55, -55, 55, -55, 55] laminate is depicted in Fig. 5. Since the transient solution does not take into account the oscillations of the external boundary conditions, we obtain a constant value on the external surfaces whatever the number of cycles (Fig. 5). The convergence of the concentration towards a constant permanent solution for 2500 cycles is ensured and the transient solution is in agreement with a finite difference calculation to valid the approach.

Case of diffusion -heterogeneous pipes

The hygroscopic properties of the resin N 5208 [START_REF] Loos | Moisture absorption of graphite -epoxy composition immersed in liquids and in humid air[END_REF] are presented in Table 2. For the heterogeneous pipes, the continuity of the chemical potential results, as it has been discuss previously, in a discontinuity of the moisture concentration at the ply interface. This discontinuity is represented by introducing a coefficient α at the ply interface (Table 3).

Internal ply in neat resin N 5208

The transient moisture concentration through the thickness of a pipe, whose internal ply is made up of resin, is plotted in Fig. 7. We notice the discontinuity of the moisture concentration at the interface between the neat resin and the composite (Fig. 7). Due to the neat resin of the internal ply, larger concentrations than in the case of the homogeneous cylinder are achieved. Moreover, the steady state, that is when the saturation concentrations of the neat resin and the composite are reached, is more quickly achieved. The transient solution that is obtained is again in agreement with the finite difference calculation.

Fig. 8 shows the moisture concentration during cycle 50. We notice large concentration gradients at points C and E where the ambient concentration undergoes abrupt changes (Fig. 8). As the composition of the external ply and the cyclic conditions do not change, the extent of the fluctuating part through the external ply is the same as in the homogeneous pipe. Since the internal ply is made up of neat resin, the extent of the fluctuating solution is more important and e 0 is about 2.5 mm. Moreover, since the saturation concentration of the neat resin is higher than for the composite, the fluctuating concentration gradients are more pronounced through the internal ply.

External ply in neat resin N 5208

The transient moisture concentration through the thickness of a pipe whose external ply is made up of resin is plotted in Fig. 9. The external ply made up of neat resin submitted to cyclic conditions presents a saturated moisture concentration higher than the composite under similar conditions. The transient solution is still in agreement with the finite difference calculation (Fig. 9). Fig. 10 shows the moisture concentration during cycle 50. The extent of the fluctuating solution for the external resin layer is about 2.5 mm. The most important fluctuating concentration gradients take place through the external ply.

Conclusion

An analytical solution is presented that allows to determine the moisture concentration field through thick laminated composite pipes, submitted to cyclic environmental conditions. The method relies on controlled assumptions, therefore the solution can be regarded as a reference solution that can be useful to calibrate FE modelling. This analytical solution leads to the time-dependent moisture concentration of structures whose composition is heterogeneous by assuming that the ply diffusion phenomena have all the same activation. It is quite valuable to perform extensive simulations to optimise ply lay-ups. This solution is only valid at a distance e 0 from the lateral surfaces of the pipe since a fluctuating solution prevails in that regions. The fluctuating solution is obtained by using a finite difference calculation and it shows that steep gradients can arise close to the lateral surfaces and therefore cause stress concentrations and hygrothermal fatigue can result from them.

Fig. 1 .

 1 Fig. 1. Amplitude of the periodic solution.

Fig. 2 .

 2 Fig. 2. Laminated hollow cylinder.

  [a] and [b]. A i and B i are determinants of matrices deduced from [a] by respectively substituting column (2i -1) and 2i by the constant vector {g}. A i and B i are similar determinants of A m and B m but considering the matrix [b]. ∆ u (ω m ) is the derivative of ∆ u with respect to ω calculated for ω m the mth positive root of ∆ u . r, β i are defined by the relations r = r/b, β i = b 2 /D i (τ ). Furthermore, the nonzero elements of [a], [b] and {g} are:

Fig. 3 .

 3 Fig. 3. Temperature cycle.

Fig. 5 .

 5 Fig. 5. Transient moisture concentration.

Fig. 7 .

 7 Fig. 7. Transient moisture concentration.

Fig. 8 .

 8 Fig. 8. Moisture concentration at cycle 50.

Fig. 9 .

 9 Fig. 9. Transient moisture concentration.

Table 1

 1 Hygroscopic properties for T300/5208 composite

	Diffusion coefficient mm 2 /s	D(t) = 0.57 exp(-5116/T (t))
		T (Kelvin)
	Ambient moisture concentration kg/m 3	c = 0.2385 H
		H (%)

Table 2

 2 Hygroscopic properties for N 5208 neat resin

	Diffusion coefficient mm 2 /s	D(t) = 2.8 exp(-5116/T (t))
		T (Kelvin)
	Ambient moisture concentration kg/m 3	c = 0.708 H
		H (%)
	Table 3
	Coefficient of concentration discontinuity