Modelling the coefficient of thermal expansion in graphite crystals: implications of lattice strain due to irradiation and pressure

Theoretical models for the coefficient of thermal expansion (CTE) first proposed in the 1970s are expanded upon, allowing them, for the first time, to be implemented over a wide temperature range. The models are of interest because they predict the effects of the changes in the crystal lattice spacing and crystallite modulus on the CTE. Hence, they can in turn be used to investigate the influence of pressure and irradiation on the CTE. To date, typographical and mathematical errors and incomplete or conflicting assumptions between the various papers had made the complex mathematical formulations difficult, if not impossible, to follow and apply. This paper has two main aims: firstly to revisit and review the CTE models, correcting the errors and compiling and updating various input data, secondly to use the revised models to investigate the effect of loading and irradiation on the CTE. In particular, the models have been applied to data for natural and highly orientated pyrolytic graphite and compared with experimental data, giving an insight into the influence of temperature, loading and irradiation on both single crystal and polycrystalline graphite. The findings lend credence to postulated microstructural mechanisms attributed to the in-reactor behaviour of nuclear graphite, which finds a wide use in predictive multiscale modelling.

Modelling the coefficient of thermal expansion in graphite crystals: implications of lattice strain due to irradiation and pressure Barry  Theoretical models for the coefficient of thermal expansion (CTE) first proposed in the 1970s are expanded upon, allowing them, for the first time, to be implemented over a wide temperature range. The models are of interest because they predict the effects of the changes in the crystal lattice spacing and crystallite modulus on the CTE. Hence, they can in turn be used to investigate the influence of pressure and irradiation on the CTE. To date, typographical and mathematical errors and incomplete or conflicting assumptions between the various papers had made the complex mathematical formulations difficult, if not impossible, to follow and apply. This paper has two main aims: firstly to revisit and review the CTE models, correcting the errors and compiling and updating various input data, secondly to use the revised models to investigate the effect of loading and irradiation on the CTE. In particular, the models have been applied to data for natural and highly orientated pyrolytic graphite and compared with experimental data, giving an insight into the influence of temperature, loading and irradiation on both single crystal and polycrystalline graphite. The findings lend credence to postulated microstructural 2018 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/ by/4.0/, which permits unrestricted use, provided the original author and source are credited.

Introduction
Significant property and dimensional changes are exhibited by polycrystalline graphite subjected to fast neutron irradiation [1] and high pressure [2]. Of particular interest are the changes to the coefficient of thermal expansion (CTE). It is postulated that changes to the CTE in polycrystalline graphite under loading and irradiation are related to a combination of changes to: the crystallites, the porosity, the crystal orientation, the microstructural stiffness or more likely a combination of these factors [3,4]. To separate out these effects, it is therefore important to understand and define the possible contribution that crystal CTE changes have on the bulk CTE. To this end, Kelly [5] derived several theoretical models capable of predicting the CTE of a graphite crystal as a function of temperature and crystal lattice spacing [6][7][8][9][10][11][12][13]. Unfortunately, to date, these theoretical models have not been widely taken up as typographical errors and incomplete or conflicting assumptions between the various papers make the complex mathematical formulations difficult, if not impossible, to follow and apply.
This paper revisits and reviews these CTE models, correcting the errors as well as compiling and updating various model input data. The models are used to investigate the effect of loading and irradiation on the CTE, and are compared with available experimental data. In particular, the models have been applied to data for natural and highly orientated pyrolytic graphite and compared with experimental data, giving an insight into the influence of temperature, loading and irradiation on both single crystal and polycrystalline graphite.

Development of a graphite semi-continuum coefficient of thermal expansion model
The relationship between crystal strains e ij and crystal stresses τ ij are defined by the Cartesian co-ordinate system, given in figure 1a, with the z-axis parallel to the hexagonal axis of the crystal c-axis as follows: The stresses τ ij are defined as the force acting on unit area parallel to the ith direction, the normal to the unit area being in the jth direction. The matrix elements, S ij , are the elastic compliance constants. This matrix can be inverted to obtain the elastic stiffness constants, C ij , as given below: 12 [14,15]; the 'equivalent cylinder' is shown dashed. (Online version in colour.) The relationships between the compliance and stiffness constants are: The stiffness constants for a graphite crystal are given in table 1 of appendix A; these were obtained using Born's long-wave method [14] and local density approximation (LDA) [ to the crystal structure in which atoms within layers are held together by strong covalent bonding, whereas bonding between layers is via weak van der Waals bonds, the values of the stiffness constants perpendicular, C 33 , and parallel, C 11 , to the basal planes are different by orders of magnitude. The shear components C 44 and C 66 = (C 11 − C 12 )/2 are also small.
The most thermodynamically stable form of graphite crystal is the ABAB hexagonal stacking, giving a theoretical density of 2.266 × 10 3 kg m −3 . The lattice a-spacing is 0.1415 nm and the c-spacing is 0.335 nm-calculated using LDA [15]; the lattice d-spacing is half that of the c-spacing.
The development of a semi-continuum model is based on the consideration of the free energy, F, in a unit cube of graphite crystallites in a thermally dilated state, assuming that there are no shear strains as follows: where the first term, U 0 , is the energy/unit volume at absolute zero temperature in the unrestrained state; the next four terms are the elastic strain energy density; and the last term is the thermal energy density of the first Brillouin zone, where v p is the frequency of the pth vibrational mode (Hz), k is Boltzmann's constant (1.38064852 × 10 −23 J K −1 ), h is Planck's constant (6.626070040 × 10 −34 J s −1 ), and T is the absolute temperature (K).
The integral is over the first Brillouin zone of the graphite crystal structure [16,17], i.e. the reciprocal lattice space between two adjacent basal planes; see figure 1b. Note that there must be a whole number of waves, denoted by the wavenumber, in the x and y planes. Therefore, the wavenumber z varies from ±1/2d = ±1.49 × 10 9 m −1 . To simplify the integration, the hexagonal prismatic Brillouin zone has been approximated as a cylinder with an equivalent radius, defined as: We use the short-hand notation d 3 σ ≡ dσ x dσ y dσ z . The integral will be performed in cylindrical co-ordinates (σ a , φ, σ z ) with σ a defined as σ a = σ 2 x + σ 2 y , tranforming to cylindrical co-ordinates and assuming cylindrical symmetry leads to It should be noted that previous authors [9] did not include the integral over the Brillouin zone at this stage of the derivation. In this paper, we explicitly separate the summation over frequencies from the summation over modes (integration over the first Brillouin zone) for clarity. This also removes confusion about the dimensionality of the different terms in equation (2.4).
The principal graphite crystallite coefficients of thermal expansion required are given by the differentials , which are the mode Grüneisen [18] parameters describing the frequency shifts due to strain which have previously been referred to as the anharmonic terms. In addition, the equations have been simplified by using the fact that, due to symmetry, the strains in the crystal basal planes are equal, i.e. e xx = e yy . It will be shown analytically later that the derivatives of the elastic constants are of greater magnitude than the elastic constants themselves; however, the elastic strains are very small, and so terms that are higher order in the elastic strains can be neglected. Thus, equation (2.7) can be simplified to Previously, the assumption was made that the higher order strain terms given above are negligible, although this is not discussed in [9]. Equations (2.8) can be solved for e xx and e yy and then differentiated with respect to temperature to give the principal crystal CTE as  [19,20], as described in the next section.

Graphite vibration model
To solve for α c and α a , we require the frequency of the vibrational modes in the longitudinal, transverse and out-of-plane principal directions and also the anharmonic modes -themselves functions of the principal vibrational frequencies. Below, we summarize the Komatsu [20] model for calculating the principal frequencies; this assumes that the graphite crystallites can be modelled as elastic extensional vibrations in a series of n thin plates [21,22]. The present authors modify Komatsu's notation [20] to match the rest of this paper. The equations for the elastic extensional vibrations in a series of n thin plates are where u n , v n , w n are the displacements at a lattice layer point (x, y) in the x-, y-and z-axis directions, differentiated with respect to time, t; ν is Poisson's ratio, ρ is the density and d is the graphite interlayer spacing; δ is related to the bending modulus of a plate. Kelly & Walker [9] give a value of δ = 6.11 × 10 −7 m 2 s −1 . In the third equation of (3.1), the second term on the r.h.s. accounts for the interaction between neighbouring plates. Komatsu [20] modified equations (3.1) to account for shear stress, obtaining the following terms: ∂w n+1 ∂y − ∂w n−1 ∂y and ∂ 2 w n ∂t 2 = −δ 2 ∂ 4 w n ∂x 4 + 2 ∂ 4 w n ∂x 2 ∂y 2 + ∂ 4 w n ∂y 4 + C 33 ρd 2 (w n+1 + w n−1 − 2w n ) The longitudinal (transverse) velocity is defined as the square root of the pre-factor in front of the x (y) second derivative term in the u n displacement equation. Assuming that Poisson's ratio  is zero, the velocity of the longitudinal wave, V L , can be simplified to Similarly, the velocity of the transverse wave V T is given by Equations (3.2) are solved using the conventional substitution in reciprocal space, i.e. in the first Brilliouin zone where σ x , σ y and σ z are the wavenumbers as discussed previously and v U , v V and v W are frequencies of the vibrational mode. Komatsu [20] solved for the principal vibrational modes, v p , which are linear and where again σ 2 a = σ 2 x + σ 2 y . It should be noted that, in the derivation, Komatsu [20] assumes that τ = C 44 /ρ is small and therefore cross terms with τ as a pre-factor are ignored.
The three vibration modes given in equations (3.6) are insensitive to σ z . For illustration, three of the vibration modes are plotted in figure    Equations (3.6) are the equations for lattice vibration used previously [9]. The anharmonic functions for the longitudinal, transverse and out-of-plane principal modes are as follows: To evaluate these, the differentials of the elastic constants ∂C ij /∂e kk and ∂δ/∂e ii are required with respect to the strain in the longitudinal and transverse in-plane x-axis and the out-of-plane z-axis. The values used in our simulations are given in table 2 (appendix A). In this paper, only the anharmonic functions [γ 3 ] zz , in their full form, are used in the derived expressions for α c . The inplane [γ 1 ] xx and [γ 2 ] xx anharmonics are simplified in the case of α a , and in effect can be regarded as constants.
The anharmonic functions [γ p ] zz are plotted as a function of σ a for different values of σ z in figure 3.
The differentials of the anharmonic functions with respect to e xx are discussed later with reference to the derivation of α a .

Coefficient of thermal expansion perpendicular to the basal plane-α c
Having obtained expressions for the orthogonal vibrations and anharmonic modes, expressions for the CTE can be derived using the semi-continuum model over a range of temperatures to compare the model with experimental data. Firstly, we consider the thermal expansion coefficient perpendicular to the basal plane.
It can be seen from table 1 (appendix A) that S 13 is two orders of magnitude smaller than S 33 , thus equation (2.11) simplifies to There are three possible vibration modes to consider when calculating α c : in-plane longitudinal mode, in-plane transverse mode and out-of-plane mode. Previously it has been argued that only the out-of-plane mode is significant. Unfortunately, the reasoning behind this assumption in [9] is unclear and there appear to be mistakes in the references and supporting data. However, due to the large difference in stiffness between the in-plane C 11 (1060 GPa) and the out-of-plane C 33 (36.5 GPa) it is not unreasonable to assume that the out-of-plane strain e zz would have less influence on the values of C 11 and C 12 than it would have on C 33 . This assumption appears to be justified by the implementation of the theory as demonstrated below. The final step is to perform the integral over a cylinder of equivalent volume to the true Brillouin zone, Data giving the changes in graphite crystal CTE as a function of temperature can be obtained by direct thermal expansion measurements or via changes to lattice spacing using X-ray diffraction (XRD). Measurements have been obtained using either naturally occurring graphite crystals or highly oriented pyrolytic graphite (HOPG). The quality of these measurements depends on the graphitic perfection of the sample used, which is usually defined using a p-factor [23,24]; p-factors range from zero to unity, representing the degree of cryptographic order. Measurements made on various graphitic structures with p-factors between 0 and 0.7, including irradiated graphite by Steward et al. [25,26], show that curves of changed interlayer d-spacing versus temperature are parallel to one another between ∼0 and ∼2600 K, indicating that the thermal expansion coefficients, which are related to the slopes of these curves, are unaffected by significant disorder of the graphite lattice. This observation agrees with the estimates of change in CTE due to strain and irradiation discussed later in this paper. Morgan [27] presents a comprehensive set of XRD measurements giving changes to α a and α c taken from Kellett & Richards [28], Nelson & Riley [29], Yang [30], Kellett et al. [31] and Baskin & Meyer [32].
Kelly [33] also presents a dataset of crystal CTE measurements, but there is confusion with regard to the source of some of the data. Figure 4 brings together various datasets for α c from Entwisle [34], Bailey & Yates [35], Yates et al. [36], Harrison [37], Nelson & Riley [29] and Morgan [27]. Where the measurement error on data are available error bars are given; most of the data originates from XRD data on graphitic materials with a p-factor of 0.2. The prediction of α c using equation (4.2) is also included. While equation (4.2) fits well at low values to medium values of temperature, it does not appear to capture the measured CTE increase at higher temperatures. Also plotted is an attempt to include the contribution to the CTE of the optical modes of the   system. It is possible to make the simplifying assumption [6] that extending the upper integration limit to √ 2σ m has the effect of including contributions from both the in-plane acoustic and the higher frequency optical mode. It is difficult to say with current experimental data whether the inclusion of optical modes produces a better fit.
The behaviour of the thermal expansion coefficient can be understood in the high and low temperature limits by examining the leading order behaviour for hv 3 kT and hv 3 kT.
For low temperatures, hv 3 kT, the leading order behaviour is For high temperatures, hv 3 kT, the leading order behaviour is We see that as the temperature tends to zero the thermal expansion coefficient exponentially approaches zero, whereas at high temperatures the thermal expansion coefficient loses all temperature dependence and becomes simply a constant. The sensitivities of the predictions of equation (4.2) to the differentials with respect to strain of δ, C 33 and C 44 are given in appendix B. Doubling ∂δ/∂e zz has a significant effect on the prediction but the effect of halving the value is much less, i.e. the influence is nonlinear. Changing ∂C 33 /∂e zz by ±1 × 10 11 Pa (a few per cent) increases/decreases the prediction by the same amount, whereas increasing the value of ∂C 44 /∂e zz by a factor of 4 increases the prediction significantly; however, reducing it by a factor 10 makes only a small difference, i.e. again the influence is nonlinear.

Coefficient of thermal expansion parallel to the basal plane-α a
The CTE, α a , is based on equation (2.12). As before, in deriving an expression for the CTE perpendicular to the basal plane, the second term in equation (2.12) can be ignored by reasoning that S 13 is relatively small, thereby making this term insignificant. Thus, the thermal strains within the basal planes are assumed to be independent of those between layers.
Thus equation (2.12) becomes Noting that the first terms inside the square brackets for the three vibrational modes given in equations (3.6) are dominant, it was therefore assumed that the other terms could be ignored. This assumption greatly simplifies the mathematics. The vibrational mode terms given by equations (3.6) become The anharmonic terms can then be easily calculated for the three modes, giving The integral on the right-hand side of equation (5.1) is accomplished using various substitutions, defining Debye characteristic temperatures for the longitudinal, transverse and out-of-plane directions, θ L , θ T and θ o , as Treating each term individually we have that with the in-plane longitudinal modes giving α a 1 = − πρk(S 11 + S 12 ) The expression for the in-plane transverse modes, α a 2 , is identical to the expression for α a 1 , with C 11 replaced with C 66 = 1 2 (C 11 − C 12 ),  For the out-of-plane transverse modes, we have It is possible to make the simplifying assumption [6] that including contributions from both the in-plane acoustic and the higher frequency optical mode can be made by extending the upper integral limit as follows: The square root arises from the particular wavenumber-frequency relationships for the in-plane mode. The contributions that each of equations (5.7), (5.8) and (5.10) makes to the total CTE, assuming that C 11 −1 ∂C 11 /∂e xx = C 66 −1 ∂C 66 /∂e xx = δ −1 ∂δ/∂e xx = 1, are given in figure 5. Solutions are included for acoustic (θ L , θ T , θ o ) and acoustic plus optical (θ M , θ Q , θ N ) modes. The resultant curves are slightly different from those previously published [6], which may be related to some errors in the previous equations. Generally, the contributions to α a from the longitudinal and transverse modes are similar in magnitude, while the magnitude of the outof-plane mode is larger, as would be expected. In all cases, the acoustic plus optical is smaller in magnitude in the lower to mid-temperature range (between 0 and ∼ 1500 K) than for the acoustic alone. This effect is larger in the case of the out-of-plane mode. These differences indicate that combining the optical and acoustic modes results in less excitation in the lower to midtemperature ranges. Above ∼ 1500 K, the contributions from the two cases converge. This has implications for the prediction of the basal plane CTE α a as discussed and shown in figure 6.
Kelly [6] .10) given in [6]. In this paper, to obtain a reasonable fit by eye, it was found necessary to use values of −8.2, −8.2 and 5.47, respectively. Data from various references for α a , as discussed in §4 of this paper, are given in figure 6 along with the prediction of α a obtained using equation (5.5).
Applying these values and summing the three expressions as equation (5.5), the following comparisons with experimental data are obtained.
The fit to the data in figure 6 is very similar to that given by Kelly [6], although there is significant scatter in the data; a much better experimental dataset would be required to justify further optimization of equation (5.5). The initial negative CTE at low temperatures arises from the interaction of the out-of-plane anharmonic modes with the two in-plane modes, whereas at temperatures above about 300 K the in-plane modes start to dominate and α a becomes positive around 675 K. With reference to figure 5, making the approximation of ignoring the optical contribution gives a better fit to the data at low to mid-temperature ranges (between 0 and 1500 K) and makes little contribution to the higher temperatures, above 1500 K.

Graphite crystal CTE model-Lennard-Jones approach to deriving a relationship for α c (a) The relationship between lattice d-spacing and temperature
To develop the methodology further, a relationship between lattice d-spacing and temperature is required. This is obtained using fits to experimental data obtained from Nelson & Riley [29], These fits are shown in figure 7. The average fit implies that the d-spacing at absolute zero and room temperature has values of 0.3328 nm and 0.3354 nm, respectively, compared with the usually accepted value of 0.335 nm at ambient room temperature [14,15]. It is of course possible to determine an equilibrium value of d from a calculation of the thermodynamically stable crystal structure. We, however, decided to take an average from experimental data because, as will become apparent later on, this lattice spacing temperature dependence is crucial for better matching to experimental data at high temperatures.

(b) Lennard-Jones approach to deriving a relationship for α c
This approach is outlined in references [10][11][12] and starts with the Lennard-Jones-type relationship for the potential energy between two atoms separated by distance d 0 , where r is the atomic separation, A = 2.43 × 10 −78 Jm 6 and r 0 is a constant. An approximation is required in which the carbon atoms are assumed to be distributed with uniform density within a set of basal planes. This is achieved by integrating equation (6.5) over the graphite crystal lattice to give the equivalent energy per unit volume of the lattice. Following a methodology suggested by Crowell [40] gives where σ = 1/q (q = 2.62 × 10 −20 m 2 is the area per atom in a basal layer [8]) and N 0 = ρN a / M = 1.13 × 10 29 atoms m −3 is the number of atoms per unit volume (graphite crystal density ρ c = 2.26 g cm −3 , Avogadro's number N a = 6.022 × 10 23 atoms or molecules per gram-mole, the molecular weight of carbon M = 12.01 atomic mass units). The radial distance in the basal planes is given by x.  The equilibrium spacing, d 0 , is given by ∂E/∂d = 0, which leads to The elastic modulus for a lattice spacing, d, perpendicular to the basal plane is given by the second derivative of equation (6.7) with respect to strain between the basal planes, i.e.
which gives an equilibrium value of Another important function is the derivative of the elastic modulus with respect to strain, with an equilibrium value of Using the average relationship between temperature and lattice spacing d, given by equation (6.4), the relationships given by equations (6.9) and (6.11) are shown in figure 8. Thus we obtain room temperature values for C 33  Minimizing equation (6.14) gives Noting that α c = (1/d)(∂d/∂T) and pre-multiplying by d before differentiating equation (6.15) with respect to T gives As previously discussed, in the case of deriving α c only the out-of-plane mode v 3 is required. It can also be argued [12] that the only significant term in the vibrational mode is the first one, whereas the only significant term in the anharmonic term is the second one. This was justified by stating that previous numerical modelling showed this to be true for the first model [11], meaning it is reasonable that this would still hold true in this new approach. It should be noted that this approximation cannot hold to zero wavenumber, as will be shown later. Applying this reasoning, equations (3.6) and (3.7) simplify to This simplifies equation (6.16) upon rearranging to exp(x) − 1 2 dx, (6.18) with an upper integration limit defined as The version of equation (6.18) derived [12] mistakenly has a denominator of 4 instead of 8; however, this mistake is compensated for by an error in the lower integration limit as discussed below. The integral in equation (6.18) does not converge for small values of x-the neglected frequency terms in the approximation (equation (6.17)) are important in the small wavenumber regime-it is therefore necessary to define a lower limit. It was originally assumed [11] that the smallest value of σ a at which the out-of-plane vibrations can be regarded as purely two dimensional is at the point where the first term starts to become significant and equal to the second term, leading to  Figure 9. Prediction of α c taking account of the change in d-spacing with increasing temperature, compared with experimental data [27,29,[34][35][36][37] and previous model (equation (4.2)). (Online version in colour.) A lower limit of x o ≈ 0.083x m , θ ≈ 93 K is required to give a good fit to available data. A smaller limit will overestimate the Debye integral whereas a higher limit will underestimate the Debye integral.
Equation (  which is plotted in figure 9 along with the previous model given by equation (4.2) for reference. As before an attempt to include optical modes is included [6], by taking an upper limit of 2θ o . We can again examine the behaviour of the thermal expansion coefficient in the high and low temperature limits by examining the leading order behaviour for θ , θ 0 T and θ , θ 0 T. It is important to remember that θ o > θ .
For low temperatures, θ , θ 0 T, the leading order behaviour is For high temperatures, θ , θ 0 T, the leading order behaviour is We again see that, as the temperature tends to zero, the thermal expansion coefficient exponentially approaches zero. However, at high temperatures, the thermal expansion coefficient only has temperature dependence through the change in lattice spacing, d, with temperature. This temperature dependence is responsible for the up-turn of the CTE at high temperatures, which appears to be in good agreement with the data. Without it, the CTE would level off at high temperatures just like the original continuum coefficient of thermal expansion model.  Figure 10. Prediction of α c , taking account of the change in d-spacing, with increasing temperature and differing amounts of crystal strain, note that the 0%, 0.015% (measured by Marrow et al. [44]) and 0.1% curves overlie each other. (Online version in colour.) Equations (5.5), (4.2) and (6.22) give methods of estimating crystal α a and α c , respectively, from first principles without resorting to atomistic modelling. Furthermore, they allow for sensitivity studies to be performed to determine the influence of crystal strain, modulus, temperature, etc. on CTE. The sensitivity of this simplified methodology is mainly dependent on the ratio C −1 33 (d 0 )(∂C 33 /∂e zz ) d 0 , which is (from equation (6.12)) equal to −17. Kelly & Duff [8] calculated this ratio using a slightly different methodology based on the work of Agranovich & Semenov [41] and Girifalco & Lad [42] to be −7.27 × 10 11 /38.6 × 10 9 = −18.8 and −6.04 × 10 11 /36.5 × 10 9 = −16.5, respectively. Thus, the first ratio of 18.8 would overpredict α c by 10% and the second ratio would underpredict by 3%.

Influence of loading polycrystalline graphite on individual crystal thermal expansion
It has been observed that both mechanical applied loading of unirradiated polycrystalline graphite [2] and irradiation creep strain [43] can induce a significant change in bulk CTE, with tensile strains reducing CTE and compressive strains increasing CTE. The equation for α c derived above can be used to investigate the effect of pressure on the thermal expansion coefficients of graphite crystals [7], assuming a change in interlayer spacing. More recently, Marrow et al. [44] measured the change in lattice d-spacing with loading using neutron diffraction and synchrotron XRD. Taking the room temperature lattice d-spacing from [44] to be 0.3353 nm implies a maximum strain of ∼0.015% associated with a bulk stain of 3000 µ . Assuming a Young's modulus of 11.9 GPa [44], this implies a bulk stress loading of 35 MPa, which is a similar order to that given in [2]. Modifying equation (6.22) to be a function of strain gives  -spacing (a) and a-spacing (b), reported by the various authors from the UKAEA [45][46][47] and from NRG in Petten [48]. (Online version in colour.)  (T) is the lattice thermal strain, which can be calculated using equation (6.4), and (P) is the lattice strain due to bulk loading. Predictions of α c are plotted in figure 10 over the range 20 to 1000 K.
From these calculations, it is clear that the change in crystal strain measured by Marrow et al. [44] is unlikely to account for the change in bulk CTE due to loading. It is also clear that crystal strains greater than approximately 1% are required to lead to significant changes in α c . Thus for practical purposes such strains are unlikely to be generated by bulk loading. However, they could result from fast neutron irradiation as discussed below.

Influence of irradiation on α a and α c
Small changes, less than a maximum of 2% at an irradiation fluence of 25 dpa, in lattice spacing due to fast neutron irradiation measured using XRD have been reported by the United Kingdom Atomic Energy Authority (UKAEA) [45][46][47] and more recently by the Nuclear Research Group (NRG) in Petten, The Netherlands [48]; see figure 11. These data have been used here in equations (5.5) and (7.1) to predict the change in CTE in the a-and c-directions as a function of fast neutron fluence and are compared with measurement of the thermal expansion measured on HOPG irradiated between 300 • C and 650 • C in figure 12. Although there is significant scatter in the measured data the predictions uphold the assumption of invariance in graphite crystal CTE for all practical purposes. However, data on irradiated polycrystalline graphite show a significant reduction in CTE with increasing fast neutron fluence. Possible reasons for this behaviour are discussed in §9.

Irradiation-induced changes to bulk polycrystalline CTE-application and some observations
To determine the bulk CTE, one can consider a simple Reuss model [49], in which all the crystals are arranged in series. For isotropic graphite, this gives a linear CTE of If all the crystallites in nuclear graphite were randomly distributed with no porosity, this would imply an unirradiated CTE of ∼ 8 × 10 −6 K −1 . However, there is porosity particularly perpendicular to the crystal c-axis, which reduces the influence of the crystal c-axis expansion, reducing the overall bulk CTE. The expansion in the c-direction is said to fill this accommodation porosity without increasing the specimen length in that direction. Thus, using the predictions of crystal CTE in equation (9.2) and normalizing the results to the virgin value a prediction of irradiated bulk CTE is compared with the measured data for well-graphitized, medium-grained nuclear graphite grade 1 in figure 13. While this does predict some reduction in CTE with increasing fluence as observed in all medium-grained graphite, the reduction is small compared with the experimental observations. The Reuss approximation effectively assumes that the phases are infinite in one direction or are at least continuously connected in one direction. However, analysis of the XRD peak broadening, based on the Scherrer equation [50], showed that the crystal size L a and L c both decreased with increasing irradiation, appearing to start to saturate at higher fluence [48,51]; this is shown in figure 14. This may suggest that the Reuss approximation breaks down in much the same way as it does in composite mechanics with discontinuous reinforcements [52][53][54].
Assuming that the reduction in crystallite size increases porosity these data can be used to try and account for the increase in porosity by weighting equation (9.1) as follows: , (9.2) where L c (0) and L a (0) are the unirradiated values of L c and L a . Figure 13 suggests a combination of the two models is required to fit the data, above and below about 5 dpa. Below 5 dpa bulk CTE is largely invariant to fast neutron fluence, whereas  Figure 15 shows the crystal coherent scattering domain size data [48,51] against irradiated CTE for EU medium-grained nuclear graphite grade 1 and grade 2, clearly showing a correlation between CTE and large coherence crystal size. These comparisons suggest that the reason for the fall in CTE above 5 dpa is caused by irradiationinduced changes in crystal perfection, leading to disruption of the crystal lattice, and resulting in a smaller scattering length.      acoustic acoustic 1/C 11 .∂ C 11 /∂e xx × 1.1 acoustic 1/C 11 .∂ C 11 /∂e xx × 0.9 acoustic and optical acoustic and optical 1/C 11 .∂ C 11 /∂e xx × 1.1 acoustic and optical 1/C 11 .∂ C 11 /∂e xx × 0.9 acoustic acoustic 1/C 66 .∂ C 66 /∂e xx × 1.  Figure 17. Sensitivity of α a to simplified anharmonics C 11 −1 ∂C 11 /∂e xx , C 66 −1 ∂C 66 /∂e xx and δ −1 ∂δ/∂e xx [27,31,[34][35][36].