Multiphysics analysis with CAD-based parametric breeding blanket creation for rapid design iteration

Breeding blankets are designed to ensure tritium self-sufficiency in deuterium–tritium fusion power plants. In addition to this, breeder blankets play a vital role in shielding key components of the reactor, and provide the main source of heat which will ultimately be used to generate electricity. Blanket design is critical to the success of fusion reactors and integral to the design process. Neutronic simulations of breeder blankets are regularly performed to ascertain the performance of a particular design. An iterative process of design improvements and parametric studies are required to optimize the design and meet performance targets. Within the EU DEMO program the breeding blanket design cycle is repeated for each new baseline design. One of the key steps is to create three-dimensional models suitable primarily for use in neutronics, but could be used in other computer-aided design (CAD)-based physics and engineering analyses. This article presents a novel blanket design tool which automates the process of producing heterogeneous 3D CAD-based geometries of the helium-cooled pebble bed, water-cooled lithium lead, helium-cooled lithium lead and dual-coolant lithium lead blanket types. The paper shows a method of integrating neutronics, thermal analysis and mechanical analysis with parametric CAD to facilitate the design process. The blanket design tool described in this paper provides parametric geometry for use in neutronics and engineering simulations. This paper explains the methodology of the design tool and demonstrates use of the design tool by generating all four EU blanket designs using the EU DEMO baseline. Neutronics and heat transfer simulations using the models have been carried out. The approach described has the potential to considerably speed up the design cycle and greatly facilitate the integration of multiphysics studies.


Introduction
Breeding blankets are designed to fulfil several high-level plant requirements, including tritium self-sufficiency, shielding non-sacrificial components from the intense neutron flux and producing heat which is ultimately used to generate electricity. Designing and engineering components for use within fusion reactors is challenging due to the high radiation fluxes and significant heat loads that they experience. Maintaining an operational and safe component within the inner vessel of a fusion reactor presents a range of difficulties; however, adding functional requirements such as tritium breeding, heat generation and heat removal further complicates the task.
Methods of design optimization such as parameter studies and a 'designing by analysis' approach are possible avenues for designing fusion reactor components that could provide solutions to this challenge. Such methods rely on human intuition and iterative analysis of models to close in on an optimal solution. Performing analysis in an isolated discipline will only find the optimal solution for performance metrics that are obtainable within that discipline. For instance neutronics optimizations may find the tritium breeding ratio (TBR) but may find unacceptable temperatures. Multiphysics analysis is required to optimize component design. To maintain data provenance it would be preferable to have a single model basis when sharing data between analysis techniques.
Traditionally, models are generated for neutronics using constructive solid geometry (CSG) and the models are suitable for use in parametric studies. Engineering analysis tends to require computer-aided design (CAD) models and CSG models are typically not compatible with engineering programs. Models for use in engineering analysis are often created via graphical user interfaces. The process of creating new engineering and neutronics models can be a time-consuming exercise. This is compounded since the models must be generated with the release of a new EU DEMO baseline design, and there are also four EU blanket designs for each iteration.
To analyze the performance of different designs within the parameter space it would be desirable to be able to produce an accurate three-dimensional (3D) geometry rapidly. Adopting a common geometry format would allow geometry to be used in multiple domains. Allowing fine details (such as cooling pipes) to be included or excluded during model generation can facilitate specific requirements of the particular analysis. Use of open source geometry-producing software such as FreeCAD [1] (used for this project), Salome [2] or PythonOCC [3] can be used to quickly generate parametric CAD geometry which can be exported into a variety of formats. CAD files in STEP format [4] are an open file standard compatible with engineering simulation software. STEP files can be easily converted into surface faceted geometry (e.g. h5m or STL files) for use in neutronics codes such as Serpent 2 [5] (used for this project) and DAG-MCNP5/6 [6]. Ideally any solution to making component models would be flexible enough to work with new DEMO baseline models and also to produce different blanket designs.

Geometry creation
It is clear that in order for any advanced method to generate any detailed geometry it must require the minimum human intervention, so it was determined that a software library should be created that allows arbitrary geometric operations to be performed, with the end goal of creating parametrically built blanket modules for the EU Demo program. This software is called the 'Breeder Blanket Model Maker' and can be found here [7]. Routines for the generation, modification and serialization of blanket envelopes were created that ultimately automatically produce detailed heterogeneous blankets for use in DEMO modelling. Demonstration neutronic and heat diffusion simulations were performed to illustrate the ease of carrying out parameter studies. Parametric models of all four EU blanket designs were generated to demonstrate the design tool. These include the helium-cooled pebble bed blanket (HCPB) [8][9][10], helium-cooled lithium lead blanket (HCLL) [11,12], water-cooled lithium lead blanket (WCLL) [13] and dual-coolant lithium lead blanket (DCLL) [14]. The process has been broken down into two parts, the first flowchart (see figure 1) summarizes the construction of all the non-breeder zone components, this includes the first wall, armour and rear plates. The second flowchart (see figure 2) summarizes the construction of the breeder zone structure.
The breeding blanket is segmented into different modules which have different shapes, orientations and positions depending upon the positioning within the reactor (see stage 1 in figure 1). The blanket designs share some common features such as: filleted corners on the toroidal HCPB, HCLL and WCLL designs or poloidal DCLL edges.

Structural components.
The procedure for the creation of the structural part of the blanket module is performed first. A number of key parameters define the structure of the blanket: first wall armour thickness, first wall thickness, rear plate count and thickness, and end cap thickness. An automated procedure regarding the automatic construction of a full detailed blanket structure was defined, and the full implementation details can be found in [7]. Additional fine detail is also included at the user's discretion including the introduction of fillets and cooling channels. The overall programmatic flow is shown in figure 2, where differing blankets follow different logical routes.

Breeder zone components.
In order to generate the full heterogeneous blanket description the internal detail of the breeder zone must be generated. This is a two-stage process: first the cooling structure is generated starting from the breeder zone envelope and generating the cooling structure from it then the cooling structure is subtracted from the breeder zone envelope, leaving the non-structural breeding material (lithium lead/lithium ceramic/neutron multiplier).

2.1.3.
Cooling structure generation. The segmentation of the cooling structure varies for each of the four breeder blanket modules, but is almost entirely a combination of poloidal-, toroidal-and radial-based segmentations. The cooling structure of the HCLL advanced plus module [12] can be represented by a series of poloidal segmentations with alternating layers of 1. Selection of example blanket module envelope from the EU DEMO baseline [15] 2. Example blanket envelope showing the front face (green), poloidal edges (blue) and toroidal edges (red) 3. Example blanket envelope with filleted poloidal edges (right) and toroidal edges (left) the fillet radius has been increased to clearly show the operation 4. Example blanket envelope with first wall armour shown in green 5. Example blanket envelope with first wall shown in green 6. First wall cooling channels (optional) 6. First wall cooling channels (optional) 7. Example blanket envelope with end caps shown in green 8. Example breeder zone envelope with back wall components shown in blue and green stiffening plates. This can be reproduced using alternate poloidal segmentation with alternating poloidal extrusion lengths, as shown in figure 2 section 2. In the case of the HCPB module there are alternating poloidal layers of lithium ceramic and neutron multiplier between the stiffening plates. The poloidal segmentation functions have been designed to allow any number of layer repetitions; this is demonstrated in stage 3 of figure 2). The wedge-shaped regions at the upper and lower extremities of the HCPB module are filled with neutron multiplier and are therefore not considered part of the cooling structure. The software is able to identify these wedge-shaped regions and group them with the other neutron multiplier regions.
Radial cuts, and thereby radial segmentation, are also implemented; these are required for both the WCLL and the DCLL blanket designs. The WCLL cooling structure can be generated with a combination of poloidal and toroidal segmentation (see stage 6 in figure 2). Both the toroidal and poloidal directions have alternating thicknesses for the structural plates and the lithium lead regions. Every other layer of the poloidal structural plate has an offset from the first wall that allows lithium lead to flow between plates. The poloidal segmentation for such a model can be carried out in a similar way to the HCLL, but the WCLL has an additional complication which requires radial segmentation. The WCLL model requires toroidal segmentation and additionally requires that the upper and lower wedge volumes should be considered to be entirely lithium lead. The resulting product of the toroidal segmentation can be seen in stage 6 in figure 2.
The DCLL cooling structure can be formed from a combination of radial and toroidal segmentation plus some detail to guide the flow of lithium lead. The procedure used was to first radially segment the blankets into three or five parts (depending upon the radial depth of the blanket). In general most of the inboard blankets accommodate three radial layers and the outboard blankets accommodate five radial layers. The addition of toroidal segmentation to the previously radial segmented breeder zone forms the first stage of the DCLL model (see stage 8 of figure 2). The DCLL blanket design allows the lithium lead to flow around the structure. An additional structural component at the upper end of each blanket module is also required by the DCLL design, the only additional complication is that a Boolean subtraction with the first radial layer is also required to obtain the desired structural plate shape (see stage 10 of figure 2).

Breeding material.
The complete description of the breeder zone comprises the description of the cooling structure and the description of the breeding material. The final stage is take the original breeder zone envelope and subtract the newly created cooling plates/stiffening plates from it, thus defining the description of the complete breeding zone. HCLL, WCLL and DCLL all require that the end regions are lithium lead and HCPB requires the end regions to be neutron multiplier.
2.1.5. Slice geometry generation. The slice geometry used in the thermal simulation can also be generated automatically. The procedure for generating the slice geometry of the HCLL is outlined in figure 3. The user specifies the blanket module from which to extract a slice. A slice envelope is created with a poloidal height equal to the poloidal height of a stiffening plate plus the poloidal height of the breeder zone. The envelope poloidal position is centred around the stiffening plate so that the slice contains half a breeder zone above and half a breeder zone bellow the stiffening plate. The cooling channel positions were found by offsetting the first wall towards the rear of the blanket, with the size of the offset progressively increasing. Surface identification for cooling surfaces was carried out by merging the coolant volumes with the structure volumes and searching for merged surfaces in the resulting geometry. The identification was manually checked after this stage to ensure a robust procedure.

Parametric geometries
As a result of the method previously described there is now an automated procedure for obtaining semi-detailed CAD geometry for the HCPB, HCLL, WCLL and DCLL blankets. The process relies on a library of common functions which can be mixed and matched to create particular blanket designs. The breeder blanket design tool is released as an open source project under the Apache 2.0 license and distributed via the UKAEA Github repository [7]. The software is subject to a test suite and the build status is updated automatically with every commit. Continuous integration practices are employed using Circle CI and Docker.
The model construction process is parametric, which allows models required for parameter studies to be generated rapidly. Currently the parameters that a user can input are: • filename of blanket envelope required for segmentation • blanket type (HCPB, HCLL, WCLL, DCLL) which also defines the geometry layout as shown in figures 1 and 2 • poloidal fillet radius for the first wall and first wall armour • toroidal fillet radius for the first wall and first wall armour • first wall armour thickness • first wall thickness • end cap thickness • thickness of each rear plate • thickness of each poloidal segmentation • thickness of each toroidal segmentation  • thickness of each radial segmentation • first wall coolant channel poloidal height • first wall coolant channel radial height • first wall coolant channel pitch • first wall coolant channel offset from the front face • output file format (STEP or STL) and tolerance.
Not all parameters are needed for each design as some are not applicable, for instance the breeder zone in the HCPB blanket has no radial segmentation option and does not require this input. Figure 4 shows each of the four blanket designs formed from a particular module from the baseline DEMO model [15]. Currently the tool requires that the first wall is a flat plane to determine the location of the internal structure.
The tool would therefore need some alterations to work with blanket envelopes with curved front surfaces (i.e. full bananashaped segments). The process of building a blanket module from an envelope typically takes a few seconds on a single core. Build time depends on the input parameters, as many very small layers would necessitate more Boolean operations than for the case of a few large layers. The process is parallelizable, and therefore a model such as the EU DEMO with 26 blanket modules typically takes less than 5 min on a quad core Intel i5 7600 CPU.

Neutronics model creation
Once the parametric CAD models have been created, one potential use is in neutronics simulations. There are several routes from CAD to neutronics models, such as conversion to CSG using conversion software such as McCad [16] or SuperMC [17]). Alternatively the use of faceted geometry is also possible. Previously, parameter studies for fusion blanket optimization have converted parametric CAD models to CSG   models using McCad and performed the simulation using MCNP [18]. This study opted to simulate using faceted geometry in the STL file format and perform the neutronics simulation using Serpent 2, which natively supports STL geometry. The process of converting from STEP to STL is quicker than from STEP to CSG and the results are easier to visually verify.
To demonstrate practical use of the parametric geometry, a series of tritium breeding simulations were obtained for the HCLL. The poloidal height of the lithium lead sections and the 6 Li enrichment of the lithium lead were varied independently. The thickness of the first wall and the stiffening plates were varied simultaneously with the poloidal height of the lithium lead sections using equations (1) and (2).
Here FWT is the firstwall thickness (m), CPT is the stiffening plate thickness (m), C D is the radial width of the first wall cooling channel (0.01 m), C p is the specific heat capacity (J K -1 ), P s is the coolant pressure (10 MPa for the helium coolant), SmD is the stress limit criterion for EUROfer (274 MPa) [19] in case of accident and LL p (m) is the poloidal height of the lithium lead sections. These equations are described in more detail in [11] and [20]. This HCLL study is a demonstration of the model-making tool developed and thermal-mechanical constraints are not taken into account. Halton sampling [21] was used as the sampling technique to select points within the parameter space. The parameter space encompassed blankets with a poloidal height of lithium lead between 0.01 m and 0.12 m and 6 Li enrichment between 0% and 100%. The requested simulation points found using the Halton sampling method are shown as red crosses in figure 6. These requested input parameters were uploaded to a online cloud-based database (Mongo Atlas [22]). Entries within the database were flagged with 'in progress', 'completed' or 'requested' to indicate their simulation status. The centralized accessible database allowed independent simulations to be carried out in parallel and the results to be coordinated. A containerized workflow was implemented using Docker [23] which held the breeder blanket model-maker software along with all the dependences required, such as neutron interaction data, FreeCAD, Python and MongoDB. Containers were then launched on EGI FedCloud resources [24] and connected to the cloud database to retrieve their simulation input parameters ( 6 Li enrichment and lithium lead poloidal height) and set the database entry for the simulation to 'in progress'. Once the simulation input parameters were received the CAD models were automatically built and combined with material The Gaussian process software used [30] was able to fit the TBR values along with their statistical errors and find the confidence values. The reference design HCLL has 90% 6 Li enrichment, 34.5 mm of poloidal lithium lead and achieves a TBR of 1.235.
properties to create the neutronics model. Upon completion of the neutronics simulation the results (TBR and heating tallies) were uploaded to the cloud-based database and the entry flag was set to 'completed'. At this point the container would either continue with the next requested simulation in the database or terminate and free up resources if all the requested simulations were complete. FENDL 3.1b [25] cross sections were used for neutron transport. Neutronics materials definitions from the EUROfusion material composition [26] were used in the model and a parametric plasma source based on [27] with plasma parameters from [28] was used. The number of starting particles run for each TBR simulation was 1 × 10 7 .

Neutronics simulations
The models generated are suitable for neutronics simulations and figure 5 shows the blanket models within a Serpent 2 geometry. Figure 6 shows the resulting TBR values from a neutronics parameter study. The TBR was found to change with the poloidal height of lithium lead. Models with a small lithium lead poloidal height contain a relatively large EUROfer fraction compared with models with a large lithium lead poloidal height. This is due to the large number of stiffening plates and it appears to have reduced tritium production. However, models with a large poloidal height can also have large quantities of EUROfer in the breeder blanket. As the poloidal height of lithium lead increases the thickness of the EUROfer first wall (see equation (1)) and the thickness of the EUROfer stiffening plates (see equation (2)) also increase.
There appears to be an optimal poloidal height which becomes more pronounced with 6 Li enrichment. Figure 5 shows the variation of TBR with 6 Li atom fraction within the lithium lead. Increasing 6 Li enrichment shows an increase in TBR as confirmed by previous studies [29]. Figure 5 shows the variation of TBR with the poloidal height of the lithium lead regions within the breeder zone of the blanket. The poloidal height has less of an effect on TBR but can be optimized. The achievable increase in TBR depends upon the 6 Li enrichment. Figure 6 shows that each different enrichment of 6 Li has a different optimal height for the poloidal lithium lead. At 90% 6 Li enrichment a TBR increase of nearly 0.1 is possible (see figure 6). The size of the 5σ confidence regions varies depending on the proximity and statistical error of nearby simulations. This is most noticeable towards the extremities of the search space where there are fewer simulations and the size of the 5σ confidence regions is larger. The variation in height of the poloidal lithium lead also has mechanical considerations as the first wall thickness and stiffening plate thickness are also increased when the lithium lead poloidal height is increased (see equations (1) and (2)). This helps explain why we observe an optimal lithium lead poloidal height for TBR values from the neutronics model.
A maximum TBR of 1.278 ± 0.010 (5σ confidence) was found using Gaussian process software [30] to fit the simulation data and statistical error (see figure 6). The highest TBR value was found for a blanket design with a lithium lead poloidal height of 0.061 m and a 6 Li enrichment of 100%. Additional constraints such as the capability of the thermalhydraulic design to cool the structure with reasonable pressure drops must also be considered, and the maximum TBR design may not meet such requirements. Thermal modes are developed in the next section.

Creation of the heat diffusion model
A slice of the HCLL blanket geometry was used to create a simplified model of the temperature field in the tungsten, EUROfer and lithium lead. Dimensions for the first wall thickness and stiffening plate thickness were calculated using equations (1) and (2) with the poloidal height of lithium lead set at 34.5 mm. The model contains a single stiffening plate with cooling channels and lithium lead either side, encased with a first wall; the material layout is shown in figure 8. A mesh with 2.3 million tetrahedra was created using Trelis [31] complete with boundary conditions and volumetric source terms. Heat diffusion simulations were carried out using FEniCS [32] which was able to apply boundary conditions to surfaces and temperature-dependent materials properties to different volumes. The heat diffusion equation used is described by equation (3).
where T is the temperature in K, λ is the thermal conductivity of the given material expressed in W m −1 K −1 and Q is the volumetric source term in W m −3 . As this equation is solved using the finite elements method, equation (3) needs to be brought to its weak formulation (or variational for mulation) as follows: where Ω is the domain on which equation (4) is solved, n is the normal direction on the external surface and v is a test function. The integration term on the boundary of Ω is determined by the boundary conditions.

The Robin boundary condition.
The Robin boundary condition allows the assignment of a convective heat flux on a boundary. In equations (5) and (6), it is shown that the heat flux depends on the temperature of the fluid and the convec- . This coefficient depends on the type of convection (natural, forced, laminar or turbulent) and the fluid in contact with the surface: Γ FWCC and Γ HSPCC are surface domains and are shown in blue and purple, respectively, in figure 7. This condition is used on the surfaces of the helium cooling channels. The coefficients h have been calculated for the first wall cooling channels (FWCC) and horizontal stiffening plate cooling channels (HSPCC). They are determined using Gnielinski correlation (equation (7)) with the parameters in table 1 in accordance with [33].
where Nu D is the Nusselt number, D h is the hydraulic diameter in m, λ He is the thermal conductivity of the fluid in W m −1 K −1 , Re D is the Reynolds number and Pr is the Prandtl number. We consider a smooth surface, and thereby the Darcy-Weisbach friction factor ξ is then given by: 3.3.2. Neumann boundary conditions. By using Neumann boundary conditions, fixed heat flux can be assigned to the front wall armour surface shown in red in figure 7, as described in equation (8): Here, J is a fixed flux density in W m −2 on a boundary. This condition is used on the surface Γ FW which corresponds to the front wall (shown in red in figure 7) with a flux of J = 0.5 MW m −2 which corresponds to the heat flux emitted by the plasma. The rest of the external surfaces are considered as insulated. This assumption is valid as long as these surfaces are part of a vacuum and are not exposed to an intense heat flux. The values of the thermal conductivity λ in W m −1 K −1 in equation (4) were found in [34].
Finally, the distribution of the volumetric source term Q in equations (3) and (4) was taken from [35] and is inputted into the finite element model using the following equations: where r is the radial distance from the front face of the breeder blanket in m. The resulting spatial distribution of volumetric heating (Q) in MW m −3 is shown in figure 9.

Heat diffusion simulations
Using the same CAD-generated models as in section 3.2, heat diffusion simulations have been performed. The steady state solution of the temperature field is shown in figure 10. Thanks to these simulations, we are able to determine the maximum temperature reached by each material and determine if the design allows the materials to stay within their maximum operating temperature limits (550 °C for EUROfer and 1300 °C for tungsten [36]). Although the design limit of 550 °C is reached in part of the stiffening plates, this design could be refined with additional cooling channels to reduce the temperature. Additional assessments involving computational fluid dynamics could be performed to check the accuracy of the temperatures predicted and the conservatism of assumptions made in this simulation. The maximum temperatures are shown in table 2.

Discussion
The simulation results show that the materials stay within their maximum operating temperature limits. However, one must be aware that some assumptions have been made that could introduce inaccuracy. The convective boundary condition on Γ HSPCC and Γ FWCC are set with constant average coolant temperatures along the channel. This assumption is reliable in term of flux balance in the module although it may introduce errors in the region near the inlets and the outlets of the channels where the coolant should be, respectively, colder and hotter than the average. One way to solve this issue would be to allow T coolant to vary along the channel. Actual cooling channel designs for HCLL might also differ from the one used in this paper, leading to variation in the temperature field. One can also notice that the simulation has been run in steady state.
Running transient simulation could allow new hot spots to be spotted before the module reaches steady state. The volumetric heat source used in this paper is a fitted exponential function mapped on a mesh using discontinuous Galerkin elements of order 0. In other words, the volumetric heat source has been discretized. Again, this might lead to inaccuracy in some big cells. This can be solved by refining the mesh, as currently the model contains 2.3 million tetrahedral elements. Finally, the lithium lead flow is not modelled in this paper for simplification purposes. Lithium lead has been considered as a solid. However, considering the very low velocity of lithium lead in HCLL modules, this assumption can be made.

Conclusion
A design tool capable of generating parametric designs for fusion breeder blankets has been demonstrated on singlemodule blanket envelopes for HCLL, HCPB, WCLL and DCLL. A wide range of design parameters can be changed to generate CAD geometry for use in parameter studies. The geometry generated is available in CAD format (STEP) and faceted geometry (STL and h5m). Conversion to CSG for neutronics simulation is achievable via existing software such as McCad or MCAM. The option of faceted geometry allows CSG geometry to be avoided in favour of more CAD-based neutronic simulation techniques such as DAGMC or Serpent 2. The provision of CAD geometry also enables manipulation to be performed with standard CAD software as opposed to CSG geometry where manipulation of the shapes is less convenient. A demonstration neutronics parameter study was performed, where poloidal lithium lead height and 6 Li enrichment were varied. This was done in order to optimize TBR for the HCLL module. A maximum TBR of 1.278 ± 0.010 (5σ confidence) was found using Gaussian processing to fit the data. The highest TBR value was found for a blanket design with a lithium lead poloidal height of 0.061 m and a 6 Li enrichment of 100%. Currently the tool allows for a large range of specific blanket module geometries to be made for use in simulations. Thermal-mechanical and other constraints are not taken into account when constructing the geometries and future research will be required to identify allowable design parameters that satisfy thermal-mechanical and thermal-hydraulic requirements. The geometry created can also be used in finite element and finite volume software to simulate heat diffusion, tritium diffusion and stress. A simplified heat diffusion problem was demonstrated in this paper. The maximum temperature within the different materials present in the midplane geometry slice of the HCLL outboard equatorial blanket module was found to be similar to previous research [33]. Maximum temper atures were 774 K for tungsten, 863 K for lithium lead and 824 K for EUROfer. The work performed here shows the value of in silico design processes, multiphysics workflows and, critically, integration with automated systems for the generation, submission and analysis of calculations. The workflow was built using modern scalable techniques to run the simulations (containerized cloud computing), store the output of the simulations (centralised cloud databases) and analyze the results (machine learning). The tool is open for extension and additional analysis such as thermal stress, mechanical stress and tritium diffusion should be the next steps for development.