1. Introduction
The PMM has found many applications in different fields of science since the introduction of Darcy’s law [1]. Seepage through soil, packed beds [2], and even biological tissues [3], flows in tube bundles and alloy solidification [4], and cooling of electronic components [5] are only a few examples. There are many analytical and numerical solution methodologies for this wide range of applications [5-8]. However the following two criteria specify the best method for a certain application.
1. The ratio of sensible length scale of the pores to length scale of the porous region or solution domain.
For many applications, such as flow through filters, this ratio can be in micro-units (10-6) or even smaller [6]. In this case, it is not practical to simulate each pore directly. Therefore all further efforts were focused on improving simulation in the interior of the porous medium by including the effects, such as thermal dispersion [9], variable porosity [10], and non-equilibrium thermal conditions [11].
But this is not the case for a nuclear reactor because the ratio is macroscopic [12]. In this case, all efforts should be focused on obtaining better results on the boundaries of porous media (the FA regions). Afterwards, a detailed simulation of a single FA can be devised and by using the same results on the boundaries of this model and all the lost data would be retrieved. This is the basic logic behind the modifications presented in this study.
2. Amount of porosity (φ) of the medium
A saturated porous medium is a combination of a fluid and solid portion. Porosity is the volume fraction of the fluid portion to the total volume [4]. For a high porosity medium, such as the core of a nuclear reactor, the solid matrix will exert force on the fluid in a complex manner. Therefore, it is not possible to relate velocity and force field with a simple analytical relation like a low porosity medium [5]. In this case, these forces should be measured by appropriate empirical or numerical means.
According to these criteria, there are three major steps that should be generally taken to fulfill the simulation of a high porosity medium with the macroscopic pore size ratio.
1. Measuring forces that are exerted on the passing fluid by the solid matrix at various mass flow rates using suitable empirical or numerical means and calculating the resistance coefficients based on these data.
2. Modifying governing equations of pure fluid to include porosity and adding the measured forces to the momentum equation as an external source.
3. Solving the modified form of governing equations on the solution domain that is elaborately divided into the porous and pure regions.
The method of performing these steps on a nuclear reactor with the proposed modifications is explained in Sec. 2. This section also includes the describing methodology of this study for the simulation of natural convection. The neutronic and thermo-hydraulic codes used and the algorithm of combining them are presented in Sec. 3. It includes introducing a novel CFD-based code that is time-dependent and three-dimensional and it has no limitation as to the details of the geometry of the solution domain.
The combined-code with the new PMM is used for the simulation of natural circulation in the Tehran Research Reactor, a typical MTR. Section 4 gives a brief description of this reactor. The specific requirements for simulation, such as defining the geometry of the solution domain are presented in Sec.5. Despite all the sophistication and details, it is still possible to solve this problem analytically. This is the reason for selecting this case to verify the modified PMM. The novel semi-analytical method for this purpose is developed in Sec. 6. The results and discussions are presented in Sec. 7. These results demonstrate the considerable capabilities of the proposed model.
2. Methodology
2-1. Measuring forces and evaluating resistance coefficients (step 1)
In a common nuclear reactor, coolant generally flows along the axis of fuel assemblies, through the core. This is the dominant direction of flow. Therefore, the porous media that are replaced with the core are not isotropic and the measured forces should be adapted to the following general form (in the Cartesian coordinate system) [5]:
where Si is the external source that should be added to the momentum equation and Dij and Cij are the tensors of viscous and inertial resistance coefficients, respectively. If the positive direction of the z-axis is assumed to be along the main direction of flow, eq. 1 can be restated as follows [12]:
where Fz is the force exerted on the fluid by the solid matrix along the z-axis and ∆p is its equivalent pressure drop. V is total volume of the porous region and L is the length of this region along the main direction of flow. The force is generally presented as a table of pressure drop data against mass flow rate or average velocity of flow. If such a table of data is available, a trend line can be plotted through them to yield a quadratic polynomial function in the following form, whereby the coefficients α and β are calculated.
By comparing eqs. 2 and 3, the unknown resistance coefficients can be calculated as follows:
This is the general process of evaluating resistance coefficients for a high porosity medium with a macroscopic pore size ratio and with one-directional flow. Ref. [12-13] used PMM for a power reactor and a research reactor, respectively. The procedures of these references for evaluating the coefficients is reproduced as the conventional approach in fig. 1(a). It is a conception of the above general process in its simplest form. The modified approach of this study is represented in fig. 1(b). The key features of the applied modifications can be declared by the following remarks.
(1) Including any variation in the cross-sectional area of the solution domain, vertical to the main direction of flow. It includes spacer grids, nozzles, headers, etc. Despite how small these variations may seem, they have significant effect on the accuracy of pressure drop calculations.
Similar to the current study, Ref. [12] used a numerical method with a CFD approach for calculating pressure loss against velocity. However, it reported a 25 percent error in comparison with the experimental results [12]. Investigations of this study revealed that the main cause of this considerable error roots in neglecting this remark. A quantitative proof of this matter is presented in Sec. 5.2.
(2) The incorporating solution of the energy conservation equation and defining all of the thermo-physical properties as functions of temperature.
In the conventional approach [12-13], the energy equation is not solved and by assuming iso-thermal conditions. Density and viscosity are evaluated at some presupposed average temperature (which is unknown prior to the solution). These values remain constant whether the problem is time-dependent or not. Whereas, in the modified approach, the energy equation is solved and all thermo-physical properties are defined as functions of temperature. This modification has a notable effect on improving the accuracy of evaluating pressure loss (see Sec. 5.2 for a quantitative proof) and also provides obvious advantages for accident analysis.
(3) Defining the resistance coefficients as functions of mass flow rate or average velocity of the flow
In the conventional approach [12-13], the resistance coefficients are always constant, whether the problem is time-dependent or not. It is due to iso-thermal assumption whereas in the modified approach, by solving the energy equation, it is possible to calculate average density and viscosity at each velocity. It results in calculating two specific coefficients at each velocity by using Eqs. (4) and (5). For many applications, such as accident analysis and simulation of pure natural convection, the mass flow rate is either unknown or changes during the solution. Furthermore, even in the case of a simple steady-state problem, the average temperature which is required to calculate average density and viscosity for Eqs. (4) and (5) is unknown prior to the solution. Therefore, in all cases, it is not accurate to choose two constant values, prior to the solution. The method used in this study fits a polynomial function directly to the data of the resistance coefficients against velocity (see Fig. 1(b)). The resulting correlations can be used for any application.
2-2. Modifying governing conservation equations (step 2)
As noted in Sec. 1, the core of a common nuclear reactor is generally categorized as a high porosity medium. When the quantity of porosity is close enough to unity, it is a common approach to disregard it in all conservation equations. It means that all of the solid obstacles in the porous medium are ignored and a pressure drop is simulated artificially by enforcing an external source in the form of Eq. (2) to the momentum equation. This approach is adopted by Ref. [12-13] for applying PMM to nuclear reactor analysis. Although it usually provides a good estimation of pressure and temperature variation this is not the case for the velocity field. Ignoring the increase in velocity due to existence of the solid matrix would also affect the boundary values of velocity that, according to Sec. 1 are curial to be calculate with the utmost precision for a nuclear reactor analysis. Therefore, all of the thermo-hydraulic conservation equations should be modified for porosity. The final form of the RANS equation [14] that is modified for porosity (φ) and used in this study is as follows (written for convenience in tensor notation).
In which, effective heat conduction coefficient (keff) is:
where, subscripts f and s refer to fluid and solid portion, respectively, Si is the external momentum source defined by eq. 2 and Sh is the heat source evaluated by neutronic equations. Definitions of other parameters can be found in Ref. [14] or any other reference that explains Reynolds-averaging approach for turbulence modeling.
2-3. Preparing reactor core and the algorithm of conducting solution (step 3)
In the conventional approach [12-13], the whole reactor core is simplified to a single porous medium, whereas in this study, the core is divided into several media that each comprises an FA. The most obvious advantage that this modification provides is as follows.
As noted in Sec. 2.1, the resistance coefficients are calculated for a single FA. These coefficients are verified by applying them to a porous region with the same scales and relevant conditions as that of the FA (see Sec. 7). By this modification, the same porous regions are placed into the core, one by one, without any changes. Therefore the same insured level of accuracy would be obtained in the main simulation. PMM with all the modifications that have been introduced so far will hereafter be referred to as the "multi-region PMM".
The general algorithm of performing multi-region PMM is represented in Fig. 2(b). The key feature of this model is that "it greatly reduces the computational resources required to conduct a large scale three-dimensional simulation without discarding any geometrical details". This virtue is achieved by imitating the common procedure for neutronic analysis of the reactor core (shown in Fig. 2(a)). The concept of cell-calculation in reactor analysis evaluates macroscopic cross-sections of an FA. Then the core is divided into several uniform regions where external dimensions of the FAs and the calculated cross-sections are applied to them [15].
There is an obvious resemblance between this procedure and the new PMM proposed in this study. In other words, the multi-region PMM is specifically devised to be the thermo-hydraulic equivalent of the neutronic procedure. This process starts with a detailed simulation of an FA and calculating correlations of resistance coefficients, as described in Sec. 2.2 (Fig. 1). Then the core is divided into proper regions. Because these regions are the same for both neutronic and thermo-hydraulic analyses, they are solved simultaneously (see Sec. 3.3 for detailed algorithm of coupling them). It is obvious that this process significantly reduces the computational cost of a large-scale simulation; nonetheless it may seem that this virtue is gained at the cost of losing all details inside FAs, because these regions are homogenized. However, as mentioned in Sec. 1 and explained in the preceding sections, all efforts in the new PMM are focused on obtaining more accurate results on the inlet and outlet boundaries of FAs, then after completing the large-scale simulation, (as Fig.2(b) suggests) these results are mapped onto the corresponding (inlet and outlet) boundaries of the detailed model of an FA and whereby all the lost data are retrieved. It makes no difference whether the problem is time-dependent or not. It should be noted that although the geometries of the inlet and outlet boundaries of the detailed model of FA and its equivalent porous medium (from the large-scale simulation) are almost always identical, but extrapolation should be anticipated for the data transfer process between them due to the fact that the applied mesh to these boundaries could be different.
2-4. Natural convection model
When heat is added to a fluid, temperature gradient causes density differences in that fluid. The buoyancy force formed by the density difference can induce a flow in the fluid. This buoyancy-induced flow is termed "natural convection" [16]. The following steps should be taken to include the effects of natural convection in a CFD simulation [14] or any similar approach that solves general conservation equations on a free geometry.
· Including the buoyancy term (ρg) in the momentum conservation equation (Eq. (7)).
· Solving the energy conservation equation (Eq. (8)).
· Defining density (ρ) as a function of temperature in a proper manner.
There are two approaches to perform the last step:
1. The Boussinesq approximation [16]. In this approach, density is assumed to be a constant value (ρ0) in all conservation equations, except for the buoyancy term in the momentum equation. In this term, density is defined as a linear function of temperature.
The constant density, ρ0 is evaluated at a presupposed average temperature, T0. Slope of the linear function is the thermal expansion coefficient (β) that is also an approximate value, due to the unknown range of temperature variation, prior to the solution.
2. If a thermohydraulic table is available for the active fluid, it is possible to extract the variation of density against temperature directly from the table at the operating pressure. Obviously, this approach is expected to have better accuracy. The method used in this study is to fit a polynomial function through these data and then use the resulting correlation instead of density in all conservation equations.
3. The applied neutronic and thermo-hydraulic codes
3-1. Introduction of the new CFD-based thermo-hydraulic code
Table 1 gives a brief review of the numerical method and special treatments implemented into the new CFD code along with the references in which they are explained. Many modifications are applied to the original methods mentioned in the table. These modifications are mainly focused on two points; reducing computational resources required and applying these methods to an unstructured and hybrid mesh grid. The most noteworthy modifications are:
| Title | Applied method | ||
|---|---|---|---|
| Numerical method | Finite volume [17] | ||
| Solver approach | Pressure-based (Segregated) [18] | ||
| Turbulence model | Simple | General | Standard k-ε [19] |
| Wall modeling | Standard wall-function [19] | ||
| Advanced | General | Realizable k-ε [20] | |
| Wall modeling | Enhanced wall-treatment [21] | ||
| Linearization shape function | Second-order linear function [17] | ||
| Method of evaluation of dependant variables on faces of a numerical cell | Convection terms (except for pressure term in the momentum equations) | Second-order upwind [22] | |
| Diffusion terms | Central-difference [23] | ||
| Pressure interpolation scheme (including the momentum equations) | Natural-convection dominated problems | Body force weighted model [16] | |
| Other problems | Central-difference [23] | ||
| Gradients evaluation method | Least squares cell-based [24] | ||
| Gradient (slope) limiter | Differential type [25] | ||
| Pressure-velocity coupling scheme (pressure-corrction equation) | Steady-state problems | SIMPLE [26] | |
| Time-dependant problems | PISO [27] | ||
| Under-relaxation method | Variables (explicit) relaxation scheme [17] | ||
| Solution accelerator | Algebraic Multi-Grid (AMG) with V-cycle recursive procedure [28] | ||
| solution method of Linear Algebraic Equations (LAE) set | LU decomposition with partial implicit pivoting [29] | ||
• Improving the method of applying the least-square scheme for gradients evaluation [24].
• Modifying the PISO scheme for pressure-velocity coupling [27].
• Using a geometric-based method for the AMG accelerator instead of a mathematic-based one [28].
As noted above, the code is developed in order to perform on unstructured and hybrid cells on a two or three dimensional numerical grid, for any free geometry. The algorithm of this code is represented in Fig. 3. The "C" programming language with double-precision variables was used for the CFD-code. Thresholds of the residuals for the convergence assessment were determined in accordance with the instructions of Ref. [14] for all of the relevant dependent variables.
3-2. Verification of the new CFD-based code
The full process of verification and validation of the new CFD code is beyond the scope of this study. Instead, the code was tested on 48 standard verification problems from Ref. [30]. The results of a few of these attempts are presented in Table 2. Figures (4), (5) and (6) are referred in the table to compare the results with the references. Other details can be found in the pertinent references introduced in Table 2. These examples should cover physical concepts of the thermo-hydraulic phenomena, which are relevant to the current study, such as turbulence simulation [31], natural convection [32], and conjugate heat transfer [33]. Note that the maximum error is less than two percent in all 48 cases. This should be reliable enough for most of the cases.
| Test case no. | Description | Ref. | Ref. result(s) | Current study result(s) | Error (maximum value) [%] | |||
|---|---|---|---|---|---|---|---|---|
| 1 | Turbulent heat transfer in a pipe expansion | [31] | Fig. 4(a) | Fig. 4(b) | 1.79 | |||
| 2 | Natural convection in a concentric annulus | [32] | Fig. 5(a) | Fig. 6(a) | Fig. 5(b) | Fig. 6(b) | 1.07 | 1.91 |
| 3 | Conjugate heat transfer in a composite solid block | [33] | The cooled wall | The Adiabatic wall | The cooled wall | The Adiabatic wall | The cooled wall | The Adiabatic wall |
| 378 [K] | 413 [K] | 378.150 [K] | 413.121 [K] | 0.04 | 0.03 | |||
3-3. Neutronic codes and algorithm of coupling them with the CFD-code
The CITVAP code from the MTR_PC package [34] is used for neutronic analysis. It solves one to three-dimensional multi-group diffusion equations in rectangular or cylindrical geometries. The code WIMSD-5B [35] is also employed to generate the macroscopic cross-sections required for these calculations. The structure of the energy groups used for the analysis is presented in Table 3.
| Energy group | 1 | 2 | 3 | 4 | 5 |
|---|---|---|---|---|---|
| Energy range | 10-0.821 [MeV] | 0.821-0.00553 [MeV] | 5530-0.625 [eV] | 0.625-0.08 [eV] | 0.08-0.00 [eV] |
The algorithm of coupling these codes with the CFD-code introduced in Sec. 3.1 is represented in Fig. 7. It should be noted that this algorithm is very straightforward. The output of each code is used to prepare the input of the next one and iteration continues until the convergence criteria are met. Similar to CFD-code (sec. 3-1) the convergence criteria are in accordance with the instructions in Ref. [14]. As noted in Sec. 2-3, by using the multi-region PMM, the same homogenized regions are used for both neutronic and thermo-hydraulic analyses. Nonetheless, as Fig. 7 suggests, data mapping (extrapolation) is anticipated in case the numerical grid is different for the two analyses. It should be noted that in the simulations of this study, a 64-bit PC platform with Intel’s core i7-980x (6 core) processor and 24 GB of RAM is used.
4. Brief description of Tehran Research Reactor (TRR)
TRR is a typical pool-type MTR. The relevant specifications of this reactor are presented in Table 4. A schematic diagram of TRR is illustrated in Fig. 8 and Fig. 9 demonstrates a symbolic view of the pool, including the core, gird plate, plenum chamber, and flapper valve. Under normal operational conditions, coolant flows downward through the core by the primary cooling circuit pump. However when shutdown or scram commences, the direction of flow reverses after a specific amount of time (hence it flows upward) and the residual heat would be removed by pure natural convection [36]. More information required for this study is presented in Sec. 5.1.
| General | |||
| Thermal power | 5 MW | ||
| Core dimensions (first operating core) | 40.5×38.54×89.7 cm | ||
| Grid plate locations for fuel assemblies | 9×6 | ||
| Each location cross-sectional dimensions | 7.71×8.1 cm | ||
| Reflector material | Graphite / Light water | ||
| Control elements | 4 Ag-In-Cd shim safety rods | ||
| 1 stainless steel regulating rod | |||
| Thermo-hydraulics | |||
| Coolant material | Light water | ||
| Cooling method | Operational condition | Forced flow primary loop | |
| Full shutdown condition | Natural circulation | ||
| Operating pressure | 0.171 Mpa | ||
| Primary coolant flow rate | 500 m3/h (2200 gpm) | ||
| Coolant inlet temperature (full power) | 37.8 ℃ (100 ℉) | ||
| Coolant outlet temperature (full power-reference core) | 46.1 ℃ (113.9 ℉) | ||
| Fuel assembly | |||
| Fuel | Plate type (MTR), LEU 20%, Al Clad | ||
| Fuel assembly dimensions (Max.) | 7.6×8.01×87.8 cm | ||
| No. of fuel plates in standard fuel assembly | 19 | ||
| No. of fuel plates in control fuel assembly | 14 | ||
| Coolant passage thickness | 0.27 cm | ||
| Fuel meat material | U3O8-Al | ||
| Fuel plate dimensions | 0.15×7.05×65.5 cm | ||
| Fuel meat dimensions | 0.07×6.00×61.5 cm | ||
5. The specific requirements of simulation
5-1. Description of the case studies
According to the SAR of TRR [36], when the shutdown commences, all of the control rods fall into the core. It causes an almost instantaneous drop in total power from the initial value of 5 MW to 400 KW, while the average temperature of the pool is 37.8 ℃. Afterwards, it takes 25 min for the holdup tank to be filled. At that moment, the coolant ceases to flow out of the pool and then the flapper valve opens, direction of flow reverses, and equilibrium condition of the pure natural convection are achieved. During this period, the total power drops to 100 KW and the average temperature of the pool reduces to 30 ℃. According to this information, two cases are selected to be studied. The conditions of these cases are presented in Table 5. Case 1 represents the actual conditions of pure natural circulation, while case 2 can be considered as the upper limit of the safety margin.
| Test cases | Case 1 | Case 2 |
|---|---|---|
| Total heat source [kW] | 100 | 400 |
| Average temperature of pool and the inlet temperature for evaluating the resistance coefficients [K] | 303.15 | 310.95 |
| The average inlet velocities for evaluating the resistance coefficients | Ten velocities in the range of 0.0 to 0.1 [m/s] (see table 8) | Ten velocities in the range of 0.0 to 0.1 [m/s] (see table 8) |
| Operational mode | Steady-state natural convection after shutdown | Steady-state natural convection after shutdown |
| Operating pressure [bar] | 1.71 | 1.71 |
| thermo-physical properties (both solids and fluid) | Same as table 6 | Same as table 6 |
5-2. Geometry and verification of the detailed model of a FA
The geometry of the solution domain for the detailed model of an FA is represented in Fig. 10. It includes a scattered view of essential consisting parts (Fig. 10(a)). The central part consists of fuel meat and clad plates as solid regions and passages between these plates and upper and lower plenum as fluid regions. The side part is a flow passage for cooling the outer side of the side plates. As Fig. 10 shows by four thin legs, this region continues to the end of the narrow passages that have been implemented in the grid plate for passing this cooling flow [36]. The top and bottom part are actually parts of the pool and plenum medium, respectively. These parts are included to simulate the effects of flow mixing and separation between the central and side parts and also for answering the question of how flow-rate divides between these two passages.
To verify this model and also to assess the importance of including all of these geometrical details, normal operational conditions of the reference core configuration (Fig. 11) are selected to test this model. Boundary conditions, thermo-physical properties, and the total mesh size that are used for this simulation are tabulated in Table 6. The results of the temperature rise and pressure drop are presented in Table 7. These results are sorted in a descending manner, from left to right, based on the numerousness of the error imposed from different geometrical and physical assumptions. The complete model is the result of discarding all of these assumptions. The results are compared to the counterpart data from the SAR of TRR [36], which are presented in the last column of Table 7 as the target results. It should be noted that every individual run of the detailed model on the system described in sec. 3.3 took almost three hours to reach the desired convergence state.
| inlet temperature (average temperature of coolant in the pool) | 310.95 [K] |
| Average inlet velocity of coolant | 0.7943 [m/s] |
| Operating Pressure (pressure at outlet) | 0.171 [Mpa] |
| Heat source | Cosine function of axial direction with the average value of the reference core configuration (see fig. 10) in full power (5MW) |
| Coolant thermo-physical properties(ρ, µ, k, Cp) | Temperature-dependant polynomials that are fitted to the data extracted from the steam table at the operating pressure. |
| Fuel and clad thermo-physical properties(ρ, k, C) | From ref. [37]. (thermal conductivity of fuel should be defined as a function of temperature) |
| Total mesh size | Approximately 4.6 million unstructured hexahedral cells |
| Applied assumptions and simplifications | Discarding the side, bottom and top parts from fig. 10 | Using the simple turbulence model (see table 1) | Discarding passages at ends of fuel plates from fig. 10 | Using constant thermo-physical properties | Complete model | Target results [36] |
|---|---|---|---|---|---|---|
| Δp [kpa] | 6.017 | 9.408 | 6.889 | 7.286 | 7.981 | ≈ 8 |
| Error [%] (relative to the target results) | 24.79 | 17.60 | 13.89 | 8.93 | 0.24 | ----- |
| ΔT [K] | 9.392 | 8.403 | 9.372 | 8.204 | 8.386 | 8.3 |
| Error [%] (relative to the target results | 13.16 | 1.24 | 12.92 | 1.16 | 1.04 | ----- |
These results show that the largest amount of error occurs when the effects of flow mixing and separation are ignored by omitting the side, top, and bottom parts. From a total flow rate of 8.521 kg/s that enters domain, share of side part is 0.0107 kg/s which means 0.126%. This amount is 27% lower than 0.0172 [kg/s] that is obtained from the fraction of cross sectional area of these two passages. Therefore, pressure loss and temperature rise in the side part are always higher than the central part. For this reason, passage of the side part is the hot channel. This matter can be observed more clearly in Fig. 12, in which counters of temperature distribution on some virtual cross-sections, vertical to the axis of the fuel assembly, are shown. The second recognized source of error is applying the simple turbulence model instead of the advanced one (see Table 1). The third source of error roots in disregarding 38 very narrow passages that exist at the both end of each plate. Cross sectional dimension of these passages is only 0.75×1.5 mm2. Comparing the results of the fifth column of Table 7 with the data of the last column prove the importance of including temperature-dependent material properties and, according to the errors reported in Table 7, using constant thermo-physical properties is recognized as the fourth source of error.
5-3. Geometry and conditions for simulation of natural convection in the whole pool of reactor
Geometry of the solution domain for the simulation of natural convection in the pool of TRR is shown in Fig. 13. It includes the core (with the configuration of the reference core (Fig .10)), plenum chamber, flapper valve, and the guide tube (see Fig. 9). This geometry took almost 11 million unstructured hexahedral meshes. The simulation is performed in two modes. First, when the whole core is simplified to a single porous medium (Fig. 13(a)) same as Ref. [12-13] and second, when the core is divided into several media that each comprise an FA (Fig. 13(b)), as the methodology of this study (Sec. 2) suggests.
6. A semi-analytical method for evaluating the equilibrium velocity of natural convection
The equilibrium velocity of pure natural convection in the steady-state condition can be evaluated by the following equation [38].
In the case of a pool-type reactor, each side can be estimated as follows.
The convention of direction for the loop integration of Eq. (13) is presented in Fig. 14. If it is possible to state pressure drop in Eq. (12) and average density (ρa) in Eq. (13), as a function of the equilibrium velocity (vE), this velocity can be easily calculated by solving Eq.(11). For this purpose, a detailed simulation of an FA (Fig. 10) is used. This model is subjected to different reverse (upward) velocities for the two cases defined in Sec. 5.1. The results are tabulated in Table 8. Then a polynomial curve is fitted to the data of the pressure drop and average density, as Fig. 15 serves as an example for case 1. If these polynomial functions are used in Eq.(11), the unknown velocity will be obtained. The results are 0.01613 and 0.02408 m/s for cases 1 and 2, respectively.
| Va[m/s] | 0.01 | 0.02 | 0.03 | 0.04 | 0.05 | 0.06 | 0.07 | 0.08 | 0.09 | 0.1 | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Case 1 | |||||||||||
| ΔP [pa] | 17.376 | 34.318 | 56.927 | 81.233 | 104.720 | 131.001 | 160.920 | 194.000 | 229.507 | 267.070 | |
| ΔT [K] | 13.574 | 6.677 | 4.482 | 3.363 | 2.691 | 2.237 | 1.913 | 1.672 | 1.485 | 1.336 | |
| ρa[kg/m3] | 993.169 | 994.549 | 994.997 | 995.226 | 995.361 | 995.449 | 995.512 | 995.558 | 995.594 | 995.622 | |
| μa×104[pa.s] | 5.9848 | 6.3575 | 6.4986 | 6.5736 | 6.6190 | 6.6495 | 6.6715 | 6.6880 | 6.7009 | 6.7113 | |
| Case 2 | |||||||||||
| ΔP [pa] | 27.868 | 35.347 | 49.198 | 67.455 | 89.812 | 116.134 | 146.947 | 181.082 | 219.087 | 261.075 | |
| ΔT [K] | 54.938 | 27.627 | 18.535 | 13.891 | 11.091 | 9.185 | 7.849 | 6.866 | 6.108 | 5.508 | |
| ρa [kg/m3] | 975.120 | 986.419 | 988.647 | 989.820 | 990.524 | 990.981 | 991.301 | 991.535 | 991.714 | 991.857 | |
| μa×104 [pa.s] | 4.2222 | 5.3462 | 5.7110 | 5.9387 | 6.0932 | 6.2022 | 6.2828 | 6.3444 | 6.3929 | 6.4325 | |
7. Results and discussions
The process of presenting the results starts with the steps of the procedure for evaluating resistance coefficients, explained in Sec. 2.1, for the two cases of this study, introduced in Sec. 5.1. Table 8 presents temperature and pressure variations and also the average density and viscosity of the coolant in the detailed model of an FA (see Sec. 5.2) when this model is subjected to ten velocities in the range of 0.01 to 0.1 m/s. An example of fitting a quadratic function to the data of pressure loss against velocity is shown in Fig. 16 for the case 1. The resulting correlation is compared with Eq. (3), whereby coefficients α and β are obtained. The results of calculating resistance coefficients using this data at each velocity via Eqs. (4) and (5) are tabulated in table 9 for the both cases. An example of fitting polynomial functions through this data and deriving the final correlations of the resistance coefficients is presented in Fig. 17 for the case 1
| Inlet velocity [m/s] | Case 1 | Case 2 | ||
|---|---|---|---|---|
| Cz [1/m] | Dz [1/m2] | Cz [1/m] | Dz [1/m2] | |
| 0.01 | 26.362 | 2435764.801 | 27.211 | 3059502.594 |
| 0.02 | 26.325 | 2283185.812 | 27.004 | 2559249.639 |
| 0.03 | 26.313 | 2230476.945 | 26.947 | 2377489.594 |
| 0.04 | 26.307 | 2202839.890 | 26.921 | 2283411.092 |
| 0.05 | 26.304 | 2186366.416 | 26.905 | 2226685.155 |
| 0.06 | 26.302 | 2175499.402 | 26.895 | 2188836.942 |
| 0.07 | 26.300 | 2167776.079 | 26.888 | 2161466.842 |
| 0.08 | 26.299 | 2162001.777 | 26.883 | 2141048.275 |
| 0.09 | 26.298 | 2157520.607 | 26.879 | 2124844.340 |
| 0.10 | 26.297 | 2153938.994 | 26.876 | 2111836.963 |
To verify the accuracy of these correlations, a porous medium in the shape of a box with the external dimensions of an FA is modeled. The correlations are applied to this medium and it subjected to the same velocities and conditions. The results of the temperature and pressure variations along with their counterparts from the detailed simulation are presented in Table 10 for the case 1. Note that the maximum error is less 0.5 percent in all cases. As noted in Sec. 2.3, such a high level of accuracy is transferred to the simulation of the whole reactor pool (or any similar large-scale simulation) because these media are placed into the core without any changes.
| Va [m/s] | 0.01 | 0.02 | 0.03 | 0.04 | 0.05 | 0.06 | 0.07 | 0.08 | 0.09 | 0.1 | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| ΔT [K] | |||||||||||
| Detailed model | 13.574 | 6.677 | 4.482 | 3.363 | 2.691 | 2.237 | 1.913 | 1.672 | 1.485 | 1.336 | |
| Replaced porous medium | 13.573 | 6.679 | 4.480 | 3.366 | 2.689 | 2.239 | 1.912 | 1.672 | 1.484 | 1.336 | |
| ΔP [pa] | |||||||||||
| Detailed model | 17.376 | 34.318 | 56.927 | 81.233 | 104.720 | 131.001 | 160.920 | 194.000 | 229.507 | 267.070 | |
| Replaced porous medium | 17.058 | 34.540 | 56.432 | 80.526 | 104.373 | 131.860 | 161.076 | 194.417 | 229.075 | 266.059 | |
Now, the results of simulating natural convection in the whole pool of TRR are presented. In this simulation, the instructions explained in Sec 2.4 by the procedures depicted in Fig. 2 and the algorithm of coupling neutronic and thermo-hydraulic analyses represented in Fig. 7 are followed. It should be noted that every individual run of the whole pool model on the system described in Sec. 3.3 took almost 17 hours to reach the desired convergence state. The average power production in all the FAs and the axial shape of the power distribution in the hot FA are demonstrated in Fig. 18. These are the final results of the neutronic analysis that are transferred to the thermo-hydraulic one for the case 1. Table 11 presents temperature and pressure variations and the average velocity in all FAs from the results of the thermo-hydraulic analysis. A general pattern of all thermo-hydraulic output variations are consistent with the distribution of the average power presented in Fig. 18(a). Furthermore, Average velocities of the coolant in the whole core, tabulated in the last row of Table 11, for the both cases, are in a very good agreement with their counterpart analytical results, presented in Sec. 6. These results are directly compared in the second row of Table 12 for case 1. Note that the calculated error is only 0.31 percent.
| Position of FA in the core | Case 1 | Case 2 | ||||
|---|---|---|---|---|---|---|
| ΔP [pa] | ΔT [K] | Vz [m/s] | ΔP [pa] | ΔT [K] | Vz [m/s] | |
| A3 | 25.899 | 7.904 | 0.01502 | 37.811 | 21.278 | 0.02256 |
| A4 | 27.463 | 8.381 | 0.01593 | 40.094 | 22.563 | 0.02392 |
| A5 | 26.570 | 8.109 | 0.01541 | 38.791 | 21.829 | 0.02314 |
| B2 | 25.835 | 7.884 | 0.01498 | 37.718 | 21.226 | 0.02250 |
| B3 | 27.213 | 8.305 | 0.01578 | 39.730 | 22.358 | 0.02370 |
| B4 | 27.676 | 8.446 | 0.01605 | 40.406 | 22.738 | 0.02410 |
| B5 | 26.678 | 8.142 | 0.01547 | 38.948 | 21.918 | 0.02323 |
| B6 | 26.124 | 7.972 | 0.01515 | 38.139 | 21.462 | 0.02275 |
| C2 | 26.867 | 8.199 | 0.01558 | 39.224 | 22.073 | 0.02340 |
| C3 | 30.419 | 9.283 | 0.01764 | 44.410 | 24.991 | 0.02649 |
| C4 | 26.761 | 8.167 | 0.01552 | 39.070 | 21.987 | 0.02331 |
| C5 | 30.896 | 9.429 | 0.01792 | 45.106 | 25.383 | 0.02691 |
| C6 | 27.177 | 8.294 | 0.01576 | 39.677 | 22.328 | 0.02367 |
| D2 | 27.757 | 8.471 | 0.01610 | 40.523 | 22.804 | 0.02417 |
| D3 | 26.462 | 8.076 | 0.01535 | 38.633 | 21.741 | 0.02305 |
| D4 | 34.218 | 10.443 | 0.01984 | 49.956 | 28.113 | 0.02980 |
| D5 | 26.947 | 8.224 | 0.01563 | 39.341 | 22.139 | 0.02347 |
| D6 | 28.031 | 8.555 | 0.01626 | 40.924 | 23.030 | 0.02441 |
| E2 | 27.671 | 8.445 | 0.01605 | 40.398 | 22.734 | 0.02410 |
| E3 | 32.224 | 9.834 | 0.01869 | 47.045 | 26.475 | 0.02806 |
| E4 | 27.122 | 8.277 | 0.01573 | 39.596 | 22.283 | 0.02362 |
| E5 | 32.463 | 9.907 | 0.01882 | 47.393 | 26.670 | 0.02827 |
| E6 | 27.120 | 8.301 | 0.01577 | 39.710 | 22.347 | 0.02369 |
| F2 | 25.558 | 7.800 | 0.01482 | 37.313 | 20.998 | 0.02226 |
| F3 | 25.744 | 7.857 | 0.01493 | 37.584 | 21.150 | 0.02242 |
| F4 | 27.477 | 8.385 | 0.01593 | 40.114 | 22.574 | 0.02393 |
| F5 | 26.124 | 7.972 | 0.01515 | 38.139 | 21.463 | 0.02275 |
| F6 | 25.846 | 7.888 | 0.01499 | 37.734 | 21.235 | 0.02251 |
| Arithmetic average | 27.729 | 8.462 | 0.01608 | 40.483 | 22.782 | 0.02415 |
| Conditions | Following the instructions of ref. [12-13] except for the simplifications of the detailed model. | Using correlations of resistance coefficients evaluated in this study (see fig. 17). | Averaged results of multi-region PMM from last row of table 11 | Semi-analytical result from sec. 6 |
|---|---|---|---|---|
| Vz [m/s] | 0.01364 | 0.01467 | 0.01608 | 0.01613 |
| Error [%] (relative to last column) | 15.44 | 9.05 | 0.31 | ---- |
| Δp [pa] | 23.525 | 25.240 | 27.729 | ---- |
| Error [%] (relative to fourth column) | 15.16 | 8.40 | ----- | ---- |
| ΔT [°C] | 9.711 | 9.213 | 8.462 | ---- |
| Error [%] (relative to fourth column) | 14.76 | 8.87 | ----- | ---- |
Contours of temperature distribution on the three vertical cross sections that cut through the center of the core are represented in Fig. 19 for the case 1. It describes how a current of natural convection transfers the residual heat form the core to the pool and dissipates this heat by mixing with the relatively cooler flow in the pool. Stream-lines in the pool of TRR colored by the velocity magnitude of flow are shown in Fig. 20. These lines are obtained by tracking some inert and massless particles released from the surface of the flapper valve (see Fig. 9). These type of results can be used to improve the quality of current and future designs. For example, it can be used to minimize pressure losses by streamlining a technique. It can also be used to achieve a better physical understanding of the general pattern of flow in a large-scale simulation. For example, for any pool-type reactor, avoiding direct flow of coolant, which just passed through the core to the free surface of the pool, is an obvious goal because it can contaminate the workspace under the containment by increasing the level of activity. Paths of stream-lines in fig. 20 prove that the design of the TRR is adequate for this purpose.
The results of the simulation using the conventional single-region PMM are presented in Table 12. This simulation is conducted using the geometry depicted in Fig. 13(a) and in the following two modes.
1. When all instructions in Ref. [12-13] are followed to the letter, excluding simplifications of the detailed model of an FA, because they clearly had a devastating effect on the accuracy of simulation (25 percent of error in initial evaluation of pressure losses is reported by Ref. [12]).
2. When the correlations of resistance coefficients calculated in this study (see Fig. 17) are applied to the single porous medium.
The results of simulation for these two modes are also compared in Table 12 with their counterparts from the multi-region PMM (the last row of table 11) and the analytical results evaluated in Sec. 4. Results presented in this table demonstrate that by applying the correlations, the accuracy of the estimation of the single-region model improves significantly but the errors are still considerable. Overall, using single-region PMM, regardless of the selected technique, would lead to losing all details inside the core and impacts of the core flow on the surrounding regions would suffer from significant inaccuracy due to unrealistic simplifications described in Sec. 2. These are the reasons for the considerable amount of errors tabulated in Table 12 for the both modes of single-region PMM and even more serious errors reported in Ref. [12] and [13].
Another noteworthy result in regards to the simulation of natural circulation is represented in Fig. 21. As noted in Sec. 2.1, subjecting the detailed model of FA to different velocities in order to estimate the relationship between the pressure (force) and velocity fields is the main step in the process of evaluating the resistance coefficients (see Fig. 1). Figure 21 demonstrates that by reducing the average inlet velocity, the position where the minimum pressure occurs gradually moves from the top of the fuel plates toward the center of it. Meanwhile, the maximum temperature occurs around the center of plates and this position does not reveal any significant movement by variation in velocity. The intriguing fact is that when velocity approaches the values pertaining to natural convection equilibrium (i.e. 0.0161 and 0.0241 m/s for the cases 1 and 2 respectively) these two positions (i.e. positions where maximum temperature and minimum pressure occur) coincide. This position also becomes the place where coolant velocity has its maximum value. Therefore, in the view of the whole circulation path (whole pool of the reactor), it is as if a hypothetical pump is located at this position which works as a "driving force" for natural convection. Furthermore, this process can be used as an alternative method for estimating equilibrium velocity of natural circulation, instead of modeling a whole reactor pool or semi-analytical method presented in sec. 6. Comparing Figs. 21(b) and (c) demonstrates that after reaching the equilibrium velocity, the position where the minimum pressure occurs does not move any further by reducing velocity but the minimum pressure becomes a larger negative value rapidly. Moreover, if the average inlet velocity is larger than the values mentioned in Fig. 21, the position where minimum pressure occurs (also known as the point of pressure gradient reversal [39]) will never appear on the fuel plates. This means that in these situations, momentum of the coolant due to forced convection is large enough to overcome all drag forces (including the effects caused by natural convection) and avoid separation of flow [39].
The result of calculating some well-known dimensionless numbers regarding natural convection in different velocities are presented in Table 13. The definition of the dimensionless numbers can be found in Ref. [16]. These numbers measure importance and strength of natural convection [16]. It should be noted that all quantities used for these calculations are a volume-weighted average of the relevant material properties, such as density, viscosity, and thermal conductivity on the volume of the coolant in the all passages between fuel plates. These results demonstrate that when velocity approaches the value corresponding to the natural convection equilibrium from larger values, a rapid growth in all of the dimensionless numbers can be detected which indicates a significant increase in the strength of the buoyancy-induced flow [16]. Figure 22 represents the variation of Rayleigh (Ra) number against average velocity for both cases in this study. Note that equilibrium velocities of natural circulation are also indicated for both cases on this diagram for comparison. These results indicate that monitoring the variation of the dimensionless numbers against velocity can be used as another alternative method for estimating equilibrium velocity of natural circulation in a qualitative manner. Although it should be obvious that this method only specifies an approximate range for equilibrium velocity and it is not as accurate as the former methods proposed for this purpose.
| Inlet velocity [m/s] | Case 1 | Case 2 | ||||
|---|---|---|---|---|---|---|
| Gr | Gr/Re2 | Ra | Gr | Gr/Re2 | Ra | |
| 0.01 | 2.433E+10 | 233.58 | 1.126E+11 | 2.240E+11 | 1651.60 | 7.507E+11 |
| 0.02 | 1.004E+10 | 27.12 | 4.998E+10 | 6.593E+10 | 171.07 | 2.692E+11 |
| 0.03 | 6.328E+09 | 7.93 | 3.234E+10 | 3.488E+10 | 46.41 | 1.547E+11 |
| 0.04 | 4.593E+09 | 3.31 | 2.381E+10 | 2.280E+10 | 18.46 | 1.058E+11 |
| 0.05 | 3.601E+09 | 1.68 | 1.882E+10 | 1.713E+10 | 9.33 | 8.178E+10 |
| 0.06 | 2.954E+09 | 0.97 | 1.553E+10 | 1.369E+10 | 5.35 | 6.664E+10 |
| 0.07 | 2.503E+09 | 0.61 | 1.321E+10 | 1.140E+10 | 3.36 | 5.628E+10 |
| 0.08 | 2.172E+09 | 0.41 | 1.150E+10 | 9.766E+09 | 2.24 | 4.873E+10 |
| 0.09 | 1.919E+09 | 0.28 | 1.018E+10 | 8.530E+09 | 1.57 | 4.293E+10 |
| 0.10 | 1.718E+09 | 0.21 | 9.136E+09 | 7.565E+09 | 1.14 | 3.834E+10 |
In the end, to practically demonstrate how the detailed results are retrieved by the multi-region PMM of this study, the values of all relevant variables on the inlet and outlet boundaries of the hot FA of case 2 from the simulation of the whole pool (located at position D4 from Fig. 11) are mapped onto their counterparts in the detailed model of a single FA. Note that the power source and the pertinent shape of distribution should be taken from neutronic part of the simulation (see Fig. 18). Contours of temperature and pressure distribution on the outer surface of this FA (the hot FA), that include the walls of the hot channel, are represented in Fig. 23. The results show that even in the intense conditions of case 2, the pure natural convection still has the capability of removing the residual heat safely, because there is a significant difference between the maximum temperature of the coolant in the hot FA (354.360 K) and the saturation temperature corresponding to the operating pressure of the system (388.478 K).
8. Conclusion
As noted in Sec. 2.3, the multi-region PMM introduced in this study was inspired by the conventional procedures for neutronic core calculations. In other words, the new PMM is specifically developed to be a thermo-hydraulic analog of the neutronic procedure in order to achieve an equivalent neutronic and thermo-hydraulic solution method. This method facilitates the use of the CFD approach in the coupled neutronic-thermohydraulic simulation of a large-scale system for an average researcher because it reduces computational cost significantly. However, the main feature of the modified method, as explained in Sec. 2-3 (see Fig. 2(b)) and practically demonstrated in Sec. 7, is the capability of retrieving all the detailed results that seem to be lost by using the PMM. In this study, the multi-region PMM was applied to simulate natural circulation in a typical MTR because it has all the requirements to test and verify the new model. Nonetheless, it should be clear that the multi-region PMM can be used for a wide range of applications in PWR, BWR, WWER, etc. This matter will be the subject of our future studies.
Fluid Flow through Packed Columns
. Chem. Eng. Prog, 1952, 48: 89-94.The role of porous media in modeling flow and heat transfer in biological tissues
. Int. J. Heat Mass Transfer, 2003, 46: 4989-5003. DOI: 10.1016/S0017-9310(03)00301-6.Thermal dispersion in a porous medium
. Int. J. Heat Mass Transfer, 1990, 33(8): 1587-1597. DOI: 10.1016/0017-9310(90)90015-M.A general two-equation macroscopic turbulence model for incompressible flow in porous media
. Int. J. Heat Mass Transfer, 1997, 40(13): 3013-3024. DOI: 10.1016/S0017-9310(96)00370-5.Detailed study of a model of heat and mass transfer during convective drying of porous media
. Int. J. Heat Mass Transfer, 1988, 31(5): 957-967. DOI: 10.1016/0017-9310(88)90084-1.Development of a Coupled CFD System-code Capability & its Applications to Simulate Current & Next Generation Reactors. PhD thesis
.CFD Simulation of a Research Reactor
. InNumerical solution of Navier-Stokes equations
". Mathematics of Computation, 1968, 22: 745-762. DOI: 10.1090/S0025-5718-1968-0242392-2.A New k-ε Eddy-Viscosity Model for High Reynolds Number Turbulent Flows - Model Development and Validation
. Computers Fluids, 1995, 24: 227-238. DOI: 10.1016/0045-7930(94)00032-T.Near-Wall Turbulence Models for Complex Flows Including Separation
. AIAA Journal, 1988, 26: 641-648. DOI: 10.2514/3.9948.The design and application of upwind schemes on unstructured meshes. Technical Report AIAA-89-0366
.The ULTIMATE conservative difference scheme applied to unsteady one-dimensional advection
. Comp. Methods Appl. Mech. Eng., 1991, 88: 17-74. DOI: 10.1016/0045-7825(91)90232-U.An Implicit Upwind Algorithm for Computing Turbulent Flows on Unstructured Grids
. Computers Fluids, 1994, 23(1): 1-21. DOI: 10.1016/0045-7930(94)90023-X.A Fast Nested Multi-grid Viscous Flow Solver for Adaptive Cartesian/Quad Grids
. International Journal for Numerical Methods in Fluids, 2000, 33: 657-680. DOI: 10.1002/1097-0363(20000715)33:5<657::AID-FLD24>3.0.CO;2-G.Numerical Study of the Turbulent Flow past an Airfoil with Trailing Edge Separation
. AIAA Journal, November 1983, 21(11): 1525-1532. DOI: 10.2514/3.8284.Solution of Implicitly Discretized Fluid Flow Equations by Operator Splitting
. Journal of Computational Physics, 1986, 62: 40-65. DOI: 10.1016/0021-9991(86)90099-9.A Multi-grid Method Based on the Additive Correction Strategy
. Numerical Heat Transfer, 1986, 9: 511-537. DOI: 10.1080/10407788608913491.Local Heat Transfer Downstream of an Abrupt Expansion in a Circular Channel with Constant Wall Heat Flux
, Journal of Heat Transfer, 1984, 106: 789-796. DOI: 10.1115/1.3246753.An Experimental Study of Natural Convection Heat Transfer in Concentric and Eccentric Horizontal Cylindrical Annuli
. Journal of Heat Transfer, 1978, 100: 635-640. DOI: 10.1115/1.3450869Safety Analysis Report for the Tehran Research Reactor (LEU)
.