Image based in silico characterisation of the effective thermal properties of a graphite foam

Functional materials' properties are influenced by microstructures which can be changed during manufacturing. A technique is presented which digitises graphite foam via X-ray tomography and converts it into image-based models to determine properties in silico. By simulating a laser flash analysis its effective thermal conductivity is predicted. Results show ~1% error in the direction the foam was ‘grown’ during manufacturing but is significantly less accurate in plane due to effective thermal conductivity resulting from both the foam's microstructure and graphite's crystalline structure. An empirical relationship is found linking these by using a law of mixtures. A case study is presented demonstrating the technique's use to simulate a heat exchanger component containing graphite foam with micro-scale accuracy using literature material properties for solid graphite. Compared against conventional finite element modelling there is no requirement to firstly experimentally measure the foam's effective bulk properties. Additionally, improved local accuracy is achieved due to exact location of contact between the foam and other parts of the component. This capability will be of interest in design and manufacture of components using graphite materials. The software used was developed by the authors and is open source for others to undertake similar studies. Crown Copyright © 2018 Published by Elsevier Ltd. All rights reserved.


Introduction
Components designed for the nuclear sector (both fission and fusion) are required to withstand challenging environments. Depending on the component's function they may be exposed to extreme environments such as high levels of radiation, temperature, pressure etc. [1]. Future generation nuclear power plants aim to increase output and efficiency over their predecessors, a byproduct of which is even more extreme physical operating conditions [2]. In addition to improving engineering concepts; a number of novel materials are proposed which have properties tailored for the specific requirements, these fall into structural [1] and functional [3] categories. In order to achieve their functional specification these materials are often highly anisotropic. Such materials include; foams [4], high entropy alloys [5], MAX phase-materials [6], composites [7], functionally graded materials [8] and nanograined materials [9].
In the most common fusion power device, the tokamak, helium, the waste product of the fusion reaction is removed at the divertor by directing plasma along magnetic field lines to strike the divertor target plates. In normal operation this region will experience thermal loads of around 10 MWm À2 , as the plasma particle kinetic energy is deposited over the target region [10]. To absorb this energy whilst remaining within operational temperature limits, the divertor design includes active cooling through pipes (built of CuCrZr) protected by tungsten (W) armour monoblocks [11]. Fig. 1 shows the location of the divertor region within a tokamak and an example monoblock. The armour is bonded to the pipe to maintain thermal conduction, but a large thermal expansion coefficient mismatch between the W and CuCrZr causes high levels of stress within the part. Therefore, a functional interlayer is used at the material interface to create a bond between the pipe and armour with improved longevity. ITER, under construction in France, uses a compliant pure Cu interlayer to relieve stress [12]. In future fusion devices it would be desirable to operate at higher thermal fluxes in the divertor to improve device efficiency, however this would further exacerbate the challenges with thermal stress. Of the various solutions being investigated several candidate designs replace the Cu interlayer with functional materials or geometric constructs [13]. As a result, the performance of these particular designs is highly influenced by micro-scale thermo-mechanical mechanisms.
There is considerable interest in developing capability to engineer micro-structures of graphite foam materials to have specific anisotropic properties capable of addressing this challenge [14]. Typically, the properties of a newly manufactured foam would be characterised by extensive experimental characterisation, requiring many samples to test the range of available orientations [15], [16]. This is usually more time consuming and expensive than the process by which the new foam structure can be created by varying the manufacturing parameters. This is a barrier for rapid development of new graphite foams when attempting to 'tune' the properties to match the desired function of the material. This work considers an alternative method of characterisation by using in silico simulation methods which have the potential to replace experimental methods to reduce time and cost. This is achieved by the use of image-based finite element methods (IBFEM) whereby an X-ray tomography image of the manufactured foam is converted directly into a high-resolution FE model which inherently accounts for the graphite foam's complex microstructure. By simulating the experimental characterisation test (typical of those conducted in the laboratory) we may apply the same conditions and observe the same material response in order to derive a measured property. This may be possible for scenarios where the bulk (or parent) material properties are already known and they are affected significantly by microstructural features caused by the manufacturing process, e.g. a graphite foam material. The additional benefit is that once a block of foam is digitised via conversion from a tomographic image, many samples, e.g. 'dog-bone' for mechanical tensile testing, can be digitally 'cut' from the same block whose volume would overlap in space. If done experimentally much more material would have to be manufactured to produce a sufficient number of samples, which is a source of variability between tests and an additional expense.
The component which forms the case study presented here is part of a heat exchanger. Because of this, the effective thermal properties of the graphite foam are of particular interest. Previous attempts have been made to estimate the effective thermal conductivity of graphite foams analytically using a unit cell model that considered parameters such as porosity, pore diameters and side length of wall [17], [18]. Although these make some steps towards estimating the effective conductivity they do not account for anisotropy. X-ray tomography has previously been used to collect statistical data about a foam's microstructure to investigate correlations with elastic modulus [19], but no simple relation was found. FE analysis has been used to include some effect from anisotropic microstructure [20], [21]. The disadvantage of these approaches is that they are largely idealised forms of graphite foams which do not account for the large variability in pore size or shape distribution and require a substantial knowledge about the foam's microstructure as model inputs. Image-based modelling to estimate thermal effective conductivity has previously been performed on carbon foams [22] and cellular aluminium [23] and was demonstrated to discover anisotropy. However, due to computational constraints, only relatively small volumes of material could be analysed in these studies. This meant that the method was only suited to materials with a relatively uniform distribution of microstructure. That is, simulations did not include a sufficiently large representative volume to exhibit the effective bulk behaviour. A statistical analysis of many instances of the modelling would be required to capture the effect of variations in microstructure. Also, because the volume of material analysed was relatively small it was not possible to replicate a laboratory experiment in silico, instead approaches based on idealised boundary conditions were used e.g. applying a temperature gradient and measuring thermal flux to calculate conductivity. Using idealised approaches not achievable experimentally can be desirable because numerical analysis will likely have less error in results compared with experimental techniques. However, there are also benefits to emulating experimental techniques through simulation. Firstly, the motivation for this research is for use of novel materials in the nuclear sector. For nuclear regulators there is currently little confidence in the modelling of such materials. Creating a workflow which can be directly related against experimental data collected according to international standards (e.g. ASME or ISO) will facilitate a wider and more rapid industrial adoption. Secondly, this research is part of a wider initiative to develop a 'virtual laboratory' whereby digital twins exist of a suite of materials characterisation apparatus. The initial vision of this approach is to have experiment and simulation running concurrently for improved analysis. Subsequently, an automated virtual laboratory could perform all characterisation tests simultaneously to destruction and still use the part in-service.
There are various experimental methods by which thermal properties can be measured. Laser flash analysis (LFA) is one characterisation test which is increasingly being adopted as the standard method for measuring thermal diffusivity (by which thermal conductivity can be inferred). Because LFA is an established and widely accepted experimental technique we will consider image-based modelling this test on foam for this work to predict its thermal properties in silico.
In this work the graphite foam, KFOAM (Koppers Inc., Pittsburgh, PA, USA), is considered because of its potential use in nuclear applications. This is because its comparatively high thermal conductivity, 240 W/mK (approximately half that of copper), and diffusivity enables rapid heat dissipation in thermally critical applications. In addition to this, its thermal properties are anisotropic, with thermal conductivity reducing by 75% in two of the orthogonal planes. If utilised correctly, the material micro-structure could be designed to favourably direct the thermal flow through the component to reduce thermally induced stresses. Finally, its inert characteristics mean that it is highly resistant to corrosion in challenging environments.
By using KFOAM as an example, this paper investigates the potential of using image-based finite element method to simulate a standard laboratory test. Thus, if the bulk properties of a solid equivalent are already known, it is shown how this technique may be used to characterise effective properties. Here we assess suitability for characterisation of thermal properties, but the same principle applies to other material properties. Initially this is intended to supplement experimental work, once there is sufficient confidence in the technique it could be used to significantly reduce and maybe replace experimental testing for similar materials where microstructure dominates the effective material properties. The paper then presents a case study of how IBFEM may be used to virtually analyse the performance of a newly manufactured material used in its desired manner. This is achieved by digitally 'cutting' the manufactured block of graphite foam to the required dimensions to be included within a component. The more conventional materials within the component are added using standard computer aided design (CAD) to produce a hybrid CAD-IBFEM model, a form of digital manufacturing. This method gives an estimate of the in-service behaviour of a novel material can be achieved without experimental work. Additionally, this should produce estimates with a higher degree of accuracy because the complex anisotropic behaviour of the material is inherent to the model instead of using homogenisation methods. The example used in this paper is that of a divertor monoblock (a heat exchange component) for a tokamak, a type of fusion energy device.

Material
The graphite foam selected for use was KFOAM P1 HD (Koppers Inc., USA). KFOAM is produced from mesophase pitch which is derived from coal tar. The foaming process was devised at Oak Ridge National Laboratory, USA and has been licensed to Koppers Inc [24e26]. The foam is produced at pressure and carbonized at 0.2 C/min to 1000 C then graphitized at 10 C/min in argon to 2800 C with a 2-h soak at temperature. Further details of the manufacturing process can be found in the paper by Klett et al. [27]. The material properties of interest to this investigation, as reported by the material manufacturer 3 are shown in Table 1. The ┴ and ═ symbols denote perpendicular and parallel orientation to the base plane from which the foam was 'grown' during manufacture. The manufacturer informed the authors that the bottom 3.2 mm of the sample would have a different foam structure as an artefact of the manufacturing process. They would normally remove this region by machining before sending to customers, but it was left as part of the sample for the purpose of this study demonstrating the advantage of in silico material 'preparation'.

Methodology
This section details the methodology for (i) three-dimensional imaging using X-ray computed tomography, (ii) processing of CT image data and conversion into digital geometry (iii) finite element mesh generation, (iv) definition of simulation boundary conditions, equation solution parameters and results analysis. The following results and discussion sections follow the same format. Fig. 2 is a schematic showing the workflow of the image-based finite element method and the software used at each step.
3.1. X-ray computerised tomography (CT) scanning X-ray tomography scanning was performed using a Nikon Metrology 225 kV system at the Manchester X-ray Imaging Facility, Research Complex at Harwell, Rutherford Appleton Laboratory, UK. Imaging and reconstruction parameters are shown in Table 2. The voltage and current parameters relate to the electron gun aimed at a target to produce X-rays. The X-ray signal is passed through the specified filter before being incident on the sample and then detector. During reconstruction, beam hardening occurs when polychromatic beams are used and the soft X-rays (lower energy) are filtered by the sample giving the false appearance of a change in attenuation through the sample. That is, the change in greyscale through the material is not caused by a change in the material (e.g. density) but rather an artefact of the imaging technique. Noise reduction is a standard image processing method that utilises smoothing algorithms. The digital filters used to correct for beam hardening and noise range from 0 to 5, with 5 being the strongest correction level.
To avoid the release of carbon dust from the sample it was scanned whilst within a sealed plastic bag. The bag is sufficiently transparent to X-rays that it will not adversely affect the imaging. Reconstruction of the 3D volume from 2D radiographs was completed using CT Pro V3.1.4785.19683 (Nikon Metrology NV, Tring, Hertfordshire, UK).

Image post-processing
CT data output essentially consists of a series of 2D slices, e.g. Fig. 11. The sliced data is made up of pixels with an associated greyscale value, giving information about the extent to which Xrays are attenuated in that point in space. In between the slices the data is interpolated to create volumetric pixels, commonly known as voxels. When the slices are stacked together these represent a 3D volume.
Before being able to create an FE mesh the data must go through a 'segmentation' process to define which regions or ranges of greyscale values belong to the various materials. To carry out this step the open source software Fiji [28] was used which is a specific distribution of ImageJ [29] which includes plugins specifically to facilitate scientific image analysis. The image was first cropped to remove the sample bag, sample holder and outer air. By use of greyscale thresholding, both foam and porous phases were defined resulting in a binarised 3D volume image.
At this stage the binarised images were analysed to investigate the graphite foam's structure. The number of voxels belonging to the foam give its volume because the volumes of the voxels are known. Using this method, it was also possible to quantify the porosity of the material. By using the graphite foam's measured mass and actual volume it was possible to calculate the density of the carbon itself, in addition to the effective density of the foam.
The lower 3.2 mm of the sample was digitally removed (cropped) as this contained the foam with a different structure due to manufacturing method. The external volume was also cropped to leave only the graphite foam (i.e. remove the external sample bag, sample holder and surrounding air). To investigate the anisotropy in the thermal paths through the graphite foam a skeletonisation process [30] was undertaken on a region of interest (200 Â 200 x 200 pixels). The skeletonisation process gradually 'erodes' (or thins) the object until its paths are only one voxel in diameter. When these points are connected they create a centreline. This makes it possible to perform analyses on the thermal paths through the foam e.g. lengths and directionality.
The tortuosity, t, is a ratio of the measured path length, L, to the Euclidian distance, C, (i.e. 'bee line') as shown in Fig. 3 and Equation (1) [31]. Euclidian distances were measured from the centrelines generated in the region of interest using the ImageJ plugin 'Analyze Skeleton' [32]. Thus, the graphite foam's tortuosity was found in a range of planes. By this method it is possible to compare how the tortuosity changes with relation to the foam's orientation as a measure of the anisotropy. The tortuosity is of interest because the effective thermal conductivity of the material is directly related to how 'quickly' heat may transfer from one bounding surface to the other. A longer path may indicate a lower effective thermal conductivity. An 'image analysis method' approach has previously been used to obtain shape factors in a metal foam [33]. This was used in a derivation of an analytical solution that included tortuosity as one element of the approach. In that instance, the foam was isotropic in its micro-structure.

FE mesh generation
To generate the FE meshes the binarised volume images were imported into ScanIP, part of the Simpleware suite of programmes, version 7 (Synopsys Inc., Mountain View, CA, USA). Two sets of meshes were created: a) cylindrical discs as used for laser flash analysis (LFA) experiments and b) fusion divertor monoblock components with a foam interlayer.
Important to the validity of any results calculated by FE is the 'quality' of the elements that make up the mesh i.e. how closely they match the ideal geometry of an equal sided tetrahedron or hexahedron. If elements are highly distorted they are known to cause numerical difficulties, often overpredicting stiffness or  resistivity in mechanical or thermal models, respectively. For image-based models this can be particularly challenging because of the non-idealised freeform geometries. Conventional metrics were used to investigate the mesh quality (e.g. ratio of longest to shortest edge length), looking at the mean results and the worst elements.

In silico experiment: laser flash analysis (LFA)
The LFA cylindrical discs (or pucks) were digitally 'cut' out of the graphite foam block to have a diameter of 12.7 mm and a thickness of 4 mm. Three discs were cut from each orientation (xz, yz and xy), one from the centre of the block and one each side, offset by approximately 1/4 and 3/4 along the axis of that orientation. Thus, there were nine virtual LFA discs in total, as shown in Fig. 4 (Xlow, Xcentre, Xhigh, Ylow, Ycentre, Yhigh, Zlow, Zcentre and Zhigh). An example is shown in Fig. 5.
In an LFA experiment a laser is pulsed on one surface of the sample while an infrared camera tracks the temperature rise on the opposite surface. If the disc was to be used in an LFA experiment in   its current form the laser would interact with the graphite foam at various penetration depths into the sample because it does not have a solid upper surface. Similarly, the infrared camera would see temperatures at a range of depths through the sample where there is a direct line of sight. In some samples there may be a direct line of sight between the laser and infra-red camera. This would make measuring thermal diffusivity of a foam via LFA extremely challenging. To mitigate this, a method used by Zhao et al. [34] was to sandwich the foam between two solid layers. To emulate this, additional thin solid layers of graphite (0.25 mm) were digitally added to the top and bottom of the discs assuming a bond which is in perfect contact with the graphite foam.

Case study: fusion energy heat exchanger component
For the monoblock, a ring shaped interlayer section (inner diameter ¼ 12 mm, thickness ¼ 1 mm) was virtually cut out of the digital foam cube and joined to CAD versions of the outer armour (22 mm Â 22 mm x 4 mm, central bore with 14 mm diameter bore) and coolant pipe (inner diameter ¼ 10 mm, thickness ¼ 1 mm) to produce a micro-scale accurate virtual monoblock, as shown in Fig. 6.
Once the desired geometries had been constructed they were converted into tetrahedral FEM meshes with sufficiently fine resolution to retain the microstructural detail, as shown in Fig. 7. Within the Simpleware software the meshing algorithm 'þFE Free' was used. Whilst discretising the surfaces and volumes, this algorithm will decimate using larger elements where possible whilst retaining geometric detail within set error limits. Reducing the total number of elements decreases the computational expense. Through trial and error, a coarseness setting level of À10 and À25 was found to produce a desirable balance between mesh 'quality' and number of elements for the LFA and monoblock meshes, respectively.

Simulation
To resolve small-scale features, IBFEM meshes can contain millions of elements compared with tens of thousands usually produced by CAD-based models. This reflects the true topology of complex surfaces (such as foams) that cannot be represented using primitive or simple CAD geometries. When topological detail is required, larger mesh sizes result in more calculations. Commercial FEM software packages cannot solve these large and complex simulations on desktop PCs. Furthermore, commercial FEM software is not suited to supercomputers because they do not make efficient use of parallelisation. Therefore, the open source solver ParaFEM, revision 2084, 4 developed by the authors [35e37], was used because it has previously been shown to scale well on parallel computing architectures i.e. the time to solve almost halves as the number of computational cores doubles [38], [39]. Visualisation and analysis of results was performed using ParaView, version 4.4 64-bit (Kitware Inc., Clifton Park, New York, USA).

Verification of thermal simulation
Verification of the parallel program developed by the authors was performed by constructing a problem whose results could be compared directly with an analytical solution. The problem chosen was a three-dimensional cuboid (confined by the domain -a < x < a, -b < y < b, -c < z < c) with unit initial temperature and zero surface temperature, as shown in Fig. 8. See Carslaw & Jaeger [40] for more details.
The analytical solution for the temperature, v, can be described by the triple Fourier series; where a l;m;n ¼ kp 2 , and the thermal diffusivity, k; is defined by its relation to thermal conductivity, density, and specific heat capacity, k ¼ K=rc p .
To represent the analytical problem in a FE form, a cuboid with dimensions a ¼ b ¼ c ¼ 2 was constructed. The main expected use of this program is image-based modelling which typically requires an unstructured mesh with large number of elements, therefore to best represent this use tetrahedral elements were chosen. The boundary conditions were such that the whole domain was given unit initial temperature and surfaces were fixed to zero for all subsequent time steps, with values being specified at the nodal locations. Material properties for oxygen free high conductivity copper (OFHC) were given. v ¼ jðx; aÞjðy; bÞjðz; cÞ v ¼ 64 Additionally, the same analytical problem was prepared in a well-known commercial FEA software package for comparison.

In silico experiment: laser flash analysis
The required boundary conditions to accurately simulate the laser profile of an LFA experiment have already been discussed in detail by Evans et al. [41], this work uses the same methodology. In brief, a thermal flux pulse is applied to one surface of the disc using a flux load distribution profile as shown in Fig. 5 to simulate the laser found in the experimental apparatus. The temperature rise on the opposite surface is tracked with respect to time to measure the half-rise time, t 1/2 . Using the exact same post-processing analysis as the laboratory experimental method, the half-rise time of the curve is used to calculate the thermal diffusivity, a, of a sample of thickness d, see Eqn. (3). This equation is a simplified version of a more complete mathematical model [42] which accounts for some real world mechanisms not included in these simulations (e.g. emissivity). When combined with density, r, and specific heat capacity, c p , this may be used in turn to calculate thermal conductivity, k, as shown in Eqn. (4). To further gain confidence in the model, the temperature rise, dT, of the sample may be verified against a calculated value by knowing the applied energy from the laser, Q, and the theoretical mass, m, from its density and volume, V, see Eqn. (5).
A time dependent heat flow analysis was carried out simulating a 0.252 s duration of the LFA experiment. For increased temporal resolution during periods of high gradients (i.e. near the start of the simulation) a variable time step size was used. Details of how the time step size was varied can be found in Table 3. The backward Euler time-stepping method was used with an iterative solver tolerance of 1.0 Â 10 À6 .

Case study: fusion energy heat exchanger component
When modelling the monoblock, boundary conditions were prescribed in such a way to emulate conditions in the divertor region of a fusion device, as shown in Fig. 6. These were a global initial temperature of 150 C, thermal flux of 10 MWm À2 on the plasma facing surface and a fixed temperature of 150 C on the inner pipe wall due to the water coolant.
A transient heat flow analysis was carried out to determine the time required to achieve a steady state temperature. The same backward Euler time-stepping method was used as in the LFA simulation. For increased temporal resolution during periods of high gradients (i.e. near the start of the simulation) a variable time step size was also used. Details of how the time step size was varied for this simulation can be found in Table 4.
The LFA simulation included a relatively short transient event, i.e. the laser pulse, and the purpose of that simulation was to accurately measure the travel time the pulse through the sample. Because the monoblock simulation only included steady-state boundary conditions, i.e. the plasma thermal load, it was possible to use fewer time steps with larger intervals to calculate the temperature profile.

Material properties
CAD based modelling typically averages localised variations, e.g. caused by microstructure, over a larger volume in a process called homogenisation, e.g. Refs. [43], [44]. In a CAD based model, the graphite foam would be represented by a solid volume representing both graphite and porous phases and be assigned effective material properties. Rather than homogenising, IBFEM modelling represents each distinct phase separately. The investigation presented in this paper uses both CAD and IBFEM simulations for comparison. To undertake the range of simulations described in section 3.4 several sets of material properties are required as described below.
For the IBFEM LFA simulations, the material properties of graphite were required which were obtained from literature values [45]. In addition to the anisotropy of the foam microstructure, graphite is a highly anisotropic material on the molecular scale. Its thermal conductivities near room temperature are 1950 W/mK and 5.7 W/mK parallel and perpendicular to the basal planes, respectively. To investigate the effect of this inherent property two variations of the IBFEM LFA model were used. Firstly, the models were given the literature values for graphite. One plane of high thermal conductivity was aligned in the direction the foam was manufactured (y axis). Then, assuming an equal distribution in the other two planes, an average of the high and low thermal conductivity values were given in the (x and z axes). Secondly, an isotropic thermal conductivity was applied in order to have a result that was purely affected by the foam microstructures and not anisotropic properties of graphite which used the same averaged value in all planes. This arrangement of material property alignment is shown in Table 5. The temperature dependent material properties of graphite were used, the values shown in Table 5 are near room temperature for comparison purposes. In addition to the IBFEM LFA models, a homogenised CAD model of a solid disc was solved for comparison. Two CAD models were created a) C foam using the graphite foam manufacturer's homogenised anisotropic material properties (i.e. effective properties measured experimentally, see Table 1) and b) using C avg_graphite i.e. where k ¼ 978 W/mK in all directions. The CAD simulations were performed to demonstrate whether simulating the LFA experiment in this manner is an appropriate method to back calculate the material's thermal conductivity. That is, to test whether results for thermal conductivity derived from simulating the experiment match the input values.
For the IBFEM monoblock model the C graphite values from literature were used. In order to compare with standard CAD practice, the same simulation was repeated twice using homogenised CAD based models. Firstly, the CAD anisotropic model used the C foam properties with the plane of higher thermal conductivity aligned parallel to the direction of the thermal load. Secondly the CAD isotropic model used the C avg_graphite material properties applied to the homogenised CAD interlayer. This was to create a comparative isotropic baseline to test the sensitivity of the model to changes in thermal conductivity of the interlayer. The values for all the material properties used are shown in Table 6. The temperature dependent material properties were used, the values shown are near room temperature for comparison purposes. For consistency and to remove any potential influence of mesh dependency all models used the same mesh geometry i.e. nodal coordinates and element structure. To homogenise the interlayer, the CAD models applied the same material properties to both foam and porous phases. The three variations of material property assignment, i.e. CAD anisotropic , CAD isotropic and IBFEM, are as labelled in Fig. 9.

X-ray computerised tomography (CT) scanning
Considering the distances between the cone beam X-ray source, sample and detector a voxel width of 35.4 mm was achieved. An annotated example tomography slice is shown in Fig. 10.

Image post-processing
An example tomography slice showing the steps from greyscale to binarised to skeletonised image is shown from the region of interest in Fig. 11. The skeletonised form is more easily viewed in 3D, which is shown in Fig. 12.
The sample's external dimensions and mass, m, were measured with a micrometre and a digital balance, respectively. The external volume, V, and effective density, r, were calculated from these values (i.e. V ¼ h x w x d and r ¼ m/V). Despite the inclusion of the 3.2 mm layer at the bottom, the calculated value of 671 kg/m 3 is relatively near the manufacturer's value of 675 kg/m 3 . These results are shown in Table 7.
The actual volume of a region of interest from the graphite foam, rather than external sample volume, was measured from the segmented image. This was found to be 2.72 Â 10 À5 m 3 . By also knowing its external dimensions it was possible to calculate the sample's porosity and the density of the graphite. These are also shown in Table 7 along with the manufacturer's value for porosity and a literature value for the density of graphite [45].
Shown in Fig. 13 is the variation in tortuosity with respect to  changing orientations in q and f, as defined in Fig. 12, as measured via the skeletonisation of a region of interest. For comparison with the orientations used by the manufacturer for measuring thermal conductivity, the tortuosity along the x, y and z axes are given in Table 8.
For the hybrid CAD-IBFEM monoblock model an additional variable is introduced in the form of converting CAD drawings of the armour and inner cooling pipe into the discretised image space i.e. curved volumes into voxels. It is useful to validate this process by comparing the calculated volumes from the prescribed dimensions with the volumes measured by counting voxels assigned to each part. This is shown in Table 9.

In silico experiment: laser flash analysis
The accuracy of the meshing to represent the graphite foam was first verified by calculating the porosity of the 'cut' LFA samples. This was accomplished by measuring the FE mesh volumes for the graphite foam samples from nodal coordinates and comparing these against the porosity of the segmented image in the locality of the surrounding region of interest. This is an informative check because the original image data is discretised into cuboids (voxels) and during meshing smoothing operations are performed to better represent the curved surfaces of the real material. These results are shown in Table 10. Table 7 Metrology of graphite foam sample as stated by the manufacturer, measured externally and by analysis of X-ray tomography image.    Table 9 for further dimensional changes from the original data introduced at the meshing stage, this is shown in Table 11. Additionally, Table 11 contains further analysis of the parts' geometries showing measured surface areas and the calculated ratio between volume and area.

Verification of thermal simulation
Fig. 14 compares the results between ParaFEM and the commercial solver using the same mesh and input parameters against the analytical solution. Even for such an 'idealised' simulation, it can be seen that the commercial solver includes a non-negligible error in its result which is greater that observed in ParaFEM's results. Therefore, it can be concluded that ParaFEM can offer sufficiently accurate results acceptable for this kind of simulation.
Firstly, it is worth considering the CAD models to validate the methodology of calculating thermal conductivity by simulating the LFA experiment. Fig. 15 a) shows the temperature rise on the rear face of the sample as would be measured by an infra-red camera in an experimental setup. To facilitate visual comparison of the halfrise time it is convenient to normalise the temperature variation between the sample's initial and maximum temperatures as shown in Fig. 15 b).
By comparing the thermal conductivity values calculated from the simulation results (see Table 12) with the original values that were input as the models' material properties (see Table 6), there is good agreement for all three samples (i.e. regardless of anisotropy in material properties). The difference between these values for the models where thermal conductivity did not vary with temperature, i.e. C foam┴ and C foam═ , was less than 1%. For the CAD model that Table 9 Dimensions and expected volume of each monoblock section according to design. OD ¼ outer diameter, ID ¼ inner diameter, l ¼ length.   included temperature dependant material properties, C avg_graphite, the difference between the room temperature thermal conductivity and that measured by the simulation is 8%. The initial and final temperatures of the sample during the test are 27 and 32.5 C. Therefore, not only does the thermal conductivity change during the simulation but also density and specific heat capacity which are used to calculate effective thermal conductivity. Because of this, the prescribed and measured thermal conductivities are not directly comparable therefore a difference of 8% is reasonable considering this. A useful additional validation is to compare the samples' final temperatures against analytical calculation. These match very closely for C foam┴ and C foam═ , and within 1% for C avg_graphite . These checks validate that it is possible to calculate a disc shaped sample's thermal conductivity by simulating the LFA experiment and processing simulation results in the same way as their experimental counterparts.

In silico experiment: laser flash analysis
For the IBFEM models there was a negligible difference in the results of the samples taken from the same plane (e.g. Xlow, Xcentre and Xhigh). For clarity in the data, only the results from the 'centre' samples of each plane are shown (i.e. Xcentre, Ycentre and Zcentre). Fig. 16 shows the normalised temperature rise on the rear surface of the samples using the material properties for a) C graphite and b) C avg_graphite . The thermal characterisation values derived from these curves are shown in Table 12. The error was calculated as the difference between the numerical result and the expected value.

Case study: fusion energy heat exchanger component
The following results demonstrate the use of IBFEM with a component simulated under in-service conditions. Fig. 17 shows a comparison of the temperature from a central plane of the monoblock once steady state had been reached for the three simulations performed (CAD anisotropic , CAD isotropic and IBFEM). The second row of images is focussed closer on the material interface region with temperature bounds rescaled to show additional information. The third row shows the surface temperature of the coolant pipe. Despite primarily providing qualitative data, it is helpful to visualise the temperature contours resultant from the thermal flux applied to the plasma facing surface and heat sink in the coolant pipe.
In order to obtain quantitative data from these models a temperature profile was taken through the monoblock. The path was drawn from the centre of the bottom surface to the centre of the plasma facing surface, as shown in Fig. 6. The results for each model are combined in Fig. 18, which also notes the change in material along the path (shown along the lower horizontal axis). Fig. 19 shows the temperature for the point at the start of this path (i.e. on the rear surface) with respect to time.

X-ray computerised tomography (CT) scanning
The manufacturer states that the average wall thickness in the graphite foam is 348 mm. Because the X-ray tomography setup used here yields a voxel width of 35.4 mm it is appropriately suited to accurately image the foam. The image provides an average of 10 voxels across the wall thickness.

Image post-processing
It can be seen that the values for porosity and density calculated from the processed image are comparable to within~2% of those   cited elsewhere (also shown in Table 7). This indicates that the digitisation of the foam on a microstructural level is an accurate representation.
When considering the variation of tortuosity with respect to changing orientation clear trends can be seen. As the azimuthal angle, q, varies from 0 to 90 , tortuosity decreases from~1.25 to 1.1 at the midpoint of 45 , then increases back to its original peak similar to 2p of a cosine curve. For variation in the polar angle, f, from À90 to þ90 a similar behaviour is observed from peaks of 1.22 to troughs of 1.18, but with three peaks emulating 4p of a cosine curve. It should be noted that the peak at 0 has a higher value of 1.26. The averaged values along the Cartesian axes show that the tortuosity in the x and z axes are similar to each other to within 1.5%. However, the average tortuosity is~6.6% higher along   the y-axis. This would indicate that the thermal conductivity would be lower along the y-axis, i.e. in the direction the foam was 'grown' during manufacture, compared to the perpendicular plane. This observation contradicts the manufacturer's measured values, Table 1. Therefore, it can be deduced that multiple mechanisms are contributing to the parent material's effective thermal conductivity. Although tortuosity is one of the factors it is an oversimplification of heat transfer, omitting mechanisms such as heat flow rates or heat losses. Image-based modelling will account for anisotropy in the local material properties, the cross-sectional areas of the flow channels through accurate topological representation and therefore better approximate the effective change.
For the image created from the CAD-IBFEM hybrid monoblock, the values for the calculated image volumes and prescribed CAD dimensions are comparable. The volumes of the inner CuCrZr pipe are identical however there is a small difference in both the armour and interlayer volumes. This is due to operations performed at the armour-interlayer interface to convert curved geometries into voxel data (cuboids) which have caused an overall change in outer diameter of the interlayer. The voxels for this image are cubes with a side length of 35.4 mm, therefore incorrectly determining the interlayer or armour dimensions by one voxel would result in a volumetric change of 3.1 mm 3 or 23 mm 3 , respectively. From the variation in expected and resultant volumes it can be calculated that the interlayer outer diameter is 13.87 mm rather than 14 mm. For this work, this level of divergence from design is not of importance as the main aim is to relatively compare the results for the CAD and IBFEM simulations (which use the same mesh and therefore have the same dimensions).

In silico experiment: laser flash analysis
One metric to indicate consistency in volume is to compare the porosity of the segmented 3D image data with that of the FE mesh. It can be seen that there is very little change with the LFA samples; the FE mesh porosity values are 4.6% lower than those for the whole block of graphite foam, shown in Table 7, which is likely due to the whole block containing the 3.2 mm layer at the bottom with a different microstructure where pores can be observed to be significantly larger, see Fig. 10.

Case study: fusion energy heat exchanger component
For the monoblock, the agreement between volumes for the segmented image and FE mesh is very good, any differences are negligible. As can be expected, the area to volume ratio is lowest for the armour and highest for the foam interlayer due to its fine microstructure. These ratios demonstrate that microstructural detail in the graphite foam have been retained during meshing.
Measuring the mesh quality shows that the meshes do contain some elements outside the desired quality thresholds for the in-out and edge length aspect ratios. However, these are only a very small percentage of the overall mesh (~0.0008%) and all elements are within acceptable limits for the other metrics. Given the fine resolution of the meshes these are unlikely to have a significant effect on the global results other than in extremely localised regions. Considering the mean and worst values for mesh quality it can be observed that the meshes are of adequate quality for FE purposes.

In silico experiment: laser flash analysis
Having verified the methodology for CAD simulations in section 4.4.1 it is now appropriate to evaluate the results of the IBFEM simulations. The fact that there was a high level of repeatability for samples digitally 'cut' from the same plane shows that there was a suitably large representative volume to exhibit the effective bulk behaviour in that plane.
It is worth first considering the results of the C avg_graphite simulations. The input thermal conductivity for these simulations is isotropic, therefore the results are not indicative of the real effective performance but rather are used to investigate effects of the anisotropy due to microstructure only. Fig. 16 shows the temperature rise on the rear surface of the sample from which the half rise time is measured. This is used to calculate thermal diffusivity and thus effective thermal conductivity, the values of which are shown in Table 12. Comparing the results of each of the planes it can be seen that the temperature rises a little slower for the Ycentre ( ┴ direction) sample compared with Xcentre and Zcentre samples (═ plane). This translates to an effective thermal conductivity that is 7.4% lower which indicates, as expected from the preliminary tortuosity analysis (where tortuosity was~6.6% higher in this direction), that the anisotropy in microstructure does impact the effective thermal conductivity and greater tortuosity has a corresponding greater reduction in effective thermal conductivity compared with the input values. As previously stated, this reduction of thermal conductivity along the y-axis contradicts results measured experimentally which show a higher effective thermal conductivity in this direction. Therefore, it is apparent that the effective thermal conductivity is dominated by the alignment of the graphitic planes rather than microstructure, but that both mechanisms have a non-negligible effect.
For the C graphite simulations, the concept was to use literature values for graphite combined with the microstructure to predict effective thermal conductivity. As seen in Fig. 16 the temperature rises more rapidly for the Ycentre sample ( ┴ direction) compared with the Xcentre and Zcentre samples (═ plane). This is opposite to the C avg_graphite simulations where the temperature rise was less rapid for the Ycentre sample but agrees with the experimental results.
For the Ycentre sample the thermal conductivity reduces in the ┴ orientation from an input of 1950 W/mK to an effective 243 W/ mK. The manufacturer's stated effective thermal conductivity is 240 W/mK which is a difference of 1%. For the Xcentre and Z centre samples the thermal conductivity reduces in the ═ orientation from an input of 978 W/mK to an effective 176 and 178 W/mK, respectively. The manufacturer's stated effective thermal conductivity is 64 W/mK which is a difference of 175%. For the ┴ direction the simulation is in good agreement with experimental results but poor agreement for the ═ plane. This suggests that the simplistic assumption made to impose the material properties of graphite on a cartesian coordinate system overlaid on the foam structure was insufficient to fully describe the anisotropic behaviour in all orientations.
Crystallographic work by Klett et al. [46] has shown that the graphitizing stage of manufacture causes a strong alignment of the basal plane along the ligament regions of the foam. The consequence of this is that there are fewer thermal barriers to cause phonon scattering along the ligament which leads to anisotropic thermal conductivity with the high values in that same direction. Observations show that there is less alignment in the junction regions therefore they exhibit a more isotropic behaviour. Fig. 20 shows a schematic describing this phenomenon. Characterisation work of graphite foams has shown a preferential alignment of the ligaments in the direction of the foam growth [46] [47], which would account for the higher effective thermal conductivity despite there also being higher tortuosity in this direction.
That is, in addition to having anisotropic microstructure, the crystallographic structure of graphitic foams is highly aligned in some regions but more amorphous in others. To better understand the complex relationship between the graphite foam's material properties and effective behaviour the improvement on the law of mixtures suggested by Markworth [48] was used to profile the change in thermal conductivity in the ═ plane, k ═ , with respect to the fractional change between a purely crystalline to purely amorphous material, f, as shown in Equation (6).
Further IBFEM simulations were performed using thermal conductivities derived from Equation (6) taken at various mix fractions. Fig. 21 shows the ratio of applied and effective thermal conductivities for each mix fraction. Using a second order fit to these results it was predicted that f ¼ 0.238 would produce the model with effective thermal conductivities nearest those stated by the manufacturer. This was tested and was found to give effective thermal conductivities of 237 W/mK in the y-axis, then 63.7 and 64.5 W/mK in the x and z axes, respectively. In this instance, this method has been used here to back calculate the relation between the anisotropy on the crystalline scale. The drawback of this approach is that it has been shown that the thermal conductivity in the ligaments and junctions is related to the graphitization rate. Therefore, this factor is specific to this grade of KFOAM and would need to be recalculated for other variants.
Experimentally, outliers are often found if a sample batch is sufficiently large. The IBFEM approach will accurately capture outlier behaviour if the deviation is caused by larger-scale foam geometry (e.g. presence of larger pores or completely fractured ligaments). IBFEM will not capture outlier behaviour if it is caused on a scale smaller than the X-ray image resolution (e.g. differences in meso-scale cracking within ligaments). Testing here was repeated by taking samples from various locations within the foam cube. There was a negligible difference in porosity of the samples (see Table 10) and consequently in the simulation results i.e. no outliers were found in this batch to demonstrate this claim.
Ideally, a better theoretical understanding of the changes in the crystallographic structure due to manufacturing processes would enable the mix fraction, f, to be identified analytically for a complete in silico characterisation of the effective thermal conductivity. An alternative approach would be to sub-divide the foam FE mesh into junction and ligament regions. In this case it would be possible to assign isotropic and anisotropic material properties to the junctions and ligaments, respectively. The material property values would be derived from the rate of graphitization with high thermal conductivity aligned along the length of the ligament. This is not trivial because the conductivity alignment would need to be calculated for each ligament by using a transformation matrix to apply it within the cartesian coordinate system. If achieved, it should provide a much-improved result.
As a final note on the image-based modelling approach to thermal characterisation, it is also possible to invert the analysis procedure to 'design' a foam with desirable material properties based on a true organic structure. That is, once an example of a manufactured foam had been digitised via X-ray CT it is possible to digitally alter the image to investigate impact on behaviour. For example, focussing on features controllable during manufacturing, the density or porosity fraction can be changed by eroding or dilating the segmented image to emulate having thinner or thicker ligaments for an optimal performance.

Case study: fusion energy heat exchanger component
In addition to investigating the characterisation of novel materials with in silico experiments it is useful to test their in-service performance within a component to identify the best candidates prior to manufacturing. The case study presented here is included as a demonstration of how the image-based modelling approach may be used instead of CAD-based modelling to investigate this without having measured the material's effective material properties experimentally to be included as input data and, despite this, provide increased accuracy through the inclusion of microstructural detail.
From images (a-i) to (c-i) in Fig. 17 it can be seen qualitatively that (a) CAD anisotropic and (c) IBFEM have similar temperature profile distributions through the cross-section whereas (b) CAD isotropic differs by generally having lower temperatures. Results for the CAD isotropic model are not included to accurately represent the foam interlayer, but rather because they are a useful baseline for comparison of other results. Therefore, CAD isotropic results will only be referred to when a noteworthy discussion point arises.
Quantitatively, one of the features of interest is the thermal gradient across the sample which is of significance because higher thermal gradients indicate higher thermally induced stresses due to non-uniform thermal expansion. An example thermal gradient can be taken by considering the distance from the peak temperature, located on the top corner, down the edge perpendicular to the plasma facing surface to the location at which the temperature reaches 220 C. For the (a) CAD anisotropic and (c) IBFEM models the peak temperatures are 916 C and 972 C giving thermal gradients of 49,360 and 51,860 C/m respectively. a difference of~5%.
An obvious difference between the CAD anisotropic and IBFEM models is the temperature field near the interlayer region. Because CAD anisotropic uses a homogenised interlayer the temperature field  is perfectly symmetrical about the centre, see Fig. 17 (a-ii). Although the temperature for the IBFEM model is similar in general, localised fluctuations can be observed see Fig. 17 (c-ii), due to specific locations where the foam is in contact with the armour or coolant pipe. This leads to a small amount of asymmetry overall, that is, the gradient is different on the left-hand side of the armour compared to the right-hand side.
The fluctuations and non-symmetry observed in the IBFEM model (not seen in either CAD model) have been introduced by the non-uniform geometries included when modelling the graphite foam on the microstructural level. The impact on the sample is demonstrated by localised 'hot-spots' on the surface of the coolant pipe as shown in Fig. 17 (c-iii). Although this only causes small global variations, the ability to consider each manufactured part individually may be invaluable for optimising the assembly of a multi-part component. This is a clear benefit over conventional CAD based modelling.
The main interest in Figs. 18 and 19 is to compare the CAD anisotropic and IBFEM models to determine whether the new approach of simulating without knowing effective material properties priori gives comparable results to the conventional FEA approach. The CAD anisotropic and IBFEM models show a relative level of agreement despite using a fundamentally different approach to represent material properties for the interlayer. That is, the CAD anisotropic values are for bulk foam properties as measured experimentally and the IBFEM values are theoretical based on combining literature values for carbon with the foam microstructure. The difference in the applied thermal conductivities in the direction of the thermal load is 88% but there's only a difference of 7% in the peak temperature in Fig. 18. As could be expected, the main difference is the temperature profile through the interlayer region where the models completely differ in their approach. The profile for the IBFEM model fluctuates significantly due to passing through foam and porous phases whereas it is continuous for the homogenised CAD anisotropic model. These are extremely localised fluctuations built into the model, and therefore, when considering the global effect of the interlayers the results are comparable. The second observable difference between the two models is a small variation in temperature on the bottom surface of the monoblock, 167 C and 174 C for the CAD anisotropic and IBFEM models respectively.
The level of agreement shown in this case study demonstrates the potential benefits and validity of the IBFEM approach with respect to the widely accepted method of homogenising complex microstructures for CAD-based modelling.
By using literature values for material properties and microstructure alone IBFEM may predict the thermal performance of components containing graphite foam without prior knowledge of its effective bulk performance and without measuring those experimentally. If foam with a new microstructure was produced it would be possible to model it from X-ray CT data without the need to perform a time-consuming range of thermal testing. This could be a powerful tool for rapid development of new functional materials. The CAD anisotropic model gives a geometrically ideal, and therefore symmetrical, result which may be globally valid but does not provide localised information within the foam. The IBFEM model may be interrogated at smaller scales to investigate local fluctuations resulting from the microstructure which could potentially be significant to the part's structural integrity. Developments underway in this wider field of materials simulation will allow modelling of the structural integrity at an even finer scale by using a multiscale cellular automata e finite element (CAFE) approach [49].
Finally, the simulation results show significant agreement; differences observed are primarily due to the increased microscale accuracy inherent in the IBFEM model. In scenarios such as this, where microstructure is a significant contributor to performance, additional confidence in predicting performance from a 'digital twin' model of a real part may be of significant operational value if the design is approaching safety limits.
It should be noted that the authors do not propose that imagebased modelling fully replace design-based modelling or lab based experimental testing but are considered a complementary tool to expedite development. Design-based modelling comes earlier in the cycle, before the manufacturing stage, to identify candidate designs based on theory and previous knowledge.
Here, a hybrid CAD-IBFEM model was presented, which assumed an ideal interface between the image-based region (foam) and the design-based regions (armour and coolant pipe). This could be used to investigate the manufacturing process e.g. selecting the optimal region from a foam block to machine for use in a component and how it should be aligned. Alternatively, simulating the bonding stages could assess impact on the ligaments at contact regions (e.g. compression leading to fracture) or optimising braze filler material by simulating material flow. Yet again, the hybrid approach could be used to investigate designing the interface area and topology to optimise for such things as heat transfer.
Full image-based modelling of a complete part can be used to reduce the number of samples that require manufacturing because the same sample can be digitally tested to destruction an infinite number of times. Because of the additional microstructural detail, image-based modelling has been shown to increase results accuracy compared with design-based modelling. However, the results will only be as accurate as the theory for the range of physical, chemical and crystallographic mechanisms included within the simulation and will not include unknown regimes. For example, the methodology demonstrated here relies on continuum mechanics and does not consider atomistic level interactions, like at material interfaces. Therefore, image-based modelling can be used to reduce testing but relies on it for verification whilst not all aspects of the material are included within the analysis.

Conclusions
Functional carbon-based materials with material properties which can be 'tuned' exist that are of interest for high-value manufacturing e.g. in nuclear applications. Some of these materials have complex anisotropic micro-structures which influence the effective (bulk) material properties. It can be simple to modify these micro-structures during manufacturing to change the material properties, but experimentally characterising the new effective performance is often time-consuming and expensive. This work presented a technique to perform characterisation experiments in silico by digitising novel materials via X-ray tomography to enable virtual testing with the image-based finite element method (IBFEM). By accounting for micro-structure to predict effective performance, the simulation only requires knowledge of the parent material which only needs characterising experimentally once. This may therefore be used to expedite the development of functional materials by assisting downselection of candidate materials. Presented here was an example that investigated the thermal conductivity of a graphite foam by simulating laser flash analysis.
In this particular example it was observed that the effective thermal conductivity was caused by a combination of the graphite foam's microstructure and crystallographic structure. There was insufficient data about the crystallographic structure. Consequently, the assumptions made about input material properties were too simplistic to predict the effective material properties in all orientations with a high degree of accuracy. Using literature values for the material property of graphite, the foam's effective thermal conductivity was predicted to be 243 and 177 W/mK in the ┴ and ═ plane, respectively. The difference between input and effective thermal conductivity was due to thermal barriers caused by the extended thermal path between one surface of the sample to the other. The graphite foam manufacturer's values, measured experimentally, were 240 and 64 W/mK. By using the law of mixtures, it was possible to adjust the input thermal conductivity such that the effective values matched those of the manufacturer. Presently, the required mix fraction was found empirically meaning increased understanding of the crystallographic structure is needed. Therefore, this study demonstrated that this use of IBFEM as a tool to predict effective material properties shows promise but that its effectiveness depends on the level of detailed knowledge about the parent material. Its use with less complex materials such as metals may be more appropriate.
This work then demonstrated how the image-based modelling technique may be used to 'digitally manufacture' a component containing such a functional material to assess its performance under in-service conditions. Doing this in silico reduces the dependency on 'real' manufacturing and potentially logistically challenging experimental testing. A case study was presented where the functional material (graphite foam) was used as part of a heat exchange component, termed monoblock, for a fusion energy device. The other parts of the component were drawn using CAD, thus a hybrid CAD-IBFEM model was created.
To test the validity of the CAD-IBFEM model its results were compared to a standard CAD models of the same component where the microstructures in the functional material region were homogenised and average material properties assigned. The CAD model used the manufacturer's anisotropic material properties for the graphite foam which were measured experimentally.
The IBFEM model was in good agreement with the CAD version which used the manufacturer's properties. This demonstrated that despite having no prior knowledge of the bulk performance of the graphite foam, faithfully modelling the geometry of the material on the microstructural level yielded results comparable to homogenisation techniques. The advantage of the IBFEM technique is that, in addition to having been shown to increase modelling accuracy for such materials, the model can be interrogated on a more localised level to provide potentially critical additional information. This capability could be used to investigate optimal interfaces between complex functional materials and their more conventional counterparts.
Although IBFEM has been used in the context of fusion engineering in this work, it could be used in a broad range of applications, particularly where materials with complex behaviours are used and bulk performance may not be easily measured. Ultimately, it is envisaged that IBFEM could be developed for quality assurance on production lines to gain confidence in component integrity, thus either increasing fabrication yields or reducing engineering reserve factors.