## Abstract

The Homogenized Heat Exchanger Thermohydraulic (HHXT) modeling environment has been developed to provide thermodynamic modeling of printed circuit heat exchangers (PCHEs). This finite element approach solves solid conduction and fluid thermohydraulics simultaneously, without the need to mesh the minuscule micro-channels of a PCHE. The model handles PCHE features such as headers, solid side walls, and channel inlet and outlet regions, in addition to the micro-channel core. The HHXT model resolves PCHE thermohydraulics using simple model definitions and minimum computational overhead, making it an ideal design tool. This work introduces the thermohydraulic model at the core of HHXT. The homogenization approach used in the model occupies a medium between simplified linear analyses of heat transfer within a PCHE and the brute force of a fully resolved finite element, or computational fluid dynamics, model. An example problem modeling an experimental PCHE is presented. The ability of the HHXT model to simulate fluid flow through a directional varying micro-channel core of two heat-exchanging streams is demonstrated. The HHXT model is being distributed for free within the research community.

## 1 Introduction

Compact and thermally efficient, printed circuit heat exchangers (PCHEs) are favored for use in next generation nuclear power plants. Containing thousands of small working fluid channels distributed in a solid 316H or 800H block, PCHEs can handle high pressures and operating temperatures required by generation IV nuclear plants. For example, the reference helium-cooled reactor of the Next Generation Nuclear Plant project (NGNP) is designed to operate with a core outlet of 800–950 °C [1]. Advanced nuclear reactors will require the certification of a nuclear service PCHE design by construction codes, such as ASME BPVC Sec. III Div. 5. Compliance with this standard requires stress to be computed from expected loading service transients.

Detailed loading transients in PCHEs are needed so that Sec. III creep–fatigue limits can be determined. A plant dynamics code modeling a sodium fast reactor driven 250 MW thermal sCO_{2} Brayton cycle plant is being used by Argonne National Lab (ANL) to determine plant transients [2,3]. These transients range from trivial load following cases to severe accidents such as a pipe break. BPVC Sec. III HBB 3100 loading classifications can be matched to these Na-sCO_{2} transients, allowing the PCHE's resistance to creep–fatigue to be certified.

Determining these transients in PCHEs requires new analysis methods. For traditional nuclear service components, such as pipping, and shell and tube heat exchangers, the stress-loading relationship can be solved analytically or with a moderately sized finite element (FE) model. In PCHEs, stresses are manifest at the micro-channel scale, while their driving loads must be resolved over the entire heat exchanger. Section III compliant FE modeling of some features can be realized at the micro-channel scale, and over slightly larger sections composing many channels, but is not practical for modeling a full-scale PCHE.

Creep–fatigue and ratcheting analysis requirements in ASME Sec. III necessitate a simplified and flexible thermohydraulic modeling approach that can be run over dozens of transients for multiple heat exchanger geometries. Using traditional FE and computational fluid dynamic (CFD) methods to model PCHE scale thermal transients yields billion degree-of-freedom (DoF) problems that are slow to solve and cumbersome to adapt to varying PCHE designs. An appropriate design tool, that bridges the gap between Sec. III loading service data and feature specific structural models, is needed.

The Heat Exchanger Thermohydraulic (HHXT) model has been developed to provide a full PCHE model needed to properly resolve Sec. III loading conditions without the complexity inherent in resolving all facets of the PCHE geometry.^{1} This is being done by expanding on previous work by Urquiza et al. in homogenizing the PCHE micro-channel core as an anisotropic media [4]. Flow distribution in the PCHE is determined using a porous media approach. Heat transfer between fluids is modeled using heat transfer coefficients and is coupled to a conduction model of the solid PCHE body.

To model thermal transients, the HHXT model is fully transient and three dimensional, and is able to account for the thermal inertia of all relevant features such as the core, walls, headers, and nozzles. Resulting transient temperature distributions in the PCHE metal can be used to calculate transient thermal stresses within the PCHE. The hydraulic, heat transfer, and thermal models all run off a common mesh, resulting in far fewer DoFs than needed in comparable CFD based solvers. This creates a model that solves transients quickly, is flexible in PCHE geometry input, and provides rich thermal transients for establishing ASME Sec. III stress transients.

While designed to analyze PCHEs, the HHXT approach is applicable to any compact heat exchanger design that consists of a singular connected solid body containing separated fluid streams, e.g., plate and fin, but not micro-tube, heat exchangers.

## 2 Review of PCHE Thermohydraulic Analysis Approaches

As with other types of the heat exchanger, it is most elementary and convenient to model the PCHE using the effectiveness number of transfer units (ε-NTU) methods. This is sufficient when the application does not require detail within the PCHE and the design of the PCHE is in a preliminary stage. Sizing and economic analysis of a PCHE design are two applications, run during the preliminary design phase of the PCHE application, that are adequately addressed using ε-NTU methods. For instance, when designing a He to sCO_{2} intermediate heat exchanger to be used in an high-temperature gas reactor, Kim and No used ε-NTU methods, combined with the CFD analysis of micro-channel heat transfer, to find a cost optimum channel design and PCHE size [5]. With sizing completed, the detailed component level design of the PCHE can begin.

Component level analyses of a PCHE require models that more accurately capture thermal performance, and resolve thermal gradients, within the PCHE body. Fluid property variation, channel geometry changes, and the transition between fluid inlets, outlets, and the counter-flow core, all effect the thermal performance in ways not captured by simple ε-NTU methods. Furthermore, they generate thermal gradation and secondary stresses, threatening the mechanical integrity of the PCHE through fatigue and ratcheting. Insufficient attention to detail within the modeling phase can overlook problematic geometry features and loading scenarios. If design margins are tight, in-service failure of the PCHE can result.

Ultimately, it is the job of the designer to ensure that all relevant aspects of PCHE geometry are addressed within the thermal and mechanical analysis. While designs and channel arrangements will vary, the process of etching, bonding, and assembling of a PCHE will result in common features which should be considered.

These features are identified in Fig. 1 for an experimental airfoil micro-channel PCHE that was constructed for CO_{2} recuperative testing at the University of Wisconsin. The features exist in distinct regions of the PCHE (outlined by magenta lines in the figure). Any model of PCHE thermodynamics needs to account for the variation of hydraulics and heat transfer from region to region. Unfortunately, due to varying levels of sophistication across PCHE model approaches, some, or most, of these regions might be neglected. Approaches to PCHE modeling can be split between linearized analyses and homogenization methods. Both approaches use estimations of the pressure drop within, and the heat transfer to the PCHE from, the fluid streams. Modeling approaches can also choose to address conduction and heat storage within the PCHE body, and sometimes will do so for the fluid streams if necessary—i.e., sodium applications.

Use of a friction factor *f*, Darcy or Fanning, is commonly used to express the pressure drop performance of a PCHE micro-channel. Heat transfer coefficients vary between the dimensional local heat transfer coefficient, *h*, and its dimensionless equivalents: Nusselt number, Nu, and Colburn factor, *j*. All these can be determined through experiment, with CFD simulations of the micro-channel, or a combination of both. Kruizenga et al. [6] and Li et al. [7] perform an evaluation of *f* and Nu for Zigzag micro-channels in sCO_{2} using both experiment and CFD, while Chen et al. [8] do the same for straight channels in He.

In linearized analysis, only the PCHE Core and Transition regions are addressed, and their behavior is reduced along the axis of the PCHE. One-dimensional discretization of fluid-to-fluid heat transfer is used, with finite-difference and ad hoc methods being the most popular. Linearized analyses capture the overall performance of the PCHE well and can account for variable fluid properties, e.g., CO_{2} operating near the critical point. Given their simplicity and potential for rapid numerical evaluation, linearized PCHE models have also found use within plant system dynamics models such as that developed at Argonne National Laboratory for the AFR-100 sodium-cooled fast reactor [9,10].

Homogenized analyses address all PCHE regions within a single 2D or 3D model of the PCHE. They do so by homogenizing the behavior of micro-channeled regions within a representative volume. The PCHE thermohydraulic problem is presented as a system of partial differential equations (PDEs) that describe the flow, conduction, heat transfer, etc. that all exist within the element. Regions of the model contain different homogenized volumes, e.g., the core region would be homogenized using a volume that captures both counter-flow streams such as that shown in Fig. 2. The PDEs can be evaluated using finite-difference or FE methods.

So far only a finite-difference homogenization method in 2D has been demonstrated by Urquiza et al. [4]. Urquiza's model could capture entrance, transition, core, exit, and solid PCHE regions. Using a porous media model for the hydraulic equations, the model captured two-dimensional flow. Urquiza's work proved the homogenization method's applicability in PCHE thermohydraulic modeling but was limited in its finite-difference implementation. Being finite difference based, meshes and geometries that could be used were limited to rectilinear space, and the assembly of boundary conditions in regions of non-linearity proved to be unstable [11].

While apparently similar to commercially available fluid-solid FE implementations (fluent, cfx), the homogenization method differs in its evaluation of the fluids and solid simultaneously within every element. Commercially available implementations evaluate each fluid and the solid on separate elements, thereby requiring each fluid channel to be covered by at least one element. Pra et al. [12] provide an example where traditional FE was used in modeling transients in an air-to-air PCHE. Here, the modelers had to suffice with modeling a single pair of hot and cold air channels, instead of the 49 pairs that composed the block.^{2} This channel pair required a 132,560 element mesh on a footprint of 196 × 896 mm, or 755,000 elements per m^{2}. In comparison, the FE implementation of the homogenization presented in this paper was able to use a 2750 element mesh over a 774 × 121 mm footprint, which at roughly 29,000 elements per m^{2}, is a marked improvement.

Computational complexity is reduced when using the homogenization method, creating a more useful design tool than is available in commercially available FE implementations. For example, ansys fluent can utilize thermal (conduction equations) or fluid (porous media or Navier–Stokes) element formulations. In fluent, the solid and fluid regions need to be meshed separately. The homogenization method utilizes a single element formulation that captures both solid conduction and porous media behavior of many fluids. Meshing of PCHE is simpler, as fluid channels do not need to be resolved using separate elements. A full PCHE can be meshed using fewer elements, thus reducing the total number of DoFs needing solution. Furthermore, an improvement in mesh element quality is possible, as small contorted elements representing each fluid, and the solid separating the fluids, within the micro-channel are not needed. This reduced complexity allows the homogenization approach to be better used as a design tool where simple setup and quick solution are necessary. Furthermore, complex analyses that require many solutions of the FE problem, such as transient simulations, can be performed more quickly using homogenized methods.

## 3 Finite Element Approach to PCHE Thermohydraulics

The FE approach used in the HHXT model evaluates a homogenized thermodynamic model of the PCHE structure. Here, the internal micro-channel regions of the PCHE are homogenized within a volume of interest, *V*, as shown in Fig. 2. The modeler must consider the different homogenized regions that would appear in the PCHEs, e.g., two stream, single stream, and no stream (solid) regions.

Within the homogenized volume, the interaction of the separate fluid streams and PCHE body is described using a system of connected non-linear PDEs. The problem can be reduced to a few degrees-of-freedom: the temperature of the core, the pressure of each fluid stream, and the temperature of each fluid stream. With multiple fluid streams, the number of DoFs is increased.

We adopt a stream notation, *ℓ*, where the first stream, *ℓ* = 1, corresponds to the solid PCHE body that occupies a fraction of physical space, *ϕ*_{1} = *V*_{1}/*V*. As a solid, this first “stream” only has a temperature DoF, *T*_{1}. All further streams, *ℓ* = 2, 3,…, correspond to different fluid channels occupying their own fraction of space *ϕ _{ℓ}* =

*V*/

_{ℓ}*V*. Each fluid stream has both pressure and temperature DoFs,

*P*, and

_{ℓ}*T*.

_{ℓ}The system of PDEs solved in the FE model consists of three different PDEs repeated to capture any number of fluid streams, see Table 1. The model contains a PDE for conduction in the solid body, along with PDEs representing the hydraulics and heat transport within each fluid stream. A non-linear Darcy porous media PDE is used for hydraulics along with another incompressible heat advection PDE for each fluid stream.

Stream index | DoF | PDE |
---|---|---|

ℓ = 1 | T_{1} | Solid conduction PDE |

ℓ = 2 | P_{2} | Darcy porous media PDE |

T_{2} | Fluid heat advection PDE | |

ℓ = 3 | P_{3} | Darcy porous media PDE |

T_{3} | Fluid heat advection PDE | |

… | … | … |

Stream index | DoF | PDE |
---|---|---|

ℓ = 1 | T_{1} | Solid conduction PDE |

ℓ = 2 | P_{2} | Darcy porous media PDE |

T_{2} | Fluid heat advection PDE | |

ℓ = 3 | P_{3} | Darcy porous media PDE |

T_{3} | Fluid heat advection PDE | |

… | … | … |

### 3.1 Solid Conduction Partial Differential Equation.

*j*coefficient

### 3.2 Darcy Porous Media Partial Differential Equation.

The Darcy Porous Media PDE provides the hydraulic solution. Instead of modeling flow through each individual micro-channel, the core is approximated as a porous medium of many fluid streams each with a porosity of *ϕ _{ℓ}*.

*k*. For orthotropic media, permeability is a diagonal matrix, $[k]perm$; giving the following form of Darcy's law [13]:

_{perm}**v**, is related to the fluid velocity in the channel

_{D}**v**by

A non-linear form of the Darcy porous media equation is used where the permeability matrix is determined from a prescribed Darcy friction factor, *f*, and is thus a function of velocity and fluid properties. In order to model the directional flow that occurs in micro-channels, the permeability matrix is also anisotropic.

**v**}

_{D}*D*is the internal hydraulic diameter of the micro-channel. This can be expressed either as the ratio of fluid volume to channel surface area, or the ratio of the homogenized volume to surface area—i.e., surface area density

_{h}If flow is laminar (*f* = 64/Re), Eq. (5) will only have one {**v _{D}**} term. Substituting this into Eq. (2) results in constant permeability; thereby returning the original laminar definition of Darcy's law.

**v**}. This is seen by combining Eqs. (2) and (5). The resulting permeability tensor components are expressed in terms of the Darcy friction factor and the pressure gradient

_{D}If the friction factor is not constant, the permeability tensor must be implicitly solved for. This is the case for most micro-channels, where *f* increases at lower velocity. Such variation is best captured with the flow's Reynolds number, Re. To be most practical, the PCHE model accepts input of Re dependent functional definitions of *f* as well as constant and laminar definitions.

### 3.3 Fluid Heat Advection Partial Differential Equation.

*P*solution

### 3.4 Heat Transfer Specification.

*h*, is the easiest to use and allows the volumetric resistance to be written in terms of the micro-channel hydraulic diameter given in Eq. (6)

*j*, or the Nusselt number, Nu, are more useful than the dimensional

*h*

### 3.5 Finite Element Implementation and Stabilization.

The PDEs are assembled using a mix of Galerkin and Streamline Upwind Petrov Galerkin (SUPG) methods, then integrated element-wise using Gaussian Quadrature. The Crank–Nicholson method is used for time integration. Currently only 2D triangular and 3D tetrahedral meshes of linear or quadratic order are supported. Standard Garlekin formulation is used for the Solid Conduction and Darcy Porous Media PDEs (Eqs. (1) and (4)). The Galerkin formulation cannot be used on the Fluid Heat Advection PDE, Eq. (9), as the heat advection term contains first-order temperature dependencies which introduce instability. Evaluating this PDE using an advective method, i.e., the SUPG method, eliminates the instability but introduces some additional complexity.

The SUPG method adds a correction by perturbing elemental approximating functions in the upwind direction of the streamline. The correction is implemented element-wise and depends on the advective term from Eq. (9) and a stabilization parameter *τ*. Many examples of the SUPG method's use in mass transport problems exist. For interested readers, a paper by Cui et al. [14] provides an example of the numerical implementation, while a book by Donea and Huerta [15] details the theory of common stabilization techniques.

The choice of the stabilization parameter, *τ*, influences the quality of the solution and has been subject to extensive research within the FE community, e.g., the review [16]. A general “optimal” definition of *τ* does not exist, resulting in multiple differing and valid interpretations of *t* to be used.

*t*in the PCHE FE model follows the methodology outlined by Knobloch [17], which involves extending a common 1D formulation into 2D and 3D elements. Here, the stabilization parameter is calculated from streamlines of length,

*s*, and Peclet number, Pe, that exist within each element of the mesh

*α*is the thermal diffusivity of the fluid in the streamline direction. Since fluid conductivity can be orthotropic in the PCHE model, the diffusivity along

*s*accounts for the directionality of

**b**

### 3.6 Model Input.

The PCHE model is implemented in matlab and integrated into its existing Partial Differential Equation Toolbox. matlab PDE tools are used in model preparation and post processing. Geometry input, meshing, setting initial and boundary condition definitions, and material property inputs are all familiar to the experienced matlab user. Likewise, postprocessing operations such as mesh-based plotting, solution interpolation, and the formatting of solution result files are similar to other matlab PDE models.

The definition of PCHE properties is not limited to simple cases. All fluid and solid body properties can be position, time, temperature, and pressure dependent. A list of all manipulatable material properties and their allowed dependencies are given in Table 2. This allows modeling of variable property problems, e.g., an sCO_{2} recuperator operating near the critical point. This also presents a way for the user to input external property routines, such as NIST real gas properties, into their heat exchanger model. Finally, there are also no geometric limitations on the PCHE model. Geometries can be 2D or 3D, can contain any number of regions, and can include non-rectilinear features.

PCHE body properties, ℓ = 1 | Units | Spatial dependence^{a} | State dependence^{a} |
---|---|---|---|

Volume fraction, ϕ_{1} | (−) | Region, position | |

Thermal conductivity^{b}, $[k]cond$ | $(WmK)$ | Region, position, time | Temperature |

Density, ρ_{1} | $(kgm3)$ | Region, position, time | Temperature |

Specific heat, $cp1$ | $(JkgK)$ | Region, position, time | Temperature |

Fluid properties, ℓ = 2,3,… | Units | Spatial dependence^{a} | State dependence^{a} |

Volume fraction, ϕ_{ℓ} | (−) | Region, position | |

Hydraulic diameter, $Dh\u2113$ | (m) | Region, position | |

Thermal conductivity^{b}, $[k]cond,\u2113$ | $(WmK)$ | Region, position, time | Temperature, pressure |

Density, ρ_{ℓ} | $(kgm3)$ | Region, position, time | Temperature, pressure |

Specific heat, $cp\u2113$ | $(JkgK)$ | Region, position, time | Temperature, pressure |

Viscosity, μ_{ℓ} | $(kgms)$ | Region, position, time | Temperature, pressure |

Darcy friction factor^{b}, $[f]\u2113$ | (−) | Region, position, time | Reynolds number |

OR permeability^{b}, $[k]perm,\u2113$ | (m^{2}) | Region, position, time | Temperature, pressure |

Colburn heat transfer coefficient, j_{1,ℓ} | (−) | Region, position, time | Reynolds number |

OR local heat transfer coefficient, h_{1,ℓ} | $(Wm2K)$ | Region, position, time | Temperature, pressure |

OR Nusselt number, Nu_{1,ℓ} | (−) | Region, position, time | Reynolds number, Prandtl number |

OR volumetric thermal resistance, $RV1,\u2113$ | $(Km3W)$ | Region, position, time | Temperature, pressure |

PCHE body properties, ℓ = 1 | Units | Spatial dependence^{a} | State dependence^{a} |
---|---|---|---|

Volume fraction, ϕ_{1} | (−) | Region, position | |

Thermal conductivity^{b}, $[k]cond$ | $(WmK)$ | Region, position, time | Temperature |

Density, ρ_{1} | $(kgm3)$ | Region, position, time | Temperature |

Specific heat, $cp1$ | $(JkgK)$ | Region, position, time | Temperature |

Fluid properties, ℓ = 2,3,… | Units | Spatial dependence^{a} | State dependence^{a} |

Volume fraction, ϕ_{ℓ} | (−) | Region, position | |

Hydraulic diameter, $Dh\u2113$ | (m) | Region, position | |

Thermal conductivity^{b}, $[k]cond,\u2113$ | $(WmK)$ | Region, position, time | Temperature, pressure |

Density, ρ_{ℓ} | $(kgm3)$ | Region, position, time | Temperature, pressure |

Specific heat, $cp\u2113$ | $(JkgK)$ | Region, position, time | Temperature, pressure |

Viscosity, μ_{ℓ} | $(kgms)$ | Region, position, time | Temperature, pressure |

Darcy friction factor^{b}, $[f]\u2113$ | (−) | Region, position, time | Reynolds number |

OR permeability^{b}, $[k]perm,\u2113$ | (m^{2}) | Region, position, time | Temperature, pressure |

Colburn heat transfer coefficient, j_{1,ℓ} | (−) | Region, position, time | Reynolds number |

OR local heat transfer coefficient, h_{1,ℓ} | $(Wm2K)$ | Region, position, time | Temperature, pressure |

OR Nusselt number, Nu_{1,ℓ} | (−) | Region, position, time | Reynolds number, Prandtl number |

OR volumetric thermal resistance, $RV1,\u2113$ | $(Km3W)$ | Region, position, time | Temperature, pressure |

Dependencies must be interpreted as having SI units.

It can be input as tensors or constants.

Hydraulic and heat transfer properties of the PCHE core channels are directly input into the model. The PCHE model accepts input of the either a porous media permeability or the Darcy friction factor for hydraulics, and a choice between the Colburn, local dimensional, or Nusselt heat transfer coefficient. These can be position, time, and Reynolds number dependent.

## 4 Modeling an Experimental Airfoil Micro-Channel PCHE

Experimental investigations into PCHE performance have been made as a means of validating the HHXT model. A scaled airfoil micro-channel CO_{2}–CO_{2} recuperative PCHE was made for this purpose. The airfoil PCHE features novel embedded fiber-optic temperature sensors, enabling experimental internal temperatures to be mapped and compared with the modeling effort.

### 4.1 Model Definition.

The airfoil PCHE consists of five cold CO_{2}, four hot CO_{2}, and one instrumented 1.5 mm thick plates. The CO_{2} plates contain a 1 mm deep etched, 94.3 mm wide, airfoil channel, containing 8.1 mm long NACA0020 airfoils spaced at 3.65 mm intervals across the flow and 6.9 mm along the flow direction. The channel pattern for the cold CO_{2} is depicted in Fig. 1. The hot side channel is simply the mirror image of this about the vertical.

CO_{2} flows through the cold and hot streams as shown in Fig. 3. Within the transition region and in parts of the entrance and exit regions, the airfoil channels curve, directing flow from the side-mounted inlets and outlets to the horizontally oriented counter-flow core.

As the airfoil PCHE is constructed out of 316 Stainless Steel, model solid properties are taken as the following constants: *k* = 17.6 (W/(m K)), *ρ* = 8030 (kg/m^{3}), and *c _{p}* = 533 (J/(kg K)).

#### 4.1.1 Flow and Heat Transfer.

To match the experimental conditions, the two fluid streams are modeled as CO_{2} at 6 MPa. Temperature-dependent functions for *k*, *ρ*, *c _{p}*, and

*μ*were input into the PCHE model. Property values were pre-determined for the isobar over a 22–250 °C temperature range as shown in Fig. 4. The Span–Wagner equation of state for CO

_{2}was used to calculate

*ρ*and

*c*[18], Vesovic's method was used in determining

_{p}*k*[19], and a correlation proposed by Fenghour et al. defined

*μ*[20].

*D*of 1.498 mm. Once bonded, the plates formed a PCHE core with a

_{h}*ϕ*of 0.264 on the cold side and 0.211 on the hot side. An average Colburn heat transfer coefficient of

*j =*0.003911 was determined from a 1D analysis of the recuperator's experimental performance. For flows of Re > 2000, the channels were found to have a friction factor fitting the following modified Colebrook equation:

*f*could be obtained through a purpose-designed experiment or a CFD study, results of which did not exist at the time the airfoil PCHE model was made. Thus, the cross-flow was arbitrarily defined as having an

*f*of 10 in the

*y*-direction of the micro-channel core. Thus, the directional friction factor tensor used in the model is

#### 4.1.2 Orthotropic Channels.

The core of the airfoil PCHE favors both flow and conduction along the axis of the recuperator (*x*-direction in all figures). Heat can be conducted in all directions within the airfoil core, but the stacked nature of the etched plates greatly favors conduction in the *x* and *y* while limiting it in *z*. Heat can be conducted within the CO_{2} but can only do so in *x* and *y*, as *z* lies perpendicular to the plane of fluid flow within the channel plates. Conduction in the cold and hot CO_{2} also differs, as there is slightly more cold fluid, with five plates of cold channels and four plates of hot.

The directional friction factor and orthotropic conductivity given in Eqs. (19) and (20) apply to the airfoil PCHE's core (region I in Fig. 5). Entrance, transition, and exit regions (regions III–XII) of the PCHE contain airfoil channels which gradually turn 90 deg. The extent of the airfoil channel's turn is evaluated using transform matrices, [Λ], which apply a rotation about one of four centers (magenta points Λ_{1−4} in Fig. 5).

The volume fraction, conductivity, and directional friction factor are defined for each region within the model. The values and transformations used at each region are given in Table 3. As the FE model accepts region and position variable inputs, the turning nature of the airfoil PCHE micro-channels can be captured.

Model regions | Solid body, ℓ = 1 | Cold CO_{2}, ℓ = 2 | Hot CO_{2}, ℓ = 3 | |||||
---|---|---|---|---|---|---|---|---|

ϕ | $[k]cond$ | ϕ | $[k]cond$ | $[f]$ | ϕ | $[k]cond$ | $[f]$ | |

I | 0.525 | $[k]s$ | 0.2644 | $[k]c$ | $[f]x$ | 0.211 | $[k]h$ | $[f]x$ |

II | 1.000 | $k[I]$ | ||||||

III | 0.789 | $49[\Lambda ]90[k]s+59k[I]$ | 0.211 | $[\Lambda ]90[k]h$ | $[\Lambda ]90[f]x$ | |||

IV | 0.789 | $49[\Lambda ]1[k]s+59k[I]$ | 0.211 | $[\Lambda ]1[k]h$ | $[\Lambda ]1[f]x$ | |||

V | 0.525 | $49[\Lambda ]1[k]s+59[\Lambda ]2[k]s$ | 0.2644 | $[\Lambda ]2[k]c$ | $[\Lambda ]2[f]x$ | 0.211 | $[\Lambda ]1[k]h$ | $[\Lambda ]1[f]x$ |

VI | 0.736 | $59[\Lambda ]2[k]s+49k[I]$ | 0.2644 | $[\Lambda ]2[k]c$ | $[\Lambda ]2[f]x$ | |||

VII | 0.736 | $59[\Lambda ]90[k]s+49k[I]$ | 0.2644 | $[\Lambda ]90[k]c$ | $[\Lambda ]90[f]x$ | |||

VIII | 0.789 | $49[\Lambda ]90[k]s+59k[I]$ | 0.211 | $[\Lambda ]90[k]h$ | $[\Lambda ]90[f]x$ | |||

IX | 0.789 | $49[\Lambda ]3[k]s+59k[I]$ | 0.211 | $[\Lambda ]3[k]h$ | $[\Lambda ]3[f]x$ | |||

X | 0.525 | $49[\Lambda ]3[k]s+59[\Lambda ]4[k]s$ | 0.2644 | $[\Lambda ]4[k]c$ | $[\Lambda ]4[f]x$ | 0.211 | $[\Lambda ]3[k]h$ | $[\Lambda ]3[f]x$ |

XI | 0.736 | $59[\Lambda ]5[k]s+49k[I]$ | 0.2644 | $[\Lambda ]4[k]c$ | $[\Lambda ]4[f]x$ | |||

XII | 0.736 | $59[\Lambda ]90[k]s+49k[I]$ | 0.2644 | $[\Lambda ]90[k]c$ | $[\Lambda ]90[f]x$ |

Model regions | Solid body, ℓ = 1 | Cold CO_{2}, ℓ = 2 | Hot CO_{2}, ℓ = 3 | |||||
---|---|---|---|---|---|---|---|---|

ϕ | $[k]cond$ | ϕ | $[k]cond$ | $[f]$ | ϕ | $[k]cond$ | $[f]$ | |

I | 0.525 | $[k]s$ | 0.2644 | $[k]c$ | $[f]x$ | 0.211 | $[k]h$ | $[f]x$ |

II | 1.000 | $k[I]$ | ||||||

III | 0.789 | $49[\Lambda ]90[k]s+59k[I]$ | 0.211 | $[\Lambda ]90[k]h$ | $[\Lambda ]90[f]x$ | |||

IV | 0.789 | $49[\Lambda ]1[k]s+59k[I]$ | 0.211 | $[\Lambda ]1[k]h$ | $[\Lambda ]1[f]x$ | |||

V | 0.525 | $49[\Lambda ]1[k]s+59[\Lambda ]2[k]s$ | 0.2644 | $[\Lambda ]2[k]c$ | $[\Lambda ]2[f]x$ | 0.211 | $[\Lambda ]1[k]h$ | $[\Lambda ]1[f]x$ |

VI | 0.736 | $59[\Lambda ]2[k]s+49k[I]$ | 0.2644 | $[\Lambda ]2[k]c$ | $[\Lambda ]2[f]x$ | |||

VII | 0.736 | $59[\Lambda ]90[k]s+49k[I]$ | 0.2644 | $[\Lambda ]90[k]c$ | $[\Lambda ]90[f]x$ | |||

VIII | 0.789 | $49[\Lambda ]90[k]s+59k[I]$ | 0.211 | $[\Lambda ]90[k]h$ | $[\Lambda ]90[f]x$ | |||

IX | 0.789 | $49[\Lambda ]3[k]s+59k[I]$ | 0.211 | $[\Lambda ]3[k]h$ | $[\Lambda ]3[f]x$ | |||

X | 0.525 | $49[\Lambda ]3[k]s+59[\Lambda ]4[k]s$ | 0.2644 | $[\Lambda ]4[k]c$ | $[\Lambda ]4[f]x$ | 0.211 | $[\Lambda ]3[k]h$ | $[\Lambda ]3[f]x$ |

XI | 0.736 | $59[\Lambda ]5[k]s+49k[I]$ | 0.2644 | $[\Lambda ]4[k]c$ | $[\Lambda ]4[f]x$ | |||

XII | 0.736 | $59[\Lambda ]90[k]s+49k[I]$ | 0.2644 | $[\Lambda ]90[k]c$ | $[\Lambda ]90[f]x$ |

Note: Transformation matrices [Λ]_{1}, [Λ]_{2}, [Λ]_{3}, and [Λ]_{4} are rotations of angle equal to that formed between the evaluated position and the vertical from the centers shown in Fig. 5. Transform [Λ]_{90} is a 90 deg rotation.

#### 4.1.3 Operating Conditions.

In experiment, the airfoil PCHE was designed to operate at CO_{2} mass flows of 0.08–0.10 (kg/s), achieving a design effectiveness of 91–93.7%. However, when operated at a lower mass flow of 0.05378 (kg/s), the airfoil PCHE effectiveness decreased to 88.9%. At this lower mass flow, the PCHE is over-sized, and heat transfer becomes unbalanced, with the heat being transferred being shifted toward one end of the PCHE and the CO_{2} streams starting to pinch at the outlet. Furthermore, turning cross-flow within the transition between the inlets and the core (regions V & X) shift heat transfer off-axis toward the hot inlet (bottom of the model). The resulting PCHE temperature distribution at steady-state was observed in experiment using a novel method of distributing fiber-optic temperature probes within the PCHE.

The fiber-optic measurements offer a point of comparison for the FE model. While multiple steady-state PCHE temperature distributions were measured, the low mass flow case is chosen for modeling, as its 2D temperature distributions are more non-linear and difficult to capture in a model. Experimental conditions for the low mass flow case are summarized in Table 4. The inlet pressures and temperatures were applied as Dirichlet boundary conditions at the inlets of the FE model. Mass flow was fixed at the inlet and outlets using Neumann boundary conditions.

Cold CO_{2} side | Hot CO_{2} side | |
---|---|---|

P_{inlet} | 5.990 (MPa) | 4.892 (MPa) |

ΔP | 4.087 (kPa) | 7.250 (kPa) |

T_{inlet} | 22.5 (°C) | 202.3 (°C) |

T_{outlet} | 144.0 (°C) | 39.1 (°C) |

$m\u02d9$ | 0.05378 (kg/s) | 0.05378 (kg/s) |

Cold CO_{2} side | Hot CO_{2} side | |
---|---|---|

P_{inlet} | 5.990 (MPa) | 4.892 (MPa) |

ΔP | 4.087 (kPa) | 7.250 (kPa) |

T_{inlet} | 22.5 (°C) | 202.3 (°C) |

T_{outlet} | 144.0 (°C) | 39.1 (°C) |

$m\u02d9$ | 0.05378 (kg/s) | 0.05378 (kg/s) |

### 4.2 Steady-State Solution.

The airfoil PCHE consists of five cold CO_{2}, four hot CO_{2}, and one instrumented 1.5 mm thick plates. The CO_{2} plates contain a 1 mm deep etched, 94.3 mm wide, airfoil channel, containing 8.1 mm long NACA0020 airfoils spaced at 3.65 mm intervals across the flow and 6.9 mm along the flow direction. The channel pattern for the cold CO_{2} is depicted in Fig. 1. The hot side channel is simply the mirror image of this about the vertical.

The FE problem was assembled on a 2750 element, 1520 node, mesh. The mesh is shown, along with the PCHE body temperature, in Fig. 7. Mesh convergence was studied by examining the total heat transferred and the average Reynolds number in the cold and hot streams for various mesh resolutions. This is shown for 884, 1352, 2750, 9722, and 61,158 element meshes in Fig. 6. Results from the 2750 element mesh are shown in this paper. At this resolution the solution is sufficiently converged without being so over meshed as to make printed figures such as Figs. 7 and 8 illegible.

The solution of the steady-state proceeds in an iterative fashion with non-linear components of the FE model reassembled at each iteration. Non-linearity results from variable CO_{2} properties, velocity, and heat transfer terms. Convergence to a residual <1 × 10^{−6} is reached within eight iterations, using 209 s of CPU time. Faster solution times are expected in future implementations, as the FE model is in the design phase and can benefit from the streamlining and pre-compilation of large portions of the toolbox.

#### 4.2.1 Comparison to Experimental Conditions.

Results from the homogenized model compare reasonably with experiment. Experimentally only the inlet and outlet conditions, mass flow, and heat transfer are well known. The model operates using *f* and *j* values that were backed out from a one-dimensional simplification of the recuperator. As a result, some play in the actual values of *f* and *j* is expected. A comparison of the model results and experiment is given in Table 5.

Experiment | Model | |
---|---|---|

T_{C,out} (°C) | 143.96 ± 0.95 | 147.92 |

T_{H,out} (°C) | 39.05 ± 0.29 | 40.58 |

Re_{C} | 11340 ± 230 | 11060 |

Re_{H} | 13780 ± 330 | 13610 |

$m\u02d9C$ (kg/s) | 0.0538 ± 0.0006 | 0.0544 |

$m\u02d9H$ (kg/s) | 0.0538 ± 0.0006 | 0.0550 |

Q_{C} (kW) | 9.60 ± 0.33 | 9.94 |

Q_{H} (kW) | −10.01 ± 0.11 | −9.95 |

ΔP_{C} (Pa) | 4090 ± 920 | 3000 |

ΔP_{H} (Pa) | 7250 ± 230 | 8150 |

Experiment | Model | |
---|---|---|

T_{C,out} (°C) | 143.96 ± 0.95 | 147.92 |

T_{H,out} (°C) | 39.05 ± 0.29 | 40.58 |

Re_{C} | 11340 ± 230 | 11060 |

Re_{H} | 13780 ± 330 | 13610 |

$m\u02d9C$ (kg/s) | 0.0538 ± 0.0006 | 0.0544 |

$m\u02d9H$ (kg/s) | 0.0538 ± 0.0006 | 0.0550 |

Q_{C} (kW) | 9.60 ± 0.33 | 9.94 |

Q_{H} (kW) | −10.01 ± 0.11 | −9.95 |

ΔP_{C} (Pa) | 4090 ± 920 | 3000 |

ΔP_{H} (Pa) | 7250 ± 230 | 8150 |

Model results for the average Reynolds number, mass flow, and total heat transfer fall within, or just on the edge of, experimental values and their uncertainties. Outlet temperatures and pressure drops are close to experimental values, but slightly outside of measured experimental uncertainty. The differences are attributed to a few experimental unknowns with require investigation, and that are discussed at the end of Sec. 4.2.4.

#### 4.2.2 Flow and Stream Heat Transfer.

The hydraulic solution is shown for both the cold CO_{2} and hot CO_{2} sides in Fig. 8. Here, the stream pressure is shown as a colormap, and vectors of the resulting velocity are superimposed. Turning of the flow through the entry, exit, and transition regions can be seen. There is also a degree of flow acceleration at the inside of the turns, as the airfoil channel allows cross-flow, and hence redistribution, to occur.

Temperatures within the CO_{2} streams are shown in Fig. 9(a). The stream temperatures are well separated at the hot side of the PCHE, but after transferring heat gradually approach each other within the PCHE core. The hot and cold stream start to pinch within region V, this can be seen in their similar temperature at left in Fig. 9(a) and the near zero heat transfer at left in Fig. 9(b). This is the behavior expected in a lower mass flow case. There is also temperature variation across the channel, as some streamlines spend longer runs within the transition regions (regions V and X).

The fluid temperatures do not change much within the entrance and exit regions (regions III, IV, VI, VII, VIII, IX, XI, and XII). This is due to the single channeled nature of these regions, where heat transfer is only between one stream and the solid body. To heat or cool the fluid in these entrance and exit regions, heat must first be conducted through the solid body from further away. The lack of much heat transfer in the entrance and exit regions is apparent when the volumetric heating of the streams is examined.

The distribution of heating and cooling within the CO_{2} streams is shown in Fig. 9(b). Here, red in the colormap indicates heating, blue cooling, and green no heating or cooling. It is apparent that the cold CO_{2} is only heated and the hot CO_{2} only cooled. Furthermore, no heating occurs in the hot entrance or cold exit region. An axial shift of more heat transfer occurring near the hot inlet side of the recuperator is apparent. Hot fluid entering the PCHE at the bottom-right must turn into the heat exchanger core. Streamlines at the bottom (−y) of the core spend less time in the transition (region X) and thus enter the core hotter, shifting heat transfer to the bottom of the core.

#### 4.2.3 Body Temperature.

The solution for the PCHE solid body temperature, *T*_{1}, is shown along with the mesh in Fig. 7. Since the model is of a lower mass flow case, the more temperature change occurs at the hot end of the PCHE (right side). Heat transfer starts in the transition region (region X), as the hot CO_{2} first encounters the cold stream. The density of heat being transferred increases, as the hot stream turns into the PCHE core (boundary between region I and X). Here, at the right side of the core, the heat transfer peaks, as can be seen in Fig. 9(b). As the hot CO_{2} continues to flow from right-to-left through the core, the streams reach a pinched condition and heat transfer is gradually reduced.

Thermal gradients within the PCHE solid body are greatest where heat transfer changes from single stream-to-body to two stream-through-body. Even in this lower mass flow case, the high heat transfer coefficient and large surface area density of the airfoil micro-channels create heating of the solid body that dominates any diffusion of temperature through the core structure. Structure temperature largely follows the temperature of the fluid it is interacting with. Thus, the body temperature within the hot CO_{2} entrance region follows the temperature of the hot inlet. When the hot entrance abruptly enters the transition region, the core body equilibrates to a temperature between the hot and cold streams.

The airfoil test recuperator was constructed to exacerbate these thermal gradients between entrance and transition region. In designing a service PCHE, these regions would be reduced as much as possible, although it would be impossible to eliminate them entirely.

#### 4.2.4 Comparison to Fiber-Optic Data.

A comparison of the model's body temperature to that measured with embedded fiber-optic temperature probes proves the model's accuracy. Figure 10 plots model recuperator body temperatures along the position of the eight fiber-optic temperature probes with the fiber's experimental reading. The eight comparisons are staggered by 100 °C to aid in visualization. It should be noted that the dips in measured temperature at the end of the PCHE (right side of plot) are due to the local cooling action of protruding capillaries used to house the fibers. These dips in temperature only exist within the fiber sensors and not within the PCHE itself.

A position-by-position comparison between the measured fiber temperatures and corresponding model temperatures is shown in Fig. 11. The linearity of the comparison is noteworthy, and most values fall within a ±20% error band. Only two lines (blue and gray) fall within the larger ±30% error band. These lines correspond to the bottom most and top most lines in Fig. 10, which are for fibers embedded within at the edges of the recuperator (along the interface between regions I and II). Heat loss through insulation and attachments at these sides are contributing factors at these locations.

A better match between the model and experiment requires a few experimental unknowns to be investigated and incorporated into the model.

The cross-flow friction factor (

*f*) in Eq. (19) needs to be properly evaluated. This could be achieved with either a CFD simulation or purpose-designed experiment looking at flow in the turning transition (regions V, VI, X, and XI)._{y}Heat loss through insulation and nozzle connections needs to be measured and implemented as boundary conditions within the model. Insulation heat loss was not quantified during the experiment. Furthermore, the extent that heat was transferred to or from external attachments such as nozzles and mounting brackets needs investigation.

An expansion of

*f*and*j*testing. This should both increase the range of Re tested for the*j*data and attempt to account for channel turning and cross-flow effects in the calculation of*f*and*j*.The thermal resistance between embedded fiber-optic instrumentation and the recuperator block needs to be quantified. In this analysis, the fiber temperature sensors were assumed to directly measure the solid PCHE temperature.

## 5 Conclusions

A system for modeling the thermohydraulic behavior inside of PCHEs, the HHXT model, has been developed to simplify the modeling of thermal loads for ratcheting and creep–fatigue analyses. The HHXT provides a multi-stream heat transfer solution in 2D or 3D for full-sized heat exchangers without needing to resolve the specific micro-channel geometries. This is achieved through a finite element solution of a system of coupled non-linear PDEs which model the following behavior: core solid thermal conduction, non-linear porous media flow in each fluid stream, thermal advection in each fluid stream, heat transfer between each fluid stream and the core solid, and thermal storage in the core and all fluids.

A 2D counter-flow sCO_{2} recuperator problem is presented as a demonstration of the HHXT model. The model is able to model directional change of the flow within a recuperator core where the cold stream enters and exits at the sides of the PCHE and the hot stream enters at the top and bottom. The HHXT model captures both the counter-flow and cross-flow heat exchange that occurs in this geometry. The model is also shown to correctly resolve the effect of the cross-flow region on the distribution of heat transfer.

The ultimate purpose of the HHXT model is to provide realistic boundary conditions for ratcheting and creep–fatigue models, by determining channel wall temperature transients throughout a full heat exchanger. On their own, mechanical models cannot scale to a full-size heat exchanger, as the resulting mesh size would be impractical. Thus mechanical modeling is limited to capturing specific components of the heat exchanger, be it the internal micro-channel structure, solid wall components, or header attachments. For component mechanical models to capture behavior over a full heat exchanger, results from the HHXT model can provide the spatial variation in channel temperatures. The HHXT model would be used to locate and quantify the largest thermal transients, then a component mechanical model could be used to evaluate the shakedown at these locations.

The HHXT model is being distributed for free within the research community. It has been developed, and is distributed, as a matlab toolbox with functionality familiar to any who have used the matlab PDE toolbox. The current version and supporting documentation can be found.^{3}

## Footnotes

Quantity of channel pairs is inferred: figures in Pra et al. show a 160 mm thick PCHE, assuming a common configuration of 1 mm deep microchannels etched in 1.5 mm thick plates and 1/4in thick top and bottom plates would give 49.1 pairs of plates.

See Note ^{1}.

See Note ^{1}.

## Acknowledgment

This work is funded as part of a DOE Integrated Research Project. Project IRP-17-14227 is focused on quantifying the mechanical integrity of printed circuit heat exchangers using experimental and modeling efforts. This is being done with the goal to build an ASME BPVC nuclear code case for compact heat exchangers.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request. The data and information that support the findings of this article are freely available.^{4} The authors attest that all data for this study are included in the paper.

## Nomenclature

**b**=heat advection coefficient

*f*=Darcy friction factor

*h*=local dimensional heat transfer coefficient

*j*=Colburn heat transfer coefficient

*k*=thermal conductivity

**k**=orthotropic thermal conductivity, or fluid permeability

*s*=streamline length within an element

*t*=time

**v**=channel velocity

**I**=identity matrix

*P*=absolute pressure

*T*=temperature

*V*=volume

*c*=_{p}specific heat

**v**=_{D}Darcy velocity

*D*=_{h}hydraulic diameter

*RV*1,*ℓ*=volumetric thermal resistance ((W m3)/K)

- Nu =
Nusselt number heat transfer coefficient

- Pe =
Peclet number

- Pr =
Prandtl number

- Re =
Reynolds number

*α*=thermal diffusivity

- Λ =
transformation matrix

*μ*=viscosity

*ρ*=density

*ϕ*=volume fraction, porosity

### Brackets

### Subscripts

## References

_{2}in Both Heating and Cooling Modes at Supercritical Pressures

_{2}-to-CO

_{2}Printed Circuit Heat Exchanger Performance Tests

_{2}Brayton Cycle Heat Exchangers