Optimization of Natural Circulation District Heating Reactor Primary Heat Exchangers

: Small modular reactors (SMRs) are gaining interest as a potential solution for cost-effective, carbon-neutral district heat (DH) production. The low pressures and temperatures permit much lighter and cheaper designs than in power plants, and efﬁciency is high as all heat generated can be sold to customers. In this work, the optimization of the primary heat exchangers in a natural-circulation 50-MW heating reactor concept was carried out to obtain an initial feasibility estimate for the concept for both baseload and load-following operation, as well as to obtain information on the characteristics of an optimized design. Studies on small natural circulation heat-only SMRs and the impact of heat exchanger design on the overall dimensions and economics have not been published before. Although a detailed heat exchanger cost model was used, the results should be considered tentative initial estimates, as much of the cost impact from the heat exchanger design comes from the effect the design has on the pressure vessel dimensions. While more detailed pressure vessel designs and cost functions are needed for ﬁnal optimization, the feasibility of the concept is shown. Optimization for different load proﬁles produced near-identical designs, with the downcomer divided approximately in half between the heat exchanger at the top and an empty space at the bottom to maximize the pressure difference available for natural circulation. Although conservative, even pessimistic estimates were used in the absence of detailed cost functions, cost prices of 30–55 EUR/MWh DH at a 10% interest rate were obtained, or only 20–40 EUR/MWh DH at a 5% interest rate. This indicates potentially good competitiveness for the considered DH SMR concept.


Background
Climate change is underway due to anthropogenic greenhouse gas emissions that are mainly caused by energy use: electricity, space heating in cold climates, and traffic fuels. In Finland, greenhouse gas emissions from power generation are already low, consistently below 70 kg CO2 /MWh [1], due to the large shares of nuclear, hydro, and biomass-based generation and the recently increasing amount of wind power. Space heating, however, is still 1/3 fossil-fueled [2]. District heating (DH) networks in the country are well developed, with a heating market share of 46% and the majority of the population within reach [2]. With a favorable public attitude toward nuclear energy [3], nuclear district heating is a potential solution for decarbonizing the remaining heating solution. In this work, the feasibility of a natural-circulation heat-only nuclear reactor is carried out by means of techno-economic optimization of the primary heat exchangers for such a concept.

District Heating Nuclear Reactors
A nuclear fission reactor liberates the nuclear binding energy of heavy elements. Macroscopically, nuclear energy manifests as heat released into the reactor fuel. Nuclear

Goals
The main goal of this study was to obtain an initial estimate for the feasibility of a small natural circulation heating reactor design for baseload DH production in a large city as well as in a smaller DH network where much of the operation would be at partial loads. The results show that the impact of load profile assumption on the optimized heat exchanger configuration is negligible, and separate designs for constant baseload and varying-load uses are not needed. From the results of this study, information is also obtained on the broad outlines of likely optimal heat exchanger design features for such a reactor concept, an approximate cost price for the generated heat, and suitability and appropriate tuning parameters of the cuckoo search algorithm for the considered optimization task. Studies on small natural circulation heat-only SMRs and the impact of primary heat exchanger design on the overall dimensions and plant economic performance have not been published before.

District Heat Demand and Temperature Levels
The primary application for the considered reactor concept is providing district heating (DH) in a Nordic context. In Finland and Scandinavia, district heating dominates the urban heating markets. The implementation is almost always a two-pipe hot-water system, with supply and return temperatures typically varying between 70-120 • C, and 35-55 • C, respectively. The temperature levels vary seasonally and also between networks, depending on the climate and the types of heat consumers in the network.
The seasonal variation of the heat load is significant: in Finland, the winter peak consumption typically exceeds the summer minimum by a factor of ten. A discretized approximation of the annual temperature and load profiles is depicted in Figure 1 as duration curves. The DH load profile is based on [23]; ambient temperatures represent 10-year averages measured in Jyväskylä, a mid-size city in central Finland [24], and the DH supply temperature represents the recommended temperature curve as a function of ambient temperature in Finnish DH networks [25]. The return temperature represents average values measured in an operating network at a given supply temperature level [26].
Two scenarios are considered, baseload (BL) and mid-load (ML). In the BL scenario, the SMR operates at full thermal power (but varying temperature levels, Figure 1) for 8200 h annually. In the mid-load scenario ML, the DH SMR system is considered in a situation where the full load corresponds to 35% of the peak load of the DH system. Table 1 summarizes each load point with which the annual operating profile is approximated. The case BL was considered in two variants, the basic 50 MW and an additional 200-MW variant for a larger city. In the second case, BL200 four reactors are placed in line with a one-meter separation of containment vessels in an elongated rectangular pool cavity. Two scenarios are considered, baseload (BL) and mid-load (ML). In the BL scenario, the SMR operates at full thermal power (but varying temperature levels, Figure 1) for 8200 h annually. In the mid-load scenario ML, the DH SMR system is considered in a situation where the full load corresponds to 35% of the peak load of the DH system. Table 1 summarizes each load point with which the annual operating profile is approximated. The case BL was considered in two variants, the basic 50 MW and an additional 200-MW variant for a larger city. In the second case, BL200 four reactors are placed in line with a one-meter separation of containment vessels in an elongated rectangular pool cavity.

Studied District Heating SMR Concept
The considered pressurized water DH reactor is based on the design proposed in [14,15]. A notable feature is that the reactor pressure vessel (RPV) is enclosed inside another pressure vessel, the containment vessel (CV), placed in a water pool. Figure 2 depicts the schematic diagram of the considered SMR concept. An intermediate secondary circuit separates the reactor's primary circuit from the DH network, connecting the two via heat exchangers. In the reactor pressure vessel, the heated primary circuit water flows up through the cylindrical riser, returning through the heat exchangers located in the annular downcomer. The primary circuit heat exchangers are 1:1 counterflow segmentedbaffle shell-and-tube heat exchangers with the hot primary water in the tubes and secondary circuit water in the shell side. The construction is essentially one-pass welded TEMA E but without tube-side entry and exit headers.

Studied District Heating SMR Concept
The considered pressurized water DH reactor is based on the design proposed in [14,15]. A notable feature is that the reactor pressure vessel (RPV) is enclosed inside another pressure vessel, the containment vessel (CV), placed in a water pool. Figure 2 depicts the schematic diagram of the considered SMR concept. An intermediate secondary circuit separates the reactor's primary circuit from the DH network, connecting the two via heat exchangers. In the reactor pressure vessel, the heated primary circuit water flows up through the cylindrical riser, returning through the heat exchangers located in the annular downcomer. The primary circuit heat exchangers are 1:1 counterflow segmented-baffle shell-and-tube heat exchangers with the hot primary water in the tubes and secondary circuit water in the shell side. The construction is essentially one-pass welded TEMA E but without tube-side entry and exit headers. The secondary circuit transfers heat to the DH network via plate heat exchangers. To lower the DH supply temperatures during warmer seasons while maintaining stable primary circuit conditions, a shunt connection, Figure 2, allows partial bypassing of the secondary heat exchangers.  The secondary circuit transfers heat to the DH network via plate heat exchangers. To lower the DH supply temperatures during warmer seasons while maintaining stable primary circuit conditions, a shunt connection, Figure 2, allows partial bypassing of the secondary heat exchangers.

Objective Function
The design of the heat exchanger affects not only the cost of the heat exchangers themselves but also the dimensions and, thereby, the cost of the reactor pressure vessel (RPV), containment vessel (CV), and, to some extent, the reactor pool. In the ML scenario, partial-load operation and full-load hours may also be affected. To consider these effects, annual net cash flow C net maximization is considered as the objective function,  Table 2, determining the pumping electricity consumption in Equation (1) is described in Section 2.4, and C TCI is calculated by incrementing the base plant costs C SMR with the free-on-board costs of the heat exchangers (HX), RPV, and CV, and the cost of the reactor pool cavity for the CV.
where the factor 3.3 [27] is used to convert the free-on-board (FOB) costs of heat exchangers and pressure vessels to TCI costs. The reactor pool cavity cost C cav represents the cost of excavating and building an extension downward from the 7-m-deep reactor pool (see Figure 3). The main pool is of fixed dimensions and cost, regardless of heat exchanger sizing, and is considered part of the C SMR . The cavity dimensions are affected by the heat exchanger sizing, however, and, thus, are evaluated as a function of the resulting dimensions. It should be noted that the C SMR figure is tentative, based on the estimated cost of a reactor similar to the Chinese DHR-400 DH SMR built in Finnish conditions [28], scaled with a capacity scaling exponent of 0.7. Although this includes the RPV and heat exchangers, for a conservative estimate, the figure was not adjusted down. While the estimate has considerable uncertainty, as there are no small DH SMRs in Europe yet, other than the costs related to RPV, CV, and excavation, the plant costs are unlikely to be substantially affected by the heat exchanger solutions. For the purposes of a natural circulation feasibility evaluation and an initial heat exchanger optimization, a detailed plant investment evaluation was considered unnecessary and beyond the scope. To find conservative upper-limit cost figures, the base C SMR cost was not adjusted for the RPV, CV, heat exchanger, and partial excavation costs. These costs, as obtained in this study, were simply added to the total scaled C SMR figure.  It should be noted that the CSMR figure is tentative, based on the estimated cost of a reactor similar to the Chinese DHR-400 DH SMR built in Finnish conditions [28], scaled with a capacity scaling exponent of 0.7. Although this includes the RPV and heat exchangers, for a conservative estimate, the figure was not adjusted down. While the estimate has considerable uncertainty, as there are no small DH SMRs in Europe yet, other than the costs related to RPV, CV, and excavation, the plant costs are unlikely to be substantially affected by the heat exchanger solutions. For the purposes of a natural circulation feasibility evaluation and an initial heat exchanger optimization, a detailed plant investment evaluation was considered unnecessary and beyond the scope. To find conservative upper-limit cost figures, the base CSMR cost was not adjusted for the RPV, CV, heat exchanger, and partial excavation costs. These costs, as obtained in this study, were simply added to the total scaled CSMR figure.
The cost estimates of the heat exchanger and pressure vessels are based on mechanical sizing, performed similarly to [29]. Precise design to any one current pressure vessel code was considered unnecessary for the purposes of this study; instead, the 1996 Finnish pressure vessel code [30] is used to estimate the dimensions and masses of the main parts of the heat exchangers and the RPV and CV. ASTM 316L stainless steel The cost estimates of the heat exchanger and pressure vessels are based on mechanical sizing, performed similarly to [29]. Precise design to any one current pressure vessel code was considered unnecessary for the purposes of this study; instead, the 1996 Finnish pressure vessel code [30] is used to estimate the dimensions and masses of the main parts of the heat exchangers and the RPV and CV. ASTM 316L stainless steel (X2CrNiMo-17-13-3) was used as material for all said components. The fixed parameter values used in mechanical and thermal calculations are summarized in Table 3 below. The sizing is based on 150 • C/5 bar maximum operating conditions in the RPV and a safety factor of 2 for the pressure, and 4 bar(g) maximum pressure in the secondary circuit. Design values of −1.0/+0.4 MPa(g) for shell pressure and 180 • C were thus considered for the heat exchangers. For the RPV, the design pressures were −0.4/+1.0 MPa(g), and for the CV −0.25/+0.4 MPa(g). The 0.4 MPa for the RPV outside and CV internal pressure is set by the possibility of leakage in secondary water pipes pressurizing the containment vessel, and the 0.25 MPa for the maximum hydrostatic pressure at the CV bottom depth. The sizing procedure is described in more detail in Appendix A.

Heat Exchanger Cost Model
The heat exchanger manufacturing cost C man,HX is estimated via a simplified version of the model of Caputo et al. [32]. The material and processing costs are evaluated separately, and the equipment FOB cost C FOB is obtained by adding the overhead costs, contingency, and manufacturer's profit, estimated at 30%, 5%, and 10% of the FOB, respectively, and a value-added tax of 24%. The total capital investment of the 16 heat exchangers is then determined using a factor of 3.3, according to Sinnott [27]: C TCI,HX = 16 · 3.3C FOB,HX .
The material costs are dominant in the manufacturing cost. From the mechanical sizing results, the component volumes are obtained and, thereby, the masses. The basic price of 316 L steel is assumed at 6.09 €/kg for all parts except the tubes; seamless drawn tubing is considerably more expensive at 29 €/kg. Material is obtained in the shape of rectangular plates, pipes, tubes, and rods. During processing, some material becomes scrap; a scrap metal price of c sh = 1.6 €/kg is assumed. The material cost calculation process of different components is described in more detail in Appendix B.
Once the dimensions are known, the processing cost can be estimated. The processing cost evaluation is a simplified version of that described in [32]. In segmentally baffled STHEs, the processing costs related to baffles represent over half of the total processing cost. The baffle-related costs (cutting, beveling, drilling, and tube bundle assembly) were, thus, modelled similarly to [32] by estimating the total processing length L pr and speed v pr to obtain the processing time t pr = L pr /v pr , specific cost per time c pr , and, from there, the total cost C pr = c pr t pr . The sum of other processing costs was estimated at 10% of the material cost C mat,HX based on the results of [32]. The total sum of processing costs of the different operations x are then calculated as where t su,x is the set-up time for each operation x, and the final term is the approximate cost of all other manufacturing processes, estimated at 10% of material cost. The equations and assumptions used for processing cost estimation are summarized in Appendix C.

Reactor Pressure Vessel and Containment Vessel Cost Models
The RPV and CV both consist of a cylindrical section, hemispherical top, and bottom heads. A detailed design of the RPV is not yet available at this stage. The costs due to the various features added to the bare pressure vessel components were estimated as representing a similar fraction of the total material cost as in the NuScale SMR concept, which utilizes a somewhat similar elongated RPV [33].
The dimensions of the RPV and the CV are found by estimating the straight section of the RPV cylinder to continue 1.5 m above the top of the heat exchanger and a distance ∆H − 0.5D RPV below (see Figure 4). A 0.5 m gap separates the CV inner surface from the RPV outer surface at the sides and bottom, i.e., D CV,i = D RPV,o + 1.0 m. The cylindrical part of the CV shell is 1 m longer than the RPV. representing a similar fraction of the total material cost as in the NuScale SMR concept, which utilizes a somewhat similar elongated RPV [33].
The dimensions of the RPV and the CV are found by estimating the straight section of the RPV cylinder to continue 1.5 m above the top of the heat exchanger and a distance Δ − 0.5 below (see Figure 4). A 0.5 m gap separates the CV inner surface from the RPV outer surface at the sides and bottom, i.e., , = , + 1.0 m. The cylindrical part of the CV shell is 1 m longer than the RPV. Using the price data in [34], a value of 13.6 €/kg is obtained for the cylindrical shell and the heads when corrected for the 2021 producer price index for the manufacture of fabricated metal products, excluding machinery [35]. This increases to 18.3 €/kg with an estimated 4.7 €/kg cost difference between carbon steel and ASTM 316L. The other vessel components, e.g., possible conical sections, flanges, and openings (but not the core and related equipment), were estimated to constitute 40% of the NuScale RPV cost [33,34]. With this assumption that the shell and heads, whose mass is obtained from the sizing, constitute 60% of the total material costs as in [34], the adjusted material cost that can be Using the price data in [34], a value of 13.6 €/kg is obtained for the cylindrical shell and the heads when corrected for the 2021 producer price index for the manufacture of fabricated metal products, excluding machinery [35]. This increases to 18.3 €/kg with an estimated 4.7 €/kg cost difference between carbon steel and ASTM 316L. The other vessel components, e.g., possible conical sections, flanges, and openings (but not the core and related equipment), were estimated to constitute 40% of the NuScale RPV cost [33,34]. With this assumption that the shell and heads, whose mass is obtained from the sizing, constitute 60% of the total material costs as in [34], the adjusted material cost that can be used to obtain the pressure vessel cost with the mass of the main parts (cylindrical shell and heads); m s+h , is obtained from as 30.54 €/kg. The manufacturing costs of the RPV and CV are then found, assuming a work-to-material cost ratio of 0.657 [34]. Although the stainless-steel vessel does not need cladding, for a conservative estimate, the work cost was not adjusted down. FOB and TCI costs are estimated using similar mark-up and installation factors as with the heat exchangers, Equations (4) and (5).

Reactor Pool Cost Model
The 7-m-deep main pool is unaffected by the heat exchanger design and is considered part of the fixed-cost C SMR . Only the cost of the cylindrical cavity extending below the main pool, C cav , varies as a function of the heat exchanger dimensions and the height DH to the RPV bottom and is, thus, evaluated in the cost model. This consists of two components: excavation C excav , and wall construction C cav,w .
The cavity has a 0.5-m-thick steel-reinforced concrete wall with 10 mm steel lining. The concrete cost is the sum of material and labor costs per m 3 . The material cost is obtained from the rebar-to-concrete mass ratio r, 0.11 for reactor buildings in Gen IV nuclear systems, the specific costs from [36], index-corrected according to [37]. Excavation cost is obtained from the estimated cavity volume, as described in Supplementary Material S1.

Decision Variables and Constraints
The decision variables and their feasible ranges are listed in Table 4. In addition to the feasibility ranges, the following additional constraints are considered: • Maximum shell-side velocity w sh,max < 1.5 m/s; • Maximum shell-side pressure drop Dp sh < 1.0 bar; RPV straight length below heat exchanger ∆H < 2.0 m.

Thermohydraulic Model
The heat exchanger model solves the heat exchanger performance as a rating problem for a given geometry. Because the operation is based on the natural circulation of the hot fluid, only the inlet state and mass flow rate of the cold flow and the total heat transfer rate are fixed and known. The heat exchanger performance (heat transfer and pressure drop) and the inlet and outlet state both affect each other, requiring an iterative solution to determine the operating point. In addition to the heat exchanger itself, the reactor core pressure drop also affects circulation. Pressure drops in the riser and downcomer outside of the core and the heat exchangers are considered negligible.
The solution process, shown as a flow chart in Algorithm 1, is based on two nested iteration loops. The outer loop adjusts the primary water inlet temperature to match the Φ tgt = 3.125 MW (50MW/16) target heat rate in a heat exchanger, while the inner adjusts the primary water flow to match the total pressure drop (reactor core plus heat exchanger), with the ∆p tgt available from the natural circulation. Damping factors α are used and adjusted during the iteration to ensure convergence. The primary water outlet temperature in step 3.1.3 is solved from tube outside area A o and overall heat transfer coefficient U using the ε-NTU method, based on dimensionless parameters NTU, ε, and C*: The overall heat transfer coefficient U is found from where R" tf,i and R" tf,o are inside and outside fouling resistances and R" w , the tube wall resistance. The fouling is estimated to be similar to a steam power plant low-pressure feed heater. Based on the Heat Exchanger Institute's recommendation [38] and a nuclear power plant feed heater diagnostics system [39], values of R" tf,i = 3.5·10 −5 m 2 K/W (feed heater fouling resistance 1 year from start-up) and R" tf,o = 5.3·10 −5 m 2 K/W (feed heater drain cooling section tube outside resistance) are used. Ideally, the tube-side flow should be fully turbulent, with Reynolds number Re > 10 4 . This ensures good, predictable heat transfer and a small area. Unfortunately, with natural circulation, this cannot be guaranteed: transition or fully laminar flow is a risk, especially at part-load when using small tubes. Small-diameter tubes can also result in a hydraulically smooth surface condition not being met. This can be advantageous, as the partially rough flow regime increases the heat transfer coefficient, although, unfortunately, the literature on the partially rough regime heat transfer is both dated and sparse. As none of the aforementioned conditions can be ruled out, all are implemented in the model.
The tube inside heat transfer coefficient h i in transition and fully turbulent regime (Re > 2300) is evaluated primarily using a variant of the Gnielinski correlation [40], where the Darcy friction factor f D for hydraulically smooth tubes is obtained from the Bhatti-Shah correlation [41], Equation (13) is said to be valid for 2300 ≤ Re ≤ 5·10 6 ; while it can be considered accurate for water at 50-150 • C range at the fully turbulent Re > 10 4 regime, yielding almost identical results to the Petukhov-Popov correlation [42] from which it is derived, it is not recommended for the 2300 < Re < 10 4 transition regime [43] where the flow intermittently switches between laminar and turbulent [44]. Based on this observation in [44], a linear interpolation between fully turbulent Re = 10 4 and fully laminar Re = 2300 heat transfer coefficient, where the intermittency factor γ is has been proposed in [45] as cited in [46]. To find a conservative estimate for the purpose of natural circulation feasibility evaluation, heat transfer in the 2300 < Re < 10 4 regime is modelled by evaluating the Nu with both (13) and (15), and using the lower result. The turbulent Nu turb,Re=10 4 in Equation (15) is obtained with Equation (13). The Nu lam,Re=2300 depends on the tube wall boundary condition. Considering constant heat flux a better approximation than the constant temperature in a counter-current heat exchanger, the laminar Nu is thus evaluated according to [46]: When surface roughness protrusions extend beyond the laminar sublayer in turbulent flow, the hydraulically-smooth surface approximation no longer applies. Roughness effects are correlated with a roughness Reynolds number e + , where e is the sand-grain roughness and C f the Fanning coefficient of friction. Although fully rough (e + > 70) flow is unlikely, partially rough (5 < e + < 60-70) [43] cannot be ruled out. In this region, the turbulent Nu is adjusted according to Norris [47] as cited in [48], where the C f is obtained from the explicit approximation of the well-known iterative Colebrook-White equation according to Chen [49],

Shell-Side Heat Transfer Coefficient
Calculating the shell-side heat transfer coefficient h o follows the Bell-Delaware method, using the Nu correlations and correction factors f according to [50] in where the factors f are for tube arrangement (f A ), turbulence development (in multi-baffled heat exchangers f N = 1), deviation of the shell-side main flow, path A in Figure 5, from pure cross flow (f G ), bypass flow B (f B ), leakage flows C and D through baffle-shell and baffletube gaps (f L ), and fluid property variation (f P ). Except for fluid property variation f P , and f B at extremely low-Re flows, these are determined from the geometry alone, unaffected by flow parameters. The equations for obtaining these are listed in [50].

Shell-Side Heat Transfer Coefficient
Calculating the shell-side heat transfer coefficient ho follows the Bell-D method, using the Nu correlations and correction factors f according to [50] in The turbulent and laminar Nu in Equation (14)  where w is the cross-flow velocity at shell centerline between two baffles without the streamed length ½πdo, the density, µ the dynamic viscosity, and ψ a void fra account for the fraction of free flow area left between the tubes as defined in [50]. The turbulent and laminar Nu in Equation (14) are obtained from Equations (22) and (23), with Re ψ,1 defined as where w is the cross-flow velocity at shell centerline between two baffles without tubes, l the streamed length 1 2 πd o , ρ the density, µ the dynamic viscosity, and ψ a void fraction to account for the fraction of free flow area left between the tubes as defined in [50].
For shell-side calculation, the number of tubes N tb is estimated according to [33] as where D OTL is the outer tube limit (see Figure 6) and P the tube pitch; of these the number of tubes in the tube window N tb,w is estimated as the same as the ratio of the bundle segment area in the window A W , based on segment height h bndl and diameter D OTL , to the area of the whole circular bundle of same diameter D OTL .

Pressure Drop, Heat Exchanger
The tube-side Δpi is found using the Darcy-Weisbach equation, where L is the tube length, the loss coefficients K are Ki,in and Ki,out into and out from the tubes, respectively. Constant , = 1 is assumed, while , is determined as a function of P and [51]. When the tube is hydraulically smooth, e + < 5, fD is obtained from Equation (14), otherwise as fD = 4Cf, with the Cf from Equation (20).
Shell-side pressure drop Δ po is the sum of Δ p in cross-flow between baffles Δ pQ, entry and exit cross flow sections ΔpQE, baffle windows ΔpW, and at the nozzles Δpnzl: where Nbf is the number of baffles; the different components are depicted in Figure 7. The calculation process follows the methodology described in [52].

Pressure Drop, Reactor Core
The reactor core consists of 37 otherwise typical 17 × 17 pressurized water reactor (PWR) fuel assemblies, but with a shorter 1-m active length (total length assumed Lcore = 1.3 m), corresponding to the LDR-50 DH reactor concept [14,15]. The core pressure drop Dpcore consists of the sum of fuel rod surface friction losses, and form losses due to the spacer grids, top and bottom nozzles, and the inlet contraction and outlet expansion losses, and the gravity loss. The gravity loss is accounted for by considering the riser height down to midpoint of the core in determining the natural circulation driving force.
The friction and form losses are evaluated with Equation (26) using the mean as velocity, and Lcore and subchannel hydraulic diameter dh as the length and diameter [53]. Friction factor fD is obtained using Equation (16) for fully turbulent flow, 64/Re for laminar. At 2300 < Re < 10 4 , linear interpolation between 64/2300 and Equation (16) at Re = 10 4 is used.
All form losses including the grid spacers, core top and bottom structures, and inlet

Pressure Drop, Heat Exchanger
The tube-side ∆p i is found using the Darcy-Weisbach equation, where L is the tube length, the loss coefficients K are K i,in and K i,out into and out from the tubes, respectively. Constant K i,out = 1 is assumed, while K i,in is determined as a function of P and θ tp [51]. When the tube is hydraulically smooth, e + < 5, f D is obtained from Equation (14), otherwise as f D = 4C f , with the C f from Equation (20). Shell-side pressure drop ∆p o is the sum of ∆p in cross-flow between baffles ∆p Q , entry and exit cross flow sections ∆p QE , baffle windows ∆p W , and at the nozzles ∆p nzl : where N bf is the number of baffles; the different components are depicted in Figure 7. The calculation process follows the methodology described in [52].

Pressure Drop, Heat Exchanger
The tube-side Δpi is found using the Darcy-Weisbach equation, where L is the tube length, the loss coefficients K are Ki,in and Ki,out into and out from the tubes, respectively. Constant , = 1 is assumed, while , is determined as a function of P and [51]. When the tube is hydraulically smooth, e + < 5, fD is obtained from Equation (14), otherwise as fD = 4Cf, with the Cf from Equation (20).
Shell-side pressure drop Δ po is the sum of Δ p in cross-flow between baffles Δ pQ, entry and exit cross flow sections ΔpQE, baffle windows ΔpW, and at the nozzles Δpnzl: where Nbf is the number of baffles; the different components are depicted in Figure 7. The calculation process follows the methodology described in [52].

Pressure Drop, Reactor Core
The reactor core consists of 37 otherwise typical 17 × 17 pressurized water reactor (PWR) fuel assemblies, but with a shorter 1-m active length (total length assumed Lcore = 1.3 m), corresponding to the LDR-50 DH reactor concept [14,15]. The core pressure drop Dpcore consists of the sum of fuel rod surface friction losses, and form losses due to the spacer grids, top and bottom nozzles, and the inlet contraction and outlet expansion losses, and the gravity loss. The gravity loss is accounted for by considering the riser height down to midpoint of the core in determining the natural circulation driving force.
The friction and form losses are evaluated with Equation (26) using the mean as velocity, and Lcore and subchannel hydraulic diameter dh as the length and diameter [53]. Friction factor fD is obtained using Equation (16)

Pressure Drop, Reactor Core
The reactor core consists of 37 otherwise typical 17 × 17 pressurized water reactor (PWR) fuel assemblies, but with a shorter 1-m active length (total length assumed L core = 1.3 m), corresponding to the LDR-50 DH reactor concept [14,15]. The core pressure drop Dp core consists of the sum of fuel rod surface friction losses, and form losses due to the spacer grids, top and bottom nozzles, and the inlet contraction and outlet expansion losses, and the gravity loss. The gravity loss is accounted for by considering the riser height down to midpoint of the core in determining the natural circulation driving force. The friction and form losses are evaluated with Equation (26) using the mean w core as velocity, and L core and subchannel hydraulic diameter d h as the length and diameter [53]. Friction factor f D is obtained using Equation (16) for fully turbulent flow, 64/Re for laminar. At 2300 < Re < 10 4 , linear interpolation between 64/2300 and Equation (16) at Re = 10 4 is used.
All form losses including the grid spacers, core top and bottom structures, and inlet (contraction) and outlet (expansion), are estimated with loss coefficients. Core inlet and outlet loss coefficients are estimated at 0.5 and 1.0, respectively [53]. For the grid spacers and top and bottom structures the loss coefficients are obtained using the correlation for PWR spacers [53], where ε is the ratio between the projected area of the obstruction (grid) to the flow area without the obstruction. The loss coefficients are determined assuming the fuel assemblies are shortened versions of the AP1000 fuel assemblies. While the details of the grid spacers and other obstructions are not known, the number and type of spacers as well as the typical pressure drop in AP1000 conditions are known [54]. From this, representative values for ε can be found, which are then used in Equation (28) for the shorter fuel design in the DH reactor conditions, considering a total of three grid spacers, and top and bottom structures.

Modified Cuckoo Search Algorithm
The Cuckoo Search (CS) is a metaheuristic global optimizer [55] metaphorically based on the well-known brood parasitism of cuckoo birds [56], as well as a search pattern known as Lévy flights observed in many foraging species in nature [56]. The CS is a populationbased method; each candidate solution to an optimization problem represents a cuckoo egg. Each egg is represented by a vector x of D decision variable values.
The literature on the CS is unfortunately somewhat contradictory. Earlier two very similar, one somewhat modified, and a fourth, significantly different variant, all described by the authors of the method as simply "Cuckoo Search", were identified. The differences were described in [22] and in further detail in [17]. The first three variants, labelled CS1 sub-variants [17], proved unsuccessful in heat exchanger optimization, but the fourth, CS2, was highly successful, and was modified further into the CS3 in [17]. Both CS2 and CS3 borrow heavily from another well-known metaheuristic, the differential evolution (DE). In this study, the CS3 was implemented.
Like all CS variants the authors are aware of, the CS3 consists of two steps. The first is based on DE [57]. The differences to DE are the use of a randomized mutation weight, and each target vector x i against which the trial vector u i competes serving also as the base vector in generating u i , rather than the latter being separately and randomly selected. All population members serve as target vector once per iteration.
In [17] the greedy local-to-best mutation strategy proved fastest, while still remaining robust. In this case, the greater number of decision variables, many of which are discrete, is likely to produce a severely discontinuous and multi-modal objective function topology. As a thorough meta-optimization was ruled beyond the scope, the more conservative CS3/rand/1/bin was selected over the greedy local-to-best variant. The rand/1/bin trial vector u i generation gives each decision variable d a value according to where the mutation weight F d is drawn from a size D-dimensioned normal-distributed random-number vector F (mean = 0.8), and ε d is a uniform-distributed random num- ber, ε d (0,1). For a rotationally invariant search effective at non-separable problems, the crossover probability CR should have a small value. The second step is a Lévy-distributed random walk. While these Lévy flights appear generally less efficient than the differential mutation at most problems, they have one important advantage: retaining a capability for long leaps in the objective function landscape even after the population has converged into a small area. This preserves some global search capability. In the CS3 variant of Cuckoo Search, the tuning parameter p a ("switching probability") determines the fraction of population on which the Lévy flights are implemented. After the differential mutation step, the population is ranked, and a fraction p a of the worst eggs are subjected to Lévy flights. Each egg serves once as a base vector x i , using its distance to a randomly chosen vector x r1 to generate a trial vector u, where α is a scaling factor, s a Lévy-distributed step size vector, n a vector of normaldistributed random numbers, and operator • the Hadamard product. Considering the role of the Lévy flights as one of global search to prevent stagnation on local optima, even when the distance ∆x best is small, a comparatively large value of α = 0.1 is used. The vector s is determined using the Mantegna's algorithm [58] where p and q are normal-distributed random-variable vectors of size D, a mean of zero, and variances of 1 and σ 2 , respectively, with Based on recommendation in [56], Lévy exponent of β = 1.5 is used. A trial vector is evaluated by competition against the original base vector; the one with better objective function value survives to the next generation. The pseudo-code of CS3/rand/1/bin is found in Appendix D (Algorithm A1), and MATLAB code of the algorithm in Supplementary Material S3. A population size of NP = 45, crossover probability CR = 0.1, and p a = 0.20 were used in all runs.
Due to the multi-constrained nature of the problem, often all members of the initial population are infeasible. Constraints are handled by penalty functions, implemented so that an infeasible solution always loses to any feasible one, no matter how poor; between two infeasible solutions, the one in breach of higher number of constraints always loses; and between two in breach of one constraint, the one in greater magnitude of violation is likely to lose. This ensures a rapid convergence towards feasible regions.

Results
The convergence behavior of the CS3 algorithm appeared robust, albeit somewhat slow, and requiring a fairly large population size of 5D. While twice greater than the 2.5D population and a greedier mutation strategy with the lower-dimensioned, less multimodal condenser optimization in the original introduction of the CS3, this is still half of what is common with the classic DE on which the CS3 is based. The convergence behavior at the three cases with 10 trial runs per case can be seen in Figure 8.
Once the number of objective function evaluations (NFE) reached 50,000, the search was terminated and re-started with halved population size and a region limited around the solutions found by the 10 runs for another 10,000 evaluations to obtain the exact solutions, listed in Table 5. The full data are available from Supplementary Material S2. Once the number of objective function evaluations (NFE) reached 50,000, the search was terminated and re-started with halved population size and a region limited around the solutions found by the 10 runs for another 10,000 evaluations to obtain the exact solutions, listed in Table 5. The full data are available from Supplementary Material S2.   The results show that the primary-circuit natural circulation design appears feasible for a DH reactor design, including partial-load operation. The designs are very similar. All designs use small 12 mm tubes-the smallest diameter, where the 3.2-mm tube sheet ligament minimum does not yet significantly increase the tube pitch ratio from its 1.25 minimum. The small tube and large shell diameters allow for placing the necessary heat transfer surface in a smaller length. There is a 2.6-m empty height in the RPV up to the bottom of the heat exchangers: the fully cooled, dense primary fluid there adds to the natural circulation force more than a similar RPV height with taller heat exchangers would.
Clear differences only appear in the flow parameters of the summer load point LP6, where the operating conditions of BL cases remain similar to the winter LP1, whereas in the case of ML, the reduced load results in clearly reduced flow rates and primary-side temperature. This results in reduced heat transfer coefficients and pressure drops (Figure 9), and transition-regime flow at the tube exit. All other cases, and most load points of ML50, are fully turbulent. The results show that the primary-circuit natural circulation design appears feasible for a DH reactor design, including partial-load operation. The designs are very similar. All designs use small 12 mm tubes-the smallest diameter, where the 3.2-mm tube sheet ligament minimum does not yet significantly increase the tube pitch ratio from its 1.25 minimum. The small tube and large shell diameters allow for placing the necessary heat transfer surface in a smaller length. There is a 2.6-m empty height in the RPV up to the bottom of the heat exchangers: the fully cooled, dense primary fluid there adds to the natural circulation force more than a similar RPV height with taller heat exchangers would.
Clear differences only appear in the flow parameters of the summer load point LP6, where the operating conditions of BL cases remain similar to the winter LP1, whereas in the case of ML, the reduced load results in clearly reduced flow rates and primary-side temperature. This results in reduced heat transfer coefficients and pressure drops ( Figure  9), and transition-regime flow at the tube exit. All other cases, and most load points of ML50, are fully turbulent. Another clear conclusion of the results is that the reactor pressure vessel and containment vessel costs far exceed the costs of the heat exchangers themselves. The cost breakdowns are shown in Figure 10 in terms of both TCIs and produced heat cost-price. Figure 9. Annual variation of main flow parameters in case ML50 (a) and BL50 (b). All data, except core Dp, refers to conditions in a single heat exchanger (one of 16). Case BL200 exhibits similar behavior to BL50.
Another clear conclusion of the results is that the reactor pressure vessel and containment vessel costs far exceed the costs of the heat exchangers themselves. The cost breakdowns are shown in Figure 10 in terms of both TCIs and produced heat cost-price. As the heat exchanger design clearly will impact the dimensions and cost of vessels (and the pool cavity), but the designs of these components are still at a very preliminary stage, the initial heat exchanger design cannot be considered final. Finally, the optimal design can only be found once the designs and cost models of said components approach the final stage as well, resulting in an iterative design process.

Discussion
The results indicate that a small district heating reactor operating in the natural circulation principle is a feasible concept under different load profiles and operating conditions. While minor differences between the configurations optimized for different scales and load profiles were found, the differences are small enough to be arguably As the heat exchanger design clearly will impact the dimensions and cost of vessels (and the pool cavity), but the designs of these components are still at a very preliminary stage, the initial heat exchanger design cannot be considered final. Finally, the optimal design can only be found once the designs and cost models of said components approach the final stage as well, resulting in an iterative design process.

Discussion
The results indicate that a small district heating reactor operating in the natural circulation principle is a feasible concept under different load profiles and operating conditions. While minor differences between the configurations optimized for different scales and load profiles were found, the differences are small enough to be arguably negligible, permitting a single standardized design of all main components: the heat exchangers, the reactor pressure vessel (RPV), and the containment vessel (CV).
A central feature of the optimal solution is the compact packaging of the heat transfer surface as high as possible in the RPV downcomer, dividing the total riser height approximately in half between the heat exchangers at the top and empty space filled with fully cooled primary water at the bottom half. This allows to maximize the pressure difference driving the natural circulation.
As the designs of the RPV and CV are not yet finalized, the results obtained here must also be treated as tentative, subject to change, as the designs and cost functions of said pressure vessels, as well as the pool cavity and the reactor core, become more detailed and accurate. The heat exchanger cost model itself, as well as some of the constraints, are currently also still on a general level. In a final design, data on the equipment available to the chosen manufacturers can be used to obtain a design tailored for such. Such final optimization should also consider the reactor pressure and temperature levels, considered as fixed constraints in this initial study, as free decision variables.
In terms of optimizer performance, the modified cuckoo search proved effective. Although considerably more conservative tuning parameters settings were needed than with the condenser optimization for which the method was originally introduced, the population size was still only half of what would typically be required with differential evolution. Finally, although analyzing the final cost competitiveness of the considered heating reactor concept was not a goal of the work, some conclusions can still be drawn. The obtained annual net cash flow rates were clearly positive for all cases, even though conservative or even pessimistic assumptions were applied throughout the process. The cost-price of the produced heat, 30-55 €/MWh, is well under the average selling price of district heat in Finland. At base-load operation, the 30-40 €/MWh cost-price is broadly similar to what could be obtained using a heat pump system with 880 €/kW th CAPEX and an average of 40 €/MWh purchased electricity price with 20 €/MWh grid and tax costs while avoiding vulnerability to the sometimes very high electricity costs during peak-load conditions. At mid-load operation, the 35-55 €/MWh heat price is broadly similar to or slightly greater than the cost price for a wood chip-fired heat-only boiler. While biomass-fired heat and power generation are likely to remain important in Scandinavian conditions for the foreseeable future, the availability of sustainably produced solid biofuels for heat generation is limited and likely insufficient to fulfill the heat demands, especially in larger population centers in, e.g., Finland or Sweden. Reactor concepts, such as the one considered in this study, thus, appear economically feasible and potentially important components of a cost-effective, clean energy system.

Conclusions
The optimization of the primary heat exchangers of a small, 50-MW district heating reactor circulation natural using a modified cuckoo search algorithm proved successful. The following main conclusions from the study were drawn: • A natural-circulation concept for a heat-only reactor appears feasible for both steady base-load and load-following operation. The optimized configurations for different scales and load profiles are near-identical; • Central to the optimized heat exchanger configurations was using small-diameter tubes at minimum allowable tube sheet ligament thickness. The resulting high density of heat transfer area per volume enables keeping the heat exchangers in the upper half of the downcomer, filling the lower half with fully cooled primary water. This maximizes the pressure difference available to drive natural circulation; • The reactor pressure vessel and containment vessel dominate the cost impact of the heat exchanger designs; • A heat-only reactor producing hot water at modest temperatures has the potential to be a highly competitive carbon-neutral DH producer. The two central reasons behind this are the lightweight construction possible for a low-temperature, low-pressure concept and the ability to convert the reactor thermal power into an energy product without the heat rejection losses inevitable in an electricity-producing power cycle; • The modified cuckoo search proved successful for the optimization task, although compared to an earlier condenser optimization study, using considerably more conservative control parameters settings emphasizes robustness over speed. The minimum acceptable shell thickness s sh,min of a cylindrical pressure vessel shell (the heat exchanger shell, or the cylindrical parts of the RPV or the CV) is the sum of required thickness s s,0 , corrosion allowance c 1 , and manufacturing tolerance c 2, s sh,min = s sh,min,0 + c 1 + c 2 (A1) where the required thickness s s,min,0 is equal to the greatest of the results obtained from sizing a cylindrical pressure vessel body against elastic buckling, s cyl,eb , non-elastic buckling, s cyl,neb , and internal pressure, s cyl,ip . For elastic buckling, an explicit equation for the thickness is not stated in [31]; a maximum acceptable pressure p max,eb for a given thickness is defined instead, where n SF is a safety factor, n SF = 3 with out-of-roundness of 0.015 or less, the auxiliary parameter C is obtained from and N bw is the number of buckling waves. Equation (A2) is solved for all integer numbers N bw ∈ {2, . . . , N bw,max }, the smallest obtained pressure defining the p max,eb . The upper limit N bw,max for the number of buckling waves is found by rounding up to the next integer the value N bw,max from Sizing against elastic buckling with the implicit Equation (A2) is performed by converting the sizing problem into a one-dimensional optimization problem to minimize the difference between the design pressure and the calculated pressure p calc obtained from Equation (A2), where s eb is the solution to min|p calc − p des | = f (s), applying Equation (A2) as the function f. The optimization is performed with Golden Section Search (GSS), as described in [59]. Two initial guesses, s a and s b , are chosen; one low and the other high enough to bracket the solution. These are then iteratively moved closer until s a -s b ≤ 0.01 mm. From the dimensions, the total volume of material can be obtained, and, with a density of r 316L = 7890 kg/m 3 , the mass. Against non-elastic buckling, the maximum pressure for wall thickness s is where a maximum acceptable shell out-of-roundness value of e = 0.015 is applied. The thickness is then obtained as with the elastic buckling by converting the problem to a onedimensional optimization problem minimize the difference between design and calculation pressures, Equation (5), and applying Equation (A6) for the function f.
Against internal pressure p, the minimum wall thickness of a cylindrical vessel s ip is

. Tubesheet Thickness
The thickness of the tubesheet s ts is obtained from where p max is the greater of the hot and cold fluid design pressures; here, 1.0 MPa, C is a constant with a value of C = 0.50 for tube sheets welded from one side to the shell and exceeding a thickness of three shell thicknesses 3s sh , and v is a correction factor, where the weld strength correction v W for austenitic stainless steel can be assumed as 1.0 when all welds are inspected, and v H is the hole field correction,

Appendix B. Heat Exchanger Material Costs
The heat exchanger material cost C mat,HX is obtained using the material volume, density, and specific costs. Table A1 lists the equations for the volumes of the main components, as well as the specific costs c. Some amount of the material, purchased as tubes, pipe, plates, and rods, is also lost as scrap in processing; a value of c sc = 1.6 €/kg is assumed for this. The baffle plate and tube sheet costs already account for the estimated amount of scrap from cutting discs or segmentally-cut baffle plates from rectangular plate. Additional scrap is also produced from drilling the holes in the tube sheet and baffle plates V sc,h , V sc,h = (N bf N h,bf s bf + 2·N tb s ts ) π 4 d 2 tbh . (A20) Scrap is also produced when standard 6-m heat transfer tubes are cut to the specified length and drilling holes. This is accounted for by determining an effective tube price c tb,eff = c tb + (c tb − c sc ) · mod(6.000; L tb,tot ) 6.000 − mod(6.000; L tb,tot ) , with which the tube material cost is evaluated. With these parameters, the material cost is found using the components volumes, specific costs, and material density,

Appendix C. Heat Exchanger Processing Cost
The costs per operation include the sum of total costs of labour, machinery amortization, and consumables and utilities consumed; assembly considers labour only. Baffle drilling is assumed to take place by placing and aligning the baffles precisely on top of each other and fastening them tightly for drilling through multiple (but no more than five) plates at a time. It should be noted that the processing speed and cost estimates are highly uncertain and can vary greatly between the manufacturers, depending on the type of equipment used. Table A2. Processing costs; baffle dimensions those depicted in Figure 6.

Process x Processing Length L x [m]
Speed v x [ m min ]

Spec. Cost c t,x [€/h]
Cutting baffles (bf) and sealing strips (ss) L cut = N bf L cut,1bf + 2N ss (2L ss + 2W ss ) v cut = 0.5 * 100 Beveling L bev = L cut (A25) v bev = 1.5 * 60 Drilling L dr = N bf N h,bf s bf + ceil N bf 5 L lt + L pt + L ot , L lt = 3.0 mm, drill head lead travel [32] L pt = 5.0 mm, drill head pre-travel [32] L ot = 5.0 mm, drill head over-travel [32] (A26) v dr = ln d o 18.2 + 1 83 § 80 Bundle assembly n/a n/a 50 * estimate based on [32] with the assumption of stainless steel processing speed being 50% of that of carbon steel. § [d o ]=mm; yields result as m/min. The function is a curve fit on a drill bit of manufacturer's recommended feed rate and RPM with carbide tips for different hole diameters on a stainless steel plate.