Robert Charles Pond
email: r.c.pond@exeter.ac.uk

John Price Hirth

Elastic and Plastic Aspects of Martensitic Transformations

Keywords: crystalline interface, deformation, dislocation theory, elasticity, phase transformations, plasticity user supplied): crystalline interface, deformation, dislocation theory

  L'archive ouverte pluridisciplinaire

The objective of this paper is to discuss elastic and plastic aspects of martensitic transformations where the interface between parent and product crystals is modelled in terms of dislocations. Bullough and Bilby [1] presented the first such treatment in 1956. Their "simple glissile interface" comprises a "surface dislocation" corresponding to a single set of dislocations which move conservatively and thereby cause the necessary lattice-invariant deformation (LID), represented in matrix notation by 2 P (for convenience, we use the later nomenclature introduced by Wayman [START_REF] Wayman | Introduction to the Crystallography of Martensitic Transformations[END_REF]). Motion of the interface also produces the lattice deformation, B , but this is not modelled explicitly in terms of defects. These two homogeneous deformations, along with a rigid body rotation, R , combine into the shape deformation, 1 P . Under these assumptions, the interface is an undistorted plane, and the shape deformation is an invariant-plane strain (IPS) given by, 2 1

RBP P =

.

(

) 1
They showed furthermore that the total Burgers vector content for such an interface is zero, consistent with the absence of stress in either crystal. These conclusions are in complete accord with the phenomenological theory of martensite crystallography (PTMC), developed contemporaneously by Wechsler, Lieberman and Read [START_REF] Wechsler | [END_REF] and Bowles and MacKenzie [4], which also considers two crystals to be homogeneous and separated by an invariant-plane (IP). In this framework, the total shape change on transformation can be written ( )

' 1 dp I P + = , ( 2
)
where I is the unit matrix, d is the displacement of all points at unit distance from the interface, and ' p is the transpose of the interface normal [START_REF] Wayman | Introduction to the Crystallography of Martensitic Transformations[END_REF][START_REF] Wechsler | [END_REF]4]. Note that, in this view, the interface plane, ' p , is determined geometrically by the affine rotation and deformations R , B , and 2 P . The former two also determine the orientation relationship (OR) between the adjacent crystals.

Since the pioneering works [1,[START_REF] Wechsler | [END_REF]4] a large number of microscopical investigations of interface structure have been undertaken (for a review see [5]).

Several studies using transmission electron microscopy showed that martensitic interfaces actually comprise two (or more) sets of discrete defects. One set of dislocations (slip or twinning) corresponds to the LID as envisaged above. The second set comprises transformation dislocations [START_REF] Christian | The Theory of Transformations in Metals and Alloys[END_REF], or disconnections [START_REF] Hirth | [END_REF], line defects with both dislocation and step character. A schematic view of such an interface is depicted in figure 1, as has been discussed qualitatively by various authors [8][9][START_REF] Olson | Dislocations in Solids[END_REF][START_REF] Dahmen | [END_REF][12]. Recently, Pond, Hirth and co-workers [13] have presented a quantitative treatment of such transformation interfaces, known as the topological model (TM). In this model, the interface is partially-coherent, the dislocation and disconnection arrays delineating coherent terraces, as in figure 1. The mechanism of transformation is lateral motion of the disconnection array along the interface, and in general must be conservative.

Figure1 here.

In this paper we demonstrate that plastic transformation strains arise in the TM through motion of the defect arrays, producing a shape deformation in agreement with the PTMC. However, the interface plane and OR are determined by superposed plastic and elastic fields of the defects that comprise the stationary interface, and we

show that small differences in the habit planes predicted by the two approaches arise as a consequence of the elastic fields.

According to the topological model, the elastic distortion field of a static interface can be written as

( )           - - - +           = Γ 0 0 0 ' , ' , '
φ φ φ φ φ φ ε ε ε ε ε ε ε ε ε . ( 3
)
The strain and rotational distortion components of ( )

' , ' , ' z y x s Γ
vary from point to point, describing both the short-and long-range fields produced by superposition of the coherency strain and the elastic field of the defect networks. (The ' , ' , ' z y x coordinate frame defines the habit plane, as opposed to the terrace plane frame z y x , , as illustrated in fig. 1). At short-range all strains are partitioned between the phases in a manner depending on their relative elastic compliances and the boundary conditions. In the case where the coherency strain is completely accommodated by the network, all component strains except ' zz ε vanish at long-range. At short-range, these

strains alternate in sign along ' x and ' y respectively, in anti-phase across the habit plane, and diminish to zero when Elastic aspects are discussed in section 2. The method of linear superposition of elastic fields is used in section 2.1, and a non-linear solution, more appropriate for larger strains, is described in section 2.2. Plastic transformation strains are discussed in section 3. A comparison of TM and PTMC predictions is presented and discussed in section 4. The topological character (Burgers vector, b, and step height, h) of disconnections is summarised in Appendix A, along with the reference states employed.

Elastic Aspects

For clarity, we consider a transformation where no LID needs to be invoked because the two crystals are naturally coherent parallel to ' , x x . The TM is not restricted to such cases; Ma and Pond [START_REF] Ma | [END_REF] for example have employed the TM to analyze habit plane variants in Fe alloys that include twist rotation in the habit plane.

The TM and PTMC apply to a planar interface between two semi-infinite crystals, a parent β phase and a product α phase. The TM method involves selecting a terrace plane of good match. Then candidate disconnections, b and h, are selected from the coherent dichromatic pattern (see appendix), where the β andα crystals are strained to give perfect matching on the terrace plane, with a coherency strain a a / 0 ∆ = ε in the y direction. This strain is partitioned to strains 2 / 0 ε ± in the two phases. Here a is the atom spacing in the y-direction in the coherent dichromatic complex and a ∆ is the difference between atom spacings in the natural dichromatic complex. The transformation proceeds by the motion of a sequence of such defects, figure 2. Once the transformation stops, the equilibrium array of disconnections is determined by adjusting their spacing so that there are no long-range elastic stresses or strains for the array. For this aspect, coordinates ( ) are coherency stresses associated with the terrace segments that can be described as arising from continuous infinitesimal coherency dislocations [16]. No loss in rigor is entailed if one represents the integral of this distribution over D d as a discrete coherency dislocation with Burgers vector cy b defined using the natural dichromatic pattern (see Appendix A) . The "overlap" step h (see appendix) is assumed to be normal to the terrace plane: modifications when this is not the case are discussed in [17].

Figure 2 here.

Linear Superposition Solution

The b z component of the disconnection has magnitude ( )

α β h h -
, so, when this is zero, the Burgers vector of the disconnections in figure 2

is ] 0 , , 0 [ y b = b
, and that for 2 is the equilibrium value that removes misfit. The net strain components are all zero, either in the (x, y, z) coordinates or when the strain matrix is transformed to any other set of coordinates, e.g. (x' ,y' ,z'). The equilibrium condition is that the strains associated with the disconnections,

the coherency dislocations is ] 0 , , 0 [ cy c b = b . Hence, cy y b b - = when D d in figure
L b y / ' ' 0 - = ε (really L b y 2 / 2 / ' ' 0 - = ε
, but the factors of 2 cancel), then leads to the expression for the equilibrium habit plane angle θ for this case, given by [13] (

)

1 2 0 tan tan - + = - h b b z y θ θ ε . ( 4
)
Figures 3 and4 here. . The tilt wall gives a rotation, corresponding to the OR,

( )       - - = - 2 cos sin 2 sin sin cos sin 2 0 1 θ θ ε θ θ θ φ h b b y z . (5)
In the small φ limit, this becomes

( ) θ φ tan / h b z = . (6)
In the rotated coordinates, the coherency strain is given by ( )

θ ε θ ν θ ε ε 2 0 2 2 0 ' 0 cos sin cos ≈ - =
, where ν is Poisson's ratio. The dislocation array has a Poisson term that nearly cancels the coherency term and the ratio

θ ν θ 2 2 sin / cos
is typically of the order 6 10 so the approximation is an excellent one.

Partitioning of elastic distortions is inherent in the above description in several ways. The coherency strain is partitioned to a value 2 / 0 ε ± in each crystal. The rotation φ was assigned to the product phase [13]. However, this was amended in [START_REF] Pond | Dislocations in Solids[END_REF] to include the partitioning of the rotation which is necessary as shown in figure 5.

Figure 5 here.

Physically, the rotation φ arises from the summed distortions (∂u i /∂x j ) of the defects. At long range the strain parts of the distortions vanish and only pure rotations are left [17]: these partition equally in the isotropic homogeneous elastic limit. The , with φ positive. The sign of φ would be negative when β h < α h .

Figure 6 here.

The results are summarized in figure 7 where habit plane inclinations are plotted as a function of φ is only significant for cases of largeφ . However, it is significant in the large φ case as a first-order correction, unlike the second-order effects discussed next. As a final illustration, the steps in figure 6a and b can be offset by half the defect spacing without changing associated long-range elastic fields (plastic shear of an equilibrium interface leaves the same elastic equilibrium interface). This is a hypothetical case, which would have higher energy than that of Figure 5, but it again demonstrates the necessity of partitioning as in the TM. The fields of part of the arrays cancel and the excess array has a net Burgers vector b z ', which partitions to b z '/2 as in the TM as described previously. We present these simple physical descriptions to emphasize that the overlap step height is the natural reference state for the superposition problem.

Figure 7 here.

As a specific example, we consider the transformation in "Ti" where the 0.86% mismatch in the x direction is neglected. For this example h α = 0.2556 nm and h β = 0.2678 nm. The values of

β β α α β η ω ω φ θ , ,

The Nonlinear Solution

The superposition solution is accurate to 3 or 4 places in predictions ofω provided that φ  < ≈ 0.8°. This suffices since it is within the experimental scatter in ω measurements, and of the order of the uncertainties in ω values in either the TM or the PTMC because of scatter in lattice parameter measurements. The superposition solution has the advantage of being described by the simple analytic equations ( 4) and (5). In this region the habit planes are those shown between the vertical bars in figure 7, and vary in a nominally linear manner with respect to the variable α β h h / . For a greater value ofφ , nonlinear effects become important and there is curvature in the lines. For the "Ti" example, the interface is described as α ω and β ω as shown in

θ φ θ ε tan )] 2 / tan( [ 0 - + = - . (7)
Equation ( 2) is similarly replaced by . 2

) 2 / cos( ) 2 / sin( 2 tan ) 2 / cos( )] 2 / sin( ) 2 / cos( [( sin 2 0 1       - - - - - - -φ θ + .
Evidently, equations [START_REF] Hirth | [END_REF] and (8) reduce to (4) and (5) in the small φ limit. Unlike the simple superposition case, equations ( 7) and ( 8) must be solved simultaneously by numerical or iterative methods. Explicitly, for 1 + i θ and i φ , one can start with the linear result 1 φ φ = i and iterate; a possible procedure is set out in Appendix B. The results for the "Ti" case are listed in Table 2, showing that the simple superposition results are satisfactory for the small φ case, i.e. for most cases analyzed or measured. Only in exceptional cases such as the α β h h / = 1.15 case, would the nonlinear corrections be needed. An example would be the martensite transformation in Pu-Ga alloys, where °≈ 6 φ and α β h h / = 1.2 [18,19]. In such a case, however, nonlinear elastic corrections would have to be considered as well as the nonlinear geometric corrections. The most obvious effect is that the coherency strains are no longer described by linear elasticity. The strains arise for stresses applied to the free surfaces normal to y to satisfy the free surface boundary condition. The strains are related to these stresses by third and higher order nonlinear elastic constants, an effect that becomes important for strains greater than about 0.5%. Anisotropic elastic effects lead to stress tensors that include more components, perhaps all six, for edges and screws, complicating defect interactions, modifying coherency relaxation and producing results analogous to those for isotropic elastic constants at an interface between elastically inhomogeneous crystals. The most important effect in the latter case is to modify partitioning. For example, in the limit that the matrix was infinitely stiff, all among the coherency stresses, the defect stresses and the tilt rotation, if present, would partition to the product phase [20]. Even in the inhomogeneous isotropic elastic case, similar changes in partitioning occur [START_REF] Dundurs | Mathematical Theory of Dislocations[END_REF]. In both cases [20,[START_REF] Dundurs | Mathematical Theory of Dislocations[END_REF], distributions of interface line forces may also be needed for the equilibrium description.

Table 2 here.

All of the effects discussed to this point apply to the analysis of the transformation defects only. The same considerations of partitioning of tilt rotations apply in the TM analysis of LID. As mentioned previously, twist rotations and their partitioning are needed when the disconnections or LID include screw components of dislocations [START_REF] Ma | [END_REF]. Interactions between disconnections and LID defects, usually neglected, can also lead to complications [START_REF] Hirth | [END_REF]: line force distributions can again arise in the interface.

We mention these possibilities to reinforce the earlier conclusion that the secondorder differences between the PTMC and TM models may not be significant compared to uncertainties arising from approximations in the models.

Plasticity Effects

In this section we show that the plastic portion of the TM description of the transformation strain agrees with the shape transformation obtained from the PTMC.

For the static equilibrium of Section 2, the elastic strain fields of the disconnections were important in relieving coherency strains. Moving disconnections also serve to produce the plastic transformation strain, converting the parent β to a In the present case, where no LID is invoked, the PTMC shape deformation, equation ( 1), is simply given by RB P = 1 . In other words, the shape deformation expresses the combined effect of the lattice deformation and the OR. The PTMC result is therefore equivalent to the plastic TM transformation deformation combined with rotation of the product byφ . As explained above, in the TM the rigid-body rotation φ is associated with the elastic distortion fields of the defects at the interface, which is unchanged on average as the train of defects moves. Thus both the rotations connected to the elastic distortions and the associated elastic distortion partitioning relative to the current interface are unchanged by the motion of the habit plane. The PTMC and the TM agree except for the elastic partitioning. That is, the two models agree on the magnitude of the components of the shape change. However the PTMC applies them relative to the interface cut on β η while the TM applies them relative toθ .

          =           = Γ δ ε γ , ( 9
Using expressions (2.7) and (2.8) of reference [START_REF] Pond | Dislocations in Solids[END_REF], one can readily confirm that both the "glide" of disconnections across the terraces and the "climb" of the coherency dislocations perpendicular to this are conservative.

Discussion

The pioneering works of Wechsler, Lieberman and Read [START_REF] Wechsler | [END_REF] and Bowles and MacKenzie [4] are phenomenological and pertain to crystals where elastic displacements are suppressed. Similarly, Bullough and Bilby's surface dislocation model [1] suppresses elastic relaxations, and is hence congruent with the PTMC. By contrast, the TM [13,[START_REF] Pond | Dislocations in Solids[END_REF] is aetiological, envisioning the interface structure as being partially coherent and transformation to occur by the motion of defects. In this approach plastic and elastic effects are distinguished. Broadly speaking, the plastic transformation strain and OR predicted by the TM for the case considered here is the same as that predicted by the PTMC. Thus, the shape transformation of the PTMC, equation (1), is equivalent to the plastic transformation strain, expression [START_REF] Dahmen | [END_REF], in combination with the OR, expression (4), of the TM. However, consideration of elastic relaxations in the TM describes the partitioned interfacial displacement field and shows that the habit plane consequently differs slightly from the PTMC prediction. One consequence is that the shape transformation according to the TM can no longer be expressed as a simple IPS like expression [START_REF] Wayman | Introduction to the Crystallography of Martensitic Transformations[END_REF]. Although the treatment presented here is for a case where no LID is invoked, it can be readily extended to the more general case. The plastic transformation strain produced when slip or twinning dislocations follow the interface can be expressed as an IPS analogous to equation [START_REF] Dahmen | [END_REF]. Moreover, their presence affects the OR and habit plane. The combination of IPSs due to disconnection and dislocation motion, along with the overall OR, can be expressed as an ILS, RB. An example of the application of the TM to a transformation involving LID is the case of ( ) 575 lath martensite in ferrous alloys [START_REF] Ma | [END_REF]. Here, the disconnections and slip dislocations in the interface have predominantly screw character, producing a partitioned twist displacement field near the interface.

The phenomenological theories cannot explain satisfactorily why habit planes sometimes vary with small composition changes of the phases (for example in Fe-C alloys) or with the mode of formation (by cold-work or cooling for example) [START_REF] Wayman | Introduction to the Crystallography of Martensitic Transformations[END_REF][START_REF] Christian | The Theory of Transformations in Metals and Alloys[END_REF]. In these approaches the lattice deformation normally chosen, B, is that corresponding to minimum strains. By contrast, in the TM this is determined primarily by the Burgers vector of the operative disconnections, as encapsulated in expression [START_REF] Olson | Dislocations in Solids[END_REF]. For a given transformation a multiplicity of disconnections may arise on a viable terrace plane, and hence distinct disconnection/dislocation arrays, leading to distinct habit planes, might arise in differing circumstances. The multiplicity of potential disconnection/dislocation combinations is particularly marked in the case of ferrous alloys [START_REF] Ma | [END_REF]. Partitioning of coherency strains has been recognised by many previous workers, e.g. [16]. However, the partitioning of tilt in phase transformations has only been recognised recently [17]. This is an irreversible thermodynamic relaxation in the strain energy of the system. Its incorporation into the TM augments the PTMC, producing an added rotation of the α β / habit planes, by 2 / φ ± when φ is small. An unrelaxed PTMC interface is depicted in figure 9, but there is no dislocation arrangement that can be superimposed on the TM structure to reproduce it. Extrinsic dislocations could be added to make, say, θ η α = when In many cases the rotation of 2 / φ ± is a small emendation but in some it is not.

We emphasize that these calculated differences apply to the idealized case of a large bicrystal with a planar interface. The partitioning arises from the superposition of surface tractions on the free surfaces normal to y in order to satisfy the free surface boundary condition. This is a good approximation for the case of lath martensite with a large aspect ratio. For lenticular plates, strain partition is resisted by compatibility forces at the plate tips. There the rotations and strains in the parent tend to be suppressed. This would be favored particularly for plates with small aspect ratios.

Compatibility strains '

xx ε would modify results for both the TM and the PTMC.

Rotations would tend to be suppressed in the TM changing the partitioning in the direction to more closely agree with the PTMC. Of course, many other factors enter when plate curvature and compatibility effects become important. We are pleased to acknowledge that the nonlinear treatment presented here was stimulated by a long interchange with P.M. Kelly.

Conclusions

tan ] ) 2 / tan( ) [( tan ] ) 2 / tan( [ 0 1 0 1 2 = + + - + + + + ε θ φ ε θ φ h b b h b b i y i z i z i y . (B1)
Inserting i φ one can solve this expression for the up-dated 1 + i θ . Then the updated 1 + i φ can be found using eq. ( 8) in the form ( ) Step height ratios and habit plane parameters for "Ti" obtained by the nonlinear method. All angles are expressed in degrees.

i i i f φ θ φ , 1 1 + + = .

  -range strains are all zero when this is satisfied. The latter result is the reason why the TM and PTMC agree exactly when b z = 0. Of course, there would be local dislocation strains within a distance shown in figure2, creates a tilt wall with angleφ . If the disconnection array is considered as the superposition of a set of pure steps and a tilt wall, the simple analytical equations of[13] apply. To appreciate this point, we can order, this is equivalent to cancelling z b , leaving just the steps h, figure4. The basic relation describing the removal of misfit,

  The extrinsic dislocations ' ez b are removed from figure3, restoring the tilt wall with the attendant partitioned distortions, which, at long-range, reduce to the partitioned rotations 2 / φ ±

φ

  in each crystal. In the original treatment, the total

  local strains are present whether the tilt components are ' z b or z b , figure 5. The differences in magnitude of strain, distortion or rotation correspond to a factor θ 2 cosand vanish to linear order. Partitioning is also nearly the same. Consider the specific case α h < β h . One can think of the interface separated into two crystals with the free surface step heights α h and β h , figure6a. These can also be represented as the superposition of pure steps of heightα h h = , figure6b, and a set of little virtual steps associated with the partitioned dislocations, figure 6c. The superposition of the two arrays shows that, relative to figure6b, the superposed tilt causes the habit plane in β to increase in this example, while that in α will decrease. The final habit planes, ω , are then given by the superposed results as 2

Figures 8 .

 8 Figures 8. A key point is that L in figure 2 becomes ) 2 / cos( / φ L L n = . Also the

δ

  product α. The dislocation component of a moving disconnection sweeping an area A produces an engineering strain per unit volume A b D = Γ . Here, n A A = with n a unit normal to A. Hence a straight disconnection of length l moving by y δ in the y for a defect of unit length. Motion of a train of disconnections as in figure 2, with D d / 1defects per unit length in the y direction, produces an engineering strain,

) with n z = 1 .

 1 figure A1(a).

  equal θ and not equal the value η β = θ + ϕ of figure7. Similarly in the other cases, one or the other habit plane could be made to agree with figure9, but both could not together so agree.

Figure 9

 9 Figure 9 here.

  Elastic and plastic strains are distinguished in the TM treatment of martensitic transformations. Thus the TM extends previous methods (surface dislocations and phenomenological theories) by including elastic relaxation by partition of the tilt wall rotations,φ , when such walls are part of the interface structure. Such relaxations lead to rotations of the habit plane by 2 / φ away from classical predictions. A nonlinear elastic/plastic analysis to determine these rotations is presented, enabling application of the TM to larger values ofφ . In the present work a simplified case is considered where no LID is invoked. When LID is involved, partitioning concepts would also apply. Plastic deformation arises in the TM through motion of defects at the interface. Transformations exclusive of LID effects are considered in detail in this paper.The plastic transformation strain, in combination with the orientation relationshipφ , predicted by the TM agrees with the previous methods.

  University of Liverpool. During that period David was not only a foremost academic but also a supportive colleague and excellent administrator.

Figure 1 .

 1 Figure 1. Schematic illustration of a parent( β )-martensite(α ) interface showing the terrace segments and defect arrays [13, 14 ]. Coherently strained terraces are reticulated by arrays of disconnections (b, h) with spacing d D and crystal slip dislocations (b, 0). The terrace and habit (primed) coordinate frames are shown and the line directions of the disconnections, ξ ξ ξ ξ D , and dislocations, ξ ξ ξ ξ L , are parallel to x and close to y' respectively.

Figure 2 .

 2 Figure 2. Schematic illustration showing the disconnection content of an interface,with Burgers vector components resolved in the terrace (upper) and habit plane (lower) frames[13,[START_REF] Pond | Dislocations in Solids[END_REF]. The terrace plane is inclined at an angle θ to the horizontal habit plane. Coherency strain is represented by the equivalent "coherency" defect content, cy b . The x-axis points out of the page.

Figure 3 .

 3 Figure 3. Schematic illustration of the tilt boundary along the equilibrium interface: the Burgers vectors of the component dislocations are ' z b . A hypothetical

Figure 4 .

 4 Figure 4. Schematic illustration of an array of "overlap" steps defining the habit plane inclination θ .

Figure 5 .

 5 Figure 5. Schematic illustration of the partitioning of distortions in both adjacent crystals arising from the dislocation components (a) z b , and (b) ' z b . The dashed lines represent bending of planes (elastic distortions) in the vicinity of the dislocations.

Figure 6 .

 6 Figure 6. Schematic illustration of free-surface step heights, (a) for the separated β and α crystals, (b) representation of the common "overlap" step height, and (c) array of virtual steps corresponding to the partitioned tilt wall, φ .

Figure 7 .

 7 Figure 7. Variation of the TM habit plane angles 2 / φ θ ω α -= (dashed line)

Figure 8 .

 8 Figure 8. (a) Schematic illustration of the habit plane angle 2 / φ θ ω α -= and (b)

Figure 9 .

 9 Figure 9. Schematic illustration of a habit plane between unrelaxed crystals according to the PTMC.

Figure

  Figure A1 (a) Illustration of natural bcc and hcp "Ti" lattices exhibiting the Burgers orientation relationship. (b) Schematic illustration of the formation of a disconnection; the b of the defect with respect to the coherent reference state is the Volterra operation required to create equivalent interfaces on either side of the crystals' surface steps. The 'overlap' step height of the defect is the smaller of the two surface steps, h(α) in this case. The sense vector ξ ξ ξ ξ points out of the page.

  Figures

Figure 7 Page

 7 Figure 7

Table 1 .

 1 Convergence is rapid since1 Step height ratios and habit plane parameters for "Ti" obtained by the linear superposition method. All angles are expressed in degrees.

	i φ	+

Table 2 .

 2

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

Acknowledgment

It was a privilege to contribute to the symposium celebrating David Bacon's outstanding career. One of us (RCP) worked with David for over thirty years at the

Appendix A

The unit cells of bcc and hcp Ti viewed along , eq. ( 7) can be written


Volker Mohles
email: mohles@imm.rwth-aachen.de

Barend Thijsse
email: b.j.thijsse@tudelft.nl

Emmanuel Jannot

Mohles I Volker

Calculation of the elastic energy of a solid solution

Keywords: computer simulation, precipitation user supplied): micro elasticity theory microscopic elasticity theory; precipitation mc, manuscriptcentral, com/pm-pml

published or not. The documents may come

Introduction

To perform a Monte Carlo (MC) simulation, a description of the energy of an atomic configuration is required. A very popular method is found in the use of a cluster expansion. An advantage of this approach remains its low computational cost for a simulation at the atomic scale. In a cluster expansion, the lattice positions are usually fixed such that they represent a perfect lattice. This rigid lattice is not adapted for the modelling of a lattice containing defects. But when modelling precipitation near a dislocation, for instance, a different approach should be considered. In a continuum medium, Hu and Chen [2] used the Macroscopic Elasticity theory of Khachaturyan to consider the dislocations' elastic stresses. At the atomic level, the Microscopic Elasticity (ME) theory of the same author [START_REF] Khachaturyan | Theory of Structural Transformations in Solids[END_REF] can be employed, according to current understanding. But in MC simulations by the present authors using the ME theory as defined by Khachaturyan, severe problems occurred.

In the quoted simulations, the precipitation of Guinier-Preston zones in AlCu alloys was to be simulated. As an example, Guinier-Preston zones of type I are single layers of Cu atoms in {100} planes. Their occurrence is known to be based on the large size mismatch between Al and Cu atoms: this mismatch prevents the formation of coherent voluminous particles because the corresponding elastic energy scales with its volume, and is too high. In contrast, the formation of plate-like precipitates is not prevented by elastic strains because these are similar to those of dislocations and hence scale more favourably with length. But in the quoted MC simulations using Khachaturyan's ME theory, always spherical precipitates emerged, even when the size mismatch between Al and Cu atoms was greatly exaggerated. This is a strong contradiction to both theoretical expectation and to well-known experimental findings. The motivation for the present work was to find the cause for the discrepancy. It lead to the development of a corrected version of the ME theory.

To assess the quality of the ME theory (starting with the original), the elastic energies were also calculated in a more straightforward way for comparison. In this approach, the energy changes in every time step of an MC simulation were derived from a relaxation simulation of "atoms" interacting by harmonic potentials. The relaxation was done using Molecular Statics (MS), i.e. a Molecular Dynamics [START_REF] Allen | Computer Simulation of Liquids[END_REF] program combined with a conjugate gradient algorithm [START_REF] Atkinson | An Introduction to Numerical Analysis[END_REF] for the atomic relaxation.

When harmonic potentials are employed for the description of the atomic interactions, the energy of a bond between two atoms of nature X and Y is given by: where E X-Y represents the ground energy, κ X-Y the bond stiffness, and d X-Y the equilibrium length of the bond between atoms of type X any Y. The total energy of a configuration is then given by:

V X-Y (∆x) = E X-Y + κ X-Y 2 ∆x -d X-Y ( ) 2 (1)
E (MS) = 1 2 r ∑ r′≠ r ∑ E n(r ),n(r′) + κ n(r ),n(r′) 2 r′ -r -d n(r ),n(r′ ) ( ) 2       (2)
where the field n(r) indicates the nature of the atom sitting at position r. Compared to the application of the ME theory, the relaxation simulation and evaluation of equation ( 2) takes much more time to compute. But the MS method and its implementation are very simple and therefore virtually guaranteed to be free of programming or interpretation errors. Therefore this method is useful to assess ME theories (and implementations thereof).

The coupled MC-MS simulations of precipitation in AlCu alloys in fact yielded the formation of Guinier-Preston-zones of type I, in full agreement with both theoretical expectation and experimental findings. Therefore the ME theory has been revisited.

Derivation of the corrected ME theory

The binary solution case

The aim of the ME theory is to provide an evaluation of the change of stored elastic energy by a configuration change, like in Fig. 1 where the solute atom 13 exchanges its position with the matrix atom 9. This atomic move modifies the local elastic strains and consequently changes the amount of elastic energy stored in the lattice.

-insert fig. 1 about here -

To calculate the change of stored elastic energy, the transformation from I to II is decomposed in three steps (Fig. 2). In a first step, the atoms are repositioned from their equilibrium positions (no stresses) to their positions in the reference lattice (symbolized by dashed lines in Fig. 2). The reference lattice corresponds to the undistorted, periodic lattice, in which all atoms experience forces from their neighbours. This step requires adding the energy ∆E I→I' to the system. Then, the solute atom 9 and the matrix atom 13 are exchanged in the reference lattice. This step also modifies the energy of the system by a value ∆E I'→II' . Finally, configuration II' is relaxed to II, where the atoms reach their new equilibrium positions. This last step reduces the energy of the system by ∆E II'→II .

-insert fig. 2 about here -

The first and the third steps are similar in nature: II'→II corresponds to a relaxation while I'→I corresponds to the inverse of a relaxation. Hence, the energy changes by these transformations are calculated the same way. When r I´ represents the atomic positions in configuration I', the positions in I can be described as r I´+ u, where u are the relaxation displacement vectors. Considering the harmonic approximation, the energy of configuration I is evaluated by a first order Taylor expansion:

E I ( ) = E I′ ( ) -F I′ ( ) ⋅ u + 1 2 u T ⋅ D ⋅ u ( ) (3)
where F(I´) represents the net forces acting on the atoms in the reference lattice positions and D is the force constants matrix [START_REF] Kaxiras | Atomic and Electronic Structure of Solids[END_REF]. Since configuration I corresponds to an equilibrium state, the displacements u minimize equation (3), leading to F(I´)=D⋅u.

By introducing the Green tensor G defined as G=D -1 , this can be transformed to u=G⋅F(I´), which in turn can be inserted for in equation [START_REF] Allen | Computer Simulation of Liquids[END_REF]

I′ →I = E I ( ) -E I′ ( ) = -1 2 F I′ ( ) T ⋅ G ⋅ F I′ ( ) ( ) (4)
The determination of (∆E) I'→I is not straightforward. The dimension of the vector F(I´) is 3N where N is the number of atoms contained in the configuration, making the size of the force constants matrix D to be 3N×3N. This large matrix size prevents from working in real space. But the determination of (∆E) I'→I becomes efficient in reciprocal space. The calculation is performed in the following for an alloy described by harmonic potentials. Using this description of the atomic interactions,

a ref = 1-c sol ( ) a mat + c sol a sol ( 5
)
with a mat and a sol being the lattice parameters of the matrix and of a pure solute lattice, respectively. These lattice parameters are naturally linked to the equilibrium lengths of the pure element bond. In case of an fcc lattice, d=2 -1/2 a. By definition, the reference lattice is not in equilibrium: the force acting on some atoms is not zero. As

f mat -mat x ( ) = κ 2 -1/ 2 a ref -d mat -mat ( ) x 0 (6) f mat -sol x ( ) = κ 2 -1/ 2 a ref -d mat -sol ( ) x 0 (7) f sol-sol x ( ) = κ 2 -1/ 2 a ref -d sol-sol ( ) x 0 (8)
Here, x 0 denotes the unit vector of x. The net force F(r) acting on an atom at position r by all neighbours is dependent on the local configuration. When all atomic interactions are limited to the nearest neighbours like with harmonic potentials, four cases are possible:

(a) A matrix atom (No. 3 in Fig. 3-a) is only surrounded by matrix atoms. Each neighbouring matrix atom acts an interatomic force f mat-mat on atom 3 (repulsive if

d mat-mat >2 -1/2 a ref as in AlCu)
. But the sum of the interatomic forces, the net force F, vanishes.

(b) A solute atom (Fig. 3-b) is only surrounded by matrix atoms. Atom 3 senses attractive interatomic forces f mat-sol coming from the neighbouring matrix atoms. In such case, the net force felt by atom 3 also vanishes.

(c) A matrix atom (Fig. 3-c) has one solute atom as a neighbour (or several ones).

Here the net force does not vanish but is f net =f mat-sol -f mat-mat , because a matrix-matrix interaction is removed and a solute-matrix interaction is added compared to the case (a).

(d) A solute atom (Fig. 3-d) has one solute atom as a neighbour (or several ones). This time, atom 3 experiences a net force f net =f sol-sol -f sol-mat because a solute-matrix interaction is replaced by a solute-solute one compared to the case (b).

-insert fig.

F(r) = - f mat -sol net (r -r′) ⋅ c(r′) ⋅ 1-c(r) ( ) r′ ∑ - f sol-sol net (r -r′) ⋅ c(r′) ⋅ c(r) r′ ∑ (9)
where The concentration field c(r) in real space has an image c ^(k ) in reciprocal space for a wave vector k restricted to the first Brillouin zone [START_REF] Kittel | Introduction to Solid State Physics[END_REF]. The transformations from one space to the other are defined by

f mat -sol net = f sol-mat -f mat -
c(r) = 1 N c ^(k ) ⋅ e i k⋅ r k ∑ (10)
and

c ^(k ) = c(r) ⋅ e -i k⋅r r ∑ (11)
Using equation ( 9), the image of the net force in reciprocal space is

F ^(k ) = - f mat -sol net (r -r′) ⋅ c(r′) r′ ∑       ⋅ e -i k⋅ r r ∑ - f sol-sol net (r -r′) -f mat -sol net (r -r′) ( ) ⋅ c(r′) ⋅ c(r) r′ ∑       ⋅ e -i k⋅r r ∑ (12)
The first term in equation ( 12) is easily derived because it represents a convolution in real space. Its expression in reciprocal space reads f ^mat -sol net (k) ⋅ c ^(k ) [START_REF] Kreyszig | Advanced Engineering Mathematics[END_REF] where ( ) X′→X from a stressed state X´ in the reference lattice to the relaxed state X is given by ∆E ( )

f ^mat -sol net (k)
X′ →X = -1 2N F ^i (k) G ^ij (k) F ^j * (k)       k′ ∑ (13)
where i and j represent the axis coordinates and F ^ * (k ) is the conjugate of F ^(k ).

Expression (13) allows the determinations of ∆E I→I' and ∆E II'→II .

The remaining quantity to evaluate is ∆E I'→II' , the change of elastic energy when passing from configuration I' to II' in the reference lattice. This change is proportional to the changing number of solute-solute bonds. When the passage from I'

to II' leads to the creation (or the annihilation) of M solute-solute bonds, the change of elastic energy from I' to II' equals

∆E I′ →II′ = M ⋅ E mat-mat ref + E sol-sol ref -2 E mat -sol ref ( ) (14)
where the E ref quantities account for the elastic energy stored in the stretching of bonds to accommodate to the lattice parameter a ref in the reference lattice. Their values are given for harmonic potentials as

E mat-mat ref = κ 2 2 -1/ 2 a ref -d mat -mat ( ) 2 (15) E mat-sol ref = κ 2 2 -1/ 2 a ref -d mat -sol ( ) 2 (16) E sol-sol ref = κ 2 2 -1/ 2 a ref -d sol-sol ( ) 2 (17)
In another ME theory proposed by Cook and de Fontaine [START_REF] Cook | [END_REF], the elastic energy change ∆E I'→II' of equations ( 14) to (17) was accounted. This helped to obtain a better agreement with MS relaxations, as is shown later on, but otherwise this theory contained the same error as the more widely known one from Khachaturyan.

Extension to multicomponent alloys

Like the previous theory developed by Khachaturyan [1], the corrected ME theory can be extended to multicomponent alloys. The central point in this extension is the definition of an expression of the net force F(r) acting on position r of the reference lattice.

When the system is composed of solute elements of n different natures, n concentration fields are defined as c p (r)=1 if the atomic position at r is occupied by a solute atom of nature p (1≤p≤n), and c p (r)=0 otherwise. Using this convention, a matrix atom is located at r when c p (r)

p=1 n ∑ = 0 (18)
The net force applied on this matrix atom surrounded by solute atoms of different nature is then defined as:

F 0 (r) = - f p 0 (r -r′) -f 0 0 (r -r′) ( ) ⋅ c p (r′) r′ ∑ p=1 n ∑ ⋅ 1-c q (r) q=1 n ∑         (19)
where f p0 indicates the interatomic force between a solute atom of nature p and a matrix atom of nature 0, by convention. In case of a solute atom only surrounded by matrix atoms, the net force F(r) applied on it vanishes. When a solute of nature q is surrounded by other solute atoms, it experiences the net force

F q (r) = - f p q (r -r′) -f q 0 (r -r′) ( ) ⋅ c p (r′) r′ ∑ p=1 n ∑ ( 20
)
where f pq represents the interatomic force between a solute atom of nature p and another atom of nature q. All the interatomic forces are defined in an fcc lattice with lattice parameter a ref as: Uniting all the cases into a single expression for F(r), one obtains from equations ( 19) and (20):

f p q r ( ) = κ 2 -1/ 2 a ref -d p q ( ) r 0
F(r) = - f p 0 (r -r′) -f 0 0 (r -r′) ( ) ⋅ c p (r′) r′ ∑ p=1 n ∑ - f p q (r -r′) + f 0 0 (r -r′) -f q 0 (r -r′) -f p 0 (r -r′) ( ) ⋅ c p (r′) r′ ∑ c q (r) p=1 n ∑ q=1 n ∑ (22)
Once the expression F(r) is known, the same method as in the binary solution case can be applied to determine the relaxation energies ∆E I→I' and ∆E II'→II . The third energy term ∆E I'→II' corresponding to the change of stored elastic energy in the reference lattice is evaluated in a similar manner as in the binary case:

∆E I′ →II′ = M p q ⋅ E p q ref + E 0 0 ref -E p 0 ref -E q 0 ref ( ) q=1 n ∑ p=1 n ∑ (23)
where M pq represents the number of bond changes between atoms of natures p and q during the transformation from configuration I' to II'. The unit energies E p q ref are defined as:

E p q ref = κ 2 2 -1/ 2 a ref -d p q ( ) 2 (24)

Accuracy of the calculation

The assessment of the accuracy of the ME theory is performed by comparing the results achieved with the ME theory with those obtained using MS relaxations. When the same harmonic potentials are used as input, the agreement should be perfect. In the calculations below, harmonic potentials were employed with the parameters: d mat- mat =2.85Å, d sol-sol =d mat-sol =2.56Å, and κ mat-mat =κ mat-sol =κ sol-sol =1.4 eV/Å 2 .

A first consistency test was conducted by looking at the energy of a configuration containing only two solute atoms as a function of the distance between the two solute atoms expressed in the number of neighbour shells NN (Fig. 4). There is a perfect agreement between MS relaxations and all ME theories when NN>1

(errors below 1meV). The situation is however different for the first neighbour shell, where the ME theory of Khachaturyan provides a value 30meV lower than the configuration energy obtained with MS relaxations. For the ME theory of Cook and

de Fontaine [START_REF] Cook | [END_REF] as well as the present one, the agreement remains excellent even for NN=1.

A second test consists in adding one more solute atom. Starting from a configuration where two solute atoms in positions r 1 and r 2 are first nearest neighbours, the third solute atom is introduced at different locations r 3 . Depending on the position of the third atom, the configuration energy is dependent on two quantities: the distance in neighbour shells NN between r 3 and r 1 , and the angle (r 2 -r 1 ,

r 3 -

-insert fig. 4 and tab. 1 about here -

The most striking difference between the previous ME theories and the version presented here is the capacity of the latter to simulate the correct shape of GP zones in AlCu alloys. As mentioned in the introduction, only spherical precipitates were obtained by MC simulations using the previous theories (both of them). On the other hand, plates located in {100} lattice planes (Fig. 5) as observed experimentally It should be noted that for certain kinds of physical mechanisms, the distinction between the various theories can be unimportant. Simulations of solute atom diffusion may be an example for this, when the solutes get in direct contact only rarely anyway. But in the present case of precipitation, where solute atoms tend to group and form nuclei for particles, it is particularly important to calculate the configuration energies of the nuclei correctly.

-insert fig. 5 about here -

Towards modeling of a real material

An additional improvement provided by the corrected ME theory is its ability to model more cases than before. In previous theories, only one misfit information could be managed: the matrix-solute bond always had the same equilibrium length as the solute-solute bond. In contrast, the new version differentiates the equilibrium length of the matrix-solute interaction from the solute-solute one. This new degree of freedom helps driving important conclusions on the effects of elastic strains on precipitation, as is shown in [START_REF] Jannot | Study of the effect of dislocations on precipitation kinetics by means of atomic simulations[END_REF].

D X-Y loc = 1 2 κ X-Y κ X-Y 0 κ X-Y κ X-Y 0 0 0 0         (25)
where X-Y represents any combination mat-mat, mat-sol or sol-sol. Using these local force constants matrixes, the net force acting on an atom at r can be expressed by

F(r) = D n(r ),n(r′) loc (r -r′) ⋅ u(r′) r′ ∑ (26)
where n(r) indicates the type of the atom at r. Since a unique small local force constants matrix does not exist, the convolution over r´ can't be simplified in reciprocal space and, hence, there is no valid formulation like F ^(k ) = D ^(k ) ⋅ u ^(k ) with a small 3×3 matrix D ^(k ) which could be afterward easily inverted to obtain a simple 3×3 Green tensor in reciprocal space. However, a difference in bond stiffness can be roughly accommodated by the following consideration. When different stiffness values are considered, the interatomic force f mat-sol defined by equation ( 7) is multiplied by the factor (κ mat-sol /κ mat-mat ). But, at the same time, the local force constants matrix is also locally multiplied by the same factor. Hence a stronger net force exists in a stiffer matrix when κ mat-sol >κ mat-mat . In such a situation, one can assume that the lattice displacements during the relaxation will be of the same amplitude. Consequently, the relaxation energy determined by equation ( 13) is only multiplied by a factor (κ mat-sol /κ mat-mat ). A simple transformation which leads to this result and remains inexpensive to implement in the relaxation energy calculation consists in the multiplication of the interatomic forces by a suitable square root, i.e. equation ( 27) may be used instead of (21): In order to model a material with a different anisotropy ratio, the description of the force constants matrix should be modified. This can be also done, as shown in [START_REF] Jannot | Study of the effect of dislocations on precipitation kinetics by means of atomic simulations[END_REF].

f p q (r) = κ p q κ 0 0 2 -1/ 2 a ref -d p q ( ) r

Conclusion

  only two parameters per bond type are necessary to evaluate the stored elastic energy of the configuration: the bond stiffness κ and the equilibrium length d of the bond. As the derivation of (4) in reciprocal space prevents from specifying different stiffness values (this aspect is discussed later on), all bonds are considered having the same stiffness κ. The first step in the evaluation of (∆E) I'→I consists in determining the net force F(r) felt by an atom located on the reference lattice I' at position r. As already indicated, the atomic positions in this lattice correspond to the ones of the undistorted lattice without macroscopic stresses. When a solid solution is considered, this reference lattice has a lattice parameter a ref which is a function of the solute concentration c sol . A simple expression for a ref valid for low solute concentrations reads

. 3 ,

 3 three different values of interatomic forces f may act on an atom in the reference lattice:

  mat and f sol-sol net = f sol-solf sol-mat according to the cases (c) and (d), respectively, and c(r) represents the solute "concentration": c(r)=1 if a solute atom is located at r, and a matrix one if c(r)=0. In Khachaturyan's version of the ME theory, the net forces were wrongly accounted. In the expression of the relaxation energy ∆F 6 relax in [1] (p. 451), they were expressed for the case of a binary phase as a sum over r´ of the term f Kanz. (r-r´)⋅c(r´), where f Kanz. are the Kanzaki forces. This definition of the Kanzaki forces is incorrect because, as outlined above, they are not a function of (r-r´) alone only but also of position r, i.e., the nature of the target atom itself (atom 3 in Fig. 3). The Kanzaki forces accept two values in case of a binary solution, depending on the nature of the atom at r: f Kanz. = f mat -sol net when the atom at r is a matrix atom, or (1-c(r))=1, and f Kanz. = f sol-sol net otherwise. The missing distinction of these two cases by c(r) in the definition of the Kanzaki forces explains the failure of the previous ME theory mentioned in the introduction. At a theoretical level, one may understand the origin of the error as a wrong definition of the stimuli in the Green's function analysis. When atoms showing a size misfit are considered, a stimulus actually corresponds to a small expansion or shrinking of an interatomic bond. Bus since a bond relates to two atoms, one has to count with care the stimuli. A simple spatial convolution involving the concentration field c(r) as in Khachaturyan's ME theory cannot achieve this task. This last statement is only valid at the atomic level.In the continuum, one can define a stimulus as the small expansion or shrinking volume. Hence, the error only appears in the ME theory but not in the macroscopic theory of Khachaturyan, which is widely used in phase field modelling for example.

  represents the image of f mat -sol net (r) in reciprocal space. On the other hand, the second term in (12) can't be analytically simplified in reciprocal space. But only elements where both c(r) and c(r´) don't vanish have to be accounted in the double sum over r and r´. These terms represent the bonding between two solute atoms. Since only few of such bonds exist in dilute alloys, the calculation of the second term of (12) is also in practice relatively simple and fast in real space.Once F ^(k ) is known, the lattice displacements u are obtained in reciprocal space as u ^(k ) = G ^(k ) F ^(k ) where G ^(k ) is the image of the Green tensor in reciprocal , the change of elastic energy ∆E

  d pq represents the equilibrium bond length between atoms of natures p and q.

  r 1 ). The calculations with three atoms (Tab. 1) show that only the new version of the ME theory provides a correct estimation (±3meV) of the configuration energies derived by MS relaxations. The Cook-de Fontaine version fails for one case, and the Khachaturyan version in several configurations.

in

  AlCu alloys nucleate during an MC simulation with the new ME theory. Details on this phenomenon are discussed in[START_REF] Jannot | Study of the effect of dislocations on precipitation kinetics by means of atomic simulations[END_REF].

  Meanwhile, the derived ME theory still possesses a serious inconvenience: all bonds have the same stiffness. The reason for this limitation relates to the definition of the force constants matrix. When atoms of two different types are present in the lattice, the large matrix D can be expressed as a combination of small local force constant matrices of size 3×3 like the following ones for the [110] first neighbour interaction in fcc an lattice:

Figure 1 .Figure 2 .Figure 3 .Figure 4 .Figure 5 .

 12345 Figure 1. Configuration change due to an exchange of position between a solute and a matrix atom

Figure 1 .

 1 Figure 1. Configuration change due to an exchange of position between a solute and a matrix atom

Figure 2 .Figure 3 .

 23 Figure 2. Decomposition of the transformation in three steps

Figure 4 .

 4 Figure 4. Configuration energies with two atoms

Figure 5 .

 5 Figure 5. Atomic configuration from an MC simulation in an Al3%at.Cu alloy using the new ME theory. only Cu-atoms are displayed.

  An analytic expression for the net force F(r) compatible with the four cases described above is required to work in reciprocal space. A useful formulation reads

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y
	3 about here -

  The change of bond stiffness also modifies the expression of the bond energies E ref in the reference lattice, where the specific stiffness has to be considered instead of a general κ value in expressions (15) to (17).Until now, it has also been assumed that the nature of the atomic interactions is purely harmonic. Such description of the bonding between atoms only allows the modelling of a material with a specific relation between elastic constants C 12 =C 44 =C 11 /2, which induces a specific anisotropy ratio χ=(C 11 -C 12 )/(2C 44 ) of 0.5.

	Performing the modifications of the model described above, the consideration
	of different bond stiffness is illustrated in an AlCu model alloy with stiffness
	κ AlAl =1.4eV/Å 2 (giving the elastic constants C 11 Al =112GPa and C 12 Al = C 44 Al =56GPa) and F κ CuCu =2.1eV/Å 2 ( C 11 Cu =187GPa and C 12 Cu = C 44 Cu =93GPa). As before, the straightforward o r MC-MS simulation is used as a benchmark. The results (Tab. 2) for a configuration
	P with two solute atoms and for different Al-Cu bond stiffness κ AlCu indicate that the
	e method described by equation (27) qualitatively reproduces the influence of a
	e r difference of bond stiffness. However, discrepancies of the order of 0.01eV are
	R observed for the first nearest neighbour interaction.
	-insert tab. 2 about here -	e
		v i e
		w
		O n l
		y
		0	(27)

  A new Microscopic Elasticity theory has been derived and validated by using slow, but unequivocal Molecular Statics simulations. The new version corrects a discrepancy found in earlier versions. Unlike the latter, the application of the present version in Monte-Carlo simulations of precipitation yields results that are in full agreement with theoretical expectations and with experimental findings in Al-Cu alloys. Moreover, the new version puts an emphasis on the atomic bonds instead of just differing atom sizes. This introduces new, independent elastic parameters (distinction between solute-solute bonds and solute-matrix ones) and thereby widens the scope in which linear elasticity can be used to describe and understand atomic behaviour.

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

Table captions Table 1 .

 captions1 Configuration energies with three atoms

Table 2 .

 2 Consideration of different bond stiffness: configuration energies with two solute atoms using the ME theory and MS relaxations

Table 1 .

 1 Configuration energies with three atoms MS relax. [eV] Error ME th. [eV] MS relax. [eV] Error ME th. [eV] MS relax. [eV] Error ME th. [eV]

		β AlCu = 1.4 eV/Å 2	β AlCu = 1.7 eV/Å 2	β AlCu = 2.1 eV/Å 2
	NN
	1	-0.066	0.005	-0.126	-0.006	-0.185	-0.025
	2	-0.028	0.000	-0.042	-0.001	-0.057	-0.003
	3	0.001	0.001	0.003	0.000	0.005	-0.001
	4	0.010	0.001	0.017	-0.001	0.026	-0.003
	5	-0.007	0.001	-0.011	0.001	-0.015	0.001
	6	0.001	0.001	0.002	0.001	0.004	0.001
	7	0.001	0.001	0.002	0.001	0.003	0.001
	8	-0.004	0.001	-0.006	0.001	-0.009	0.001

Table 2 .

 2 Consideration of different bond stiffness: configuration energies with two solute atoms using the ME theory and MS relaxations

Page 22 of 22 http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

Acknowledgements

This research was carried out under project number MC4.03167 in the framework of the Strategic Research program of the Netherlands Institute for Metals Research in the Netherlands (www.nimr.nl).


Anna Serra

D J Bacon

Interaction of a moving twin boundary with perfect dislocations and loops in an HCP metal

Keywords: atomistic simulation, crystal defects, crystalline interface, deformation twinning, dislocation interactions, topological theory of defects, twinning atomistic simulation, crystal defects, crystalline interface

Interaction of a moving 101-2 twin boundary with perfect dislocations and loops in an HCP metal

Introduction

Slip in the HCP metals is dominated by dislocations with Burgers vector b=1/3 that glide on either the (0001) basal plane or { } first-order prism planes.

However, the low crystal symmetry limits the number of independent slip systems available for glide of these perfect dislocations and, as a result, plastic deformation by twinning is common. A twin grows into the adjoining crystal by the diffusion-less motion of the interface between them and is therefore influenced, in principle, by crystal defects encountered by the moving interface. This interaction between a twin boundary and other defects is modelled in the present paper by atomic-scale computer simulation.

We consider a model for α-zirconium (Zr), which is an important metal for nuclear technology, for example. It twins by several different modes, of which twinning is the most common, i.e. the twin boundary is a plane in both crystals. Such twins are broad due to the high mobility of their boundaries [1].

Earlier computer simulation studies of the interaction between a twin boundary and crystal defects were restricted by the model size and rigid boundary conditions to cases where the twin interface was static and the defects were created sufficiently close to it to react with it spontaneously. This restriction has now been overcome by the development of a boundary condition technique that enables the twin interface to move freely. It allows one crystal in a twin-related pair to expand at the expense of the other by glide of twinning dislocations [2,3], unrestricted by the boundary conditions (see section 2).

The method has already been used to investigate the interaction of a twin boundary with single and small clusters of point defects [4], a situation relevant to neutron-irradiated Zr alloys. However, the defects considered in that paper were small (containing from 2 to 35 point defects) and did not have true dislocation structure.

There are two reasons for broadening the study in the present work to straight dislocations and larger defect clusters. First, the earlier simulations with rigid boundary conditions showed that a straight, perfect dislocation near a twin boundary may be absorbed by it spontaneously, thereby creating an interfacial step with dislocation character [5]. (Such an interfacial line defect is known as a 'disconnection' [START_REF] Pond | Dislocations in Solids[END_REF].) Under the application of a suitable shear strain, the disconnection may be a source for creation of twinning dislocations, which, by their motion in the interface, cause it to move.* In simulation of some dislocation-interface reactions, the product disconnection was seen to glide unchanged with the interface, although the Burgers vector did not lie in the plane of the step riser of the defect . However, due to the small number of atoms that could be simulated in the earlier studies, motion of the disconnection was restricted by pinning of the twin interface at the rigid boundaries of the model. This difficulty is avoided in the present work and so the process can be studied in full. The second reason for re-examining the interaction of a twin boundary with clusters of point defects is that clusters containing up to only 35 point defects, as in [4], are barely resolvable by transmission electron microscopy (TEM), whereas the clusters observed by TEM in irradiated Zr can be sufficiently large to have distinct dislocation character [START_REF] Griffiths | [END_REF]. Thus, by simulating dislocation loops formed by hundreds of point defects, it should be possible to see if the mechanisms revealed earlier and here for straight, infinitely-long dislocations occur for the more complicated geometry when a twin boundary meets the side of a dislocation loop. In this context, we note that although Higashiguchi et al. [8] observed twins in neutron-irradiated Zr deformed at either 77 or 300K, Onimus et al. [9] found that twinning did not play a role in the deformation of Zr irradiated with neutrons at 350 o C. Instead, dislocation channelling occurred in narrow bands, the orientation of which depended on the conditions of loading.

Before proceeding to the model, a brief introduction to the crystallography of the defects studied is required. The admissible interfacial disconnections separating regions of identical boundary structure are best identified using the dichromatic pattern associated with the twin. It is shown in Fig. 1 and consists of the lattice sites of two crystals in twin-related orientation, projected with the sites of a lattice plane of one crystal coincident with those of a plane of the other. The two crystals are drawn as white (λ) and black (µ) in Fig. 1. In addition to the crystal coordinate system it is convenient to define an interfacial coordinate system with x parallel to the common tilt axis z parallel to the normal to the interface and y parallel to the twinning disconnections is obtained from the difference in the same reference frame between translation vectors t(λ) and t(µ) of the two lattices, defined in the following sense:

b = t(λ) -P t(µ), ( 1
)
where P is the transformation matrix relating the coordinate frames of the two crystals.

The sign conventions used here are chosen so that the positive line sense of the disconnection is taken to be out of the paper (+x direction) and b is given by the RH/FS rule, e.g. [START_REF] Hull | Introduction to Dislocations[END_REF]. Furthermore, the sense of step 'up' or 'down' refers to a boundary being traversed from left to right, i.e. in the +y direction. Burgers vectors are written as b ±m/±n , where the magnitude of the subscripts m and n indicates the number of lattice planes, spacing d, at the step in the upper and lower crystals, respectively, and the subscript sign indicates the step sense. Thus, a twinning dislocation has Burgers vector b 2/2 (step up) or b -2/-2 (step down). The latter is shown in Fig. 1 and is the difference between t λ = λ defining a step down by two planes on the lower surface of the λ crystal, and t µ = µ defining a step down by two planes on the upper surface of the µ crystal. The same notation is adopted for crystal dislocations, i.e. b of a dislocation of the type 1/3 at 60 o to the x axis in the λ crystal is labelled b 1/0 in Fig. 1 and is defined by t λ that steps up by one lattice plane in λ, whereas b -1/0 corresponds to t λ that steps down by one plane. (With the two-subscript notation for Burgers vectors, the magnitude and sign of each subscript is conserved in dislocation and disconnection reactions, as will be seen in section 3.)

Model

The computational model had about one million atoms (dimensions L x , L y , L z = 23.4, 35.5, 28.9 nm) and consisted of two half-crystals in twin orientation separated by either a planar twin boundary or one containing a twinning disconnection in edge orientation with line direction parallel to the x-axis λ and Burgers vector b -2/-2 = α : α is approximately 0.0645 for the lattice parameter ratio c/a = 1.5952 (see [1] for details).

This bi-crystal model had periodic boundary conditions in the interfacial glide plane x-y of the twinning disconnection and fixed outer boundaries in the form of rigid blocks of defect has its extra half-plane in the upper (λ) crystal and glides to the right (+y) under positive applied shear strain.) The atomic system was relaxed to minimise the potential energy after each shear displacement using the interatomic potential for Zr from [START_REF] Ackland | [END_REF] that we refer to as A95. The applied resolved shear stress for each level of strain was calculated from the resultant force on the upper rigid block, as described in [12]. To check the stability of some configurations, molecular dynamic simulations for T = 100K

were performed using two different interatomic potentials, namely A95 and the potential described in [13], which we denote as M07. In these cases, strain was applied at a constant rate of 25x10 6 s -1 by applying a constant velocity to the upper rigid block of atoms. The principal difference between the two potentials arises from the value of the intrinsic stacking fault energy on the basal and first-order prism planes. The fault energy is relatively low on the basal plane in the A95 model, resulting in a preference for dissociation into two Shockley partials. Dissociation also occurs in the M07 model, but the partial spacing is much smaller [14]. shown in [15] that the screw dislocation gliding under stress in the λ crystal is either transmitted through the twin boundary into the µ crystal or is reflected back into the λ crystal with a change of glide plane, e.g. basal-to-prism or vice versa. Decomposition into products that change the structure of the boundary does not occur. Second, it was found in [5] that a 60°-mixed dislocation decomposes when entering the interface, thereby changing its structure and creating interfacial disconnections. More detail follows below. Note that perfect dislocation loops have the same 1/3 λ Burgers vectors, and for these defects it is necessary to consider all three orientations of b, as will be seen in section 4.

Interaction of a straight dislocation with the twin boundary

The interaction between the twin boundary and a 60° crystal dislocation gliding on the basal plane depends on the sign of the dislocation and the order of the two Shockley partials relative to the boundary. (If the basal glide plane of the crystal dislocation lies between planes A and B, say, in the ...ABAB... stacking sequence, the order of the Shockley partials is reversed from that when the glide plane lies between B and A.) The partials have Burgers vectors b 30 and b 90 at 30 0 and 90 0 to the line direction, respectively. The interaction of 60° dislocations of opposite line sense with the twin interface was reported in [16] The dichromatic pattern in Fig. 1 shows a set of possible reactions represented by the following Burgers vector balance:

b 1/0 = b -(2n-1)/-2n + nb 2/2 . ( 2
)
The magnitude of b -(2n-1)/-2n depends on the height of the step (nd) and the c/a ratio: when c/a just less than the ideal value (8/3) 1/2 , as in Zr, the smallest magnitudes correspond to n = 3 and 4 (see Fig. 1). The dislocation is attracted to, and absorbed spontaneously by, the twin boundary in this reaction when less than about 8 nm from it. In [5] it decomposed into a b -5/-6 disconnection line and three b 2/2 twinning disconnections (Fig.

-1/0 = (b 4/4 + b L ) + b T + 2b -2/-2 , (3)
where b L and b T are the Burgers vectors of the leading and trailing partials respectively.

The outcome is represented schematically in Fig. 2(b).

The reaction of eq. ( 3) and Fig. 2(b) is the same as that found with the smaller atomic model in [16] when the leading partial has b L = b 90 . However, it was reported in [16] that when the leading partial has b L = b 30 the following reaction occurs: The difference in the mechanisms described in the previous two paragraphs is believed to arise from the atomic structure of the twin boundary. As discussed in [15], rows of atoms in the boundary alternate in spacing, as a result of which the hydrostatic stress alternates between positive and negative along a row. Consequently, the atomic configuration of the interstitial atom at the lowest corner of a loop is not the same for the two situations described above. However, in neither case are any point defects removed from or added to the loop and it is simply transformed from the λ crystal to the µ crystal with the same shape, size and Burgers vector. twinning dislocation is less than that of the interstitial loop and τ c is 40MPa.

Discussion and Conclusions

The development of a method that allows an interface containing a disconnection line, i.e. a defect with step as well as dislocation character, to be simulated with periodic boundary conditions in the interface plane has allowed us to investigate the interaction of a moving twin boundary with two types of crystal defect, namely a straight, infinitely long dislocation line in a basal plane and a periodic row of perfect dislocation loops.

The interface considered is that of the twin, which is important and common in the deformation of the HCP metals. The interaction between straight crystal dislocations of the 1/3 slip system and the twin interface have been modelled previously, but by simulation where the number of atoms in the model was small and fixed boundary conditions were employed, as a consequence of which the interface was unable to be translated through the model by glide of twinning dislocations. With this restriction overcome, it has been possible to revisit the

2(a)

). The current simulations confirm the earlier finding [5] that this disconnection moves with the twin boundary along its riser (the basal plane of one of the crystals) because twinning dislocations can move past it by glide and self-compensated climb.

Furthermore, this disconnection acts as a source of twinning dislocations when stress is applied in the absence of pre-existing twinning dislocations. Thus, this decomposition product provides a simple mechanism for twin growth/shrinkage, as predicted in [5].

The decomposition product of a 60° dislocation of opposite sign (b = b -1/0 ) with its extra half-plane directed away from the interface is more complex (Fig. 2

  Burgers vectors of admissible interface

  ±z directions[2,3]. With this construction, when a twinning disconnection moving in the ±y direction crosses one periodic boundary, it re-enters at the opposite side of the simulation box and, due to the step character of the disconnection, the twin interface is displaced in the z direction by the step height ±2d. Crystal dislocations or dislocation loops were introduced on either side of the twin boundary. To assist understanding of the reactions and mechanisms found, atoms with irregular coordination number were recorded, as will be shown in the figures that follow. Most of the modelling used molecular statics simulations, i.e. temperature T = 0K. In these, the b -2/-2 disconnection was induced to glide by applying shear strain e yz in increments of 10 -5 via displacement of the upper rigid block relative to the lower one. (The b -2/-2

3. 1 .

 1 Dislocations considered A straight perfect dislocation with line direction along the x axis λ can be pure screw, i.e. b = 1/3 λ , or have 60 °-mixed character, i.e. b = ±1/3 λ or b = ±1/3 λ . The latter vectors are represented by b ±1/0 in the dichromatic pattern of Fig. 1. However, only the 60 ° cases are considered here, for two reasons. First, it was

  Fig.2is a schematic summary to be discussed below of the results of the earlier simulations for the model with a small number of atoms and fixed boundary conditions[16]. There are three different possibilities for the reaction, depending on the sign of the dislocation line and the stacking sequence of the basal plane on which it glides. In the Burgers vector notation used, b of the reacting crystal dislocation is either b 1/0 or b -1/0 , as indicated. The outcomes are as follows.

2

  (a)). In the larger, periodic model of the present work, it formed a b -7/-8 defect and four b 2/2 defects. This outcome was found to be independent of the leading/trailing order of the partials of the crystal dislocation. If the boundary initially contains a b -2/-2 disconnection, it annihilates with one of the b 2/2 defects, leaving only three. The periodicity of the simulated box allows each b 2/2 twinning dislocation to move freely under the influence of the others and they spread equidistantly along the interface.Under applied shear strain e yz > 0, the b 2/2 disconnections glide to the left (-y direction).The critical resolved shear stress, τ c , necessary for the upper block to move a distance |b 2/2 |, i.e. for the total distance moved by the four (or three) defects to equal L y , is value represents the obstacle stress when a b 2/2 twinning dislocation encounters a b -7/-8 disconnection and overcomes it by conservative climb and subsequent detachment on the other side. (It was shown in[16] that climb of a b ±2/±2 defect from one side of a b -7/-8 disconnection to the other is compensated by climb of the latter along its riser plane as the twin boundary translates by ±2d.) Furthermore, it has been shown that if the b -7/-8 disconnection exists in the interface in the absence of pre-existing twinning dislocations, it can act as a source of dipoles containing pairs of equal and opposite b ±2/±2 defects when τ reaches a value of 110MPa[14].Thus, the disconnections such as b ±7/±8 and b ±5/±6 can be created by the decomposition of crystal dislocations at the twin boundary and are effective sources of b ±2/±2 twinning dislocations. They move with the twin boundary by conservative climb along their riser[5,16], thereby supplying the twinning disconnections necessary to keep the boundary moving. 3.3. Dislocation with b = b -1/0 The reaction is more complex since although three disconnections form in the boundary, one with b 4/4 and two with b -2/-2 , the b -1/0 dislocation is not fully absorbed. Instead, only the leading partial of the b -1/0 line is attached to the b 4/4 disconnection. The trailing partial stands off from the boundary, separated from it by an intrinsic basal stacking fault. The leading Shockley partial stabilizes the b 4/4 disconnection because this defect on its own is not stable and decomposes into two b 2/2 twinning dislocations, in accordance with Frank's rule. We can therefore describe the reaction as b

b 4 . 3 . 4 . 2 .

 4342 Fig. 4 shows visualisations of the loop and boundary containing a b -2/-2 disconnection at different stages of the process. Fig. 4(a) is a near projection as the b -2/-2 defect approaches the loop from the left and the twin boundary is about to meet the lower corner of the loop. The conversion of the loop from λ to µ is half-complete in Fig. 4(b), where the x-axis projection shows clearly the transformation of the loop sides from basal to prism and vice versa. The loop acts as an obstacle to the repeated glide of the b - 2/-2 twinning disconnection, as can be seen in Fig. 4(c) by the z-axis projection of the configuration shown in (b). The b -2/-2 dislocation is seen to be bent strongly as it breaks away from the loop. The critical stress for this is 60MPa.

4 . 3 .

 43 Fig. 6(b). The other three sides of the loop remain separate from the twin boundary. As the boundary moves upwards into the λ crystal by glide of the pre-existing b -2/-2 twinning dislocation, the loop moves upwards along its glide prism with side BS1 incorporated in the twin boundary. The repeated emission of the two b -2/-2 defects assists the motion of the twin boundary and the maximum τ for the process is only 20MPa. When the centre of gravity of the loop is shifted by ½[0001] λ so that the 90° partial leads the entry of BS1 into the twin boundary, all four sides of the loop are absorbed into the twin boundary, as shown in Fig. 6(c). As the twin boundary moves upwards under applied stress, it drags the absorbed interstitial cluster along the 60° direction of the glide prism of the original loop. This process provides a strong obstacle to the glide of the b -2/-2 twinning dislocation and τ c is 55MPa. MD simulations of T = 100K showed the same features.

  straight dislocation in the twin interface and now investigate its effect on boundary mobility. Furthermore, by simulating the interaction between a moving twin boundary and either interstitial or vacancy clusters that are large enough to have clear dislocation character, it has been possible to interpret some of the effects observed in terms of the results for straight dislocations. We have considered straight, 1/3 crystal dislocations lying parallel to the interface and perpendicular to the direction of twinning shear. The screw dislocation parallel to the interface does not decompose when it meets the moving interface and has no discernible effect on its motion. The 60° mixed dislocation is attracted by the boundary in all situations, but the form of the decomposition product and its effect on twin boundary motion depend on the sign of the dislocation (b ±1/0 in the notation here). The dislocation with b = b 1/0 has its extra half-plane intersecting the interface and fully decomposes in it to form a b -7/-8 disconnection and four b 2/2 twinning dislocations, independently of the nature of the leading/trailing Shockley partials (Fig.

  (b)). Three disconnections form in the boundary, one with b 4/4 and two with b -2/-2 , but the b -1/0 dislocation is not fully absorbed. Instead, only the leading partial of the b -1/0 line is held in the core of the b 4/4 disconnection (Fig. 2(b)). The (b 4/4 + b L ) complex thus formed is stable under positive and negative shear strain e yz and acts as an obstacle to the glide of b ±2/±2 twinning dislocations. The critical resolved shear stress for twinning dislocations to glide over the interface rises to 70-80MPa. The b 4/4 disconnection resulting from the partial absorption of the b -1/0 crystal dislocation has a riser lying perpendicular to the basal planes of the upper (λ) crystal.

Fig. 1 .Fig. 2 .

 12 Figure captions

Fig. 3 .

 3 Fig. 3. Near-x-axis visualisation of the twin boundary when first a 60°-mixed b -1/0 crystal dislocation is absorbed to form a (b 4/4 + b L ) complex and then a second dislocation of opposite sign (b 1/0 ) is absorbed to form a b -5/-6 disconnection. The trailing 30° partial of the first dislocation is seen to be standing away from the interface. The nature of the defects involved is indicated by the usual dislocation symbols. Atoms in the interface are visualised and coloured by their coordination number.

Fig. 4 .

 4 Fig. 4. (a, b) Visualisations of interface and dislocation atoms showing two stages when a twin boundary meets and passes through an interstitial dislocation loop with b parallel to the x-axis. (c) The configuration in (b) in z-axis projection: the curved shape of the twinning dislocation gliding from left to right demonstrates the obstacle strength the loop provides to twinning. (d) Alternative to (a) when the order of the Shockley partials making up the basal-plane sides of the loop is reversed. The Shockley partials are seen in (c) and (d).

Fig. 5 .

 5 Fig. 5. (a,b) Visualisation of interface and dislocation atoms showing two stages in the interaction between a twin boundary and a vacancy dislocation loop with b parallel to the x-axis. The loop is seen to be smaller in (b) after the interaction due to the detachment of vacancies on the left and right. (c) Change in potential energy of the

Fig. 6 .

 6 Fig. 6. Visualisations of the interaction between the twin boundary and interstitial loops with b at 60° to the x axis. (a) The loop in the model with the M07 interatomic potential is pushed on its glide prism ahead of the boundary. In the model with the A95 potential, either (b) one side of the loop is trapped in the boundary or (c) all sides are absorbed, depending on the order of the basal-plane Shockley partials. In both cases, the interstitials are transported along the glide prism of the original loop as the interface moves up or down.

Fig. 7 .

 7 Fig. 7. (a) Near-x-axis visualisation after absorption of a vacancy loop with b at 60° to the x axis in the twin boundary. A b -2/-2 twinning dislocation is approaching from the left. (b) z-axis visualisation showing the asymmetric shape of the b -2/-2 twinning dislocation as it passes through the periodic row of loops. Two model cells are set sideby-side to make the dislocation shape clearer.

  FIGURE 1

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

Acknowledgements

Discussions with Dr. Yu N Osetsky are acknowledged. A.S. was supported by the Spanish Ministry of Science and Innovation (FIS2006-12436-C02-02) and the Catalan Government (CIRIT 2005 SGR -00779). The computing was carried out in CESCA (www.cesca.es).


Mark Robert

Keywords: computer simulation, dislocation theory user supplied): core structure, bcc iron, Frenkel-Kontorova

published or not. The documents may come

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters