Multi-scale modelling of nuclear graphite tensile strength using the site-bond lattice model

Failure behaviour of graphite is non-linear with global failure occurring when local micro-failures, initiated at stress-raising pores, coalesce into a critically sized crack. This behaviour can be reproduced by discrete lattices that simulate larger scale constitutive responses, derived from knowledge of microstructure features and failure mechanisms. A multi-scale modelling methodology is presented using a 3D Site-Bond lattice model. Microstructure-informed lattices of both ﬁ ller and matrix constitu- ents or ‘ phases ’ in Gilsocarbon nuclear graphite are used to derive their individual responses. These are based on common elastic modulus of “ pore-free ” graphite, with individual responses emerging from pore distributions in the two phases. The obtained strains compare well with experimentally obtained data and the stress-strain behaviour give insight into the deformation and damage behaviour of each phase. The responses of the ﬁ ller and matrix are used as inputs to a larger scale composite lattice model of the macroscopic graphite. The calculated stress-strain composite behaviour, including modulus of elasticity and tensile strength, is in acceptable agreement with experimental data reported in the literature, considering the limited microstructure data used for model's construction. The outcome supports the applicability of the proposed deductive approach to the derivation of macroscopic properties. © 2016 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Synthetic graphite has been used in the nuclear industry as a fast neutron moderator since the first demonstration of a chain nuclear reaction in the 1940s. Its retention of strength at elevated temperatures has allowed it to be used as a structural component in high temperature gasecooled reactor designs [1]. Nuclear graphite can be regarded as having a three-phase microstructure, which depends on the raw materials and the manufacturing process used. Filler particles, which derive from calcined petroleum or pitch coke, are dispersed within a matrix of binder material, usually consisting of graphitised coal-tar pitch mixed with finely ground filler particles. Both of these solid phases host populations of pores, with sizes that cover the length scale from a few nm upwards to mm [2]. Graphite belongs to the class of quasi-brittle materials, which exhibit limited non-linear stress/strain response prior to maximum or peak stress, with a macroscopic effect akin to plasticity [3]. This non-linearity is partly attributed to the generation and growth of micro-cracks, which occur at length scales dictated by the prominent microstructure features. As tensile strain is applied, microcracks initiate around the larger pores due to local stress amplification, and their effect leads to a reduction in graphite's stiffness [4]. Ultimate tensile failure results from crack growth and coalescence into a flaw of critical size [4]. The continued evolution of these processes determines the behaviour beyond the peak stress; graphite may exhibit a limited post-peak softening or fail at peakstress, depending on its microstructure and the loading conditions [3]. Therefore, the structural integrity of nuclear graphite is controlled by the organisation of the three phases (filler, matrix and pores) and the component's service conditions.
Conventional modelling strategies, such as the finite element method [5], assume the behaviour within model elements is homogeneous and scale-independent, which inherently fails to account for the effects of microstructure failure mechanisms in the material response. This might be inappropriate in certain cases, particularly when modelling the graphite responses at length scales close to its microstructure features. If fracture mechanics is to be used in the structural integrity assessment of nuclear graphite, the fracture process zone ahead of a macroscopic crack is of specific interest, as this is influenced by microstructure features and specimen size [3], [6]. Integrity assessment requires local information in such cases.
Local effects have been incorporated into various modelling approaches. This can be done either statistically, using the weakest link assumption [7], or within a continuum framework extended by cohesive zone models [8]. However, these approaches are restricted by a dependence on phenomenological calibration against macroscale data. Lattice models are a branch of discrete models whereby nodes are connected by elements into a statistically parallel network. Such models have been developed for quasi-brittle materials including graphite [7] [9], and concrete [10e11], incorporating microstructure information into the element properties to allow the simulation of macro-scale behaviour as a consequence of micro-failure mechanisms. Lattice models can be constructed as irregular, representing specific, usually imaged, microstructures, or regular. Irregular models present a substantial problem, as the calibration of lattice element properties with measured continuum properties becomes a trial and error process. Regular models can be calibrated in many cases analytically and allow for up-scaling to potentially representative volume elements, although one issue with most 3D regular lattices has been that they could not be tuned to reproduce desired Poisson's ratio values [12]. The regular 3D sitebond model proposed by Jivkov and Yates [13] is capable of reproducing the range of Poisson's ratio required for quasi-brittle materials, and this model has been further extended in this research.
The original model, using beams as lattice elements, was firstly applied to study damage evolution in concrete in tension [14] and under complex loads [15]. One problem with that formulation is that the structural beams between nodes introduce local micropolar behaviour in the lattice, i.e. rotation-dependent energy conservation. That behaviour has not been confirmed experimentally and consequently the calibration of beam properties remains incomplete. To avoid this, the site-bond model has been reformulated with spring bundles, the stiffness coefficients of which have been calibrated analytically [16]. The work presented here makes a substantial step in the development of the site-bond model with a new approach to represent lattice elements and their behaviour. Experimentally-derived distributions of microstructure features [17] are used to inform separate site-bond lattice models for filler and matrix constituents (i.e. phases) in the microstructure of nearisotropic Gilsocarbon graphite. The results are compared with experimental data of the tensile deformation of the individual constituents, obtained recently by strain imaging of the microstructure during a mechanical test of the same graphite [18]. Brief details of that experiment and the data obtained are provided here. The experimentally measured response of the separate microstructure constituents is then used to inform a multi-scale sitebond methodology. Such a model may complement and scientifically underpin the conservatism of structural integrity assessment methodologies for graphite, providing size estimates for areas of significant local damage ahead of macroscopic cracks, or damage evolution laws for use in continuum scale models.

Theory and method
Within the site-bond methodology [13], material volume is represented with a discrete assembly of truncated octahedral cells illustrated in Fig. 1(a). This choice of shape is in accordance with statistical studies that demonstrate its suitability for representing a topologically averaged microstructure [19]. The computational counterpart of the regular assembly is a 3D lattice, or mathematical graph, consisting of sites at cell centres, connected by bonds to 14 neighbouring sites. This yields two distinct bond types, B 1 and B 2 , which emanate from sites in the 6 principal (normal to square faces) and 8 octahedral (normal to hexagonal faces) directions, respectively, as illustrated in Fig. 1(b). The bond types have lengths L and ffiffiffi 3 p L/2, with L representing the cell extension in the principal directions. The bonds' behaviour is associated with an inter-cell volume, called the support volume, formed by the two pyramids with common base at the face normal to the bond; Fig. 1(c) shows the support volume of a principal bond. An example of a site-bond model is given in Fig. 1(d).
Previously the site-bond model has been applied to the threephase graphite microstructure, i.e. matrix, filler particles and pores, with filler particles considered to be located at sites and to occupy fractions of cell volumes, pores considered to be located in some support volumes, and the remaining volume occupied by matrix [20]. In this manner, the network of bonds, modelled as bundles of independent springs, one axial and two transversal, represent the potential micro-failures both within and between particles. With the mapping of particles to sites, the model length scale, L, has been determined from experimentally measured particle size distribution and volume fraction. With known L, the springs' stiffness coefficients have been calculated from the energy equivalence between discrete and continuum cells under homogeneous strain fields according to the procedure outlined by Zhang et al. [16]. This methodology, derived specifically for the site-bond geometry, allows for accurate representation of isotropic elastic materials with Poisson's ratio ranging from À1 to 0.5, an improvement on previous lattice arrangements where only zero Poisson's ratio has been allowed [12].
However, when the observable microstructure cannot be considered as a three-phase composite, but rather a two-phase composite with pores dispersed in a solid, the model length scale cannot be calculated in the same manner. In such a case (e.g. Ref. [21]), the model length scale is arbitrary, and similarly to the finite element analysis, improved accuracy is achieved by reducing the scale, i.e. increasing the number of cells in the assembly representing the specimen. Further in Ref. [21] the bonds have been represented by 1D connector elements in ABAQUS [22], rather than spring bundles. The combination of connector elements with only axial non-linear and dissipative physical response, and a geometrically non-linear formulation, i.e. finite deformation analysis, makes the representation more physically realistic. Firstly, springs represent only conservative behaviour, similar to inter-atomic potentials in molecular dynamics, while connectors allow for energy dissipation. Secondly, the finite deformation analysis ensures conservation of angular momentum at sites in accordance with recent advances in geometric theory of solids [23].
The model presented in this paper combines and extends the developments presented in Refs. [20] and [21]. The graphite is considered as a three-phase composite and the approach used in Ref. [20] is applied at the composite level. However, the composite level properties are derived from separate models of filler particles and matrix at the constituent level, where the approach used in Ref. [21] is applied. The procedure for bond calibration follows [16], with the exception that only the axial stiffness coefficients are used for the connector elements. The axial stiffness coefficients of principal and octahedral bonds are given by Equations (1) and (2) respectively, where E and n are macroscopic elastic modulus and Poisson's ratio.
In the absence of transversal springs, the analytical results of [16] dictate an initial macroscopic Poisson's ratio of 0.25 must be used in order to maintain energetic equivalence in the calibration procedure. However, the finite deformation analysis reduces the emergent ratio to the prescribed value used in Eqns. (1) and (2) such that the actual value for graphite of 0.2 may be used. Further, local heterogeneity of graphite solid phases due to arrangements of crystals and the presence of unresolved porosity is represented by variable E for different bonds, following uniformly random distribution within ±10% of a nominal value. This is in addition to the stiffness changes resulting from resolved pores within each phase, the introduction of which is described in Section 2.3.
In summary, the present work introduces a two-scale approach for graphite modelling, where filler particles and matrix are modelled separately as two-phase materials (pores dispersed in solids) and their responses are used to inform a larger scale twophase model of graphite (filler particles and matrix).

Pore-free bond behaviour
The bond response follows a linear relationship in compression and a linear-softening relationship in tension as shown in Fig. 2. The tensile behaviour encapsulates both the deformation of the bond support volume, V, by storing elastic energy, and the failure of the face between the two cells with area A, by dissipating energy in surface generation. The energy released upon bond failure, G C , is the sum of the elastic energy stored at face failure initiation, G E , and the dissipated energy in full face separation, G D . The force, F d , and the displacement, U d , at face failure initiation can be calculated from known G E and the bond stiffness. The failure displacement, U f , can be calculated from known G D and the failure initiation point.
In a previous work on graphite grades IG110 and PGX (both nearly isotropic) [20], the total released energy was equated exclusively to the face separation energy, gA, where g is the enthalpy for creation of two surfaces in graphite, derived by atomic scale calculations to be 9.7 J/m 2 [24]. This did not allow for determination of F d , U d , and U f , from bond stiffness and separation energy alone and required an assumption for the ratio G D /G E , alternatively U f /U d . Irrespective of the selected ratio, the model was not able to predict correctly the relative tensile strengths across different graphite grades without grade-dependent factors (i.e. microstructure dependent factors) to increase the released energy. The average filler particle size provides a significant three-fold difference between the two grades [25]. Therefore, it can be deduced that the grade-dependent factor should be related to a volumetric term, whereby the additional released energy is associated with the stored energy in support volumes.
Hence, here a separate scaling for G E and G D , by the support volume and face area, respectively, is proposed via volumetric and surface constants, U and g: While the value of g is the same as before, the value of the volumetric constant, U, is not as easily derived from experimental or atomic scale calculations. It accounts physically for the volumetric deformation of the bond support volume. This arises not only from the change in bond length but also from the necessity to maintain solid unbroken cells, i.e. support volumes from the same cell remain in contact. The procedure for calibrating this constant is described in Section 2.4.

Material and microstructure
The material considered is moulded IM1-24 Gilsocarbon (GCMB grade) polygranular nuclear graphite, manufactured by Graftech (formally UCAR). The bulk material has weakly-anisotropic properties; depending on orientation, the elastic Young's modulus is between approximately 11.6 and 11.9 GPa, with a Poisson ratio of 0.2 and a tensile strength between 19 and 20 MPa [26]. It is one of the graphite grades used in the nuclear cores of the UK Advanced Gasecooled Reactor fleet. The same grade, from different billets, has been studied in previous work by some of the authors [27e29].
High resolution computed X-ray tomography data were obtained with a voxel size of 1.8 mm in experiment EE9036 at the Diamond Light Source (I12 beamline). Full details of the experimental conditions and standard back-projection tomographic reconstruction from radiographs are reported elsewhere [18]. The imaged volume discussed here (4.32 Â 4.32 Â 4.81 mm) contains filler, matrix and pores. The pores have lower X-ray attenuation; there is also some incidental phase contrast due to the imaging conditions, which aids the detection of pores. The regions of filler and matrix can then be identified by the morphology of their pores; filler particles exhibit a characteristic onion-skin structure of lenticular pores, and the matrix has a less organised structure, Fig. 3. A total of 55 filler particles and 25 matrix sub-volumes of different sizes were extracted from the dataset. The smallest volume of a selected filler particle is 0.05 mm 3 , and the largest volume is 4.4 mm 3 . The matrix volumes vary from 0.16 to 1.16 mm 3 . The results for both phases are shown in Fig. 4.
The microstructures may be segmented using an image intensity threshold to define the pores and solid graphite. All X-ray tomography images were converted to 8-bit datasets before the segmentation. It was not possible to apply a single threshold for the grey-scale dataset, so the segmentation procedure was performed in ImageJ software using a multi-step thresholding with the corresponding smoothing and binarisation steps for pore boundary determination and large pore filling. The thresholds were verified visually by comparison with the original grey-scale image. The filler particle boundaries have been manually identified by using the visible matrix pores, which surrounded a particle. These pores are quite large and have a well distinguishable structure that is different from the lenticular pores of the particle. Parts of the particle boundary connected with solid matrix were restored then assuming the ellipsoidal shape of the particle. The shape of the unbroken Gilsocarbon filler particles is typically ellipsoidal, often close to a spherical shape. Analyses of tomographic data from the same graphite billet found that the fraction of filler particles varies within 14e29%; the fraction in the volume from Ref. [18] is 29%.
Within both phases the smallest pore volume that could be resolved was restricted by the resolution of the tomography data. The mean filler pore volume throughout the imaged volume is 15080 mm 3 (standard deviation 6320 mm 3 ). The pores in the matrix have very different shapes and cover a wide range of volumes from 6 to 10 5 mm 3 ; the most frequent pore volume is approximately 100 mm 3 . The largest individual pores observed occupy a volume of about 10 6 mm 3 in the filler and 10 5 mm 3 in the matrix; more than 90% of the total pore volume in a subset may be spatially combined into one large pore. It is important to note that the selection of matrix volumes excluded regions that contained larger pores (>100 mm), which occur due to gas porosity. These pores have been quantified using laboratory tomography data of lower resolution (Skyscan 1272) (10 mm/voxel), their fraction has been estimated as 6.3% of the total volume of the sample.
The porosity fraction in the different phases was extracted for a randomly selected subset of 20 filler particles and 7 matrix volumes from tomographed volumes of the same graphite billet [17]. The observed porosity fraction in the filler is typically lower than in the matrix; the mean porosity of the filler subsets is 12.2% with a standard deviation of 3.6%, while the matrix has average porosity of 16% and a standard deviation of 3.1%. The cumulative probability for porosity observed in both filler and matrix phases are shown in Fig. 5(a). A region was chosen at random for each phase, so that its pore size distribution could be used for the subsequent simulations. The pore size distribution from the selected region for both phases is shown in Fig. 5(b) with the filler and matrix samples containing 1204 and 24394 pores respectively.
The deformation of filler and matrix, up to an applied tensile stress of 7.5 MPa, has been studied by digital volume correlation (DVC) of X-ray computed tomography images, obtained during a tensile test. Full details of test and the image correlation analysis are reported elsewhere [18]. Briefly, the DVC analysis of subvolumes that contained filler and matrix was used to calculate the axial strains in each xy-plane as the gradient of the average vertical displacements in the z-direction, which corresponded to the tensile axis. Only those displacements contained inside the ellipsoidal volume that defined the filler particle have been considered for the filler. In the studied volumes, the axial strain at 250 MPa applied stress varies from 0.0003 to 0.001 in the filler and from 0.0005 to 0.0013 in the matrix; the average axial strain is larger in the matrix (774 mε ± 178 mε) than for the filler particles (667 mε ± 197 mε), which suggests the elastic modulus of the matrix may be lower. The mean strain of the tomographed volume of the tensile sample at 250 MPa was measured to be 730 mε.

Pore-affected bond behaviour
The mapping of microstructure to the model follows the procedure outlined in a previous work [21] whereby micro-cracks are considered to initiate at pores. Pores, with sizes selected at random from an experimental pore size distribution, are assigned to faces of  cells until the desired porosity is reached. The presence of porosity is reflected in changes of the tensile response of corresponding bonds, Fig. 6. The peak force changes from the pore-free value, F d , to a new value, F' d , according to where V is the support volume of the bond and V 0 is the support volume remaining after the corresponding pore volume is removed, i.e.
In the same manner, the displacement at face failure initiation, U d , changes to a new value, U' d , according to Equations (4) and (6) represent pore-corrected force and displacement parameters via pore-corrected (or effective) areas and lengths, respectively.
The stored elastic energy at face failure initiation, G E , scales with the change of support volume: Differently, the damage energy, G D , scales with the face area: where A 0 is the face area remaining after the corresponding pore area is removed: In this manner, the failure displacement, Uf', reduces proportionally less from the pore-free value than the initiation displacement, Ud', so the amount of softening is reduced as the pore size is increased, resulting in an increasingly brittle response. The bond stiffness is also reduced.
Similarly to tensile behaviour, the presence of porosity also alters the compressive response of the bond. Compressive stiffness is reduced by the same factor as tensile stiffness according to pore size, although this decrease is only maintained for a relative displacement equivalent to the diameter of the pore present, after which time the stiffness increases back to the original value.
In the process of random pore allocation to faces, some pore volumes may exceed the corresponding bond support volumes. Such bonds are removed from the model and the excess pore volume, i.e. the difference between the allocated pore volume and the removed support volume, is distributed to neighbouring bonds in the same manner until all the volume is allocated. In this manner the size distribution of pores distributed to faces will be a representative sample of the experimental distribution, although the spatial distribution of pores will be entirely random.

Calibration of the volumetric constant
Preliminary studies using the pore representation outlined in Sub-section 2.3 were undertaken to calibrate the volumetric constant. Four different grades of nuclear graphite, IG110, NBG-18, PGX and Gilsocarbon, were simulated with models scaled according to the size and volume fraction of the corresponding filler particles following [20]. The microstructure information used for grades IG110, NBG-18 and PGX, including pore size distributions, porosity and filler particle sizes was taken from the microscopy studies by Kane et al. [25]. The filler particle size distribution used for the IM1-24 Gilsocarbon was that shown in Fig. 4. The microstructure data used for all grades and the corresponding references are  Table 1. A filler particle volume fraction may vary within graphite grades. The value used was 0.2, an average of the two values obtained of 0.144 and 0.252 presented in Ref. [17]. Although these values were specific to Gilsocarbon, the same filler particle volume fraction of 0.2 was used for all grades. The pore distributions taken from Ref. [25] do not differentiate between pores found in the matrix and filler phases. Hence, for simplicity only the pore size distribution for the matrix phase, shown in Fig. 5, was used to calibrate the volumetric constant of the IM1-24 Gilsocarbon. It is not necessary to repeat the process for the filler phase, because for a prescribed porosity its larger pores can be considered as represented in the model by the coalescence of numerous smaller pores assigned to one and the same lattice bond.
Measured Young's moduli, tensile strength and typical porosities for each grade in its virgin (i.e. as supplied, without any effects of fast neutron irradiation or radiolytic oxidation) are shown in Table 1. To calculate the pore-free, axial stiffness coefficients of bonds, pore-free values of the Young's moduli are required. These were calculated with a series of normalised simulations with Young's modulus equal to one, without porosity and with virgin state porosity. Several realizations with the latter were analysed, differing in the spatial distribution of pores but identical size distributions. The simulations were performed without failure of bonds, i.e. in the elastic regime of bond behaviour. From these simulations the ratios between the (unit) pore-free and the simulated average virgin-porosity moduli were calculated. The ratios were used to scale the experimental virgin-porosity moduli to pore-free moduli, reported in row 3, from where pore-free axial stiffness coefficients were calculated by Eqns. (1) and (2). The close proximity of the calibrated pore-free modulus of Gilsocarbon, 14955 MPa, to a pore-free value of 15 GPa, derived from nanoindentation experiments [30], gives confidence in the calibration procedure. In the absence of equivalent data (to the knowledge of the authors') for the other grades this is considered an adequate validation.
Models with calibrated coefficients were subject to displacement controlled uniaxial tension until failure, using the assumption that energy released at bond failure equals the energy of face separation, gA, specifically assuming as previously [20] [21].
Failure was considered to be the point at which the simulation failed to find equilibrium using a time increment size less than a threshold value (1 Â 10 À25 was deemed suitably small). Several control parameters were adapted to improve convergence in the highly non-linear model. Viscous regularization was utilized to improve the dissipation of energy from damaged bonds to the surrounding bonds. A damping coefficient of 1 Â 10 À5 was chosen after initial studies showed it gave the best compromise between aiding convergence and producing a consistent peak stress. The quasi-Newton method was used for the analysis and the method of extrapolation was suppressed, preventing excess iterations.
The difference in simulated and experimental values of tensile strength, s Sim T and s Exp T , as a function of model cell size (volume), for a model that does not include the proposed volumetric correction, is shown in Fig. 7 with red marks. Here, the cell size, shown in Table 1, reflects the different structures of the four graphite grades in terms of particle size distribution. It is apparent that, when using the cell size as a representation of the average filler volume, the relationship between strength discrepancy and size is approximately linear, i.e. the larger the average filler particle size the larger the difference between simulated and expected tensile strength.
This relation between strength discrepancy and cell size occurs from using the energy of face separation as a sole measure of energy released upon bond failure, failing to account for the energetic separation into area, G D ¼ gA, and volumetric terms, G E ¼ UV, as proposed in Section 2.1. As such the volumetric constant can be calibrated from the linear trend.
The peak bond force in the preliminary studies (without volumetric term) and the proposed model can be expressed as: with stiffness coefficient, K, being the same for both models. The linear trend shown in Fig. 7 can be expressed as: where m is the gradient of the linear trend and S is the cell size. From comparison between simulations and experiment, one can write: Table 1 Comparisons between grades for model inputs (Young's modulus, virgin porosity), the resulting model cell size and literature values of tensile strength and mean particle size.
As such, by substituting Equations (11) and (12) into Equation (13) the volumetric term can be expressed as: The volumetric term is therefore a function of cell size (volume) and differs for bonds B 1 and B 2 , according to their support volumes, V, and face areas, A. The gradient m, was calculated from Fig. 7 as 0.0022. The units of the volumetric term are J/m 3 . Rerunning the same simulations with the calibrated volumetric term produced results with considerably less discrepancy from experimental values, as shown in Fig. 7 with blue marks.

Single phase modelling and results
In this work a two-scale methodology is introduced in order to build up the composite response of graphite directly from the mechanical response of the individual phases. In this section the single phase procedure used for both filler and matrix phases is outlined; specifically, the modelling of filler particles and matrix incorporates the experimentally measured pore size distributions of the individual phases.
Five site-bond models for each phase were generated with porosities randomly selected from the measured porosity distributions shown in Fig. 5(a). Within each model, pores were randomly assigned to faces with sizes from the corresponding measured pore size distributions, Fig. 5(b), until the required porosity for the particular model was achieved. All models were constructed as lattices occupying cubic regions of 10-cell sides for computational efficiency. Fig. 8 and Fig. 8(b) show the stress-strain response obtained from the filler models and matrix models, respectively. The response is visibly different for both phases with significantly more energy dissipation and nonlinearity exhibited by the filler phase simulations as a result of bonds entering the softening region of the constitutive behaviour. The matrix phase shows less pre-peak nonlinearity with sudden "avalanche" failure shortly after peak stress. One of the five matrix simulations failed to run (results not shown), presumably as a result of multiple failures occurring in the initial solution increment. The stress-strain curves illustrate the effect of variable porosity on the responses of different phases, which will be used as input to the composite level model in Section 4. There are significant variations of elastic modulus within the models for each phase, which do not relate simply to the total porosity. This suggests that the response results from both the porosity value and the different spatial distributions of pores across samples. It appears that porosity increase generally leads to elastic modulus reduction with this trend more prominent within the filler phase. However, this is not a comprehensive trend with samples of comparable porosities showing different moduli, which is attributed to the spatial arrangements of the pores. Furthermore, it has been shown that pore shape affects modulus [36] [37] although this phenomenon is not yet represented in the current model.
The behaviour of each phase can be understood when the initial model states are considered. The initial porosity present on a filler and matrix model are shown in Fig. 9. The brittle response of the matrix phase results from high proportion of bonds that are both removed pre-simulation due to porosity, on average 10.5%, and damaged but not yet failed, on average 59.9%. The high number of damaged bonds explains the catastrophic failure, with a large number of damaged bonds reaching a critical load at the same simultaneously. The more "graceful" failure of the filler phase emerges from lower proportions of the same values, 7.1% and 30.9% respectively, which allow damage to evolve. Table 2 lists the Young's modulus calculated in each simulation for the initial load increment, which was sufficiently small so as no failures occurred. It should be noted that some values of E for the filler phase are higher than the pore-free value used to calibrate the models. This is because of the introduced random distribution of E to different bonds with 10% standard deviation. The average value is within 1.5% of the pore-free value of 14995 MPa, suggesting little effect of the filler pores on its stiffness. In contrast, the matrix porosity has a substantial effect on its stiffness, reducing the porefree value by more than 20%.
Total strains are composed of an elastic part and permanent part arising from the generation and growth of micro-cracks. The stressstrain curves in Fig. 8 were used to extract the total strain for each simulation at a stress value of 7.5 MPa, which can be compared to the experimentally measured strains in filler particles and matrix. The experimental values were measured at global tensile stress of 7.5 MPa. Fig. 10 shows the cumulative probability of measured axial strain in both filler and matrix samples, together with the cumulative probability of strains obtained by simulations. The cumulative probability of simulated strains arises from the five model realizations with different porosities per phase. Both experimental and simulated results for the two phases indicate that for the same global stress, matrix strains are higher than filler strains. The difference in the pore systems of the matrix and the filler, a consequence of the distributions of porosity and pore volumes shown in Fig. 5, results in lower stiffness and lower strength of the matrix. The lower stiffness can be attributed primarily to the higher matrix pore volume fraction, while the lower strength e to the higher propensity to micro-crack generation.
Care should be taken with direct comparison of the experimental and simulated results, since the simulations result from a local stress of 7.5 MPa as opposed to the global 7.5 MPa in the experiment; stiffer regions within a heterogeneous bulk specimen attract an increased amount of the load, leading to stress partitioning (similar to that observed in particulate composites) according to the phase properties and position within the specimen. The comparison does however yield some interesting discussion regarding the stress-state of the samples. Both simulations of filler and matrix phases predict lower axial strain at 50% probability than the experimental data; the filler phase model has a larger discrepancy (25% difference from experiment as opposed to 15% for matrix phase). However, the lower compliance of the filler phase would suggest that the filler particles experience a larger stress than the matrix under the same applied global stress in experiments. The variability in porosity and pore size distribution in the microstructure is larger than that used to construct the models analysed, which may account for the smaller variation seen in the simulations. The modelling approach is judged to be promising; both phases were calibrated using the same porefree value of E, leaving the resulting responses of the two phases to emerge from the porosity of the microstructures alone.

Composite modelling and results
The lattice for the matrix-filler composite was based on a cubic cellular structure of length 10 cells, giving C ¼ 1729 cells. Particle sizes from the experimental distribution, Fig. 4, were assigned at random to each site. The cell size was calculated from the model   volume as a function of cell length L, Equation (16), the cumulative assigned filler particle volume and the desired filler volume fraction, q F , according to Equation (17). The value of q F was taken as 0.2, an average of the two values obtained of 0.144 and 0.252 presented in Ref. [17].
Bonds, chosen at random, were assigned random properties of the filler samples, derived in Section 4. This process continued until the cumulative total of the bond support volumes matched the desired particle volume fraction. The remaining bonds were assigned random properties of the matrix samples, also derived in Section 4. Specifically, the stress-strain behaviour of the different filler and matrix samples (i.e. Fig. 8) were used to inform bond behaviour. The local gradient of the non-linear stress-strain behaviour (i.e. from Fig. 8) was used as input to the stiffness calibration used for the single phase models, Equations (1) and (2). This stiffness and the corresponding strain from the sample behaviour were used to calculate the force-displacement relationship of each bond.
The results of the composite simulations are shown in Fig. 11. The initial Young's modulus value of 12900 MPa is higher than typical literature values for virgin Gilsocarbon. This may be a result of failing to simulate the largest pores in the matrix phase (with size >100 mm); these occupy 6.3% of the total volume and were excluded in the distributions of porosity in filler and matrix, reported in Vertyagina and Marrow [17], that were used in the simulations. This extra porosity, which could be included in future models either in the matrix phase models themselves, or as additional porosity applied to the composite model, would slightly decrease the elastic modulus. The tensile strength value obtained of 13.6 MPa is lower than the value of 19e20 MPa quoted in the literature [26]. The simulated tensile strength would be expected to decrease with inclusion of the larger pores, so this suggests that the model does need further refinement, particularly with regards to the peak stress values. Nonetheless, the comparison between simulation and experiment is encouraging at this stage of the model development.

General discussion
It was found that the multi-scale model was more numerically stable and less computationally expensive than the single phase models. This is because the multi-scale model does not require the pores to be modelled explicitly. Instead the effects of pores are homogenised within the models of the individual phases, which are represented by a small number of possible stress-strain behaviours that are reasonably consistent with experimental behaviour (i.e. Fig. 10). In doing so there is no need to remove bonds prior to simulation as a result of porosity. This makes finding an initial equilibrium less expensive and hence improves the numerical stability of the model.
Despite these promising results there are still limitations with the model, mainly linked to limited data and numerical controls. Firstly, in addition to the large pores that are currently neglected, there is finer scale porosity in the graphite microstructure [38] that is unresolved using the experimental techniques described in Section 2.2. Its effect should be included in the pore-free modulus, however. Furthermore no consideration has been taken of residual stresses, which exist following manufacturing [39] and so affect on the stress-state of each phase. With respect to numerical stability, the controls required to obtain convergence, including viscous regularization, and their impact on the model response need to be further understood to increase confidence in the results and hence allow for reduced conservatism in safety assessments of graphite components.
Further work includes an improved calibration procedure compatible with theories of discrete elasticity [23]. Following this, more rigorous studies may be undertaken using the Site-Bond methodology, deriving damage evolution laws and characterising the size of the fracture process zone for use in continuum scale analyses for structural integrity assessment of graphite.

Conclusions
A multi-scale modelling methodology is presented, whereby microstructure-informed Site-Bond lattice models of both filler and matrix phases are used to construct a larger scale composite Site-Bond model of the nuclear grade graphite Gilsocarbon. A key feature of the proposed model is that it requires a single calibration of the elastic properties of "pore-free" graphite, from where it can predict the elastic properties of real graphite from the knowledge of microstructure characteristics, such as particle and pore density and size distribution. The single-phase model results suggest that the evolution of damage is more prevalent in the filler phase than the matrix phase with filler models demonstrating "graceful" stress-strain behaviour resulting from micro-failures as opposed to the brittle failure seen in the matrix phase. Filler particles are shown to be stiffer than the matrix phase.
The modulus and tensile strength value calculated from a multiscale composite model, 12.9 GPa and 13.6 MPa respectively, informed with the responses of the single-phase models, are encouraging when compared to the values found in the literature, 11.6 GPa and 19e20 MPa respectively. The proposed semiempirical approach, using a calibrated pore-free stiffness and deriving longer scale behaviour from microstructure information, has the potential to develop into a deductive methodology for calculating emergent behaviour, when combined with richer microstructure information and validated by damage characterisation experiments. Fig. 11. The stress-strain response of the multi-scale simulation