Journal of Nuclear Energy Science & Power Generation TechnologyISSN: 2325-9809

All submissions of the EM system will be redirected to Online Manuscript Submission System. Authors are requested to submit articles directly to Online Manuscript Submission System of respective journal.

Research Article, J Nucl Ene Sci Power Generat Technol Vol: 8 Issue: 1

Mesh Size Sensitivity of the Simple Cheap Vessel Benchmark Problem

Djillali Saad1* and Ahmed Hadjam2

1Department of Nuclear Materials and Mechanical Test, Nuclear Research Center of Birine, P.O. BOX 180, Ain-Oussera 17200, Algeria

2Thermal-Hydraulic Laboratory, Nuclear Research Center of Birine, P.O. BOX 180, Ain-Oussera 17200, Algeria

*Corresponding Author :Djillali Saad
Department of Nuclear Materials and Mechanical Test, Nuclear Research Center of Birine, P.O. BOX 180, Ain-Oussera 17200, Algeria
Tel: +2130550884021
E-mail: [email protected] or [email protected]

Received: December 11, 2018 Accepted: February 18, 2019 Published: February 22, 2019

Citation: Saad D, Hadjam A (2019) Mesh Size Sensitivity of the Simple Cheap Vessel Benchmark Problem. J Nucl Ene Sci Power Generat Technol 8:1.


A mesh sensitization study was applied to the Simple Cheap Vessel Problem provided by Innovation Systems Software (ISS). The generations of the different meshes used in the present work are realized using the CoupleMeshing program, locally developed. This program allowed us to reduce the input of the COUPLE mesh part composed of several hundreds of lines (formatted or unformatted) in only five unformatted lines. A new RELAP5/SCDAP execution interface of the code that integrates the CoupleMeshing program has been proposed. A call is being made through this article to Mr. C. Allison to support this new interface. The results found clearly show the influence of the mesh size on the different parameters studied.

Keywords: RELAP5/SCDAP; Couple mesh; Mesh sensitivity; Lower head vessel; Relocated debris; Molten pool


The RELAP5/SCDAPSIM code is designed to predict the behavior of reactor systems under normal and accident conditions. This code is in continuous development as part of an international program of development of nuclear technologies called SDTP (SCDAP Development and Training Program) [1]. The COUPLE model in the SCDAP/RELAP5 code is designed to calculate the heat up of a reactor vessel lower head in response to contact with reactor core material that slumps to the lower head [2].

The COUPLE model is largely inspired from the standalone COUPLE program developed at INEL. Numerical computation, of transient thermal conduction and associated heat transfer in two-dimensional structures, in the COUPLE model is based on the combination of triangular and quadrilateral finite elements. The COUPLE model is two-dimensional and its application is quite general, one of these famous applications is for the RELAP5 code. The primary application of the COUPLE model is the computation of the thermal behavior of the lower part of the primary pressure vessel and the molten material or debris that may have fallen to this region from the core region [3].

The COUPLE model does not use an automatic mesh generator and the finite element mesh must be defined manually (input data). The input data can contain a very general mesh but it would require a considerable input. Hence the need for the use of a semiautomatic mesh generation procedure allowing simpler input with a loss of generality is used.

Two related meshes, the simple mesh and the transformed mesh, are considered in the semi-automatic mesh generation. The simple mesh is a two-dimensional mesh of JM lines of nodes, with IM nodes in each line, and with JM IM nodes arranged in a square pitch. Finally, the simple mesh is transformed into a mesh of finite elements to represent the real geometry of the application.

The axisymmetric r-z model is applied to a cross section of the pressure vessel. The mesh takes into consideration the outer section of carbon steel with an inner lining of stainless steel. The transformed mesh will not, probably, have the square pitch of the simple mesh composed only of the vertical and horizontal straight lines connecting the nodes, but also with circular arcs positioned as needed to connect the nodes. But the node connections remain the same. That is, the four nodes forming the square in the simple square mesh are the same nodes forming the quadrilateral finite element in the transformed mesh [3].

Although it is called, in this simple mesh, to introduce the minimum of geometrical data necessary for the constitution of the transformed mesh, but this minimum remains always painful. Note here that this simple mesh can easily contain a hundred lines just for its part of the geometric data (part of the simple mesh is shown in Figure 1).

Figure 1: Left: Format currently used (old style) for introduction of the COUPLE mesh model; Right: Proposed format.

The problem raised in the preparation of the input file for the code RELAP5/SCDAP Mod 3.x and earlier versions, especially in its part SCDAP is undoubtedly the part of the COUPLE model. This part of the input deck can easily contain several hundred lines (depends on the density of the mesh). In addition to the formatted format of the data (RELAP5/SCDAP Mod 3.x and earlier versions), the simple mesh of the lower part of the reactor vessel is very difficult. The manual preparation of this simple mesh can easily contain errors in the position of the nodes without the possibility to identify them. This difficulty is increased when one has to make several times a different simple mesh to see for example the effect of the sensitivity of the mesh in the calculations.

Although in the RELAP5/SCDAP Mod4.0 version is able to use the new style of COUPLE format. The new format increases user’s flexibility and reduces the input errors when users write input decks. Certainly, in this new format, the formatted style of data has become unformatted and therefore more flexibility, but the problem of simple mesh preparation remains unresolved. In this new format, the style has become unformatted instead of formatted. The beginning of each line of data is started by the 5M0xxxxx card. There is the possibility of inserting comments. All this does not solve the main problem which is the large number of lines (hundreds) that must be written for each single mesh.

To solve the equations, the discretization (the mesh) is obligatory. The main goal of a mesh of finite elements is to bring the geometry of the modeled body closer. Space steps are therefore of great importance for the accuracy of the numerical solution. However, the more the mesh will be refined, the more computation time will be important. We must therefore find a compromise between accuracy and speed of calculation. For that, tests have been made with several meshes to see from what moment the refinement brings no more notable improvements.

In most cases, its tests are performed on applications (software) that support automatic mesh. However this task becomes painful when the preparation of the mesh is done manually or it is semi-automatic as in the case of the COUPLE model. Note here, that his tests are done several times until the desired result. The number of tests depends in general on the type of problem considered (tens of tests in most cases).

In order to overcome this difficulty and make this study possible, we thought about automating the simple mesh preparation and all the necessary data for the COUPLE model part. A program in Fortran 90 language “CoupleMeshing” is realized. This application allowed us to transform the part of the input deck of the COUPLE model consisting of several hundreds of formatted (or unformatted) lines to five (05) unformatted lines only (Figure 1).

In-Vessel Severe Accident Phenomena

Various scenarios cause a severe accident; it is characterized by core degradation i.e. loss of geometrical integrity by melting or debris formation [4]. During a station blackout (SBO), the core is insufficiently cooled, therefore the residual heat is not removed completely and the water continues to boil which leads to the discovery of the core. At the beginning of the accident, the excessive core heat up of the core is indicated by the axial and radial power profiles of the decay power and by the instructions imposed on the minimum level of water in the core [5]. As a result, the core temperature continued to increase thus threatening the integrity of the structural materials and eventually the fuel, subsequently causing core damage. Two different phases characterize the core degradation:

An early phase that covers, the beginning of the core uncovery, the heat up process and the melting of the reactor materials with a relatively low melting points (Zircaloy cladding, B4C absorber material). However, the geometry of the core remains essentially intact.

And a late phase that leads to the melting of ceramic materials (UO2), the formation of molten pools, the relocation of debris to the lower plenum and eventually to the rupture of the lower head [6].

In general, the early phase of core degradation is dominated by thermohydraulic phenomena. Hence, the importance of choosing a better thermohydraulic modeling for the reactor response and the timing of the event. The late phase is dominated, in addition to the thermohydraulic phenomena, by additional phenomena such as the oxidation and the interaction of the materials of the core. Severe accidents codes are more diversified, less validated and hence have more uncertainties in the late phase [7].

Core degradation and relocation to lower plenum

The mechanical failure of the ZrO2 oxide layer of the cladding is the main cause of the initiation of debris relocation. This oxide layer acts as a protective wall and prevents slumping of fused materials into the core. For low-pressure accidents, cladding in zircaloy begins to swell and rupture occurs when the core temperature reaches 1000- 1200 K [8]. The release of volatile fission product gases is a direct consequence of the fuel cladding rupture. These products will have serious consequences if they are released into the environment. The candling phenomenon is formed by the simultaneous melting of control blades and fuel rods that move downwards in colder areas and freeze to form a crust [5]. The oxidation processes are affected when the re-solidification of the liquid materials is triggered thereby causing flow blockages that reduce the local fluid flow. Hence, steam starvation zones are created that cannot be oxidized and are likely to produce Zircaloy rich melt that will not be oxidized until relocation to another part of the core [8]. Local blockages also prevent the melt of new relocations, which allows the formation of local molten pools. The potential for molten metal relocation through the core support plate when it is still intact was identified in the XR2-1 experiment performed at SNL [9]. Before the central support plate is completely plugged, some of the fusion control rods will be able to move into the lower plenum. The central core support plate, whose function is the separation of the core region to the lower plenum, finally yields due to the accumulation of the debris charge followed by the ejection of a large amount of debris into the lower head.

Debris characteristics in lower plenum

Numerical simulation of the debris bed in the lower plenum is very difficult due to the complex physical phenomena involved. At the present state of knowledge, the relative quantities, the location of the phases and the layers in the lower part of the head cannot be described deterministically [8]. Figure 2 illustrates some important phenomena involved in the configuration of debris in the lower head.

Figure 2: Phenomena affecting the lower plenum debris configuration [8].

The amount of water available in the lower plenum plays a key role in the consequence of the relocation of any mass to the latter. Thus hot debris comes into contact with the remaining water, undergoes coarse fragmentation (particulate debris) that can lead to non-energy interaction (steam spike) or cause an energetic steam explosion and severely damage the lower head when remaining metals becomes oxidized [10]. The composition of the debris quenched in the lower part of the lower head will vary with height because of the different properties of the material as the melt builds up [11].

Porosity is the key parameter in particle bed behavior. The cooling ability of the debris depends inter alia on the space between the particles. Experimental tests and analytical calculations have confirmed that the porosity of randomly compacted spheres is about 40% independent of particle size [12]. Based on data from the TMI-2 accident, the size of the entrained particles varies between 1 and 5 mm [13] and the porosity shows only a small dependence on the variance of the size distribution [12]. After a complete failure of the lower core support, the interaction with the water of the lower head is limited in the case where the melt enters the latter in the form of a large single jet (diameter of about 50 cm). The relocated melt will remain largely liquid and will immediately form a molten pool surrounded by crusts in the lower head, possibly with an overlying water layer on top [5]. The stratification of the molten material essentially depends on the relative timing and nature of the relocation process. The early relocation of ceramic materials promotes the formation of several layers of ceramic and metal layers in the lower plenum, while late relocation promotes mixing and reducing the number of layers [8]. The failure of the vessel and the progression of the ex-vessel accident (e.g. steam explosion and coolability) are conditioned by the differences in composition and the amount of molten liquid due to the different thermal conductivities [14].

The integrity of the lower head during the TMI-2 accident was explained by an inherent cooling mechanism due to the gap formation between the debris and the lower head wall [8,13]. Microscopic cracks, which develop in the wall of the lower head, trap water and during its evaporation, when the debris move towards the lower plenum, it produces a contact resistance gap with a thickness typically between 1 and 2 mm [15]. This result is confirmed by calculations carried out on SCDAP/RELAP5-3D, taking into account the difference between the corium-to-vessel-gap compared to a non-existing gap, these calculations show a significant decrease in the temperature of the lower surface [15].

In-vessel hydrogen generation

The oxidation process of the core materials (zircaloy, stainless steel and B4C) is responsible for the production of hydrogen. Oxidation accelerates the core temperature rise during the progression of severe accident. In terms of hydrogen production and core degradation, zircaloy and steam oxidation is considered the most important process [8]. In the 1500 to 1700 K temperature range, the exothermic heat of zircaloy oxidation exceeds the disintegration heat, and the cladding heat up rate increases significantly (1 K/s to 10-20 K/s) [16]. The production rate of hydrogen is mainly conditioned by the participating masses, temperatures, surfaces areas and the availability of water and steam [6]. The factors that dominate hydrogen production are the temperature at which the ZrO2 layer breaks and the underlying molten oxidizing metal is released [17]. Due to the complexity of Zircaloy surface estimation and the effect of blockages, a major uncertainty in the determination of hydrogen generation is the timing of the cladding failure and loss of core geometry [6,10]. However, oxidation of B4C is highly exothermic and generates 6 to 7 times more hydrogen than oxidation of the same amount of Zircaloy [18]. The contribution of the oxidation of B4C is less important than that of Zircaloy because of the amount of B4C much lower than the mass of Zircaloy present in the reactor.

One of the lessons to be learned from the Fukushima accident is that the accumulation and subsequent detonation of hydrogen gas produced by rapid oxidation can break down the containment structures and lead to widespread radioactive contamination [19]. The hydrogen mitigation system used in the inert containment significantly reduces the risk of hydrogen combustion.

The Simple Cheap Vessel Problem

The Simple Cheap Vessel Problem is the most rigorous of the benchmark problems, in that the core is modeled with two groups of fuel rod and two groups of Ag/In/Cd control rod components, and a COUPLE mesh has been added, as shown in Figure 3. Each group of fuel rods represents 18,408 rods and each group of control rods represents 118 rods. These groups of fuel and control rods are representative of an entire reactor core. This problem is intended to test the meltdown of the core, and the transfer of core molten material to the lower head [20,21]. The core volumes are initialized with steam, and allowed to heat up, meltdown, and relocate to the lower head.

Figure 3: Nodalization diagram for Simple Cheap Vessel Benchmark Problem [20].

The nodalization of the problem studied is illustrated in Figure 3. Through the time-dependent junction 50, the flow of fluid begins in the time-dependent volume 10. The lower head of the container is represented by the single volume 60 which is connected, in turn, to a bypass volume 80. The six pipe volumes 100, 101 and 105 model the two center channels and the bypass channel respectively. Each volume of these pipes is connected by a crossed junction. Timedependent single volumes 110 and 125 model the upper plenum of the reactor system [22].

As for the SCDAP model, it includes four components two for fuel and two others for control rods. The core model is PWR and the power history used is based on the TMI-2 reactor [23]. In the SCDAP model, the fuel and control rod components are subdivided into six thermal structure nodes. The fuel nodes of the SCDAP model have flow boundaries from volumes in pipes 100 and 101. Simulation heat transfer between fuel and control rods components is supported by the SCDAP radiation model. The creep rupture of the vessel is also modeled. The power for the SCDAP heat structures uses a general table 900, and has a value of 934 MW.

The CoupleMeshing program uses unformatted input. The axisymmetric model is used to model the lower head of the vessel. Two configurations of discretization were used in this study; 11 and 22 nodes in the r-direction and 16 and 27 nodes in the z-direction. Modeled materials include relocated debris, stainless steel, Inconel, carbon steel and coolant. RELAP5 single volume 60 is used for the lower head flow boundary volume. It can hold slumped material during the transient.

Advantage of the CoupleMeshing program

The “CoupleMeshing” program is made using the Fortran 90 language. The GUI is implemented using the DISLIN library [24]. The interface of the “CoupleMeshing” program is illustrated in Figure 4. Since this application is capable of generating the mesh of the lower part of the reactor vessel necessary for the execution of the COUPLE model, it is also capable of introducing, to the same input file, the thermohydraulic data of the analyzed problem. The input file thus generated is obtained from a simple text file composed of only five unformatted lines.

Figure 4: RELAP5-SCDAP-COUPLE Interface.

The first line of this input file contains the format style for generating the COUPLE mesh Model and the title of the problem being processed. If the first word is “old”, the file is generated according to the old style format (formatted). If the first word is “new”, the file is generated in the new style format (unformatted) (Figure 1).

The second line of the input file describes the control data specific to each COUPLE mesh. Namely, the COUPLE flag, the number of the RELAP5 volume to receive COUPLE debris, the debris source indicator for COUPLE mesh, the flag for breakup of COUPLE debris, the maximum number of RELAP5-SCDAP hydraulic time steps to use per COUPLE time steps and the maximum COUPLE time step.

The third line of the input file groups the geometric data and the instructions necessary to execute the mesh. Figure 5 schematizes the lower portion of the pressure vessel and the geometric data necessary for the graphical representation of the model. These data are:

Figure 5: Schematization of the lower portion of the pressure vessel and the geometrical data needed for mesh generation

The center of the radius of the curvature for the three arcs of the structures constituting the lower portion of the pressure vessel. The first arc represents the outer radius of the vessel for its carbon steel part. The second intermediate arc is at the same time the inner radius of the carbon steel portion of the vessel and the outer radius of the stainless steel portion. The third arc represents the inner radius of the lower portion of the pressure vessel for its stainless steel part.

The thickness, at the cross-section, of the lower portion of the pressure vessel for the two materials: carbon steel and stainless steel.

The height of the lower portion of the pressure vessel.

We adopted the same methodology used by the COUPLE code. According to this methodology, two indices (i and j) are used to define each node of the mesh. The index i is in the radial direction and the index j in the axial direction. The proposed program “CoupleMeshing” greatly simplifies the task and transforms the twodimensional model into one-dimensional. It suffices to choose only the desired number of divisions in the radial direction, the remainder for the axial direction is calculated automatically (Figure 1, the division number is chosen equal to 13).

As for the discretization in the axial direction, we propose through this program a sophisticated method, but simple, to fully exploit the capabilities of this application. To discretize the thickness of the part of the carbon steel vessel, for example, it is sufficient to choose the number of subdivision and a multiplication factor (factdivcs: multiplication factor for carbone steel) (Figure 2). The multiplication factor allows us to refine the mesh near the outer wall thickness (0 < factdivcs < 0.5) or inner wall thickness (0.5 ≤ factdivcs < 1) or simply subdividing the thickness equitably (factdivcs=1) (Figure 6).

Figure 6: The different possibilities of discretizations through the thickness (carbon steel or stainless steel).

The last two lines of the input file contain the mechanical data of the materials used in the simulation. Namely, the density of the material, the horizontal thermal conductivity, the axial thermal conductivity, the specific heat capacity, the initial temperature of finite element mesh and the gap heat transfer coefficient for case of materials in both sides of the gap being solid state.

The rest of the data required to create the COUPLE mesh model file is taken by default as recommended by the RELAP5/ SCDAPSIM code manual [3].

“CoupleMeshing” is a light program that can be easily integrated into the RELAP5 execution interface code (Figure 7). Instead of adding hundreds of lines of data at the end of each input file (sample.i) to represent the COUPLE mesh model, just read, with this new interface, only the five lines (stored in a text file , sample.txt) and launch the mesh generator through the “Generate Couple Mesh” button. The operation will create a mesh file with all other data for the COUPLE mesh model and insert it at the end of the input file (sample.i).

Figure 7: Proposed new RELAP5-SCDAP Interface.

Of course, this new interface is not yet realized without the agreement and support of the ISS. Let us note here that the tools necessary for the realization of this new interface are carried out namely the mesh generation program “CoupleMeshing”.

Results and Discussion

After more than twenty different meshes tested on our problem, using CoupleMeshing program, we concluded that the results found are divided into two large families. From 1 to 11 radial divisions made on the lower portion of the pressure vessel give almost the same result. The second group of results is for divisions ranging from 12 to 22. For reasons of limitation of our version of the code RELAP5/ SCDAP Mod. 3.4, we could not test other more refined mesh. For this purpose, the results presented below are for two types of mesh: a coarse mesh (150 elements and 176 nodes), which corresponds to 11 radial divisions, and a refined mesh (546 elements and 576 nodes), which corresponds to 22 radial divisions (Figures 8 and 9).

Figure 8: Coarse mesh with 176 nodes and 150 elements.

Figure 9: Refined mesh with 594 nodes and 546 elements.

Figure 10 shows the evolution, as a function of the height of the lower head of the vessel, of the debris material indicator. The type of material is identified according to an indicator that varies between 0.3 and 1 (Table 1). The debris material indicator is defined by the AFBULK variable which is not taken into account by default. To take it into account, the card 208 is used.

AFBULK Type of material
0.3 Mostly Ag-In-Cd
0.4 Mostly stainless steel
0.5 Mostly Zr
0.6 Mostly ZrO2
0.7 More than 50% UO2
1.0 More than 70% UO2

Table 1: Indicator of type of material in element [3].

Figure 10: Debris relocated at the center of the lower head.

The bottom of the lower head of the vessel (less than 0.2 m) contains mostly Ag-In-Cd, while the rest contains mostly stainless steel. This result is justified because the mass of stainless steel relocated to the lower head of the vessel is much larger compared to the masses of UO2 or ZrO2 (Figures 11- 13). The two types of mesh give the same results for a height that varies between 0.3 and 1.2 m. The results are distinct for the bottom and the top of the lower head of the vessel. The coarse mesh shows a sudden change for the bottom of the lower head of the vessel, while the refined mesh shows a gradual change. This result is justified by the mesh step used in the two cases of the mesh. For the top of the lower head of the vessel, the results are quite different. The first type of mesh indicates that the top of the lower head of the vessel contains in the majority of the stainless steel. While the top of the lower head of the vessel is empty (filled with water) for the second type of mesh.

Figure 11: Total liquefied steel mass.

Figure 12: Total liquefied UO2 mass.

Figure 13: Total liquefied ZrO2 mass.

Figure 14 shows the evolution of the average temperature, in both cases of mesh, debris as a function of time. The endothermic peak that corresponds to 55 s, in the case of the coarse mesh, from which the temperature tends to zero, is explained by the quality of the mesh used in the simulation. Indeed, when the temperature, in some meshes of discretization, is badly estimated and when the number of these meshes is small (coarse mesh), a great influence on the average calculation of the temperatures is expected. This influence can be underestimated as around 55 s, as it can be overestimated (55 to 90 s). It is clear that mesh quality plays a very important role in any finite element calculation.

Figure 14: Average debris temperature in COUPLE mesh.

Figures 15 and 16 shows the evolution as function of time at the maximum temperatures of the debris and the lower head structure respectively. It is obvious that the temperature of the debris proportionally influences the temperature of the surrounding structure. This proportionality is well illustrated in the two figures cited above. Figure 17 shows the evolution, as a function of time, of the temperature of the outer surface cladding of the Zircaloy. This figure is calculated according to the SCDAP model of the code. In this figure, we notice a significant increase in temperature around 90 s, probably this sudden increase is due to the cladding rupture. This result is confirmed by the level of damage given in Figure 18.

Figure 15: Maximum debris bed temperature in COUPLE mesh.

Figure 16: Maximum temperature of structural material in COUPLE mesh.

Figure 17: Cladding outside surface temperature.

Figure 18: Level of damage at fuel rod component (channel #1).

In this figure, we notice the total rupture of the fuel rod in the axial volumes 2, 3 and 4. The level of damage in these axial volumes is reached 1 which corresponds to a molten pool (Table 2). Figure 18 also confirms the moment of rupture, which is close to 90 s. The two cases of mesh (coarse and refined) take this phenomenon into consideration. The behavior of the system (debris+the structure of the lower head) is related to what happens above. In the case of the coarse mesh, the peak of temperature is overestimated compared to the case of the refined mesh ( Δt=260 K).

0.0 Intact geometry
0.1 Rupture due to ballooning
0.2 Rubble (fragmented)
0.4 Cohesive debris
1.0 Molten pool

Table 2: DAMLEV Damage state [3].

When core meltdown is caused by a severe accident, the energy transferred from the melt to the coolant is considerable (Figures 19 and 20). An accurate estimate of the temperature and pressure in this part of the reactor helps immensely in any design calculation of this important part of the reactor, in order to manage the accident in good conditions. The accuracy of the numerical solution is automatically linked to the step in space used. The more the mesh is refined, the more precise the results are. However, beyond a certain limit of the step used, the results will converge towards a single limit. We must therefore find a compromise between precision and calculation costs.

Figure 19: Total energy transfer from debris and structural material to fluid at boundaries of debris and structural.

Figure 20: Molten pool power density.


A mesh sensitization study, on parameters of the molten pool, is carried out. The study is applied to the Simple Cheap Vessel Problem provided by ISS. A program in FORTRAN 90 language (CoupleMeshing) is realized to facilitate mesh generation of the COUPLE mesh model. This program is easily integrated into a new executing interface of the RELAP5/SCDAP code. A new form of this interface is proposed in this work. The program is tested in several mesh configuration cases. The program is able to automatically support the radial division number and mesh refinement in areas that are close to the outer (or inner) surface of the lower head vessel. This program allowed us to gain a considerable time in the constitution of the simple mesh for the COUPLE mesh model, which allowed us to test tens of mesh in a reasonable time. The results found for the calculated parameters, show the considerable influence of the mesh quality.


Track Your Manuscript

Recommended Conferences

Share This Page

Media Partners