Liangyu Chen
email: lychen@sei.ecnu.edu.cn

Svyatoslav Covanov
email: svyatoslav.covanov@loria.fr

Davood Mohajerani

Marc Moreno Maza

Big Prime Field FFT on the GPU

Keywords: I.1.2 [Symbolic and Algebraic Manipulation]: Algorithms-Algebraic algorithms, Analysis of algorithms Algorithms, Experimentation Fast Fourier transforms, Finite fields, Generalized Fermat numbers, Graphics processing units, CUDA

We consider prime fields of large characteristic, typically fitting on k machine words, where k is a power of 2. When the characteristic of these fields is restricted to a subclass of the generalized Fermat numbers, we show that arithmetic operations in such fields offer attractive performance both in terms of algebraic complexity and parallelism. In particular, these operations can be vectorized, leading to efficient implementation of fast Fourier transforms on graphics processing units.

INTRODUCTION

Prime field arithmetic plays a central role in computer algebra and supports computation in Galois fields which are essential to coding theory and cryptography algorithms.

The prime fields that are used in computer algebra systems, in particular in the implementation of modular methods, are often of small characteristic, that is, based on prime numbers that fit in a machine word. Increasing precision beyond the machine word size can be done via the Chinese Remainder Theorem (CRT) or Hensel Lemma.

However, using machine-word size, thus small, prime numbers has also serious inconveniences in certain modular methods, in particular for solving systems of non-linear equations. Indeed, in such circumstances, the so-called unlucky primes are to be avoided, see for instance [START_REF] Arnold | Modular algorithms for computing Gröbner bases[END_REF][START_REF] Dahan | Lifting techniques for triangular decompositions[END_REF].

In this paper, we consider prime fields of large characteristic, typically fitting on k machine words, where k is a power of 2. When the characteristic of these fields is restricted to a subclass of the generalized Fermat numbers, we show that arithmetic operations in such fields offer attractive performance both in terms of algebraic complexity and parallelism. In particular, these operations can be vectorized, leading to efficient implementation of fast Fourier transforms on graphics processing units.

We present algorithms for arithmetic operations in a "big" prime field Z pZ, where p is a generalized Fermat number of the form p = r k + 1 where r fits a machine-word and k is a power of 2. We report on a GPU (Graphics Processing Units) implementation of those algorithms as well as a GPU implementation of a Fast Fourier Transform (FFT) over such big prime field. Our experimental results show that 1. computing an FFT of size N , over a big prime field for p fitting on k 64-bit machine-words, and 2. computing 2k FFTs of size N , over a small prime field (that is, where the prime fits a 32-bit half-machineword) followed by a combination (i.e. CRT-like) of those FFTs are two competitive approaches in terms of running time.

Since the former approach has the advantage of reducing the occurrence of unlucky primes when applying modular methods (in particular in the area of polynomial system solving), we view this experimental observation as a promising result.

The reasons for a GPU implementation are as follows. First, the model of computations and the hardware performance provide interesting opportunities to implement big prime field arithmetic, in particular in terms of vectorization of the program code. Secondly, highly optimized FFTs over small prime fields have been implemented on GPUs by Wei Pan [START_REF] Moreno Maza | Fast polynomial arithmetic on a GPU[END_REF][START_REF] Moreno Maza | Solving bivariate polynomial systems on a GPU[END_REF] in the CUMODP library www.cumodp.org we use them in our experimental comparison.

Section 7 reports on various comparative experimentations. First, a comparison of the above two approaches implemented on GPU, exhibiting an advantage for the FFT over a big prime field. Second, a comparison between the two same approaches implemented on a single-core CPU, exhibiting an advantage for the CRT-based FFT over small prime fields. Third, from the two previous comparisons, one deduces a comparison of the FFT over a big prime field (resp. the CRT-based FFT over small prime fields) implemented on GPU and CPU, exhibiting a clear advantage for the GPU implementations. Overall, the big prime field FFT on the GPU is the best approach.

A discrete Fourier transform (DFT) over Z pZ, when p is a generalized Fermat prime, can be seen as a generalization of the FNT (Fermat number transform), which is a specific case of the NTT (number theoretic transform). However, the computation of a DFT over Z pZ implies additional considerations, which are not taken into account in the literature describing the computation of a NTT, or a FNT [START_REF] Burrus | Fast convolution using fermat number transforms with applications to digital filtering[END_REF][START_REF] Dimitrov | Generalized fermat-mersenne number theoretic transform[END_REF].

The computation of a NTT can be done via various methods used for a DFT, among them is the radix-2 Cooley-Tukey, for example. However, the final complexity depends on the way a given DFT is computed. It appears that, in the context of generalized Fermat primes, there is a better choice than the radix-2 Cooley-Tukey. The method used in the present paper is related to the article [START_REF] Covanov | Fast arithmetic for faster integer multiplication[END_REF], which is derived from Fürer's algorithm [START_REF] Fürer | Faster integer multiplication[END_REF] for the multiplication of large integers. The practicality of this latter algorithm is an open question. And, in fact, the work reported in our paper is a practical contribution responding to this open question.

The paper [START_REF] Burrus | Fast convolution using fermat number transforms with applications to digital filtering[END_REF] discusses the idea of using Fermat number transform for computing convolutions, thus working modulo numbers of the form F = 2 b + 1, where b is a power of 2. This is an effective way to avoid round-off error caused by twiddle factor multiplication in computing DFT over the field of complex numbers. The paper [START_REF] Dimitrov | Generalized fermat-mersenne number theoretic transform[END_REF] considers generalized Fermat Mersenne (GFM) prime numbers that prime of the form (q pn -1) (q -1) where, typically, q is 2 and both p and n are small. These numbers are different from the primes used in our paper, which have the form r k + 1 where r is typically a machine-word long and k is a power of 2 so that r is a 2k-th primitive root of unity, see Section 3.

COMPLEXITY ANALYSIS

Consider a prime field Z pZ and N , a power of 2, dividing p -1. Then, the finite field Z pZ admits an N -th primitive root of unity 1 ; We denote by ω such an element. Let f ∈ Z pZ[x] be a polynomial of degree at most N -1. Then, computing the DFT of f at ω via an FFT amounts to:

1. N log(N ) additions in Z pZ, 2. (N 2) log(N ) multiplications by a power of ω in Z pZ. If the bit-size of p is k machine words, then 1. each addition in Z pZ costs O(k) machine-word operations, 2. each multiplication by a power of ω costs O(M(k)) machine-word operations, where n → M(n) is a multiplication time as defined in [START_REF] Zur Gathen | Modern Computer Algebra[END_REF]. Therefore, multiplication by a power of ω becomes a bottleneck as k grows. To overcome this difficulty, we consider the following trick proposed by Martin Fürer in [START_REF] Fürer | Faster integer multiplication[END_REF][START_REF] Fürer | Faster integer multiplication[END_REF]. We assume that N = K e holds for some "small" K, say K = 32 and an integer e ≥ 2. Further, we define η = ω N K , with J = K e-1 and assume that multiplying an arbitrary element of Z pZ by η i , for any i = 0, . . . , K -1, can be done within O(k) machine-word operations. Consequently, every arithmetic operation (addition, multiplication) involved in a DFT of size K, using η as a primitive root, amounts to O(k) machine-word operations. Therefore, such DFT of size K can be performed with O(K log(K) k) machine-word operations. As we shall see in Section 3, this latter result holds whenever p is a so called generalized Fermat number.

Returning to the DFT of size N at ω and using the factorization formula of Cooley and Tukey, we have

DFT JK = (DFT J ⊗ I K )D J,K (I J ⊗ DFT K )L JK J , (1)
see Section 4. Hence, the DFT of f at ω is essentially performed by: 1. K e-1 DFT's of size K (that is, DFT's on polynomials of degree at most K -1), 1 See Section 4 for this notion.

2. N multiplications by a power of ω (coming from the diagonal matrix D J,K ) and 3. K DFT's of size K e-1 . Unrolling Formula (1) so as to replace DFT J by DFT K and the other linear operators involved (the diagonal matrix D and the permutation matrix L) one can deduce that a DFT of size N = K e reduces to:

1. e K e-1 DFT's of size K, and 2. (e -1) N multiplication by a power of ω. Recall that the assumption on the cost of a multiplication by η i , for 0 ≤ i < K, makes the cost for one DFT of size K to O(K log 2 (K) k) machine-word operations. Hence, all the DFT's of size K together amount to O(e N log 2 (K)k) machine-word operations, that is, O(N log 2 (N ) k) machineword operations. Meanwhile, the total cost of the multiplication by a power of ω is O(e N M(k)) machine-word operations, that is, O(N log K (N ) M(k)) machine-word operations. Indeed, multiplying an arbitrary element of Z pZ by an arbitrary power of ω requires O(M(k)) machine-word operations. Therefore, under our assumption, a DFT of size N at ω amounts to

O(N log 2 (N ) k + N log K (N ) M(k)) (2)
machine-word operations. When using generalized Fermat primes, we have K = 2k and the above estimate becomes

O(N log 2 (N ) k + N log k (N ) M(k)) (3)
The second term in the big-O notation dominates the first one. However, we keep both terms for reasons that will appear shortly. Without our assumption, as discussed earlier, the same DFT would run in O(N log 2 (N ) M(k)) machine-word operations. Therefore, using generalized Fermat primes brings a speedup factor of log(K) w.r.t. the direct approach using arbitrary prime numbers.

At this point, it is natural to ask what would be the cost of a comparable computation using small primes and the CRT. To be precise, let us consider the following problem. Let p1, . . . , p k pairwise different prime numbers of machine-word size and let m be their product. Assume that N divides each of p1 -1, . . . , p k -1 such that each of fields Z p1Z, . . . , Z p k Z admits an N -th primitive roots of unity, ω1, . . . , ω k . Then, ω = (ω1, . . . , ω k ) is an N -th primitive root of Z mZ. Indeed, the ring Z p1Z ⊗ ⋯ ⊗ Z p k Z is a direct product of fields. Let f ∈ Z mZ[x] be a polynomial of degree N -1. One can compute the DFT of f at ω in three steps:

1. Compute the images f1

∈ Z p1Z[x], . . . , f k ∈ Z p k Z[x] of f . 2. Compute the DFT of fi at ωi in Z piZ[x], for i = 1, . . . , k, 3.
Combine the results using CRT so as to obtain a DFT of f at ω. The first and the third above steps will run within O(N × M(k) log 2 (k)) machine-word operations meanwhile the second one amount to O(k × N log(N )) machine-word operations, yielding a total of

O(N log 2 (N ) k + N M(k) log 2 (k)) (4)
These estimates yield a running-time ratio between the two approaches of log(N ) log 2 2 (k), which suggests that for k large enough the big prime field approach may outperform the CRT-based approach. We believe that this analysis is part of the explanation for the observation that the two approaches are, in fact, competitive in practice, as we shall see in Section 7.

SPARSE RADIX GENERALIZED FERMAT NUMBERS

The n-th Fermat number, denoted by Fn, is given by Fn = 2 2 n + 1. This sequence plays an important role in number theory and, as mentioned in the introduction, in the development of asymptotically fast algorithms for integer multiplication [START_REF] Schönhage | Schnelle multiplikation großer zahlen[END_REF][START_REF] Fürer | Faster integer multiplication[END_REF].

Arithmetic operations modulo a Fermat number are simpler than modulo an arbitrary positive integer. In particular 2 is a 2 n+1 -th primitive root of unity modulo Fn. Unfortunately, F4 is the largest Fermat number which is known to be prime. Hence, when computations require the coefficient ring be a field, Fermat numbers are no longer interesting. This motivates the introduction of other family of Fermatlike numbers, see, for instance, Chapter 2 in the text book Guide to elliptic curve cryptography [START_REF] Hankerson | Guide to Elliptic Curve Cryptography[END_REF].

Numbers of the form a 2 n +b 2 n where a > 1, b ≥ 0 and n ≥ 0 are called generalized Fermat numbers. An odd prime p is a generalized Fermat number if and only if p is congruent to 1 modulo 4. The case b = 1 is of particular interest and, by analogy with the ordinary Fermat numbers, it is common to denote the generalized Fermat number a 2 n + 1 by Fn(a). So 3 is F0(2). We call a the radix of Fn(a). Note that, Landau's fourth problem asks if there are infinitely many generalized Fermat primes Fn(a) with n > 0.

In the finite ring Z Fn(a)Z, the element a is a 2 n+1 -th primitive root of unity. However, when using binary representation for integers on a computer, arithmetic operations in Z Fn(a)Z may not be as easy to perform as in Z FnZ. This motivates the following. Definition 1. We call sparse radix generalized Fermat number, any integer of the form Fn(r) where r is either 2 w + 2 u or 2 w -2 u , for some integers w > u ≥ 0. In the former case, we denote Fn(r) by F + n (w, u) and in the latter by F - n (w, u).

Table 1 lists sparse radix generalized Fermat numbers (SRGFNs) that are prime. For each such number p, we give the largest power of 2 dividing p -1, that is, the maximum length N of a vector to which a radix-K FFT algorithm where K is an appropriate power of 2. Notation 1. In the sequel, we consider p = Fn(r), a fixed SRGFN. We denote by 2 e the largest power of 2 dividing p-1 and we define k = 2 n , so that p = r k + 1 holds.

As we shall see in the sequel of this section, for any positive integer N which is a power of 2 such that N divides p -1, one can find an N -th primitive root of unity ω ∈ Z pZ such that multiplying an element x ∈ Z pZ by ω i(N 2k) for 0 ≤ i < 2k can be done in linear time w.r.t. the bit size of x. Combining this observation with an appropriate factorization of the DFT transform on N points over Z pZ, we obtain an efficient FFT algorithm over Z pZ.

Representation of Z pZ

We represent each element x ∈ Z pZ as a vector ⃗

x = (x k-1 , x k-2 , . . . , x0) of length k and with non-negative integer coefficients such that we have

x ≡ x k-1 r k-1 + x k-2 r k-2 + ⋯ + x0 mod p. ( 5
)
This representation is made unique by imposing the following constraints 1. either x k-1 = r and x k-2 = ⋯ = x1 = 0, 2. or 0 ≤ xi < r for all i = 0, . . . , (k -1). We also map x to a univariate integer polynomial fx ∈ Z[T ] defined by fx = ∑ k-1 i=0 xit i such that x ≡ fx(r) mod p. Now, given a non-negative integer x < p, we explain how the representation ⃗

x can be computed. The case x = r k is trivially handled, hence we assume x < r k . For a nonnegative integer z such that z < r 2 i holds for some positive integer i ≤ n = log 2 (k), we denote by vec(z, i) the unique sequence of 2 i non-negative integers (z 2 i -1 , . . . , z0) such that we have 0 ≤ zj < r and z = z 2 i -1 r 2 i -1 + ⋯ + z0. The sequence vec(z, i) is obtained as follows:

1

. if i = 1, we have vec(z, i) = (q, s), 2. if i > 1, then vec(z, i)
is the concatenation of vec(q, i -1) followed by vec(s, i -1), where q and s are the quotient and the remainder in the Euclidean division of z by r 2 i-1

. Clearly, vec(x, n) = ⃗ x holds. We observe that the sparse binary representation of r facilitates the Euclidean division of a non-negative integer z by r, when performed on a computer. Referring to the notations in Definition 1, let us assume that r is 2 w + 2 u , for some integers w > u ≥ 0. (The case 2 w -2 u would be handled in a similar way.) Let z high and z low be the quotient and the remainder in the Euclidean division of z by 2 w . Then, we have

z = 2 w z high + z low = r z high + z low -2 u z high . (6)
Let s = z low -2 u z high and q = z high . Three cases arise: (S1) if 0 ≤ s < r, then q and s are the quotient and remainder of z by r, (S2) if r ≤ s, then we perform the Euclidean division of s by r and deduce the desired quotient and remainder, (S3) if s < 0, then (q, s) is replaced by (q + 1, s + r) and we go back to Step (S1). Since the binary representations of r 2 can still be regarded as sparse, a similar procedure can be done for the Euclidean division of a non-negative integer z by r 2 . For higher powers of r, we believe that Montgomery algorithm is the way to go, though this remains to be explored.

Finding primitive roots of unity in Z pZ

Notation 2. Let N be a power of 2, say 2 , dividing p-1 and let g ∈ Z pZ be an N -th primitive root of unity.

Recall that such an N -th primitive root of unity can be obtained by a simple probabilistic procedure. Write p = qN + 1. Pick a random α ∈ Z pZ and let ω = α q . Little Fermat theorem implies that either ω N 2 = 1 or ω N 2 = -1 holds. In the latter case, ω is an N -th primitive root of unity. In the former, another random α ∈ Z pZ should be considered. In our various software implementation of finite field arithmetic [START_REF] Li | The modpn library: Bringing fast polynomial arithmetic into maple[END_REF][START_REF] Chen | The basic polynomial algebra subprograms[END_REF][START_REF] Haque | Dense arithmetic over finite fields with the CU-MODP library[END_REF], this procedure finds an N -th primitive root of unity after a few tries and has never been a performance bottleneck.

In the following, we consider the problem of finding an Nth primitive root of unity ω such that ω N 2k = r holds. The intention is to speed up the portion of FFT computation that requires to multiply elements of Z pZ by powers of ω.

Proposition 1. In Z pZ, the element r is a 2k-th primitive root of unity. Moreover, the following algorithm computes an N -th primitive root of unity ω ∈ Z pZ such that we have

ω N 2k = r in Z pZ. Algorithm 1 Primitive N -th root ω ∈ Z pZ s.t. ω N 2k = r procedure PrimitiveRootAsRootOf(N, r, k, g) α ∶= g N 2k β ∶= α j ∶= 1 while β ≠ r do β ∶= αβ j ∶= j + 1 end while ω ∶= g j return (ω) end procedure
Proof Since g N 2k is a 2k-th root of unity, it is equal to r i 0 (modulo p) for some 0 ≤ i0 < 2k where i0 is odd. Let j be a non-negative integer. Observe that we have

g j2 2k = (g i g 2 k q ) 2 2k = g i2 2k = r i i 0 , (7)
where q and i are quotient and the remainder of j in the Euclidean division by 2k. By definition of g, the powers g i2 2k , for 0 ≤ i < 2k, are pairwise different. It follows from Formula (7) that the elements r i i 0 are pairwise different as well, for 0 ≤ i < 2k. Therefore, one of those latter elements is r itself. Hence, we have j1 with 0 ≤ j1 < 2k such that g j 1 N 2k = r. Then, ω = g j 1 is as desired and Algorithm 1 computes it. ◻

Addition and subtraction in Z pZ

Let x, y ∈ Z pZ represented by ⃗ x, ⃗ y, see Section 3.1 for this latter notation. Algorithm 2 computes the representation → x + y of the element (x + y) mod p. Proof At Step (1), ⃗

x and ⃗ y, regarded as vectors over Z, are added component-wise. At Steps (2) and (3), the carry, if any, is propagated. At Step (4), there is no carry beyond the leading digit z k-1 , hence (z k-1 , . . . , z0) represents x + y.

Step (5) handles the special case where x + y = p -1 holds.

Step (6) is the overflow case which is handled by subtracting 1 mod p to (z k-1 , . . . , z0), finally producing → x + y. ◻ A similar procedure computes the vector → x -y representing the element (x -y) ∈ Z pZ. Recall that we explained Algorithm 2 Computing x + y ∈ Z pZ for x, y ∈ Z pZ procedure BigPrimeFieldAddition(⃗ x, ⃗ y, r, k) 1: compute zi = xi + yi in Z, for i = 0, . . . , k -1, 2: let z k = 0, 3: for i = 0, . . . , k -1, compute the quotient qi and the remainder si in the Euclidean division of zi by r, then replace (zi+1, zi) by (zi+1 + qi, si), 4: if z k = 0 then return (z k-1 , . . . , z0), 5: if z k = 1 and z k-1 = ⋯ = z0 = 0, then let z k-1 = r and return (z k-1 , . . . , z0), 6: let i0 be the smallest index, 0 ≤ i0 ≤ k, such that zi 0 ≠ 0, then let zi 0 = zi 0 -1, let z0 = ⋯ = zi 0 -1 = r -1 and return (z k-1 , . . . , z0). end procedure in Section 3.1 how to perform the Euclidean divisions at Step (S3) in a way that exploits the sparsity of the binary representation of r.

In practice, the binary representation of the radix r fits a machine word, see Table 1. Consequently, so does each of the "digit" in the representation ⃗

x of every element x ∈ Z pZ. This allows us to exploit machine arithmetic in a sharper way. In particular, the Euclidean divisions at Step (S3) can be further optimized.

Multiplication by a power of r in Z pZ

Before considering the multiplication of two arbitrary elements x, y ∈ Z pZ, we assume that one of them, say y, is a power of r, say y = r i for some 0 < i < 2k. Note that the cases i = 0 = 2k are trivial. Indeed, recall that r is a 2k-th primitive root of unity in Z pZ. In particular, r k = -1 in Z pZ. Hence, for 0 < i < k, we have r k+i = -r i in Z pZ. Thus, let us consider first the case where 0 < i < k holds. We also assume 0 ≤ x < r k holds in Z, since the case x = r k is easy to handle. From Equation (5) we have:

xr i ≡ (x k-1 r k-1+i + ⋯ + x0 r i ) mod p ≡ j=k-1 ∑ j=0 xjr j+i mod p ≡ h=k-1+i ∑ h=i x h-i r h mod p ≡ ( h=k-1 ∑ h=i x h-i r h - h=k-1+i ∑ h=k x h-i r h-k ) mod p
The case k < i < 2k can be handled similarly. Also, in the case i = k we have xr i = -x in Z pZ. It follows, that for all 0 < i < 2k, computing the product x r i simply reduces to computing a subtraction. This fact, combined with Proposition 1, motivates the development of FFT algorithms over Z pZ.

Multiplication in Z pZ

Let again x, y ∈ Z pZ represented by ⃗

x, ⃗ y and consider the univariate polynomials fx, fy ∈ Z[T ] associated with x, y; see Section 3.1 for this notation. To compute the product x y in Z pZ, we proceed as follows.

For large values of k, fxfy mod T k + 1 in Z[T ] can be computed by asymptotically fast algorithms (see the paper [START_REF] Chen | Parallel integer polynomial multiplication[END_REF][START_REF] Covanov | Fast arithmetic for faster integer multiplication[END_REF]). However, for small values of k (say k ≤ 8), using plain multiplication is reasonable. uiT i , we observe that for all 0 ≤ i ≤ k -1 we have 0 ≤ ui ≤ kr 2 and compute a representation → ui of ui in Z pZ using the method explained in Section 3.1. 3: We compute uir i in Z pZ using the method of Section 3.4.

4: Finally, we compute the sum

k-1 ∑ i=0
uir i in Z pZ using Algorithm 2. end procedure

FFT BASICS

We review the Discrete Fourier Transform over a finite field, and its related concepts. See [START_REF] Fürer | Faster integer multiplication[END_REF] for details.

Primitive and principal roots of unity. Let R be a commutative ring with units. Let N > 1 be an integer. An element

ω ∈ R is a primitive N -th root of unity if for 1 < k ≤ N we have ω k = 1 ⇐⇒ k = N . The element ω ∈ R is a principal N -th root of unity if ω N = 1 and for all 1 ≤ k < N we have N -1 j=0 ω jk = 0. ( 8
)
In particular, if N is a power of 2 and ω N 2 = -1, then ω is a principal N -th root of unity. The two notions coincide in fields of characteristic 0. For integral domains every primitive root of unity is also a principal root of unity. For non-integral domains, a principal N -th root of unity is also a primitive N -th root of unity unless the characteristic of the ring R is a divisor of N .

The discrete Fourier transform (DFT). Let ω ∈ R be a principal N -th root of unity. The N -point DFT at ω is the linear function, mapping the vector ⃗ a = (a0, . . . , a N -1 ) T to ⃗ b = (b0, . . . , b N -1 ) T by ⃗ b = Ω⃗ a, where Ω = (ω jk ) 0≤j,k≤N -1 . If N is invertible in R, then the N -point DFT at ω has an inverse which is 1 N times the N -point DFT at ω -1 .

The fast Fourier transform. Let ω ∈ R be a principal N -th root of unity. Assume that N can be factorized to JK with J, K > 1. Recall Cooley-Tukey factorization formula [START_REF] Cooley | An algorithm for the machine calculation of complex fourier series[END_REF] DFT JK = (DFT

J ⊗ I K )D J,K (I J ⊗ DFT K )L JK J , (9)
where, for two matrices A, B over R with respective formats m × n and q × s, we denote by A ⊗ B an mq × ns matrix over R called the tensor product of A by B and defined by

A ⊗ B = [a k B] k, with A = [a k ] k, . (10)
In the above formula, DFT JK , DFT J and DFT K are respectively the N -point DFT at ω, the J-point DFT at ω K and the K-point DFT at ω J . The stride permutation matrix L JK J permutes an input vector x of length JK as follows

x[iJ + j] ↦ x[jJ + i], (11)
for all 0 ≤ j < J, 0 ≤ i < K. If x is viewed as a K × J matrix, then L JK J performs a transposition of this matrix.

The diagonal twiddle matrix D J,K is defined as

D J,K = J-1 ⊕ j=0 diag(1, ω j , . . . , ω j(K-1) ), (12)
Formula ( 9) implies various divide-and-conquer algorithms for computing DFTs efficiently, often referred to as fast Fourier transforms (FFTs). See the seminal papers [START_REF] Püschel | SPIRAL: Code generation for DSP transforms[END_REF] and [START_REF] Frigo | The design and implementation of FFTW3[END_REF] by the authors of the SPIRAL and FFTW projects, respectively. This formula also implies that, if K divides J, then all involved multiplications are by powers of ω K .

BLOCKED FFT ON THE GPU

In the sequel of this section, let ω ∈ R be a principal N -th root of unity. In the factorization of the matrix DFT JK , viewing the size K as a base case and assuming that J is a power of K, Formula (9) translates into a recursive algorithm.

This recursive formulation is, however not appropriate for generating code targeting many-core GPU-like architectures for which, formulating algorithms iteratively facilitates the division of the work into kernel calls and thread-blocks. To this end, we shall unroll Formula (9). Notation 3. Assuming c = 0, that is, N = K e , we define the following linear operators, for i = 0, . . . , e -1:

Ui(ω) = I K i ⊗ DFT K (ω K e-1 ) ⊗ I K e-i-1 ⋅ I K i ⊗ D K,K e-i-1 (ω K i ) , Vi(ω) = I K i ⊗ L K e-i K , Wi(ω) = I K i ⊗ L K e-i K e-i-1 ⋅ D K,K e-i-1 (ω K i ) . (13)
Remark 1. We recall two classical formulas for tensor products of matrices. If A and B are square matrices over R with respective orders a and b, then we have

A ⊗ B = L ab a ⋅ (B ⊗ A) L ab b . ( 14
)
If C and D are two other square matrices over R with respective orders a and b, then we have

(A ⊗ B) ⋅ (C ⊗ D) = (A ⋅ C) ⊗ (B ⋅ D). ( 15
)
Our GPU implementation reported in Section 6 is based on the following two results. We omit the proofs, which can easily be derived from Remark 1 and the Cooley-Tukey factorization formula; see [START_REF] Moreno Maza | Fast polynomial arithmetic on a GPU[END_REF]. Proposition 2. For i = 0, . . . , e -1, we have

Ui(ω) = Vi(ω) I K e-1 ⊗ DFT K (ω K e-1
) Wi(ω) [START_REF] Li | The modpn library: Bringing fast polynomial arithmetic into maple[END_REF] The following formula reduces the computation of a DFT on K e points to computing e DFT's on K points. Proposition 3. The following factorization of DFT K e (ω) holds:

DFT K e (ω) = U0(ω)⋯Ue-1(ω)Ve-1(ω)⋯V0(ω).

(17)

IMPLEMENTATION

We have realized a GPU implementation in the CUDA language of the algorithms presented in Sections 3 and 5. We have used the third and the fourth Generalized Fermat primes from Table 1, namely P3 ∶= (2 63 + 2 34 ) 8 + 1 and P4 ∶= (2 62 + 2 36 ) 16 + 1. We have tested our code and collected the experimental data on three different GPU cards.

In this section, we discuss implementation techniques. Our experimental results are reported in Section 7.

Parallelization. Performing arithmetic operations on vectors of elements of Z pZ has inherent data parallelism, which is ideal for implementation on GPUs. In our implementation, each arithmetic operation is computed by one thread. An alternative would be to use multiple threads for computing one operation. However, it would not improve the performance mostly due to overhead of handling propagation of carry (in case of addition and subtraction), or increased latency because of frequent accesses to global memory (in case of twiddle factor multiplications). Memory-bound kernels. Performance of our GPU kernels are limited by frequent accesses to memory. Therefore, we have considered solutions for minimizing memory latency, maximizing occupancy (i.e. number of active warps on each streaming multiprocessor) to hide latency, and maximizing IPC (instructions per clock cycle). Location of data. At execution time, each thread needs to perform computation on at least one element of Z pZ, meaning that it will read/write at least k digits of machine-word size. Often, in such a scenario, shared memory is utilized as an auxiliary memory, but this approach has two shortcomings. First, on a GPU, each streaming multiprocessor has a limited amount of shared memory which might not be large enough for allowing each thread to keep at least one element of Z pZ (which depending on the value of k, can be quite large). Second, using a huge amount of shared memory will reduce the occupancy. At the same time, there is no use for texture memory and constant memory for computing over Z pZ. Conclusively, the only remaining solution is to keep all data on global memory. Maximizing global memory efficiency. Assume that for a vector of N elements of Z pZ, consecutive digits of each element are stored in adjacent memory addresses. Therefore, such a vector can be considered as the row-major layout of a matrix with N rows and k columns. In practice, this data structure will hurt performance due to increased memory overhead that is caused by non-coalesced accesses to global memory. In this case, an effective solution is to apply a stride permutation L kN k on all input vectors (if data is stored in a row-major layout, this permutation is equivalent to transposing the input to a matrix of k rows and N columns). Therefore, all kernels are written with the assumption that consecutive digits of the same element are N steps away from each other in the memory. As a result, accesses to global memory will be coalesced, increasing memory load and store efficiency, and lowering the memory overhead. Decomposing computation into multiple kernels. Inside a kernel, consuming too many registers per thread can lower the occupancy, or even worse, can lead to register spilling. In order to prevent register spilling, register intensive kernels are broken into multiple smaller kernels. Size of thread blocks. Our GPU kernels do not depend on the size of a thread block. So, we choose a configuration for a thread block that will maximize the percentage of occupancy, the value of IPC (instruction per clock cycle), and bandwidth-related performance metrics such as the load and store throughput. We have achieved the best experimental results for thread blocks of 128 threads, or 256 threads. Effect of GPU instructions on performance. Our current implementation is optimized for the primes P3 ∶= (2 63 + 2 34 ) 8 + 1 and P4 ∶= (2 62 + 2 36 ) 16 + 1. Therefore, we rely on 64-bit instructions on GPUs. As it is explained in [6], even though 64-bit integer instructions are supported on NVIDIA GPUs, at compile time, all arithmetic and memory instructions will first be converted to a sequence of 32-bit equivalents. This might have a negative impact on the overall performance of our implementation. Specially, compared to addition and subtraction, 64-bit multiplication is computed through a longer sequence of 32-bit instructions. Finally, using 32-bit arithmetic provides more opportunities for optimization such as instruction level parallelism.

EXPERIMENTATION

We compare our implementation of FFT over a big prime field against a comparable approach based on FFTs over small prime fields. To be precise, we implement the two approaches discussed in Section 2. Recall that the first approach computes an FFT of size N over a big prime field of the form Z pZ where p is a SRGFN of size k machine words. The second approach uses s = 2 k half-machine word primes p1, . . . , ps and proceeds as follows:

1. projection:

compute the image fi of f in Z p1Z[x], . . . , Z p k Z[x], for i = 1, . . . , k, 2. images: compute the DFT of fi at ωi in Z piZ[x], for i = 1, .
. . , k (using the CUMODP library [START_REF] Moreno Maza | Fast polynomial arithmetic on a GPU[END_REF]), 3. combination: combine the results using CRT so as to obtain a DFT of f at ω. We use half-machine word primes (instead of machine-word primes as discussed in Section 2) because the small prime field FFTs of the CUMODP library impose this choice. Experimental results are gathered in Section 7.1.

We also have implemented and tested a sequential, CPU version of both approaches. For the small prime field approach, we use the NTL library [START_REF] Shoup | Number theory library (ntl) for c++[END_REF], supporting FFT modulo machine-word size primes of 60 bits. However, for the big prime field approach, we have implemented our own arithmetic in a sequential C++ program. Experimental results are gathered in Section 7.2.

Big prime vs small prime on the GPU

The output of the two approaches is the DFT of a vector of size N over a ring R which is either a prime field or a direct product of prime fields, and for which each element spans k machine-words. Hence these two approaches are equivalent building blocks in a modular method. For realizing the benchmark, first, we perform the reduction step, followed by computing s = 2k FFTs of size N over small prime fields. In the small field case, we use the highly optimized implementation of the following FFT algorithms from the CUMODP library (see [START_REF] Moreno Maza | Fast polynomial arithmetic on a GPU[END_REF][START_REF] Moreno Maza | Solving bivariate polynomial systems on a GPU[END_REF] and [START_REF] Haque | Dense arithmetic over finite fields with the CU-MODP library[END_REF]): the Cooley-Tukey FFT algorithm (CT), the Cooley-Tukey FFT algorithm with precomputed powers of the primitive root (CT-pow), and the Stockham FFT algorithm. The above codes compute DFTs for input vectors of 2 n elements, where 20 ≤ n ≤ 26 is typical.

Our CUDA implementation of the big prime field approach computes DFT over Z pZ, for P3 ∶= (2 63 + 2 34 ) 8 + 1 and P4 ∶= (2 62 + 2 36 ) 16 + 1, and input vectors of size N = K e where K = 16 for P3, and K = 32 for P4. Furthermore, for P3, we have 2 ≤ e ≤ 5, while for P4 (due to the limited size of global memory on a GPU card), we have 2 ≤ e ≤ 4.

The benchmark is computed on an NVIDIA Geforce GTX 760M (CC 3.0), an NVIDIA Tesla C2075 (CC 2.0), and an NVIDIA Tesla M2050 (CC 2.0). The first card has effective bandwidth of 48 GB/s, with 4 streaming multiprocessor, and the total number of 768 CUDA cores. Figures 1 and2 show the speedup of the big prime field FFT compared to the small prime field approach, measured on the first GPU card. Moreover, Table 2 presents the running times of computing the benchmark on the mentioned GPU cards. In each table, the first three columns present the running time of computing the small prime field FFT based on the Cooley-Tukey algorithm, the Cooley-Tukey FFT algorithm with precomputed powers of the primitive root, and the Stockham algorithm, respectively. Also, the last column presents the running time of computing the big prime field FFT.

As it is reported in [START_REF] Moreno Maza | Fast polynomial arithmetic on a GPU[END_REF], the FFT algorithms of the CU-MODP library gain speedup factors for vectors of the size 2 16 and larger, therefore, the input vector should be large enough to keep the GPU device busy, and thus, provide a high percentage of occupancy. This explains the results displayed on Figures 1 and2; for both primes P3 and P4, when N = K 2 and N = K 3 , our big prime field FFT approach significantly outperforms the small prime field FFT approach.

More importantly, for both primes P3 and P4, and with vectors of size N = K 4 , our experimental results demonstrate that computing the big prime field FFT is competitive with the small prime field approach in terms of running time. For both primes P3 and P4, we can compute FFT for an input vector of size N = K 4 , which is equivalent of 2 16 and 2 20 elements, respectively, and is large enough to cover most of the practical applications.

Eventually, for P3, and for a vector of size N = K 5 , the Cooley-Tukey (with precomputation) and Stockham FFT codes are slightly faster than the big prime field FFT. Nevertheless, for each of the tested big primes, there is a bit size range of input vectors over which the big prime field approach outperforms the small prime approach, which is coherent with the analysis of Section 2. For P3 ∶= ( 2 ) for computing the benchmark for N = K e for P3 and P4 (timings in milliseconds). Figure 3 shows the percentage of time spent in each operation in order to compute the big prime field FFT on a randomly generated input vector of size N = K 4 (measured for both primes and on the first mentioned GPU card). As it illustrated, for both primes, computation follows a similar pattern, with multiplication by twiddle factors as the main bottleneck. Finally, Table 5 presents the profiling data for computing the base-case DFT K on a GTX 760M GPU.

CPU vs GPU implementations

Table 3 presents the sequential running time of computing the benchmark for both primes on three different CPUs (measured in milliseconds). In addition, Table 4 shows the speedup range for computing the small and the big prime field approaches on CPU and GPU. For each prime, the first and the second column show the lowest and the highest running time of the same approach on CPU and GPU, respectively. Also, the last column contains the lowest and the highest speedup ratio of computing the same approach on CPU to its counterpart on GPU. .

CONCLUSION

Our results show the advantage of the big prime field approach. To be precise, for a size range of vectors, one can find a suitable large prime modulo which FFTs outperform the CRT-based approach. Our current implementation is preliminary and much optimization is possible, in particular for the multiplication in the big prime field. We anticipate better recent GPU cards with larger memory bandwidth and better support for 64-bit arithmetic. The CUDA code presented with this article is part of the CUMODP library freely available at www.cumodp.org.

Algorithm 3 2 : 1 ∑

 321 Computing xy ∈ Z pZ for x, y ∈ Z pZ procedure BigPrimeFieldMultiplication(fx, fy, r, k) 1: We compute the polynomial product fu = fxfy in Z[T ] modulo T k + 1. Writing fu = k-i=0

Figure 2 :

 2 Figure 2: Speedup diagram for computing the benchmark for a vector of size N = K e (K = 32) for P4 ∶= (2 62 +2 36 ) 16 +1.

Figure 3 :

 3 Figure 3: Running time of computing DFT N with N = K 4 on a GTX 760M GPU.

Table 1 :

 1 SRGFNs of practical interest.

	(2 63	+ 2 53	) 2	+ 1	2 106
	(2 64	-2 50	) 4	+ 1	2 200
	(2 63	+ 2 34	) 8	+ 1	2 272
	(2 62	+ 2 36	) 16	+ 1	2 576
	(2 62	+ 2 56	) 32	+ 1	2 1792
	(2 63	-2 40	) 64	+ 1	2 2560
	(2 64	-2 28	) 128	+ 1	2 3584

p max{2 e s.t. 2 e p -1}

  Figure 1: Speedup diagram for computing the benchmark for a vector of size N = K e (K = 16) for P3 ∶= (2 63 + 2 34 ) 8 + 1.

		150.00 148.32		K 2		2.00	1.93		K 4
					K 3					K 5
							1.50	1.47
	Speedup	100.00			94.50	Speedup	1.00		1.16	0.79	1.31	0.87
		50.00		48.83
							0.50
			8.83	5.22	6.89
		0.00					0.00
			CT	CT-pow	Stockham			CT	CT-pow	Stockham
		50.00	49.34		K 2		2.00	1.82		K 4
					K 3
		40.00					1.50
					33.06
	Speedup	20.00 30.00		25.15		Speedup	1.00		1.00	1.07
		10.00					0.50
			3.08	1.94	2.29
		0.00					0.00
			CT	CT-pow	Stockham			CT	CT-pow	Stockham

Table 2 :

 2 63 + 2 34 ) 8 + 1, this range is [2 12 , 2 16 ] while P4 ∶= (2 62 + 2 36 ) 16 + 1, this range is [2 15 , 2 20 ]. Our GPU implementation of the big prime field arithmetic is generic and thus can support larger SRGFNs, see Table1. Running time of computing the benchmark for N = K e on GPU (timings in milliseconds).

		Computing the benchmark for N = K e
		for P 3 ∶= (2 63 + 2 34 ) 8 + 1 (K = 16)
		Measured on a NVIDIA GTX-760M GPU
	e	CT	CT-pow	Stockham	Big FFT
	2	8.30	2.73	5.29	0.05
	3	10.96	6.49	8.55	1.24
	4	50.49	30.29	34.37	26.06
	5	820.82	444.07	490.72	558.22
		Measured on a NVIDIA Tesla C2075 GPU
	e	CT	CT-pow	Stockham	Big FFT
	2	9.44	2.93	5.16	0.03
	3	11.72	6.27	7.54	0.89
	4	31.85	15.57	19.07	17.71
	5	418.58	191.57	205.13	371.48
		Measured on a NVIDIA Tesla M2050 GPU
	e	CT	CT-pow	Stockham	Big FFT
	2	12.92	3.12	5.35	0.03
	3	15.35	6.66	8.00	0.88
	4	35.59	15.93	19.62	17.41
	5	424.98	198.46	206.71	364.88
		Computing the benchmark for N = K e
		for P 4 ∶= (2 62 + 2 36 ) 16 + 1 (K = 32)
		Measured on NVIDIA GTX-760M GPU
	e	CT	CT-pow	Stockham	Big FFT
	2	18.30	9.33	12.26	0.37
	3	62.98	39.74	46.72	20.40
	4	1772.9	974.01	1042.62	971.28
		Measured on NVIDIA Tesla C2075 GPU
	e	CT	CT-pow	Stockham	Big FFT
	2	19.82	9.56	11.56	0.27
	3	44.50	23.39	27.98	15.16
	4	891.35	437.29	464.69	695.02
		Measured on NVIDIA Tesla M2050 GPU
	e	CT	CT-pow	Stockham	Big FFT
	2	27.22	9.91	11.62	0.27
	3	51.81	23.93	28.60	14.80
	4	902.35	449.53	465.51	678.34

Table 3 :

 3 Running time of computing the benchmark for N = K e using NTL library on CPU (timings in milliseconds).

		Measured on		Measured on		Measured on
		Intel Xeon X5650 @ 2.67GHz CPU	AMD FX(tm)-8350 @ 2.40GHz CPU	Intel Core i7-4700HQ @ 2.40GHz CPU
	Computing the benchmark for N = K e for P3 ∶= (2 63 + 2 34 ) 8 + 1 (K = 16)	Computing the benchmark for N = K e for P3 ∶= (2 63 + 2 34 ) 8 + 1 (K = 16)	Computing the benchmark for N = K e for P3 ∶= (2 63 + 2 34 ) 8 + 1 (K = 16)
	e NTL Small FFT	NTL Big FFT	e NTL Small FFT	NTL Big FFT	e NTL Small FFT	NTL Big FFT
	2	2.51	1.85	2	4.06	4.13	2	3.12	0.73
	3	23.19	35.08	3	16.06	20.01	3	14.19	21.06
	4	372.19	750.40	4	296.00	528.00	4	232.76	505.96
	Computing the benchmark for N = K e	Computing the benchmark for N = K e	Computing the benchmark for N = K e
		for P4 ∶= (2 62 + 2 36 ) 16 + 1 (K = 32)		for P4 ∶= (2 62 + 2 36 ) 16 + 1 (K = 32)		for P4 ∶= (2 62 + 2 36 ) 16 + 1 (K = 32)
	e NTL Small FFT	NTL Big FFT	e NTL Small FFT	NTL Big FFT	e NTL Small FFT	NTL Big FFT
	2	14.94	12.91	2	12.00	8.00	2	12.48	9.79
	3	384.10	692.16	3	296.00	396.00	3	233.26	496.03
	4	11303.76	33351.29	4	10128.00	22992.00	4	7573.65	26089.53

Table 4 :

 4 Speedup ratio ( T CPU

	T GPU

  Computing the benchmark for N = K e for P3 ∶= (2 63 + 2 34 ) 8 + 1 (K = 16) (timings in milliseconds)Computing the benchmark for N = K e for P4 ∶= (2 62 + 2 36 ) 16 + 1 (K = 32) (timings in milliseconds)

	e SmallFFT CPU SmallFFT GPU	Speed-up	e	SmallFFT CPU	SmallFFT GPU	Speed-up
	2	2.51 -4.06	2.73 -12.92	0.19X -1.48X	2	12.00 -14.94	9.33 -27.22	0.44X -1.60X
	3	14.19 -23.19	6.27 -15.35	0.92X -3.69X	3	233.26 -384.10	23.39 -62.98	3.70X -16.42X
	4 232.76 -372.19	15.57 -50.49	4.61X -23.90X	4 7573.65 -11303.76 437.29 -1772.92 4.27X -25.84X
	e BigFFT CPU	BigFFT GPU	Speed-up	e	BigFFT CPU	BigFFT GPU	Speed-up
	2	0.73-4.13	0.03-0.05	14.6X -137.6X	2	8.00 -12.91	0.27 -0.37	21.62X -47.81X
	3	20.01-35.08	0.88-1.24	16.13X -39.86X	3	396.00 -692.16	14.80 -20.80	19.03X -46.76X
	4 505.96 -750.40	17.41-26.06	19.41X -43.10X	4 22992.00 -33351.29 695.02 -971.28 23.67X -479.62X

Table 5 :

 5 Profiling results for computing base-case DFT K on a GTX 760M GPU (collected using NVIDIA nvprof).

	Measured on a GTX760M GPU	P 3 = (2 63 + 2 34 ) 8 (K = 16)	P 4 = (2 62 + 2 36 ) 16 (K = 32)
	Metric Description	Mult by r	Add/Sub	Mult by r	Add/Sub
	Achieved Occupancy	74%	45%	62%	46%
	Device Memory Read Throughput	22.0 GB/s	20.7 GB/s	22.50 GB/s	10.00 GB/s
	Device Memory Write Throughput	22.0 GB/s	20.7 GB/s	23.19 GB/s	10.06 GB/s
	Global Memory Replay Overhead	0.28	0.03	0.25	0.01
	Executed IPC	0.72	2.41	0.78	4.56
	Instruction Replay Overhead	0.47	0.13	0.53	0.028
	Global Store Throughput	22.08 GB/s	20.74 GB/s	44.42 GB/s	9.91 GB/s
	Global Load Throughput	24.57 GB/s	20.91 GB/s	46.39 GB/s	10.22 GB/s
	Global Memory Load Efficiency	90.44%	43.60%	48.70%	98.86%
	Global Memory Store Efficiency	94.95%	43.75%	49.35%	99.99%

Acknowledgements

This work is supported by IBM Canada Ltd (CAS project 880) and NSERC of Canada (CRD grant CRDPJ500717-16).


H Ramezani Dana

P Casari

R El Abdi

S Fréour

F Jacquemin

CHARACTERISATION OF THE HYGRO-THERMO-MECHANICAL BEHAVIOUR OF ORGANIC MATRIX COMPOSITES INSTRUMENTED WITH OPTICAL FIBRES: A STUDY OF INTERFACIAL BONDING

Keywords: hygroscopic ageing, mechanical characterisation, interfacial bonding, optical fibre

Moisture diffusion can decrease the mechanical stiffness and strength of organic matrix composites. Recently, Fibre Bragg Grating (FBG) sensors have been used in order to study the hygroscopic ageing of polyester/glass composites at room temperature. In the present study, the hygroscopic strain measurement of polyester /glass composite has been achieved at room temperature. Furthermore, measurements have been done at higher temperatures in order to better understand the combination of both hygroscopic ageing and varying temperatures on the mechanical properties of these composite samples. The Bragg wavelength ( ) was found to shift linearly over a temperature range from room temperature to 35 °C. Beyond 35 °C, the Bragg wavelength does not linearly vary as a function of the temperature. A strong variation of the Bragg wavelength above a specific temperature threshold was found. This could be explained according to two mechanisms. Firstly, hygroscopic ageing could result in a decrease of the glass transition temperature (T g ) of the polymer matrix. Therefore viscoelastic behavior may appear beyond this temperature. Hygroscopic ageing could also degrade the interfacial shear strength between the fibre and the resin. Differential Scanning Calorimetry (DSC) analyses showed that hygroscopic aging does not affect the Tg of the polyester resin. Furthermore, it has been found that the hygroscopic aging degrades the interfacial adhesion of the optical fibre/polyester according to the considerable decreases of the interfacial shear strength observed in practice.

Introduction

The organic matrix composites (OMC's) have been used widely in structural applications over the past three to four decades [START_REF] Hollaway | A review of the present and future utilization of FRP composites in the civil infrastructure with reference to their important in-service properties[END_REF]. The service life of such structural components is strongly associated to their durability, which itself depends on the environmental conditions they are subjected to. Recently, the effects of humidity on the mechanical properties of polymer matrix composites have been reported by several authors [START_REF] Post | Investigation of 3D Moisture Diffusion Coefficients and Damage in a Pultruded E-glass/Polyester Structural Composite[END_REF][START_REF] Lin | Moisture sorption-desorption-resorption characteristics and its effect on the mechanical behavior of the epoxy system[END_REF][START_REF] Adams | The effect of immersion in sea water on the dynamic properties of fibrereinforced flexibilised epoxy composites[END_REF][START_REF] Adams | The dynamic properties of fibre-reinforced polymers exposed to hot, wet conditions[END_REF]. Glass fibers and polyester resin, the two constituents of the composite samples used in this study, present heterogeneous hygroscopic properties. For example, the coefficient of moisture expansion (CME) and the maximum moisture absorption capacity (Ms) are entirely different for these two materials. As a result of this, moisture absorption gives rise to the hygroscopic swelling of the matrix, which, coupled with the strong heterogeneity between the properties of both the matrix and the hydrophobic E-glass reinforcements, induces heterogeneous local stresses within the composite [START_REF] Jacquemin | Analytical modeling of transient hygro-elastic stress concentration-Application to embedded optical fiber in a non-uniform transient strain field[END_REF]. Considerable efforts have been made by researchers to develop analytical models that can predict the mechanical states occurring during both the transient stage and the permanent regime of the moisture diffusion process of organic matrix composites submitted to hygro-mechanical loads [START_REF] Neumann | Free-volume dependent moisture diffusion under stress in composite materials[END_REF][START_REF] Sar | Coupling moisture diffusion and internal mechanical states in polymers-A thermodynamical approach[END_REF]. Several studies highlight that the resulting mechanical states bring along structural damage [START_REF] Perreux | A study of the coupling between the phenomena of water, absorption and damage in glass/epoxy composite pipes[END_REF][START_REF] Weitsman | Effects of fluids on the deformation, strength and durability of polymeric composites-An overview[END_REF][START_REF] Weitsman | A correlation between fluid-induced damage and anomalous fluid sorption in polymeric composites[END_REF]. Therefore, the determination of internal mechanical states of polymer matrix composite submitted to the hygro-mechanical loads is critical for reliably predicting the long term behavior and durability of these composite structures. Different techniques have been developed in order to measure the hygro mechanical strain [START_REF] Karalekas | Monitoring of hygrothermal ageing effects in an epoxy resin using FBG sensor: A methodological study[END_REF][START_REF] Mihailov | Fiber Bragg grating sensors for harsh environments[END_REF]. Recently, optical fibre sensors have been used to achieve the measurement of strain and temperature [START_REF] Kang | Simultaneous Measurement of Strain and Temperature of Structures Using Fiber Optic Sensor[END_REF]. Optical fibre sensors provide many advantages compared to other strain measurement techniques, such as electrical and piezoelectric sensors, due to their ability to investigate the internal mechanical states within the bulk of the specimen [START_REF] Giaccari | On a direct determination of non-uniform internal strain fields using fibre Bragg gratings[END_REF][START_REF] Lai | Studies of hygrothermal degradation of a single fiber composite: An iterative approach with embedded optical sensors and numerical analysis[END_REF][START_REF] Lai | An experimental-numerical study of moisture absorption in an epoxy[END_REF]. The FBG is an optical filtering device that reflects light of a specific wavelength and is present within the core of an optical fiber waveguide. The reflected wavelength depends on the spacing of a periodic variation or modulation of the refractive index that is present. The sensing function of an FBG originates from the sensitivity of both the refractive index of the optical fiber and the grating period within the fiber to externally applied load and temperature. In the authors' previous study, the hygroscopic properties of glass/polyester composites have been identified by using FBG sensors [START_REF] Dana | Hygroscopic strain measurement by fibre Bragg gratings sensors in organic matrix compositesapplication to monitoring of a composite structure[END_REF]. Authors used FBG sensors in order to characterize the evolution of the hygroscopic strain field in unidirectional fibre reinforced composites, immersed in deionized water, during the moisture diffusion process. Materials properties such as the Coefficients of Moisture Expansion (CME) as well as the transverse diffusivity D 2 were identified from the treatment of the collected data [START_REF] Dana | Identification of moisture diffusion parameters in organic matrix composites[END_REF]. Additionally, the effect of temperature change on the Bragg wavelength ( ), and consequently on the hygroscopic strain, was presented. In this study, the shift in the fiber Bragg wavelength ( ) due to the influence of temperature has been examined for the glass/polyester composite specimen. It was observed that the response of Bragg wavelength shifting with temperature is linear over a range from 23° C to 35° C. On the contrary, the Bragg wavelength does not exhibit same behavior beyond 35° C. This could be explained by several phenomena. Among them, the two main factors are, firstly, the influence of the hygroscopic ageing of composite samples which could lead to a fall in T g , and secondly, a viscoelastic behavior appearing at higher temperatures. Moreover, the absorbed water can induce the debonding of the fibre-resin interface.

For that purpose, the characterisation of the glass transition temperatures of elaborated polyester resin and acrylate coating before and after hygroscopic aging has been achieved. In addition, the determination of the interfacial debonding stress of the silica optical fibre/polyester interface before and after hygroscopic aging has been carried out using the single-fibre pull-out test.

Materials and testing methods

Materials and specimens

An ortho-phthalic polyester resin (POLYLITE 420-731) has been used in this study. The gelling point of this resin is relatively short (about 20 min). Methyl Ethyl Ketone Peroxide (MEKP) has been used as a catalyst for initiating the polymerization of this polyester resin. The employed optical fibres are standard SMF-28e single mode fibres with a 250 µm diameter. This fibre is composed of a silica core and an acrylate coating. The acrylate coating is often used in standard optical fibres in order to protect the silica fibre core. Despite having a hydrophilic behavior, it ensures the protection of fibre against micro-bending and external physical stresses. The optical fibres used in this study were supplied by IX-fiber (SMF28-N/A-0.4 m). The polyester resin/optical fibres samples have been manufactured by a fairly simple method. In this method, a metallic mold, which consists of two steel plates, has been employed. The optical fibres have been placed between two plastic chocks -themselves covered by rubber gaskets (Figure 1A). The gaskets have the role of protecting the optical fibre against external stress, and consequently prevent fibre breakage. The metal mold was placed into an environmental chamber and heated up close to 50° C in order to reduce the temperature gradient between the mold and polyester resin. The prepared resins were then transferred into the mold and kept in an environmental chamber for 4 hours at ambient temperature (23° C). After polymerization, the optical fibre/resin specimens were withdrawn from environmental chamber and cut in order to carry out the pull-out test. After cutting and polishing, the specimen surfaces were cleaned with ethanol in order to remove any residual oil and dirt resulting from milling operation. The final dimensions of the specimens prepared for the pull-out test stated below were (25 × 10 × 2) mm 3 . A series of specimens, some containing optical fibres, and others without the fibres, were used for the ageing tests.

Figure1. Preparation of mold plate (A) and mold position in the heating chamber (B).

The tested unidirectional composite samples were made of E glass fibres embedded in an ortho-phthalic polyester (POLYLITE 420-731). The specimens were manufactured through the vacuum assisted resin infusion (VARI) method. Several composites, as well as neat resin square plates (1.3 mm thick) were fabricated. The fibre volume fraction (v f ) of the elaborated composites varied from 17 % to 22 %. The initial weights and dimensions of the specimens were recorded. The FBGs employed in this study were printed on standard SMF-28e single mode optical fibres with a 250 μm diameter. Bragg gratings (10 mm long) were uniform, and centered on a 1555 ± 0.2 nm wavelength. The use of FBGs in composite samples is limited by the ability to insert these sensors in these materials (so as to retrieve the optical signal) whilst minimizing their effect on the composites mechanical behavior. In order to address this issue, the optical fibre containing Bragg grating was inserted between two plates, bonded by polyester resin. The FBGs were either aligned with the reinforcing fibres or set perpendicularly to them. In the area of the grating, the acrylate coating was removed from the FBG in order to get the most direct strain provided by the surrounding material. The final size of the obtained instrumented samples is 90 × 90 × 3 mm 3 . In addition, this kind of architecture enables us to limit the stresses experienced by optical fibres during the fabrication process. Two series of specimens were manufactured. The first subset of samples were equipped with FBGs sensors. This group of samples was intended to provide the time-dependent evolution of the internal strain states throughout the moisture diffusion process. The second subset of specimens, which did not contain any Fibre Bragg Grating, were devoted to moisture uptake characterization by means of periodic mass measurements. These samples have the same final size as the instrumented specimens (90*90*3 mm 3 ). The experimental conditions that these samples were subjected to (with regards to the hygroscopic ageing process) is detailed in the next section.

Hygroscopic ageing process 2.2.1 Moisture absorption experiments

The hygroscopic ageing tests were carried out on both the composite and neat resin samples in order to identify their diffusive behavior. The initial weights of the samples were recorded. Thereafter, the specimens were immediately placed into deionized water maintained at controlled ambient temperature (20 °C). The change in mass was measured using an analytical electronic scale, marketed by SARTORIUS ® , with an accuracy of 0.1 mg. The weight gain versus the square root of time (√ ) curves for the composite and neat resin samples were then determined in order to follow their moisture absorption kinetics. The results (detailed in paragraph 3.1) show that the studied samples exhibit a Fickian diffusive behavior.

Identification of moisture diffusion parameters

The identification method used in this study for determining the diffusion parameters was explained in a previous work [START_REF] Dana | Identification of moisture diffusion parameters in organic matrix composites[END_REF]. This numerical method was intended for identifying both the macroscopic transverse moisture diffusion coefficient (D 2 ) and the maximum moisture absorption capacity (M s ) of the samples. The numerical method was based on the comparison between the 3D Fick's solution and the experimental measurement of the weight gain occurring during the diffusion process [START_REF] Pierron | A Novel procedure for identification of 3D moisture diffusion parameters on thick composites: theory, validation and experimental results[END_REF]. The method consists of finding out the unknown values of the problem by minimizing the standard deviation (Equation 1) using a Gauss-Newton algorithm:

∑ [ ( ) ] Equation 1
where ( ) is the moisture content calculated at time from the analytical solution to Fick's model, whereas is the corresponding experimental point.

Hygroscopic strain measurement method

FBG sensors have been used in order to achieve local strain measurements [START_REF] Colpo | Characterization of residual stresses in a single fibre composite with FBG sensor[END_REF][START_REF] Bosia | Deformation characteristics of composite laminates-part I: speckle interferometry and embedded Bragg grating sensor measurements[END_REF][START_REF] Guemes | Response of Bragg grating fiber-optic sensors when embedded in composite laminates[END_REF]. As mentioned previously, in the area close to the grating, the FBGs acrylate coating was removed in order to directly get the strain experienced by the surrounding composite. The acrylate coating could actually induce many side problems such as fibre matrix debonding. Moreover, it interferes with the local strain measurement due to the heterogeneous properties of acrylate compared to the composite. The diffusivity, maximum moisture absorption and coefficients of moisture expansion are different from one material to the other, requiring careful posttreatment of the experimental data gathered whenever the acrylate coating cannot be removed. However, the instrumented samples prepared in the present work enable us to directly determine the time-dependent evolution of the local internal axial strain states experienced by the optical fibre throughout the moisture diffusion process. Since the internal strains of the optical fibre are strongly related to the distribution of the mechanical states in the embedding material, this sensor provides the data necessary to characterize the local strain experienced by the composite surrounding the optical fibre. The principle of local strains measurement through FBGs has been presented in various articles [START_REF] Mihailov | Fiber Bragg grating sensors for harsh environments[END_REF][START_REF] Giaccari | On a direct determination of non-uniform internal strain fields using fibre Bragg gratings[END_REF][START_REF] Helan | Principle of Fiber Bragg Gratings Measurement[END_REF]. An FBG consists of a series of grating slices formed along the fibre axis. If the Bragg condition is fulfilled, a light signal propagating into the device can interfere constructively with the waves reflected by each of the slices. Consequently, a back reflected signal with a center wavelength commonly known as the Bragg wavelength ( ) is formed. The Bragg wavelength ( ) depends on the effective index of refraction ( ) and on the Bragg period ( ) of the grating according to the wellknown Bragg equation: Equation 2In the case where the FBG is submitted to a homogeneous axial strain and uniform temperature change , the Bragg wavelength experiences a deviation ( ) from the reference value ( ) corresponding to the unloaded state ( ; ): Equation 3Coefficients and depend on the nature of the optical fibre and the FBG printing parameters. Some preparatory tests have been carried out on the optical fibres (before the instrumentation of the composite specimens) at various controlled environmental conditions, thus permitting the characterisation of the FBGs sensitivity parameters, namely a =7.8*10 -6 /°C and b = 0.78*10 -6 /(μm/m), respectively. According to relation (3), under isothermal conditions, the Bragg wavelength shift is proportional to the axial strain ( ). According to Eq (3), one obvious major limitation of FBG sensors in composites is their simultaneous sensitivity to both deformation and temperature. As a result, in order to measure moisture absorption through strain variation, the temperature must be separately measured. This process allows thermo-optical effects, and thus thermomechanical deformations, to be deduced as a result the Bragg wavelength shifts caused by strains due to moisture absorption only. In the practical case considered in this work, the measured Bragg wavelength shift enables the strain experienced in the center of the instrumented samples to be determined from relation 3 relative to the value of interest; namely, the product . The Bragg wavelength shifting is recorded at ambient temperature (23°C) in order to investigate the effect of hygroscopic ageing. The measured Bragg wavelength is plotted as a function of water uptake for the composite specimens. As an example, the measured Bragg wavelength shift ( ) for the samples containing 17 % fibres equipped with FBGs positioned perpendicularly (90°) and/or parallel (0°) to the reinforcing fibres are 5.03 nm and 0.27 nm respectively. From relation [START_REF] Lin | Moisture sorption-desorption-resorption characteristics and its effect on the mechanical behavior of the epoxy system[END_REF], assuming that the thermal contribution (

) is negligible, hygroscopic strain values close to 4160 × 10 -6 and 260 × 10 - 6 have been found in the two principal directions.

Measurement of the Glass Transition Temperature T g

In order to determine the glass transition temperature (T g ) of elaborated samples before and after hygroscopic ageing, we applied the Differential Scanning Calorimetry (DSC) method. In this research work, a METTLER TOLEDO-STARe-DSC was used to measure the temperature and heat flows associated with thermal transitions in the elaborated resin samples. All prepared DSC samples had a mass of less than 20 mg. The samples were heated from room temperature (25° C) to 300° C at a rate of 20° C min -1 .

Characterization of the interfacial bonding stress 2.5.1 Interfacial shear strength measurement

Pull-out tests have been achieved in order to characterize the interfacial bonding stress of the optical fibre/ resin samples. It is necessary to recall that the pull out tests were done on the bare fibres. The pull out tests were performed at the IPR laboratory (Département Mécanique & Verres) at "Université de Rennes 1". The pull out tests were carried out on a LLOYD LR 50 K universal testing machine equipped with a force sensor of 100 N. Prior to the test, the base of the optical fibre/resin samples were fixed by lower jaw clamp. The optical fibre that comes out of the resin sample was wrapped around a pulley with a radius of 0.2 m. The upper pulley is mobile and subject to a constant imposed strain rate (here the selected strain rate is equal to 0.04 min -1 ). An automatic acquisition system enables the force and displacement data to be backed up in real time. At the end of the test, the matrix and fibres are collected in order to measure the embedded length (L e ) and study the microstructure by Scanning Electron Microscopy (SEM).

The principles of the pull out test

Initially, this test was developed on fibres surrounded by a polymer matrix [START_REF] Broutman | Measurement of the fiber-polymer matrix interfacial strength, Interfaces in composites[END_REF], before being used in Organic Matrix Composites. The principle of this test consists of submitting a fibre to mechanical loading. One end of the fibre is embedded in a matrix, and the other end is exposed to an increasing uniaxial tensile stress (Figure 2a) until decohesion at the interface occurs and, subsequently, the complete extraction of the fibre is achieved. Often, the curve of applied force versus displacement is plotted in order to analyze the experimental results (Figure 2b). The obtained curve can be divided into three areas. The first zone (A) corresponds to the elastic strain of the free part of the fibre and to interfacial debonding. Here, the maximum force ( ) corresponds to the entire fibre debonding from the matrix. Decohesion is followed by a sudden drop of the debonding force from to (B). Subsequently, the fibre extraction process starts from the force , involving friction between the fibre and the matrix [domain (C)] [START_REF] Delobelle | Embedment of optical fibre into a composite structure[END_REF]. The value of the interfacial shear strength can be determined from the following equation: Equation 4According to Equation 4, ( ) is expressed as a function of the embedded length ( ), the debonding force ( ) and the radius of the fibre ( ). The value of the debonding force ( ) can be determined from the applied load-displacement curve. In the present study, the embedded length ( ) of inserted fibre in the matrix has been measured through optical microscopy.

Results and discussion

Moisture uptake

Figure 3 shows the evolution of the moisture uptake as a function of the square root of time, obtained for the three glass fibre-polyester composite samples containing a volume fraction of the fibre equal to 17 %, 21 % and 22 % respectively. As one can see from figure 3(a), the three samples exhibit an almost linear evolution of the moisture uptake for several months until a pseudo-plateau -indicating that the saturation of the diffusion process has been reached. The curve of weight gain versus the square root of time obtained for the neat polyester resin sample is shown on Figure 3 (b). It can be noted that these samples present a Fickian diffusion behavior. According to figure 3, the maximum moisture absorption capacity of the neat polyester is nearly twice more that of the composite specimens (1.5 % versus 0.78 % respectively).

The diffusion kinetics of the polyester resin (with and without optical fibres) have been determined using the same gravimetric test. This was intended to assess the possible effect of the acrylate coating on the moisture absorption process in the studied samples. The computed evolutions obtained through the identified diffusive parameters found using Fick's model were plotted on the same graph as the experimental data collected on the resin samples (Figure 4). Figure 4 enables us to evaluate the reliability of the identification method. The moisture diffusion parameters of composites, neat polyester resin samples are summarised up in Table 1

Hygroscopic strain measurement

The measured Bragg's wavelength shift enables the strain experienced in the center of the instrumented samples to be determined from relation [START_REF] Lin | Moisture sorption-desorption-resorption characteristics and its effect on the mechanical behavior of the epoxy system[END_REF], assuming that the thermal contribution, namely, the product is negligible by comparison to the quantity of interest: . The Bragg wavelengths from before and after ageing and the corresponding strains for all tested samples were gathered in Table 2 It has been shown that the effect of temperature variations on the measured strain is significant. Therefore, the determination of purely hygroscopic strain requires a correction, particularly in the direction parallel to the reinforcing fibres. Various studies have been achieved in order to separate the purely hygroscopic strain from the temperature effect on the FBG measurement [START_REF] James | Simultaneous independent temperature and strain measurement using in-fibre Bragg grating sensors[END_REF][START_REF] Xu | Discrimination between strain and temperature effects using dual-wavelength fiber grating sensors[END_REF][START_REF] Rao | Simultaneous strain and temperature measurement of advanced 3-D Braided composite materials using an improved EFPI/FBG system[END_REF]. We propose a reliable method in order to correct the temperature effects before estimating the hygroscopic strains. This method is based on the analysis of the strain change experienced by the material due to a thermal loading once the aged specimens have reached a saturated state. On the basis of this technique, the effects induced by the temperature change from those of the hygroscopic strain on the measured Bragg wavelength shifting ( ) have been separated by plotting the Bragg wavelength as a function of temperature (Fig. 5). It should be remembered that the response of Bragg wavelength shifting to the temperature is linear over a range from 23° C to 35° C.

Figure 5. Bragg wavelength as a function of temperature for the composite specimens instrumented perpendicularly to the reinforcing fibres.

In the present study, the Bragg wavelength shifting is recorded at higher temperature in order to investigate the effect of hygroscopic ageing, simultaneously with temperature. The measured Bragg wavelength is plotted as a function of temperature on Figure 6 for the composite specimens instrumented by FBG positioned perpendicularly (Figure 6a) or parallel (Figure 6b) to the reinforcing fibres, respectively. We could find that beyond 35° C, the Bragg wavelength does not linearly vary as a function of the temperature anymore. The latter event could be due to multiple phenomena. This reaction could be explained in two ways. Firstly, by the influence of the hygroscopic ageing of composite samples which could have led to a fall of T g and secondly, by a viscoelastic behavior appearing a higher temperatures. Moreover, the absorbed water can induce the debonding of the fibre-resin interface [START_REF] Gu | Correlation of fiber pull-out strength and interfacial vibration damping techniques by micromechanical analysis[END_REF]. Thus, another possible explanation could come from the hygroscopic stress induced by absorbed moisture. In this case, when the mechanical states exceed the interfacial debonding stress of silica optical fibre-resin system, such nonlinearities as observed on figure 6 stress could be coupled with the fall of the glass-transition temperature. According to figure 6, this effect occurs in both types of composite samples where FBGs were located perpendicular and parallel to the reinforcing fibres. As a result, we propose the characterization of the glass transition temperature of elaborated polyester resin and acrylate coating before and after hygroscopic aging. On the other hand, the determination of interfacial debonding stress of silica optical fibre/polyester seems to be relevant for properly interpreting the obtained results. The characterization of that property will be handled by using the single-fibre pull-out test before and after hygroscopic ageing.

The glass-transition temperature measurement

Differential scanning calorimetry (DSC) analysis reveals information on the glass transition temperature (T g ) and degree of cure. For the thermosetting resins, the T g could be determined from the inflection point of the curve plotting the first derivative of the heat capacity as a function of temperature. Figure 7 (a) represents the DSC curves of the glass polyester composite specimens before hygroscopic ageing. The thermogram shows that the T g of composite specimens with different fibre contents is almost identical, with an average value of 50° C. An exothermic peak appears on the DSC thermograms of figure 7a. This means that the whole crosslinking, and therefore the entire hardening of the resin, was not fully achieved for these composite samples. It is important to note that in cross-linked systems, the glass transition temperature depends on the reticulation degree: it shifts to higher temperatures when the reticulation degree increases. The thermograms of polyester resin and acrylate coating in un-aged state are presented in figure 7(b). The thermograms highlight that the neat polyester resin has a glass transition temperature near to the Tg of glass polyester composite specimens (about 55° C). The thermograms indicate that the polyester resin has a lower glass transition temperature than the acrylate resin (55° C versus 160° C, respectively). We could also observe the same exothermic peak on the polyester resin thermogram, as on the other composites thermograms. We can also indicate that the acrylate glass transition temperature is considerably more important than in the polyester resin. Figure 8 represents the thermograms of polyester samples before and after hygroscopic ageing. Knowing that the determination of the T g could be difficult after ageing, we studied the effect of hygroscopic ageing on the three polyester specimens. The latter helps to improve the reliability of the obtained results. We found a very good reproducibility of the obtained results in non-aged and aged states for polyester resins. It can therefore be observed in figure 8 that hygroscopic ageing did not significantly affect the T g of polyesters (with an average value of 56°C after ageing). Some authors have reported a fall of T g during hygroscopic ageing which could vary from 10°C up to 80°C approximately [START_REF] Ghorbel | Hygrothermal effects on the physico-chemical properties of pure and glass fiber reinforced polyester and vinylester resins[END_REF][START_REF] Struik | Physical ageing of amorphous polymers and other Materials[END_REF][START_REF] Davies | Sea water aging of glass reinforced composites: shear behaviour and damage modelling[END_REF][START_REF] Castaing | Effects of water absorption and osmotic degradation on long term behavior of glass fiber reinforced polyester[END_REF][START_REF] Maurin | Seawater ageing of low styrene emission resins for marine composites: Mechanical behaviour and nano-indentation studies[END_REF]. For the unsaturated polyester resins, T g is expected to decrease considerably beyond 1 % of absorbed water. However, this event is not observed in our study. Nevertheless, peak height of the DSC thermogram decreases when the polyester experiences humid aging. This phenomenon is due to the increasing of the degree of polymerisation during the hygroscopic ageing, but it doesn't seem to significantly affect the T g in the present case. The transition temperature of polyester after aging, measured in this study, has been compared with the Tg of aged polyester predicted by the Fox law. We could note that in this case the two values (54.1°C, predicted by the Fox equation and 56.0°C, according to DSC measurements) are in reasonable agreement. It has been shown that the presence of fibres, and consequently the presence of interfaces, enhances the plasticisation of the composite matrix, in comparison with pure resin samples, and slows down the T g kinetics [START_REF] Ghorbel | Hygrothermal effects on the physico-chemical properties of pure and glass fiber reinforced polyester and vinylester resins[END_REF]. Unfortunately the obtained results on the aged acrylate sample are not exploitable. In fact, the thermogram corresponding to this sample shows just an intense peak representing an unusual reaction occurring in the capsule during the heating phase. A further study is planned, involving a new method of preparation for the aged acrylate coating samples. In view of these results, we observe that the glass transition temperature remains the same regardless of its aged state or/and non-aged state.

Evaluation of interfacial stress via the pull-out test

We studied interfacial debonding stress of silica optical fibre/polyester (O.F./polyester) before and after hygroscopic aging by using the single-fibre pull out test. Nevertheless, several samples broke either at the level of the free part of the optical fibre or at the fibre/sample junction (fibre output). We have not included these measures in the interfacial debonding analysis. Figure 9 shows the force-displacement curves of the O.F./polyester sample obtained via pull out-test before ageing.

The force-displacement curves of the O.F./polyester system exhibit a linear elastic portion up to F = 2 N, followed by a continuous force decrease. In this specific case, the decline of force might be controlled by the friction on the total embedding length of the optical fibre (Figure 9). Figure 10 shows the force-displacement curves of the aged O.F./polyester samples obtained via pull-out tests. Five samples have undergone the characterisation process. All obtained curves exhibit a linear elastic portion, evolving to a maximum, followed by a continuous or discontinuous decrease of the measured force (Figure 10 b). Contrary to the other specimens, sample number 2 deforms elastically up to F = 3 N, thus 3 times higher than the maximum value attained in the reproducibility tests (Figure 10 a). The force-displacement curves obtained on samples 3, 7 and 9 exhibit the classical 3 parts described previously in paragraph 2.5.2. As a consequence, for these specimens the force decrease appears to be controlled by the friction that takes place on the total embedded length of the optical fibre. Unfortunately, the force value of the sample number 2 is so great that it minimises the scale for the other specimens and does not enable the clear visualisation of the three steps of the forcedisplacement curve. For this purpose, the force-displacement curve of sample number 2 was shown in a separate graph (Figure 10. a). Force-displacement curves of sample number 5 don't highlight any friction between the fibre and resin during the extraction. Most of the time, the maximum force value undergone by the studied samples (except sample n°2) is weaker than 1 N. Compared to the un-aged samples (figure 9), which could reproducibly withstand a maximum force of around 2 N, the ultimate strength of the optical fibrepolyester resin interface is strongly affected by hygroscopic aging. Moreover, the effect induced by hygroscopic aging yields a strongly wide overall response of the material, according to the variability of the force-displacement curves drawn on figure 9. Such a strongly variable macroscopic behavior might result from lower scale effects, such as interface weakening or fibre-matrix partial debonding having occurred during the hygroscopic aging process. Since damage occurs at the microscopic scale, a lesser reproducibility of the results is expected. The values of F.O./polyester samples interfacial shear strength ( ) determined according to Eq. ( 4) from the experiment performed on un-aged samples are reported in Table 3. We find that the values of interfacial shear strength ( ) for these specimens are scattered but significantly weaker after immersion than the corresponding values obtained on un-aged samples. This is more obvious with the graph of figure 11: outcome could be explained by the physico-chemical modification of the interface in the presence water. Actually, plasticisation of the polymer matrix interface may significantly impair both its mechanical stiffness and ultimate strength.

According to Table 3, the hygroscopic aging has an important effect on the strength of the optical fibre / polyester interface. In the view of the results of paragraph 3.3, where we evaluated the influence of hygroscopic aging on the glass transition temperature of the polyester resin, we could explain this phenomenon by the relatively low T g of polyester. In fact, the polyester debonding stress in its aged state is strongly reduced, because the T g of the polyester is relatively low. These results are in agreement with the results described in paragraph 3.2. Actually, beyond 40°C the temperature is very close to the T g of polyester.

Consequently, a decrease of the mechanical properties of the fibre/resin interface does occur in the case of the polyester based system, after exposure to an ambient fluid. This was confirmed by the lower value of the maximum force required to pull out the optical fibre from the O.F./polyester system in its aged state compared to the un-aged samples.

Conclusions and perspectives

In this work, the hygroscopic strains experienced by composite samples immersed in deionized water were followed from the transient to the permanent stage of the diffusion process. The measurements were carried out using Fibre Bragg Grating (FBG) sensors. The thermal strains, due to the temperature change, occurring during the investigation have been determined in order to achieve a more realistic estimation of the hygroscopic strains. In this study we found that the hygroscopic aging does not affect the glass transition temperature of polyester resin, according to Differential Scanning Calorimetry (DSC) analysis. However, the hygroscopic aging influences the interfacial adhesion of the O.F./Polyester system. In fact the influence of hygroscopic aging on the interfacial adhesion of O.F./resin reveals that the moisture diffusion degrades this interface from a physico-chemical standpoint. We found that the interfacial shear strength of aged O.F./Polyester system considerably decreases because of the lower glass transition temperature of the polyester. We can conclude that the drop in the Bragg wavelength ( ) of glass/polyester composites in their saturated state should rather be associated with a decrease in their debonding strength ( ) due to moisture absorption.

Figure 2 .

 2 Figure 2.Tensile testing device (a) and typical curve of pull-out test (b)

FiberFigure 3 .

 3 Figure 3. Moisture absorption curves for (a) glass fibre -polyester composite samples and (b) neat polyester resin sample.

Figure 4 .

 4 Figure 4. Experimental and numerical (Fick model) curves of moisture uptake for neat polyester specimens with (a) and without (b) optical fibres.

Figure 6 .

 6 Figure 6. Bragg wavelength as a function of temperature for the polyester based composite specimens instrumented perpendicularly (a) and parallel (b) to the reinforcing fibres.

  could happen. It should be noted that the hygroscopic

Figure 7 .

 7 Figure 7. Glass transition temperature of composite samples (a) and resin samples (b) before hygroscopic aging.

Figure 8 .

 8 Figure 8. Hygroscopic ageing effect on the glass transition temperature.

Figure 9 .

 9 Figure 9. Force-displacement curves of the un-aged O.F./Polyester sample obtained owing to pull-out tests.

Figure 10 .

 10 Figure 10. Force-displacement curves of the aged O.F./Polyester samples obtained during pull-out tests performed under reproducibility conditions.

Figure. 11 .

 11 Figure. 11. Interfacial shear stress values of O.F./polyester samples before and after ageing.

Table 1 .

 1 : Moisture diffusion parameters of studied samples from the collected experimental gravimetric study.

		Transverse moisture	Maximum moisture
		diffusion parameter	absorption capacity
	Composite 17 %	D 2 ( mm 2 /s) 2.07*10 -7	M S (%) 0.85
	Composite 21 %	1.74*10 -7	0.78
	Composite 22 %	1.97*10 -7	0.79
	Neat polyester resin with optical fibre	2.31*10 -7	1.51
	Neat polyester resin	2.07*10 -7	1.31

Table 2 .

 2 . Hygroscopic strain for all instrumented specimens.

	v f (%) Angle (°)	(nm) Before	(nm) After		ε (10 -6 )	CME
			immersion	immersion	(nm)		(β)
	17	0	1554.68	1554.95	0.27	255	0.021
	17	90	1550.98	1556.01	5.03	4157	0.471
	21	0	1554.68	1555.14	0.46	379	0.042
	21	90	1551.48	1555.93	4.45	3677	0.397
	22	0	1554.68	1554.97	0.29	239	0.026
	22	90	1550.81	1555.18	4.37	3611	0.397

Table 3 .

 3 Interfacial shear stress values of O.F./polyester samples before and after ageing.

		Un-aged		Aged
	Sample #	Force (N)	(MPa)	Force (N)	(MPa)
	1	2.00	2.96	3.03	3.86
	2	3.99	2.15	0.04	0.10
	3	3.40	2.86	0.96	1.46
	4	5.91	3.66	0.34	0.55
	5	9.04	3.13	0.18	0.31
	Average		2.95		1.26
	Std		0.54		1.54
	deviation