Anthony R Bunsell
email: anthony.bunsell@mines-paristech.fr

Alain Thionnet
email: alain.thionnet@mines-paristech.fr

Life prediction for carbon fibre filament wound composite structures

Keywords: MinesParistech, Materials carbon fibres, composite materials, computer simulation, experimental, finite-element modelling, hydrogen, Monte-Carlo, multiscale modelling pressure vessels, failure, lifetime prediction Life prediction, composite filament wound structures, fibre failures, viscoelasticity, debonding, pressure vessels

The processes governing the failure of filament wound composite structures have been examined. It is shown that the fibres controlling the failure of such a structure, when it is internally pressurised, can be considered to be subjected only to tensile loads. A multi-scale model has been developed which considers the effects on the scale of the elastic fibres and includes the effects of the viscoelastic matrix as well as debonding around fibre failures. The intact fibres neighbouring fibre breaks are subjected to an increase in stress and a higher probability of failure than elsewhere in the composite. During a monotonic failure test the initially random fibre failures are seen to begin to coalesce in a way governed by the stochastic nature of the fibre breaks and this eventually leads to failure. Under prolonged loading the relaxation of the matrix around fibre breaks causes the overloads in the neighbouring fibres to evolve and to induce delayed fibre breaks which can eventually lead to instability in the structure. The model allows these processes to be taken into account while calculating the behaviour of the whole composite structure.

Introduction

High performance carbon fibres were first produced commercially in the second half of the 1960s in the UK [START_REF] Phillips | Carbon fibre reinforced plastics. The first fifteen years[END_REF]. This was the result of work at the Royal Aircraft Establishment, Farnborough. The carbon fibres were made from polyacrylonitrile precursors and over a period of forty years their composites have become established as one of the most important materials used in advanced aerospace structures. The presently emerging generation of civil aircraft is being made with half their weight made up of carbon fibre reinforced epoxy resin. Today, however, the largest markets for these materials are no longer confined to the aerospace industry but wider industrial sectors are finding uses for these remarkable materials. This success has to no small degree been aided by researchers, such as Anthony Kelly and his colleagues, who have allowed an in-depth understanding of the physical mechanisms governing the behaviour of composites [START_REF] Kelly | Strong Solids[END_REF].

One of the major challenges, however, is in the prediction of long-term behaviour and reliability of composite structures. This is particularly important as advanced composites are often used in structures the failure of which would be catastrophic. It is only through extending the legacy that an earlier generation of researchers has left to us that techniques based on a true understanding of damage processes in advanced composites can be developed and reliability in assessing composite structures as well as optimising designs can be achieved. One of the most challenging applications is in pressure vessels for the storage of gas at high pressures. Pressure vessels have been in use since the beginning of the industrial revolution and a large body of empirical knowledge has been accumulated over the years, often at the expense of human lives lost in explosions. Traditionally these pressure vessels have been made of steel and the experience which has been obtained through accidents is now largely understood in terms of metallurgical mechanisms determining crack propagation. The standard procedure for assessing the reliability of a metal pressure vessel is the hydraulic test in which the vessel is subjected to one and a half times the maximum in-service pressure. If it survives, it is deemed suitable for continued use in service. This is instinctively satisfying and is justified, as at the higher pressure, any incipient crack will either propagate, due to the increased stress in the vicinity of the crack or a plastic zone will develop around the crack tip, which, at lower pressures, will impede further crack propagation. However steel is heavy and such pressure vessels are unsuitable for the storage of gases used as a fuel for vehicles. The finite availability of oil is leading to increasing use of natural gas for such applications and this can be seen as a forerunner of the use of hydrogen. The difficulty is that these gases have to be pressurised to very high pressures for although hydrogen has three times the energy content of petrol, for the same mass, it is much lighter and has to be compressed to 70MPa (700 bars) to become an efficient and useable alternative fuel for the future.

Advanced composites, usually made of high strength carbon fibres in an epoxy matrix, used for pressure vessels provide an obvious solution to the needs for transport and they are being increasingly used for onboard applications for vehicles such as buses, utility vehicles and private cars [START_REF] Bunsell | Composite pressure vessels supply the answer to transport problems[END_REF]. Figure 1 shows carbon fibre composites pressure vessels on the roof of a bus running on natural gas. Composites are ideally suited for this application. This is due to their light weight, high strength and rigidity. They are presently used for storing natural gas at pressures usually up to 20MPa and are actively being investigated for the storage of hydrogen at pressures of 70MPa. However, legislation and common sense require that pressure vessels be tested periodically to ensure their safe continuing use in service. At present, there are no proof testing techniques or in-service reliability assessments techniques which are mentioned in standards that are suitable or based on the failure processes known to control lifetimes of composite structures. The standards which are described are based on the behaviour of metals which fail by processes which do not occur in composites. Composite materials do not fail by simple crack propagation and in these materials, plastic deformation is not a dominant mechanism determining fracture toughness. A hydraulic test on a composite pressure vessel has only one clear outcome, which is that fibres, which ordinarily would not have broken, fail during the test and the vessel is closer to failure than it was before the test, as is illustrated in Figure 2. Visual inspection is not an alternative technique as damage, due to, say an impact, produces delamination on the reverse side of the composite and this is not available for inspection. There are techniques, such as ultrasonic inspection and others, that permit localised damage due to impact to be identified, however, the pressure vessels will be used for decades, possibly twenty or thirty years and the general assumption that carbon fibre reinforced epoxy resin composites do not age, due to the combined effects of pressure and time is easily shown to be a fallacy. This paper will demonstrate how such composites evolve with pressure and time and propose means to determine residual lifetimes.

Filament wound composite pressure vessels, of the type considered in this study, are made of fine carbon fibres, typically T700 fibres or equivalent, made from polyacrylonitrile precursors. The fibres have diameters of 7µm, they are perfectly elastic with a Young's modulus of 235 GPa and in common with all populations of fibres, show considerable scatter in their strengths, giving a Weibull modulus of typically 6. The fibres are placed on geodesic paths around a mandrel, which later serves as a liner to ensure gas tightness. This means that when the vessel is pressurised the fibres are subjected only to tensile forces and at the level of the fibres an analogy can be made with the behaviour of a unidirectional composite loaded in the fibre direction. As with a unidirectional composite, failure of the material means that the fibres have to break. In addition the burst strengths of the vessels are controlled by the winding of a fretting layer of fibres over the cylindrical part of the pressure vessel, at 90° to its major axis. In the case of a carbon fibre composite, the fibres support around 99% of the load. The failure of carbon fibres is therefore clearly the mechanism which determines the failure of the composite pressure vessel. It is not crack propagation, as in a metal pressure vessel, as is assumed in the standards.

The consequences of failure of fibres in a composite material have been examined in a number of studies, beginning with Cox who developed a 2-D analytical analysis of load sharing between elastic fibres embedded in an elastic matrix (4). This study revealed the fundamental mechanism controlling composite behaviour to be load transfer from the matrix, which undergoes shear near the fibre, to which it is well bonded, to the fibre which experiences tensile force loads which increase from its end. By writing equations of equilibrium between the shear forces in the matrix and the tensile forces in the fibre, Cox was able to provide analytical solutions for the stress states in and around short or discontinuous fibres in a composite. Subsequent studies considered the effect of perfectly plastic matrices and were shown to simplify the analysis [START_REF] Kelly | Strong Solids[END_REF]. Other studies considered the effects of fibre failure on the other unbroken fibres in a composite [START_REF] Rosen | Tensile failure of fibrous composites[END_REF][START_REF] Zweben | Tensile failure of fiber composites[END_REF]. It was shown that the effect of a fibre break in a composite was confined to a small volume of material around the break and that only the fibres immediately neighbouring the fibre break experienced any change of stress. With the availability of computers and increasing computational capacity this situation has been the subject of many subsequent two and three dimensional numerical studies, including a few which consider the effects of a viscoelastic matrix [START_REF] Nedele | Three dimensional finite analysis of the stress concentration at a single fibre break[END_REF][START_REF] Landis | Micromechanical simulation of the failure of fiber reinforced composites[END_REF][START_REF] Lifschitz | Time-dependent longitudinal strength of unidirectional fibrous composites[END_REF].

Studies in our group have examined this problem with the aim of improving the models and applying the results to various filament wound structures, including pressure vessels. In order to tackle this problem an immediate practical problem has to be considered. In a unidirectional cfrp specimen with a cross section of 20mm x 1mm there are about 136,000 In any cross section of a CNG tank there are around 1.25 million fibres. It is clear that any attempt to model the behaviour of a complete structure with a discretisation at the level of the individual fibres would encounter unacceptable computational times. Nevertheless the mechanisms which govern failure occur at the scale of the individual fibres. The FE 2 multiscale approach which has been adopted can be seen to share points in common described in references (10 and 11).

Failure process on the microscopic scale

Studies of the effect of the stochastic nature of fibre failures in composites have underlined the controlling nature of this type of damage [START_REF] Argon | Fracture: strength and toughness mechanisms[END_REF][START_REF] Curtin | Stochastic damage evolution and failure in fiber-reinforced composites[END_REF][START_REF] Phoenix | Statistical strength theory for fibrous composite materials[END_REF].

The micromechanics study considers a three-dimensional analysis of a representative volume element (RVE) comprising all the processes to be considered. These are fibre strength and its stochastic nature, fibre tensile elastic modulus and the properties of the matrix. Several levels of refinement are considered so that at the simplest level both fibres and matrix are considered to be elastic. Interfacial debonding between the broken fibre and the matrix is considered at a higher level of refinement. Finally, at the highest level of sophistication, the matrix is considered to be viscoelastic, so that its properties are time dependent, as are the effects of load transfer between neighbouring fibres, which are due to shear of the matrix.

In this study both the fibres and the matrix are considered to be homogeneous continua. The orthorhombic spatial frame used at the microscopic level is R loc = (O, y , x , z ) for which (x, y, z) describe the coordinates of a point M. The vector z is aligned parallel to the fibres. Earlier studies have allowed the Representative Volume Element (RVE) to be defined for the undamaged composite [START_REF] Blassiau | Micromechanisms of load transfer in a unidirectional carbon-epoxy composite due to fibre failures, Part 1: Micromechanisms and 3D analysis of load transfer, the elastic case[END_REF][START_REF] Blassiau | Micromechanisms of load transfer in a unidirectional carbon-epoxy composite due to fibre failures, Part 2: Influence of viscoelastic and plastic matrices on the mechanism of load transfer[END_REF][START_REF] Blassiau | Micromechanisms of load transfer in a unidirectional carbon-epoxy composite due to fibre failures, Part 3[END_REF][START_REF] Baxevanakis | Comportement statistique à rupture des composites stratifiés[END_REF]. The original 2D model showed that a constant average unidirectional composite failure stress could be calculated using a RVE consisting of six parallel fibres of 8mm length [START_REF] Baxevanakis | Comportement statistique à rupture des composites stratifiés[END_REF]. This model has been extended to the 3D case by first considering RVEs consisting of six by six square arrays of fibres with lengths of 8mm. This has been refined for geometrical considerations to arrays consisting of 32 fibres. The degrees of damage were modelled at six levels of failure in the RVE, written as E-N, in which N represents the state of damage in the group of thirty two fibres. The values of N represent the fraction of single broken fibres to unbroken fibres and were N=32, 16, 8, 4 and 2, so that, for example, when N equals 8, one in eight fibres is considered broken. The associated periodic cell was given the abbreviation C-N. The failed state E-0 is represented by C-0. The state which represents isolated fibre breaks, so that no interaction occurs between breaks, is described by the damage state E-32 and gives unvarying failure stress for the composite. This is the justification for choosing the RVE consisting of 32 unbroken fibres as representing three-dimensional state of the undamaged composite structure.

It follows that the composite is considered to be periodic with a hexagonal fibre arrangement in the ( y , x ) plane, the periodic cell, which is representative of the RVE and called CS32, consisting of cells of 32 fibres. Its geometry is considered to be a parallelepiped with the ---z axis parallel to the fibres and the section is a square with sides of length c (Fig. 3). From reference 15, the length in the z direction was taken to be between the planes z = 0 and z = L = 8mm. The reference origin O is the geometrical centre of the section in the plane z = 0. At this scale, at point M, the stress tensor is described as σ and the strain as ε. The fibre is always considered to be linearly elastic, homogeneous and isotropic. In contrast, two types of behaviour have been considered for the matrix. The most reductive behaviour considered the matrix to be linearly elastic, homogeneous and isotropic. A more realistic model considers the matrix to be linearly viscoelastic with behaviour described by the following equation:

ε ε ε ε γ γ γ γ + + + + ε ε ε ε = = = = γσ γσ γσ γσ + + + + σ σ σ σ ∞ ∞ ∞ ∞ C C . 0 .
in which γ is the fourth order relaxation tensor; C 0 is the initial rigidity tensor and C ∞ the final rigidity tensor after an infinite time. As the matrix is istoropic, each of the tensors is isotropic. In this way the notion of time can be included in the model of the behaviour of the composite and permits delayed composite failure to be studied.

Earlier studies have shown that this approach accurately describes what happens in the vicinity of a fibre break [START_REF] Blassiau | Micromechanisms of load transfer in a unidirectional carbon-epoxy composite due to fibre failures, Part 1: Micromechanisms and 3D analysis of load transfer, the elastic case[END_REF][START_REF] Blassiau | Micromechanisms of load transfer in a unidirectional carbon-epoxy composite due to fibre failures, Part 2: Influence of viscoelastic and plastic matrices on the mechanism of load transfer[END_REF][START_REF] Blassiau | Micromechanisms of load transfer in a unidirectional carbon-epoxy composite due to fibre failures, Part 3[END_REF] and can account for:

• the random nature of fibre breaks and the points of failure along the fibre length; • the number of breaks in a RVE at five levels of damage going from the undamaged to the failed state. These states are : C32, C16, C8, C4, C2 (Fig. 2) and contain respectively 1, 2, 4, 8, 16 broken fibres amongst the 32 present in the CS32 cell; • the load transfer to and along the intact fibres due to the fibre failures.

Debonding is modelled by the undoing of the finite element boundary nodes of the matrix and the broken fibre. Outside this zone the fibre-matrix bond is considered to be perfect. The debonded length is considered to be twice the distance from the broken end of the fibre to the point where the fibre and matrix are again adhered. In this way the model accounts for :

• the effects of interfacial debonding around the fibre break on load transfer to intact fibres. Ten debonding lengths have been considered from 3.5 to 35µm in steps of 3.5µm; From the above considerations the model considers • the effects of the viscoelastic nature of the matrix, which is considered to be linear, on load transfer to neighbouring intact fibres. In addition the present study also takes into account the fibre volume fraction.

The goal of these calculations is to simulate the failure of fibres and to observe the evolution of the axial stress in neighbouring intact fibres as a function of the damage state, time, distance from the plane in which the break lies and as a function of any debonding of the fibre-matrix interface. In order to do this we define a coefficient of load transfer, k r . This is obtained as the mean of two axial stresses, in the fibre considered, calculated between two consecutive normal sections of the fibre in the FE mesh compared to the same quantity in the undamaged state.

The coefficient of load transfer along the lengths of the fibres is given by: (

)dxdydz

z y, x, , V O, t 0, d CS32, dxdydz z y x , t d Z , V t, d, C, f f r k = = = = = = = = = = = = ∫ ∫ ∫ ∫ ∫ ∫ ∫ ∫ σ σ σ σ ∫ ∫ ∫ ∫ σ σ σ σ ∫ ∫ ∫ ∫ + + + + + + + + Z Z S S Z Z i i F F i i zz f zz V C 1 1 , , , , ,
In which C is the cell representing the state of damage considered d is the debonded length t is the time after fibre failure V f is the fibre volume fraction z is the coordinate along the fibre from the plane of failure (z=0) Z i+1 and Z i are the abscissa of the plane sections between which k r is calculated

2 z z i i Z + + + + = = = = + + + +1
S F indicates the cross section of the fibre considered x, y are the coordinates of the plane section of the cell σ zz is the axial stress in the fibre considered It has been verified that the discretisation employed is that which results in convergence. The mesh is so fine that the stresses are nearly constant in each element.

Whether or not there is debonding at the fibre/matrix interface or if the matrix is viscoelastic, or not, the definition for k r remains valid. However, in order to separate the different effects we can write

υ υ υ υ + + + + + + + + = = = = r d r r r K K K k
where r K is the load transfer exclusively due to a fibre break; if this effect is neglected.

A data base at the microscopic scale

Comparisons have been made with the self consistent model in which one broken fibre is considered embedded in a homogeneous medium having the macroscopic properties of the intact composite [START_REF] Nedele | Three dimensional finite analysis of the stress concentration at a single fibre break[END_REF]. The results obtained by the present model, for the case in which interfacial bonding is maintained after failure, are in good agreement with such self consistent models. However, such models inevitably underestimate the stresses on intact fibres near the plane of failure as the homogeneous equivalent material is more rigid than the matrix material [START_REF] Blassiau | Damage accumulation processes and life prediction in unidirectional composites[END_REF]. This leads to difficulties in understanding the subsequent behaviour of these fibres in a model that takes the effects of time into account. These difficulties are removed in the present analysis which, in the case of a viscoelastic matrix and with data on the stochastic nature of fibre strength, allows the effects of matrix relaxation around fibre breaks to be evaluated. The result is seen as an increase in time, over increasing lengths, of the stresses on intact fibres neighbouring a fibre break. This is the cause of delayed fibre breaks and the evolving accumulation of damage in the composite.

If only one cell is considered and a total debonded length of 70 µm taken, as the effect of the fibre break become insignificant beyond this length: the calculation time to obtain geometrical convergence for a periodic homogenisation requires approximately a discretisation for 250,000 nodes and two hours of calculation time for an all elastic model. If this approach is applied to a macroscopic structure, for which the number of Gauss points is in the thousands, the calculation time would be unacceptably long. A less time consuming approach is therefore necessary.

The load transfer function has been smoothed and only the load transfer in the fibre nearest to the fibre break calculated. In order to separate the effects of interfacial debonding and those due to the viscoelastic shear of the matrix, we put the multiscale FE 2 model. Using the overall stress applied on the macroscopic scale, the duration of the load and the state of damage, this data base allows, by linearity, to know the stress state in the components at the microscopic scale. The tensile stress in the fibres can therefore be computed, without resolving the problem at the microscopic scale by the finite element method and the evolution of the damage, consisting of increasing numbers of fibre breaks in the RVE can be evaluated. It is assumed that the failure of a fibre induces a debond length of 35µm which has been shown to result in the greatest stress concentration in neighbouring intact fibres.

In order to compare the relative influences of different physical phenomena, for different fibre volume fractions, six levels of refinement in the model must be considered : the first taking into account only fibre breaks; the second adds the effect of load transfer, through the matrix on intact fibres neighbouring the fibre break; the third includes the effects of fibre/matrix debonding from the point of failure and then the previous three degrees of sophistication are again considered but the viscoelastic nature of the matrix is included.

The use of parallel computation to describe the macroscopic scale

The carbon fibres used in the specimens and in the pressure vessels had radii of 3.5µm so that the dimensions of the CS32 cell, representing a volume fraction 64%, were 0.044mm x 0.044mm x 8mm. The volume of the CS32 cell was V 64% =0.0155mm 3 .

Calculations carried out for a laboratory plate specimen have shown that the mesh necessary for convergence can be carried out with finite elements (of the c3d8 type) containing 4 RVEs.

It can be seen that if the cylindrical part of the pressure vessel is 500mm in length, has an internal radius of 100mm and a composite and wall thickness of 10mm, assuming a fibre volume fraction of 64%, the finite element mesh would consist approximately of 70 million elements each consisting of four RVEs and as many nodes. This would result in approximately 500 million Gauss points and 200 million of degrees of freedom. A pure FE 2 calculation would involve 500 million x 2 = 600 million hours of calculation time. The use of parallel computing can be seen to be indispensible.

So as to reduce the calculation time at the macroscopic scale certain assumptions have been made so as to further simplify the problem when considering the long-term macroscopic behaviour of a unidirectional composite. These hypotheses are based on experimental results and do not affect the quality of the final result.

At the macroscopic scale, the composite can be seen as an equivalent homogeneous material, which can be described within orthorhombic coordinates, b loc = ( ( ( ( ) ) ) )

x x x 3 2 1 r r r , ,
. The x1 r vector is aligned parallel to the fibre direction. At this scale at the point M, the stress tensor is and Σ and the strain tensor is Ε.

The first hypothesis is that the fibre failures are the only mechanism determining the failure of the structure. That is to say that the principal mechanism controlling the failure of a unidirectional composite is fibre failures and that non linear behaviour, seen when such composites are loaded off axis, such as intralaminar cracking and viscoelastic behaviour, are not considered. Secondly it is assumed that, unlike the microscopic behaviour, the viscoelastic effects of the matrix, on the macroscopic scale, are negligible, as is seen experimentally. The The variation of the elastic modulus of the equivalent homogeneous material as a function of the density of fibre breaks could be solved as a problem of periodic homogenisation. However, an approximation, based on the law of mixtures, which gives an upper bound, gives an acceptable result. Finally at the macroscopic scale and in the b loc base the behaviour of the unidirectional composite the behaviour of the unidirectional composites is given by

( ( ( ( ) ) ) ) ( ( ( ( ) ) ) ) ( ( ( ( ) ) )
)

M E M a M ∑ ∑ ∑ ∑ = = = = :
with a(M) = (M) a 0 , where (M) a 0 describes the rigidity tensor of the undamaged material at the point M, except

( ( ( ( ) ) ) )                         - -- - = = = = T R N M N a M a 1 ) ( 0 1111 1111
in which N R , N T (32 in our case) are respectively the number of broken fibres and the total number of fibres in the RVE.

As a result of the second assumption it can be seen that it is only the increase in load in the intact fibres neighbouring a fibre break at the microscopic level which is important in determining behaviour.

Also the second hypothesis justifies only the increase in loading parallel to the fibre axis being considered and not in other directions.

Multiscale calculation

The input data for this calculation are:

at the macroscopic scale: the characteristics of the undamaged composite and the law describing the fall in longitudinal rigidity of the equivalent homogonous material; -at the microscopic scale, for each Gauss point: the random choice of 5 failure stresses according to the Weibull statistics, shown experimentally to describe the stochastic failure strengths of the fibres; the data base allowing the axial stresses supported by the fibres in the RVE to be known and the evolution of the damage from the undamaged state in CS32 to complete failure represented in C2.

The multiscale simplified FE 2 calculation is iterative and calculates the state of damage (number of fibres broken), at iterative step n, having determined the state of damage at step n-1 in the following manner:

a time increment is considered together with or without a change to the macroscopic stress depending on whether an increase in load is made or whether the test is at constant load; -the finite element calculation at the macroscopic scale determines the macroscopic stress on the composite in the direction of the fibres; -at each Gauss point for the microscopic scale:

• a localisation step -Knowing the macroscopic stress and the damage state of the material, the data base for load transfer is applied so as to calculate the axial stress in the fibres contained in the RVEs. This stress allows the number of broken fibres in the • a homogenisation step -so as to calculate the macroscopic behaviour of the homogenised equivalent material, taking into account the new state of damage; -the next incremental step in the iteration process as a function of time can be calculated.

Model applied to unidirectional plates

The approach described above to model damage accumulation has been used to predict the failure of unidirectional plate specimens loaded in the fibre direction. It has been shown that the model accurately predicts tensile failure to within ±2% and that by increasing the complexity of the model to include fibre-matrix debonding the scatter of composite properties can be simulated [START_REF] Blassiau | Micromechanisms of load transfer in a unidirectional carbon-epoxy composite due to fibre failures, Part 3[END_REF]. The same reference explores the effects of the viscoelastic properties of the matrix on the long term behaviour of the composite under a steady load. It is shown that the model is capable of predicting the scatter of rates of fibre failures when compared to experimental curves obtained by acoustic emission. The present study has taken into account the variation of Weibull modulus which occurs when selecting just thirty two fibres from a very large fibre population. It has been shown that, even though the large population will have a determinant value of Weibull modulus, the random choice of thirty two fibres will lead to a wide distribution of the value [START_REF] Berger | Statistical analysis of the failure stresses of ceramic fibres: Dependence of the Weibull parameters on the gauge length, diameter variation and fluctuation of defect density[END_REF]. The results of including this effect can be seen in Figure 4 (21). These results confirm that damage increases with time and that the spread of the increase of damage, when different specimens are subjected to the same applied load, is increased over and above that obtained due to the scatter in fibre strengths because of the effects of considering a limited number of fibres in the RVE. When these results are compared to experimentally obtained curves, as shown in Figure 5 there is very good agreement. Although, at the microscopic scale and near fibre breaks the compliance of the material is altered there is no discernable effect on macroscopic strain as the breaks are isolated in the body of the composite. This means that although breaks can be detected by acoustic emission there is no measureable increase in macroscopic strain.

Application to the behaviour of a composite pressure vessel

The above model can be applied to filament wound structures, such as tubes; however our immediate interest is in understanding the behaviour of composite pressure vessels. So as to extend the above approach to a pressure vessel, the structure has been modelled using finite elements, as shown in Figure 6. Calculating the stress fields in the pressure vessel and applying the above multi-scale model allows damage accumulation to be calculated as a function of time and pressure.

Figure 7 shows a simulation of a monotonic burst test on a pressure vessel with a fibre volume fraction of 64%.The dotted line shows the initial linear curve of pressure related to radial expansion of the vessel. The open circles represent points calculated by the model of pressure versus radial expansion showing the deviation from linearity corresponding to the beginning of mechanical instability. The filled circles represent the average number of fibre breaks in the RVEs, taken over the entire pressure vessel. Very few fibres are broken in plies which are not at 90° with respect to the long axis of the vessel and even in the fretting layer, the low number emphasises that failure is due to localised coalescence of breaks. The point at which the curve begins to deviate significantly from linearity is taken to show when fibre breaks begin to interact and coalesce, resulting in a physical instability of the pressure vessel, as has been shown experimentally by high resolution tomography [START_REF] Scott | Advanced micromechanical analysis of highly loaded hybrid composite structures[END_REF].

Figure 8 shows a simulation of the accumulation of average number of fibre breaks in each RVE taken over the whole pressure vessel when it is held at 50MPa, The pressure vessel is designed to have a burst pressure of approximately 80 MPa,. The rate of damage can be seen to fall with time. Again, the average numbers are low as ultimate failure is due to clustering of fibre breaks rather than the overall failure rate reaching a critical level. The number of fibres in the fretting plies at 90° to the major axis of the pressure vessel is approximately 100x10 6 . Fibres in off-axis plies do not significantly contribute to numbers of fibre failures or ultimate burst strength. The result reflects the experimental curves obtained experimentally on pressure vessels and published elsewhere [START_REF] Blassiau | Micromechanisms of load transfer in a unidirectional carbon-epoxy composite due to fibre failures, Part 2: Influence of viscoelastic and plastic matrices on the mechanism of load transfer[END_REF]. As with the plate specimens, isolated breaks do not produce a measureable overall change in dimensions of the pressure vessel. A physical instability arises when the local density of fibre breaks, at some point or points in the structure, reaches a critical threshold and the regions of damage begin to interfere. This is shown in Figure 9. The point of physical instability is shown by the point of inflexion in the curve of damage accumulation for the simulated pressure vessel subjected to 70 MPa (700 bars). Although this point does not represent the burst of the pressure vessel it does represent a limit for its reliable use.

The inflection point, representing the limit of reliability, occurs after increasingly long periods of loading as the applied pressure is reduced and so are not shown for lower pressures.

Proof testing of composite pressure vessels

The above analysis, together with the experimental results can be used to assess the state of damage of a pressure vessel even if it has not been previously monitored, providing that its loading history is known. Acoustic emission has been used by other authors to detect damage in carbon fibre composite pressure vessels [START_REF] Downs | Acoustic emission from depressurization to detect/evaluate significance of impact damage to graphite/epoxy pressure vessels[END_REF] however usually to detect accidental damage rather than to predict future behaviour. Figures 8 and9 reveal the state of damage as a function of time for a pressure vessel held at a steady load. If the loading pattern to which the pressure vessel has been subjected and over what period of time, is known, the results of this and previous studies have shown that the acoustic emission produced by the failure of the carbon fibres will follow a similar curve. The manner in which random fibre breaks have been induced in service is not of concern as it is the fibre break density which governs the acoustic activity when the pressure vessel is subjected to a constant pressure test. A master curve can therefore be drawn, using the above model, for an ideal pressure vessel which accumulates damage in this way and arrives at its inflexion point, or state of instability, exactly at the end of the period required; say twenty years. Any pressure vessel subjected to a similar load history will lie on such a curve and the question which needs to be answered is, will its curve reach the point of instability before or after that of the master curve. The proof test consists of subjecting the pressure vessel to a constant pressure equal to the maximum in-service pressure and monitoring the acoustic emission rate over a given period, of say twenty four hours. If the inflection point of the pressure vessel is greater than the required service life then the curve on which the pressure vessel lies will be under that of the master curve. A comparison of the numbers of emissions over the same period, of say twenty four hours will give a smaller number of emissions for the pressure vessel under test than that predicted by the master curve.

In this case the vessel can remain in service. If not it should be removed. This comparison of curves is shown schematically in Figure 10 (3). The proposed test does not require knowledge of the above model, which however will allow the master curve to be determined for a given series of pressure vessels. The model can also be extended to predict the critical damage level corresponding to the point of inflexion.

Conclusion

The accumulation of fibre breaks in unidirectional composites has been accurately modelled.

A simplified multiscale finite element approach has been employed in which the RVE takes into account the stochastic characteristics of fibre breaks, the variation of Weibull moduli due to finite numbers of fibres being considered, the load transfer through the matrix and its viscoelastic behaviour, as well as interfacial debonding. It has been shown that this model accurately describes the experimental results obtained in monotonic tensile tests and also in steady load tests. The model allows the scatter to be expected in experimental tests to be determined.

The use of parallel computation has allowed a multiscale model to be developed that allows the behaviour of filament wound structures and in the present study the whole of the pressure vessels to be described. The analogy between the behaviour of unidirectional composites and that of the reinforcements in the pressure vessels has been justified by a comparison between experimental and theoretical results. It has been shown that when the damage density is such that random fibre breaks are no longer isolated in the matrix but that an interaction begins to occur between damaged regions, the rate of damage begins to accelerate. This represents a point of physical instability which determines the limit of the safe use of a pressure vessel.

The proof testing technique which is proposed is based on an understanding of the real failure processes controlling the stability of the pressure vessel, rather than relying on intuition or models based on the behaviour of metals. The technique is simple to implement and therefore of low cost. It should allow this technology to progress at a faster rate than hitherto by optimising composite pressure design and allowing a reliable evaluation of residual lifetimes to be made.

Figure 1 :

 1 Figure 1: Natural gas cylinders made from carbon fibre reinforced epoxy resin on the roof of a bus.

3 Figure 2 :

 32 Figure 2: A hydraulic test would greatly increase the number of broken fibres in the pressure vessel.

5 Figure 3 :

 53 Figure 3: Model of the failure process on the scale of the carbon fibres.

Formatted

KK

  is the additional part of the load transfer exclusively due to the debonding, is exclusively due to the viscoelastic properties of the matrix,

  data base describing the microscopic scale which needs to be taken into account in Formatted: Lowered by 5 pt Formatted: Lowered by 5 pt Formatted: Lowered by 5 pt Formatted: Lowered by 5 pt Formatted: Lowered by 5 pt Formatted: Lowered by 5 pt Deleted: ¶ Deleted: ¶ Deleted: O Deleted: in which Deleted: ¶ ¶ Deleted: is the part of the load transfer exclusively due to the failure of the fibre ¶ Deleted: ¶ Deleted: is due to the debond with 0 is not considered ¶ ¶ Deleted: is due exclusively to the viscoelastic behaviour of the matrix with 0

  of the homogenised material is considered only to be due to fibre failure and reflects the density of broken fibres.

9 RVE

 9 to develop by comparison with the failure stresses selected from the Weibull analysis;

Figure 4 :

 4 Figure 4: Simulated fibre failures in unidirectional composites loaded to 80% of failure load showing scatter due to the stochastic nature of fibre strengths and also variations in Weibull moduli between the A and B curves.

Figure 5 :

 5 Figure 5: Comparison of simulated and experimental results of constant load tests, at 80% failure stress, with unidirectional specimens loaded in the fibre direction.

10 Figure 6 :

 106 Figure 6: A finite element representation of a composite pressure vessel

Figure 7 :

 7 Figure 7: Simulated burst tests a carbon fibre epoxy resin composite pressure vessel, with a fibre volume fraction of 64%. The open dotted line shows the initial linear curve of pressure related to radial expansion of the vessel. The open circles represent points calculated by the model of pressure versus radial expansion showing the deviation from linearity corresponding to the beginning of mechanical instability. The filled circles represent the average number of fibre breaks in the RVEs, taken over the entire pressure vessel. The low number emphasises that failure is due to localised coalescence of breaks.

Figure 8

 8 Figure 8 Simulated fibre failures in a pressure vessel, designed to fail at 80MPa, pressurised to 50MPa and then maintained at that constant pressure showing the average number of fibre breaks in the RVEs as a function of time. The rate of damage can be seen to reduce with time. The number of fibres in the 90° plies, which control the burst pressure, is approximately one hundred million. Off-axis plies do not contribute significantly to the number of fibre failures or burst pressure.

Figure 8

 8 Figure8shows that the damage accumulated in the composite pressure vessel increases with time, however, this damage consists of regions of fibre breaks isolated by the effect of load

Figure 9 :

 9 Figure 9: Fibre breaks occurring in a pressure vessel during pressurisation and then during the maintenance of the pressure constant. The curve for a pressure of 50 MPa (500 bars) shown in the previous figure is reduced because of the reduction in scale. The point of inflection seen on the curve for 70 MPa (700 bars) reveals the point at which regions of damage begin to interact and the pressure vessel is no longer reliable (21).

Figure 10 :

 10 Figure 10: The comparison of the acoustic activity, representing fibre failures, between a pressure vessel in service and the master curve of an ideal vessel which takes exactly the service life required to reach its point of instability.

Figure 1 :

 1 Figure 1: Natural gas cylinders made from carbon fibre reinforced epoxy resin on the roof of a bus.

Figure 2 :Figure 3 :

 23 Figure 2: A hydraulic test would greatly increase the number of broken fibres in the pressure vessel.

Figure 4 :Figure 5 :

 45 Figure 4: Simulated fibre failures in unidirectional composites loaded to 80% of failure load showing scatter due to the stochastic nature of fibre strengths and also variations in Weibull moduli between the A and B curves.

Figure 6 :Figure 7 :Figure 8 Figure 9 :Figure 10 :

 678910 Figure 6: A finite element representation of a composite pressure vessel

  The test would involve simply monitoring damage by technique for a predetermined time and a comparison with the results predicted for the master curve. This could be automated to make the test simpler.

	the acoustic emission
	F o r	Deleted: 20
	P
	e
	e r	Deleted: 21
	R	Deleted: ,
	e
	v i e
	w
	O n l
	y
	12

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

Acknowledgement

The authors wish to thank Dr. D.H. Allen, Dean of Engineering at the University of Nebraska, Lincoln, USA for his support and valuable suggestions as well as Mr.S.Camara, of the same university, for help with the study and for the supply of some of the figures.


Flyura Djurabekova

Lorenzo Malerba

Roberto C Pasianot

Pär Olsson

Kai Nordlund

Keywords: diffusion, molecular dynamic simulations, Monte-Carlo, vacancies Kinetic Monte Carlo, di-vacancy, vacancy cluster, diffusion, metals mc.manuscriptcentral

Introduction

The microstructure and microchemical evolution in materials subjected to irradiation depends largely on the kinetics of thermally activated processes occurring at the atomic scale. Kinetic Monte Carlo (KMC) numerical methods are especially suitable to model the evolution of complex systems, when a large number of such processes are possible [?, ]]Chatt07,Fich91. They are therefore nowadays widely used, especially to treat problems of phase separation or damage evolution in materials subjected to annealing or irradiation [?, ]]Bouar02,Fu05,Krasn07,Domain99,Schmau02,Vincent08,Sois08. However, they require the precise knowledge of all involved activation energies, for all possible processes. The correct, on-the-fly evaluation of these energies, in such a way that the computational model is accurate, while remaining applicable in practice, is far from being a trivial problem. Even in the long studied case of atomic species redistribution in alloys via the vacancy diffusion mechanism [?, ]]Bouar02, Domain99, Schmau02,Vincent08, Young66 the energy barrier not only depends on the type of atom, but is also a complicated and a priori unknown function of the local atomic configuration [?, ]]Djurabek07a. Many heuristic approaches based on total energy calculations have been traditionally used [?, ]]Young66, Kang89 to allow at least partly for this dependence. The total energy of the system is calculated using either an empirical interatomic potential (EIP) [?, ]]Domain99, or pair interaction energies fitted either to an EIP [?, ]]Schmau02, Bouar02 or to density functional theory (DFT) calculations [?, ]]Vincent08, Sois08, often using a broken bond scheme [?, ]]Schmau02,Bouar02,Sois08. The main concern of these approaches is to ensure that the transition rates respect the detailed balance [?, ]]Fich91, so as to be thermodynamically consistent. Often, for these conditions to be fulfilled, the characteristic energies of the relevant thermally-activated processes are explicitly correlated to energy differences between initial and final local equilibrium states [?, ]]Young66, Domain99,Vincent08.

However, the fact that a certain state is thermodynamically favoured does not necessarily mean that it is going to appear faster or more often than less favoured ones. The broken bond scheme, in which the activation energy is obtained as the difference between saddle point and initial energies, both expressed in terms of pair interactions [?, ]]Schmau02, Bouar02, Vincent08, Sois08, does attempt to be more realistic. Nonetheless, being based on a rigid lattice pair energy scheme, it remains a strongly simplified approach [?, ]]Djurabek07a.

One way of improving the accuracy of an atomistic KMC (AKMC) model in the sense discussed here is to feed it with activation energies rigorously calculated as functions of the local atomic configuration using, for example, an EIP, or DFT methods. If the number of involved energies is limited, this can be done by pre-tabulating all possibilities [?, ]]Djurabek07b, Barash05; if it is combinatorially large, advanced regression methods with good predictive capability (e.g. artificial intelligence) may be used [?, ]]Djurabek07a,Shastry05,Castin08c. In models of this type, thermodynamics does not appear explicitly in terms of binding energies, energy differences between states or favoured kinetic paths towards lower energy states. It is, instead, "hidden" in the difference between activation energies of direct and inverse processes, i.e., the energy difference between an initial state A and a final state B can be obtained as difference between the barrier from A to B and the barrier from B to A.

In this Letter we present an application of such an AKMC model for the apparently simple case of the migration of a di-vacancy in iron to illustrate the fact that energy barriers dont always correlate with energy differences between initial and final state, as previously pointed out for the simpler case of a monovacancy [?, ]]Djurabek07a,Bocquet02. We choose to study this particular case since di-vacancies in iron (and probably in all or most bcc metals) are most stable in 2nn configuration, as consequence positron annihilation does not allow di-vacancies to be distinguished from single vacancies [?, ]]Kuri06. Thus, there is a lack of experimental data on di-vacancies in iron and this is in itself a motivation for our work. At the same time, iron is a material of primary importance for applications, within and without the nuclear industry, which is why a number a of EIPs to be compared exist. We show that the activation energies do not necessarily correlate with the energy difference between the involved local equilibrium states and may increase the frequency of appearance of states that are not necessarily the thermodynamically most favoured ones. Nonetheless, the model does respect thermodynamics, as is demonstrated by solving the corresponding master equations, which can be explicitly shown to embody Boltzmann's statistics (its stationary solution is thermodynamic equilibrium) [?, ]]VanKampen03, and comparing it to the AKMC simulation results. We also show, however, that even in a DFT framework relatively small differences between activation energies may appear, depending on the used approximation, with consequences that may be non-negligible a priori. In particular, our study reveals that the migration mechanism for the di-vacancy in iron is more complex than one would intuitively deduce from simple di-vacancy stability considerations. Such a mechanism cannot be predicted by a model that explicitly couples activation energies and energy differences between local equilibrium states and requires an atomic-level study to be clearly identified.

F

Methodology

KMC methods determine the evolution of a system by stochastically choosing events whose probabilities are given in terms of frequencies. The latter are written using the classical transition-state-theory expression for the frequency of thermally activated phenomena: Γ = νe -Ea/kB T , where k B is Boltzmann's constant, T the absolute temperature, ν the attempt frequency (here considered constant and equal to 6 ×10 12 s -1 ), and E a the activation energy. The associated time is estimated using the residence time algorithm [?, ]]Bortz75. In our case, E a corresponds to the energy barrier for the exchange of a vacancy with a neighbouring atom. The possible states of the cluster were defined by the mutual distance between the two vacancies, in terms of nearest neighbour (nn) shells. There are five possible states of the bound configurations of a di-vacancy: 1nn, 2nn, 3nn, 4nn, and 5nn. The possible states and transitions are pictorially represented in Fig. 1. Three transitions, namely 3nn → ∞, 4nn → ∞, and 5nn → ∞, lead to di-vacancy splitting. The energy barriers for all possible transitions were calculated by the drag method [?, ]]Henkel00 using (for a check of trends) four widely used EIPs, denoted here as: AB [?, ]]Ackland97, MH [?, ]]Mend03, AM [?, ]]Ackland04, and DD [?, ]]Dudarev05. For these simple transitions, the drag method has been checked to provide the same results as the nudged-elastic band method (NEB) [?, ]] NEBb. All the classical potential calculations were carried out with the DYMOKA code [?, ]]DYMOKA. The corresponding DFT calculations were performed in two different ways, namely using pseudopotentials generated both within the projector augmented wave (PAW) approach [?, ]]Kresse99 and in the form of ultrasoft pseudopotentials (USPP) [?, ]]Kresse99,Vander90,Kresse96a, as implemented in the Vienna ab initio simulation package (VASP) [?, ]]Kresse96b. The USPP formalism can be derived from the PAW formalism by linearization. In general, the USPP method does not describe transition metals and magnetic systems as well as the PAW method does. However, a large amount of DFT data found in the literature concerning Fe, and impurities in it, have been in fact obtained in the USPP approximation[?, ]]Do-main01,Becquart03, Vincent05,Vincent06b and used to parametrise AKMC models [?, ]]Vincent08. All calculations were spin polarized and the exchange correlation functional used was the generalised gradient approximation [?, ]]PW91 with the Vosko-Wilk-Nusair correlation interpolation [?, ]]VWN for the PAW approach. The barriers were calculated using the nudged elastic band method [?, ]]NEBb implemented in VASP on supercells of 128 bcc sites at constant volume conditions. The reciprocal space was sampled with 27 k-points according to the method by [?, ]]MonkPack. The plane wave energy cutoff was 300 eV for PAW and 240 eV for USPP. The convergence using these parameters has been reported earlier [?, ]]Domain01,Olsson.

Results and Discussion

Fig. 2 shows the schematic energy landscape, according to the AKMC model, of di-vacancy states and transitions along the path connecting two dissociated states through, in order, 5nn, 1nn, 2nn, and 4nn states. The dissociated state is the reference, i.e. has zero energy. According to all used model Hamiltonians, there are three bound states, namely (in order of increasing binding) 4nn, 1nn, and 2nn. DFT predicts also the 5nn state to be bound (no EIP grasps this feature). According to USPP calculations 1nn and 4nn states are almost equal in energy, while all EIPs and PAW calculations predict 1nn to be more strongly bound. Yet, all EIPs favour the transition from the 2nn ground-state to 4nn over the transition to 1nn. Interestingly, the same result, though unnoticed, was also obtained with a pair potential for iron in a work published almost 40 years ago [?, ]]Tsai70. We also obtained the same trend using two other EIPs for Fe, from Ref. with [?, ]]Djurabek07a,Bocquet02. Of all EIPs, the one globally providing the closest results to DFT is AM (or, almost indifferently, MH) and was therefore used for the comparative studies.

The AKMC simulations were performed with pre-tabulated energy barrier values [?, ]]Djurabek07b, calculated as described above. The di-vacancy was always introduced in the simulation box in 2nn state, which is the most stable according to all used model Hamiltonians. The cluster was considered to exist so long as the two vacancies did not part beyond 5nn distance. The simulation was repeated 50 times at different temperatures between 250 K and 500 K. The migration properties of the di-vacancy were studied in detail, by monitoring the visited configurations, their frequency, the time spent in each of them and the total lifetime of the complex, as well as its mean free path before dissociation. The effective migration and dissociation energies of the di-vacancy were obtained as slopes of the corresponding Arrhenius plots for the diffusion coefficients and lifetimes calculated with barriers from different EIPs, as well as from USPP and PAW DFT calculations (Table 1). In addition, for comparison the same type of simulation was also performed by estimating the energy barriers with the heuristic equation proposed in [?, ]]Domain99 and used, e.g., in [?, ]]Kang89:

E b = E 0 + Ef -Ei 2 (1)
where E i and E f are, respectively, the total energy before and after the vacancy jump, calculated on rigid lattice at each AKMC step using the AM potential. E 0 is a constant barrier, customarily assumed to equal the single vacancy migration energy in bcc-iron. The latter was calculated to be 0.63 eV with the drag method using the AM and MH potentials. The barriers calculated from the different EIPs and DFT varieties provide different mean free paths of the di-vacancy, consistently with the different dissociation and migration energies predicted. Despite the fact that the DD EIP gives the highest effective migration energy for the di-vacancy, the even higher barrier for dissociation from the 4nn state (Fig. 2) allows the di-vacancy to ramble the longest distance during its lifetime (Fig. 3).

Fig. 4 shows the occupation statistics of di-vacancy states N i /(N 2 (n AKM C = 0)) versus AKMC steps according to AM-calculated barriers (Fig. 4 from eq. (1) (Fig. 4(b)), obtained as average over 50 different simulations at 300 K.

Here N i is the number of simulations in which the di-vacancy was found in the inn state at each AKMC step, n AKM C , and N 2 (n AKM C = 0) = 50 the total number of simulations (every di-vacancy starts from the 2nn state). The curves tend to zero due to the eventual dissociation of all di-vacancies. Consistently with the energy landscape of Fig. This means that the di-vacancy migrates by oscillating between 2nn and 4nn states rather than 2nn and 1nn (Fig. 4(a)). What would be intuitively expected, knowing that the 1nn state is more strongly bound than 4nn is that, on the contrary, while migrating the di-vacancy should oscillate between 2nn and 1nn states. The latter mechanism, which is the only one that has been considered even in recent DFT investigations [?, ]]Fu05, is, not surprisingly, predicted by an AKMC model in which eq. ( 1) is used (Fig. 4(b)). Here the presence of 2nn and 1nn states is very clear while the 4nn state hardly ever appears.

A different migration mechanism should imply, in principle, also different effective migration and dissociation energies and, therefore, a different effect on mass transportation in metal alloys of interest due to the thermal diffusion. Although we do observe differences between the diffusion parameters assessed from AKMC simulations with tabulated barriers and with barriers obtained from eq. ( 1), this difference is surprisingly small. The involved energy values are rather similar as there is an accidental compensation between predicted migration energies and binding energies in rigid lattice (eq. ( 1)) with the AM potential. However, we emphasise that in different circumstances (no rigid lattice, different EIP) this similarity will most likely vanish, and hence it is important to know the true mechanism of cluster migration.

The observed counterintuitive divacancy migration mechanism is a direct consequence of the values of the energy barriers. The question to be answered is whether the 2nn-4nn mechanism is consistent with thermodynamics, because Fig. 4(a) may wrongly lead to the conclusion that Boltzmann's statistics is violated, if a higher energy state is occupied much more often than a lower energy state. To answer this question it is necessary to trace the probability density of states as a function of time, instead of AKMC steps, since the residence time in the different states is different, regardless of how often they appear. This is obtained in a natural way by solving the master equation system that governs the process [?, ]]Fich91, Chatt07. Note that we are solving here the master equation for an open system, consistent with the AKMC simulations (where we start the simulations from a di-vacancy and we stop collecting data when the di-vacancy dissociates):

                               dC1 dt = K 21 • C 2 + K 31 • C 3 + K 51 • C 5 -(K 12 + K 13 + K 15 ) • C 1 dC2 dt = K 12 • C 1 + K 42 • C 4 -(K 21 + K 24 ) • C 2 dC3 dt = K 13 • C 1 + K 43 • C 4 -(K 31 + K 34 + K 3∞ ) • C 3 dC4 dt = K 24 • C 2 + K 34 • C 3 + K 54 • C 5 -(K 42 + K 43 + K 45 + K 4∞ ) • C 4 dC5 dt = K 15 • C 1 + K 45 • C 4 -(K 51 + K 54 + K 5∞ ) • C 5 (2)
Here C i (i=1,2,. . .,5) denote the probability density of di-vacancies in the inn state at a given time. The solution of the equations was obtained with the initial condition:

C 2 =1 and C 1 = C 3 = C 4 = C 5 = 0 at 300 K. The K ij coefficients
K ij = k ij • Γ ij = k ij • ν 0 e -E ij m /kBT (3)
i.e. they are obtained by multiplying the rate Γ ij of the i → j transition, times the number of possible ways in which this transition may occur, k ij (multiplicity of equivalent transitions). These multiplicities solely depend on crystallographic structure and geometry of the cluster:

k 12 = k 13 = k 54 = 6; k 15 = k 42 = k 45 = k 51 = 2; k 21 = k 24 = k 34 = k 4∞ = k 5∞ = 8; k 31 = k 3∞ = k 43 = 4.
All the multiplicities are doubled taking into account the equal probabilities of both vacancies to jump for each transition. Multiplicities of not listed transitions k ij are zero as they are impossible in the bcc-structure. . Probability density of di-vacancy states C i = N i /(N 2 (t = 0)) at 300 K as a function of time for the same bound states as in Fig. 4, using barriers from the AM potential. N i is the number of simulations in which the di-vacancy is found in the inn state within the considered ∆t < 10 -5 ; N 2 (t = 0) = 500 is the total amount of the simulations, each starting from a di-vacancy in 2nn state. The dark grey lines are "horizontal" averages of the light grey ones over larger time intervals (∆t ∼ 3 × 10 -3 s). The solid black lines are the solutions of eqs. (2). The inset (left bottom corner) zooms in the beginning of the 1nn and 4nn distributions.

In Fig. 5 we compare the results for the probability density of states corresponding to the AKMC simulations of Fig. 4(a) , plotted versus time, with the numerical solution of eqs. (2) at 300 K (solid black lines). We increased the statistics to 500 AKMC simulations to reduce the statistical fluctuations. The comparison requires a proper treatment of the AKMC simulation data, because time advances at a different pace in each simulation according to residence time algorithm. This happens because the ∆t associated with each AKMC step is not constant, thus the same number of AKMC steps corresponds to a different time. When the statistical averaging is applied (as described in the figure caption), the AKMC results are directly comparable, and in fact in excellent agreement, with the solution of the master equations (2).

The probability that, in an infinitesimal time interval, the di-vacancy is found in a 2nn state is always the highest, followed by 1nn (about one order of magnitude lower) and, finally, 4nn (about two orders of magnitude lower). This result is fully consistent with the binding energies of Fig. 2,i embodied by the EIP, even though the migration mechanism favours the appearance of the 4nn state, as shown in Fig. 4(a). Analytically, under a few simplifying assumptions, it is possible to show that the solution of the master equations corresponds indeed to classical Boltzmann's statistics. For example, if we do not allow dissociation in eqs. (2), the ratio between the concentrations of 1nn and 4nn divacancies after infinite time (obtained imposing that the derivatives on the left hand are zero), C ∞ 4 /C ∞ 1 , after neglecting C ∞ 3 and C ∞ 5 (which have been found to be very small) is:

C ∞ 1 C ∞ 4 ≈ K 42 K 21 K 12 K 24 ∝ exp[ (E 1 b -E 4 b ) kBT ] (4)
where we have equated

E 24 m -E 42 m = E 2 b -E 4 b and E 21 m -E 12 m = E 2 b -E 1 b , the binding energies E i
b being positive values (depth of the wells in Fig. 2). This ratio, using AM binding energies, will be significantly larger than 1 for all temperatures of interest. On the other hand, it is clear also from Fig. 5 that the C 1 /C 4 ratio remains essentially constant, on average, at a value much larger than 1 (about 28, as has been verified), independently of the fact that the total di-vacancy population is dwindling. This means that the steady-state in terms of visited configurations and also time spent in them is reached almost instantaneously and that, correspondingly, the time spent in the 1nn state is always larger than in 4nn. Nonetheless, so long as the barriers are those predicted by the model Hamiltonian we have used, the mechanism of migration of the di-vacancy remains 2nn-4nn, i.e. the number of times the 4nn configuration is visited is significantly larger than for 1nn, even though the time spent in the former is shorter. This is implicitly demonstrated by the inset in the lower left corner of Fig. 5, which blows up the beginning of the curves. Since at the outset all di-vacancies are in the 2nn state, both 1nn and 4nn curves start from zero but, of the two, the 4nn one grows faster and reaches earlier the steady-state (before di-vacancies start to dissociate), in agreement with the faster kinetics dictated by the AM energy barriers of Fig. 2.

We also performed AKMC simulations with the barriers that should be considered the most reliable ones, i.e. those obtained from DFT(PAW) calculations. In this case, the migration mechanism was somewhat different, namely a mixture of 2nn-1nn and 2nn-4nn exchanges.

Conclusions

In summary, we have shown how a kinetic path can be "thermodynamically counterintuitive", as a consequence of the relative value of the involved activation energies, while in fact fully respecting the laws of thermodynamics. If a proper calculation of all relevant activation energies is performed, the mechanism of migration of the di-vacancy in iron consists in an oscillation between not only 2nn and 1nn states, but also, largely, between 2nn and 4nn states, despite the fact that the 1nn state is the closest in energy to the fundamental 2nn state. Such a "thermodynamically counterintuitive" mechanism could not be predicted by a model linking explicitly the activation energies to the energy difference between initial and final state. The use of models in which the kinetics is explicitly linked to these energy differences [?, ]]Young66, Domain99,Vincent08 may therefore suggest mechanisms that, though thermodynamically acceptable, do not in fact occur. This emphasises the importance of devoting great care and effort to a proper evaluation of activation energies, to ensure the accuracy of models describing the evolution of physical systems driven by thermally activated processes, possibly avoiding to resort to traditionally used heuristic equations.

This work was performed in the framework of the bilateral collaboration agreement between SCK•CEN and CNEA, sponsored by the Belgian Scientific Policy Office, under contract BL/52/A01. The DFT calculations were performed on the supercomputers at Centre de Calcul Recherche et Technologie (CCRT) in the framework of an EDF-CEA contract.

Figure 1 .Figure 2 .

 12 Figure 1. Different possible di-vacancy states and transitions. 0 is the position of the first vacancy; different states are obtained depending on the location of the second one. ∞ positions correspond to a dissociated di-vacancy. Transitions are possible only via the exchange of the given vacancy with a nearest neighbour atom (e.g. 1→2).

Figure 3 .

 3 Figure 3. Mean free paths of di-vacancy against inverse temperature. The lines are guides to the eye.

Figure 4 .

 4 Figure 4. Occupation statistics of three bound di-vacancy states (1nn, 2nn, and 4nn) visited by one di-vacancy at each AKMC step during its diffusion migration. The di-vacancy always starts from 2nn state. The statistics was defined by 50 simulations at 300 K using (a) pre-tabulated barriers and (b) barriers from eq.(1), both with the same potential.

Figure 5

 5 Figure 5. Probability density of di-vacancy states C i = N i /(N 2 (t = 0)) at 300 K as a function of time for the same bound states as in Fig.4, using barriers from the AM potential. N i is the number of simulations in which the di-vacancy is found in the inn state within the considered ∆t < 10 -5 ; N 2 (t = 0) = 500 is the total amount of the simulations, each starting from a di-vacancy in 2nn state. The dark grey lines are "horizontal" averages of the light grey ones over larger time intervals (∆t ∼ 3 × 10 -3 s). The solid black lines are the solutions of eqs.(2). The inset (left bottom corner) zooms in the beginning of the 1nn and 4nn distributions.

Table 1 .

 1 [?, ]]Simon93 and from Ref.[?, ]]Chak02. Finally, this result is also provided by DFT(USPP) calculations. Conversely, according to DFT(PAW) both transitions are almost equally probable. Thus, all tested EIPs favour a transition to a less strongly bound, and therefore thermodynamically less favoured, state. In DFT the situation is smoothed: with USPP two energetically equivalent states are reached with different frequency, while with PAW two states of different stability are reached with essentially the same frequency. Regardless of the method differences, most important for the conclusions of the current work is that in no case is the energy barrier correlated with the energy difference between states, in agreement Migration Em and dissociation E d energies as well as their intercepts with the y-axis D0 and τ0 of di-vacancy in α-Fe obtained by AKMC method using pre-tabulated barriers from the different EIP's and DFT varieties (see text). AM,eq.(1) denotes the results, where the barriers were assessed at each AKMC step (rigid lattice) from eq. (1).

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

  2 , 2nn and 4nn states are by far the most often observed ones.

	January 22, 2010 11:14 Philosophical Magazine Letters orig-Djurabekova-Malerba-resubmitted
	8	Submitted for publication
	F o r
	P
	e
	e r
		R
		e
		v i e
		w
		O n l
		y

  .e. with the thermodynamics

	January 22, 2010 11:14 Philosophical Magazine Letters orig-Djurabekova-Malerba-resubmitted
	10	Submitted for publication
	F o r
	P
	e
	e r
		R
		e
		v i e
		w
		O n l
		y

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