Hanany Tolba

Nouha Dkhili

Julien Nou

Julien Eynard

Stéphane Thil

Stéphane Grieu

GHI forecasting using Gaussian process regression

Keywords: Global horizontal irradiance, forecasts, Gaussian processes, non-parametric regression, machine learning

In this paper, online Gaussian process regression (GPR) is used to model and forecast Global Horizontal Irradiance, at forecast horizons ranging from 30 min to 5 h. It is shown that the covariance function (or kernel) is a key element, deeply influencing forecast results. As a consequence, Gaussian processes with simple kernels and with more complex kernels have been tested and compared to the classic persistence model. Using two datasets of 45 days, it is shown that online GPR models based on quasiperiodic kernels outperform both the persistence model and GPR models based on simple kernels, including the widely used squared exponential kernel.

INTRODUCTION

Integration of fluctuating power generation (from renewable energy sources) in the electricity grid has always represented a major challenge. A grid operator has to ensure balance between electricity supply and demand. Unfortunately, many difficulties oppose the fulfillment of this goal, even more so in the context of smart grids due to the deployment of distributed power generators. The intermittency of solar power aggravates the problem of voltage regulation in distribution grids [START_REF] Moreno-Munoz | Very short term forecasting of solar radiation[END_REF]. Therefore, the stability of the network becomes dependent upon its successful compensation of supplydemand variations. As a result, maintaining the continuity and quality of service requires the ability to account for changes in the state of the network in real-time. To answer this question, statistical models, have been extensively used in recent years to forecast Global Horizontal Irradiance (GHI). These predictive models include the persistence model and auto-regressive models (AR) [START_REF] Boland | Time series modelling of solar radiation[END_REF]. Different models based on Artificial Neural Networks (ANN) can be found in the literature to forecast solar irradiance from meteorological and geographical data [START_REF] Mellit | A 24-h forecast of solar irradiance using artificial neural network: Application for performance prediction of a grid-connected PV plant at Trieste, Italy[END_REF][START_REF] Marquez | Forecasting of global and direct solar irradiance using stochastic learning methods, ground experiments and the NWS database[END_REF]. A comparative analysis of these methods can be found in [START_REF] Lauret | A benchmarking of machine learning techniques for solar radiation forecasting in an insular context[END_REF][START_REF] Reikard | Predicting solar radiation at high resolutions: A comparison of time series forecasts[END_REF], among others. However, the use of sophisticated methods remains scarce in the literature of solar irradiance. Especially, the relatively recent development in non-linear modeling, such as Gaussian Process Regression (GPR), remains very limited. In [START_REF] Lauret | A benchmarking of machine learning techniques for solar radiation forecasting in an insular context[END_REF] and [START_REF] Voyant | Machine learning methods for solar radiation forecasting: A review[END_REF], a GPR model is used to forecast GHI and is compared with other models. The conclusion is that using GPR models lead most of the time to the best results. However, in these studies, no procedure has been conducted to select an optimal kernel for the GPR model. In fact, only the Squared Exponential (SE) kernel, i.e. the default choice, has been considered. The present paper will show that this approach can be drastically improved by choosing an adapted kernel.

GAUSSIAN PROCESSES

Definition

A Gaussian Process (GP) is a collection of random variables, any finite number of which have a joint Gaussian distribution [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. A Gaussian process defines a prior over functions, which can be converted into a posterior over functions once some data has been observed. To indicate that a random function f (x) follows a Gaussian process, we write f (x) ∼ GP(µ(x), k(x, x )), where x and x are arbitrary input variables, µ(x) = E [f (x)] is the mean function (usually assumed to be zero) and

k(x, x ) = E (f (x) -µ(x))(f (x ) -µ(x ) T ) is the covari- ance function or kernel.

Covariance function (kernel)

A covariance function encodes our assumptions about the function which we wish to learn. This initial belief could be how smooth the function is or whether the function is periodic. The covariance function is also known as the kernel. Any function could be a covariance function as long as the resulting covariance matrix is positive semi-definite. For scalar-valued inputs, we briefly describe commonly used kernels (see [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] for an exhaustive list of kernels).

The Squared Exponential (SE) kernel is given by:

k SE (x, x ) = σ 2 exp -(x-x ) 2 /2 2
(1) where σ > 0 is the amplitude and > 0 is the correlation length parameter or characteristic length-scale. Intuitively, controls how fast the functions sampled from your GP oscillate.

The Rational Quadratic (RQ) kernel is given by:

k RQ (x, x ) = σ 2 1 + (x-x ) 2 /2α 2 -α
(2) with α > 0 and > 0. This kernel is equivalent to a scale mixture of SE kernels with different correlation length distributed according to a Beta distribution with parameters σ, α and [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF].

The Matérn class of kernels is defined by:

k Mν (x, x ) = σ 2 1 2 ν-1 Γ(ν) √ 2ν |x -x | ν • B ν √ 2ν |x -x | (3)
where σ is the amplitude, is the correlation length parameter, Γ is the standard Gamma function and B ν is the modified Bessel function of second kind of order ν. The parameter ν controls the degree of regularity (differentiability) of the resultant GP. If ν = 1 /2, the exponential kernel is obtained:

k E (x, x ) = σ 2 exp (-|x-x | / )
(4) This kernel corresponds to the continuous version of a classical discrete autoregressive model AR(1) (also known as Ornstein-Uhlenbeck process). Two other kernels of Matérn class are widely exploited in the literature. They correspond to the cases ν = 3 /2 and ν = 5 /2, respectively, and are defined by:

k M 3/2 (x, x ) = σ 2 1 + √ 3 |x -x | • exp - √ 3 |x -x | (5) k M 5/2 (x, x ) = σ 2 1+ √ 5 |x -x | + √ 5 (x -x ) 2 3 2 • exp - √ 5 |x -x | (6)
The periodic kernel is given by:

k Per (x, x ) = σ 2 exp - 2 sin 2 ( π(x-x ) /P ) 2 (7)
A periodic kernel assumes a globally periodic structure of period P in the function we wish to learn. The parameters σ and have the same interpretation found in the kernels previously defined.

Kernel composition

It is possible to combine several kernel functions to obtain a more complex kernel function: the only constraint is that the resulting covariance matrix must be a positive semi-definite matrix. Addition and multiplication are two ways of combining covariance functions while keeping the positive semi-definite property:

k(x, x ) = k 1 (x, x ) + k 2 (x, x ) (8) k(x, x ) = k 1 (x, x ) × k 2 (x, x ) (9)
For example, we can model a quasiperiodic GP by multiplying a periodic kernel by a non periodic kernel: this gives a way to transform a global periodic structure into a local periodic structure. In this paper, a wide variety of kernel structures are built to model GHI by adding or multiplying two kernels from the kernel families discussed in subsection 2.2.

GAUSSIAN PROCESS REGRESSION (GPR)

Standard Gaussian process regression

Consider the standard regression model:

y = f (x) + ε
(10) where x ∈ R D×1 is the input vector, f is the regression function, y is the observed value and ε ∼ N (0, σ 2 ε ) is an independent, identically distributed Gaussian noise. Gaussian process regression is a Bayesian nonparametric regression which assumes a GP prior over the regression functions [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]

: it consists in approximating f (x) ∼ GP(µ(x), k(x, x )) using a training set of n observations D = {(x i , y i ), 1 i n}.
As shorthand notation, we merge all the input vectors x i into a matrix X ∈ R n×D and all corresponding outputs y i into a vector y ∈ R n×1 , so that the training set can be written as (X, y).

From equation ( 10), it can be seen that:

y ∼ N µ(X), K + σ 2 ε I (11) where K = k(X, X) ∈ R n×n .
In this setting, the joint distribution of the observed data y and the latent noisefree function on the test points f * = f (X * ) is given by:

y f * ∼ N µ(X) µ(X * ) , K + σ 2 ε I K * K T * K * * (12)
where

K * = k(X, X * ) ∈ R n×n * and K * * = k(X * , X * ) ∈ R n * ×n * .
It can be shown that the posterior predictive density is then also Gaussian (see [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]):

f * | X, y, X * ∼ N (µ * , σ 2 * ) (13
) where:

µ * = µ(X * ) + K T * K + σ 2 ε I -1 (y -µ(X)) (14)
σ 2 * = K * * -K T * K + σ 2 ε I -1 K * ( 15
k * = [k(x * , x 1 ) . . . k(x * , x n )] T
(16) Equation ( 14) becomes:

µ * = µ(x * ) + k T * K + σ 2 ε I -1 (y -µ(x * )) (17
) Thus µ * , the mean prediction for f (x * ), can be written as a linear combination of kernel functions, each one centered on a training point:

µ * = µ(x * ) + k T * α = µ(x * ) + n i=1 α i k(x i , x * ) (18) where α = (K + σ 2 ε I) -1 (y -µ(x * )).
In the sequel, these coefficients α i will be referred to as parameters; they are updated each time a new observation is made (as opposed to the parameters of the kernel, referred to as hyperparameters, which are not updated once training is over (see subsection 3.3).

Online Gaussian Process Regression (OGPR)

An obstacle that could be met when using the standard GPR approach is the difficulty with which to incorporate a new training point or set of points. When a forecasting algorithm is run in situ, we do not have a fixed data set D = {(x i , y i ), 1 i n}. Rather, data is received one or a few observations at a time, and the total amount of information keeps growing. Moreover, updates are often required every hour, minute or even second. Using the above-mentioned methods in this case would rapidly incur prohibitive computational overheads, since the matrix (K + σ 2 ε I) has to be inverted, with a complexity O(n 3 ). A solution is to use the so-called Online Gaussian Process Regression (OGPR) approaches [START_REF] Csató | Gaussian Processes -Iterative Sparse Approximations[END_REF][START_REF] Huber | Recursive gaussian process: On-line regression and learning[END_REF][START_REF] Ranganathan | Online sparse gaussian process regression and its applications[END_REF][START_REF] Kou | Sparse online warped Gaussian process for wind power probabilistic forecasting[END_REF][START_REF] Bijl | Online sparse gaussian process regression using FITC and PITC approximations[END_REF]. Formally, let us suppose that we know the distribution of the GP given the first n training points (X, y) and for simplicity, let us consider µ = 0 (without loss of generality). Suppose that we make n + new observations regrouped in the matrix X + . Then:

K K + + σ 2 ε I n+n+ -1 = k(X, X) + σ 2 ε I n k(X, X + ) k(X + , X) k(X + , X + ) + σ 2 ε I n+ -1 (19) =      K + σ 2 ε I n A K + B K T + B T K ++ + σ 2 ε I n+ D      -1 (20)
Using the block matrix inversion theorem, we get:

K K + + σ 2 ε I n+n+ -1 = A -1 + A -1 B∆ -1 B T A -1 -A -1 B∆ -1 ∆ -1 B T A -1 ∆ -1 (21)
where

∆ = D -B T A -1 B ∈ R n+×n+ .
It can be seen that the inversion of the (n + n + ) × (n + n + ) matrix now only requires the inverse of A, which is already known, and the inversion of a n + × n + matrix: the computational cost is thus O(n 3 + n 2 ) rather than O((n + n + ) 3 ) when performing direct matrix inversion.

Training a GPR model

As mentioned in subsection 2.2, there is freedom in choosing the kernel function of a GP. Common choices, such as the kernel previously covered, include hyperparameters. The hyperparameters of the GPR model (10), denoted θ, group those of its kernel and the variance of the noise. For instance, the squared exponential k SE kernel defined by 1 contains two hyperparameters: the amplitude σ and the correlation length , to which is added the noise variance σ 2 ε : thus, in this case θ = [σ σ 2 ε ] T . These hyperparameters have to be estimated from the data. To do so, the probability of the data given the aforesaid hyperparameters is computed. Assuming a GPR model with Gaussian noise (10), the log marginal likelihood is given by [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]:

L(θ) = log [P(y|X, θ)] (22) = - 1 2 y T K + σ 2 ε I -1 y + log det K + σ 2 ε I + n log(2π) (23)
Estimation of the hyperparameters θ is usually achieved by maximizing the log marginal likelihood ( 22). The maximization process may be accelerated if the gradient of the log likelihood is known and a gradient-based algorithm, such as a conjugate gradient method, can be used [START_REF] Blum | Optimization of gaussian process hyperparameters using Rprop[END_REF]. The interested reader can turn to [START_REF] Moore | Fast methods for training gaussian processes on large datasets[END_REF], [START_REF] Ambikasaran | Fast direct methods for gaussian processes[END_REF] or [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] for more information.

DATA DESCRIPTION

The dataset is derived from measurements taken in Perpignan (southern France), at the PROMES-CNRS laboratory, using a Rotating Shadowband Irradiometer (RSI). Typical uncertainties of the RSI are about ±5%. The station is located approximately 20 km west of the Mediterranean Sea. Winter is mild and summer is hot and dry. In addition, there is a lot of wind, often resulting in a cloudless sky.

Evaluation of the GPR models included in the comparative study has been made using a dataset covering a period of 45 days, from June 5 th , 2015 to July 18 th , 2015, containing GHI data with a 30 min time step. If not specified otherwise, the first 30 days of the dataset are used for training (see subsection 5.1), whereas the last 15 days are used to evaluate models' forecasting skill.

MODELING AND FORECASTING RESULTS

Initialization and estimation of hyperparameters

For modelling GHI using a GPR model, we first need to identify the kernel function most adapted to this type of time series. The procedure of kernel identification adopted in this study is inspired by the work developed in [START_REF] Duvenaud | Structure discovery in nonparametric regression through compositional kernel search[END_REF]. The idea is to construct a wide variety of kernel structures by adding or multiplying kernels.

In particular, the kernel families discussed in subsection 2.2, as well as their sum and product are considered. While all the possible combinations have been evaluated, only the following kernels have been included in this study:

• simple kernels: k E , k SE , k RQ , k M 3/2 , k M 5/2 , k Periodic ; • quasiperiodic kernels; products: k Per×E , k Per×SE , k Per×RQ , k Per×M 3/2 , k Per×M 5/2 ; sums: k Per+E , k Per+SE , k Per+RQ , k Per+M 3/2 , k Per+M 5/2 .
Indeed, the other combinations of non-periodic kernels exhibited the same behavior as simple kernels. For clarity's sake, their results have not been displayed.

The hyperparameters of each kernel have been estimated from GHI training data via the minimization of the log marginal likelihood (22). Due to the non-convexity of L(θ), the optimization may not converge to the global maximum. The classical approach to tackle this issue is to use multiple starting points randomly selected from a specific prior distribution. In [START_REF] Chen | How priors of initial hyperparameters affect gaussian process regression models[END_REF], the authors consider different types of priors for the initial values of hyperparameters and investigate the influence of the priors on the predictability of GPR models. The results reveal that the estimates for k SE are robust regardless of the prior distributions, whilst they are very different using different priors for k Per which implies that the prior distributions have a considerable impact on the estimates of the parameters in the case of a periodic kernel.

Regarding the initial values of the hyperparameters θ, the following choices have been made.

• The correlation length has been chosen to be equal to the standard deviation of the training data. • When the tested kernels involve a periodic kernel, an initial value of one day has been chosen for the period P , in order to overcome problems that could be encountered in its estimation.

• The initial values of remaining hyperparameters, if any, are randomly drawn from a uniform distribution: θ i ∼ Uniform(0, 1).

Forecasting results using online GPR

The evaluation metric considered in this study is the normalized Root Mean Square Error (nRMSE):

nRMSE = 1 n * n * i=1 y test (i) -y forecast (i) 2 1 n * n * i=1 y test (i) (24)
where y test ∈ R n * ×1 is the test data and y forecast ∈ R n * ×1 is the forecast given by the considered models.

Let us assess the models' performance globally by comparing their forecasts at different time horizons. In Figure 1, nRMSE is plotted versus forecast horizon, for all GPR models listed in subsection 5.1, as well as the persistence model, used for reference. Detailed numerical results can be found in Table 1. First, it be seen considered models can be divided into three classes: persistence, which gives the worst results; GPR models based on simple kernels, with improved performance; and GPR models based on quasiperiodic kernels, which perform considerably better, especially at higher forecast horizons.

P e r s is . k E k R Q k M 3 / 2 k M 5 / 2 k S E k P e r × S E k P e r × M 5 / 2 k P e r × M 3 / 2 k P e r × E k P e r k P e r + S E k P e r + M 5 / 2 k P e r + M 3 / 2 k P e r + E k P e r + R Q k P e r × R Q 0.5 1 1.5 2
Even at the lowest forecast horizon (30 min), GPR models based on simple kernels give results comparable to the Regarding GPR models based on simple kernels, there does not seem to be a clear best-performing kernel. Depending on the forecast horizon, k M 3/2 , k M 3/2 , k RQ , and k SE all alternatively give the best results, while k E lags behind. It should also be noted that k SE -the kernel usually chosen in previous work on GHI forecasting -does not always provide the best results among simple kernels.

A realization of a GPR model based on a periodic kernel is a periodic signal. As a consequence, this model gives rather constant results in terms of nRMSE as the horizon increases. This GPR model gives good forecasting results for clear days, but it cannot explain the variability in GHI data due to atmospheric disturbances. GPR models based on quasiperiodic kernels, however, are being composed of both a periodic and a non-periodic kernel. As a consequence, they have the advantage of the periodic kernel, while still being able to explain GHI rapid changes during the day. Among quasiperiodic kernels, k Per×RQ outperforms other kernels, but followed closely by k Per+M 3/2 , k Per+M 5/2 and k Per+E (see Table 1).

To get further insight into models' performance, let us study the temporal evolution of GHI during both training and testing. Figures 2, 3 and 4 show 30 min, 4 h and 48 h forecasts obtained using the persistence model, displayed for reference, and three selected GPR models: the classic k SE kernel and two of the best-performing quasiperiodic kernels, k Per+SE and k Per×RQ . Here, a dataset of nine days is used (seven days for training and two days for testing).

Recall that during training each data sample is used, while during testing a new observation is incorporated each whole multiple of the forecast horizon. So, a new observation is taken into account every 30 min in Figure 2, every 4 h in Figure 3 and every 48 h in Figure 4 (which results in no update at all, since the testing dataset lasts two days). It also implies that training is identical in these three figures; only testing differs.

When looking at the training stage, it can be seen that each GPR model fits the data quite well. There are very few differences between k SE and k Per+SE , because they both generate rather smooth signals; the k Per×RQ kernel, however, allows for more irregular signals and thus follows GHI more closely during atmospheric disturbances. Now turn to the testing stage. When there is a new observation each 30 min (Figure 2), the persistence model is excellent and all GPR models give good results, with a slight advantage to quasiperiodic kernels, especially k Per×RQ . But when the horizon increases the difference between the models become apparent. When there is a new observation each 4 h, the simple kernel k SE struggles to explain GHI variation, as demonstrated by the large confidence interval obtained (Figure 3a). When no update is available it simply converges to its mean value (around 280 W m -2 ). This phenomenon is accentuated in Figure 4a: no updates are made during the entire testing and since this GPR model does not have additional information on the GHI behavior, it simply gives a constant value equal to its mean. On the contrary, quasiperiodic kernels have an additional information on GHI: it has a daily pattern. Therefore, when no update is available they converge to a periodic value. This is the advantage of quasiperiodic kernels: when no additional information is available, they reproduce the periodic behavior learned during training (see Figures 4b and4d). When comparing the quasiperiodic kernels, k Per×RQ is preferred because it allows for more drastic changes in the modeled signal, as can be seen in Figures 3b and3d.

CONCLUSION

The main objective of this paper was the elaboration of a Gaussian process regression model adapted to multihorizon GHI forecasting. In previous research on GHI forecasting using GPR models, the squared exponential kernel k SE was used as the default choice (see e.g. [START_REF] Lauret | A benchmarking of machine learning techniques for solar radiation forecasting in an insular context[END_REF][START_REF] Voyant | Machine learning methods for solar radiation forecasting: A review[END_REF]). However, as it has been explained throughout the present document, the kernel is a key element of GPR models. A comparison between several GPR models based on simple kernels and more complex kernels has shown that, while GPR models globally outperform the persistence model, GPR models based on quasiperiodic kernels are better suited to model and forecast GHI. The proposed interpretation is that the structure of GHI is better modeled through an omnipresent periodic component representing its global structure and a random latent component explaining rapid variations due to atmospheric disturbances.

  ) Let us examine what happens in the case of a single test point x * . Let k * be the vector of covariances between the test point and the n training points:

Fig. 1 .

 1 Fig. 1. nRMSE versus forecast horizon for all models.

  GPR model based on k Per×RQ kernel.

Fig. 2 .

 2 Fig. 2. 30 min forecasts, for the persistence model and OGPR models based on three different kernels: k SE , k Per×RQ and k Per+SE (dataset of nine days, with a training dataset of seven days and a testing dataset of two days).

Fig. 3

 3 Fig. 3. 4 h forecasts, for the persistence model and OGPR models based on three different kernels: k SE , k Per×RQ and k Per+SE (dataset of nine days, with a training dataset of seven days and a testing dataset of two days).

Fig. 4 .

 4 Fig. 4. 48 h forecasts, for the persistence model and OGPR models based on three different kernels: k SE , k Per×RQ and k Per+SE (dataset of nine days, with a training dataset of seven days and a testing dataset of two days).

Table 1 .

 1 Values of nRMSE obtained for all models. The best results are highlighted.

	F. horizon 30 min	1 h	2 h	3 h	4 h	5 h
	Persistence 0.2177 0.3485 0.5833 0.7935 0.9749 1.1130
	k E k M 3/2 k M 5/2 k RQ k SE	0.2183 0.2818 0.4028 0.4737 0.5998 0.6382 0.2125 0.2652 0.3574 0.3815 0.5075 0.5263 0.2126 0.2651 0.3558 0.3792 0.5042 0.5235 0.2179 0.2743 0.3739 0.4144 0.5377 0.5532 0.2153 0.2666 0.3387 0.3707 0.4660 0.5427
	k Per	0.2340 0.2441 0.2596 0.2720 0.2772 0.2779
	k Per+E k Per+M 3/2 k Per+M 5/2 k Per+RQ k Per+SE	0.1689 0.1931 0.2209 0.2292 0.2538 0.2518 0.1703 0.1954 0.2223 0.2308 0.2508 0.2548 0.1705 0.1956 0.2226 0.2312 0.2510 0.2553 0.1692 0.1929 0.2200 0.2355 0.2526 0.2529 0.1714 0.1986 0.2268 0.2345 0.2548 0.2588
	k Per×E k Per×M 3/2 k Per×M 5/2 k Per×RQ k Per×SE	0.1779 0.2098 0.2585 0.2480 0.3118 0.3192 0.1922 0.2280 0.2660 0.2920 0.3211 0.3188 0.1927 0.2289 0.2677 0.2947 0.3237 0.3220 0.1675 0.1899 0.2176 0.2197 0.2484 0.2424 0.2045 0.2506 0.3061 0.3462 0.3825 0.3886
	persistence model (nRMSE kernels already give better results (nRMSE 0.22), while quasiperiodic 0.18). As the horizon increases, performance of the persistence model degrades quickly, while those of GPR models degrade more slowly. At the highest horizon (5 h), persistence gives nRMSE 1.11; for simple kernels nRMSE 0.55; for quasiperiodic kernels nRMSE 0.30.


Stéphane Andrieux

Amel Ben Abda

Huy Duong Bui

Reciprocity principle and crack identification

In this paper we are concerned with the planar crack identification problem defined by a unique complete elastostatic overdetermined boundary datum. Based on the reciprocity gap principle, we give a direct process for locating the host plane and we establish a new constuctive identifiability result for 3D planar cracks.

Introduction

Let be a smooth open domain of R 2 or R 3 , and a curve of R 2 or a surface of R 3 smooth enough but not necessarily simply connected (cracks). We consider the following elastostatic problem (P):

(P) ⎧ ⎨ ⎩ div σ = 0 in \ σ • n = F m on ∂ σ • N = 0 on (1)
with:

• the linear elastic stress-strain law

σ = Aε (u) , (2)
• the relation between displacement and linearized strain

ε(u) = (∇u) s , ( 3
)
where A is the stiffness tensor, ∇ is the gradient operator and (.) s denotes the symmetric part of the tensor (.). Equations ( 4) are added to ensure the existence of a solution. Since the problem involves only the 'Neumann' boundary conditions, the solution is known up to a rigid displacement. The conditions [START_REF] Brésis | Analyse fonctionnelle[END_REF] are added to select a solution:

∂ F m = 0 and ∂ OM ∧ F m = 0, (4)
∂ u = 0 and

∂ OM ∧ u = 0, ( 5
)
where O is an arbitrarily fixed point of the space. We denote by u the unique solution of the considered problem. The aim of this paper is to give results on the identifiability of cracks lying inside the domain , when the overspecified boundary data of the problem (P) are available on the whole exterior boundary ∂ of the domain (prescribed surface force F m and displacement response on the boundary, denoted by u m which is supposed to be known). To our knowledge, theoretical results on the identifiability of cracks are sparse and are concerned mainly with the thermal conductivity equation in two dimensions [START_REF] Friedman | Determining cracks by boundary measurements Indiana[END_REF][START_REF] Alessandrini | Unique determination of multiple cracks by two measurements Quad[END_REF][START_REF] Bryan | A uniqueness result concerning the identification of a collection of cracks from finitely many elastostatic boundary measurements[END_REF].

For the scalar conductivity equation and in the particular case of plane cracks (or straight cracks in 2D), with given overspecified data on the whole boundary, Andrieux and Ben Abda [START_REF] Andrieux | Identification de fissures planes par une donnée de bord unique: un procédé direct de localisation et d'identification[END_REF] have proved a theorem on crack identifiability based on the concept of reciprocity gap. Andrieux [START_REF] Andrieux | Fonctionnelles d'écart à la réciprocité généralisée et identification de fissures par des mesures surabondantes de surface C[END_REF] has extended the result to other types of operators (non-symmetrical). The proof of the theorem given in [START_REF] Andrieux | Identification de fissures planes par une donnée de bord unique: un procédé direct de localisation et d'identification[END_REF] is direct: after showing the identity between the support of the field discontinuity and the domain occupied by the crack, one identifies the temperature discontinuity via its Fourier series. This result has been numerically checked by Bannour et al [START_REF] Bannour | A semi-explicit algorithm for the reconstruction of 3D planar cracks Inverse[END_REF] in two and three spatial dimensionals. Here we propose to study the elastostatic case. Consider a three-dimensional domain , containing one or many cracks lying on the same plane of R 3 , ⊂ . The space is referred to the direct orthonormal frame of reference (O, e 1 , e 2 , e 3 ), with x 1 , x 2 , x 3 denoting the associated coordinates. The plane is defined by the equation: N • x + c = 0, where N is the unit normal . In section 2, we shall give explicit formulae for localizing the plane . Section 3 is devoted to the complete identification of the crack .

Localization of the crack plane

Let V be the set of regular fields, such that the associated stress fields, via the elastic law, are in equilibrium inside :

V = {v ∈ [H 1 ( )] 3 , div Aε(v) = 0 in }. (6)
The reciprocity functional gap for considered problem is defined on the space V by:

RG(v) = ∂ {F m • v -u m • Aε(v) • n}.
Remark that whenever the domain is uncracked ( = ∅) the functional RG is identically zero on V , which is nothing but the reciprocity theorem of Maxwell-Betti in elastostatics. This remark justifies the name 'reciprocity gap functional'. The following lemma expresses the value of RG, defined with the boundary data by an integral over the crack surfaces.

Lemma 1. For any v of V , RG(v) = [[u]] • Aε(v) • N, ( 7
)
where

[[u]] is the jump of u across .
Proof. The result is obtained by applying Green's formula.

Remarks. Denoting by [[u]

] the extension of u to the whole plane , by letting it equal zero outside , we obtain immediately:

RG(v) = ∩ [[u]] • Aε(v) • • N.
Let us consider the particular case of isotropic elasticity, with E and ν being the Young modulus and the Poisson ratio respectively. We introduce the following family of displacement field (v ij ), with the parameters (or indexes) (ij ), the (k) component of which is defined by:

v ij k = A -1 kl mn E ij mn x l , (8) with
• (E ij ) canonical basis of second-order symmetrical tensors of R 3 :

E ij mn = 1 2 (δ i m δ j n + δ i n δ j m ). • A -1 compliance

tensor of the medium (inverse of A).

The indices take values in the set {1, 2, 3} and δ is the Kronecker delta. We use the convention of summation for repeated indexes. The stress fields σ ij associated to v ij are constant ones (zero divergence and hence belonging to the space V ):

σ ij = Aε(v ij ) = E ij . ( 9
)
Using lemma 1, it is easy to show that:

RG(v ij ) = N ⊗ [[u]] s ij ≡ R ij . ( 10
)
Denoting the tensor R with components R ij by:

R = N ⊗ [[u]] s ( 11
)
one can establish the following two propositions.

Proposition 1. The norm of the mean value of the displacement discontinuity, across the crack, is given explicitly by:

[[u]] = 2R 2 -(Tr R) 2 , ( 12
)
moreover, denoting the tangential discontinuity by:

[[u t ]] = [[u]] -[[u]] • N N, one gets ⎧ ⎨ ⎩ [[u t ]] = 2R 2 -(Tr R) 2 [[u]] • N = Tr R. ( 13
)
We assume now that the applied loadings result in a tangential displacement discontinuity [[u t ]], with a non-zero mean value. The selection of such a loading is still an open problem. We only observe that formula (12) allows us to check easily if the hypothesis is satisfied by a given loading system. Define the unit vector

U = [[u]] [[u]]
.

Equation (10) shows that the identification of the normal to the crack plane, by the RG functional, consists of the identification of two unit vectors N and U using the data R:

R = R 2R 2 -(Tr R) 2 .
Proposition 2. The two non-vanishing principal values of the tensor R are:

λ 1 = 1 + Tr R 2 and λ 2 = 1 -Tr R 2 , Tr R 1.
All possible vectors U and N can be chosen by the permutation and the change of signs in the following pairs of vectors, in the basis of eigenvectors

( 1 , 2 , 3 ) of R: λ 1 , -λ 2 , 0 , λ 1 , --λ 2 , 0 . ( 14
)
Since the normal N can be restricted without loss of generality to pointing into an a priori chosen hemisphere, two systems of applied loads are necessary for the selection of the normal.

If two systems of loadings 1 and 2 lead to the tensors R 1 and R 2 , one then obtains:

N = 1 3 ∧ 2 3 | 1 3 ∧ 2 3 | , ( 15
)
with the condition that these systems of loadings must be discriminative for the crack detection (in the sense that the vectors 1 3 and 2 3 are not parallel). Having determined the normal N by the previous proposition, let (T , V , N) be an orthonormal direct basis. With X k , k = 1, 2, 3 as the new coordinates in the basis (O, T , V , N), the equation of the crack plane is given by X 3 + C = 0. Under the condition U t = 0, we have the following proposition.

Proposition 3 (identification of the crack plane). The constant C is given by:

|C| = 1 2R 2 -2 T r(R) 2 α=T ,V RG 2 (v α ) , ( 16
)
where the fields v α are defined by:

v 1 T = -X 2 1 /2E -νX 2 2 /2E + (2 + ν)X 2 3 /2E, v 2 T = νX 1 X 2 /E, v 3 T = νX 1 X 3 /E, (17a) v 2 V = -X 2 2 /2E -νX 2 1 /2E + (2 + ν)X 2 3 /2E, v 1 V = νX 1 X 2 /E, v 3 V = νX 1 X 3 /E. ( 17b
)
Remark. Results from the lemma and the above propositions are valid for edge cracks, when intersects ∂ . It must be emphasized that unilateral contact problems, with or without friction on the crack surfaces, fall without restriction in the framework of the above results. The only additional difficulty lies in the a priori choice of external loadings giving rise to non-negative normal crack displacement discontinuity and to tangential displacement blockaded by friction.

Complete identification

In this section, we give a constructive method for proving the complete identification of the crack, once the crack plane has been determined. Again, the right tool is the linear form RG. Based on lemma 1, the identification of is obtained by the interpretation of [[u]] as a continuous linear form of L 2 ( ) (Riesz's representation theorem [START_REF] Brésis | Analyse fonctionnelle[END_REF]). The proof consists of two steps. In the first step, we show that the support of the displacement discontinuity coincides with the domain occupied by the crack (in other words, the displacement discontinuity cannot vanish in a null open set of ). In the second step, we construct two new families of fields v which enables us to identify the Fourier transform of the displacement discontinuity [[u]] .

Lemma 2. The support of [[u]

] coincides with the domain occupied by the crack:

supp[[u]] = .
Proof. Clearly one has supp[[u]] ⊂ . Suppose that the inclusion is strict. The proof is based on the statement that the discontinuity must be identically equal to zero, so that there is a contradiction with what has been assumed, that is that the loading leads to a jump of the tangential displacement with a non-vanishing mean value.

Suppose that supp[[u]] ⊂ . Then there is an open subset I in where:

[[u]] = 0 σ • N = 0.
Step 1. Denote by S 1 the intersection of and , and consider u 0 , the trace of the displacement u on S 1 :

u 0 = u S 0 + u A 0
where u A 0 is the normal component of u, u S 0 its tangential component and S splits into two subsets denoted by + and -.

Consider, now, the following elastostatic problem defined on + :

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ div σ + 1 = 0 in + σ + 1 • n = 0 on ∂ + ∩ ∂ u + 1 = u S 0 on S 1 σ + 1 = Aε(u + 1 ) and ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ div σ + 2 = 0 in + σ + 2 • n = F m on ∂ + ∩ ∂ u + 2 = u A 0 on S 1 σ + 2 = Aε(u + 2 )
. By the superposition principle, one has:

σ + = σ | + = σ + 1 + σ + 2 and u + = u | + = u + 1 + u + 2 .
Let S be the symmetry with respect to the plane . One extends σ + 1 and σ + 2 to S( + ) by: u - 1 (SX) = Su + 1 (X), u + 2 (SX) = -Su + 2 (X). Denote by ( σi , ũi ) the fields:

σi = σ + i in + σ - i in S( + ),
where

σ - i = Au - i and ũi = u + i in + u - i in S( + ).
Therefore ( σi , ũi ) satisfies the equilibrium equations in + ∪ S( + ) ∪ .

Step 2. Let ( σ , ũ) be equal to ( σ1 + σ2 , ũ1 + ũ2 ) by the Almansi lemma (i.e. unique continuation property for the elastic system), ( σ , ũ) coincides with (σ, u) on

( + 1 ∪ - 1 ∪ ) ∩ , therefore [u] = 0 on .
Theorem. If contains a plane crack not necessarily connected and if the loads F m are chosen in such a manner that [[u]] = 0, then is determined uniquely by the boundary data (F m , u m ).

Proof.After lemma 2, identifying is equivalent to identifying the support of [[u]]. As in [START_REF] Andrieux | Identification de fissures planes par une donnée de bord unique: un procédé direct de localisation et d'identification[END_REF] the reconstruction of [[u]] is based on the linear form RG(v) (with appropriate fields v of V ). The elastic case is more complicated and one has to identify the Fourier transform of [[u]] using two families of fields of V . After a new change of frame of reference (O , T , V , N) with O in the plane , one introduces the index ξ := (ξ 1 , ξ 2 , 0) running over R 2 and the family of vector fields Z ξ and Z * ξ of the complex space R 3 + iR 3 , with the family parameter ξ :

Z ξ = (ξ + i|ξ |N), Z * ξ = (ξ -i|ξ |N). (18)
Following a method introduced by Calderon [7], we construct the families of fields:

w + (X, ξ ) = ∇ exp(-iZ ξ • X) + ∇ exp(-iZ * ξ • X), (19a) w -(X, ξ ) = ∇ exp(-iZ ξ • X) -∇ exp(-iZ * ξ • X). ( 19b
)
These fields belong to the space V since the divergences of the associated stress fields σ + (X, ξ ) and σ -(X, ξ )) vanish. Applying lemma 1, we get:

RG(w + (ξ )) = 2E (1 + v) |ξ | 2 [[ û]] N (ξ ), RG(w -(ξ )) = - 2E (1 + v) i|ξ |ξ • [[ û]] t (ξ ),
where

  [[ û]] N and [[ û]] t denote respectively the Fourier transforms of the normal and the tangential components of the crack displacement discontinuity [[u]], with respect to the plane . Hence, we are let to the following equations:[[ û]] N = (1 + v) 2E|ξ | 2 RG(w + (ξ )), (20a)|ξ |ξ • [[ û]] t (ξ ) = (1 + v) 2iE RG(w -(ξ )),(20b)which show that the Fourier transform[[ û]] can be obtained uniquely by the boundary data (F m , u m ) since the parameter (index) ξ runs over R 2 . For example, if we take ξ = (ξ 1 , 0, 0) (respectively ξ = (0, ξ 2 , 0)) equation (20b) gives us the partial Fourier transform of [[ û]] t with respect to ξ 1 (respectively ξ 2 ). The support of the discontinuity and hence the crack itself, can be identified by the boundary data.


Teddy Fen-Chong

E ´veline

Andre ´dubault

Jean Louis Halary

Viscoelastic Characteristics of Pentane-Swollen Polystyrene Beads

Keywords: expanded polystyrene, swollen beads, pentane, plasticization, expansion, viscoelasticity

Expansion of pentane-swollen polystyrene beads was followed at and above the glass transition by using a very simple optical technique. Measurements allowed a suitable determination of the glass transition temperature and of the viscoelastic characteristics of the pentane-plasticized polystyrene chains inside the beads, at different pentane contents. This was achieved by combining the well-known Williams, Landel, and Ferry approach and some calculations derived from a recent micromechanical modeling of the expanded polystyrene microstructure. The results, analyzed in terms of free volume characteristics, revealed the peculiar plasticizing character of pentane, as compared to the usual polymer diluents. A value of 42°C was found as the lowest temperature for which the dimensional variation of the beads could result from glass transition motions. The relevance of these conclusions was discussed, not just in terms of the changes in dimensions of the individual beads at zero expansion time but also in terms of the dimensional evolution observed in the expanded polystyrene structures.

INTRODUCTION

For many years, solid polymeric foams have been the subject of extensive studies, dealing with both fundamental and technological aspects of synthesis, properties and applications (see, for example, ref. 1). One among the pending issues concerns the changes in dimensions that the foams can undergo after manufacturing.

The later point is especially crucial when dealing with expanded polystyrene (EPS), made from beads of pentane-swollen polystyrene. [START_REF] Cigna | [END_REF][3][4][5] In that case, after temperature has been raised above both the pentane boiling point and pentane-swol-len polystyrene glass transition temperature, the beads expand and acquire a cellular microstructure, as shown in Figure 1. Following this preexpansion, EPS slabs are produced by further expansion of a great number of preexpanded beads in a mold at temperatures between 100 and 115°C. After demolding, a temperature gradient exists across the cross-section of the slab until it uniformly reaches room temperature. At the end of the processing cycle, EPS is reported to still contain 2-3 wt % of blowing agent. [START_REF] Horak | Paper Three-Day Conf[END_REF] The concentration of pentane in the EPS decreases with time as the blowing agent diffuses and is replaced by air. From this moment, EPS slabs dimensions are reported to increase slightly for a short period before decreasing (shrinkage). This phenomenon goes on even after the temperature has reached constant room temperature. Until six weeks after complete cooling, it can be high enough to prevent EPS panels from being stuck onto plaster panels for building applications, [START_REF] Horak | Paper Three-Day Conf[END_REF][START_REF] Ja ¨rvela | To ¨rma ¨la ¨, P.; Ja ¨rvela ¨[END_REF] such that foamed polystyrene needs to be cured for reaching stabilization.

In order to reduce this costly and long-lasting cure time, it is necessary to know more on the influence of processing parameters on the dimensional evolution of EPS. To this end, it can be anticipated that two factors play a prominent role in the process: (1) the intracellular pressure, governed by the exchange between residual pentane and air across the cell walls; and (2) the mechanical resistance of the walls, associated with the modulus and glass transition temperature of the pentane-plasticized polystyrene.

In this paper, we present a novel way of determining the thermomechanical behavior of EPS beads as a function of temperature and pentane content, based on the optical inspection of the expansion of individual beads. As will be shown, analysis of the experimental data requires the combined use of the viscoelastic equation of Williams, Landel, and Ferry, 9 the so-called WLF equation, and of a recent micromechanical modeling of EPS cellular microstructure. [START_REF] Fen-Chong | [END_REF]

EXPERIMENTAL Materials Characterization

As-Received Beads

The pentane-swollen polystyrene beads were provided by Shell Co. (Research and Technology Center, Louvain-La-Neuve, Belgium). The diameter of these beads ranges from 0.9 to 1.3 mm. Beads of similar diameter were selected by inspection of the individual samples under a microscope equipped with a micrometric reticle. This timeconsuming procedure also allowed the elimination of some broken beads.

With the aim of characterizing the polystyrene that had been used for the bead preparation, the polymer was extracted from the as-received material. For this purpose, expandable beads were first dissolved in tetrahydrofurane to an amount of 5 wt % of beads. After 24 h under stirring, the solution was evaporated at 100°C under vacuum for 72 h in order to remove completely the volatile components (pentane and diluent). The dry solid residue was further dissolved in tetrahydrofurane at 0.1 wt % for gel permeation chromatography (GPC) analysis. It was found to consist of polystyrene (weight-average molecular weight around 185,000 g ⅐ mol Ϫ1 , polydispersity close to 2), and of a small amount (1%) of low-molecular weight additives (400 -500 g ⅐ mol Ϫ1 ). These unidentified additives could not be separated from polystyrene by precipitation: indeed, the latter done by putting a 0.1 wt % beads toluene solution in a large excess of methanol led to a solid exhibiting the same GPC chromatogram as the material obtained by evaporation. The glass transition of this unplasticized polystyrene T g 0 was determined by differential scanning calorimetry (DSC) at a heating rate of 10°C/min. T g 0 , taken at the onset of the change in heat capacity, was found to be 97 Ϯ 1°C.

Determination of the Amount of Pentane in the Beads

It is worth mentioning that the blowing agent, which we call pentane to be consistent with the industrial usage, is not actually pure n-pentane but a n-pentane-rich mixture of different isomers. As a first approximation, we will consider in this study that the properties of the blowing agent are identical to those of pure n-pentane.

Firstly, the amount of pentane was determined in the as-received beads, so-called beads A, by using the GPC technique. For this purpose, the expandable beads were dissolved in tetrahydrofurane (THF) at about 5 wt %. Direct evaluation of the pentane/polystyrene ratio would be imprecise due to the huge difference between the areas of the relevant peaks, as shown in Figure 2. This difficulty was overcome by analyzing ternary systems of bead-PMMA-THF (PMMA: polymethylmethacrylate) instead of the bead-THF mixtures (Fig. 2). For this purpose, a weighted amount of a reference material, a PMMA of 1830 g ⅐ mol Ϫ1 , was added to a weighted amount of beads in the solution, in order to yield comparable GPC peak areas for PMMA and pentane. The amount of pentane in the beads was derived using a calibration curve, relating the ratio Mass (Pentane)/ Mass (PMMA) to the ratio Peak Area (Pentane)/ Peak Area (PMMA), and determined from GPC measurements on additional n-pentane-PMMA-THF solutions of known composition (Fig. 3).

Following this procedure, the weight fraction of pentane in the as-received beads (beads A) was found to be around 8 wt %. Two additional series of beads, exhibiting a lower pentane ratio, were prepared from the as-received beads. In the first procedure, the original beads were placed in a vacuum oven (1 mm Hg) for 65 h at room temperature before testing (beads B). The weight fraction of pentane in these beads was found to be reduced to 6 wt %. A more drastic pentane removal was achieved by following another procedure. The original beads were placed in the vacuum oven for 120 h at room temperature and then for 40 h at 40°C, and finally they were left for equilibration for 90 h at 67°C and atmospheric pressure (beads C). The weight fraction of pentane in these beads was shown to drop to 1 wt %.

Bead Expansion Visualization Setup

Figure 4 shows the schematic diagram of the experimental setup. It consists of a thermally controlled Mettler hot-stage whose specimen compartment was modified to include a very small bath of silicone oil. One bead at a time was then plunged into this bath by means of a needle attached to a transparent glass slide. Above this slide was a microscope with a micrometric reticle that permits measurement of the diameter as the bead expands. Measurements were performed by eye in real time or by attaching to the binocular either a film or a video camera. In the former case, the diameter was measured after each 15 s in the first minute of expansion, then after each minute in the following 9 min, and finally after each 2 min until 20 min had gone by.

Figure 5 shows the expansion of one bead: it is worth noting that the bead remained spherical at any time, which means that the expansion proceeded, as desired, in a highly isotropic fashion.

RESULTS AND DATA ANALYSIS

The primary purpose of the experimental collection of the plots of bead diameter vs. expansion time at different temperatures (Fig. 6) was to determine the glass transition temperature T g (), of the pentane-swollen polystyrene, as a function of initial pentane weight fraction . Due to the plasticizing effect of pentane, one may anticipate that T g () should be lower than T g 0 and should decrease with increasing initial pentane content in the beads. In addition, the rate of bead expan- sion decreases with decreasing expansion temperature. For each series of pentane-swollen polystyrene beads, T g () is conventionally defined as the highest temperature for which no expansion was observed for 20 min under our setup. As shown in the second entry of Table I, the value of glass transition determined in this way decreases markedly with increasing pentane amount. By the way, it is worth pointing out that this procedure of the determination of T g is much more accurate in expandable materials than the familiar DSC technique. The latter, indeed, suffers from severe limitations in this case, because the heat flow signal is disturbed at and above T g by the changes in contact area between the specimen and the detector, as the result of expansion.

Second, the experimental data (Fig. 6) can be used to extract the expansion from the slopes of the isothermal curves of bead diameter vs. time. Analysis of the curves of expansion rate vs. time yields the initial expansion rates v 0,T . Such quantities depend on the expansion temperature and also on the initial pentane content of the beads.

Analysis of v 0,T values requires some theoretical assumptions with the aim of accounting for the bead mechanical behavior at the beginning of expansion. Let us first consider a foam with closed cells like an expanded bead as shown in Figure 1. This body is assumed to be at constant and uniform temperature with a traction-free boundary. The porosity is also supposed to be constant and uniform in the whole foam. As indicated in Figure 7, we consider the mechanical problem of predicting its changes in dimensions due to the presence of an intracellular pressure that is assumed to be identi-cal in each cell. The mechanical behavior of the cell membranes (the matrix) is assumed to be linear, isotropic, and elastic. It can be shown [START_REF] Fen-Chong | [END_REF] that the elastic strain of such an isolated bead is given by

,T ϭ c 1 Ϫ c ⅐ 1 Ϫ ,T 2 ⅐ ͑1 Ϫ 2 ⅐ ,T ͒ ⅐ 1 k ,T ⅐ P ,T (1)
where ,T is the elastic strain undergone by an isolated bead at temperature T, c is the porosity of the foam (volume fraction of cells), k ,T and , T are the bulk modulus and the Poisson ratio of the matrix, respectively, and P ,T is the intracellular pressure in excess of the atmospheric pressure.

All of these quantities depend on the temperature T and on the pentane content of the bead. If the cell membranes are supposed to present a linear and isotropic viscoelastic behavior, and if the Poisson ratio and the pentane content are assumed to be constant, then by means of the correspondence principle, 11 eq. ( 1) can be converted to

,T ͑t͒ ϭ c 1 Ϫ c ⅐ 1 Ϫ 2 ⅐ ͑1 Ϫ 2 ⅐ ͒ ⅐ J ,T ͑t͒ d dt P ,T ͑t͒ ϭ c 1 Ϫ c ⅐ 1 Ϫ 2 ⅐ ͑1 Ϫ 2 ⅐ ͒ ⅐ ͭ ͵ Ϫϱ t J ,T ͑t Ϫ u͒ ⅐ d du P ,T ͑u͒ ⅐ du ϩ i J ,T ͑t Ϫ t i ͒ ⅐P ,T i ͮ (2)
in which t is time, R is the Riemann-convolution product, and J ,T (t) is the creep matrix bulk com- pliance [convolution-inverse of the pentane-swollen polystyrene bulk relaxation modulus k ,T (t)], and ʈP ,T ʈ i is any pressure discontinuity occurring at time t i . Equations ( 1) and ( 2) were derived to predict the strain of an expanded bead between two relatively close states. This can be extended for an initial state corresponding to the beginning of the expansion. In this way, is the initial pentane wt % content of the bead and the intracellular pressure can be considered equal to the vapor pressure of pentane at temperature T, P vs (T). It is assumed that the temperature of each expandable bead instantaneously reaches the temperature T of the silicone oil bath, such that the intracellular pressure is proportional to the Unit-Step function H(t): P ,T ͑t͒ ϭ P vs ͑T͒ ⅐ H͑t͒.

(

By deriving eq. ( 2) with respect to time, and by substitution of eq. ( 3) for P ,T , we obtain the expression of the bead initial expansion strain rate v 0,T at t ϭ 0:

v 0,T ϭ ͫ d dt Tͬ tϭ0 ϭ c 1 Ϫ c ⅐ 1 Ϫ 2 ⅐ ͑1 Ϫ 2 ⅐ ͒ ⅐ ͫ d dt J ,Tͬ tϭ0 ⅐ P vs (T) (4)
It is worth noting that eq. ( 4) is not applicable to large times of beads expansion because in this case the intracellular pressure is no longer constant and the matrix behavior becomes elastovisco-plastic. As shown by eq. ( 4), consideration of the initial expansion strain rates at different temperatures for identical series of beads with a given wt % pentane content gives a means of deducing the type of (thermo)mechanical behavior of the cell membranes as a function of initial pentane content. For instance, if the time-temperature equivalence principle is assumed above the glass transition of pentane-swollen polystyrene, the translation factor 9 can be easily calculated: where a T g ( ) is the shift factor of the WLF equation between the temperature T and T g , taken as the reference temperature. Comparing eq. ( 5) with the experimental results needs the relative vapor pressure to be known. This is done by using literature data for the absolute vapor pressure of pentane [START_REF] Brandup | Polymer Handbook[END_REF] and by subtracting the atmospheric pressure (about 0.1 MPa). In such calculations, the hydrostatic pressure of silicone oil (about 44 Pa) is negligible. The temperature dependence of the intracellular pressure of pentane is given by the phenomenological Dupre ´'s law [START_REF] Brandup | Polymer Handbook[END_REF] :

v 0,T P vs ͑T͒ ϭ ͓͑d/dt͒J ,T ͔ tϭ0 ͓͑d/dt͒J ,Tg͑͒ ͔ tϭ0 ⅐ v 0,Tg͑͒ P vs ͓T g ͔͑͒ ϭ a Tg͑͒ ͑T͒ ⅐ v 0,Tg͑͒ P vs ͓T g ͔͑͒ , (5)
log͓P vs ͑T͔͒ Ϸ Ϫ6.5 ϩ 30 T ϩ 3 ⅐ log͑T͒, ( 6
)
where T holds for the temperature in Celsius and log is the decimal logarithm.

The next step of our data analysis was to consider the temperature dependence of the initial strain rates v 0,T at a given pentane content. The temperature dependence of the translation factor can be inferred by fitting the quantity log[v 0,T / P vs (T)] Ϫ log{v 0,T g () /P vs [T g ()]},which equals log[a T g () (T)] according to eq. ( 5). Since T g () has been defined as the highest temperature for which no expansion was observed for 20 min, it can be determined according to the precision of our optical measurements (half a micrometer graduation, i.e., 10 m) and to the mean expandable beads diameter (about 1 mm). It turns out that v 0,T g () is equal to about 4.8 10 Ϫ4 min Ϫ1 whatever the class of beads (A, B, or C) under consideration. Figures 8(a-c) show the temperature dependence of the normalized bead initial expansion rate for beads A, B, and C, respectively. If the capability of a bead to expand is assumed to be controlled, in the initial stages, by the ability of the polymer chains to develop glass transition motions, the temperature dependence of the initial expansion rate, or more precisely, of the quantity log[a T g (T)] can be analyzed by using the well-known WLF equation 9 :

log͓a Tg ͑T͔͒ ϭ log ͫ v 0,T P vs ͑T͒ ͬ Ϫ log ͫ v 0,Tg P vs ͑T g ͒ ͬ ϭ C 1 g ⅐ T Ϫ T g C 2 g ϩ ͑T Ϫ T g ͒ . (7)
As shown in Figure 8(a-c), the fits of the experimental data by eq. ( 7) are satisfactory in all cases. The values of the viscoelastic coefficients of the WLF equation, C 1 g and C 2 g taken at the reference temperature T g , are given in Table I, together with the critical glass transition temperature of the pentane-swollen polystyrene at infinitely long observation time T ϱ (), defined as

T ϱ ͑͒ ϭ T g ͑͒ Ϫ C 2 g ͑͒. (8)
Below this temperature, the glass transitionlike molecular motions cannot be observed even if the experimental time window were be infinite. T ϱ () does not seem to vary with the pentane content of the beads within the error bar on the determinations (about Ϯ5°C).

At this stage of our investigations, a pending question concerns the possibility of reducing the curves of Figure 8(a-c) to a unique master curve that should account for the initial pentane content of the beads. Simple replacement in eq. ( 5) of c/(1 Ϫ c)b y PS / C5 ⅐ /(1 Ϫ ), where PS ϭ 1050 g ⅐ dm Ϫ3 and C5 ϭ 626.2 g ⅐ dm Ϫ3 are the densities of pure polystyrene and of pentane, [START_REF] Brandup | Polymer Handbook[END_REF] respectively, would be misleading because it supposes that the whole pentane population is only in the cells; this assumption is in contradiction with the plasticizing effect that we have evidenced for pentane in the form of a polystyrene T g depression.

The viscoelastic coefficients C 1 g and C 2 g can be discussed in the framework of the free volume models by assuming the validity of two simple concepts 9 :

1. At T Ն T g , the viscosity depends on the fractional free volume according to the Doolittle equation:

ln ͑T͒ ϭ ln A ϩ B f͑T͒ , (9)
where f(T) is the fractional free volume available at temperature T, A and B are constants, the latter being usually taken to be unity.

2. The fractional free volume above T g , f(T), increases with temperature according to

f͑T͒ ϭ f͑T g ͒ ϩ ␣ f ⅐ ͑T Ϫ T g ͒, (10)
where f(T g ) is the fractional free volume available at T g , and ␣ f the free volume thermal expansion coefficient above T g . These quantities are related to the WLF coefficients by

f ͓T g ͔͑͒ ϭ B 2.3 ⅐ C 1 g ͑͒ and ␣ f ͑͒ ϭ f ͓T g ͔͑͒ C 2 g ͑͒ . ( 11
)
The values of f [T g ()] and ␣ f () are also given in Table I. It turns out that both fractional free volume available at T g and thermal expansion coefficient of the free volume above T g are found to drasti-

DISCUSSION

The above results need to be discussed in two separate ways, depending on whether emphasis is laid on the phenomenon of polystyrene plasticization or on the issue of the dimensional variations of EPS after manufacturing.

As far as plasticization of polystyrene in the beads is concerned, let us focus on the variations of T g (), f [T g ()], and ␣ f () as a function of the pentane weight fraction (Table I). Even though there is no proof that all the pentane molecules are in interaction with the polystyrene segments and actually act as a plasticizer, the weakly scattered plots shown in Figure 9 demonstrate that the observed evolutions make sense, at least from a qualitative viewpoint. The increase by a factor of 2.4 of the free volume available at T g when T g is depressed from 97 down to 65°C looks reasonable, and in agreement with the order of magnitude of f(T g ). [START_REF] Ferry | Viscoelastic Properties of Polymers[END_REF] More striking is the observed huge increase in ␣ f even for as small as 6 wt %. To our knowledge, such an effect has never been reported in the literature. However, the usual diluents (as those used to modify the mechanical characteristics of polyvinylchloride) are bulky and low vapor pressure molecules. One may suggest that the peculiarity observed here would be the result of the numerous motional degrees of freedom available in the flexible molecule of pentane over the temperature range of interest. As a direct consequence of the large sensitivity of ␣ f on , C 2 g () drastically decreases and T ϱ () becomes closer and closer to T g () when is increased. At the same time, T g () decreases markedly; hence, according to eq. ( 8), T ϱ () is roughly independent of .

The next point of interest deals with the value found for T ϱ , i.e., around 42°C. This result means that the relaxation associated with the glass transition motions of polystyrene chains in the beads is unlikely to occur below 42°C. This conclusion holds whatever the actual expansion degree of the beads is, since T ϱ was shown not to depend on the residual pentane amount. On the other hand, it does not mean that at room temperature molecular relaxation is impossible: it might occur on a very large time scale, provided some short spatial scale motions, analogous to those evidenced in low temperature creep experiments, are liable to develop.

Transposition of the above conclusions, relative to individual EPS beads, to the case (of industrial interest) of manufactured structures could be achieved only with prudence. In particular, one should bear in mind that the thickness of the cell walls is very much smaller than the diameter of the beads. As a consequence, the viscoelastic behavior of the plasticized polystyrene chains might be affected by the passage from volume to surface properties, leading to possible changes in the values of T g and T ϱ .

CONCLUSION

It is clear that the present study was not able to settle definitely the question of dimensional evolution in expanded polystyrene foams, a problem which has had much concentrated effort in the literature, especially from a technological viewpoint (see, for instance, refs. 1 and 4). As a matter of fact, the present contribution aimed at providing an insight in the relevant physicochemical phenomena. With this respect, three main objec-tives have been reached: (1) a novel methodology has been developed for studying the blowing of individual EPS beads at and above their glass transition temperature; (2) an original method for the data analysis has been derived, based on the combination of a very recent micromechanical modeling with the well-known WLF equation; and (3) quantitative results have been obtained in this way regarding the plasticization of polystyrene by the residual blowing agent.

This approach may be of interest for future work on polymeric foams, including the detailed behavior of EPS industrial structures and, tentatively, that of extruded polystyrene foams.

Figure 1

 1 Figure 1 Cellular microstructure of a section of an isolated expanded polystyrene bead (scanning electron micrograph): the cell membranes appear in white or grey.

Figure 3

 3 Figure3GPC calibration curve for pentane with PMMA as reference. The linear slope is found to be 0.98.

Figure 2

 2 Figure 2 GPC chromatogram of dissolved expandable beads with PMMA as the reference material.

Figure 4

 4 Figure 4 Schematic diagram of the experimental apparatus for the bead expansion experiments.

Figure 5

 5 Figure 5 Observation of the expansion of an isolated bead with an optical microscope. The initial bead is at the top and the expanded bead at the bottom. The end of the thermocouple can be distinguished in the right corner of the cell.

Figure 6

 6 Figure 6 Bead diameter vs expansion time. Example of 8 wt % pentane-swollen beads at various temperatures above T g .

Figure 7

 7 Figure 7Micromechanical analysis of an intracellular-pressurized foam with closed cells. Medium 1 is a cell with its membrane, medium 2 is made up of an infinity of medium 1. The considered problem consists of deriving the overall strain of the foam with linear, isotropic, and elastic cells membranes, subjected to a traction-free outer boundary.

Figure 8

 8 Figure 8 Experimental data treatment by WLF analysis: (a) 8 wt % pentane swollen PS, (b) 6 wt % pentane-swollen PS, and (c) 1 wt % pentane-swollen PS.

Figure 8 (Figure 9

 89 Figure 8 (Continued from previous page)

Table I WLF Parameters of Pentane-Swollen Polystyrene (The WLF Coefficients for Pure Polystyrene are Extracted from Ref. 13)

 I

	(%)	T g (°C)	C 1 g	C 2 g (°C)	f(T g )/B (%)	␣ f /B (10 Ϫ4 °CϪ1 )	T ϱ (°C)
	09 7 Ϯ 1	13.4 Ϯ 0.2	55 Ϯ 2	3.24 Ϯ 0.05	5.9 Ϯ 0.3	42 Ϯ 3
	1 Ϯ 0.5	88 Ϯ 1	8.3 Ϯ 0.2	50 Ϯ 2	5.24 Ϯ 0.13	10.5 Ϯ 0.7	38 Ϯ 3
	6 Ϯ 0.5	72 Ϯ 1	6.1 Ϯ 0.2	32 Ϯ 2	7.13 Ϯ 0.23	22.3 Ϯ 2.1	40 Ϯ 3
	8 Ϯ 0.5	65 Ϯ 1	5.6 Ϯ 0.2	18 Ϯ 2	7.76 Ϯ 0.28	43.1 Ϯ 6.3	47 Ϯ 3

We thank Thavarajah Shanmuganathan (a student from the University Paris XI-Orsay) for his contribution to the experiments and Dr. T. Lalot (Laboratoire de Synthe `se Macromole ´culaire, University Paris VI) for his help in the GPC measurements. We are also indebted to Prof. A. Zaoui (E ´cole Polytechnique) and Prof. L. Monnerie (ESPCI) for kind and stimulating discussions.


V D

A FIRST STAGE IN THE DEVELOPMENT OF MICROMECHANICAL SIMULATIONS OF THE CRYSTALLOGRAPHIC PROPAGATION OF FATIGUE CRACKS UNDER MULTIAXIAL LOADING

Keywords: Stage I crack, Fatigue, Mode II, Multiaxial, Simulations, Dislocations

Simulations of the nucleation of dislocations, glide and annihilation ahead of a fatigue crack growing along a localized slip band (a 'long' Stage I crack or a Stage II crack with a K value close to the threshold) are performed for the case of push-pull or reversed torsion loadings, ignoring, in a first approach, the effect of grain boundaries. The crack growth rates are deduced from the dislocation flux at the crack tip. An influence of the normal stress on the friction between the crack flanks as well as on the condition for dislocation emission is introduced. A slower Stage I growth rate is then predicted for reversed torsion, consistent with experimental data.

INTRODUCTION

It is widely acknowledged that reversed shear strain is the driving force for Stage I fatigue crack growth, but that the presence of an opening stress on highly sheared facets can assist the development of microcracks [START_REF] Tong | Kinetic study on the propagation of fatigue crack in pure aluminium single crystals under different modes of loading[END_REF][START_REF] Gell | The characteristics of Stage I fatigue fracture in a high strength nickel alloy[END_REF][START_REF] Tanaka | Fatigue crack growth along planar slip bands[END_REF]. The absence of such an opening stress in the case of torsional loading is often considered responsible for slower Stage I growth and longer fatigue life compared with equivalent strain ranges in push-pull. As a consequence, most, if not all, multiaxial fatigue criteria [START_REF] Brown | A theory for fatigue failure under multiaxial stress-strain conditions[END_REF][START_REF] Fatemi | A critical plane approach to multiaxial fatigue damage including outof-phase loading[END_REF] incorporate both the maximum shear strain range and the stress (or strain) normal to the facets that undergo the maximum shear strain range. The relative magnitude of the second term has to be adjusted through materials constants so as to unify the tensile and torsional data. As was underlined in a previous paper [START_REF] Doquet | Crack initiation mechanisms in torsional fatigue[END_REF], these material constants are extremely variable from one material to another, but we do not have any satisfactory explanation for this variability, because the effect the normal stress has during Stage I is not clear, from the point of view of the physical mechanism.

The aim of the present work is thus to suggest mechanisms by which the normal stress could influence Stage I fatigue crack development or the crystallographic propagation of a Stage II crack close to the threshold, and to represent them by micromechanical simulations of dislocation glide in the plastic zone, and then to try to derive some of the microstructural or mechanical parameters that could determine the influence of the normal stress.

The simulations are similar to those developed by Pippan [START_REF] Pippan | Dislocation emission and fatigue crack growth threshold[END_REF], or Wilkinson and Roberts [START_REF] Wilkinson | A dislocation model for the two critical stress intensities required for threshold fatigue crack propagation[END_REF]. The difference is that the present study incorporates the effect of the normal stress on friction forces along crack flanks as well as on the condition for dislocation nucleation.

PHILOSOPHY OF THE SIMULATIONS

Stage I fatigue cracks grow along localized slip-bands. They are submitted to either pure shear loading (mixed mode II+III) in torsion, or combined opening and shear in push-pull. But plasticity, i.e. dislocation nucleation and glide, is only related to the shear components and is restricted to the slip band co-linear to the crack, otherwise there is a transition towards Stage II propagation (if non-coplanar slip is activated) [START_REF] Koss | Fracture along planar slip bands[END_REF]. As a consequence, and unlike the shear stresses, the singular opening stress is not relaxed (or shielded) by the dislocations stress field [START_REF] Lin | Cleavage, dislocation emission, and shielding for cracks under general loading[END_REF].

In the present simulations, only mixed mode I+II loading has been considered for the sake of simplicity. This means that the dislocations that are emitted and glide along the coplanar slip plane have a pure edge character (instead of mixed edge and screw when a mode III component is considered). They are supposed to be perfect dislocations. However, the possibility of cross-slip is not envisaged. This means that the present simulations are pertinent essentially for medium stacking-fault-energy F.C.C. metals and alloys.

Dislocations are supposed to be emitted from the crack tip as has often been observed during in situ tests performed in a transmission electron microscope [START_REF] Ohr | An electron microscope study of crack tip deformation and its impact on the dislocation theory of fracture[END_REF].

No major microstructural obstacle to dislocation glide is considered as a first approach. This will be introduced in a forthcoming paper, since the interactions of the crack with the microstructure are an essential feature of Stage I growth. It is assumed that the tip of the plastic zone associated with the crack tip does not encounter any grain or twin boundary.

The geometry considered here (a straight, two-dimensional crack in an infinite half space) is thus very simple compared to the real geometry of Stage I cracks, but it is believed that such a simple model can, however, provide some insight into the crack growth mechanisms.

The author is aware of the limitations of Linear Elastic Fracture Mechanics in the case of short cracks. But, as in the work of Tanaka and Mura [START_REF] Tanaka | Fatigue crack growth along planar slip bands[END_REF] who consider 'long stage I cracks', though still in the size range influenced by the microstructure, stress intensity factors will be used in the following. However, a 'correction for plasticity' will be introduced, as shown below, based on the stress field of dislocations [START_REF] Lin | Cleavage, dislocation emission, and shielding for cracks under general loading[END_REF]. This may be acceptable for torsional fatigue, where Stage I cracks can, in many instances, even in high-cycle fatigue, reach lengths that justify the use of LEFM before branching [START_REF] Pinna | Etude de la propagation des fissures de fatigue sous chargement de cisaillement plan[END_REF]. This can also be advocated for cyclically softening materials, e.g. Inconel 718, where Stage I cracks also reach a consequent size, even in push-pull [START_REF] Socie | Fatigue damage maps[END_REF], or for high strength steels with a fine microstructure (e.g. 4340 steel ) where, according to Lankford [START_REF] Lankford | On the small crack fracture mechanics problem[END_REF] or Smith [START_REF] Smith | Short fatigue cracks. In: Fatigue Mechanisms: Advances in Quantitative Measurement of Physical Damage[END_REF], the minimum crack length can be less than 100 mm, and yet careful use of LEFM parameters may yield reasonable predictions. In any case, this should not be a problem in describing crystallographic propagation of Stage II cracks when DK is close to the threshold, and as concerns Stage I cracks, this can be considered as an assumption whose relevance has to be evaluated with regard to the physical soundness of the resulting predictions.

Concerning the mode II stress intensity factor, allowance has to be made for the friction forces distributed along the crack flanks that make the effective stress intensity factor, Keff II , different from its nominal value, Knom II [START_REF] Pinna | Quantification of friction effects during mode II fatigue crack propagation in a maraging steel[END_REF]. These friction forces exist, even in pure mode II, because of the asperities of the crack flanks that are formed either when the crack adopts a zig-zag path between several parallel slip bands in the same grain, or when it crosses a grain boundary and is thus slightly tilted because of a difference in crystallographic orientation. The reality of these friction forces is demonstrated in Fig. 1 by the wear debris coming out of the mouth of a Stage I crack photographed during a combined tension and torsion test carried out in a SEM on a tubular specimen of Co45Ni. Lynch and Ryder [START_REF] Lynch | The fatigue behaviour of two aluminium-zinc-magnesium alloys[END_REF] have shown how much Stage I growth rate in an aluminium alloy submitted to torsional fatigue is affected by a modification of the tribological conditions: in the presence of an inert fluid that tends to remove the fretting debris from the crack flanks, a four-fold increase in Stage I crack growth rate is observed, whereas Stage II crack growth rate is decreased.

Friction forces will of course be intensified in the presence of a normal compressive stress and reduced if there is an opening stress. Let us now consider that n edge dislocations with Burgers vector b i =±b (depending on their sign) are present in the plastic zone. Considering the shear stress produced by each dislocation along the plane of the crack, it can be shown [START_REF] Lin | Cleavage, dislocation emission, and shielding for cracks under general loading[END_REF] that they have a shielding effect on the stress distribution at the crack tip, characterized by

Ktip II =Keff II -∑ n i=1 mb i (1-n)ǰ2px i (1)
where x i stands for the distance of the i-th dislocation from the crack tip. m and n are the shear elastic modulus and Poisson's ratio, respectively.

The criterion for dislocation nucleation at a crack tip was recently analysed by Rice [START_REF] Rice | Dislocation nucleation from a crack tip: an analysis based on the Peierls concept[END_REF] on the basis of the Peierls concept. A tension-shear coupling equation, based on atomic models, was later established by Sun et al. [START_REF] Sun | Estimates from atomic models of tension-shear coupling in dislocation nucleation from a crack tip[END_REF]. According to them, nucleation occurs if Ktip II reaches a critical value, Knucl II , with:

Knucl II = S 2m (1-n) [cr us -a(cu us -cr us ) A p 2 -y BD (2)
where y=arctan(K II /K I ) is the phase angle between mode I and II loadings. (In torsion, a pure shear case, y=P/2, whereas for push-pull, at any time, the normal and shear stresses are equal in magnitude, but since K I is zero when the normal stress is compressive, y=P/4 during the tensile phase and P/2 in the compressive phase, see Fig. 2). a is a positive constant that depends on the crystallographic and electronic structure of the material. c us is termed 'the unstable stacking energy'. It is defined as the energy increase per unit area of slip plane when the lattice on one side of the plane is shifted relative to the lattice on the other side, to the unstable equilibrium position at or near to a sliding displacement of b/2. Now, this can be achieved in two ways, either letting a certain amount of opening displacement occur (this is the relaxed case) or without any opening displacement (this is the unrelaxed case). The latter case corresponds to a value cu us of c us , which is higher than cr us in the previous case. Equation [START_REF] Gell | The characteristics of Stage I fatigue fracture in a high strength nickel alloy[END_REF] thus predicts a lower threshold stress intensity factor for dislocation nucleation when y decreases, i.e. when an opening stress is present.

In this paper, only fully reversed loadings are envisaged (the effect of the K min /K max ratio will be treated in a forthcoming paper), so that K II changes sign. As a consequence, negative as well as positive edge dislocations can be emitted when, respectively:

Ktip II ∏-Knucl II or Ktip II Knucl II (3)
Once emitted, these dislocations glide or rest according to the law:

n i =n 0 sign(b i )Ωsign(t i )Ω|t i |-t f m where x is zero if x∏0 and x otherwise (4)
where n i is the velocity of the i-th dislocation, t f is the lattice friction on the dislocations, n 0 and m, two constants. t i is the shear stress on the i-th dislocation, i.e. according to Ohr [START_REF] Ohr | An electron microscope study of crack tip deformation and its impact on the dislocation theory of fracture[END_REF]:

t i = Keff II ǰ2Px i - mb i 4P(1-n)x i -∑ j≠i mb j 2P(1-n) S x j x i Ω 1 x j -x i (5)
The second term on the right-hand-side of Eq. ( 5) represents the image stress on the dislocation due to the free surface of the crack, the third term is the shear stress due to the other dislocations (and their own image dislocations, hence the ǰx j /x i term). Annihilation between a positive and negative dislocation is allowed when their mutual distance is less than a critical value y e . A dislocation will return to the crack tip and then annihilate if it approaches closer than a critical distance y r . The cyclic loading path is followed by incremental time steps, Dt, small enough for the velocity of each dislocation to be considered constant over Dt. In practice, Dt is continuously adjusted so that, during one time step, dislocations move at most by a set fraction of the distance to their nearest neighbour's previous position.

Two successive cycles only are simulated because the second cycle is representative of the steady state.

The crack is considered to grow by one Burgers vector each time a pair of positive-negative dislocations has been emitted at the crack tip, or when a positive (or negative) dislocation returns to the crack tip. In the latter case, it is assumed that even though the crack tip geometry before the dislocation nucleation is, in principle, recovered when this dislocation comes back, the free surface increment created at nucleation, that has been exposed to the environment and hence gas adsorption, cannot be rewelded. Both events (positive-negative pair emission, or emission and return of a dislocation) correspond to some cyclic plastic flow at the crack tip and should thus contribute to its growth.

IMPLEMENTATION OF THE SIMULATIONS

The nominal stress intensity factors at the crack tip are calculated at each time-step as:

K I =s ncp Ωǰpa when s ncp 0 Knom II =t cp Ωǰpa (6)
where s ncp and t cp are the current values of the normal stress and shear stress on the critical plane (i.e. the plane where the shear stress range is maximum), respectively.

Experimental information concerning crack flanks frictional interactions for a long crack loaded in mode II has been obtained, for a maraging steel, in the framework of other studies [START_REF] Pinna | Etude de la propagation des fissures de fatigue sous chargement de cisaillement plan[END_REF][START_REF] Pinna | Quantification of friction effects during mode II fatigue crack propagation in a maraging steel[END_REF]. Precracked specimens, covered with square-shaped microgrids with a 5 mm pitch were loaded in shear inside the chamber of the scanning electron microscope, and the displacement profiles measured using microgrids as reference marks. It was shown that the measured displacements resulted from the superposition of the nominal load and a uniform friction stress, ranging from 20 to 70 MPa (depending on crack flanks roughness) in the absence of a normal stress, and reaching 110 MPa when a 150 MPa compressive stress was applied. In addition, a thin tubular specimen containing a transversal crack has been progressively loaded in torsion, after the application of a tensile or compressive stress, inside the SEM, and the threshold stress intensity factor for which slip was observable at the crack tip was measured. The results are plotted in Fig. 3 as a function of the normal stress, s n . It can be seen that the influence of a compressive normal stress on crack flanks friction is more significant than that of an opening stress: a 80 MPa opening stress seems sufficient to reduce the frictional resistance to crack flanks sliding to zero, but #230 MPa is necessary, in compression, to double this resistance.

In the present study, an attempt is made to reproduce these effects qualitatively through empirical equations.

The influence of the normal stress, s n , on the friction shear stress, c, distributed along the whole crack length, 2a, is expressed by:

c=c 0 exp(-k + Ωs n ) if s n 0 c=c 0 exp(-k - Ωs n ) if s n ∏0 ( 7
)
c 0 is a constant which characterizes the friction stress in the absence of any normal stress, and is thus related to the tortuosity of the crack path, which is generally less pronounced for a Stage I crack than for the long cracks involved in the study mentioned above. The order of magnitude of c 0 will thus be 5-10 MPa in the simulations. k + and k are two constants, the latter connected with the friction coefficient of the material. Most of the simulations will be performed with k + = 0.057, k -=0.014. (This means that c is divided by 100 for s n =80 MPa, but it is only multiplied by three for s n =80 MPa.) The value chosen for k leads to a stronger influence of s n on frictional effects than that measured in maraging steel, since in many FCC materials, crack flanks asperities must impinge more easily, crashed against one another by a compressive force, than in this high strength steel, thus leading to enhanced frictional effects. It follows:

Keff II =Knom II ±cǰPa (8)
the sign affecting the second term depending on the loading (-) or unloading (+) situation. Of course, when the loading changes direction, there will be periods where Knom II changes, but Keff II remains constant, until the frictional resistance to reverse crack flanks slip is overcome.

Figure 4 shows how Knom II and Keff II vary in time during a reversed torsion and push-pull cycle described, respectively, by: t=t 0 Ωsin vt, s=s 0 Ωsin vt, with s 0 =2t 0 (equivalent stress ranges in the sense of Tresca), assuming that the crack lies along the maximum shear strain plane in each case (i.e. longitudinal/transversal or at 45°to the tensile axis, respectively). Note that DKnom II is the same in both cases.

Sun et al. [START_REF] Sun | Estimates from atomic models of tension-shear coupling in dislocation nucleation from a crack tip[END_REF] 2) gives a nucleation stress intensity factor in the tensile phase of the push-pull cycle that is 92% of its value in torsion or in the compressive part of the cycle. (Note that the effect would be more pronounced for 90°out-of-phase tension and torsion loading-a case for which y varies continuously between 0 and P/2, and consequently, Knucl II can be as low as 84% of its value in pure shear. This case thus deserves a special study.) Knucl II for pure mode II was set to 0.3 MPaǰm in accordance with the order of magnitude given in Ref. [START_REF] Ohr | An electron microscope study of crack tip deformation and its impact on the dislocation theory of fracture[END_REF].

The critical distance for positive-negative dislocations annihilation, y e , was taken as 16 nm, which is the annihilation distance for edge dipoles in copper given by Essmann and Differt [START_REF] Essmann | The nature of the wall structure in persistent slip bands of fatigued metals[END_REF].

The critical distance for dislocation annihilation at the crack tip, y r , was chosen as 10-10 m, which is the order of magnitude of the core radius.

A lattice friction force t f of 20 MPa was chosen, which is close to the local shear stress measured in the channels of persistent slip bands in copper by Lepinoux and Kubin [START_REF] Lepinoux | In situ TEM observations of the cyclic dislocations behaviour in persistent slip bands of copper single crystals[END_REF].

The constants in the expression of dislocations velocity [Eq. ( 4)] are chosen as n 0 =13 ms-1 and m=0.88, which are typical experimental values for FCC metals [START_REF] Wilkinson | A dislocation model for the two critical stress intensities required for threshold fatigue crack propagation[END_REF].

RESULTS

Pure shear case (reversed torsion)

The sequence of events and the positions of the dislocations during two loading cycles in reversed torsion with Knom II,max =0.394 MPaǰm and a friction stress c 0 =5 MPa along the crack flanks are sketched in Figs 5 and6.

Nothing happens initially while Ktip II is less than 0.3 MPaǰm. Then, five positive dislocations are successively nucleated, producing, each time, a shielding effect that decreases quickly as the dislocation glides away from the crack tip [see the serrated evolution of Ktip II in Fig. 5( b)]. Note [Fig. 6( b)] the existence of a dislocation-free zone (D.F.Z.) between the crack-tip and the last emitted dislocation. This fact, consistent with experimental observations [START_REF] Ohr | An electron microscope study of crack tip deformation and its impact on the dislocation theory of fracture[END_REF], does not follow naturally from approaches based on continuously distributed dislocations that ignore image forces, as in the work of Bilby, Cottrell and Swinden [START_REF] Bilby | The spread of plastic yield from a notch[END_REF], but they can be introduced, as in the work of Weertman et al. [START_REF] Weertman | Double slip plane crack model[END_REF], who consider dislocations distributed along two slip planes parallel to the crack, one above and one below the crack plane, and 'assume a dislocation-free zone in front of the crack to correspond to the slip plane separation'.

Nothing occurs then during a large portion of the unloading phase because of the stationary value of Keff II on one hand [see Fig. 5(a)] and the lattice friction that delays the reversed motion of the dislocations on the other hand. The last emitted dislocation is the first to be unlocked and, as Keff II grows negative, it travels back to the crack tip, where its sudden antishielding effect [evidenced by a serration on the Ktip II curve of Fig. 5(c)] leads to the emission of a negative dislocation. Annihilation of the positive and negative dislocations follows, and Ktip II recovers its 'regular' value until the following positive dislocation stacked in the plastic zone travels back to the crack tip in its turn.

This process is repeated until the third emitted positive dislocation has been annihilated. The last two positive dislocations remaining in the plastic zone, however, do not travel back to the crack tip for two reasons: first, the back-stress due to preceding dislocations, that would favour a reverse motion, is weak for the second one and null for the first one; and second, the decrease of Keff II is slower and slower (sinusoidal nominal waveform), and perhaps not sufficient to overcome the lattice frictional resistance to dislocation glide.

The next negative dislocation nucleated comes to meet the second positive one and annihilate it. Then, two additional negative dislocations are nucleated before the first one approaches the last remaining positive one close enough to annihilate it. Two additional negative dislocations are nucleated before 'unloading', during which four of the five negative dislocations present in the plastic zone return to the crack tip where they promote the emission of positive nucleations and annihilate them, whereas the last one will be annihilated later by a positive dislocation coming to meet it.

After that, a similar sequence of events repeats itself and finally, the computed crack growth during the second period is 10 Burgers vectors.

Similar simulations have been carried out for various loading ranges and various friction stresses c 0 . The results are drawn in a bilogarithmic da/dN versus DKnom II plot in Fig. 7. It can be seen that the existence of a threshold stress intensity range below which the crack does not propagate for lack of cyclic plasticity at the tip, as well as the well-known shape of da/dN versus DK curves in the vicinity of this threshold are qualitatively reproduced by the simulations, as in the work of Wilkinson and Roberts [START_REF] Wilkinson | A dislocation model for the two critical stress intensities required for threshold fatigue crack propagation[END_REF], and Pippan [START_REF] Pippan | Dislocation emission and fatigue crack growth threshold[END_REF]. by the opening stress, in push-pull, whereas for higher loads, nucleation is no longer critical, but reverse dislocation glide, which is more natural under reversed torsion, is. According to this trend, the ratio between torsional and tensile fatigue lives should thus increase as the loading range decreases. This corresponds to experimental observations and is thus very encouraging [START_REF] Doquet | Crack initiation mechanisms in torsional fatigue[END_REF].

At this stage of the study, the mere effect of the normal stress on the nucleation threshold does not seem to be a sufficient explanation for the differences in fatigue lives usually observed between push-pull and reversed torsion. But this conclusion might change with the introduction of grain boundaries against which the dislocations emitted by the crack will pile up, thus generating stress concentrations in the neighbouring grain. The higher number of dislocations emitted in push-pull for an equivalent loading range should generate higher stress concentrations, thus promoting slip activation in the next grain and thus, probably, an easier crossing of the grain boundary.

Besides, when the influence of the normal stress on crack flanks friction is introduced, the difference in Stage I kinetics between the two loading modes is increased, and this difference is amplified when the friction stress c 0 increases (compare curves C and E, or D and F in Fig. 9).

CONCLUSIONS

As a preliminary stage towards a simulation of the crystallographic propagation of fatigue cracks, the crack development process by the emission/annihilation of dislocations at the crack tip has been simulated, ignoring, in a first approach, the effect of grain boundaries. The originality of these simulations lies in two aspects. First, in the distinction made between nominal and effective loading, allowance is made for the friction stresses that resist the sliding of crack flanks, these stresses being either decreased or enhanced by a normal stress; and second, in the incorporation of tension-shear coupling with respect to the critical stress intensity for dislocation nucleation, following Rice and co-workers.

The crack speed (da/dN) versus DK curves were obtained for push-pull and reversed torsion. The existence of a threshold stress intensity factor for crack growth, as well as the shape of the da/dN versus DK curves in the vicinity of the threshold are qualitatively reproduced. The influence of friction stresses along the crack flanks is also clearly illustrated: it shifts the da/dN versus DKnom II curve to higher DK values compared to the 'intrinsic' curve corresponding to zero friction. Due to the effect of a normal stress on the dislocation nucleation threshold, crystallographic growth is predicted to be slower under reversed torsion than under push-pull for equivalent (Tresca) stress (or strain) ranges, and the ratio between torsional and tensile fatigue lives is predicted to increase as the loading range decreases, which is consistent with experimental data on many materials. The difference in the crack growth rates calculated under the two loading modes is amplified when normal stress-dependent frictional effects are introduced.

If the mechanisms proposed here to explain the influence of the normal stress on Stage I fatigue crack growth were pertinent, the difference between tensile and torsional fatigue lives should thus depend on the roughness of Stage I cracks and the friction coefficient of the material, both affecting the intensity of the friction stress along the crack flanks. It should also depend on the crystallographic and electronic structure of the material, which determine the degree of tension-shear coupling in relation to dislocations emission. Very specific experiments have to be designed to check these points.

Some results of the simulations suggest that the interaction between a Stage I crack and a grain boundary should be different under push-pull and reversed torsion. Work is thus now in progress to complete the model by taking into account the existence of microstructural obstacles to dislocation glide (essentially, grain and twin boundaries). The possibility of a non-zero mean stress should also be introduced. Stage I crack growth under multiaxial non-proportional loadings deserves special investigation.

Fig. 1 .

 1 Fig. 1. Aspect of a Stage I crack formed in a Fig. 2. Nominal loading paths, in terms of stress intensity factors, Co45Ni alloy under 90°out-of-phase tension and for a crack lying along a maximum shear strain plane under reversed torsion loading. torsion and push-pull.

Fig. 4 .

 4 Fig. 4. Variation of nominal and effective mode II stress intensity factors for a crack lying along the maximum shear Fig. 3. Variation of the mode II stress intensity factor at stress plane under (Tresca) reversed torsion and push-pull the onset of glide, measured at the tip of a crack progress-loadings of equivalent amplitude. c 0 =5 MPa, k + =0.057, ively loaded in plane shear, in a maraging steel. k -=0.014.

Fig. 5 .

 5 Fig. 5. Results of K II values from the simulation in the pure shear case (reversed torsion). c 0 =5 MPa, t f =20 MPa, Knucl II =0.3 MPaǰm. (a) Variation of Knom II , Keff II and Ktip II with time. ( b), (c): details of (a).

Fig. 6 .

 6 Fig. 6. Results of the simulation in the pure shear case (reversed torsion) c 0 =5 MPa, t f =20 MPa, Knucl II =0.3 MPaǰm. (a) Number of dislocations in the plastic zone. ( b) Positions of the dislocations. (The dislocations are numbered in the chronological order of their emission. Those which annihilate others soon after their emission are mentioned only by number.)

Fig. 7 .

 7 Fig. 7. Crack growth rates calculated various friction stresses versus nominal DK II in the pure shear case (reversed torsion).

Fig. 9 .

 9 Fig. 9. Comparison of the calculated crack growth rates versus nominal DK II curves. (A) push-pull, no friction; (B) reversed torsion, no friction; (C) push-pull, c 0 =5 MPa; (D) push-pull, c 0 =10 MPa; (E) reversed torsion, c 0 =5 MPa; (F) reversed torsion, c 0 =10 MPa.

In addition, the influence of the friction stress along the crack flanks is illustrated: quite similarly to the closure effects for mode I propagation, the friction effects shift the da/dN versus DKnom II curve compared to the 'intrinsic' curve corresponding to zero friction. In the absence of quantitative data on the friction stresses along the flanks of a Stage I crack, a discussion on the calculated values of the threshold stress intensity factors would be premature, and what is more, experimental values of mode II thresholds in single crystals, that could be compared to the calculated values, seem to be rare. This threshold for cyclic plasticity at the crack tip is clearly less relevant for engineering applications than the threshold for grain boundaries crossing which is, in many instances, intimately connected with the fatigue limit.

In phase opening and shear case (push-pull) Figure 8 shows the results of a simulation performed for Knom II,max =0.394 MPaǰm and c 0 = 5 MPa (as in the previous case of Figs 5 and6, i.e. reversed torsion) in the case of push-pull loading.

It can be seen that the positive dislocations nucleated in the tensile phase are much more numerous than in the pure shear case (77 instead of five), because of a higher maximum effective stress (0.369 MPaǰm instead of 0.306 MPaǰm) and a lower nucleation threshold.

The truncation of the Keff II cycle in the compressive part of the cycle (Keff II,min =0.273 MPaǰm only, because of enhanced frictional effects) does not, however, allow all these dislocations to return to the crack tip, and 28 of them stay in the plastic zone. [There again, those which return to the crack tip are responsible for a sudden shielding that disappears when the dislocation annihilates, hence the serrated evolution of Ktip II in Fig. 8(a).] The crack growth rate per cycle (49 Burgers vectors) is nevertheless much higher than in the pure shear case.

Several such simulations were performed for push-pull with various loading ranges and friction stresses, c 0 . The results are compared to torsional results in Fig. 9 in the form of da/dN versus DKnom II curves. Comparing the curves labelled A and B, it appears that even without friction, the crack growth rates are higher in push-pull than in reversed torsion, but the difference, appreciable for small loading ranges, vanishes as the loading range increases. This is probably because close to the threshold, the critical stage for crack propagation is dislocation nucleation, which is made easier