Lev Truskinovsky

Anna Vainchtein

The origin of nucleation peak in transformational plasticity

Keywords: Martensitic phase transitions, Lattice models, Nonlocal interactions, Peierls-Nabarro landscape, Nucleation

A typical stress-strain relation for martensitic materials exhibits a mismatch between the nucleation and propagation thresholds leading to the formation of the nucleation peak. We develop an analytical model of this phenomenon and obtain speciÿc relations between the macroscopic parameters of the peak and the microscopic characteristics of the material. Although the nucleation peak appears in the model as an interplay between discreteness and nonlocality, it does not disappear in the continuum limit. We verify the quantitative predictions of the model by comparison with experimental data for cubic to monoclinic phase transformation in NiTi.

Introduction

In displacement-controlled experiments shape memory alloys and other martensitic materials display a nucleation peak: prior to nucleation of a new phase the load reaches a maximum but then drops to a distinctly lower value. The subsequent plateau is associated with phase boundaries propagating along the specimen at essentially constant stress [START_REF] Lexcellent | Internal loops in pseudoelastic behaviour of ti-ni shape memory alloys: experiment and modelling[END_REF]Shaw andKyriakides, 1995, 1997b;[START_REF] Sun | An inclusion theory for the propagation of martensite band in NiTi shape memory alloy wires under tension[END_REF]. The nucleation peak is not unique to transformational plasticity and is also observed during initiation of a conventional plastic deformation in mild steels [START_REF] Butler | L uders front propagation in low carbon steels[END_REF][START_REF] Hall | Yield Point Phenomena in Metals and Alloys[END_REF][START_REF] Kyriakides | On the propagation of luders bands in steel strips[END_REF][START_REF] Froli | A mechanical model for the elastic-plastic behavior of metallic bars[END_REF], where it has been attributed to the fact that the stress required to release the trapped dislocations is higher than the stress needed to sustain their motion [START_REF] Cottrell | Dislocation theory of yielding and strain ageing of iron[END_REF][START_REF] Johnston | Dislocation velocities, dislocation densities and plastic ow in lithium uoride crystals[END_REF]. For martensitic materials the nucleation-induced load drop has been observed in 3D numerical simulations based on various plasticity-like phenomenological models (Shaw and Kyriakides, 1997a;[START_REF] Kyriakides | On the propagation of luders bands in steel strips[END_REF][START_REF] Sun | An inclusion theory for the propagation of martensite band in NiTi shape memory alloy wires under tension[END_REF] but due to the complexity of these models, the physical parameters responsible for the size of the peak have not been identiÿed. At a qualitative level, the nucleation peak in these materials has been associated with the presence of su ciently ÿne grains and heuristically linked to the strong locking of phase boundaries and the relative ease of their glide upon release (Shaw andKyriakides, 1995, 1997b).

In this paper we develop an analytical model of the nucleation peak phenomenon in martensites: our model supports the intuition developed in plasticity theory and adapts it to the case when the principal carriers of inelastic deformation are phase boundaries. Speciÿcally, we consider a prototypical mass-spring system consisting of rigid elements (crystal planes) connected by bi-stable elastic springs representing transforming shear layers. To mimic the three-dimensional nature of the actual problem, we complement the up-down-up interactions between the nearest neighbors (NN) by a harmonic interaction of the next-to-nearest neighbors (NNN). The bi-stable discrete models without NNN interactions (e.g. M uller [START_REF] Villaggio | A model for an elastic-plastic body[END_REF][START_REF] Fedelich | Hysteresis in discrete systems of possibly interacting elements with a two well energy[END_REF]Puglisi andTruskinovsky, 2000, 2002a, b) capture many important features of transformational plasticity but fail to predict the peak phenomenon. Recent numerical studies of the models incorporating NNN interactions showed that the nucleation peak can be recovered (e.g. [START_REF] Ye | E ect of phonon anomalies on the shear response of martensitic crystals[END_REF][START_REF] Triantafyllidis | On higher order gradient continuum theories in 1-D nonlinear elasticity. Derivation from and comparison to the corresponding discrete models[END_REF][START_REF] Rogers | Discretization and hysteresis[END_REF][START_REF] Froli | A mechanical model for the elastic-plastic behavior of metallic bars[END_REF][START_REF] Pagano | A simple model for phase transformations: from discrete to the continuum problem[END_REF]; none of these models, however, have been developed analytically to the extent that they could explain the necessity of the peak phenomenon and identify the microparameters controlling the size of the stress drop.

We begin by ÿnding the limits of instability of a homogeneous state and determine an analytical expression for the nucleation threshold in the most general case. In order to obtain an analytical characterization of the propagation threshold we use a piecewise linear approximation for the NN interactions. This simpliÿcation allows us to reconstruct the non-equilibrium Peierls-Nabarro landscape for the propagating phase boundaries and compute the martensitic analog of the Peierls stress. We show that the presence of nonlocal interactions in the discrete model makes nucleation and propagation thresholds di erent and then prove that the nucleation peak does not disappear in the continuum limit. An important question is whether the quasi-continuum strain-gradient approximation of the discrete model [START_REF] Mindlin | Second gradient of strain and surface-tension in linear elasticity[END_REF][START_REF] Triantafyllidis | On higher order gradient continuum theories in 1-D nonlinear elasticity. Derivation from and comparison to the corresponding discrete models[END_REF] is capable of capturing the peak phenomenon. We show that the nucleation threshold is approximated well only in the case of long-wave instability (see also [START_REF] Triantafyllidis | The in uence of scale size on the stability of periodic solids and the role of associated higher order gradient continuum models[END_REF] and that although the peak is captured, the stress drop is grossly exaggerated. Finally, we verify the quantitative prediction of the model by using the experimental data for cubic to monoclinic phase transformation in NiTi wires. The comparison with experiment leads to the bounds for the measure of nonlocality which are in a good agreement with the independent estimates based on a realistic interatomic potential.

The proposed model stays in the same prototypical relation to transformational plasticity as the well known Frenkel-Kontorova model to classical plasticity. A nontrivial formal correspondence between the two models in the case of inÿnite systems has been established by [START_REF] Truskinovsky | Peierls-Nabarro landscape for martensitic phase transitions[END_REF].

The structure of the paper is as follows. In Section 2 we formulate the discrete problem for a ÿnite system and specify the boundary conditions. In Section 3 we analyze stability of the single-phase equilibrium in the general setting discovering the possibility of both macro and microinstabilities. In Section 4 we introduce a piecewise linear model for the local interactions and obtain an explicit representation of a generic metastable equilibrium state. In Section 5 we study the energy barriers between the neighboring metastable states and show that propagation always reduces to a succession of transformations inside individual elements, while nucleation can involve transformation of several elements at once. Section 6 contains the derivation of the explicit formulae for the size of the nucleation peak in both discrete and continuum problems. In Section 7 we establish a correspondence between the microparameters of the lattice and the experimental measurements of both the stress drop and the size of the nucleation band. Finally, Section 8 contains comparison of the discrete and strain-gradient models showing that the agreement is at most qualitative. The conclusions are summarized in Section 9.

The model

Consider a system of N + 1 particles linked to their nearest and next-to-nearest neighbors by elastic springs (see Fig. 1). Let u k , 0 6 k 6 N , be the displacements of the particles with respect to a load-free homogeneous reference conÿguration with spacing . Denote the strain in the kth NN spring by w k = (u k -u k-1 )= . Then the total energy of the system can be written as

= N k=1 1 (w k ) + 2 N -1 k=1 2 w k+1 + w k 2 + B (w 1 ; w N ); (1)
where 1 (w) and 2 (w) are the energy densities of the NN and NNN interactions, respectively. The term B corresponds to the energy of the boundary elements. We assume that the chain is placed in a hard device with the total displacement d:

i i-1 i+1 0 1 N-1 N ε 2ε NNN spring NN spring
u N -u 0 = N k=1 w k = d:
(2)

Due to the nonlocality of the model, the boundary condition (2) must be complemented by other constraints. For instance, to mimic an "extra hard device" one can additionally impose conditions w 1 = 0 and w N = 0, making the term B irrelevant. Alternatively, one may consider a "zero-moment device" by assuming that B = 0. Let f i (w) = i (w), i = 1; 2, denote the forces in NN and NNN springs, respectively. The equations governing the equilibrium of the interior particles with 2 6 k 6 N -1 have the form

f 1 (w k ) + f 2 w k+1 + w k 2 + f 2 w k + w k-1 2 = F; ( 3
)
where F is the total force in the system. The natural boundary conditions read

f 1 (w 1 ) + f 2 w 2 + w 1 2 + 1 @ B @w 1 = F; f 1 (w N ) + f 2 w N + w N -1 2 + 1 @ B @w N = F; (4)
An additional assumption adopted in what follows, B = 2 (w 1 ) + 2 (w N ) ;

(5) means that the boundary NNN springs are cut in half and reconnected parallel to the NN springs (see the dashed springs in Fig. 1). The formal advantage of this choice of B is that the boundary equations ( 4) can be included into the bulk equations (3) if we additionally assume the existence of ÿctitious 0th and (N + 1)th springs satisfying

w 0 = w 1 ; w N +1 = w N : (6)
The real advantage of (4) and ( 5), however, is that the corresponding boundary conditions ensure the existence of a trivial solution with the uniform strain distribution [START_REF] Triantafyllidis | On higher order gradient continuum theories in 1-D nonlinear elasticity. Derivation from and comparison to the corresponding discrete models[END_REF][START_REF] Charlotte | Linear elastic chain with a hyper-pre-stress[END_REF]. While in the rest of the paper we will be using mostly (5), the e ect of switching to B = 0 is brie y discussed in Section 6.

Nucleation threshold

It is easy to see that for arbitrary spring potentials the trivial solution of the problem (3), (4) and ( 5) is given by w k = d=L, where L = N . The homogeneous response of the system is then characterized by the formulae F = f 1 (d=L) + 2f 2 (d=L) and = L( 1 (d=L) + 2 2 (d=L)). To analyze the stability of this solution we introduce the tangential moduli of NN and NNN springs: K(w) = 1 (w) and (w) = 2 (w)=2.

The homogeneous conÿguration is stable if and only if the quadratic form Bv • v with

B =              K + 3 0 : : : 0 K+ 2 . . . 0 0 . . . . . . . . . 0 . . . 0 K + 2 0 : : : 0 K+ 3 :              (7)
is positive deÿnite for all v = 0 such that

N k=1 v k = 0: (8)
Finding the boundaries of stability in the space of parameters is important for the subject of the paper because they coincide with the limits of a barrierless nucleation. Before addressing the stability problem systematically, we observe that K + 4 ¿ 0 and 6 0 are su cient for stability. Indeed, in this case all terms in the quadratic form

Bv • v = (K + 4 ) N k=1 v 2 k - N -1 k=1 (v k+1 -v k ) 2 : (9)
are nonnegative. Similarly, by writing the quadratic form as

Bv • v = K N k=1 v 2 k + N -1 k=1 (v k+1 + v k ) 2 + 2 (v 2 1 + v 2 N ):
we obtain that conditions K ¿ 0 and ¿ 0 are also su cient for stability. In terms of the main nondimensional parameter of the problem,

= K 4 ; (10)
these stability intervals can be written as -∞ ¡ ¡ -1 and 0 ¡ ¡ ∞ and hence the instability limits are located in the interval

-1 6 6 0: (11)
To ÿnd the exact locations of the stability boundaries, consider the (zero) eigenvalue problem Bv = 0. In the bulk of the chain (1 6 k 6 N -1) this means

(K + 2 )v k + v k-1 + v k+1 = 0: (12)
On the boundaries we obtain

(K + 3 )v 1 + v 2 = 0; (13) (K + 3 )v N + v N -1 = 0: (14)
The eigenvector v must be nonzero and satisfy (8). We seek solution in the form v k = k and obtain the characteristic polynomial 2 + (4 + 2) + 1 = 0: (15)

This means that 1; 2 = -1 -2 ± 2 ( + 1), and according to (11), we have to consider the following cases:

Special case = -1: Here 1 = 2 = 1 and v k = A 1 k + A 2 . The boundary equations ( 13) and ( 14) reduce to v 1 = v 2 and v N = v N -1 , implying that A 1 = 0 and thus v k = A 2 . The constraint (8) implies that A 2 = 0. Since there are no nontrivial solutions, there is also no stability change.

Special case = 0: Here 1 = 2 = -1 and v k = (A 1 k + A 2 )(-1) k . Applying ( 13) and ( 14), we obtain A 1 = A 2 = 0 which again means that there is no change of stability.

Generic case -1 ¡ ¡ 0. We ÿrst observe that in this interval 1; 2 can be written as 1; 2 = e ±i! ; where ! is deÿned by

= -cos 2 ! 2 : (16)
Applying the boundary equation ( 13), we obtain (up to a multiple)

v k = cos[(k -1=2)!]:
The second boundary equation ( 14) implies 2 sin(!=2) sin(N!)=0, and since 0 ¡ ! ¡ , we obtain sin(N!) = 0 and

! = n N ; (17)
for 1 6 n 6 N -1. All these solutions are nontrivial and satisfy (8). Therefore the instability of the trivial solution can take place at any of the bifurcation points

K + 4 cos 2 n 2N = 0:
The corresponding unstable modes are

v k = cos[(k -1 2 ) n N ]: (18)
To locate the stability boundary, we begin with the case ¡ 0. Then

K + 4 cos 2 n 2N ¿ K + 4 cos 2 2N and hence K + 4 cos 2 2N ¿ 0 (19)
is both necessary and su cient for stability. The instability develops through the growth of the long-wave mode

v k = cos[(k -1 2 ) N ]:
As N tends to inÿnity, the unstable wave length also becomes inÿnite and (19) reduces to K + 4 ¿ 0. At ÿnite N we have K + 4 cos 2 2N ¿ K + 4 , meaning that the discrete homogeneous conÿguration may still be stable when the macroscopic modulus E = K + 4 is already negative. One can say that in the discrete problem the macroscopic instability is delayed due to the ÿnite size of the system (see also [START_REF] Triantafyllidis | On higher order gradient continuum theories in 1-D nonlinear elasticity. Derivation from and comparison to the corresponding discrete models[END_REF].

The situation is di erent when ¿ 0. In this case

K + 4 cos 2 n 2N ¿ K + 4 cos 2 (N -1) 2N = K + 4 sin 2 2N ;
and hence the necessary and su cient condition for stability is

K + 4 sin 2 2N ¿ 0: (20)
The instability develops through the growth of the short-wave mode

v k = (-1) k sin k -1 2 N :
In the limit of inÿnite N , the wave length approaches the interatomic distance. Also (20) becomes K ¿ 0, and since now K + 4 sin 2 ( =2N ) ¡ K + 4 , this microscopic instability develops before the macroscopic (or homogenized) system becomes unstable (see also [START_REF] Triantafyllidis | The in uence of scale size on the stability of periodic solids and the role of associated higher order gradient continuum models[END_REF]. In the context of martensitic phase transitions, this e ect may be linked to the observation of the pre-martensitic tweed microstructures [START_REF] Kartha | Disorder driven pretransitional tweed pattern in martensitic transformations[END_REF].

We now summarize the necessary and su cient conditions for stability of the trivial solution. In terms of elastic moduli K and we obtain the following stability intervals

K + 4 cos 2 2N ¿ 0; K + 4 sin 2 2N ¿ 0: (21)
In the limit of inÿnite N the inequalities (21) reduce to the known conditions K ¿ 0, K + 4 ¿ 0 (e.g. [START_REF] Mindlin | Second gradient of strain and surface-tension in linear elasticity[END_REF]. In terms of the nondimensional parameter the stability intervals take the form

¡ -cos 2 2N ( ¡ 0); ¿-sin 2 2N ( ¿ 0); (22)
which in the limit of inÿnite N gives ¡-1 ( ¡ 0) and ¿ 0 ( ¿ 0). The structure of the stability domain in the plane K -is illustrated in Fig. 2a for both ÿnite and inÿnite N ; in Fig. 2b, showing the plane -N -1 , we illustrate the dependence of the stability boundaries on the size of the system. For the NNN system with the "moment free" boundary conditions ( B = 0) similar analysis of the stability for trivial equilibria (which now have boundary layers) can be found in [START_REF] Charlotte | Linear elastic chain with a hyper-pre-stress[END_REF].

Nontrivial solutions and metastability

To model martensitic phase transitions, we assume that each NN spring has a double-well energy generating a non-monotone "up-down-up" force-strain relation depicted in Fig. 3. If the negative slope of the force-strain relation for the NN spring is su ciently steep, it is easy to show that the total strain can reach the threshold where the homogeneous phase is absolutely unstable. To study the resulting nucleation and the subsequent growth of the new phase, we need to describe the postcritical behavior of the system. More precisely, in the domain of loadings where trivial solution is unstable, we need to ÿnd the nontrivial solutions of the equilibrium equations (3) corresponding to either local or global minima of the energy.

Although the instability limits for the trivial solution can be found in the general case, we succeeded in obtaining a complete analytical description of the postcritical behavior only in the case when the functions f 1 (w) and f 2 (w) are either linear or piecewise linear. Speciÿcally, we considered

f 1 (w) = Kw for w ¡ w c K(w -a) for w ¿ w c (23)
and

f 2 (w) = 2 w: (24)
The parameters of the NN potential are the critical strain w c , transformation strain a and elastic modulus K (see Fig. 3). The linear NNN spring is characterized by the elastic modulus which will also be used as the measure of nonlocality. We assume that K ¿ 0 and ¡ 0, which is suggested by the linearization of the Lennard-Jones potential (see [START_REF] Charlotte | Linear elastic chain with a hyper-pre-stress[END_REF]; to ensure that independently of N the homogeneous states are stable in their domain of deÿnition, we also assume that K +4 ¿ 0. In terms of the nondimensional parameter , the assumptions on the moduli can be summarized as

-∞ ¡ ¡ -1: (25)
It is not hard to see that under these assumptions the microinstability is excluded while the macroinstability always takes place at w = w c .

For the piecewise linear model ( 23), ( 24) the total energy (1) reduces to

= 1 2 Bw • w -q • (w -w c ) ; (26)
where the vector q k = KaÂ(w k -w ck ) prescribes distribution of phases, Â(x) is a unit step function and w ck = w c , k = 1; : : : ; N . The equilibrium equations (3) together with the boundary conditions (6) can then be rewritten as

Bw = F + q; ( 27
)
where F k = F, k = 1; : : : ; N . To eliminate the redundant parameters we rescale the variables. By selecting L = N as the length scale and K as the scale of force, we deÿne

u k = u k N ; d = d N ; F = F K ; = KN ; B = 1 K B: (28)
Unless specially mentioned, in what follows, we will be using only rescaled variables with the bars dropped. The dimensionless problem depends on the two main parameters: and N .

To ÿnd the metastable states we can simply prescribe the distribution of phases q and solve the linear problem ( 27). The resulting equilibrium conÿgurations are automatically local minimizers of the energy (26), because the sti ness matrix B is positive deÿnite. To obtain the global minima of the energy, we need additionally to minimize with respect to the phase geometry q.

The solution of the linear system ( 27) can be written as

w = B -1 F + B -1 q = w 0 + w 1 : (29)
The ÿrst term in (29) corresponds to the uniform conÿguration with all springs in the ÿrst phase. Without explicitly computing B -1 one can see that

w 0 k = F 1 + : (30)
It is convenient to express w 1 not in terms of the variables q, but in terms of their discrete derivatives p deÿned by

p i = q i+1 -q i a : (31)
Notice that these relations can always be explicitly inverted, yielding

q k = q 1 + a N i=1 p i Â(k -i -1): (32)
The physical meaning of the variables p i is clear from the representation (33)

p i =        1 if w i ¡ w c
To simplify the subsequent formulae we also set p N = 0.

In the interior nodes (2 6 k 6 N -1) w 1 k must satisfy the di erence equation

1 + 1 2 w 1 k + 1 4 (w 1 k+1 + w 1 k-1 ) = q 1 + a i p i Â(k -i -1): (34)
We seek the general solution of (34) in the form

w 1 k = w h k + w in k ; (35)
where w h k satisÿes the homogeneous equation and w in k is a particular solution. The general solution of the homogeneous problem can be represented [START_REF] Mickens | Di erence Equations: Theory and Applications[END_REF] as a linear combination of k 1 and k 2 , where 1; 2 are the roots of the characteristic polynomial (15). By writing 1; 2 = e ± , where = 2 arccosh | |;

(36) we obtain

w h k = C 1 e k + C 2 e -k : (37)
The constants C 1 and C 2 are to be found from the boundary conditions (6).

To obtain a particular solution of (34) with q 1 = 0 and p j = 0 for j = i, we observe that in this case the right hand side of (34) is zero for k 6 i and constant for k ¿ i +1. Therefore we can write

w pi k =    w p- k = A 1 e k + A 2 e -k for k 6 i w p+ k = A 3 e k + A 4 e -k + a p i 1 + for k ¿ i; (38)
and then "glue" the two sides together by requiring that

w p- i = w p+ i ; w p- i+1 = w p+ i+1 : (39)
Since we are looking for a particular solution, we may always let A 1 = A 2 = 0, and thus consider w p- k = 0. Solving (39) for A 3 and A 4 , we obtain

w pi k = Â(k -i -1=2) 1 - cosh[(k -i -1=2) ] cosh( =2) ; where = a + 1 (40)
is the macroscopic transformation strain. Finally, by superposition, we obtain the particular solution in the form

w in k = q 1 1 + + N i=1 p i Â(k -i -1=2) 1 - cosh[(k -i -1=2) ] cosh( =2) : (41)
The general solution of ( 34) is now given by ( 35), ( 37) and (41). Applying the boundary conditions (6), we obtain

w 1 k = q 1 1 + + N i=1 p i sinh[(N -i) ] cosh[(k -1=2) ] cosh( =2) sinh(N ) + Â(k -i -1=2) 1 - cosh[(k -i -1=2) ] cosh( =2) : (42)
Before combining ( 30) and ( 42), we notice that the relation between the force F and the total displacement d can be written in the form:

F = 1 + d - l N ; ( 43
)
where l is the number of springs in phase II related to q 1 and p via

l = Nq 1 a + N i=1 (N -i)p i : (44)
Finally, substituting (43), ( 44) into ( 29), ( 30) and ( 42), we obtain the representation for a generic metastable conÿguration

w k = d + N i=1 p i sinh[(N -i) ] cosh[(k -1=2) ] cosh( =2) sinh(N ) + i N -1 + Â(k -i -1=2) 1 - cosh[(k -i -1=2) ] cosh( =2) : (45)
To compute the energy we introduce a new variable w with components wk = w k -d representing the p-dependent part of the strain ÿeld w. Then ( 26) can be rewritten as

= 1 + 2 d 2 -a(d -w c ) l N + ˜ ; (46)
where

˜ = 1 N 1 2 B w -q • w (47)
is the contribution due to phase changes. Observing that in equilibrium ˜ = -1 2N q • w and using (32) to eliminate q, we obtain

= 1 + 1 2 d 2 -(d -w c ) l N + 2 2N Jp • p ; (48)
where the kernel J is given by

J ki = sinh[(N -i) ] sinh[ k] sinh(N ) sinh + k i N -1 + Â(k -i -1) k -i - sinh[(k -i) ] sinh : (49)
For a given distribution of phases the loading parameter d cannot take arbitrary values since the strains must satisfy the constraint

p i (w i -w c ) ¡ 0; p i (w i+1 -w c ) ¿ 0 for i : p i = 0: (50)
Conditions (50) generate bounds on d(i 1 ; i 2 ; : : : ; i n ), where i 1 ; i 2 ; : : : ; i n are locations of the phase boundaries, parameterizing a particular metastable branch. For example, if n = 1 we obtain the limits d

-(i) ¡ d(i) ¡ d + (i), where d ± (i) = w c + sinh[(N -i) ] cosh[(i ± 1=2) ] cosh 2 sinh(N ) + i N -1 : (51)
Formulae ( 43), ( 48) and (51) describing single-interface solutions are illustrated in Fig. 4. An equilibrium branch parametrized by i begins at d = d -(i), where w i = w c , and terminates at d = d + (i), where w i+1 = w c . In the interval d -(i) ¡ d ¡ d + (i) the force is linear, while the energy is quadratic. Due to the symmetry of the problem, the ith branch is indistinguishable from the (N -i)th branch. Equilibria with two interfaces are illustrated in Fig. 5 for a chain with N = 6. Notice that unlike the case without NNN interactions (e.g. Puglisi and Truskinovsky, 2002b), here not only the volume fractions but also the actual locations of the interfaces distinguish the branches. For example, the blow-up in Fig. 5b shows that the energy of the branch p=(0 -1 0 1 0 0) is higher than the energy of a branch p=(-1 0 1 0 0 0) although both branches have the same fraction of phase II.

d b (-1 1 0 0 0 0) (0 -1 1 0 0 0) (0 0 -1 1 0 0) (0 0 0 -1 1 0) (-1 0 1 0 0 0) (0 -1 0 1 0 0) (0 0 -1 0 1 0) (1 0 0 0 -1 0) (-1 0 0 1 0 0) (0 -1 0 0 1 0) (1 0 0 -1 0 0) (0 1 0 0 -1 0) p = (-1 0 0 0 1 0) (1 0 -1 0 0 0) (0 1 0 -1 0 0) (0 0 1 0 -1 0) (1 -1 0 0 0 0) (0 1 -1 0 0 0) (0 0 1 -1 0 0) (0 0 0 1 -1 0) (0 -1 0 1 0 0) (-1 0 1 0 0 0) Ψ Ψ (a) (b)
To ÿnd the global minimum of the energy we can ÿrst minimize among metastable solutions with a given number of interfaces n. The resulting lower envelopes ˆ (n; d) are shown in Fig. 6 for N =6 and n 6 4. One can see that the single-interface solutions have the lowest energy, which is expected since our choice of nonlocal interactions leads to the penalization of interfaces. Similar observation has been made by [START_REF] Rogers | Discretization and hysteresis[END_REF] for a discrete model with long-range forces.

Peierls-Nabarro landscape

To understand the degree of relative stability of various metastable states we need to evaluate the energy barriers. Since conÿgurations with at most one interface are expected to have the lowest energy, it is instructive to start the analysis with the corresponding section of the energy landscape. By focusing on the conÿgurations with a single phase boundary we obtain an analog of the Peierls-Nabarro landscape known in the theory of dislocations; it is also consistent with experimental observations for martensitic materials (e.g. [START_REF] Krishnan | Stress induced martensitic transformations[END_REF].

Consider a generic metastable conÿguration with a single phase boundary at k = i. To evaluate the barrier between this and the neighboring local minimum we need to choose a path connecting the conÿguration w k (i) to the conÿguration w k (i + 1) with one extra spring in phase II. Since the (i + 1)th spring must change phase, it is natural to choose the strain w i+1 as a parameter and minimize the total energy with respect to all w k with k 6 i and k ¿ i + 2. We obtain the system of equations

1 + 1 2 w k + 1 4 (w k+1 + w k-1 ) = F + a for k 6 i F for k ¿ i + 2; ( 52
)
which must be supplemented by the boundary conditions ( 2) and ( 6). The conÿgurations satisfying ( 52), ( 2) and ( 6) have the lowest energy among all states with a given w i+1 . This follows from the positive deÿniteness of the matrix B (i+1) obtained from (7) by deleting the (i + 1)th row and the (i + 1)th column. It is convenient to write the explicit representation for these conÿgurations by using instead of w i+1 another parameter :

w k ( ) =      d - N -1 + C 1 (sinh( k) -sinh( (k -1))) ; k6 [ ] + 1 d - N + C 2 (e k -e (k-1) -e (2N -k) + e (2N -k+1) ) ; k ¿ [ ] + 1; (53)
where [ ] denotes the integer part of . It can be directly checked that (53) solves ( 52) for k 6 [ ] and k ¿ [ ] + 2, with [ ] = i and the total force F given by

F = 1 + d -N : (54)
By matching the strains at k = i + 1 = [ ] + 1 and imposing (2), we compute the values of the constants

C 1 = 2{e (2N -[ ]) [ -[ ] -e ( [ ] -1)] -e ([ ]+1) [ -[ ] -1 -e ( -[ ])]} (e 2 -1)(e 2N -1)
and

C 2 = 2{( -[ ]) sinh[([ ] + 1) ] -( -[ ] -1) sinh[[ ] ]
} (e 2 -1)(e 2N -1) :

In particular, we obtain an explicit relation between the parameters w i+1 and :

w [ ]+1 ( ) = d - N -1 + C 1 (sinh[([ ] + 1) ] -sinh[[ ] ]): (55)
One can see that parameter w i+1 oscillates as the function of i with period 1, while increases monotonically. Since the integer values of correspond exactly to metastable conÿgurations, the function ( ) obtained by substituting (53) in ( 46) is exactly the Peierls-Nabarro (PN) potential of our system. In order to move from one valley of this potential at = i to the neighboring one at = i + 1 the system must overcome the Peierls barrier. It is located at = i deÿned by

w i+1 ( i ) = w c : (56)
The height of the Peierls barrier i→i+1 = ( i ) -(i) can be explicitly computed from ( 46) and (47).

The typical structure of the PN landscape is illustrated in Fig. 7 for a chain with N = 10. One can see that at d = 0:5 the metastable conÿgurations correspond to the integer values = 2 (or 3), = 1 (or 4), = 0, = 5 and = 6. The corresponding strain proÿles are shown in Fig. 8 for the metastable states = 0, = 1 and = 2 and for the saddle point conÿgurations with = 0 and = 1 .

Propagation threshold and nucleation peak

It is realistic to assume that the system driven by external loading d remains in a local minimum conÿguration until the minimal energy barrier around this state reaches a critical threshold H determined by the level of uctuations (imperfections). In particular, the maximal delay strategy, associated with the gradient-ow dynamics, requires that the system stays on a given metastable branch until it becomes absolutely unstable (H = 0) and then evolves towards the nearest local minimum along the path of steepest descent. Various resulting force-strain paths for the NN system are discussed in Puglisi and Truskinovsky (2002a, b). For the NNN system the two characteristic paths with and without a threshold are shown in Fig. 9. One can see that as the elongation increases the system initially stays in the trivial conÿguration but eventually reaches the state where the smallest energy barrier becomes equal to H . Then nucleation takes place and the system escapes from the local minimum through the ÿrst saddle point with a subcritical height. After the nucleation event the phase boundary propagates along the chain in a stick-slip fashion, with the system getting temporarily trapped in each of the single interface metastable equilibria (parametrized by i). The resulting graphs of the force F(d) exhibit characteristic serrations. A visible nucleation peak originates from the fact that the nucleation event involves two springs transforming at once while the propagation involves only one spring changing phase at a time.

To understand why the ÿrst two springs change phase simultaneously, which never happens in the NN system (Puglisi and Truskinovsky, 2002a, b), we need to compare the height of the energy barrier for the transition 0 → 1 with the height of the barriers for several subsequent transitions. The height of the barriers is shown in Fig. 10 for the chain with N = 10 and = -2. Notice that the barrier for the transition 1 → 2 is lower than the barrier for the transition 0 → 1, moreover, the barrier 1 → 2 vanishes at d + (1) and beyond this point the branch i = 1 does not exist any more (transition 0 → 1 deteriorates into 0 → 2). A similar calculation for N = 20 and = -2 shows that the barrier for the transition 0 → 1 is higher than two subsequent barriers for the transitions 1 → 2 and 2 → 3; as the last two barriers vanish at su ciently large d, the transition 0 → 1 ÿrst deteriorates into 0 → 2 and then into 0 → 3. One can see that in this case nucleation event involves simultaneous transformation of three springs.

Our computations show that "massive" nucleation and the associated nucleation peak phenomenon occur only when NNN interactions are su ciently strong (large | -1 |). Physically, it is the consequence of the nonlocal character of the model. The nonlocality does not have to be of the NNN type: a di erent nonlocal model with long-range interactions also exhibits the nucleation peak [START_REF] Rogers | Discretization and hysteresis[END_REF].

To understand what is going on, we observe that before the nucleation all springs are stretched uniformly, whereas after the nucleation, due to the presence of the internal boundary layers, the springs that are closer to the phase boundary have higher strain and hence are closer to the critical threshold than the springs far away. This facilitates the subsequent switching events and results in the smaller force required for the propagation of a phase boundary comparing to nucleation. In the system without nonlocality ( = -∞), all springs outside the interface are stretched uniformly, and therefore the propagation of the interface does not take place until the critical strain is reached in all springs simultaneously. In this case the nucleation peak is absent: phase propagation e ectively reduces to successive nucleation events in the shorter and shorter chains which requires the same critical force.

To obtain the upper bound for the size of the nucleation peak it is enough to consider the maximum delay strategy. In our model barrierless nucleation takes place when the force reaches the spinodal limit F max = (1 + -1 )w c . On the other hand, the advance of the interface from k = i to k = i + 1 takes place at F(i; N ) = F M + F P (i; N ), where

F M = 1 + w c - a 2 (57)
is the Maxwell force and F P (i; N ) is the Peierls force given by

F P (i; N ) = F(d + (i); i) -F M = a cosh[( 1 2 + i) ] sinh[ (N -i)] cosh 2 sinh(N ) - 1 2 : (58)
One can show that the function F P (i; N ) depends on i weakly away from the narrow boundary layers near i = 1 and i = N -1. In the limit of inÿnite N the Peierls force (58) approaches the constant value [START_REF] Truskinovsky | Peierls-Nabarro landscape for martensitic phase transitions[END_REF] lim

N →∞ F P (i; N ) = F P = a 2 1 + : (59)
As the NNN interactions get weaker ( → -∞), the Peierls force tends to the spinodal limit F max -F M and the nucleation peak disappears. The conÿguration of the resulting hysteresis loop in the continuum limit is shown in Fig. 11. Although the serrations disappeared, the nucleation peak remains ÿnite with the force dropping by the amount

= 1 + w c -F P = 1 2 a 1 - 1 + : ( 60
)
This quantity is positive as long as -∞ ¡ ¡ -1 and is always less than a=2-the di erence between the spinodal and Maxwell forces. The half-height of the narrow part of the hysteresis loop is given in the continuum limit by (59).

To estimate the number of springs involved in the nucleation event, we recall that the nucleation always takes place at d = w c . Therefore we must ÿnd the metastable branch with the smallest nonzero i which is deÿned at this value of d. Setting d ± from (51) equal to w c , we obtain the equation for the number of participating springs i sinh

[(N -i) ] cosh[(i ± 1=2) ] cosh 2 sinh(N ) = 1 - i N : (61)
For su ciently large N and generic the solution of this equation can be approximated by

i nuc = N 1 + e : (62)
In particular, (62) implies that in the continuum limit the number of springs involved in the nucleation event is inÿnite. The dimensionless size of the transformed portion of the chain (martensite band) remains ÿnite:

l 0 = 1 1 + e : (63)
When nonlocal interactions are absent ( = -∞), the nucleus contains only one spring and l 0 = 0.

To illustrate the e ect of the outside boundary layers on the nucleation phenomenon we replace the special boundary conditions (6), which suppress boundary layers, by the "zero-moment" conditions (4) with B = 0. In this case the zero-interface solution is no longer trivial because the strain decreases exponentially near the boundaries (see Fig. 12b).

Suppose that the chain is originally in phase I. When the critical value d = d + (0) is reached, the strains in the middle of the chain pass the threshold w = w c and the one-phase solution becomes unstable (point A in Fig. 12). This leads to the formation of either two symmetric interfaces in the center (point C in Fig. 12) or of a single interface near one of the boundaries (point B in Fig. 12). Our computations show that the single-interface conÿguration B has a lower energy than the two-interface conÿguration C. However, during transition from A to C only two springs (7th and 9th) transform into phase II, while transition from A to B requires transformation of at least three springs which are also initially farther below the threshold. As a result, transition from A to B may encounter a higher barrier than transition from A to C, and then the nucleation will take place in the interior of the chain. While this possibility 0.6 0.7 0.8 0.9 can be investigated rigorously, for the subject of this paper it is enough to mention that the nucleation peak survives in both cases (see Fig. 12a).

Veriÿcation of the model

The explicit formulae from the previous section can be used to obtain bounds on the value of the nonlocality measure . The expression for from ( 10) and ( 60) has the form

= E 4 1 -1 - 2 E -2 ; ( 64
)
where E = K + 4 is the macroscopic elastic modulus of the homogeneous chain and is the macroscopic transformation strain. The magnitude of the stress drop at the peak and the transformation strain are available from the experimental data of [START_REF] Shaw | Thermomechanical aspects of NiTi[END_REF] on NiTi wires. For instance, in the experiment conducted at 70 • and the loading rate 4 × 10 -5 s -1 , the measurements gave = 0:039 GPa, and = 3:97%. The Young's moduli of the two participating phases are di erent, E = 56:7 and 27:5 GPa. By using two separate values we estimate to be between -1:1 and -1 GPa, which implies that is in the range -14:7 ¡ ¡ -7:3. Similar comparison of experiment and theory for CuAlNi yields = -33:4 [START_REF] Truskinovsky | Peierls-Nabarro landscape for martensitic phase transitions[END_REF].

To obtain an independent estimate of we can assume that the interactions between particles are governed by the Lennard-Jones potential. In this case we have 1 (r= -1) = 2 2 (r=2 -1) = U (r), where 1 (w) and 2 (w) are the energies of the NN and NNN springs, respectively, and U (r) has the form

U (r) = K 72 r 12 -2 r 6 : (65)
The coe cients in (65) are chosen to ensure that the elastic modulus at equilibrium r = equals K. Linearizing around the unstretched homogeneous state with the spacings r = and r = 2 , we obtain [START_REF] Charlotte | Linear elastic chain with a hyper-pre-stress[END_REF] = U ( ) 4U ( 2) :

This yields = -56:5. Taking E = 56:7 GPa as the value of the Young's modulus we obtain = -0:255 GPa which despite the rather rigid form of the potential ( 65) is within a reasonable range from the values obtained above. Finally, we can use the above estimates of to predict the initial size of the martensite band L 0 . By taking -14:7 ¡ ¡ -7:3 we obtain 0:001 ¡ L 0 =L ¡ 0:005, where L is the size of the specimen. While we could not ÿnd direct experimental measurements for L 0 in the literature, this parameter is usually estimated to be of the order of the specimen's diameter D (e.g. [START_REF] Sun | An inclusion theory for the propagation of martensite band in NiTi shape memory alloy wires under tension[END_REF]. This would be in agreement with our estimate since for thin NiTi wires used in experiment D=L∼0:0025 [START_REF] Leo | Transient heat transfer e ects on the pseudoelastic behavior of shape-memory wires[END_REF].

Strain-gradient approximation

Now we can compare the exact results obtained in the discrete model with the predictions of the quasi-continuum strain gradient approximation which neglects discreteness but retains the original internal length scale. To obtain the approximate model we temporarily reintroduce dimensional variables and perform the Taylor expansion in the small parameter = L=N . By preserving the ÿrst nonlocal term in the energy functional we obtain

= L 0 (w) + 1 2 A 2 (w ) 2 d x; (66)
where we omitted nonessential null Lagrangian contributions and introduced

A = - 1 12 (K + 16 ): (67)
The functional (66) must be minimized subject to the constraint

L 0 w(x) dx = d (68)
and the clamping boundary conditions

w (0) = w (L) = 0 (69)
which represent the continuum analog of (6). The Euler-Lagrange equations can be written in the form (w)-A 2 w =F (analog of (3)), where (w)=Ew-KaÂ(w-w c ).

Recall that E = K + 4 is the macroscopic modulus.

The trivial solution of the Euler-Lagrange equations w(x) = d=L is stable whenever the second variation of energy where n ¿ 1 is an integer. The solution exists for

E + n L 2 A = 0: (72)
If A = 0, the stability condition obviously reduces to E ¿ 0. When A ¿ 0, we obtain E + ( n =L) 2 A ¿ E + ( =L) 2 A and the trivial state is stable if and only if

E + L 2 A ¿ 0: (73)
The instability takes place through the growth of a long-wave mode v = cos( x=L).

When A ¡ 0, the trivial solution is always unstable because one can always ÿnd large enough n at which E +( n =L) 2 A ¡ 0; the corresponding instability is of the short-wave type.

To compare the stability limits in the strain gradient model with the ones obtained in the discrete model, we shall ÿrst rewrite them in terms of dimensionless , and replace L= by N . In the case when A ¿ 0, we obtain

-4 ¡ ¡ - 1 -2 =3N 2 1 -2 =12N 2 ( ¡ 0): (74)
Recall that in the discrete case the corresponding interval is -∞ ¡ ¡ -cos 2 =2N ( ¡ 0). It is easy to see that for large N both models predict the same upper boundary = -1 (E = 0), which means that the strain-gradient model captures the onset of macroscopic instability in the discrete model rather well.

As we showed above, when A ¡ 0 the stability range does not exist due to the short-wave instability. However, in view of the underlying discrete structure, a shortwave instability is unphysical if its wave length is less than the length scale of the lattice . If we assume that the mode number cannot exceed n = N , we obtain the following stability condition:

E + 2 A ¿ 0; (75)
which in terms of can be written as a combination of two intervals,

¡ -4 ( ¡ 0) (76)
and While the former stability range obviously complements (74), the latter has its own analog in the discrete model, -sin 2 =2N ¡ ¡ ∞ ( ¿ 0). Observe that the lower limit in the strain-gradient model is signiÿcantly higher than in the discrete model suggesting that even with the restriction on the minimum wave length, the strain-gradient approximation grossly exaggerates the instability domain on the side of the short-wave instability. The general trend, however, is predicted correctly: as in the discrete model, the short-wave instability in the gradient model develops before the strain enters the macroscopic spinodal region (see also [START_REF] Bardenhagen | Derivation of higher order gradient continuum theories in 2,3-D non-linear elasticity from periodic lattice models[END_REF]. Fig. 13 summarizes the stability analysis in both the K -and -N -1 planes. It should be compared with the analogous stability diagrams for the discrete model (Fig. 2). Now consider the nontrivial solutions of the variational problem (66), (68). As shown in [START_REF] Carr | Structured phase transitions on a ÿnite interval[END_REF], the strain-gradient model allows for local minima with at most one interface. For the piecewise linear model all these metastable conÿgurations can be obtained analytically [START_REF] Truskinovsky | Ericksen's bar revisited: Energy wiggles[END_REF]. A representative picture is shown in Fig. 14, with the parameters chosen to match those for the discrete model with N = 20. We observe that while the strain-gradient approximation captures the structure of the absolute minimizers of the discrete problem, it fails to reproduce the rich structure of the metastable equilibria which in the quasi-continuum model all collapse into a single branch. This is the reason why the strain-gradient approximation does not generate a realistic hysteresis. Also, while all multiple-interface equilibria of the discrete model are metastable, conÿgurations with more than one interface are absolutely unstable in the quasi-continuum approximation. This di erence is due to the fact that in the discrete model one cannot vary the volume fraction of phases continuously, and thus certain perturbations that make multiple-interface solutions unstable in the continuum case are impossible in the discrete model. Finally, we remark that the estimate -14:7 ¡ ¡ -7:3 obtained in the previous section from the experimental data implies that A ¡ 0, contrary to the typical assumption in the strain-gradient models that the coe cient A is positive.

- 1 -2 =3 1 -2 =12 ¡ ¡ ∞ ( ¿ 0): (77

Conclusions

Our study clariÿes why the behavior of the local and nonlocal systems with bi-stable elements is di erent. In the case of nonlocal interactions the elements around the phase-boundary are "pre-conditioned" due to the presence of the boundary layers. This results in smaller energy barriers for the propagation of the phase boundary compared to the local model and gives rise to the ÿnite mismatch between the spinodal stress and the Peierls stress. Another di erence is that in the nonlocal models the nucleation event typically involves transformation of more than one element. This produces nucleation peak which has been long known from experiment. We showed that the nucleation peak persists in the continuum limit where it manifests itself through an instantaneous formation of a ÿnite band scaled with the length of the sample. Although the nucleation peak is captured by the strain-gradient approximation, the resulting shape of the hysteresis loop is unrealistic. This suggests that the nonlocal coupling in the actual physical system is weaker than what is required for the gradient model to be applicable.

Fig. 1 .

 1 Fig. 1. A ÿnite chain of particles with nearest-neighbor (NN) and next-to-nearest-neighbor (NNN) interactions.

FigFig. 3 .

 3 Fig. 2. (a) Domain of stability for the trivial solution in K -plane: solid lines-discrete model, N = 3; dashed lines-continuum limit, N → ∞. (b) The dependence of the stability limits on the size of the system. In both ÿgures stability domain is in gray.

  and w i+1 ¿ w c (I to II phase switch) 0 if sign(w i -w c ) = sign(w i+1 -w c ) (no phase switch) -1 if w i ¿ w c and w i+1 ¡ w c (II to I phase switch):

Fig. 4 .

 4 Fig.4. The overall energy-strain and force-strain relations along the trivial (i = 0, N ) and the single-interface (1 ¡ i ¡ N ) metastable solutions. Here N = 10, = -2, wc = a = 1.

Fig. 5 .

 5 Fig. 5. (a) Energy-strain relation for the homogeneous solution (thin lines) and the two-interface solutions with various locations of the interfaces (thick lines). (b) The blow-up of (a) around conÿgurations with p = (0 -1 0 1 0 0) and p = (-1 0 1 0 0 0). Parameters: = -2, wc = a = 1, N = 6.

Fig. 6 .

 6 Fig.6. Energy-strain relation for the energy minimizers in the families of solutions with up to four interfaces.Here N = 6, = -2, wc = a = 1.

Fig. 8 .

 8 Fig. 7. (a) Energies of the single-and zero-interface branches of equilibria available at d = 0:5. (b) Peierls-Nabarro energy landscape along the path w k ( ) connecting metastable equilibria at d = 0:5. Parameters: N = 10, = -2, wc = a = 1.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Branch-switching sequence for the path of maximal delay, H = 0, and for the path with a ÿnite critical energy barrier H = 0:003. Parameters: = -2, N = 10, wc = a = 1.

Fig. 11 .

 11 Fig. 11. The maximum hysteresis loop in the continuum limit for the nonlocal (NNN) model with -∞ ¡ ¡ -1.

Fig. 12 .

 12 Fig. 12. Nucleation behavior of the system with external boundary layers: (a) single-interface (solid lines) and two-interface (dashed lines) metastable states; (b) the corresponding strain proÿles. The numbers in brackets indicate locations of phase the boundaries for conÿgurations with phase I on the left end. Parameters: N =15, = -2, wc = a = 1.

  problem determining the stability boundaries reduces to solving linear equation Ev -2 Av = 0. The general nontrivial solution of this equation compatible with (71) takes the form v(x) = cos nx L ;

  Fig. 13. (a) Domain of stability for the trivial solution in K -plane in the strain-gradient model with L= = N = 3 and in the continuum limit (N → ∞, dashed lines). (b) The dependence of the stability boundaries on the size of the system. In both ÿgures stability domain is in gray.

Fig. 14 .

 14 Fig. 14. The rescaled energy-strain and force-strain relations for the absolute energy minimizers (thick line) and the metastable single-interface solutions (thin lines) for (a) the discrete chain with N = 20 and (b) the strain-gradient approximation with = 1 20 . Other parameters: = -2, wc = a = 1. Unstable single-interface solutions are shown by dashed lines.

Acknowledgements

This work was supported by the NSF grants DMS-0102841 (L.T.) and DMS-0137634 (A.V.).


Brice Lecampion
email: brice.lecampion@csiro.au

Andrei Constantinescu
email: andrei.constantinescu@lms.polytechnique.fr

Laurent Malinski

Laurent Malinsky
email: malinsky@lms.polytechnique.fr

Identification of Poroelastic Constants of "Tight" Rocks from Laboratory Tests

Keywords: Poroelasticity, Permeability, Experimentation, Rocks, Laboratory tests, Hydraulic conductivity

come

Introduction

This paper presents a method for the identification of porolastic constants of rocks with very low permeability, also denoted as "tight rocks." As the concept of very low permeability often depends on the intended applications, we specify here that "tight" rocks refer to materials having an intrinsic permeability k under 10 -20 m 2 . Typical applications are hydro-mechanical responses of rock masses for geological modeling, design of underground waste repositories, or others rock engineering projects involving large time scales.

These extreme values of permeability reach actually the physical possibilities of direct measurements. In these ranges, it is not possible to accurately measure the flow output of the sample or to control the inside pore pressure. It is known ͑Hart and Wang 2001͒, that classical steady state tests are not appropriate for this type of rocks. Therefore, several experimental methods and interpretation ͑Escoffier et al. 2001;[START_REF] Hart | A single test method for determination of poroelastic constants and flow parameters in rocks with low hydraulic conductivities[END_REF] of the corresponding results have been developed to solve this problem.

Our main interest in this work lies, on the one hand, in interpreting the hydromechanical coupling during drained isotropic confinement and pulse tests performed on cylindrical specimens and, on the other hand, in determining the poroelastic parameters from these experiments.

Both tests can be performed on rocks with extremely low hydraulic conductivity and are of transient nature spanning several days. As presented in the next sections, the measurements of radial and/or axial displacements during a drained isotropic confinement test permit us to identify the drained and undrained bulk moduli as well as a diffusivity coefficient D lumping the hydraulic conductivity and the poroelastic parameters. The decrease of the reservoir pressure recorded during a pulse test is governed by two constitutive parameters: a lumped dimensionless coefficient of the poroelastic parameters and the same diffusivity coefficient D. A combined interpretation of these two tests should permit us to completely identify the poroelastic constants.

The main difficulty in the interpretation is the lack of simple closed-form solutions for such experiments ͑Adachi and Detournay 1997͒. Even in the framework of linear poroelasticity ͑Biot 1941; [START_REF] Detournay | Fundamentals of poroelasticity, in comprehensive rocks engineering[END_REF][START_REF] Coussy | Poromechanics[END_REF], no such solutions are known. In order to overcome this type of difficulty, it is common to interpret such transient tests with tabulated charts ͑Neuzil et al. 1981;[START_REF] Escoffier | Mesure de la perméabilité et du coefficient de Biot dans les argilites de l'est ͑in French͒[END_REF], sometimes neglecting the poroelastic coupling. As both tests are fully two dimensional, we shall first show that one-dimensional closed form solutions can be derived and applied providing the. specimen is slender enough. The validity of the one-dimensional solution for the twodimensional field is checked using fully coupled finite element simulations.

A complete identification method for the pulse and drained isotropic confining test is proposed. The parameter identification is here formulated as an inverse problem. The poroelastic constants are obtained by minimization of a cost functional measuring the discrepancy between computed predictions and measurements.

Finally, experiments performed on deep argillaceous rocks are analyzed using the proposed method and several details related to the modeling and the identification procedure are discussed.

Tests Procedure

Both the pulsed and the drained isotropic confining test are performed on a cylindrical core ͓͑r , z͒ ͓0,R͔ ϫ ͓0,L͔ with R = 1.8 cm, L = 7.2 cm͔ under axisymmetric conditions.

The specimen is set between two porous end caps in a triaxial cell. The cell allows us to apply a confining pressure and the end caps permit us to impose the pore pressure at the top and bottom of the sample. The complete experimental setup is based on an extremely rigid servo-controlled frame. Precise details of the experimental setup can be found in Malinsky et al. ͑2002͒. However, let us simply highlight here several important aspects of the test facility: 1. The cylindrical surface of the specimen is lined with an impermeable jacket. The validity of this system for extremely low hydraulic conductivities rocks has been assessed by tests performed on steel cores with a rough cylindrical surface. 2. The compressibility of the pore pressure injection system has been reduced to a minimum. 3. The effect of thermal variations are corrected dynamically using a reference circuit. 4. The injected water is a reconstructed formation water. 5. Radial and axial displacements are measured via extensometers, while pressure-volume controllers are placed in the top and bottom hydraulic circuit. Prior to the tests, the sample is resaturated under isotropic confinement for several weeks. After this preparation, the sample is in an homogeneous and isotropic hydromechanical state which can be considered as a reference initial state with regards to the linear theory of poroelasticity.

The drained isotropic confinement test consists of the sudden application at t = 0 of an isotropic confining normal stress on all the surfaces of the specimen. The surface tractions are given by The recorded axial and radial displacements measure the evolution of the specimen from the undrained to the drained values.

The pulse test consists of the instantaneous injection of a water volume V 0 in a deformable reservoir of stiffness C res connected to the porous sample ͑see Fig. 2͒. This injection produces an instantaneous pressure increment p 0 in the reservoir. After a transient phase, the system reservoir specimen reaches an equilibrium state where the specimen pore pressure and the reservoir pressure are equal. The decrease of the reservoir pressure from p 0 toward its equilibrium value depends on both the poroelastic coefficients and the permeability of the specimen, as well as on the reservoir stiffness.

This where q = fluid flux vector flowing in the specimen. The confining stress is kept constant during the test and is equal to the one of the reference state. However, as the pore pressure cannot be applied on a porous material without applying a mechanical pressure, the top and bottom sufaces ͑z =0,L͒ of the samples are also subject to a compression ͑t͒ • n =p res ͑t͒n If this last boundary condition is ignored, then the Skempton effect will be equally ignored. As a consequence the estimation of the poroelastic parameters will be incorrect.

Modeling and Parameter Identification

Closed Form Solutions

The poroelastic fields describing the tests are fully twodimensional. As an illustration of this fact, we have displayed spatial distribution of the shear stress component during a drained isotropic confinement test in Fig. 3.

As explained by Adachi and Detournay ͑1997͒, the poroelastic coupling does not allow to obtain a closed form solution for these particular configurations. This stems essentially from the hydromechanical boundary conditions on the cylindrical surface: no hydraulic flux and an applied radial stress. In spite of the apparent simplicity of these boundary conditions, the theory of linear poroelasticity yields to nonhomogeneous hydromechanical fields in the core during the transient phase. It has also been shown ͑Adachi and Detournay 1997͒, that for the limiting case of a very slender specimen ͑m = R / L Ӷ 1͒, the continuity equation governing the pore pressure evolution in the sample, and the balance equation uncouple. In the limiting case, stress and pore pressure fields are homogeneous in the sample, which finally translates into a one-dimensional closed-form solution. The expressions of the closed form solutions for both tests are given in the Appendix.

We have first addressed the validity of the slender assumption in the case of our specimen, which has a given aspect ratio m = R / L = 0.25. A series of axisymmetric finite element simulations using the object oriented finite element code Cast3M ͑2005͒ have been performed for different aspect ratios in order to investigate the hypothesis of slender specimen ͑m = R / L Ӷ 1͒.

We have analyzed the relative error between the axisymmetric FEM and one-dimensional solution for different quantities measured during the experiments. As a typical example in Fig. 4, we have plotted for the case of the drained isotropic confinement test the relative error on the vertical displacements for a series of different aspect ratios. We remark that the error is only significant during the transient phase and decreases with the value of the aspect ratio m = R / L. Similar patterns have also been observed for the pulse test. We also note that computations for different values of the shear modulus did not show a major influence on this relative error.

As a consequence, we can conclude that for the aspect ratio of the specimen used in our experiments ͑m = 0.25͒, the coupled effect is negligible. The relative error ⑀ on the axial strain is at most of 1% and will therefore permit the usage of the one-dimensional closed-form solution in the sequel. This fact simplifies the computional effort for the identification process, saving both programming and computation time compared to a finite element solution.

Parameter Identification

The question of the parameter identification problem is: "How to invert the displacements and the reservoir pressure measured during the tests in order to obtain an estimation of some poroelastic parameters?" A direct interpretation of the solutions of both test shows that different parameters are related to one or the other solutions ͑see Table 1͒. The drained K, the undrained K u bulk moduli, and the diffusivity coefficient D ͓expression ͑3͔͒ can be identified from the recorded evolution of the axial ͓u z ͑z = L , t͔͒ and/or radial ͓u r ͑z = L /2,t͔͒ displacements during the drained isotropic confinement test ͓see Eq. ͑5͒ in the Appendix͔. From the pulse test, the dimensionless parameters ͓expression ͑6͔͒ and the diffusivity coefficient D can be identified from the recorded reservoir pressure decay p res ͑t͒ ͓see Eq. ͑5͒ in the Appendix͔.

A close inspection of the expressions for displacements and pressure show that they involve infinite series and Laplace transforms. As a consequence, we do not possess a simple expression relating directly the unknown parameters to the measurements. The proposed solution of the identification problem will be therefore based on the minimization of a least square cost functional measuring the distance between predictions and measurements. The minimization of this functional is performed using a Levenberg-Marquardt gradient based algorithm ͑see Gill et al.

for details͒

One of the important technical requirements of gradient based optimization algorithms is the accuracy of the gradient component. In the case where the direct problems are solved numerically, special numerical techniques such as the adjoint state method and the direct differentiation needs to be used in order to obtain precise estimation of the gradient ͑see Lecampion and Constantinescu 2005 for a derivation of the poroelastic case͒. However, here the knowledge of the closed form solutions permitted the computation of the exact gradients of the solution in the directions of the unknown parameters, the so-called sensitivities

du dK du dK u , ¯dp dM , dp d
The sensitivities of the axial displacement with respect to different parameters in the case of the isotropic confinement test are displayed in Fig. 5. We should note that the experimental response contains information on the undrained modulus K at early times, but the information fades later. At large times, the displacement depends, as expected, mainly on the drained modulus. The sensitivity on the diffusivity coefficient is maximum during the transient phase. Similar patterns can be observed for the pulse test and could be used for fine tuning of the identification algorithm. We recall that the estimation of unknowns parameters from measurements is an inverse problem and that it may not have a unique solution ͑Vogel 2002͒. However, for both tests considered here, algebraic manipulations of the closed form solutions permit us to show that uniqueness and stability are assured at least for the ideal mathematical problem ͑see Lecampion 2002 for details͒.

In the case of real data, where the poroelastic model may not be adequate, the identification can prove to be difficult and unstable. The validity of the results can be simply assessed by running the minimization algorithm with several different starting points for the initial guess of the parameters.

Results and Discussion

The drained isotropic confinement and the pulse test have been performed on several argillaceous rock samples. The tested argilite belongs to the Callovo-Oxfordian formation of Meuse Haute-Marne ͑Bure, France͒ ͑Malinsky et al. 2002͒. The samples are taken from cores at a depth of about 500 m, and its properties of this rock are extensively described in ANDRA, Direction Scientifique ͑1999͒. In this section, we present some of the identification results for the poroelastic constants obtained using the tests and the identification procedure described previously.

Drained Isotropic Confinement

Fig. 6 and7, display the measured and fitted responses of an argilite sample submitted to a drained isotropic confinement. As a first remark, one can note the high level of anisotropy between the axial and radial deformation. The identification performed using only the axial displacement of the isotropic drained consolidation on an argilite sample gives a perfect fit ͑Fig. 6͒. For different initial guesses the same set of identified parameters is obtained K = 1.04 GPa; K u = 9.82 GPa; D = 4.02 ϫ 10 -9 m 2 s -1

The identification performed using only the radial displacement recorded at the middle of the specimen yields the following values of the parameters: K = 2.12 GPa; K u = 16.9 GPa; D = 1.74 ϫ 10 -9 m 2 s -1

Comparing the experimental and computed responses one can observe a certain misfit ͑see Fig. 7͒. A close inspection of this figure shows on the one hand that the measured radial displacement has not reached its asymptotic value after 15 days, whether the computed one already tends to its asymptotic value. On the other hand, at early times the measured initial delay due to the diffusion process is rather small when compared to the computed one. In this case, using the knowledge of the uniqueness result for the identification problem one can state that the mismatch is not an artifact of the identification. The linear poroelastic model is probably not suited to represent the material behavior. This has later been confirmed by swelling observed on the same sample.

͓ML -1 T -2 ͔ Drained isotropic confinement D = k f MK K n diffusivity ͓L 2 T -1 ͔ Drained isotropic confinement/pulse test = MK C res K u R 2 L ͑ 1-B 3 ͒ ͓•͔
The anisotropy is clearly straightforward when looking at the parameters estimated with either axial or radial displacement. This anisotropy can be explained by the high level of damage of the sample, which most probably occurred during sample preparation.

Knowing K, K u , and D, it is possible to compute the intrinsic permeability assuming reasonable value for Biot's modulus M. Table 2 presents the results for different values of M. An intrinsic permeability on the order of magnitude of 1 ϫ 10 -21 m 2 matches the expected values for this argilite.

Pulse Test

The identification method has also been applied to pulse tests. The experimental as well as the identified response are displayed in Fig. 8. Table 3 presents the identified parameters as well as the value of the cost functional for different initial guesses. The best fit has been obtained for the following set of values: = 0.653; D = 1.65 ϫ 10 -9 m 2 s -1 However, the match between the measured and computed responses is far from perfect. As before, taking into account the stated uniqueness of the identification procedure, one can conclude that this relative misfit is attributed to the poroelastic model, which is not able to capture the complex behavior of this rock sample.

At this point it is important to note that the exact value of the reservoir stiffness C res intervenes directly in the identification process and has to be known with precision in order to extract accurate values of the poroelastic parameters. The reservoir stiffness can be estimated by additional tests on the reservoir. Unfortunately, in our case, this estimation was far from perfect.

In Table 5, we have displayed the value of the coefficient MK / K u and that of the intrinsic permeability k for a range of acceptable values of the reservoir stiffness C res . We should note that all the obtained values are in a standard range for an argilite.

Discussion

As a first conclusion, we note that all the identified values of the parameters lies in the physical range known from the literature ͑see Table 4͒ Despite the reported misfit between model predictions and measurements we can conclude that the first order response of the One should not be mislead by the obtained experimental results. The technical difficulties of performing tests on such argillaceous materials are enormous and the values of the parameters should be taken with precaution as preliminary estimations. Several experimental and modeling issues regarding the chemically reactive and swelling behavior of this argilite should also be investigated in more details.

An important advantage has been obtained by performing the two tests, drained isotropic confinement and pulse test, together. By combining the identified coefficients from the two tests ͑K , K u , D, and ͒, it is possible to estimate all the poroelastic parameters and the permeability with only the exception of the shear modulus. Moreover, as the diffusivity coefficient D is estimated twice in two different ways, it ensures the robustness of the estimation. Complete results of this procedure could not be displayed here due to confidentiality agreements.

Finally, we note that the closed form solution used for the pulse test is similar to the one proposed by Hsieh et al. ͑1981͒ for the uncoupled case. The difference is that we explicitly take into account the presence of a mechanical pressure on the reservoir surface. The presence of this boundary condition largely affects the identified parameters in the coupled case.

Conclusions

In this paper, we have discussed the identification of several poroelastic parameters from drained isotropic confinement and pulse tests. We have shown that for specimens with an aspect ratio R / L Ӷ 0.25, we can use the one-dimensional closed form solutions obtained for a slender specimen. This permits us to speed up the identification procedure. The proposed identification method has been shown to provide robust and efficient solutions to this inverse problem. Note that the induced mechanical loading during a pulse test greatly affects the estimated value of the poroelastic parameter as it induces a Skempton's effect.

The identified values of the constitutive parameters for the argilite of Meuse Haute-Marne are within the expected range. The partial mismatches between the poroelastic model and the experiment show that probably other physical phenomena ͑swelling, chemical reaction . . .͒ should also be taken into account.

We can conclude that the combination of the drained isotropic confinement and the pulse test permits a complete identification of the poroelastic parameters and hydraulic conductivity with the exception of the shear modulus. The methodology presented here can be extended to other porous materials such as tight sandstone or granite. From the linear poroelastic constitutive equations in cylindrical coordinates ͑Biot 1941; [START_REF] Detournay | Fundamentals of poroelasticity, in comprehensive rocks engineering[END_REF][START_REF] Wang | Theory of linear poroelasticity with applications to geomechanics and hydrogeology[END_REF], one obtains the overall axial, and the radial displacement at the middle of the specimen. These displacements are recorded during an experimental test

u z ͑z * = 1,t * ͒ L = 3K ͫ 1 -K u -K K u I͑1,t * ͒ ͬ u r ͑z * = 1/2,t * ͒ R = 3K ͫ 1 -K u -K K u K͑1/2,t * ͒ ͬ
where I͑l , t * ͒ = ͐ 0 1 K͑ , t * ͒d is given by the following series:

I͑1,t * ͒ = 2 2 ͚ n=1 ϱ exp͑-n 2 2 t * ͒ͩ 1 n - cos n n ͪ 2
Series ͑4͒ and ͑5͒ are precisely estimated even for early times ͑t * Ϸ 0.01͒ when keeping the first ten terms of the series.

Pulse Test

This solution is similar to the classical one given by Hsieh et al.

͑1981͒ with the difference that we take into account the effect of the mechanical loading on the specimen due to the reservoir pressure.

The initial and boundary conditions ͑1͒ give the following additional equations to the pressure equation ͑2͒: Another dimensionless parameter related to the constitutive parameters and the reservoir compressibility is introduced

␥ = MK C r K n R 2 L
The Laplace transform of the solution ⌸ ¯is given in the Laplace domain by

⌸ ¯͑z * ,s͒ = cosh͓ ͱ s͑1/2 -z * ͔͒ s coshͩ ͱ s 2 ͪ+ ͱ s sinhͩ ͱ s 2 ͪ ͑5͒ with = ␥ͩ1 - B 3 ͪ= KM K u C res R 2 L ͫ1-ͱ M͑K u -K͒ 3K u ͬ ͑6͒
Note that if we do not take into account the mechanical pressure, the solution has a similar form but with = ␥: the Skempton's effect is ignored. No simple analytical expression of the inverse Laplace transform exists for expression ͑5͒. The inversion can be performed numerically using Stehfest algorithm, for example ͑Cheng et al. 1994͒.

•

  n = n H͑t͒ where ϭstress tensor; and n=outside unit normal and H͑t͒ denotes the Heavyside step function. The cylindrical surface ͑r = R͒ of the sample remains undrained, i.e., the fluid flux through the surface is zero q • n = 0 while the top and bottom surfaces ͑z =0,L͒ are drained ͑see Fig. 1͒

  p = pore pressure.

  transient method has been proposed first by Brace et al. ͑1968͒ for rocks with a negligible storage coefficient. The method has been extended afterward to other types of rocks by Hsieh et al. ͑1981͒ and is actually largely used ͑Zeynaldi-Andabily and Rahman 1995; Escoffier et al. 2001͒. However, in the case of a poroelastic coupling, special care has to be taken with respect to the real boundary condition. The cylindrical surface ͑r = R͒ of the specimen is kept undrained during the test. The pore pressure on the boundary between the reservoir and the specimen, ⌫ res , at the top and bottom surfaces, equals the fluid pressure in the reservoir p res ͑t͒ at all time instants. The evolution of the reservoir pressure

Fig. 1 .

 1 Fig. 1. Boundary conditions ͑left͒ and evolution of radial and axial displacements ͑right͒ for drained isotropic confinement test

Fig. 3 .

 3 Fig. 3. Isovalue contours of shear component of stress tensor during transient phase of isotropic confining test. Finite-element model Cast3M, axisymmetry, poroelastic coupling.

Fig. 4 .

 4 Fig. 4. Relative error on axial strain between closed form and coupled finite-element solution for different aspect ratios m = R / L

Fig. 5 .

 5 Fig. 5. Dimensionless total axial displacement u z / L and associated sensitivity coefficient on K , K u , and D during drained isotropic confinement test

  . Previous laboratory experiments ͑Vincké et al. 1998; Escoffier et al. 2001͒ suggest a value around 10 -21 m 2 for the intrinsic permeability of these rocks. The value that we can deduce ͑Tables 3 and 5͒ tends to be slightly lower: 10 -22 m 2 . The values of the elastic undrained moduli found are around known values. The drained values obtained here are lower than the one reported previously ͑Escofier et al. 2001͒. This may be partly due to the high level of damage to the sample observed after the test.

Table 2 .Fig. 6 .Fig. 7 .

 267 Fig. 6. Experimental and identified response ͑axial strain͒ to isotropic confinement test performed on argilite of Meuse Haute-Marne

  The boundary conditions are homogeneous⌸͑z * = L,t͒ = ⌸͑z * = 0,t͒ = 0 ⌸͑z,t * = 0͒ = 0The dimensionless pore pressure is easily obtained from the solution of the diffusion equation in a bounded one-dimensional domain with an instantaneous volumetric source term at t * =0 ͑Carslaw and Jaeger 1959͒ ⌸͑z * ,t * ͒ = -BK͑z * ,t * ͒ with K͑z , t͒ K͑z * ,t * ͒ = 2

Table 1 .

 1 Identifiable Coefficients and Corresponding Tests

	Coefficient	Dimensions	Test
	K , K u drained and
	undrained bulk moduli

Acknowledgments

The authors would like to thank Professor E. Detournay and Dr. J. I. Adachi for the benefit of useful discussions. This research was funded by the French Agency for Radioactive Waste Management ͑ANDRA͒ through a Ph.D. Grant for B.L. The experimental campaign was run for ANDRA under Contract No. ANDRA 017079-Contract-cadre AQL. The experimental tests were performed at G3S, Ecole Poly Y technique, France by Serge Chanchole and François Coste.

Appendix. Closed-Form Solutions for Slender Case

For the slender case ͑R / L Ӷ 1͒, the problems are one dimensional and the continuity and balance equations uncouple. The diffusion equation governing the pore pressure now reduces to

where tr͑͒ = kk = volumetric part of the stress tensor; and B͓•͔ = Skempton dimensionless coefficient ͑between brackets ͓•͔ we indicate the physical dimensions of the parameters͒. The diffusivity D͓L 2 T -1 ͔ is given by ͑Detournay and Cheng 1993; Wang 2000͒

where


M Meyer

Yu L Raikher

O Sandre

A Bée

V Cabuil

V Dupuis

P Licinio

R Perzynski

Magneto-orientational properties of ionically stabilized aqueous dispersions of Ni(OH) 2 nanoplatelets

Keywords: . 82.70.Dd Colloids -83.80.Hj Suspensions, dispersions, pastes, slurries, colloids -75.50.Tt Fineparticle systems; nanocrystalline materials -75.20.Ck Nonmetals

Magnetic and orientational behavior of nickel hydroxide nanoplatelets ionically stabilized in a liquid matrix is studied. Under an applied field the platelets orient their faces normal to its direction. For characterization of the individual behavior of dispersed and non-interacting particles three techniques are used: SAXS, SQUID and magneto-optics. Analysis reveals that nickel hydroxide in a platelet phase is paramagnetic with a pronounced anisotropy of the intrinsic susceptibility, the major component of which (in the direction normal to platelet face) exceeds the minor one by about 25%.

Due to their abilities to self-organize [1,2], concentrated dispersions of nanorods or nanoplatelets raise large scientific and applicational interest. Self-organization of these anisometric (and often intrinsically anisotropic) nano-objects, which are usually sterically stabilized [3][4][5], results from competition between the orientational entropy favoring the isotropic state, and the excludedvolume entropy favoring a liquid-crystalline order [6]. Aqueous dispersions of clay [7,8] and hexagonal gibbsite [9] are among the very few examples of electrostatically stabilized dispersions of nanoplatelets reported in the literature. In such systems, the effect of shape anisotropy (anisometry) of individual objects on the interparticle interaction is strongly modified by Debye electric screening [10]. Tuning of the latter mechanism allows one to change the interparticle repulsion from strongly anisotropic to roughly isotropic. Accordingly, typical phase diagrams evidence a complicated interplay between nematic and dense gel phases [7][8][9].

The system described in this paper belongs to this class of ionic colloids. Namely, we study an aqueous dispersion of Ni(OH) 2 (nickel hydroxide) nanoplatelets, which in a dilute regime displays very strong field-induced orientational properties. Ni(OH) 2 platelets are getting extensively involved as parent systems for manufacturing new a e-mail: regine.perzynski@upmc.fr types of nanostructured high-discharge capacity battery electrodes [11][12][13]. An insight to the state of the particles, a deeper understanding of their properties and a possibility to control the internal structure of their dispersions thus seem to be quite timely.

In some aspect, the effects encountered are similar to those reported for lyotropic colloids of mineral particles: goethite nanorods [14], Fe-montmorillonite [15] or gibbsite [16,17]. However, the orientational ordering that takes place in nickel hydroxide nanodispersions is significantly different: under a magnetic field the platelets orient their planes perpendicular to the direction of the latter as if they were striving to build up a diskotic nematic phase.

In what follows, we begin with a brief description of the colloid preparation and the analysis of the morphology and size of the suspended platelets. Then we focus on the magneto-orientational properties of the system which are characterized at room temperature by small-angle X-ray scattering (SAXS) and by magneto-optical birefringence methods. As the particle paramagnetic susceptibility is an issue of our interest, the dilute dispersions of Ni(OH) 2 platelets are explored by magnetic SQUID measurements in the temperature range 5-300 K. Finally, the experimental results obtained by different techniques are compared and discussed to identify the individual response (averaged over an ensemble) of Ni(OH) 2 platelets to an applied field. Nickel hydroxide particles are chemically synthesized by adding ammonium hydroxide NH 4 OH to a Ni(NO 3 ) 2 solution and heating at 80 • C under vigorous stirring. This synthesis is extensively described elsewhere [18]. Electron microscopy and X-ray diffraction are used to investigate the structure, size and morphology of the obtained objects [18]. Figures 1 and2 show typical electron transmission micrographs and diffraction patterns of β-Ni(OH) 2 platelets, evidencing lamellar structure of spacing ≈ 4.7 Å (close to the bulk value 4.6 Å) and in-plane hexagonal packing. Powder X-ray diffraction (PXRD) of disordered samples, see [18], shows diffraction peaks characteristic of β-Ni(OH) 2 structure, enlarged because of the small size of the crystals. The diffraction patterns contain several peaks [hkl] of various widths, which be indexed. Assuming that the distortion can be neglected, we use Scherrer formalism to deduce the diameter D = 53 ± 3 nm of the platelets from peak [100] (±3 nm being the uncertainty of the measurement) and their thickness b = 8 ± 1 nm from peak [001]. These determinations are fairly compatible with the typical log-normal diameter distribution D 0 = exp( ln D ) = 55 nm and width σ = 0.2 deduced from TEM micrographs, and with Figure 1.

To properly disperse the platelets in water, a sufficiently high surface density of charge is necessary [18,19]. For this purpose, coating with citrate ligands attained by adsorption equilibrium is performed. The volume fraction Φ of the platelets in a stable aqueous dispersion can reach up to 0.15 [19] at 8 < pH < 12, while the concentrations of unadsorbed citrate [cit] free typically range between 5 • 10 -3 mol/L and 3 • 10 -2 mol/L. The phase diagram as a function of Φ and of the ionic strength of the solution can be found in [19]. In qualitative terms, the dispersions are flowing at low Φ and they are "soft solids" at larger concentrations, the boundary between these two states being ionic strength dependent. As we are interested here in the behavior of individual platelets (averaged in the dispersion over an ensemble of non-interacting particles), we consider very dilute systems with pH = 10.5, where the osmotic pressure is experimentally almost independent of the ionic strength. In this regime the surface charge on the platelets is of the order of -15 μC/cm 2 (meaning roughly one electron charge per nm 2 ) under the concentration of unadsorbed citrate [cit] free ≈ 5 • 10 -3 -3 • 10 -2 mol/L. Quasielastic light scattering experiments, performed on dilute dispersions at Φ = 6 • 10 -4 lead to a typical hydrodynamic diameter of 60 nm fully compatible with the powder analysis by TEM and PXRD. It shows that there is no particle aggregation in these dilute dispersions.

Small-angle X-ray scattering measurements under applied magnetic field are performed on D24 beam line in LURE (Laboratoire d'Utilisation du Rayonnement Electromagnétique, Orsay, France). The configuration of the apparatus is tuned to wavelength λ = 0.148 nm selected by a germanium monochromator with the sample-detector distance 2.76 m. This allows to cover the wave vector interval 6 • 10 -2 nm -1 < q < 1 nm -1 . Scattering patterns of the samples at 0.4 • 10 -2 < Φ < 2.9 • 10 -2 and [cit] free = 10 -2 mol/L are recorded on a 2D detector with a magnetic field ranging between 0 and 640 kA/m. The signal of the solvent cell is used to subtract the experimental background from the scattering patterns. Being isotropic in zero field, the scattering patterns become anisotropic in its presence. The zero-field measurement allows to determine the form factor of the dispersed objects (data not shown). Because of their polydispersity (σ = 0.2 from TEM micrographs), this form factor is not oscillatory at large q. We deduce the gyration radius R g = 20 ± 1 nm of the nanoplatelets from the slope -R 2 g /3 of a semi-logarithmic plot of the intensity as a function of q 2 in the range qR g ≤ 2. Their thickness b = 6.1 ± 0.3 nm is measured in the intermediate q regime π/R g < q < π/b from the slope -b 2 /12 of a Kratky-Porod plot ln(Iq 2 ) vs. q 2 , which (a) the scattering pattern elongated along the field direction (to better reveal the anisotropy, a greyscale with only 8 levels is used); α is the azimuth angle defined in the text; (b) the angular dependence of the scattered intensity I at a given scattering vector q = 0.29 nm -1 fitted with the aid of formula (1). Its experimental value at H = 0 is 680 a.u. confirms the platelet shape of the dispersed objects. We then calculate the nanoplatelet diameter D = 56 ± 3 nm using the expression R 2 g = D 2 /8 + b 2 /12 assuming axial symmetry of the dispersed objects. This value of D agrees fairly well with other our measurements. Figure 3a illustrates the anisotropy of the scattering pattern for a sample of Φ = 2.9 • 10 -2 under a field H = 640 kA/m. The pattern is elongated along the direction of the applied field h = H/H, meaning that the platelets orient their major dimension perpendicular to the field. The scattering intensity at a given q value varies along the azimuthal angle according to [17,20]

H ν ν ν H H S 2 = 1 S 2 = -½ S 2 = 0
I(q, H, α) = I(q, H = 0) 1 -(1 -3 cos 2 α) S 2 (H) , ( 1
)
where I(q, H = 0) denotes the isotropic value of I in the absence of field and

S 2 = 3 2 [ (νh) 2 -1 3 ]
is the orientation order parameter built on unit vector ν normal to the platelet face and defined as a statistical average over the particle ensemble; the azimuthal angle is α = arccos[(qh)/q]. The orientational structures corresponding to the limiting values of the order parameter S 2 are shown in Figure 4.

Figure 3b shows a fit with equation (1) of the data of Figure 3a recorded on the 2D detector along a full circle at q = 2.9 nm -1 under H = 640 kA/m. The value of the orientation order parameter rendered by this adjustment is S 2 = 0.096. At constant H, the parameter S 2 with the accuracy ±8% depends neither on q nor on Φ within the respective ranges investigated. Figure 5a presents the field dependence of S 2 ; the data scale as H 2 , confirming that we deal with a polarizable system. As seen from Figure 5a, the observed value of S 2 is always positive, in contrast to what was obtained with goethite nanorods in [10], where a nonmonotonic behavior was observed. Moreover, as S 2 (H) is found to be independent of both q and Φ, one concludes that the orienting entities indeed possess shape anisotropy and behave independently.

The individual field-induced behavior is confirmed by another experiment, which is also sensitive to the parameter S 2 , but which is based on a completely different technique and range of wavelengths. Namely, with a set-up detailed in [21], a magneto-optical measurement is performed. We find that under an applied magnetic field varying between 0 and 900 kA/m, dilute dispersions of Ni(OH) 2 platelets display an optical birefringence Δn which is negative. In the experimental range of H and Φ, the induced optical anisotropy Δn is proportional to the square of the applied field, see Figure 5b, and linear in the volume fraction Φ of the dispersed particles. Birefringence of such a type implies that under field a progressive orienting of the objects with a negative intrinsic optical anisotropy takes place. Given that, the attained optical anisotropy takes the form

Δn = -|δn o |Φ S 2 (H), (2)
where the intrinsic anisotropy δn o is a function of the particle shape and of its intrinsic refraction index. Substituting in equation ( 2) the obtained SAXS results, one obtains for the mean value of the intrinsic anisotropy the estimate δn o = -0.05 ± 0.005. We note that the field-induced orientation reported here is essentially different from the one obtained in [15,16], where the platelets positioned their flat faces parallel to the field. The explanation is that the mineral microdisks in [16] are diamagnetic; hence, due to the shape effect, their field-induced dipoles align with the major dimension of the object.

In principle, the experimentally established facts of S 2 (H) being positive and δn o negative, could be accounted for with the aid of two alternative models. One assumes that the particles have permanent magnetic moments directed along ν. The second one is based on the platelets being paramagnetic with the intrinsic susceptibility tensor possessing a pronounced anisotropy. However we have not found any indications of permanent magnetic moments neither in our measurements at room temperature (see just below) nor in the literature. According to [22][23][24][25], where magnetic properties of disperse Ni(OH) 2 were investigated quite thoroughly at 4-300 K, the particles lose global antiferromagnetism at about 20 K and any signs of metamagnetic behavior at about 100 K. This eliminates any hypothesis of existence of a spontaneous magnetic moment, and makes the anisotropic paramagnetic model with the maximum component χ of intrinsic susceptibility along the normal ν to the platelet the only plausible one.

Magnetization measurements performed with a SQUID (Superconducting Quantum Interference Device, Quantum Design, MPMS) confirm the paramagnetic nature of the response of dilute dispersions of Ni(OH) 2 to an applied field at room temperature. As follows from Figure 6, the magnetization of a colloid with volume fraction Φ = 1.4% is proportional to the applied field up to H = 4 • 10 3 kA/m. Using the data presented in this figure, we obtain for the susceptibility of the suspension in the direction of the field χ s = M/H = 4.2 • 10 -5 . The equilibrium susceptibility of a suspension of uniaxially anisotropic paramagnetic particles subjected to Brownian rotary motion depends on the orientation order parameter in the following way [26]:

χ s = Φ χ ⊥ + 2 3 S 2 + 1 3 Δχ , (3)
where χ ⊥ is the susceptibility per unit volume of the platelet in the direction perpendicular to ν, and Δχ = χχ ⊥ . Since the orientation order parameter S 2 grows with H, the susceptibility (3) is weakly field dependent. Assuming proportionality of the field-induced magnetization to the particle volume fraction, one obtains the paramagnetic susceptibility of the particle substance

χ = χ s /Φ = (3.0 ± 0.1) • 10 -3 ; (4)
here a possible correction due to shape anisotropy is negligible. The estimate given by equation ( 4) agrees fairly well with the value χ = 2.6 • 10 -3 found in magnetic measurements of [24] on dense Ni(OH) 2 gels at 300 K under much weaker field H = 280 kA/m. For a particle with an anisotropic paramagnetic susceptibility embedded in a liquid matrix, the orientation-dependent part of Zeeman energy is

U = -1 2 μ 0 ΔχV H 2 (νh) 2 , ( 5
)
where V is the particle volume. Performing statistical averaging over the particle orientation angle with the Boltzmann distribution W ∝ exp(-U/kT ) and expanding the result for ΔχV H 2 /kT < 1, one finds for the orientation order parameter

S 2 (H) = μ 0 ΔχV H 2 15kT . ( 6
)
Modeling our objects as identical round cylinders and using the particle size measurements, one finds V ∼ 2 × 10 -23 m 3 . Further, from S 2 (640 kA/m) = 0.1 with the aid of equation ( 6) one arrives at the estimate Δχ ≈ 6 × 10 -4 . This means that the intrinsic magnetic susceptibility tensor for the platelet phase of Ni(OH) 2 is of the form

χ ik = χ ⊥ δ ik + (χ -χ ⊥ ) ν i ν k , ( 7
)
with χ ⊥ ≈ 2.7•10 -3 and χ /χ ⊥ ≈ 1.25. Under thermalization in the absence of the field the net value of the initial susceptibility of a colloid should be χ s = 1 3 Φ(χ + 2χ ⊥ ), that is about 2.9 • 10 -3 Φ. Experimentally, with the given accuracy this value is hardly distinguishable from the one of equation ( 4). This is understandable as the degree of orientation in the suspension of paramagnetic particles is not high: ∼ 10%. However SAXS and birefringence experiments, being sensitive exclusively to S 2 , unambiguously prove its existence.

Crystallographic data confirms that in the platelets the axis ν, i.e., the normal to the particle face, is singled out. Indeed, it is well known, see [24,25] for example, that the platelet crystal is a layered structure with a hexagonal inlayer packing of Ni atoms. In the antiferromagnetic state at low temperatures, the c-axis becomes the easy axis of magnetic anisotropy. Although at room temperature there is no spontaneous spin ordering, it is reasonable to suppose that the spin-orbit interaction is still present, and the field-induced paramagnetic polarization is enhanced in the direction of the crystal field.

The temperature dependence of the field-cooled (FC) and zero-field-cooled (ZFC) susceptibilities in the temperature range 5-300 K in the field H = 4.1 kA/m are also investigated. The results are presented in the inset of Figure 6. The magnetic behavior at high temperatures is typically paramagnetic: FC and ZFC magnetizations coincide and behave as 1/T . At low temperatures, the FC magnetization continuously increases as the temperature is decreased which most probably indicates metamagnetic behavior. On the other hand, the ZFC magnetization presents a peak at a temperature 25 K, which matches closely the Néel point T N for the low-temperature antiferromagnetic phase which does not differ greatly for the platelet and bulk phases, see [22][23][24][25]. However, already in [24] some indications of appearance of particle spontaneous magnetic moments below T N were mentioned. This was at the time when Néel was preparing his concept [27] according to which fine grains of an antiferromagnet always acquire permanent magnetic moments of the order μ A N p , where μ A is the magnetic moment of a sublattice atom, N the total number of such atoms in the grain, and 1 3 < p < 2 3 . By now this fact is well established for ferric hydroxides like ferrihydrite and mineral cores of ferritin [28]. It seems that in nickel hydroxide the electron magnetic resonance could be an appropriate tool to provide the evidence for the existence of such states and to monitor their behavior as a function of temperature, as it is for ferrihydrite [29].

Conclusions

Individual magnetic and orientational behavior of nickel hydroxide nanoplatelets is studied by means of stabilizing particles in a liquid matrix. For characterization of the particles three independent techniques are used: SAXS, SQUID and colloid magneto-optics. Under an applied field the suspended platelets orient their faces normal to the field direction; the order parameter attained under H = 640 kA/m at room temperature is about 0.1. Analysis of this effect is performed by combining the results of all the aforementioned measurements. It reveals that nickel hydroxide in a platelet phase is paramagnetic with a pronounced anisotropy of the intrinsic susceptibility. The field-induced orientation in a dilute dispersion of Ni(OH) 2 has the same symmetry as the liquid-crystalline ordering arising spontaneously in dense assemblies of same particles. This infers that a magnetic field is an important controlling factor which can either reduce the critical concentration (at given temperature) or increase the reference temperature (at given concentration) of the liquidcrystalline transition in Ni(OH) 2 platelet systems. Lowtemperature magnetic measurements indicate that below 100 K the particles develop metamagnetic and antiferromagnetic properties. Further studies of these phenomena seem to be very promising.

Fig. 1 .

 1 Fig. 1. High-resolution transmission electronic microscopy picture of a side view of a platelet; the measured interplane distance of the crystal lattice is 4.7 Å; the white bar corresponds to 10 nm.

Fig. 2 .

 2 Fig. 2. Transmission electronic microscopy: diffraction (a, b) and pictures (c, d); in (a, c) the beam goes through an ensemble of platelets, in (b, d) the beam is focused on one platelet.

Fig. 3 .

 3 Fig. 3. SAXS scattering for a Ni(OH)2 platelet suspension at Φ = 2.9 • 10 -2 and [cit] free = 10 -2 mol/L in a field H = 640 kA/m;(a) the scattering pattern elongated along the field direction (to better reveal the anisotropy, a greyscale with only 8 levels is used); α is the azimuth angle defined in the text; (b) the angular dependence of the scattered intensity I at a given scattering vector q = 0.29 nm -1 fitted with the aid of formula (1). Its experimental value at H = 0 is 680 a.u.

Fig. 4 .

 4 Fig. 4. Schematic view of orientational distributions of the platelets in a liquid matrix corresponding to the limiting values of the orientational order parameter (left to right): easy-axis, easy-plane, disordered.

Fig. 5 .

 5 Fig. 5. Field dependences: (a) of the orientation order parameter S2.For each value of the field, S2 is averaged over 0.15 nm -1 < q < 0.3 nm -1 and 0.4 • 10 -2 < Φ < 2.9 • 10 -2 ; measurements are done at [cit] free = 10 -2 mol/L. The full line corresponds to S 2 = 2.9 • 10 -7 H 2 , where H is taken in kA/m; (b) of the reduced magnetic birefringence Δn/Φ for three different samples at Φ = 0.5% ± 0.1% with 5 • 10 -3 mol/L < [cit] free < 10 -2 mol/L.

Fig. 6 .

 6 Fig. 6. Magnetization of a dispersion of platelets at Φ = 1.4% and [cit] free = 2.6 • 10 -2 mol/L at room temperature. Inset: temperature dependence of the initial susceptibility of the same sample measured at H = 4.1 kA/m under field-cooled (circles) and zero-field-cooled (squares) regimes.

  Partial support of the work by ECO-NET Program No. 16274QL is gratefully acknowledged. Y.R. acknowledges also support on the part of UrB/SB RAS Collaboration Program and RFBR project 07-02-96026. P.L. acknowledges support from CAPES 1973/01-0, all the authors thank C. Bourgaux for providing help on D24 beam line in LURE and N. Malikova for careful re-reading of the manuscript.