V Aleshin

W Desadeleer

M Wevers

K Van

Den Abeele

K Van Den Abeele
email: koen.vandenabeele@kuleuven-kortrijk.be

Characterization of hysteretic stress-strain behavior using the integrated Preisach density

Keywords: hysteresis, rock mechanics, non-destructive testing, nonlinear constitutive behavior

In this paper, we introduce the concept of Integrated Preisach-Mayergoyz (IPM) density to analyze static uniaxial compression tests at values well below the critical strength, and to characterize the elasticity of materials with hysteresis in their stress-strain relationship. The IPM density can be deduced from a particular force protocol following basic data treatment. The advantage of the IPM density over prior approaches is that no second order differentiation of the data is required which reduces the errors and uncertainties typical for past practice in the specific context of rock elasticity using scanning curves and PM density analysis. The characterization of the elasticity of the material is established in terms of a non-hysteretic strain contribution in the form of a nonlinear but reversible equation of state, and a hysteretic contribution represented by the IPM density. The IPM inversion procedure is tested for simulated stress-strain data subjected to additive noise, and the results are compared to the traditional methodology. In addition, we analyze the hysteretic and non-hysteretic characteristics of five natural building stones, and show evidence for a classification based on the inferred properties.

Introduction

The class of materials for which hysteresis is believed to be responsible for the nonlinear elastic and acoustical properties, is quite extensive (rocks, minerals, construction materials, composites, ceramics, polycrystalline metals, martensites). An enlightening exemplar for observing the separation of loading and unloading curves in static stress-strain tests is the family of sandstones, where cyclic static tests under pressures of a few tens of MPa's reveal a difference between the ascending and descending branches (loop thickness) which is comparable to the elastic strain itself. A similar branch separation can be observed in the stress-strain behavior of fatigued metals. This makes hysteresis a material state property of primary importance and, consequently, the problem of its characterization and its inversion from measurements arises.

A traditional model for the description of rate-independent hysteresis is the Preisach space approach [START_REF] Preisach | Über die magnetische Nachwirkung[END_REF] which was initially developed for ferromagnetics. Later, this methodology was generalized in the perspective of control theory, including the proofs of several important theorems by Krasnosel'skii and Pokrovskii [START_REF] Krasnosel | Systems with Hysteresis[END_REF], whose results are discussed in a classical paper by Mayergoyz [START_REF] Mayergoyz | Hysteresis models from the mathematical and control theory points of view[END_REF]. In the 1990's, this formalism was applied for the elasticity of martensites [START_REF] Ortin | Preisach modeling of hysteresis for a pseudoelastic Cu-Zn-Al single crystal[END_REF] and geomaterials [START_REF] Mccall | Equation of state and wave propagation in hysteretic nonlinear elastic materials[END_REF], and in this context the term PM model or PKM model (or space) was adopted. Key player in this model is the PM-density distribution which represents the statistical distribution of microscopic units with a bistable state within a representative volume.

As far as material characterization is concerned, the basic methodology for inverting the hysteretic material characteristic (i.e. the PM density) from experimentally obtained stress-strain data was published by Guyer et al. [START_REF] Guyer | Quantitative implementation of Preisach-Mayergoyz space to find static and dynamic elastic moduli in rock[END_REF]. It was noted that the PMdensity inversion is a highly underdefined problem: the number of constraints that can be obtained from a few hysteretic loops is typically much smaller than the number of points defining the PM density (2D surface). To avoid this difficulty, one generally uses regularization (or "smoothening") methods (see for instance Tikhonov and Arsenin [START_REF] Tikhonov | Solutions of Ill-posed Problems[END_REF] for a fundamental description of the underlying theory). Guyer et al. [START_REF] Guyer | Quantitative implementation of Preisach-Mayergoyz space to find static and dynamic elastic moduli in rock[END_REF] proposed three different ways to realize the regularization, resulting in three inversion algorithms:

"simulated annealing" (direct minimization of a regularized objective function), "normal modes" (approximation of the PM density with an orthonormal set of functions) and "exponential decay" (imposing an explicit form of the PM density characterized by a limited number of model parameters). Visually, the agreement between the three obtained density surfaces from the same set of data is acceptable;

however, an extensive quantitative comparison and error calculation has not been performed. Nonetheless, using numerical tests with simulated data, it is straightforward to calculate the actual PM density reconstruction error for any of these approaches directly.

In this paper, we look at the problem from a slightly different point. On one hand, the fact that the PM density inversion is underdefined can be compensated to a certain extent by operating a special stress protocol that consists of an extended sequence of loops with different upper and lower stresses. Even though the finite measurement precision limits the amount of experimental stress-strain data points, it can be significantly increased. On the other hand, we note that the problem of reconstructing the PM density is highly ill-posed: similar stress-strain responses may correspond to highly different PM distributions. Indeed, as shown by Mayergoyz [START_REF] Mayergoyz | Hysteresis models from the mathematical and control theory points of view[END_REF], the PM density can be obtained as the second derivative of the hysteresis information contained in the measured data. However, it is well known that each differentiation procedure is an extremely bad conditioned operation, especially on noisy data, and tends to blow up small errors on the measurements. Furthermore, once the PM density has been inferred from the data by double differentiation, it is then necessary to perform a 2D integral calculation of the PM density over a specific area to predict the stress-strain curves in response to an arbitrary stress protocol.

The above mentioned traditional practice opens a simple question: is it possible to avoid the double differentiation and the subsequent double integration, and characterize the material by a property that relates to the data more directly? In other words, can we use for instance an integral characteristic instead of a differential one?

Sure enough, there is not much scientific progress in using an integral property instead of a differential one. However, from a practical point, the approach that will be outlined in the next paragraphs completely avoids solving an ill-posed problem and performing a badly conditioned operation. Indeed, we will show that the characterization of both hysteretic and non-hysteretic material elasticity can be solved with a precision of a few percents even for very noisy data, if we are not obliged to infer the traditional representation of the PM density.

Ideally, when solving the hysteretic elasticity characterization problem, the inferred material property describing the mechanical hysteresis should not depend on the experimental conditions; be able to predict results for new experiments; and be obtained by means of a well-posed inversion.

Such a method is proposed in the present communication and is discussed in detail below for the treatment of simulated noisy data and real experimental data.

PM density and IPM density

The traditional treatment of hysteresis using the Preisach space (or PM space) consists in the evaluation of the following expression for the uniaxial stress(σ)-strain(ε) relationship:

( ) ( ) , NH H history of ε ε σ ε σ σ = + . (1)
Here, ε NH (σ) is a one-variable function representing the non-hysteretic contribution (sometimes referred to as classical nonlinearity since it is generally approximated by a simple polynomial expansion [START_REF] Landau | Theory of Elasticity[END_REF]). The second term accounts for the hysteresis itself and is history-dependent. The dependence on the stress history appears in the following way. We suppose that any representative (physically small) volume of the material contains a large number of hysteretic mechanical units, which can be found in one of two states: "open" or "closed". A single mechanical element is assumed to be characterized by two stress values σ o and σ c , so that it opens if σ < σ o when the actual stress σ decreases, closes if σ > σ c when the stress σ increases. In all other cases it keeps its previous strain state. Here we adopt the conventions (convenient for compression tests) that σ and ε are both positive for compression. Even though there is no experimental evidence for this, we can associate the hysteretic elements 13]. Experimental evidence for the connection between crack density and macroscopic effects of hysteresis can be found in Van Den Abeele et al. [START_REF] Van Den Abeele | Nonlinearity and crack density in thermally loaded CFRP[END_REF].

We assume that the material (or better, each representative volume) is characterized by a distribution ρ(σ c ,σ o ) of hysteretic elements (called PM density) with characteristic stress values σ o and σ c , so that ρ(σ c ,σ o )dσ c dσ o is the probability to find an element in the rectangle (σ c ,σ

c +dσ c )×(σ o ,σ o +dσ o ) and
( )

, 1 c o c o d d ρ σ σ σ σ +∞ +∞ -∞ -∞ = ∫ ∫
. We also accept that open units do not contribute to the stress-strain relationship given in Eq. ( 1) and that the total hysteretic strain contribution is given by γ when all elements are closed. Under these assumptions, the hysteretic addition ε H to the total strain simply equals the fraction of closed elements multiplied by γ:

( ) ( ) , , c c H c o c o c o c o d d d d ε γ ρ σ σ σ σ ρ σ σ σ σ Ω Ω = = ∫∫ ∫∫ % , (2)
where Ω c denotes the area in the (σ c ,σ o )-or PM-space containing the closed elements, and ρ γρ = %

. An example is displayed in Figure 1

The area of integration is a right-angled triangle defined by the corners (σ c ,σ c ), (σ o ,σ o ) and (σ c ,σ o ) (as for instance indicated in Figure 4 which will be discussed later). It is easy to prove that Eq. ( 2) can be rewritten as ( )

1 2 3 1, ... 1 i H n i i n H H H H H ε = = - + -± =- - ∑ , (4)
where the values H i correspond to the IPM densities at the (σ c ,σ o ) points i=1,2..n on the staircase-line as illustrated in Figure 1. The last sign is "+" if n is odd and "-" otherwise.

Here, we assume for simplicity that the section 1-2 is vertical (in practice, this is easily achievable by introducing an additional point very closely to the first one). Otherwise, Eq. ( 4) should read:

2 3 4 ... H H H H ε = - + -.
The relation between the PM density and the IPM density can also be rewritten as follows:

( ) ( ) 2 , , c o c o c o H σ σ ρ σ σ σ σ ∂ = -∂ ∂ % . ( 5
)
As it will be shown later, the advantage of the IPM density is that its distribution can be found directly from the stress-strain data. In addition, we avoid the ill-posed operation of calculating the second derivative of potentially noisy data. As a result, the IPM density will appear to be a more relevant material characteristic than the conventional differential one.

Data treatment procedure for IPM density reconstruction

In this Section we formulate the essential algorithm for retrieving the IPM density from obtained stress-strain data. The proposed method requires a specific stress protocol as displayed in Figure 2. This protocol contains N L =N(N+1)/2 loading-unloading cycles,

with N an integer number (e.g. N=8 in Figure 2), and is constructed in such a way that each compression-decompression cycle in the left hand side of Figure 2 corresponds to one of the possible triangles in the PM space, formed by any two orthogonal grid lines in the figure on the right hand side. Inversely, each of the N L =N(N+1)/2 triangles that can be found on the N×N grid, corresponds to a compression-decompression cycle in the stress protocol (a single loop in the stress-strain relationship). In the rest of this paper, we will refer to this particular protocol as the N-level stress protocol.

Starting from the basic equations of the PM model, i.e. Eqs. (1-2), it is easy to verify that the thickness ∆ε of a hysteretic loop at a constant stress value σ equals the integral of the PM density over a rectangle with corners (σ,σ), (σ max ,σ), (σ max ,σ min ) and

(σ,σ min ), where σ min and σ max denote the lowest and highest stresses for this particular loop (as shown in Figure 3):

( ) ( ) ( ) max min , c c o o d d σ σ σ σ ε ε σ ε σ σ ρ σ σ σ ↑ ↓ ∆ = - = ∫ ∫ % . ( 6
2 max max , ε ρ σ σ σ σ ∂ ∆ = ∂ ∂ % . ( 7
)
The conventional practice of inferring the PM density thus implies the second derivative of the loop thickness information in the experimental stress-strain curves. Since experimental stress-strain curves are usually subject to substantial noise, this differential operation is badly conditioned. We can avoid the differential density reconstruction by working with the integral formulation.

In the integral approach, it suffices to measure the thicknesses ∆ε from the available stress-strain data for all σ-levels at the grid points in the PM space. For an N×N grid with N L =N(N+1)/2 possible triangles (i.e. loops in the stress-strain relation)

and the corresponding N-level stress protocol, the maximum number of possible "nonzero" thickness measurements is N M =N(N 2 -1)/6. This means that each N-level protocol gives us the PM density integrals over N M rectangles, overlapping as well as nonoverlapping. Once these values are known, it is possible to infer the integrals over the individual 1×1 squares in the PM space (which we call blocks or bins). The number of such bins is N B =N(N-1)/2. Remark that the number of unknown bin integrals N B is always smaller then the number N M of known data points (i.e. loop thickness measurements, or rectangles (bins combinations)), except in the trivial case for N=2, when both numbers equal 1. As a consequence, we have an over-determined system of linear equations that can be solved by minimizing the corresponding objective function: ,

M N i k l i k l i F I ε γ = ∈ ⎛ ⎞ = ∆ - ⎜ ⎟ ⎝ ⎠ ∑ ∑ , ( 8
)
where ∆ε i is the result of the i-th thickness measurement and , ,

k l k l i I ∈ ∑ represents the
summation of all bin integrals (i.e., PM density integrals over the (k,l) th bin) belonging to the corresponding i-th rectangle. The residual value of the objective function can be used as a criterion for deciding whether we are dealing with a good PM system or a system that significantly deviates from the ideal one. (Note that an ideal PM system is defined by Eqs. [START_REF] Preisach | Über die magnetische Nachwirkung[END_REF][START_REF] Krasnosel | Systems with Hysteresis[END_REF] or [START_REF] Preisach | Über die magnetische Nachwirkung[END_REF][START_REF] Ortin | Preisach modeling of hysteresis for a pseudoelastic Cu-Zn-Al single crystal[END_REF], and its properties can by found in the classical article by Mayergoyz [START_REF] Mayergoyz | Hysteresis models from the mathematical and control theory points of view[END_REF] or in the fundamental book by Krasnosel'skii and Pokrovskii [START_REF] Krasnosel | Systems with Hysteresis[END_REF].

The problem of calculating the individual bins values from their combinations basically requires the subtraction of two closely valued positive numbers, each known with some precision. In that case, the relative error of the result is much higher than the relative errors of the individual components. In other words, the calculation of the individual bins values is a coarse-grid equivalent of the continuous differentiation given by Eq. ( 7). However, in our current approach, we don't really need the bin values themselves, but only a specific combination of them corresponding to the grid points where the IPM density is calculated. To find the IPM density in a given grid point (σ c ,σ o ) we must collect all bins located at the upper-left hand side from that point (see Figure 4). By summing bin integrals, we avoid the above mentioned instability problem since the sum of positive numbers always has a small relative error. In other words, we calculate a combination of items from another combination of items, that's why the requirements to error on each item are not too high. Note that, despite the fact that we combine bins in the IPM reconstruction algorithm, it is advisable to assure the stability of the process of deriving the bins themselves. A typical regularization technique [START_REF] Tikhonov | Solutions of Ill-posed Problems[END_REF] consists in considering a stabilizing addition to the objective function (Eq. 8) of the following kind

( ) ( ) 2 2 , 1 , , , 1 , ,
'

k l k l k l k l k l k l F F I I I I α α + + = + - + - ∑ ∑ , (9)
which penalizes large differences between neighboring bins and smoothens the coarsegrained PM density. The same consideration was made by Guyer et al. [START_REF] Guyer | Quantitative implementation of Preisach-Mayergoyz space to find static and dynamic elastic moduli in rock[END_REF], however, with one major substantial difference. In this article the proposed smoothening is of vital importance for finding the differential density, while here it represents just a small correction that is not critical at all. The stabilizing parameter α must be predetermined depending on the noise value (as will be discussed later).

As presented in Figure 4, the IPM density H(σ c ,σ o ) consists of elementary square bins and small triangular bins which are located at the diagonal σ o =σ c . To estimate the integrals in the triangular bins, we assume that the PM density in each triangular bin has a constant value (N in total). This is, of course, not exactly true and may lead to errors which will be discussed later. From Eq. ( 6), it follows that the area inside any closed stress-strain loop, extending from σ min to σ max under the assumption of

a constant PM density 0 ρ % , is ( ) max min 3 0 m a x m i n 6 A d σ σ ε σ ρ σ σ = ∆ = - ∫ % . ( 10
)
A c c e p t e d m a n u s c r i p t 13 13 Thus, by measuring the areas A for all the smallest stress-strain loops, one easily obtains the values of the triangular bin integrals:

( ) ( ) 2 max min 0 max min 3 2 A σ σ ρ σ σ - = - % (11)
Here 0 ρ % , σ max and σ min are, of course, different for different triangular bins.

With this, we have all information for obtaining the IPM density H(σ c ,σ o ) at the discrete grid points. Finally, we interpolate the set of IPM points with a smooth function to get the IPM density distribution everywhere in the (σ c ,σ o ) space. Inversely, this distribution can then again be used to predict the hysteretic contribution to the total strain following any arbitrary stress protocol.

We now return to the nonhysteretic contribution to the total strain. The component ε NH (σ) in Eq. ( 1) can only be reconstructed once the IPM density is found.

To do this, we consider the loading curve ( )

ε σ ↑
for the largest loop of the protocol and we note that

( ) ( ) ( ) min , NH H ε σ ε σ σ σ ↑ = + . ( 12
)
Alternatively, we can also consider the unloading curve for which

( ) ( ) ( ) ( ) max min max , , NH H H ε σ ε σ σ σ σ σ ↓ = + - . ( 13
)
Using both the measured loading and unloading branches and with the known IPM density values at hand, we can retrieve ε NH (σ) twice. Certainly, for an ideal PM system The performed numerical analysis enabled us to choose the number N from the following compromise consideration: on one hand, to coarse grid hides local features of the density, but, on the other hand, too fine grid makes it difficult to measure thickness of the smallest internal loops. N=8 was established as a optimum in our situation.

In conclusion, the algorithm for the IPM density computation and nonlinear elasticity calculation can be summarized in the following seven steps:

1. Choose N and execute the corresponding N-level stress protocol depicted in Figure 2.

2. Measure the thickness of all loops at the grid points.

3. Minimize the objective function (Eq. 9) to get all square bins.

4. Measure the areas of the smallest loops to get all triangular bins.

5. Collect all bins related to the IPM density at the grid points.

6. Interpolate the discrete set of IPM points to obtain the IPM density distribution everywhere in the (σ c ,σ o ) space.

7. Determine the non-hysteretic contribution.

In the next Section, we will illustrate the use of this method for simulated data, and discuss how to choose the optimal values for certain parameters used in the algorithm.

Testing the algorithm with synthetic data

In order to test the proposed method and illustrate the IPM density reconstruction, including errors, we performed the following actions: 1. We assume a particular analytical model form for the PM (or IPM) density.

2. We compute the response (stress-strain loops) for an 8-level stress protocol by means of Eq. (4).

3. We accept a particular noise model and contaminate the stress-strain response with noise of given amplitudes.

4. We approximate the noisy stress-strain data by piecewise polynomials of a given degree.

5. We reconstruct the IPM density using the algorithm discussed in section 3, and calculate the global reconstruction error by comparing the result with the assumed analytical IPM density.

6. We vary certain parameters of the algorithm in order to minimize this error.

7. We repeat actions 3-6 for different noise amplitudes.

The above detailed plan makes our method and its results fully reproducible by the reader.

The particular model form we assumed for the verification of the reconstruction model is a two-dimensional Gaussian distribution of the PM density:

( ) ( ) ( ) 2 2 2 2 || , e x p e x p 4 4 c o c o c o m d d σ σ σ σ ρ σ σ ρ ⊥ ⎛ ⎞ ⎛ ⎞ - + = - - ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠ % % (14)
where the parameters ⊥ d and || d account for the widths of the distributions perpendicular and parallel to the diagonal σ o =σ c , and m ρ % is a normalization constant whose value is not critical for the discussion here. For this particular formulation, we numerically computed the IPM density H(σ c ,σ o ) given by Eq. 3, and obtained, using Eq. 4, the hysteretic component of the stress-strain curve for an 8-level stress protocol In order to simulate real experimental conditions, we introduced noise contamination on the numerical stress-strain data. In our simulation, we used a model of purely additive noise:

ηξ ε ε + = ' , ( 16
)
where ε is the numerically calculated value of the strain response, ' ε is the strain value contaminated with noise, η is the noise amplitude (an instrumental constant), and ξ is a random number with Gaussian distribution of zero mean and unit variance. The choice of the noise model resulted from a comparison of experimental and simulated stressstrain data. Using the additive noise model, the behavior of the polynomial approximation error versus loop number for real data and for simulated data subjected to Eq. ( 16) were found to be quite similar (we omit the graphs for brevity). In an additive noise model, smaller loops are more drastically affected by the noise than larger loops. Figure 5 illustrates a typical case for η=3 10 -5 . The loop corresponding to a small excursion with stress amplitude σ max -σ min equal to 12.5% of the total stress excursion (50 MPa) almost completely loses its shape, while a big loop with σ max -σ min =50 MPa is still clearly discernible and reproducible. As a result, it is reasonable to use a different (smaller) degree of polynomials to approximate the smallest loops.

In order to reconstruct the IPM density, we first approximated the noisy stressstrain data by piecewise polynomial functions of degree N p1 for the smallest loops and N p2 for the other loops. Next, we determined the thicknesses of all loops except for the smallest ones. For the smallest ones, we measured their areas. If due to the presence of noise, some lengths or areas are found to be negative, we assign a zero value.

The obtained values are then transformed into the IPM density using the aforementioned procedure. The reconstruction algorithm uses several parameters: the degrees of the polynomials N p1 and N p2 for the approximation of the stress-strain loops (the smallest and the others, respectively), the regularization parameter α (see Eq. ( 9)) for the bins, and the spline order N s for the final 2D IPM approximation of the gridded data. Their With respect to the reversible elasticity parameters, as listed in Eq.(15), we obtained an error on the static Young's modulus E o of 0.9% for the largest noise level, and an error of 4% for the first order nonlinear coefficient β.

Numerous numerical experiments lead us to the following general conclusions about the optimal parameter values.

A c c e p t e d m a n u s c r i p t

• The choice of the spline order N s does not influence the bin and gridded IPM errors. For N s =3 the additional error introduced by the splines is minimal, at all noise amplitudes.

• For polynomial degrees N p1 and N p2 we recommend the values N p1 =3 and N p2 =6, again for all noise amplitudes.

• The optimum regularization parameter α depends on the noise amplitude η (in accordance with the general theory by Tikhonov and Arsenin [START_REF] Tikhonov | Solutions of Ill-posed Problems[END_REF]). Nevertheless, the global error is not crucially sensitive to α, since α serves as an internal regularization only, which mainly affects the bin values.

Table 1 lists the recommended values of α for different η. The use of an α value which deviates by 50 or even 100% from the recommended value is often not critical.

For instance, for η = 4•10 -5 the value α=1.5 (optimal) leads to a global reconstruction error of 4.0% while at α=0.5 the error equals 4.2%. However, if one takes α=10 then the global reconstruction error grows to 7.1%. So, we may conclude that only orders of α-values are essential.

A proper choice of the degree N p1 is more important: for example, for N p1 =4 instead of 3 (optimal), the global error easily increases by 0.3-0.8% for noisy data (η > 1•10 the assumption of a constant PM density inside each of the triangular bins then is violated more easily.

IPM density reconstruction from experimental data

We now apply the reconstruction method to real experimental data resulting from uniaxial compression tests on natural building stones. We have analyzed the following rock materials: Serena sandstone, Sander sandstone, Vicenza limestone, Lecce limestone and Portland limestone. Some physical properties of these rocks are listed in Table 2. Among other natural stones, these materials were extensively tested in laboratory conditions within the EU-FP5 grant DIAS [15] for comparison with the results of a combination of three prototype devices (micro-drilling, micro-indentation and Rayleigh wave acoustics) developed for in-situ characterization of their mechanical and acoustical parameters in historical buildings.

After having determined the uniaxial compression strength (UCS) for each rock type, we instrumented cylindrical samples (with a length of about 50mm, and a diameter of 25mm) with 2 dynamic strain gauge extensometers (Instron 2620-603 with gauge length 10 mm (Ext1) and 25 mm (Ext2) and full-scale range ± 1mm) and placed them between the grips of a MTS 810 TestStar load frame (Figure 8). This system is capable of loading up to 100kN. The MTS software allows to translate custom designed command signals, such as our 8-level stress protocol. The experiments are load-driven and can be controlled very accurately (error between command and actual loading is less than 0.025kN). The tests on the different samples are performed with the same speed (0.125kN/s). For each rock type, we define the loading maximum σ 60% for the 8- In the first part, we initialize the complete system by performing ten cycles of 60% of the maximum stress loading and unloading. The initialization assures that the stress-strain behavior becomes repeatable. In the second part, the 8-level stress protocol is executed. Finally, we continue with additional loading-unloading commands which are used to test the IPM reconstruction. The first "conditioning" loops and the final "additional" loops are not used in the reconstruction algorithm. Only the second part between 5000 and 13000 seconds is used for the inversion.

Figure 10, similarly to Figure 5 for the synthetic data, illustrates two loops of the stress-strain curve for Serena sandstone (σ 60% =58MPa). The large loop (0→σ 60% →0) switches all the hysteretic elements contained in the covered pressure range. This loop evidently shows the largest difference between loading and unloading curves. The inset represents one of the smallest loops obtained at the highest ambient pressure (58MPa→51MPa→58MPa). At these pressures almost all hysteretic movement (in compression or shear) of crack surfaces inside the rock is disabled [START_REF] Mavko | Frictional attenuation: An inherent amplitude dependence[END_REF], and the loop becomes considerably less hysteretic than the others. Correspondingly, the two loops in coincidence between the inverted data points (large and small crosses) confirms the validity of the PM model for this material.

In Figure 12 we have plotted two IPM portraits, one for a sandstone (Sander)

and one for a bioclastic limestone (Lecce). The petrographic parameters (type and size of grains, pores and binding material) of these two rocks are considerably different, as well as their physical properties. In particular, as can be seen from Table 2, the porosity, grain size, and UCS values all differ about an order of magnitude. The reconstructed hysteretic portraits of both rocks deduced from their corresponding experimental stressstrain curves are found to be considerably different too: the IPM densities have distinctive shapes. We observe that the contour lines for the IPM density of Sander sandstone are more inclined with respect to the diagonal than for Lecce limestone. This means that Serena sandstone has relatively more hysteretic features with low opening pressures.

Figure 13 summarizes all reconstructed IPM portraits for the rocks considered in the present study. We clearly see that the two sandstones display a similar structure of the IPM distribution with rather different characteristics compared to the three limestones. Some characteristic parameters of the IPM-distributions are listed in Table 3. The maximum IPM density H max corresponds to the hysteretic contribution to the strain ε Η,max at 60% of the ultimate compression strength (σ 60% ). We notice that these values are comparable for rocks of the same types: around 8•10 -4 for the class of the 'harder' sandstones, and from 0.6 to 2•10 -4 for the class of the 'softer' limestones.

Interestingly enough, the normalized parameter H max /σ 60% 2 which characterizes the total "amount" of hysteresis per unit of pressure-pressure space does not show the same classification tendency. This means that apparently similar portraits for rocks of the same class, despite the visual similarity, can have significantly varying quantitative parameters. We have also found that, if properly prepared, two samples of the same rock produce the same IPM densities within their error intervals.

The 2D information in the IPM density distribution makes it hard to directly compare quantities between the rocks. On the other hand, simple numerical parameters such as the ones listed in Table 3 carry severely limited information on the full reconstructed portraits. For a quick and meaningful understanding of the differences between the rocks, we compared particular 1D profiles expressed in normalized variables in Figure 14. We therefore considered the two sections H(σ c ,σ o =0) and

H(σ c =σ 60% ,σ o ).
Note that, despite the reduction of the number of arguments, both characteristics are not "local" but rather integral quantities, since they cover all hysteretic elements contained in particular triangles with only one varying side. In fact the two profiles simply correspond to the hysteretic strain contributions on the loading and unloading branches of the largest loop. In order to compare the comportment for different rocks, we first normalized these dependencies, H(σ,σ o =0)/H max and H(σ c =σ 60% ,σ)/H max , and plotted them as function of σ/σ 60% in Figure 14 for all the rocks under test. Certainly, the values σ 60% and H max used in the definition of the section and in the normalization bring in a certain influence of the arbitrary parameter of 60% UCS. Nevertheless, in a similar manner as in the IPM portraits, rocks of the same types tend to group around similar profiles. In conclusion, we can say that the normalized 1D curves, although suffering from the point of generality, offer a simple means for the comparison of the hysteretic properties.

In an analogous manner, we analyzed the non-hysteretic contribution ( ) NH ε σ to the total strain for the different rocks using Eq. ( 12) and/or [START_REF] Aleshin | Friction in unconforming grain contacts as a mechanism for tensorial stress-strain hysteresis[END_REF]. When normalized to its maximum at σ 60% and expressed in terms of the normalized stress value σ c /σ 60% , we can observe the same clustering of the two groups, as illustrated in Figure 15. In addition,

A c c e p t e d m a n u s c r i p t

24 24 we used the first order nonlinear expression of the non-hysteretic stress-strain relation (Eq. 15) to quantify the values of the quasi-static Young's Modulus E 0 and the first order nonlinearity coefficient β. The values are listed in Table 4 and, for the Young's modulus, they are compared to the values obtained from dynamic velocity measurements which were performed on similar cylindrical samples using resonance spectroscopy. The observation is that the static modulus is generally smaller than the dynamic modulus and that the limestones have a smaller 'classical' first order nonlinearity coefficient compared to the sandstones.

As a final point in the experimental analysis, we also investigated the predictive power of the considered methodology. As mentioned before, the final part of the executed protocol (Figure 9, gray line for times later than 13000 seconds) was not used for the construction of the system's IPM portrait. It was added to the protocol in order to serve as a test for the stress-strain prediction. Figure 16 illustrates the actual measurement of the stress-strain curve together with the predicted behavior that was calculated using Eq. (1), Eq. ( 4) and the reconstructed IPM density for Sander sandstone (Figure 13). The agreement is most satisfactory. The experimental and predicted shapes only weakly deviate from each other: the relative error of the deviation

p i i e x p i RMS RMS ε ε ε - ) is 2 %. ( , e x p , , ( )/ ( )
The major difference is a constant shift, which can be attributed to the fact that successive cycles are never perfectly repeatable. Even after a long initialization process, a small shift of the successive loops remains present in the experiment.

A c c e p t e d m a n u s c r i p t 6. Conclusions

We introduced the IPM density as a means for characterizing the hysteretic mechanical properties of solid materials and showed that it is a robust concept, even in the presence of noise. The reconstruction of the IPM density is based on a particular force protocol and involves simple mathematical data treatment. The methodology has been illustrated

for synthetic data and applied for the characterization of several classes of natural building stones. The reconstruction of the IPM density and the analysis of the nonhysteretic stress-strain contribution revealed that sandstones are more hysteretic and more nonlinear than limestones in the compressive stress range of 0 to 60% of the limiting strength. The physical reason for this difference is to be found in the different composition and binding properties of the two classes.

Even though a limiting stress value of 60% of the UCS was used, the protocol and methodology can be executed for any arbitrary stress range. The resulting IPM density then enables us to predict any stress-strain relation within this range with great accuracy.

In the future, we intend to characterize the influence of damage mechanisms on the IPM density of various materials by investigating the alterations in the distribution characteristics as the result of progressive thermal and/or mechanical fatigue loading.

This may for instance lead to a very interesting application in the quantitative evaluation of the effect of fire damage on concrete walls, housing, tunnels, etc.

In addition we plan to perform an intense comparative study between the static nonlinear and hysteretic elasticity, and the dynamic nonlinear properties of the same material. This will inform us about the link, or the absence of a link, between the two strain ranges. Figure 11. Non-hysteretic contribution ε NH (σ) (discrete data points) and the largest hysteretic loop (solid lines) in the stress-strain behavior of Serena Sandstone. The small crosses correspond to the inversion of ε NH (σ) from the loading curve (Eq.12) and the large ones are obtained from the unloading curve (Eq. 13). 3 and defines the color scale of the plot. Lecce limestone. The values σ c and σ o range from 0 to σ 60% (60% of UCS) for each rock. The maximum IPM density H max corresponding to this interval is listed in Table 3 and defines the color scale of the plot.

A c c e p t e d m a n u s c r i p t 4 4

 4

A c c e p t e d m a n u s c r i p t 5 5

 5

A c c e p t e d m a n u s c r i p t 6 6

 6

A c c e p t e d m a n u s c r i p t 7 7

 7 hypothetically with internal cracks, which are closed at higher compression and open at smaller values. It is then obvious that σ o <σ c for any element. The conjectural association of hysteretic units with microcracks and partially contacting surfaces has been studied in recent work by Pecorari[START_REF] Pecorari | Nonlinear interaction of plane ultrasonic waves with an interface between rough surfaces in contact[END_REF][START_REF] Pecorari | Adhesion and nonlinear scattering by rough surfaces in contact: Beyond the phenomenology of the Preisach-Mayergoyz framework[END_REF], and Aleshin and Van Den Abeele[11-

  in which the population of closed units is represented by the dark gray area. Taking into account the above accepted rules for individual element switching, any stress loading history leads to a configuration in the PM space in which the closed (Ω c ) and open (Ω o ) areas are divided by a staircase line (Figure 1) defined by the loading history and, in particular, carrying information about the relevant stress extrema (e.g. Krasnosel'skii and Pokrovskii [2]. Starting from the classical representation of the PM-density, ρ(σ c ,σ o ), we can also define the integral PM density (IPM density) H(σ c ,σ o ), which is a dimensionless quantity:

A c c e p t e d m a n u s c r i p t 9 9

 9

)

  strain value at a compression σ when loading or unloading the sample within the particular hysteretic loop. From the above equation it readily follows that ( )

A c c e p t e d m a n u s c r i p t 12 12

 12

A c c e p t e d m a n u s c r i p t 14 14 these

 14 two values must coincide. For realistic cases, they deviate with a discrepancy characterizing the validity of the discussed PM model.

A c c e p t e d m a n u s c r i p t 16 16 (-

 16 see Figure2). For the non-hysteretic component we choose a simple first order nonlinear σ-ε NH dependence where E o denotes the static Young's modulus and β reflects a first order elastic nonlinearity parameter. In the simulation we have set 5000MPa , with σ expressed in MPa. This expression is based on a fit of real experimental data, but is not critical to the present study since the IPM reconstruction algorithm works with loop thicknesses, thereby eliminating the reversible contributions.

A c c e p t e d m a n u s c r i p t 17 17

 17

Figure 6 .

 6 illustrates the resulting reconstructed IPM density, H rec (σ c ,σ o ), smoothed by a 2D spline interpolation, corresponding to the case with a noise amplitude η=3•10 -5 , and the corresponding local reconstruction error In some regions close to the diagonal σ o =σ c , the noisy data may lead to very small negative IPM densities in the spline approximated distribution.The underlying reason for the negative values, which are obviously not physical, is purely of mathematical nature. Even though we assume an anti-symmetric extension of the IPM density into the half-plane σ o >σ c , i.e. H rec (σ c ,σ o ) = -H rec (σ o ,σ c ), small negative values may appear as approximation errors, since the anti-symmetry axis σ o =σ c is angled 45° with respect to the IPM space grid lines. In any case, we have found that these approximation errors increase the reconstruction error only slightly. An alternative assumption which sets H rec (σ c ,σ o )=0 for σ o >σ c does not greatly affect the situation.

A c c e p t e d m a n u s c r i p t 18 18

 18

  role has been investigated in terms of the optimization of the reconstruction. The main criterion in the parameter's optimization is the final global IPM reconstruction error for the inverted and spline approximated distribution in comparison with the exact IPM distribution, i.e. squares reconstruction error was computed on a 700×700 mesh in the (σ c ,σ o )-space, restricted to the half space σ c <σ o . (As the final form of the reconstructed IPM density is given by a spline, it is possible to calculate the analytical function at any (σ c ,σ o )-point.) In addition, we also calculated the least squares reconstruction errors for the individual bin integrals and for the gridded IPM data (again, in comparison with the exact values).The values shown in Figure7represent the three above mentioned errors for different noise amplitudes and for optimized algorithm parameters. One clearly observes that the global reconstruction error on the bin values is much higher compared to the gridded and spline interpolated IPM reconstruction errors. We expected this beforehand, since the reconstruction of the bin values is a coarse grid equivalent solution of the ill-posed problem of the PM density reconstruction. This graph confirms that the IPM inversion is well-posed.

A c c e p t e d m a n u s c r i p t 20 20

 20

  of the material's strength at failure. During the experiment, both extensometers and the actual loading force are registered every 0.2s.

Figure 9

 9 Figure 9 illustrates the complete stress protocol which is loaded on the

Figure 10 are

 10 Figure 10 are the most and the least hysteretic ones within the protocol.

A c c e p t e d m a n u s c r i p t 23 23

 23

Figure 1 .

 1 Figure captions

Figure 2 .

 2 Figure 2. Specific N-level stress protocol (with N=8) used in the reconstruction

Figure 3 .

 3 Figure 3. Illustration of the strain thickness measure ∆ε at constant stress value σ for a

Figure 4 .

 4 Figure 4. The IPM density at the marked grid point (σ c ,σ o ) is the sum of the PM density

Figure 5 .

 5 Figure 5. The largest loop simulation of the 8-level stress protocol (0→50MPa→0) and

A c c e p t e d m a n u s c r i p t 30 30 Figure 6 .

 306 Figure 6. The reconstructed IPM density H rec (σ c ,σ o ), left, and the local reconstruction error, right, for synthetic data contaminated with additive noise (η=3•10 -5 ) based on the

Figure 7 .

 7 Figure 7. Global reconstruction errors for different noise amplitudes η using the

Figure 8 .

 8 Figure 8. Experimental set-up for cyclic uniaxial compression tests on cylindrical rock

Figure 9 .

 9 Figure 9. Typical experimental stress command protocol, including initialization cycles

Figure 10 .

 10 Figure 10. Experimental data corresponding to the largest loop of the 8-level stress

Figure 12 .

 12 Figure 12. Reconstructed IPM density distributions for two different rocks (Serena

Figure 13 .

 13 Figure 13. Reconstructed IPM portraits for various rocks. In the right column: Serena

Figure 14 .

 14 Figure 14. Comparison of the normalized IPM densities at σ o =0 and at σ c =σ 60% for

Figure 15 .

 15 Figure 15. Comparison of the normalized non-hysteretic contribution for different

A c c e p t e d m a n u s c r i p t 32 32 Figure 16 .

 3216 Figure 16. Prediction of additional loading-unloading cycles made for Sander

A c c e p t e d m a n u s c r i p t 33 33 Table 1 :

 331 Recommended values of the regularization parameter α for different noise amplitudes.

Figure 1 .

 1 Figure 1. Schematic representation of areas of closed (dark gray, Ω c ) and open (light gray, Ω o ) hysteretic units in a PM space (σ c ,σ o ), separated by a staircase line (1, 2, 3,

Figure 2 .

 2 Figure 2. Specific N-level stress protocol (with N=8) used in the reconstruction

Figure 3 .Figure 4 .

 34 Figure 3. Illustration of the strain thickness measure ∆ε at constant stress value σ for a

A c c e p t e d m a n u s c r i p t 41 41 Figure 5 .

 415 Figure 5. The largest loop simulation of the 8-level stress protocol (0→50MPa→0) and

Figure 6 .

 6 Figure 6. The reconstructed IPM density H rec (σ c ,σ o ), left, and the local reconstruction error, right, for synthetic data contaminated with additive noise (η=3•10 -5 ) based on the

Figure 7 .A c c e p t e d m a n u s c r i p t 44 44 Figure 8 .Figure 9 .

 74489 Figure 7. Global reconstruction errors for different noise amplitudes η using the

Figure 13 .

 13 Figure 13. Reconstructed IPM portraits for various rocks. In the right column: Serena

Figure 14 .σFigure 15 .Figure 16 .

 141516 Figure 14. Comparison of the normalized IPM densities at σ o =0 and at σ c =σ 60% for

  -5 ).

	Strictly speaking, the magnitude of the reconstruction error may vary for
	different values of the PM density widths ⊥ d and || d (see Eq. 14), but our numerical
	tests have shown that the general recommendations hold. Here we used ⊥ d =30 MPa and

|| d =50 MPa, which correspond to the differential PM density shown in the 2D grayscale plot of Figures 2 to 4. Smaller value of the widths will generally increase the error since

Table captions Table 1 .

 captions1 Recommended values of the regularization parameter α for different noise amplitudes.

Table 2 .

 2 Material properties

Table 3 .

 3 Maximum IPM density H max , maximum pressure σ 60% of the protocol, and relative H max related to the total pressure area σ 60% 2 for different rocks.

Table 4 .

 4 Linear and classical nonlinear moduli for different rocks. Comparison of the

Quasi-Static (QS) Young's modulus (deduced from a first order polynomial expansion of the nonhysteretic contribution in the quasi-static stress strain relation) and the dynamic (Dyn) Young's modulus (deduced from velocity measurements at zero stress).

Table 2 :

 2 Material properties.

		Density (kg/m³)	Total Porosity (%)	Mean Grain Size (µm)	UCS (MPa)
	Serena sandstone	2560	6.5	200-800	95.8
	Sander sandstone	2200	19.7	150-200	58.3
	Portland limestone	2030	20.1	300	31.4
	Vicenza limestone	1970	29.3	400-2000	29.3
	Lecce limestone	1460	47.4	80-100	11.0

Table 3 :

 3 Maximum IPM density H max , maximum pressure σ 60% of the protocol, and relative H max related to the total pressure area σ 60% 2 for different rocks.

		H max (-)	σ 60% , MPa	H max /σ 60% 2 , MPa -2
	Serena sandstone	8.48•10 -4	58.0	2.56•10 -7
	Sander sandstone	7.55•10 -4	34.8	6.34•10 -7
	Portland limestone	1.86•10 -4	18.7	5.40•10 -7
	Vicenza limestone	5.92•10 -5	17.6	1.97•10 -7
	Lecce limestone	1.06•10 -4	6.50	2.65•10 -6

Table 4 .

 4 Linear and classical nonlinear moduli for different rocks. Comparison of the Quasi-Static (QS) Young's modulus (deduced from a first order polynomial expansion of the nonhysteretic contribution in the quasi-static stress strain relation) and the dynamic (Dyn) Young's modulus (deduced from velocity measurements at zero stress).

		E 0 -Dyn	E 0 -QS	1 st order QS-
		(GPa)	(GPa)	Nonlinearity β (-)
	Serena sandstone	25.0	20.5	414
	Sander sandstone	15.7	14.0	1070
	Portland limestone	20.7	18.1	63
	Vicenza limestone	20.0	18.4	76
	Lecce limestone	8.6	9.1	125

Acknowledgments

The authors gratefully acknowledge the support of the Flemish Fund for Scientific Research. (G.0206.02, G.0257.02, and G.0554.06), the provisions of the European Science Foundation Program NATEMIS and the European FP5 and FP6 Grants DIAS (EVK4-CT-2002-00080) and AERONEWS (AST3-CT-2003-502927).


U Polimeno

M Meo

Understanding the Effect of Boundary Conditions on Damage Identification process when using nonlinear elastic wave spectroscopy methods

Understanding the effect of boundary conditions on damage identification process when using nonlinear elastic wave spectroscopy methods U. Polimeno, M.

A c c e p t e d m a n u s c r i p t

Introduction

Due to economic reasons, there is a need to extend the life of commercial and military aircraft well past their original design life. However, this gives rise to the urgent need for the development of cost effective NDT techniques to monitor the condition of aircraft structures. The ultimate goal to increase aircraft safety and reduce accidents is to incorporate sensors, actuators and electronics in a structural health monitoring system to monitor and report on the health of the aircraft structure, and to actuate in order to prevent damage or failure of the aircraft structures. This will allow the aircraft operator to move away from costly time-based maintenance procedures toward real-time health condition monitoring allowing timely decisions on preventative and scheduled maintenance before failure occurs. Currently, a new class of promising NDT techniques, called Nonlinear Elastic Wave Spectroscopy (NEWS), which could be incorporated in a structural health monitoring system, has being developed and it monitors the integrity of structures by analysing the material nonlinear elastic behaviour [START_REF] Guyer | Nonlinear mesoscopic elasticity: evidence for a new class of materials[END_REF][START_REF] Ostrovsky | Dynamic nonlinear elasticity in geomaterials[END_REF][START_REF] Van Den Abeele | Nonlinear elastic wave spectroscopy (NEWS) techniques to discern material damage, part I: nonlinear wave modulation spectroscopy[END_REF][START_REF] Van Den Abeele | Nonlinear elastic wave spectroscopy (NEWS) techniques to discern material damage, part II: Single-Mode Nonlinear Resonance Acoustic Spectroscopy[END_REF][START_REF] Van Den Abeele | Micro-damage diagnostic using nonlinear elastic wave spectroscopy (NEWS)[END_REF][START_REF] Meo | Nonlinear elastic wave spectroscopy identification of impact damage on sandwich plate[END_REF][START_REF] Landau | Theory of Elasticity[END_REF][START_REF] Mf | Fundamental and applications of nonlinear acoustic of nonlinear wave propagation in mechanics[END_REF][START_REF] Van Den Abeele | Elastic pulsed wave propagation in media with second -or higher-order nonlinearity. Part I: theoretical framework[END_REF][START_REF] Meo | Barely visible Impact damage detection in aircraft composites[END_REF] caused by the presence of damage.

The NEWS techniques originally were developed to estimate the nonlinear elastic properties of heterogeneous geomaterial like rocks or sandstones. It has been shown that some damaged materials have a complex compliance and have a nonlinear behavior that cannot be explained with classical nonlinear models [START_REF] Guyer | Nonlinear mesoscopic elasticity: evidence for a new class of materials[END_REF][START_REF] Ostrovsky | Dynamic nonlinear elasticity in geomaterials[END_REF]. The stress strain curve of damaged materials contains terms that describe classical nonlinearity, as well as hysteresis, and discrete memory [START_REF] Ostrovsky | Dynamic nonlinear elasticity in geomaterials[END_REF]. Considering stress and strain

A c c e p t e d m a n u s c r i p t

for one dimensional problem, the constitutive relation for wave propagation analysis in quasi static condition can be written as follow:

( )

, K d σ ε ε ε = ∫ & ( 1
)
where σ is the stress, ε is the strain , ε& is the time derivative strain and K is the elastic modulus.

According to the experimental evidence [START_REF] Van Den Abeele | Nonlinear elastic wave spectroscopy (NEWS) techniques to discern material damage, part I: nonlinear wave modulation spectroscopy[END_REF][START_REF] Van Den Abeele | Nonlinear elastic wave spectroscopy (NEWS) techniques to discern material damage, part II: Single-Mode Nonlinear Resonance Acoustic Spectroscopy[END_REF][START_REF] Van Den Abeele | Micro-damage diagnostic using nonlinear elastic wave spectroscopy (NEWS)[END_REF][START_REF] Meo | Nonlinear elastic wave spectroscopy identification of impact damage on sandwich plate[END_REF] one of the possible formulations for the non linear non classical hysteretic elastic modulus is shown in equation [START_REF] Ostrovsky | Dynamic nonlinear elasticity in geomaterials[END_REF].

( ) ( ) ( ) [ ] { } , K K t sign ε ε βε δε α ε ε ε = -- - ∆ + + 2 0 1 & & K ( 2 )
where K 0 is the linear modulus, , β and δ are classical nonlinear coefficients of Taylor expansion of stress versus strain [START_REF] Landau | Theory of Elasticity[END_REF][START_REF] Mf | Fundamental and applications of nonlinear acoustic of nonlinear wave propagation in mechanics[END_REF][START_REF] Van Den Abeele | Elastic pulsed wave propagation in media with second -or higher-order nonlinearity. Part I: theoretical framework[END_REF], α is a material hysteresis measure that takes into account also the history of stress through the terms ( ) ε& sign , ε ∆ is the local strain amplitude, represented in equation 3(expression for simple continuous sine excitation [START_REF] Van Den Abeele | Micro-damage diagnostic using nonlinear elastic wave spectroscopy (NEWS)[END_REF]) and equation 4:

2 min max ε ε ε - = ∆ (3) 1 ) ( 1 1 ) ( 1
- = ⇒ < = = ⇒ > = ε ε ε ε ε ε sign dt d sign dt d (4)
The nonlinearity induced by the damage and simulated by this model can be estimated by two groups of methods available for NEWS application: resonance and non-resonance. Resonance methods analyse the dependence of the resonance frequency on the strain amplitude while exciting the sample at relatively low amplitudes. By observing the relative frequency shift it is possible to have a measure of the internal degradation of the microstructural properties of the material. In Fig. 1 a schematic representation of the method is shown; according with experimental evidence [START_REF] Guyer | Nonlinear mesoscopic elasticity: evidence for a new class of materials[END_REF][START_REF] Ostrovsky | Dynamic nonlinear elasticity in geomaterials[END_REF][START_REF] Van Den Abeele | Nonlinear elastic wave spectroscopy (NEWS) techniques to discern material damage, part I: nonlinear wave modulation spectroscopy[END_REF][START_REF] Van Den Abeele | Nonlinear elastic wave spectroscopy (NEWS) techniques to discern material damage, part II: Single-Mode Nonlinear Resonance Acoustic Spectroscopy[END_REF][START_REF] Van Den Abeele | Micro-damage diagnostic using nonlinear elastic wave spectroscopy (NEWS)[END_REF] the sample under a sweep load centred on one of its natural mode show different behaviour of structural response depending on the status of sample: intact or damaged. These methods, shown in Fig. 1 [START_REF] Ostrovsky | Dynamic nonlinear elasticity in geomaterials[END_REF][START_REF] Van Den Abeele | Nonlinear elastic wave spectroscopy (NEWS) techniques to discern material damage, part II: Single-Mode Nonlinear Resonance Acoustic Spectroscopy[END_REF], are usually referred to as Nonlinear Resonance Ultrasound Spectroscopy (NRUS) [START_REF] Van Den Abeele | Nonlinear elastic wave spectroscopy (NEWS) techniques to discern material damage, part II: Single-Mode Nonlinear Resonance Acoustic Spectroscopy[END_REF]. The resonance methods can be rather labour-intensive and require many readings to be taken (frequency response curves for various excitation levels). Therefore, they may not be suitable for fast in-process measurements.

A c c e p t e d m a n u s c r i p t

Non-resonance methods are used to study nonlinear effects on elastic waves propagating through the medium. According to the previous experimental work [START_REF] Van Den Abeele | Nonlinear elastic wave spectroscopy (NEWS) techniques to discern material damage, part I: nonlinear wave modulation spectroscopy[END_REF][START_REF] Meo | Nonlinear elastic wave spectroscopy identification of impact damage on sandwich plate[END_REF] and the mathematical model in equations ( 1) and ( 2), intermodulation phenomena can be observed in the structural response of damaged material when an excitation as shown in Fig. 2 (a) is used. In particular, the excitation is made by the superposition of two sine waves with fundamental frequencies f 1 and f 2 respectively. The spectrum response of the sample is characterized by two effects: the presence of harmonics of both fundamental excitations frequency (nf ; n=1,3,5…; for example: 3f 1 , 5f 1 ...nf 1 and 3f 2 , 5f 2 ... .nf 2 ), and the generation of sidebands around f 2 ( f 2 ±nf 1 ; n=2,4,…; for example:

, f 2 -2f 1 , f 2 +2f 1 …..f 2 -nf 1 , f 2 +nf 1 )
as shown in Fig. 2. While these two methods have been successfully used to detect impact damage [START_REF] Meo | Nonlinear elastic wave spectroscopy identification of impact damage on sandwich plate[END_REF][START_REF] Meo | Barely visible Impact damage detection in aircraft composites[END_REF], the effect of geometrical discontinuities, boundary conditions etc…on the nonlinear response of structures have not been assessed.

In particular, in this paper the effect of the clamp acting on the sample was studied to analyse how boundary conditions may affect the damage detection measures and introduce hysteretical nonlinear effects.

Experimental setup

The test pieces investigated were carbon fibre epoxy composite plates with unknown material properties. The size of the plate size was 210x100x4mm 3 and it was inspected with a C-scan to assess the presence of damage or defects. During this preliminary inspection no damage was found and the specimens appeared in pristine condition. The experimental setup used is shown in Fig. 3.

The composite plate was fixed on one side and left free to move on the other end as showed in Fig. 3. Two different values of torque were applied on the bolts (measured with a torquemeter), T 1

Hz

Amplitude

f 1 2f 1 3f 1 f 2 -2f 1 f 2 -f 1 f 2 f 2 +f 1 f 2 +2f 1 Hz (c). Input load f 1 f 2
Nonlinear modulation effects Hz Amplitude Amplitude and T 2 , in order to analyse the effect of the clamp different on the material nonlinear elastic behaviour. The torque T 1 corresponds to a firmly fixed sample edge, and T 2 was 20% greater then T 1 . The plate excitation was obtained by using a low distortion speaker, operating in low and medium acoustic frequency range. The speaker was driven by a signal generator through the power amplifier, and the data was acquired by using an accelerometer and a Brüel & Kjaer Pulse system. The accelerometer (ENDEVCO 27A12) was located at 2 cm from the free short edge and at 5cm from the long edges, on the front face of the plate. Two different types of NEWS tests were carried out: NRUS and NWMS [START_REF] Guyer | Nonlinear mesoscopic elasticity: evidence for a new class of materials[END_REF][START_REF] Meo | Nonlinear elastic wave spectroscopy identification of impact damage on sandwich plate[END_REF]. The first resonance frequency (flexural mode) for the plates considered was around 200 Hz. The experiments were performed using different drive amplitudes in order to measure the resonance shift versus external amplitude. For the NWMS tests the input signal was the superposition of two sine wave, where the fundamental frequencies (f 1 and f 2 ) are the first and the third resonance modes of the sample (199<f 1 <201 Hz, 2230<f 2 <2260 Hz). Different external load amplitudes were tested in order to correlate the non linear phenomena with the external amplitudes.

f 2 f 1 (a). Input load

Results discussion

Two different torque values were applied , (T 2 >T 1 ) and the NRUS and NWMS techniques were used. The results of the NRUS test for the two clamps values are shown in Fig. 4. The results in Fig. 4 show that in both cases there is a significant shift in the resonance frequency with increasing load amplitude. However, since the same external load was used for both test, when using T 2 as the clamping torque, the value of acceleration amplitude is higher than when using T 1 . This behaviour is caused by an increase of the stiffness due to the higher torque used. By measuring the nonlinear hysteretic parameter α (equation 5), an indication of the level of hysteretical nonlinear effects introduced by the contact between the sample and the clamp could be obtained. In a previous work [START_REF] Meo | Barely visible Impact damage detection in aircraft composites[END_REF] the parameter α was used to estimate relative damage severity, while in this work, α was used to study the effect of the boundary conditions on the damage detection process. The nonlinear hysteretic parameter α was evaluated as follows:

0 0 α ε - ∆ = i f f f ( 5
)
where ∆ε is the average strain amplitude, f 0 is the natural frequency of intact sample or the lowest resonance mode measured, f i is the natural mode measured for each drive amplitude, represented by the maximum value of each curve shown in Fig. 4. In Fig. 5, the behaviour of parameter α as a function of the acceleration is showed. An increase of the clamp torque causes an increase of the nonlinear parameter, showing that measuring only this parameter may give false data for the damage detection process. Various other amounts of clamping torque were studied and similar conclusions were drawn. In Fig. 6 and Fig. 7 the spectra of the accelerations measured with the NWMS test are shown for both clamp conditions. Different amplitudes of external loads were used in order to correlate the amplitude of nonlinear features versus the structure linear response, represented by fundamental amplitude. The generation of nonclassical nonlinear hysteretical features is clearly observable as highlighted by the presence of harmonics and sidebands. Increasing Clamp Torque

f 1 2 f 1 3 f 1 Frequency [Hz] f 2 -2f 1 f 2 -f 1 f 2 f 2 +f 1 f 2 +2f 1 f 1 2 f 1 3 f 1 f 2 -2f 1 f 2 -f 1 f 2 f 2 +f 1 f 2 +2f 1 Frequency [Hz]
A c c e p t e d m a n u s c r i p t The results show that by varying the boundary conditions (torque value) a variation of the structural nonlinear effects was observed. In particular, a decrease of the second and third harmonic, and the first right and left sidebands was observed as shown in Fig. 8 and Fig. 9. Increasing Torque Increasing Torque

f 2 -2f 1 f 2 -f 1 f 2 f 2 +f 1 f 2 +2f 1 f 1 2 f 1 3 f 1 f 1 2 f 1 3 f 1 f 2 -2f 1 f 2 -f 1 f 2 f 2 +f 1 f 2 +2f 1 Frequency [Hz]
Frequency [Hz] In Fig. 10 sidebands amplitude are displayed against amplitude of f 2 . Tests were conducted by varying the amplitude of frequency f 2 and keeping constant the amplitude of frequency f 1. In this case, the effect of a rising of clamp torque is a decrease of the amplitudes of non linear properties. The results of these tests show some of the limitations of the NDT based on the monitoring the nonlinear elastic effects. Non-classical nonlinear effects can also be generated also by the presence of contact regions, i.e. boundary conditions. There is, therefore, the need to optimise the NEWS methods and to further study these effects with other geometrical configurations, for example structures in pure contact (kissing bonds, joints etc…), and to develop methods to discern the nonclassical nonlinear effects caused by geometrical configurations from that generated by the presence of impact damage.

Increasing Torque

Increasing Torque

Increasing Torque Increasing Torque

A c c e p t e d m a n u s c r i p t Conclusion

This paper investigates the generation of nonclassical nonlinear effects due to changes in boundary conditions and analyse the limitations of NEWS techniques in detecting damage in the presence of unknown boundary conditions. In this study variable boundary conditions were studied by changing the torque of clamp system used to hold the structure examined.

Two different NEWS techniques were used and both techniques, NRUS and NWMS, highlighted the presence of nonclassical non linear features in the spectrum responses. In the first case a considerable shift of first natural mode was clearly observable in the undamaged sample.

Using NWMS, the intermodulation effects were evident, highlighted by the presence of sidebands and harmonics, induced by the boundary conditions. However, the non linear properties showed a decrease in amplitude when the clamp torque was increased. The results of these tests demonstrate that boundary condition can cause linear material to behave in a non linear elastic manner, and therefore providing false-positive indications of damage.

The results shows that more work is needed to demonstrate the effectiveness of the NEWS methods and further research studies and methodology development are needed to discern nonlinear effects generated by contacts between mating parts from those generated by the presence of damage.

Fig. 1 .

 1 Fig. 1. Resonance Frequency versus drive amplitude. (a) intact sample (b) damaged sample.

Fig. 2 .

 2 Fig. 2 . Linear and non linear modulation effect. (a) Excitation frequencies. (b) of a linear undamaged. (c) Excitation frequencies. (d) Nonlinear response of a damaged material : harmonics and sidebands.

Fig. 3

 3 Fig. 3 Schematic representation of equipment configuration.

Fig. 4 .

 4 Fig. 4. NRUS tests on one sample with two clamping torque: T 1 and T 2 , where T 2 >T 1 .

Fig. 5 .

 5 Fig. 5. Non linear parameter versus acceleration. The marker "■" indicates data for torque level T 1 and the marker "▲" indicates data for torque level T 2 , where T 1 <T 2 .

Fig. 6 .

 6 Fig. 6. Fourier spectra of accelerations measured with NWMS approach relative to the clamp condition T 1 . (a) Increasing f 1 amplitude. (b) Increasing f 2 amplitude.

Fig. 7 .

 7 Fig. 7. Fourier spectra of accelerations measured with NWMS approach relative to the clamp condition T 2 . (a) Increasing f 1 amplitude. (b) Increasing f 2 amplitude.

Fig. 8 .

 8 Fig. 8. Effect of the clamping torque on the observed second and third harmonics. The marker "■" indicates data for torque level T 1 and the marker "▲" indicates data for torque level T 2 , where T 1 <T 2 .

Fig. 9 .

 9 Fig. 9. Effect of the clamping torque on the observed sidebands against F1 amplitude. The marker "■" indicates data for torque level T 1 and the marker "▲" indicates data for torque level T 2 , where T 1 <T 2 .

Fig. 10 .

 10 Fig. 10. Effect of the clamping torque on the observed sidebands against F2 amplitude. The marker "■" indicates data for torque level T 1 and the marker "▲" indicates data for torque level T 2 , where T 1 <T 2 .

m.meo@bath.ac.uk

& & & &


Corina Fetecau
email: fetecau@math.tuiasi.ro

D Vieru

A note on the second problem of Stokes for Newtonian fluids

Keywords: Stokes' second problem, Newtonian fluids, Exact solutions

 3 and 4). The amplitudes of the wall shear stresses corresponding to the cosine and sine oscillations are almost identical, except for a small initial time interval. The time required to attain the steady-state for the cosine oscillations of the boundary is smaller than that for the sine oscillations of the boundary. This time decreases if the frequency of the velocity of the boundary increases.

Introduction

The study on the flow of a viscous fluid over an oscillating plate is not only of fundamental theoretical interest but it also occurs in many applied problems. In the literature this motion is termed as Stokes' second problem or Rayleigh problem [START_REF] Schlichting | Boundary Layer Theory[END_REF]. The first exact solutions corresponding to this problem for non-Newtonian fluids seem to be those of Rajagopal [START_REF]A note on unsteady unidirectional flows of a non-Newtonian fluid[END_REF]. Recently, Erdogan [START_REF] Erdogan | A note on an unsteady flow of a viscous fluid due to an oscillating plane wall[END_REF] established exact transient solutions for the motion of a Newtonian fluid due to the cosine and sine oscillations of the plate. For large times these solutions tend to the steady-state solutions. However, the Erdogan's solutions are complicated enough and they are not presented as sum of the steady-state and transient solutions. The purpose of this note is to present new exact solutions for the Stokes' second problem describing the flow for small and large times. In order to obtain these solutions the Laplace transform method is again used, but our solutions are directly presented as sum of the steady-state and transient solutions and are simpler than the previous ones.

Direct computations show that the steady-state solutions can be reduced to the classical forms [START_REF]A note on unsteady unidirectional flows of a non-Newtonian fluid[END_REF][START_REF] Erdogan | A note on an unsteady flow of a viscous fluid due to an oscillating plane wall[END_REF], while the diagrams of the transient solutions, as it was proved by graphical illustrations, are in close proximity of those obtained in [START_REF] Fetecau | Starting solutions for some unsteady unidirectional flows of a second grade fluid[END_REF] by a different technique. The required time to reach the steady-state for cosine oscillations of the wall is smaller than that for the sine oscillations of the wall. This time decreases if the frequency of the velocity of the boundary increases.

Statement of the problem

Consider a Newtonian fluid at rest over an infinitely extended plate which is situated in the (x, z) -plane. At time

+ = 0 t
the plate begins to oscillate in its plane and transmits the motion into the fluid. Its velocity is of the form

, ) t , y ( u ) t , y ( i v v = = (1)
where i denotes the unit vector in the x -direction.

The governing equation, in the absence of a pressure gradient along the x -axis and

neglecting body forces, is [2, 3] , 0 t , y ; y ) t , y ( u t ) t , y ( u 2 2 > ∂ ∂ ν = ∂ ∂ ( 2
)
where ν is the kinematic viscosity of the fluid.

)
where ω is the frequency of the velocity of the wall.

Solution of the problem

In order to determine exact solutions to the previous problems, describing the motion of the fluid at small and large times after the start of the boundary, the Laplace transform method is used. Consequently, applying the Laplace transform to Eq. ( 2) and using the initial and boundary conditions ( 3) and ( 4), we find that . q y exp q

U ) q , y ( u , ly respective , q y exp q Uq ) q , y ( u 2 2 2 2 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - ω + ω = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - ω + = (5)
Regarding the right parts of Eqs. (5) as a product of two functions ) q ( U 1 and ) q , y ( U 2 and using the known results [5]

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - πν = ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - ω ω = ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ω + ω ω + - - t 4 y exp t t 2 y q y exp L )}; t sin( ), t {cos( q ; q q L 2 1 2 2 2 2 1 (6) and { } ∫ ∫ - = - = * = ⋅ - t 0 2 1 t 0 2 1 2 1 2 1 1 ds, s) t (y, U (s) U ds s) (y, U s) (t U t) (y, U (t) U q) (y, U (q) U L ( 7
)
where * denotes the convolution product of two functions,

{ } ) q ( U L ) t ( U 1 1 1 - = and { } ) q , y ( U L ) t , y ( U 2 1 2 - =
, we find for the velocity field ) t , y ( u the simple expressions , ds s 4

y exp s s )] s t ( cos[ 2 Uy ) t , y ( u 2 t 0 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - ⎮ ⌡ ⌠ - ω πν = (8)
A c c e p t e d m a n u s c r i p t

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - ⎮ ⌡ ⌠ - ω πν = (9)
The solutions ( 8) and ( 9) can be also written into equivalent forms , ds s 4

y exp s s )] s t ( cos[ 2 Uy ds s 4 y exp s s )] s t ( cos[ 2 Uy ) t , y ( u 2 t 2 0 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - ⎮ ⌡ ⌠ - ω πν - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - ⎮ ⌡ ⌠ - ω πν = ∞ ∞ (10)
⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - ⎮ ⌡ ⌠ - ω πν - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - ⎮ ⌡ ⌠ - ω πν = ∞ ∞ (11
⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - ⎮ ⌡ ⌠ - ω πν - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν ω - ω = ∞ ν ω - (
⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - ⎮ ⌡ ⌠ - ω πν - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν ω - ω = ∞ ν ω - (13)
As soon as the velocity fields have been determined, the corresponding shear stresses can be immediately obtained by means of the known relation . y

) t , y ( u ) t , y ( ∂ ∂ µ = τ (14)
Consequently, introducing (12) and ( 13

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - ⎮ ⌡ ⌠ - ω πν ν µ + + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - ⎮ ⌡ ⌠ - ω πν µ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ π - ν ω - ω ν ω µ = τ ∞ ∞ ν ω - (15
⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - ⎮ ⌡ ⌠ - ω πν ν µ + + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - ⎮ ⌡ ⌠ - ω πν µ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ π - ν ω - ω ν ω µ - = τ ∞ ∞ ν ω - (16)
By letting now ∞ → t into these last four relations, we attain to the classical steadystate solutions (cf. with [START_REF]A note on unsteady unidirectional flows of a non-Newtonian fluid[END_REF] or [START_REF] Erdogan | A note on an unsteady flow of a viscous fluid due to an oscillating plane wall[END_REF], Eqs. ( 12) and ( 17))

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν ω - ω = ν ω - 2 y t cos Ue ) t , y ( u 2 y sc and ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν ω - ω = ν ω - 2 y t sin Ue ) t , y ( u 2 y ss (17) and ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ π - ν ω - ω ν ω µ = τ ν ω - 4 2 y t sin e U ) t , y ( 2 y sc , ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ π - ν ω - ω ν ω µ - = τ ν ω - 4 2 y t cos e U ) t , y (
)] s t ( cos[ 2 Uy ) t , y ( u 2 t tc ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - ⎮ ⌡ ⌠ - ω πν - = ∞ , (19)
⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - ⎮ ⌡ ⌠ - ω πν ν µ + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - ⎮ ⌡ ⌠ - ω πν µ - = τ ∞ ∞ (
⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - ⎮ ⌡ ⌠ - ω πν ν µ + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - ⎮ ⌡ ⌠ - ω πν µ - = τ ∞ ∞ (22)
as it results from Figs. 1, 2, 3 and 4, they are identical to those obtained in [START_REF] Fetecau | Starting solutions for some unsteady unidirectional flows of a second grade fluid[END_REF] 15) and ( 16) we can immediately determine the shear stress at the wall

ξ ξ ⎮ ⌡ ⌠ ν ω + ξ ξ π - = νξ - ∞ (23) , d ) y sin( e ) / ( U 2 ) t , y ( u t 2 0 2 4 ts ξ ξ ⎮ ⌡ ⌠ ν ω + ξ ξ ν ω π = νξ - ∞ (24) , d ) y cos( e ) / ( U 2 ) t , y (
, ds s s )] s t ( cos[ 2 U 4 t sin U ) t ( t wc ⎮ ⌡ ⌠ - ω πν µ - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ π - ω ν ω µ = τ ∞ (27)
for the cosine oscillations of the boundary, and

, ds s s )] s t ( sin[ 2 U 4 t cos U ) t ( t ws ⎮ ⌡ ⌠ - ω πν µ - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ π - ω ν ω µ - = τ ∞ (28)
for the sine oscillations of the boundary. The graphs of the wall shear stresses, corresponding to the cosine and sine oscillations of the boundary, are presented in Figs. 5 and6. For a better comparison their plots are presented together as well as separately.

The two oscillations have similar amplitudes and a phase shift that persists for all times.

Note that the improper integral in Eq. ( 27) is divergent at 0 t = . The some behavior is observed by analyzing Eq. (25) for 0 y = which strengthens the excellent agreement (as demonstrated by figure 3) between formulas (21) and (25). On the other hand formula (28), for the sine oscillations of the boundary, contains an improper integral which is convergent at 0 t = in agreement with the expression (26) obtained by the Fourier transform technique.

Finally, as it was to be expected, it is worthy pointing out that making 0 → ω into Eqs.

(12) and ( 15) we obtain the classical solutions , t 2

y erfc U t 2 y erf 1 U ) t , y ( u ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - = , t 4 y exp t U ) t , y ( 2 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - πν µ - = τ (29)
corresponding to the first problem of Stokes [START_REF] Zierep | Similarity Laws and Modeling[END_REF]. In order to obtain these results, we have

⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ν νπ ν + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - ν ν ν - = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - ⎮ ⌡ ⌠ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ν νπ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - ⎮ ⌡ ⌠ ∞ ∞ (30)
In the above relations, ) ( erf ⋅ and ) ( erfc ⋅ are the error function and the complementary error function of Gauss.

Conclusions

In this note we have established new exact solutions corresponding to the second problem of Stokes for Newtonian fluids. These solutions, unlike those obtained by Erdogan [START_REF] Erdogan | A note on an unsteady flow of a viscous fluid due to an oscillating plane wall[END_REF], are simpler and directly presented as a sum of steady-state and transient solutions. The steady-state solutions can be brought to the classical forms, while the graphs of the transient solutions, as depicted in Figs.1, 2, 3 and 4, are in close proximity to those corresponding to the similar solutions obtained in [START_REF] Fetecau | Starting solutions for some unsteady unidirectional flows of a second grade fluid[END_REF] by a different technique.

For large values of t, when the transients disappear, the motion of the fluid is described by the steady-state solutions which are periodic in time and independent of the initial conditions. However, they satisfy the governing equation and the associated boundary conditions.

The variations of the wall shear stresses for the cosine and sine oscillations of the boundary are presented in Figs. 5 and6 for several time intervals. The two oscillations have similar amplitudes and a phase shift that persists for all times. The two sets of formulas (one obtained by Laplace transform, the other by Fourier sine transform) are in perfect agreement.

A c c e p t e d m a n u s c r i p t 2 Corina

 2 Fetecau, D. Vieru and C. Fetecau

u s c r i p t 6 Corina

 6 Fetecau, D. Vieru and C. Fetecau which are periodic in time and independent of the initial conditions. As regards the

Finally, we would

  like to know the required time to attain the steady-state for the cosine and sine oscillations of the boundary. The variations with the distance from the wall of the corresponding starting and steady-state velocities for the oscillations of in Figs.7 and 8. It is clearly seen from these figures that the required time to reach the steady-state for sine oscillations of the wall is greater than that for the cosine oscillations of the boundary. Furthermore, this time decreases if the frequency ω of the velocity increases.

A c c e p t e d m a n u s c r i p t

c c e p t e d m a n u s c r i p t

	, 0 , 0 , y , y ( ( u ( u ( u u	, 0 0 > > A , t all for ) t cos( ω t all for ) t sin( ω , 0 y for 0 U U > = = = ) ) ) 0 t t , y as 0 ) t ∞ → →	The boundary and initial conditions are	(3a) (4a) (4b) (4c
	( u	, y	0	)	=	0	for	y	>	, 0	(3b)
	( u	, y	t	)	→	0	as	y	→	, ∞	(3c)
	or

Acknowledgement

The authors acknowledge the support from the Ministry of Education and Research, trough CEEX-program, Grant: CERES-2-Cex06-11-56/25-07-2006, no.1/24.07.2006. We also thank the referees for their valuable comments and suggestions which helped us improve the quality of the present paper.