Abstract

Printed circuit heat exchangers (PCHEs) have rapidly gained popularity since being introduced nearly three decades ago, and they are currently widely deployed in the petrochemical and aviation industry. Their compactness, thermohydraulic efficiency, inherent suitability for high temperature/pressure fluids containment, and demonstrated durability are some of the reasons the nuclear industry is seeking to adopt this technology as well. However, the relatively strict nuclear-related regulatory design code, especially when classified as critical to the safety of the reactors, is posing challenges to adopting the technology. From stress analysis point of view, one undesirable feature of PCHEs is their geometrical complexity, which is implied by their multilength-scale features. As a result, a full-scale model of a utility-scale exchanger cannot simply be solved on a computer because meshing such components results in a vast number of degrees-of-freedom. This work seeks to address the challenge of stress analyses to PCHEs by presenting a method to simplify the geometry of PCHE designs. The models proposed by this work can be practically analyzed on a standard computer and provide a path for implementing ASME design rules. The analyses presented herein are divided into five separate investigations. Each is carried out to incrementally simplify the analyzed model by addressing features such as the shapes of the flow passages, the complex distribution of stress in large components, the three-dimensionality of the stress and strain, the thermal stresses caused by thermohydraulic operation observed experimentally and more.

1 Introduction

Diffusion bonded heat exchangers (DBEs) are attractive components for use in the nuclear industry due to their compact size, effectiveness, low pressure drop, and demonstrated durability in the petrochemical industry. Developers of next-generation nuclear power plants (NGNPs) are especially interested in this technology because of the inherent suitability of the microchannel for high pressure (HP) and high temperature pressure containment, and the implied capital cost reduction. Furthermore, in a recent report prepared by Nestell and Sham for the U.S. Department of Energy, diffusion bonded printed circuit heat exchangers (PCHEs), a subset of DBEs, were identified as a particularly suitable technology for NGNP and applications of nonphase changing fluids, such as found in supercritical water, helium Brayton, or the recompression supercritical carbon-dioxide (sCO2) power cycles [1].

However, developers of NGNP who want to use this technology as the primary core heat exchanger [2,3] are facing challenges in following ASME guidelines. The current set of ASME rules for design, fabrication, and inspection of high temperature nuclear use are predicated on traditional geometries and manufacturing technologies, and arguably not well suited to DBEs. In addition, due to the little practical experience of operating DBEs at temperatures, pressures, and with working fluids found in NGNP designs, it is unknown how well will they perform over time.

To address some of the knowledge gaps and the state of the ASME guidelines, the U.S. Department of Energy funded an Integrated Research Project (IRP) that includes team members from five universities, DBE manufacturers, national laboratories, and engineering consulting firms. The goals of the IRP are twofold: to identify issues related to operating PCHEs at prototypical NGNP conditions, and to draft an ASME code case for diffusion-bonded components that includes the mechanical, welding, and in-service inspection design rules [4]. To support these efforts, several laboratory-scale exchangers (small but complete DBEs) with heat duty ranges of 5–50 kW and made of 316 L SS and Incoloy 800H, were purchased from Vacuum Process Engineering, Inc (Sacramento, CA). These exchangers were used for qualification of nondestructive examination techniques, testing with high-temperature working fluids and for producing data required for the validation of numerical models simulating their mechanical and thermohydraulic characteristics.

This work is supported by the IRP in order to demonstrate a method of modeling the stress and strains developing due to operational loads, as required by the current ASME Section III, Division 5 rules covering high temperature reactor components. More specifically, for a component operating at elevated temperature and classified as critical to the safety of a reactor (Class A), ASME guidelines require full structural analysis of the component to determine whether the geometry and loads are acceptable. However, performing structural analysis of a complete PCHE is not practical because the numerous small features and the multilength-scale geometry result in a finite element mesh too complex to be solved directly. This paper describes several steps of simplifications to the geometry of the component that reduce the size of the analyzed model without significantly compromising its accuracy. These simplifications make the application of design by analysis (DBA) to PCHEs feasible.

2 Background

Hundreds of studies have been published in recent years on the topic of PCHEs, mostly focusing on characterizing the pressure drop and heat transfer of different geometries and with different working fluids. The investigated geometries can generally be divided into two categories: continuous and discontinuous (examples are shown in Fig. 1). In the continuous designs, the fluid is distributed into microchannels at the heat exchanger inlet (inlet header) and is prevented from mixing between channels until exiting the heat exchanger (outlet header), where the individual streams merge. The individual channels can be straight or take a tortuous path, but always include a continuous wall separating adjacent channels. In the discontinuous designs, there are no individual channels. Rather, there is a fin structure that enhances fluid mixing and periodically disrupts the formation of boundary layers. Common examples of fin structures include airfoil, S-shaped and staggered rectangular fins [57].

To the best of our knowledge, of the two different fundamental geometrical designs, only the continuous is currently commercially available. Vendors typically apply the ASME Code Section VIII (nonnuclear) to size the heat exchanger. Numerous studies describe the application of the Section VIII rules to DBEs [811]. The Section VIII rules, in short, treat each channel as a standalone pressure vessel and simplify the geometry by treating the semicircular channel cross section as a rectangular stayed plate. This design approach is considered “design by rule” (DBR) because analytic formulations for geometrical sizing are provided in the standard, and numerical simulation is generally not required. The efficacy of this design approach is demonstrated by many components currently in service.

Although used successfully, DBR is not considered in this work. ASME guidelines for components classified as critical to the safety of the reactor and operating at elevated temperatures (found in Section III, Division 5, Subsection HB, Subpart B for Class A components) currently require a full structural analysis of the component. The code and associated Code Cases offer three paths for DBA: elastic analysis, elastic-perfectly plastic (EPP) analysis (available as two code cases [12,13]), and inelastic analysis. Following any of the DBA approaches requires performing structural analysis on an accurate geometrical representation of the component. However, a complete and accurate model of a utility-scale PCHE cannot be solved without a major simplification of the geometry. For example, according to a study published in 2003 by Heatric, a PCHE manufacturer, one of the largest industrial PCHEs deployed was in an offshore gas platform. The heat duty of this unit was only rated at 24.4 MW (relatively small when compared to the units found in utility powerplants), but it still includes millions of microchannels [8]. As shown in Sec. 3.1, the mesh required to numerically model a three-dimensional (3D) unit cell (two channels and 3 mm deep) included approximately 400,000 elements. Thus, the mesh of millions of microchannels, approximately 300 mm long, will likely include hundreds of billions of elements. The desire to perform DBA on PCHEs and the need to simplify their geometry has prompted this investigation.

Because of the challenges associated with modeling the geometrical complexity of DBEs, only a handful of studies have been published on the topic of structural analysis in recent years. Lilo et al. [14] investigated the stresses developing between semicircular channels due to different pressure loads in an isothermal cross section and using plε formulation. Lee and Lee [15] investigated the stress intensity developing as a result of temperature gradient between channels on a two-dimensional (2D) model using plε formulation. Mahajan et al. [11] have shown that maximum creep strain and its location are dependent on the PCHE core size when the side walls are accounted for. This study also investigated the von Mises stress (σvm) in varying sizes of models, and corner tip radii of semicircular channels. Pra et al. [16] used a 2D plane stress plσ model to investigate the transient stresses developing during severe thermal cycling. The study also physically tested a PCHE unit under severe thermal cycling and reported that the tested specimen survived 100 loading cycles without a breach of the pressure boundary. Simanjuntak and Lee [17] investigated the reduction of stress concentration at the channel tip as chemical etching is applied to round the flat face of the semicircle. Under the assumption, the diffusion bonded plates are perfectly aligned, the study has shown a significant reduction in peak stress due to the rounded channel tip.

Of the studies above, only two have modeled the stresses developing as a result of both temperature gradient and pressure loads simultaneously. Both studies used 2D planar formulation to simulate the core of the PCHE: one used plε formulation and the other used plσ. While some authors argue the usage of simplified models leads to accurate simulate of certain regions of actual PCHEs, to the best of our knowledge, no analyses in support of this claim were published. The goal of this work is to close the gap between the simplified models analyzed and the actual PCHE geometry. We arrive at the simplified model by incrementally simplifying the complex 3D model, and discuss the limitation and applicability of each simplification step. The simplified model derived in this work will be used in subsequent studies to perform DBA, in accordance with ASME Section III, Division 5 requirements on an actual PCHE geometry, using the EPP and full inelastic paths.

3 Simplification of the Geometry of the PCHE Core

The following work is broken into five distinct analyses. In Sec. 3.1, the model geometry is investigated to determine whether the thermal and pressure induced stresses experienced by a corrugated heat exchanger can be modeled by a simplified, straight geometry. In Sec. 3.2, the efficacy of reducing from a 3D to a 2D model for stress analysis is examined for three different planar stress–strain relations. In Sec. 3.3, the relative significance of thermally induced stresses induced by thermal gradients in the in-plane and out-of-plane directions is discussed in the context of actual operating conditions observed in an experimental loop. In Sec. 3.4, we compare results of a heat exchanger analysis to those of a subset with greatly reduced channel count. Finally, in Sec. 3.5, the sensitivity of stresses to channel misalignment, which is commonly observed in real components, is investigated.

All structural calculations in this work were performed with MOOSE (Multiphysics Object Oriented Simulation Environment), an open-source finite element code developed by Idaho National Laboratory. The mesh, boundary conditions, constraints, and formulations are discussed in detail in each subsequent section. Realistic temperature-coupled thermophysical and mechanical properties of 316H SS were assigned to all models that were loaded thermally. These properties were obtained from ASME Section II, Subpart D (for UNS S31609 and austenitic stainless steel group 3A when applicable). Other than temperature dependency, all analyses in this section assumed the following: small strain, small displacement, linear stress–strain relations, and stress-free temperature (SFT) at 25 °C. This work also assumes the diffusion bonded material maintains the characteristics of the base metal. Specifically, the material remains approximately homogeneous, isotropic, and maintains the temperature dependency shown in Fig. 2. The authors acknowledge that some of the base metal properties may be altered in the diffusion bonded material, but more research is needed to investigate the effects of the real diffusion bonded properties on the models presented in this study.

3.1 Influence of Corrugation on the Stress.

For the analyses in this section, two corrugated models were constructed to cover the range of exchangers found in the literature: a model with a channel corrugation angle of 18 deg (3DZZ18) and corrugation angle of 37 deg (3DZZ37). 3DZZ18 was created to represent an exchanger with a moderate channel corrugation angles that is often found in the literature [18,19], whereas 3DZZ37 was construct based on the geometry of a PCHE that operated at Georgia Institute of Technology (GT). The exchanger tested at GT was manufactured in accordance with DBR depicted in ASME BPVC Sec-VIII to meet the pressure boundary sizing requirements when subjected to the operating conditions shown in Table 1. The PCHE was made of dual certified 316/316 L, included a total of 17 etched layers, 8 for helium (hot-side), and 9 for sCO2 (cold-side) in an alternating configuration and shown in Fig. 3. Geometrical reconstruction of this exchanger was performed by using a high-resolution (90 μm) 3D computed tomography (CT) system. The tomography scans were performed at the University of Wisconsin-Madison in accordance with the ASTM 1695-E CT resolution qualification procedure. More details on the method for geometrical reconstruction were documented in Ref. [20]. The channel features relevant for the analyses presented by this work are shown in Table 2.

To investigate the effect the flow-path corrugation has on the stress developing in the core region, three models were created. Each model simulated a PCHE made of two diffusion bonded plates, each with one channel in the middle. The reference models: 3DZZ37 and 3DZZ18 has a flow path either mimicking the geometry obtained from the X-ray reconstruction (shown in Fig. 3) or close to the geometries found in other studies. The simplified model (3DS) has straight channels constructed with dimensions obtained from a cross section of both reference models taken perpendicular to the flow direction and has the channel dimensions shown in Table 2. The geometries of all models were created in CAD and imported for meshing with ten-noded tetrahedral elements in coreformcubit, a commercial meshing software. Special attention was given to the stress risers near the sharp channel tips by refining the mesh. Each mesh ultimately included more than 400,000 quadratic hex-20 elements in both models.

The pressure and temperature loads shown in Table 3 were applied to the walls of the semicircular flow-paths in the model, where the HP channel was the lower of the two. The models (shown in Fig. 4) were constrained with one plane of symmetry on the face perpendicular to the overall flow direction to simulate the continuous core of the PCHE. Two additional nodal constraints were assigned to prevent rigid body motion.

Two lines of interest were created on both corrugated models for the extraction of high-resolution discretized σvm as shown in Fig. 4 on the top. The lines marked as “Front” were located on the front plane of symmetry while the lines marked as “Center” were at the center of the body and perpendicular to the local flow directions. The stress in 3DS was uniform through its thickness, therefore only one location for extraction of data was used. All lines on the three models were oriented horizontally and were located on the upper wall of the pressurized channel, where the highest stresses in the model are observed (see Fig. 5 on the top).

Von Mises stress was sampled at 1000 equally spaced locations along all lines for extraction of von Mises data and composed into the curves shown in Fig. 5 at the bottom. These data show good agreement between all lines and locations. The small variation that is observed is attributed to the variation in the mesh that could not be perfectly uniform between the two models because of the use of an automatic mesher and use of tetrahedral elements. As expected, a sharp increase in σvm is observed at the extreme ends of all curves, and associated with the vicinity to the singularity caused by the geometrical discontinuity at the channel tips. All curves show a similar rate of increase and overall magnitude of stress.

3.2 Efficacy of Two-Dimensional Models for PCHE Simulation.

In this section, we compare the σvm developing on a 2D model, which was solved with three different planar formulations, and compare the results with a reference 3D model representing the core of a freely supported 3D PCHE (see Fig. 6). The meshes of both models were refined in areas where high stress gradient was expected and the material properties depicted in Fig. 2 were assigned. The 3D model was constrained by three planes of symmetry, to simulate the center of a long 4 ×4 channel PCHE, and two planes of symmetry were used for the 2D models. The 2D model was solved with three following planar formulations: plε,plσ, and generalized Plane Strain (Gplε) formulation. The generalization of the plane strain formulation was introduced by Cheng et al. [21]. In the formulation, the strain tensor components are assumed to be functions of only the in-plane coordinates x and y
ε=ε(x,y)
(1)
Under the small-strain assumption, the strain tensor components are related to the displacement field as
ε=12(u+uT)
(2)
Equation (2) has the following 3D solution
u=[u(x,y)az22v(x,y)bz22(ax+by+c)z]
(3)
In the equation above, a, b, and c are constant scalars. In this study, a and b are set to 0 such that the strain tensor is limited to the form presented in the following equation:
εij=[εxx(x,y)εxy(x,y)0εyx(x,y)εyy(x,y)000εzz(scalar)]
(4)

To investigate the mechanical response of the models to both pressure and temperature loads, the loading conditions in Table 3 were applied to all models. The analysis uses a perfect plasticity constitutive model to approximate the time-independent inelastic response of the material. Von Mises stress data were extracted along three vertical lines as shown in Fig. 6. Figure 7 shows excellent agreement between the Gplε and the reference 3D model. This plot also shows that the plσ and plε formulations differ significantly from the full 3D results. Since Gplε produces such a good agreement, in this work, it is assumed that models solved using this formulation provide a good approximation of freely supported models of PCHE. Additional discussion on this topic is given in Sec. 4.

The temperature load on the 3D model resulted in a uniform temperature distribution along ẑ direction. Since neither of the 2D models can compute stresses caused by “out-of-plane” temperature distribution, that may lead to additional stresses, the use of Gplε should be limited to cases that the thermal stresses arising from streamwise temperature variations can either be neglected or superimposed on the model. In the next set of analyses, we consider whether the streamwise temperature gradient plays an important role in the overall stresses developing under normal operating conditions.

3.3 Stresses Induced by Cross-Channel and Streamwise Temperature Distribution.

We first consider the range of temperature fields that are expected to develop in PCHEs during normal operation using thermohydraulic testing data obtained from a physical test loop. The measured temperature profiles are then divided into temperature gradients: in the flow (streamwise) direction and the transverse (cross-channel) plane, which are aligned with ẑ and x̂,ŷ directions in our models, respectively. A limiting scenario is then selected from the experimental data and applied to a number of different models, to compare in isolation the magnitude of the resulting von Mises stress due to the thermal stresses in the cross-channel and streamwise directions.

3.3.1 Experimental Data.

A 316 L counterflow PCHE that is shown in Fig. 3 was tested between high temperature helium and sCO2 loops at GT shown in Fig. 8. Eight instrumentation ports were distributed along the exchanger core and intercepted the flow of both fluid circuits. These ports allowed a direct fluid temperature and pressure measurements during operation. Third-order polynomial functions were fitted to the helium and sCO2 temperature and pressure measurements (noted as T3(z) and P3(z)), respectively, and used in conjunction with gas property data available from the NIST REFPROP database [22] to derive continuous functions representing the fluid's properties, and heat transfer characteristics. Details on the experimental system, the PCHE instrumentation, methodology of derivation of thermohydraulic figures of merit, and range of operating conditions were published in Ref. [20].

In this work, using the calculated coefficient of convection of the different fluids, and a simplified heat exchanger model with an average wall thickness δ = 1 mm separating the channels of different streams, the temperature gradients through the steel were computed for all operating conditions as shown in Fig. 9. The conditions analyzed here were limited to those where the fluid heat capacitances were well balanced between different fluid circuits, and the system operated at a steady-state. Only normal operating conditions are analyzed in this work because they represent the overwhelming majority time of the operation of the exchanger in service. In all analyzed cases, helium was used as the hot fluid and entered the exchanger at approximately 550 °C and 3.3 MPa, and the sCO2 was used as the cold fluid, and entered at approximately 250 or 350 °C and 20 MPa.

The limiting case was chosen out the range of operating conditions as the case with the highest local thermal gradient in the ẑ direction. This case was chosen because the high thermal gradient is expected to yield the highest thermal stress, thus show the highest possible discrepancy between the simplified model and a full 3D model, discussed in Sec. 3.2. This scenario is highlighted in red on Fig. 9 and has the following approximated streamwise functional temperature profile: T3(z)=532453Z281z21723z3. For this limiting case, the mean difference in wall temperature between the hot and cold channels wall was calculated as 7.51 °C.

3.3.2 Stresses Associated With the Different Temperature Gradients.

Four models were created to investigate the effects the limiting temperature distribution has on the σvm developing in the exchanger. A 2D model (2DTG) with semicircular channels representing the thermal stresses in the cross-channel direction, and three simplified 3D beam models representing the thermal stresses developing due to the streamwise temperature distribution. The cross-channel temperature difference of 7.5 °C was applied to the channel walls in 2DTG, which was solved with Gplε to represent the cross-channel stress at the hot end of the PCHE core. All 3D models had a length of 306.4 mm, simulating the length of the core of the exchanger shown in Fig. 3. These models had a square cross-sectional area that was varied between the three models to represent different length-to-cross section aspect ratios. The 3D model notations used are R1/1 for the square block, R5/1 for a model with streamwise length five times larger than the side of its square cross section, and R10/1 for a model with streamwise length ten times larger than the side of its square cross section. The model's dimensions and mesh are shown in Fig. 10. The 3D models were constrained by two planes of symmetry, which were used to decrease the size of the analyzed models and an additional nodal constraint to prevent rigid body motion.

The solution of the four analyzed models is shown in Fig. 11. The top two plots show the temperature distribution and corresponding σvm on 2DTG. The maximum stress in the model is observed on the channel walls and equals approximately 18 MPa. The bottom right plot shows the streamwise temperature profile that resulted in a temperature difference of approximately 230 °C across the simulated exchanger core. The maximum temperature gradient is approximately 1110 °C/m and takes place near the cold end of the core. The image at the bottom left shows, qualitatively, that the location of maximum thermal stress in all 3D model is at the corner furthest away from the central axis of the exchanger (noted as ẑ by the triad). The bottom left image also shows the location for extraction of stress data which is plotted in Fig. 12.

Results extracted from the dotted lines on the 3D models are shown in Fig. 12. A logarithmic scale was used for the stress axis in order to allow a comparison of the stress data side by side. A significant difference of both the magnitude and location of maximum stress is observed between the three models. In R10/1, the maximum stress is approximately 0.5 MPa, and observed at the location of maximum temperature gradient, which is at the cold end of the core. As the aspect ratio of the beam approaches the cubical shape of 1/1, maximum stress increases up to approximately 20 MPa. The location of the maximum stress shifts away from the location of the maximum temperature gradient and approaches the center of the beam. It should be noted that when the shape of the beam is approaches a cube, the maximal streamwise thermal stress, for the specific temperature distribution analyzed herein, is approximately three times greater than that of the cross-channel thermal stress.

3.4 Number of Channels Needed to Be Included in a Design by Analysis.

The use of 2D Gplε models to simulate the PCHE leads to a significant reduction in the computational efforts required to simulate a complete model. However, when applied to a utility-scale PCHE that is expected to have millions of microchannels, will still result in a model requiring the use of a super-computer. For this reason, studies published on the topic of PCHEs either focused on small lab-scale PCHE or simply modeled unit cells, but without discussing the relevance of the simulation to an entire exchanger that includes a large number of channels confined by a solid envelope. This investigation was performed with two goals in mind: demonstrate the load distribution within the cross section of a PCHE that includes side walls, and to propose a method to simplify the geometry of the cross section for the purpose of performing DBA.

3.4.1 Distribution of Stress in the Core of PCHEs.

The model analyzed in this section (51X50) simulates the cross section of a 60 kW heat-duty exchanger with 6 mm solid envelope as shown in Fig. 13. The size of 51 plates with 50 channels each was chosen because it was the largest model that could undergo a complete thermo-elastic simulation with a sufficiently fine mesh on a workstation with 48 GB of RAM. Model 51X50 was constructed of a single meshed unit cell that was replicated and merged to form a larger model with a perfectly periodic mesh. The dimensions of this cell are identical to the one shown in Fig. 6 and included approximately 1370 elements. 51X50 included a total of 1275 cells, 650 LP, and 625 HP cells and ultimately included nearly 2 × 106 elements. The loading conditions applied to the two types of cells are shown in Table 3 and a symmetry plane was used to constrain and decrease the size of the analyzed model.

The goal of this analysis is to demonstrate the variation in stress in the largest possible model that includes a solid envelope. As noted above, the model included only two types of cells depending on the boundary conditions: LP and HP. Taking advantage of the periodicity of the geometry and mesh, and using a postprocessing algorithm, the solved model was decomposed into the two types of cells. The coordinates of the nodes were then translated and stacked such that the nodes from the different cells shared the exact same location. Aligning all the nodes of each type of cell allowed a direct and detailed comparison between the stresses developing in cells of different regions in the model. Plots of the σvm of all data-points and the two types of cells are shown in Fig. 14 on the top. The data-points shown in black were obtained from the reference cells (5a and 5b), which were located adjacent to the plane of symmetry, and the center of the exchanger. To demonstrate the error, one should expect by analyzing a single cell to represent the σvm in the entire model, node by node data from the reference cells was subtracted from all other cells of each type, and plotted in Fig. 14 on the bottom.

Statistical analysis of the 650 LP cells (shown in Fig. 14 on the top right) shows the mean stress through the cell is 18.1 MPa with a standard deviation of 7.08 MPa. Analysis on the difference between the reference 5a cell and all other LP cells show a maximum difference in σvm of 28.4 MPa and a minimum difference of −16.5 MPa. The mean absolute difference of 5a and the rest of the LP cells is 3.86 MPa. A similar analysis on the HP cells shows the mean stress is 33.1 MPa with a standard deviation of 17.9 MPa. Analysis of the difference between 5b cell and the other 624 HP cells is shown and −23.9 and 9.6 are the minimum and maximum differences, respectively. The mean absolute difference between 5b and the rest of the HP cells is 3.89 MPa.

As noted in Sec. 2, the goal of this work is to produce a method to perform DBA on PCHEs in accordance with ASME Section III, Division 5 requirements to ensure proper sizing based on the expected operational loads. Because of the variation in stress experienced by cells at different location, no single and simplified model can perfectly simulate the complete spectrum of stresses shown above. The combination of the ASME requirement to perform a detailed DBA, the complexity of the stress distribution in PCHEs, and the size of the mesh required to model the cross section of a utility-scale PCHE present challenges that this work cannot overcome. Therefore, instead of accurately modeling the distribution of stress in the model, we seek to bound the stress.

3.4.2 Bounding the Stress in the Mode.

The use of periodic mesh facilitated the identification the data-points with the highest stress value for each particular node location. These data were used for the creation of new dataset, representing the bounding worst scenario of stress for 51X50. This bounding case dataset, along with all other points analyzed in this section, is plotted in Fig. 15. The points highlighted in blue represent the new bounding scenario dataset, which was sampled from all other cells shown in black. Although the bounding bounding case dataset is an artificial dataset that may include stress results from many different cells of different regions, the resulting dataset is smooth and continuous across the HP and LP cells, which is required by specific ASME acceptance criteria (notice the continuity of the data on both sides of the line Y =0).

Two additional models were created to compare bounding worst scenario experienced by models of different sizes. 17X16 included 17 plates with 16 channels each and a solid envelope of 1.92 mm, and 7X6 included and 7 plates with six channels each and a solid envelope of 0.72 mm. The thickness of the envelope was adjusted proportionally to the 6 mm envelope of the already analyzed 51X50, based on the linear number of cells in each model. 17X16 and 7X6 were constrained and loaded with identical boundary conditions and HP and LP cell configuration as 51X50. Similarly to 51X50, the solved models were divided into HP and LP cells, the nodes of each type of cell were translated and stacked, and a dataset of the bounding worst scenario was created for all models, and is shown in Fig. 16.

Statistical analysis of the difference between the bounding stresses obtained from the three models shows the mean absolute difference in stress between 51X50 and 17X16 is 0.8 MPa with a standard deviation of 1 MPa. Of all 2738 analyzed data-points, the maximum difference between the two models is observed near the singularity at the channel corner tip of the HP cell, where model 17X16 predicts a stress 7.2 MPa lower than 51X50. Similar analysis between 51X50 and 7X6 show an excellent agreement as well. The mean absolute difference between the two datasets is 0.74 MPa with a standard deviation of 2.6 MPa. The maximum difference between the two models, again, is observed near the point of singularity at the tip where model 7X6 predicts a stress 9.2 MPa lower than 51X50.

3.5 Sensitivity of the Stress to Channel Vertical Misalignment.

In this analysis, we investigate the sensitivity of the stress to misalignment of the PCHE channels. To consider the range of channel arrangement possibilities in a real PCHE, the unit shown in Fig. 3 was milled using a wire electrical discharge machining. A total of 20 specimens were extracted from the cross section of the PCHE and scanned using a high-resolution digital microscope. Example of one scan is shown in Fig. 17 and demonstrates that channels may be stacked in or out of vertical alignment within one exchanger. For this particular unit, the commonly made assumption of perfectly align HP and LP channels is inaccurate.

Three models were created to simulate the range of offset possibilities. Model 2X20 is the reference model with perfectly vertically aligned channels. The model includes four channels with the dimensions shown in Fig. 10 on the right. Model 2X21/3 and 2X21/2 are the two investigated models where the channels were offset by 1/3 and 1/2 of the cell width, respectively. All models included a relatively thick envelope as the comparison in stress distribution is focused on the ligaments between the channels. The models were solved with Gplε formulation, constrained to prevent rigid body motion and loaded similarly to all other models in this work (see Table 1). The three models are shown in Fig. 18.

The sensitivity of the models to vertical misalignment was compared using the von Mises stress. Data were extracted along a line, passing through areas of high stress in all models and plotted in Fig. 19 at the bottom. Comparison between the three models shows that 2X21/3 and 2X21/2 are experiencing higher peak stress than 2X20 by approximately 20%. However, the model with the vertically aligned channels experiences a highest stress over the center of the channel curve locations 0.7–1.6 and 3–4. When comparing the stress developing across the ligament connecting the two high pressure channel (curve locations 1.9–2.5), the models with the vertically misaligned channel are experiencing a significantly higher stress than 2X20.

4 Discussions

This study performed five analyses to investigate if simplified models can be used to simulate the stress experienced by 3D PCHE geometries under thermohydraulic loads. Pressure and temperature loads obtained from experimental observations on NGNP-relevant fluids were applied to the analyzed models. The models were constrained to simulate a freely supported body in order to focus on the elastic and thermal stresses caused by the interaction between the pressure boundary and the working fluids. Stresses caused by nozzle connection, restraints that may prevent the exchanger from thermally expanding and other external factors were not addressed. These system-specific loads would need to be included in the design analysis of an operating exchanger, but cannot be addressed in this general study. In addition, the channel profile used for representing the flow-path cross section was idealized as semi-elliptical and assumed perfectly uniform through the cross section of the exchanger. The authors acknowledge that for an accurate assessment of fatigue characteristics of such exchanger, a more accurate representation of the channel profile and channel-tip corner radii will be needed. Furthermore, the assumption of perfect uniformity of the channel profile should be re-evaluated when the simplification methodology is applied to the utility-scale PCHE.

To ensure the numerical models analyzed by this work converged properly, specific convergence studies were performed on the models presented in Secs. 3.1, 3.2, and 3.5. These studies indicated that the results reported herein are independent of the mesh, but omitted from this paper for conciseness. It should be noted that numerical convergence at the vicinity to the channel tip, that was modeled as perfectly sharp, could not be achieved due to the geometrical discontinuity and associated numerical singularity. Further studies should focus on accurately modeling the initiation of damage at these stress concentrations. Models developed in this work will be used in subsequent studies for performing DBA on PCHEs operating under a variety of conditions. The methodology for simplifying the geometry presented herein will be used for exploring other PCHE and DBE geometries that can be simplified for the application of the DBA design rules. The discussion below is divided based on the analyses presented in Sec. 3.

4.1 Influence of Corrugation on Stress and Sensitivity to Plate Miss-Alignment.

Analyses presented in this work have demonstrated that the corrugation of the exchanger can be represented using a simplified straight channel for the purpose of stress analysis when both pressure and temperature loads are present. The reference model in this investigation was an exchanger with corrugated channels based on the channel geometry of the component tested at GT, which had a relatively steep flow-path zigzag angle of 37 deg. The cross section of the channel was uniform when measured perpendicular to the local flow direction. The simplified model had a straight channel that was constructed based on the cross section of the corrugated channel taken perpendicular to the flow direction. Von Mises stress was extracted from both models, above the high pressure channel, where the highest stresses were observed. The trends and magnitude of the von Mises stress from the simplified model mimicked closely those obtained from the reference model. Because good agreement was observed between the two models when solved with realistic material properties and using realistic loads, we argue the simplified model can be used for the simulation of the stress and strain experienced by the freely supported corrugated PCHE.

4.2 Two-Dimensional Planar Formulation.

The analysis comparing the results obtained from three planar formulations to a full 3D model under pressure and thermal loads has shown that the plε is significantly overpredicting the von Mises stress when compared with the full 3D model. The large regions of plasticity in the model can be explained by the selection of the SFT at 25 °C and the artificial restriction of the model from thermally expanding in the normal direction to the plane of analysis. Although a more intelligent selection of SFT, such as 537 °C for the particular loading scenario used in the investigation could mitigate the plasticity in the model, the same artificial mechanism of restraining the model from expanding would still exist and induce artificial stresses. As for the plσ model, as expected, it underpredicted the σvm stress as all stress components associated with the out of plane direction are completely neglected. The only formulation that accurately simulated a freely supported 3D model subjected to thermal and pressure loads was the Gplε. Results obtained from the Gplε model were in excellent agreement with the 3D model of an exchanger with straight channels.

4.3 Stresses Induced by Cross-Channel and Streamwise Temperature Distribution.

Four models were created in this section to analyze the magnitude of the thermal stress caused by realistic temperature distribution obtained from experiments. One 2D model was used to investigate the cross-channel stresses induced by the temperature difference between channels of different fluid circuits. Three 3D models were used to investigate the thermal stresses induced by the streamwise temperature distribution on different exchanger geometries. The 2D model was thermally loaded by channel temperature difference of ΔT=7.5 °C and solved using Gplε formulation. This load resulted in large regions of the model experiencing a von Mises stress of 16 MPa, and stress in excess of 20 MPa was observed in the vicinity of the sharp channel tips. The maximum von Mises stress observed on the 3D models varied significantly with the thickness of the model in the transverse to streamwise directions. Geometries simulating narrow (lab scale) heat exchangers have shown that the streamwise temperature gradient may lead to maximum von Mises stresses of less than 1–6 MPa, depending on the thickness of the exchanger, much less than the thermal stresses arising from cross-stream gradients. Models with cubical shape, simulating the utility scale exchanger, will experience much larger stresses of approximately 55 MPa, approximately twice the cross-stream thermal stresses. It should be noted that the maximum von Mises stress observed in the cubical model is approximately 1/2 the yield strength for 316H at the range of temperatures used in the analyses.

The results presented above should be considered bounding the stresses during normal operation for the following reasons: the particular streamwise temperature profile applied to the 3D models was selected because it showed the highest local thermal gradient of all analyzed experimental data. Other experimental temperature profiles had a lower maximum local thermal gradient and thus are expected to lead to lower maximum thermal stresses. In addition, the nonporous solid blocks used to simulate the exchanger in the 3D models are expected to yield higher thermal stresses than a real exchanger with a significantly softened core due to the presence of the channels. To demonstrate the influence, the softened core has over the maximum local stresses developing due to the particular temperature profile, one more model was created. This model included a cavity, rather than a solid core, simulating the limiting case of low stiffness in the highly porous core in the cubical model. The model had an identical length and width as R1/1, loaded with identical temperature profile, was also constrained by two planes of symmetry, and included only a 6-mm envelope simulating the side walls of a highly porous PCHE. The softened model results show a decrease of approximately 50% in the peak von Mises stress, which developed approximately at the same location as in R1/1. The softened model mesh, von Mises plot, and a comparison with R1/1 are shown in Fig. 20.

These results suggest a Gplε model for the thermal stress may not be sufficient to accurately resolve the critical thermal stresses in a full-scale exchanger, though a Gplε model is likely sufficient for the lab scale exchanger. The R1/1 solid analysis is bounding, as the presence of the channels will reduce the effective stiffness of the exchanger, leading to reduced thermal stresses. A design model superimposing a Gplε analysis for the cross-flow thermal stresses with a simplified 3D analysis, potentially accounting for the effective stiffness reduction caused by the channels, may be feasible, at least for design by elastic analysis where stresses can be superimposed.

4.4 Number of Channels Needed to Be Included in the Analysis.

This work aimed at analyzing the largest possible PCHE cross section with the goal of investigating whether a perfectly periodic state of stress develops in the utility scale PCHE, as hypothesized by other studies. The investigated model included 51 plates with 50 channels each and was reduced by half using a vertical plane of symmetry. The model was loaded by both internal pressure and temperature and solved using Gplε formulation. Due to the size of the model, using an engineering workstation with 48 GB of RAM memory was required, and it took approximately 1 h to solve. Analysis on the results presented in Fig. 14 shows an average variation in the von Mises stress, depending on the location of the cell, of approximately 20% within the LP cells and approximately 10% within the HP cells. This variation is likely caused by the interaction between the pressurized porous core and the solid envelope, which is expected to be present in both lab-scale and utility-scale PCHEs. Because of the variation observed in the large model, no simplified model with periodic boundary conditions can simulate the entire range of load possibilities observed in the large model.

As noted in Sec. 2, the goal of this study is to produce a method for performing DBA on PCHEs using the approved design rules for safety-critical components operating at high temperatures in nuclear reactors. For the purpose of following the design rules, a model that bounds the damage may be used to ensure proper sizing of the exchanger. In Sec. 3.4.2, a method for generating a bounding case for the von Mises stress was demonstrated. The method requires generating a PCHE model with a perfectly periodic mesh representing the core of the exchanger. This was completed by meshing a single cell then replicating and offsetting the mesh to form any size of larger model. A solid envelope was then added to the core and the model was then solved and geometrically deconstructed back into the individual cells for postprocessing. A filtering code was then applied to all the cells, to identify the highest local von Mises stress for every given node in all replicated cells. The highest von Mises stress data were used to form a new dataset representing the worst scenario state of stress. This upper bounding state of stress dataset was created to three models of different sizes and has shown to be in a very good agreement for all three. It should be noted that the thickness of the envelope in the reduced size models was decreased proportionally to the model's linear thickness. In the authors' opinion, the method of generating a worst scenario state of stress on a small model allows to effectively ensure proper sizing of any size of heat exchanger, under the limitations discussed in Sec. 3.3.

4.5 Sensitivity to Plate Misalignment.

Optical scans of specimens obtained from a real PCHE component have shown that the common assumption of vertically aligned channels is not necessarily accurate. The channel arrangement shown in Fig. 17 suggests that within one exchanger, multiple possibilities of channel arrangement are possible. Although the arrangement shown in the image could have been intentional by the manufacturer, in the authors' opinion, the misalignment could have been caused by other factors, such as relative movement of the chemically etched plates during diffusion bonding or low tolerance of the chemically etched plates. The goal of the analysis was to investigate whether the different possible channel arrangements may lead to a significant variation in stress under identical boundary conditions. As shown in Fig. 18, three models were created to simulate the range of channel alignment possibilities. Comparison of results between the three models shows an increase of approximately 20% in peak stress in the misaligned models. In addition, regions of the ligament connecting between the two high pressure channels (curve location 1.9–2.5 in Fig. 19) show approximately a twofold increase in the von Mises stress in the nonvertically aligned models. The model with the perfectly aligned channels showed a higher stress only in locations close to the center of the channel (curve location 0.7–1.6 and 3–4 in Fig. 19). Between the three models analyzed, the stress appears to be bounded over the entire analyzed range by the models with the perfectly aligned channels and the one with the 1/2 cell width offset. In the authors' opinion, for DBA analysis, both of these limiting cases should be considered.

5 Conclusions

This work presents several methods to derive a simplified model that can be used for bounding the stress and strain experienced by a PCHE with the popular zigzag channel geometry subjected to thermohydraulic loads. Sensitivity analysis has shown that the stress developing in a corrugated exchanger with a uniform channel cross section can be accurately represented by a straight channel model constructed using the channel profile taken perpendicular to the local flow direction of the corrugated model. This work also showed that 2D models solved with Gplε formulation can be used to simulate the stresses developing in a 3D freely supported exchanger. Examination of the results obtained from plε and plσ show that the former may overpredict the stresses, depending on the selection of the SFT, and the latter may underpredict the von-Mises stress because it neglects all stress components in the streamwise direction. Analysis comparing the cross-channel and streamwise stresses induced by realistic temperature profiles shows that the streamwise thermal stresses may be significantly greater than the cross-channel stresses, depending on the geometry of the component. For lab-scale (long and narrow) components, the cross-channel stress may be 3–20 times greater than the streamwise stress while for the utility-scale exchangers (cubical), the cross-channel stress may be 1/3 of the streamwise stress.

Analysis of the stress distribution within the cross section of a 60 kW PCHE demonstrates the complexity of the distribution of stress caused by the confinement of the porous and pressurized core by a solid envelope. Because channels in different regions within the cross section experience different stress distribution, we argue that a simplified model, such as a cell with a periodic boundary condition, cannot accurately simulate the mechanical response of the larger model. Instead of accurately representing the distribution of stress within the model, this work presents a method to bound the stress experienced by a model, and has shown that a simplified model can be used to reliably bound the stress of a large model. Such dataset may be used for the sizing of the component in accordance with ASME rules, since it will result in a conservative design from a structural perspective.

Finally, this work investigated the sensitivity of the stress distribution to the vertical alignment of the PCHE channels. We have shown that the assumption of perfectly aligned channels may be inaccurate using optical scans of specimens extracted from a real component. To investigate the effect the misalignment has on the stress distribution, three models simulating the range of possible channel arrangement were created. Results extracted near the high pressure channel from the solved models have shown an increase in peak stress of approximately 20% due to channel misalignment. Comparison between the three models also showed that two data-sets, the perfectly aligned and 1/2 offset datasets, may be used for bounding the stress on all channel alignment possibilities.

Acknowledgment

This research was performed using funding from the U.S. Department of Energy Office of Nuclear Energy's under an Integrated Research Program entitled “Advancements toward ASME Nuclear code case for compact heat exchangers (IRP-17-14227)” with Award No. DE-NE0008714 and Nuclear Energy University Programs entitled “Thermal Hydraulic and Structural Testing and Modeling of Compact Diffusion-Bonded Heat Exchangers for Supercritical CO2 Brayton Cycles (NEUP-16-10578)” with Award No. DE-NE0008589. The first author would also like to extend special thanks to Stephen R. Johnston for his help in editing this paper. Messner's contribution was sponsored by the U.S. Department of Energy, under Contract No. DEAC02-06CH11357 with Argonne National Laboratory, managed and operated by UChicago Argonne LLC.

Funding Data

  • Nuclear Energy University Programs (Grant No. 0008714; Funder ID: 10.13039/100006999).

Nomenclature

Gplε =

general plane strain formulation

plε =

plane strain formulation

plσ =

plane stress formulation

Mechanical Properties
SFT =

stress-free temperature

σvm =

von Mises stress

σy =

yield stress

Thermal Properties
Cp =

heat capacity

K =

conductivity

α =

thermal expansion coefficient

Acronyms
ASME =

American Society of Mechanical Engineers

DBA =

design by analysis

DBE =

diffusion bonded heat exchangers

DBR =

design by rule

GT =

Georgia Institute of Technology

IRP =

integrated research project

NGNP =

next generation nuclear power

PCHE =

diffusion bonded printed circuit heat exchangers

References

1.
Nestell
,
J.
, and
Sham
,
T. L.
,
2015
, “
ASME Code Considerations for the Compact Heat Exchanger
,” Oak Ridge National Laboratory, Oak Ridge, TN, Report No.
ORNL/TM-2015/401
.https://www.osti.gov/biblio/1214514-asme-code-considerations-compact-heat-exchanger
2.
Liao
,
J.
, and
Utley
,
D.
,
2019
, “
Study on Reactor Vessel Air Cooling for Westinghouse Lead Fast Reactor
,”
Nucl. Technol.
, 206(2), pp.
191
205
.10.1080/00295450.2019.1599614
3.
Liao
,
J.
,
Ferroni
,
P.
,
Wright
,
R. F.
,
Bachrach
,
U.
,
Scobel
,
J. H.
,
Sofu
,
T.
,
Tentner
,
A. M.
,
Lee
,
S. J.
,
Epstein
,
M.
,
Frignani
,
M.
, and
Tarantinoe
,
M.
,
2021
, “
Development of Phenomena Identification and Ranking Table for Westinghouse Lead Fast Reactor's Safety
,”
Prog. Nucl. Energy
, 131, p.
103577
.10.1016/j.pnucene.2020.103577
4.
Keating
,
R.
,
Nestell
,
J.
,
McKillop
,
S.
,
Allen
,
T.
, and
Anderson
,
M.
,
2019
, “
Potential ASME Code Case for Construction of Compact Heat Exchangers in High Temperature Reactors
,” Proceedings of the International Conference on Nuclear Engineering (
ICONE
),
The Japan Society of Mechanical Engineers
, May 19–24, p.
1138
.10.1299/jsmeicone.2019.27.1138
5.
Ngo
,
T. L.
,
Kato
,
Y.
,
Nikitin
,
K.
, and
Ishizuka
,
T.
,
2007
, “
Heat Transfer and Pressure Drop Correlations of Microchannel Heat Exchangers With s-Shaped and Zigzag Fins for Carbon Dioxide Cycles
,”
Exp. Therm. Fluid Sci.
,
32
(
2
), pp.
560
570
.10.1016/j.expthermflusci.2007.06.006
6.
Pidaparti
,
S. R.
,
Anderson
,
M. H.
, and
Ranjan
,
D.
,
2019
, “
Experimental Investigation of Thermal-Hydraulic Performance of Discontinuous Fin Printed Circuit Heat Exchangers for Supercritical CO2 Power Cycles
,”
Exp. Therm. Fluid Sci.
,
106
, pp.
119
129
.10.1016/j.expthermflusci.2019.04.025
7.
Jentz
,
I. W.
, and
Anderson
,
M. H.
,
2021
, “
Coupled Heat Transfer and Hydraulic Modeling of an Experimental Printed Circuit Heat Exchanger Using Finite Element Methods
,”
ASME J. Therm. Sci. Eng. Appl.
, 13(3), p. 031017
.10.1115/1.4048312
8.
Pua
,
L. M.
, and
Rumbold
,
S. O.
,
2003
, “
Industrial Microchannel Devices: Where Are We Today?
,”
ASME
Paper No. ICMM2003-1101.10.1115/ICMM2003-1101
9.
Le Pierres
,
R.
,
Southall
,
D.
, and
Osborne
,
S.
,
2011
, “
Impact of Mechanical Design Issues on Printed Circuit Heat Exchangers
,”
Proceedings of sCO2 Power Cycle Symposium
, University of Colorado Boulder, Boulder, CO, May 24–25, pp.
24
25
.https://www.heatric.com/app/uploads/2018/04/Impact-of-mechanical-design-issues-on-P CHE.pdf
10.
Carlson
,
M. D.
,
2016
, “
The Selection Evaluation and Rating of Compact Heat Exchangers (Search) Software Suite? code Capabilities and Experimental Comparison
,” Sandia National Lab. (SNL-NM), Albuquerque, NM, Report No.
SAND2016-2257C
.https://www.osti.gov/biblio/1347637-selection-evaluation-rating-compact-heat-exchangers-search-software-suite-code-capabilitiesexperimental-comparison
11.
Mahajan
,
H. P.
,
Devi
,
U.
, and
Hassan
,
T.
,
2018
, “
Finite Element Analysis of Printed Circuit Heat Exchanger Core for High Temperature Creep and Burst Responses
,”
ASME
Paper No. PVP2018-84748.10.1115/P VP2018-84748
12.
ASME
,
2015
, “
Case N-861: Satisfaction of Strain Limits for Division 5 Class a Components at Elevated Temperature Service Using Elastic-Perfectly Plastic Analysis
,”
The American Society of Mechanical Engineers
, New York.
13.
ASME
,
2015
, “
Case N-862: Calculation of Creep-Fatigue for Division 5 Class a Components at Elevated Temperature Service Using Elastic-Perfectly Plastic Analysis
,”
The American Society of Mechanical Engineers
, New York.
14.
Lillo
,
T. M.
,
Williamson
,
R. L.
,
Reed
,
T. R.
,
Davis
,
C. B.
, and
Ginosar
,
D. M.
,
2005
, “
Engineering Analysis of Intermediate Loop and Process Heat Exchanger Requirements to Include Configuration Analysis and Materials Needs
,”
Idaho National Laboratory (INL)
, Idaho Falls, ID, Report No.
INL/EXT-05-00690
.10.2172/911251
15.
Lee
,
Y.
, and
Lee
,
J. I.
,
2014
, “
Structural Assessment of Intermediate Printed Circuit Heat Exchanger for Sodium-Cooled Fast Reactor With Supercritical CO2 Cycle
,”
Ann. Nucl. Energy
,
73
, pp.
84
95
.10.1016/j.anucene.2014.06.022
16.
Pra
,
F.
,
Tochon
,
P.
,
Mauget
,
C.
,
Fokkens
,
J.
, and
Willemsen
,
S.
,
2008
, “
Promising Designs of Compact Heat Exchangers for Modular HTRs Using the Brayton Cycle
,”
Nucl. Eng. Des.
,
238
(
11
), pp.
3160
3173
.10.1016/j.nucengdes.2007.12.024
17.
Simanjuntak
,
A. P.
, and
Lee
,
J.-Y.
,
2020
, “
Mechanical Integrity Assessment of Two-Side Etched Type Printed Circuit Heat Exchanger With Additional Elliptical Channel
,”
Energies
,
13
(
18
), p.
4711
.10.3390/en13184711
18.
Chen
,
M.
,
Sun
,
X.
,
Christensen
,
R. N.
,
Skavdahl
,
I.
,
Utgikar
,
V.
, and
Sabharwall
,
P.
,
2016
, “
Pressure Drop and Heat Transfer Characteristics of a High-Temperature Printed Circuit Heat Exchanger
,”
Appl. Therm. Eng.
,
108
, pp.
1409
1417
.10.1016/j.applthermaleng.2016.07.149
19.
Kim
,
I. H.
, and
No
,
H. C.
,
2013
, “
Thermal–Hydraulic Physical Models for a Printed Circuit Heat Exchanger Covering He, He–CO2 Mixture, and Water Fluids Using Experimental Data and CFD
,”
Exp. Therm. Fluid Sci.
,
48
, pp.
213
221
.10.1016/j.expthermflusci.2013.03.003
20.
Katz
,
A.
,
Aakre
,
S. R.
,
Anderson
,
M. H.
, and
Ranjan
,
D.
,
2021
, “
Experimental Investigation of Pressure Drop and Heat Transfer in High Temperature Supercritical CO2 and Helium in a Printed-Circuit Heat Exchanger
,”
Int. J. Heat Mass Transfer
,
171
, p.
121089
.10.1016/j.ijheatmasstransfer.2021.121089
21.
Cheng
,
A. D.
,
Rencis
,
J. J.
, and
Abousleiman
,
Y.
,
1970
, “
Generalized Plane Strain Elasticity Problems
,”
WIT Trans. Modell. Simul.
,
10
, pp.
167
174
.https://www.witpress.com/Secure/elibrary/papers/BE95/BE95020FU.pdf
22.
Lemmon
,
E. W.
,
Bell
,
I.
,
Huber
,
M. L.
, and
McLinden
,
M. O.
,
2018
, “
NIST Standard Reference Database 23: Reference Fluid Thermodynamic and Transport Properties-REFPROP
,” National Institute of Standards and Technology,
Gaithersburg, MD
.