Normal and Accidental Scenarios Analyses with MELCOR 1.8.2 and MELCOR 2.1 for the DEMO Helium-Cooled Pebble Bed Blanket Concept

As for Light Water Reactors (LWRs), one of the most challenging accidents for the future DEMOnstration power plant is the Loss of Coolant Accident, which can trigger the pressurization of the confinement structures and components. Hence, careful analyses have to be executed to demonstrate that the confinement barriers are able to withstand the pressure peak within design limits and the residual cooling capabilities of the Primary Heat Transfer System are sufficient to remove the decay heat. To do so, severe accident codes, as MELCOR, can be employed. In detail, the MELCOR code has been developed to cope also with fusion reactors, but unfortunately, these fusion versions are based on the old 1.8.x source code. On the contrary, for LWRs, the newest 2.1.x versions are continuously updated. Thanks to the new features introduced in these latest 2.1.x versions, the main phenomena occurring in the helium-cooled blanket concepts of DEMO can be simulated in a basic manner. For this purpose, several analyses during normal and accidental DEMO conditions have been executed. The aim of these analyses is to compare the results obtained with MELCOR 1.8.2 and MELCOR 2.1 in order to highlight the differences among the results of the main thermal-hydraulic parameters.


Introduction
The exploitation of fusion as energy source requires the demonstration of a limited impact in terms of risk to the staff, the public, and the surrounding environment, well below the limits established by the national safety authorities.Hence, a systematic safety analysis has to follow the design development, to demonstrate that the safety objectives are met for each solution proposed.
In Light Water Reactors (LWRs), these analyses are executed employing the so-called "severe accident codes."These codes are able to simulate accidents from the initiating event until the release of radioactive materials outside the containment building [1].The same codes may be also employed for fusion applications [2][3][4][5][6] but the intrinsic and deeper differences between the two reactor types cannot assure the quality and the correctness of the obtained results.For this purpose, some severe accident codes have been adapted to cope also with specific fusion phenomena.
One of the codes that have been adapted is MELCOR, in particular, the version based on the old, and no longer maintained, 1.8.2 version (M 1.8.2) [3].On the contrary, newer MELCOR versions are continuously released and developed for LWRs, and, in the latest versions, an extension to cope also with helium-cooled reactor has been implemented [7,8].Hence, these latest versions (M 2.1.x)should be also capable of treating the main phenomena occurring in the DEMO Helium-Cooled Pebble Bed (HCPB) concept in a basic manner.
DEMO is the first fusion reactor designed to prove the fusion plant capability to produce electrical power in a safe and commercially acceptable way.Several design solutions have been proposed for such reactor, based on different blanket concepts [9].Parallel to the development of the blanket concepts, also several different Primary Heat Transfer System (PHTS) designs have been developed, and the HCPB blanket concept and its PHTS seem to be one of the most promising solutions proposed [9,10].
2 Science and Technology of Nuclear Installations Therefore, the aim of the present work is to check the thermal-hydraulic capabilities of the two MELCOR versions (1.8.2 and 2.1.6342)in support of future researches and analyses on DEMO reactor.Three scenarios have been analysed: a normal operational scenario and two off-normal operational scenarios, characterized by the release of helium inside the Vacuum Vessel (VV).The present paper can be then considered as a parallel one to [11], but more focused on a MELCOR version-to-version comparison instead of on design analysis.

The Helium-Cooled Pebble Bed (HCPB) Blanket Concept.
The HCPB blanket concept is under development at the Karlsruhe Institute of Technology (KIT), as part of the European efforts on the fusion technology researches [10].This HCPB concept consists in a solid breeder material (Li 4 SiO 4 or Li 2 TiO 3 ) in form of flat pebbles encapsulated into boxes cooled by helium.Several boxes compose the blanket structure, and each box can be subdivided into three radial zones: (i) The First Wall made of EUROFER and protected through a Tungsten layer (W layer).The aim of this First Wall is to protect the outer regions of the machine and to provide a suitable heat transfer area for the plasma heat.
(ii) The Breeding Zone needed for the production of Tritium.The Breeding Zone is composed as a sandwich (in tangential direction) made of pebble beds, cooling plates/pipes, and Be pebbles for neutron multiplication.
(iii) The Back Supporting Structure (BSS) able to hold the box in position and to feed the box coolant circuits thanks to 4 manifolds.
The blanket is subdivided into 16 sectors each composed by three Out-Board (OB) and two In-Board (IB) segments.In total, 48 OB and 32 IB segments compose the blanket itself, and each sector covers an angle of 22.5 ∘ [12,13].Each segment is in turn composed by six poloidally arranged boxes.In Figure 1, a sketch of the boxes arrangement is shown.

The HCPB Primary Heat Transfer System (PHTS).
The conceptual cooling strategy for the HCPB Blanket is based on the use of four independent cooling circuits: two for the OB segments and two for the IB ones.This design increases the overall blanket heat transfer performances thanks to the cooling of the FW channels in counter flow.The two cooling circuits are connected to the blanket segments via the manifolds in the BSS.Each cooling circuit consists of five Cooling Trains (CTs) plus one spare for the OB and two CTs plus one spare for the IB.Each CT has a Steam Generator (SG) and a compressor (Figure 2).The PHTS design data reported in Table 1 are taken from [12] except for the temperature range which is taken from [13].
The connection among the blanket OB segments and the PHTS is made through 48 pipes crossing the VV upper port.The pipes size is DN 250 inside the VV upper port and DN  350 outside the VV upper port to reduce the velocity of the coolant and the pressure drops.The coolant is then collected by five hot headers (200 m 3 each), which in turn redistribute the coolant to the CTs (main pipes size DN 1200).Finally, five cold headers collect and redistribute the coolant to the blanket segments (Figure 2).The size of piping, compressors, SGs, and valves are the same for both the OB and the IB coolant loops and the overall piping layout relies on quite good standardization.Further details on this design can be found in [12][13][14][15].

Employed Models
3.1.Normal Operational Scenario.The PHTS spatial nodalization for the DEMO stationary scenario has been built based on the data reported in [12] and consists in the following: (i) 12 Control Volumes (CVs) employed to simulate the PHTS.The inlet and the outlet blanket pipes are simulated as a single lumped CV, as well as the main coolant pipes connecting the headers with the SGs.
(ii) 12 Flow Junctions (FJs) employed to simulate the connections among the CVs.
(iii) 2 Heat Structures (HSs) employed to simulate the heat exchanges between the FW and the box coolant pipes and the primary and the secondary SG sides.The total heat exchange area is 2,500 m 2 and 15,000 m 2 for the blanket and the SG HSs, respectively.As boundary conditions, injection of 910.5 MW and removal of A single CV with a free volume of 70,000 m 3 at 40 ∘ C and 0.1 MPa is employed to simulate the Expansion Volume (EV), which is the "containment" of a fusion reactor.These boundary data have been taken as in the previous Power Plant Conceptual Studies [16].Three HSs, characterized by an inner layer of steel (1 cm) and an outer layer of concrete (1 m), simulate the boundary walls of the EV.The sum of the heat exchange areas of these three HSs is 7,850 m 2 .In Figure 3, a sketch of the employed model is reported.

Off-Normal Operational Scenarios.
In the off-normal operation scenarios, the PHTS has been simulated as in the normal operational scenario, and a valve has been introduced to simulate the break across PHTS and the VV.The VV has been in turn connected with the EV through a rupture disk.The break is supposed to occur 100 s after the beginning of the simulations (at 0.0 s in the graphs), and two break sizes have been investigated: 0.12 m 2 and 0.01 m 2 .The rupture disk has a total flow area of 1 m 2 and an opening set-point of 0.15 MPa [17].The VV has been simulated with a free volume of 1,860 m 3 at 210 Pa and 200 ∘ C. The free volume value is an assumption based on the discussions performed with the DEMO VV designers.
The initial conditions of the system are the same as the normal operation scenario and, after the break, the following plasma decay power trend has been modelled [18]: (i) In 1 s, the plasma power falls down from 100% to 5%.(ii) From 1 s to 3,600 s, the plasma power remains constant at 5%.
(iii) From 3,600 s to 7,200 s, the plasma power decreases from 5% to 1%.
(iv) After 7,200 s, the plasma power remains constant at 1%.
In Figure 3, a sketch of the employed models is reported.

Normal Operational Scenario Results
As shown in Table 2, both MELCOR versions are able to simulate the main PHTS parameters.The slight differences shown in the blanket and SG temperatures and pressures are probably due to the improvements and modifications introduced in the M 2.1 version.Further details about these improvements are provided in the following paragraphs.

Off-Normal Operational Scenario Results
5.1.Break Size of 0.12 m 2 .The first off-normal operational scenario investigated is characterized by a break size of 0.12 m 2 , that is, double guillotine break of one of the four  manifold pipes passing through the BSS.The break is supposed to occur 100 s after the beginning of the simulation, that is, at 0.0 s in the graphs.The break involves only this pipe and the helium released is supposed to flow directly into the VV without any resistance provided by the surrounding blanket boxes.
As shown in Figure 4, the depressurization rate of the PHTS is quite fast, and in less than one minute an equilibrium condition is reached with the VV and EV.Both MELCOR versions provide similar results for the PHTS and EV, but some differences exist in the VV pressure trends.The VV behaviour can be subdivided into two phases: an initial fast pressurization phase followed by a depressurization phase, triggered by the opening of the rupture disk.The pressurization phase lasts for about 2 s in both MELCOR versions, and similar values in terms of maximum pressure and temperature are shown (0.66 MPa and 1050 ∘ C).The subsequent depressurization phase lasts for about 35 s in both MELCOR versions, but in M 2.1 a steeper trend is shown.This difference is due to the mass flow rates across the rupture disk (Figure 5) and to the VV atmospheric temperatures between 10 s and 20 s (Figure 6).Although it should be noted that several modifications were introduced in the latest MELCOR versions, such modifications do not provide any remarkable difference in terms of total pressure and mass flow rates.
Figure 7 shows the blanket, the SG, the EV, and the VV atmosphere temperature trends during the transient.For the first 500 s, the blanket temperature suffers cyclical increase and decrease due to the combined action of the plasma decay heat, the helium expansion, and the residual mass flow rate in the circuit.However, after 500 s, the influence of the plasma decay heat became the most important driving phenomenon, and the blanket temperature starts to increase monotonically.The two MELCOR versions begin to provide different results after about 300 s for the PHTS and about 100 s for the VV, that is, during the phase when the plasma decay heat starts to become the only important parameter.Inside the documents assessing the improvements and modifications introduced to M 1.8.2 for the ITER safety assessment [4,5], no mention is made regarding the employment of helium as coolant; therefore, the helium modelling approach should be the same in both MELCOR versions.Although minor differences in the thermal-hydraulic modules are present for the two code versions, their influence cannot be excluded: (i) In M 1.8.2, specific relations for nucleate boiling, critical flux, and film boiling were introduced [4].Of these three modifications, only the critical heat flux model can modify the results in the present scenario, because the other two models are employed only if "water" is employed as coolant.
(ii) In M 2.1, ad hoc modifications were introduced to simulate HTGRs, and the gas properties were reassessed [8].These new modifications introduced  are not specifically related to the thermal-hydraulic aspects, but the modification of the gas properties could have an impact on the results.
(iii) M 2.1.xare completely recodified versions, so it cannot be assured that the results of the two versions are the same even if the same models are employed.
Taken one by one, the previous causes seem to be not sufficient to explain all the differences in the gas temperature behaviours, so a combined effect of whole causes seems to be the most appropriate explanation.Moreover, these modifications have also a minor influence on the mass flow rates and thus on the total pressure trends, as shown in Figures 4  and 5.
Regarding the SG atmospheric temperature (Figure 7), except for initial oscillation due to the helium expansion, the overall temperature remains quite constant during the entire transient, but a difference of about 30 ∘ C among the results of the two MELCOR versions is shown.The reasons of this behaviour should be the same discussed for the blanket temperature behaviour.On the contrary, the EV temperature results are quite similar thanks to EV size, which dumps the influence of the differences between the two MELCOR versions.In the early instants, the EV temperature increases due to the helium ingress, and then it decreases thanks to the heat exchanges with the cold structures and, in the long term, with the outer environment.The maximum pressures and temperatures reached inside the EV are ∼0.27MPa (Figure 4) and ∼270.0 ∘ C (Figure 7), respectively.
A special highlight should be made for the VV temperature, because it is mainly influenced by the heat transfer coefficients (HTCs) among its free volume and the surrounding HSs (Figure 8).For the first 30 s, the blanket temperatures provided by both versions are quite similar; then M 2.1 starts to predict higher temperatures.This behaviour is due to the lower HTC values calculated by M 2.1.A similar behaviour was also shown in [19], in which a comparison between MELCOR 1.8.2 and MELCOR 1.8.6 against an In-Vessel LOCA in the ITER facility was performed.The reasons found behind this different behaviour were mainly related to the improvement on numerical methods and modelling approaches introduced in the 1.8.6 version [20].Therefore, considering that MELCOR 2.1 is successive evolution of the  code, the differences highlighted on the temperature and HTC behaviours can be explained in a similar manner.

Break Size of 0.01 m 2 .
The second off-normal operational scenario is characterized by a break size of 0.01 m 2 at the same location of the previous scenario.This case has been investigated to highlight the differences between the two MELCOR versions during a slow transient.As for the previous case, the helium released flows directly into the VV without any resistance provided by the blanket boxes.
Figures 9 and 10 show both pressure trends of the PHTS, the EV, and the VV with two different time scales.For this small break size scenario, the PHTS depressurization phase lasts for about 10 min instead of about 1 min as in the previous larger case (Figure 4).The rupture disk between VV and EV opens in ∼5.0 s in both MELCOR versions, but the subsequent PHTS depressurization is slightly delayed in time for M 1.8.2 due to a different prediction of the mass flow rate across the rupture disk (Figure 11).This difference leads also to different VV atmospheric temperature and pressure behaviours (Figures 10 and 12), but no influences  are highlighted in the EV behaviour (Figures 13 and 14).The different mass flow rate across the disk can be explained considering the notion that M 2.1 executes "double precision calculations," while M 1.8.2 executes only "single precision calculations."Quite similar results were obtained also in [19], where a comparison among M 1.8.2 and M 1.8.6 (both introducing the same models for fusion applications) was reported.
Figures 13 and 14 show the atmospheric temperature trends inside the PHTS and in the EV for MELCOR 2.1 and MELCOR 1.8.2, respectively.Both MELCOR versions show that about 400 s after the break the temperature of the whole system remains quite low, except for the blanket and the outlet pipes.This behaviour is due inversion of the mass flow rate across the blanket (BL) and its outlet pipes (OP) (Figures 15 and 16).In the early instants, the mass flow rate remains above 10 kg/s and the temperatures in the BL and its OP increase thanks to the plasma decay heat, while in the hot header (HH), the SG, the cold header (CH), and the blanket inlet pipes (IP), the temperatures decrease due to the helium expansion.In the following instants, the temperatures of the blanket and its outlet pipes abruptly decrease due to the decrease of the mass flow rates to very low values (∼0.0 kg/s).However, at about 450 s for M 2.1 and at about 425 s for M 1.8.2, the mass flow rate across the BL and its OP falls below 0.0 kg/s meaning that flow inversion is occurring.During this inversion phase, a cold fluid stream flows from the SG to the OP causing an abrupt decrease of their temperatures, while the blanket continues to remain hot thanks to the plasma decay heat.Once the mass flow rates across the loop fall to 0.0 kg/s, the BL temperature continues to increase due to the decay heat, while in the other plant zones it slightly decreases due to the termination of the helium expansion.Indeed, the total pressure equilibrium is reached at about 600 s, as well as the termination of the helium circulation in the loop, meaning that the total pressure and the mass flow rates are closely connected.
Although the SG shows a quite surprising behaviour due to its almost constant trend along the entire incidental transient, this behaviour is probably due to the presence of an HS, which reduces the effects of the helium expansion through the release of the heat accumulated during the stationary phase.For this purpose, to check whether the same effect can be reproduced also in the other zones of the loop, improved plant nodalisation has been built, characterized by the introduction of an HS in each CV.The thicknesses of these HSs were taken from [12] and an external layer made of mineral wool was considered along the entire coolant loop, to obtain a total heat loss toward the EV below 0.1%.As shown in Figures 15 and 16, these HSs prevent the mass flow rate inversion through the BL and the OP between 400 and 600 s.Similarly, the atmospheric temperatures in various CVs of the loop, instead of decreasing to about 100 ∘ C, increase to 300 ∘ C (IP and CH), to 500 ∘ C (HH), and to 600 ∘ C (OP) as shown in Figures 17 and 18.
The results obtained show clearly the influence of such HSs, but their effect is not strong as for the SG HS.This difference is probably due to greater heat transfer area of the SG HS compared to the other HSs and to the different boundary conditions imposed.The SG HS has a total heat transfer area of 15,000.0m 2 , while the heat transfer area of the biggest HS attached to a CV of the loop is about 800 m 2 , that is, 1/20 of the SG HS size.The boundary conditions imposed were as follows: (i) For the SG HS, the power removal rate for the entire transient has been user imposed.
(ii) For the other HSs, only the inner and outer surfaces were linked to their CV and to the EV, respectively, but no other conditions were imposed.
So, for future analysis, it seems appropriate to always add HSs to each CV in order to obtain more reliable results.

Conclusions
In the present work, the results of a comparison study performed with two versions of the MELCOR code for the DEMO reactor behaviour have been presented.The code versions employed are M 2.1 for LWRs and M 1.8.2 for fusion reactors [5,21].This second version (M 1.8.2 for fusion reactors) has been widely employed in the past for ITER and DEMO safety analyses [4,5,22,23], while M 2.1.xare extensively used for LWR applications, but in the newer versions their capabilities have been extended to cope also with HTGRs [7,8].Therefore, the newer M 2.1 version can be also employed to treat the main phenomena occurring in the DEMO HCPB blanket concept in a basic manner.
For this purpose, it has been found appropriate to check the differences among the prediction provided by the two MEL-COR versions.To do so, the results from these two versions during normal and two off-normal operational scenarios have been reported and discussed in the paper.An overall good agreement has been found during the analysis of the normal DEMO operational scenario, but different pressure, temperature, and mass flow rates trends have been highlighted for the two off-normal scenarios, especially for the small break one.However, it should be also highlighted that specific user assumptions could have influenced the results obtained, but if user effects are present, then they should be the same for both MELCOR versions because the same identical DEMO plant nodalisations were employed.The obtained results are not surprising considering the deep evolution of the MELCOR code in terms of thermal-hydraulic modelling and numerical methods.However, at this stage, M 2.1 version can be applied only to simple thermal-hydraulic analysis for the future fusion reactors due to the lack of the specific models, as the Plasma Facing Components (PFCs) oxidation under steam atmosphere, the water freezing, and the improvements on the aerosol deposition model, which on the contrary characterize the M 1.8.2 fusion version [4].
Nevertheless, M 2.1 provides more reliable thermalhydraulic results thanks to its numerical double precision.Therefore, in the future, it could be interesting to see the introduction of the specific fusion models also into the M 2.1 version in order to obtain an updated and more reliable integral code for fusion safety analysis.Another tip that can be also proposed for future analysis is to introduce HSs in each CV.

Figure 4 :
Figure 4: Total pressure trends of the PHTS, the VV, and the EV.

Figure 5 :
Figure 5: Mass flow rates across the break and the rupture disk.

Figure 6 :
Figure 6: Atmospheric temperature trends in the VV.

Figure 7 :Figure 8 :
Figure 7: Atmospheric temperature trends in the blanket, the SG, the EV, and the VV.

Figure 9 :
Figure 9: Total pressure trends of the PHTS, the VV, and the EV.

Figure 10 :
Figure 10: Total pressure trends in the EV and the VV during the first 25.0 s after the break.

Figure 11 :
Figure 11: Mass flow rate across the break and the rupture disk.

Figure 15 :Figure 16 :
Figure 15: Mass flow rates through the inlet and the outlet blanket pipes in M 2.1 for the base and the modified cases.