Two Modeling Approaches for sCO₂-PCHEs Transients

Abstract

Printed circuit heat exchangers (PCHEs) are ideal for sCO₂ power cycles due to their compactness, effectiveness, and high-pressure capability. However, their unique architecture complicates the modeling of their dynamic behavior in power cycles, which experience rapid transients. Both the significant computational resources required and the high investment costs of experiments limit their widespread application. To model the component- and system-level transients of sCO₂-PCHEs, this study presents two 1D modeling approaches for different purposes: one for component-level simulation based on local properties of sCO₂, and the other for system-level simulation based on transfer functions. Given the significant discrepancies observed when using a fixed time constant in the latter approach, this paper introduces the concept of an optimal time constant to model the transient behavior of PCHEs as a first-order system with minimal prediction error. This optimal time constant varies with operating conditions, contrary to what the name might suggest. These results demonstrate the potential of PCHEs in advanced power cycles and provide valuable insights for accurate system-level sCO₂ power cycle control studies.

1. Introduction

Supercritical carbon dioxide (sCO₂) power cycles are promising solutions for transforming traditional thermal powers to renewable energies [1, 2]. The sCO₂ cycles find massive applications across vast working scenarios thanks to their high conversation efficiency but low carbon emissions. Among various types of heat exchangers, printed circuit heat exchangers (PCHEs) are favored as primary heaters and recuperators in sCO₂ cycles because of their outstanding heat transfer performance and remarkable resistance to high temperature and pressure. Massive experimental and numerical studies have been conducted to investigate sCO₂-PCHEs’ dynamic performance under varying scenarios, such as nuclear power [3, 4], concentrated solar power [1], wasted heat recovery [5, 6] and hypersonic aero-engines [7]. For these investigations, the most challenging task is to handle the difficulties brought by sCO₂ fluid’s steep property variations.

However, the temperature-dependent properties violate the assumption of ‘constant property’ employed by most dynamic modeling approaches [8, 9], resulting in significant and lagging prediction errors, for which both three- and one-dimensional simulation techniques have paid special efforts to address the issue. The high-fidelity three-dimensional modeling approaches, such as Computational Fluid Dynamics (CFD) [10, 11] and Direct Numerical Simulation (DNS), capture the dramatic property variations at the price of magnificent computation cost [12]. The one-dimensional modeling methods prefer the discretized approach to capture local properties while requiring remarkable computational costs [13]. On the other hand, AI-assisted technologies have shown promise in enhancing power plant efficiency, but their application in predicting system-level off-design performance and optimizing control strategies remains underexplored [14, 15]. This limitation is primarily due to the scarcity of operational data needed to train AI models effectively. The data, which are derived from experimental tests, require significant investment costs, making them accessible only to a few well-equipped laboratories [16].

Utilizing inappropriate empirical correlations might induce unexpected convergence and stability issues to dynamic modeling [17, 18]. The dynamic simulation, aiming to provide insightful knowledge of sCO₂-PCHEs’ transients for control studies, prefers the lower-order system to proxy PCHEs’ transients with acceptable errors. While experimental observations and analytical solutions have confirmed the feasibility of using the first-order system approximation for shell-and-tube and plate heat exchangers [19, 20], the same modeling technique has not been reported for sCO₂-PCHEs, as we know. Therefore, this chapter will propose a component-level 1D model based on sCO₂ local properties and a first-order approximation for sCO₂-PCHEs. The feasibility of both models under different operating conditions are compared. A data-fitting approach will also be presented to compute the key parameter, time constant, of the proposed low-order approximation. As will be seen, a varying time constant, determined by both geometrical parameters and operating conditions, can significantly reduce the approximation’s errors. By the end of this study, the transients of sCO₂ PCHEs can be approximated with minimal cost both experimentally and numerically, with good accuracy.

2. Modeling methodologies

2.1 1D numerical model for sCO${}_2$-PCHEs

A sCO₂ Recompression Brayton Cycle (Figure 1a) usually chooses PCHEs as its precooler, high-temperature and low-temperature recuperator. The sCO₂ fluid experiences dramatic thermal-hydraulic variations in these components, as shown in the changes of Prandtl numbers in Figure 1b and c.

Figure 1.

sCO₂’s property variations per components of recompression Brayton cycles.

To investigate the dynamics of sCO${}_2$ power cycles, 1-D numerical modeling [2] usually discretizes a straight- or zigzag-channel PCHE into $N-1$ heat transfer units (Figure 2), where each unit consists of one hot volume, one cold volume, and one separating metal wall volume. For the $j_{\mathit{th}}$ unit $j=\left(1,\dots ,N-1\right)$, the governing conservation equations for mass, energy, and momentum for either the cold or hot side are,

Figure 2.

Schematic of discretized PCHE model.

$$\frac{\partial M_j}{\partial t}=A\cdot l\cdot \left[\frac{\partial {\rho }_{:j}}{\partial p_{:j}}\cdot {\dot{p}}_{:j}+\frac{\partial {\rho }_{:j}}{\partial T_{:j}}\cdot {\dot{T}}_{:j}+\frac{\partial {\rho }_{:j+1}}{\partial T_{:j+1}}\cdot {\dot{T}}_{:j+1}\right]+{\dot{w}}_{:j}-{\dot{w}}_{:j+1}$$E1

$$A\cdot l\cdot {\rho }_j\cdot c_{v,j}\cdot {\dot{T}}_j=w_{:j+1}\cdot h_{:j+1}-w_{:j}\cdot h_{:j}+Q_{x,j}$$E2

$$\frac{l}{A}\frac{\partial w_j}{\partial t}=p_{:j}-p_{:j+1}-\mathrm{\Delta }{p}_{\mathit{\text{fric}},j}$$E3

where the mass $M_j$, averaged density ${\rho }_j$, averaged pressure $p_j$, averaged temperature $T_j$, averaged specific heat capacity $c_{v,j}$, and mass flow rate $w_j$ are volume-average quantities of volume $j$, $A$ is the cross-sectional area (identical for node and volume), $l=L/\left(N-1\right)$ is the volume length, ${\rho }_{:j},p_{:j},T_{:j},h_{:j},w_{:j+1}$ are averages of density, pressure, temperature, specific enthalpy, and mass flow rate over the cross-section area of node $:j$, and $\mathrm{\Delta }{p}_{\mathit{\text{fric}},j}$ is the local pressure drop through the volume caused by the friction.

For the energy conservation equation (Eq. (2)), $Q_{x,j}$ denotes the heat transferred between a flow volume $j$ and its adjacent wall volume $j$, where the subscript $x$ is either $h$ for hot side or $c$ for cold side, such as,

$$Q_{x,j}=\left(T_{w,j}-T_{x,j}\right)\cdot \omega \cdot l\cdot {\mathit{hc}}_{x,j}$$E4

$${\mathit{hc}}_{x,j}=\frac{{\mathrm{Nu}}_{x,j}\cdot k_{x,j}}{D_{\mathit{hyd}}}$$E5

where $T_{w,j},T_{x,j}$ are temperatures of wall volume $j$ and flow volume $j$_, respectively, $\omega$ is the wetted perimeter, and ${\mathit{hc}}_{x,j}$ is the local heat transfer coefficient for convective heat transfer, which is evaluated via the local Nusselt number ${\mathrm{Nu}}_{x,j}$. Thus, for this heat transfer unit $j$, the heat transferred from hot volume $Q_{h,j}$ to cold volume $Q_{c,N-j}$ is balanced at the in-between wall volume,

$$l\cdot {\overline{\rho }}_w\cdot A_w\cdot {\dot{T}}_{w,j}=Q_{h,j}+Q_{c,N-j}$$E6

Lastly, a critical task for the above modeling is to select a proper $\mathit{Nu}$ number correlation to determine the Nusselt numbers Nu and then the local heat transfer coefficients ${\mathit{hc}}_{x,j}$ (Eq. (5)). A comparison of multiple correlation candidates, as well as more discussions for the above models, are presented in the literature [2].

2.2 A first-order approximation for sCO₂-PCHEs

As reported in several published studies, the exit temperature of heat exchangers can be expressed as a first-order approximation [19, 20]. In specific, if one inlet is experiencing a step change in mass flow rate or temperature, the exit temperature is a function of time constant $\tau$ and time delay $t_r$, such as,

$$\left\{\begin{array}{lll}t\le t_r& :& T\left(t\right)=T_0\\ t>t_r& :& \frac{T\left(t\right)-T_{\infty }}{T_0-T_{\infty }}=exp\left(-\frac{t-t_r}{\tau }\right)\end{array}\right.$$E7

Using this approximation, the dynamic behavior of the heat exchanger can be represented as a first-order system described by a transfer function (TF) [17, 21], such as,

$$G{\left(s\right)}_{1\mathit{st}}=\frac{Y\left(s\right)}{F\left(s\right)}=\frac{K}{\mathit{\tau s}+1}exp\left(-t_{r}s\right)$$E8

where $Y\left(s\right),F\left(s\right)$ are the Laplace transformation of system’s response $y\left(t\right)$ and forcing input $f\left(t\right)$, respectively, $\tau$, time delay $t_r$ and gain $K$ are to-be-determined parameters. For sCO₂-PCHEs, the time delay $t_r$ usually equals 0 regarding PCHEs’ short flow path, and the gain $K$ normally equals 1 for the symmetry of two flows, i.e., a step change in one flow would lead to a change of the same scale in the other flow. Therefore, this chapter will demonstrate how to calculate the time constant $\tau$ with great attention.

2.3 Time constant calculation

The time-domain response of the first-order system (Eq. (8)), when taking a step-down input $f\left(t\right)=\mathit{Au}\left(t\right)$, can be described as,

$$T\left(t\right)=T_0{e}^{-t/\tau }+\mathit{KAu}\left(t-t_r\right)\left(1-e^{-\left(t-t_r\right)/\tau }\right)$$E9

where $T_0$ is the initial value of $T\left(t\right)$ for $t\le t_r$. Therefore, we can locate the best guess of time constant through the data-fitting approach. The retrieved time constant $\tau$ as well as gain $K$ and delayed time $t_r$ should fit the system with the experimental or simulation results with minimum errors [22, 23], such as,

$$\underset{\tau ,K,t_r\in R}{argmin}\frac{1}{n}\sum _{i=1}^n{\left[T_i-y_i\right]}^2$$E10

where $T_i=T\left(t=t_i\right)$ is the time-domain response of the system (Eq. (9)) at time $t_i$, $y_i$ is the PCHE’s outlet temperature at the same moment, $i=1,\dots ,n$ indexes $n$ moments of a dynamic experiment or simulation process. In the above optimization, mean square errors ($\mathit{MSE}$s) are employed to evaluate prediction errors, for which maximum absolute error ${\sigma }_{\mathit{max}}$ is an alternative metric, such as,

$${\sigma }_{\mathit{max}}=max\left(|T_i-y_i|\right),i=1,\dots ,n$$E11

3. Validation case

3.1 Set-ups

This section sets up a parameterized, general-purpose Modelica program to simulate sCO₂ PCHEs’ dynamic behavior under various experimental scenarios, with which we can calculate time constant $\tau$ accordingly. The program, as published in our open source project (SCOPE) [2], adopts two source mass flow components and two sink pressure components. By setting the temperature, mass flow rate, and pressure of the sources to the experimental scenarios, the program can simulate the PCHE’s dynamic behavior under various geometrical parameters and operating conditions. The graphical user interface of the simulation program is illustrated in Figure 3.

Figure 3.

Modelica test program GUI.

3.2 First-order system validation

We select the configuration (Table 1) of a HTR of a 10 WMe sCO₂ Recompressed Brayton Cycle [24] as a demonstration. When taking a step-down change in hot inlet mass flow rate, the hot outlet temperature exponentially drops as defined in Eq. (9). The time constant $\tau$ identified by Eq. (10) is $251.92s$, which yields a good match between SCOPE’s simulation output and the first-order system, with an MSE of ${0.6761}^{°}{C}^2$ and a ${\sigma }_{\mathit{max}}={1.4188}^{°}C$ (Figure 4a, Table 2).

Name	Symbol	Value
Hot/Cold inlet temperature	$T_{h,\mathit{\text{in}}}/T_{c,\mathit{\text{in}}}$ (°C)	578.15 / 181.95
Hot/Cold outlet temperature	$T_{h,\mathit{\text{out}}}/T_{c,\mathit{\text{out}}}$ (°C)	192.0 / 532.1
Hot/Cold inlet pressure	$P_{h,\mathit{\text{in}}}/P_{c,\mathit{\text{in}}}$ (bar)	89.62 / 237.48
Hot/Cold mass flow rate	${\dot{m}}_h/{\dot{m}}_c$ (kg/s)	100.17 / 100.17
Weight / Height	$W/H$ (mm)	360 / 3710
Channel length	$L_{\mathit{fp}}$ (mm)	800
Channel diameter	$D_{\mathit{ch}}$ (mm)	1.5
Metal wall thickness	$T_{\mathit{\text{wall}}}$ (mm)	0.51
Channel number	$N_{\mathit{ch}}$ ([1])	700,000

Table 1.

Geometry and operation parameters for time constant evaluation tests.

Figure 4.

Comparison between local property-based simulation and LTI system prediction.

Name case	Standard	Off-Des. case A	Off-Des. case B
${\dot{m}}_{\mathit{hot},\mathit{\text{in}}}$	100.17 kg/s	95 kg/s	105 kg/s
${\tau }_{\mathit{op}}$	251.92	197.66	320.08
$\mathit{MSE}$	0.6761	0.3209	1.2938
${\sigma }_{\mathit{max}}$	1.4189	2.1662	1.8018
${\mathit{MSE}}_{\mathit{std}}$	–	206.7400	200.5752
${\sigma }_{\mathit{max},\mathit{std}}$	–	17.2889	17.3021

Table 2.

Detailed discrepancies between standard case and two off-design cases.

However, the same value of $\tau$ will produce significant errors when applying for other starting mass flow rates (Figure 4b and c). For a below-design-point operation where ${\dot{m}}_{\mathit{hot},\mathit{\text{in}}}$ starts from $95\mathit{kg}/s$ (Figure 4b), the hot outlet temperature predicted by the first-order system, with $\tau =251.92$, diverges from the SCOPE’s result by an $\mathit{MSE}$ of ${206.74}^{°}{C}^2$ and an ${\sigma }_{\mathit{max},\mathit{std}}$ of ${17.2889}^{\circ }C$ at $t=1000s$. Similarly, for an above-design-point operation where ${\dot{m}}_{\mathit{hot},\mathit{\text{in}}}$ starts from $105\mathit{kg}/s$, the discrepancies are $\mathit{MSE}={320.08}^{°}{C}^2$ and ${\sigma }_{\mathit{max},\mathit{std}}={17.3021}^{°}C$ at $t=1000s$. The remarkable discrepancies reflect the drawback of the assumption of constant heat transfer coefficients held by the low-order system approximation. Under the assumption, the time constant is a system instinct property determined by heat exchangers’ geometrical parameters and should remain constant during a system’s transient. However, for sCO₂-PCHEs where the working fluid experiences steep non-linear variations, the time constant’s value varies as the initial operating condition changes. As a consequence, employing the $\tau$ identified at the on-design condition will produce significant prediction discrepancies if the PCHE works at alternative operating conditions (Figure 4b and c). A variable time constant, on the contrary, can achieve minimal prediction errors under various operating conditions.

To prove the above point, we identify the value of time constants per operating condition, which produces good agreements between the first-order system and SCOPE’s outputs. In specific, a time constant $\tau$ valued at 197.66 produces a prediction with $\mathit{MSE}={0.3209}^{°}{C}^2$ and ${\sigma }_{\mathit{max}}={2.1662}^{°}C$ for the below-design-point operation (Figure 4b) and a time constant $\tau$ valued at 320.08 generates an $\mathit{MSE}={1.2938}^{°}{C}^2$ and ${\sigma }_{\mathit{max}}={1.8018}^{°}C$ for the above-design-point case (Figure 4c). Therefore, it is reasonable to introduce a notation named optimal time constant, ${\tau }_{\mathit{op}}$, as a refined key parameter to generate minimal prediction errors for a sCO${}_2$ PCHE applications under varying configurations toward dynamic and control strategic studies.

4. Conclusions

For sCO₂-PCHEs, a variable time constant determined by the system’s operating conditions is favorable to approximating the component’s dynamic behavior with a first-order system. Moreover, PCHEs are critical components for sCO₂ power blocks considering their substantial impacts on the heat transformation and flexibility of the entire power block. A better quantitative understanding of its flexibility through analyzing time constant can effectively serve the investigation of the control strategy of sCO₂ power cycles. Understanding the ${\tau }_{\mathit{op}}$‘s distribution under varying operating conditions provides informative knowledge for designing and optimizing control strategies. For example, it is reasonable to maintain the working fluids at relatively low initial mass flow rates and alter the mass flow rate with the highest possible magnitude (within the permitted limit) to achieve the lowest time constant and highest flexibility. Therefore, for power block control studies, it is essential to take the trade-off between $\tau$ and $\eta$ into account and utilize the advantages to achieve the control objective with the highest speed but lowest cost.