Coupled Heat Transfer and Hydraulic Modeling of an Experimental Printed Circuit Heat Exchanger Using Finite Element Methods

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₂ 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₂ 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.¹ 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.

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₂ 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₂ 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₂ 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₂ 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.² This channel pair required a 132,560 element mesh on a footprint of 196 × 896 mm, or 755,000 elements per m². 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², 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.

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, ϕ₁ = V₁/V. As a solid, this first “stream” only has a temperature DoF, T₁. 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.

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 ϕ_ℓ.

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.

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.

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.

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.

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₂ 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.

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.

Experimental investigations into PCHE performance have been made as a means of validating the HHXT model. A scaled airfoil micro-channel CO₂–CO₂ 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.

The airfoil PCHE consists of five cold CO₂, four hot CO₂, and one instrumented 1.5 mm thick plates. The CO₂ 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₂ is depicted in Fig. 1. The hot side channel is simply the mirror image of this about the vertical.

CO₂ 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³), and c_p = 533 (J/(kg K)).

To match the experimental conditions, the two fluid streams are modeled as CO₂ 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₂ was used to calculate ρ and c_p [18], Vesovic's method was used in determining k [19], and a correlation proposed by Fenghour et al. defined μ [20].

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₂ 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₂ 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 Λ₁₋₄ 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.

In experiment, the airfoil PCHE was designed to operate at CO₂ 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₂ 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.

The airfoil PCHE consists of five cold CO₂, four hot CO₂, and one instrumented 1.5 mm thick plates. The CO₂ 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₂ 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₂ properties, velocity, and heat transfer terms. Convergence to a residual <1 × 10⁻⁶ 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.

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.

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.

The hydraulic solution is shown for both the cold CO₂ and hot CO₂ 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₂ 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₂ 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₂ is only heated and the hot CO₂ 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.

The solution for the PCHE solid body temperature, T₁, 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₂ 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₂ 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₂ 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.

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.

1. The cross-flow friction factor (f_y) 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).

2. 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.

3. 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.

4. 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.

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₂ 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.³

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.

There are no conflicts of interest.

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.⁴ The authors attest that all data for this study are included in the paper.

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

RV1,ℓ =

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

{} =

column vector row vector

[ ] =

matrix

⌈⌋ =

diagonal matrix

| | =

magnitude of vector, absolute value of constant

cond =

anisotropic conductivity matrix

ℓ =

stream index

perm =

anisotropic permeability matrix