A multi-scale conjugate heat transfer modelling approach for corrugated heat exchangers International Journal of Heat and Mass Transfer

The paper compares two serrated plate-ﬁn Heat Exchanger (HE) corrugation modelling methods using Computational Fluid Dynamics (CFD). The ﬁrst method follows closely recent literature studies and models a ﬁnite length single channel of a corrugation layer inside the HE core. The second method utilises the conjugate heat transfer methodology and models a section of the HE core with both cold and hot ﬂuid streams separated by a solid conducting wall (HE corrugation). The results of latter model are then extrapolated for the full dimensions of a HE core layer to obtain ﬂow and heat transfer characteristics. The conjugate heat transfer analysis methodology presented is novel and eliminates the need for analyt- ical/empirical modelling currently widely used within industry. Furthermore, it provides more detailed information about the ﬂow and heat transfer inside the HE core enabling potential for more efﬁcient HE designs. Predictions at the corrugation level were carried out at 88 6 Re corrug 6 2957 with mesh independence studies completed for all the computational domains. The results obtained in the HE corruga- tion predictions were then implemented to the multi-scale HE unit model where the ﬂow inside the HE core was modelled using two porous media simpliﬁcations whilst the heat transfer was simpliﬁed using the effectiveness source term. The HE unit predictions were validated against industrial experimental data with good agreement found between the numerical and experimental results. All the simulations were completed using the open-source CFD package OpenFOAM.


Introduction
Heat Exchangers (HE) are the devices used widely in aerospace, automotive and other industries in which heat is transferred between two or more fluid streams, provided there is a significant temperature gradient between them [1]. Compact HE have a high heat transfer area to flow volume ratio, which is typically the case where flow hydraulic diameter (d h ) is of % 8 À 9 mm or lower [2]. Those are often used in aerospace due to the small packaging size and are typically of a platular type where the flows are separated by the plates. Kays and London [3] pioneered the design of compact HE, experimentally assessed a wide range of HE corrugation designs and developed an analytical methodology (NTU À model) enabling estimation of the HE core performance. However, such methodology provides only a limited understanding about the flow and heat transfer inside the HE.
Thus, in recent years there has been increasing number of studies using CFD aiming to improve the efficiency of flow and heat transfer inside the HE. However, because of the high surface area density to the flow volume the numerical analysis of compact HE is challenging and typically requires splitting the simulations into two parts to make the computations feasible. Firstly, a small section of HE core is analysed in detail to obtain flow and thermal performance data of the HE corrugation. This data is then used in a macro HE unit model where the flow and heat transfer inside the HE core are simplified by employing porous media and heat transfer effectiveness models to obtain pressure drop and thermal characteristics. The flow in the HE headers, however, can be modelled using conventional CFD with recent examples of similar related studies for compact HE in [4][5][6][7]. The two step approach also highlights the need for accurate detailed modelling of HE corrugation in order to predict performance of the complete unit. A typical approach to the HE corrugation analysis in the literature is to consider only a single period of a certain HE corrugation. It is typically modelled by employing a periodic flow and heat transfer assumption with either a constant temperature or constant heat transfer rate condition. The approach was introduced by Patankar [8] in 1977 and is still observed in the recent studies [9,10]. Whilst this approach is computationally inexpensive it has two disadvantages. Firstly, it does not capture complex entry/exit effects of the HE core shown in [11]. Secondly, estimating heat transfer requires special treatment as described in [8] since with the periodic flow assumption standard heat transfer conditions would result in the perpetual addition of energy. As alternatives to this method, researchers have developed finite length corrugation models such as in [12,13] and introduced conjugate heat transfer [14,15] aiming to improve solution accuracy.
In this paper a more recent computational approach using a single finite length channel similar to [12,13] is compared to a novel method utilising a conjugate heat transfer approach by modelling a slice of the HE core. The flow is simulated at 88 6 Re corrug 6 2957 using laminar and RANS k À xSST modelling assumptions [16]. Data obtained in these detailed simulations is applied to a HE unit model where two porous flow models are used to simplify the flow. A heat transfer effectiveness () model, measuring the actual heat transfer/maximum heat transfer possible at certain conditions is used to simplify the heat transfer inside the HE core. The HE unit data is then validated against industrial experimental data.

Background & methodology
A compact plate-fin HE unit shown in Fig. 4 (assembled form) with an asymmetric serrated corrugation (Figs. 1 and 3) was selected and used for both numerical and experimental testing. The corrugation has a hydraulic diameter of d h ¼ 1:134 mm with a channel height of % 2:5 mm and a mean width between the serrations of % 0:86 mm. The HE core is composed of 10 cold (JET A-1 fuel) and 9 hot (BP Turbo Oil 2380) flow layers with the identical corrugation in cross-flow with the solid (aluminium) separating the flows. Fig. 3c shows the schematic of the core with Table 1 providing dimensions of the HE core components.
Prior to the more computationally extensive work, the complex asymmetric corrugation (Fig. 3a) was simplified (Fig. 3b) modelling a single period of the geometry in order to aid the meshing process using the fully periodic flow assumption [17]. The two CFD domains used for modelling the HE corrugation in this study are shown in Fig. 2. A finite length single column model is shown in the Fig. 2a. Using a single column domain has been established in the literature, but has significant disadvantage because of the limitations in terms of only considering heat transfer via a constant temperature or heat flux boundary conditions through the walls of the domain. These boundary conditions are not directly suitable in this case because the computations cannot account for the crossflow nature of the heat transfer between two fluid streams. These effects are of particular importance in many cross-flow HE applications where one of the streams does not remain at a constant temperature throughout the HE core. It should be noted that estimation of the effectiveness inside the HE core using the single column model data is feasible, but would rely on applying corrections using the empirical/analytical data, provided by Kays and London [3], thus, not done in this study. The predicted pressure drop from the single column computational domain were useful for verifying the HE section model described below.
A proposed improved numerical approach that enables the direct capture of the cross-flow heat transfer effects occurring in many HE is to take a section of a HE core. Such a computational domain consists of both hot and cold fluid streams which are separated by a conducting solid. A schematic diagram of such section is given in Fig. 2b with the output of the solid domain with a size of the 1=16th of the layer area of the HE core in Fig. 7. The result from this approach can then be extrapolated in both flow-wise and cross-flow directions to obtain the overall HE core layer performance for a particular flow condition.
The CFD predictions were completed using OpenFOAM with both laminar and k À xSST RANS modelling assumptions with the governing equations used [18]: where u -velocity vector, p -pressure, q -density, T -temperature, S and Q -source terms to the momentum and energy equa- where l t is the dynamic turbulent viscosity and calculated based on the turbulence model [16] and a eff ¼ a þ l t =ðqPr t Þ, where a ¼ k=qC p and Pr t is the turbulent Prandtl number taken to be 0.85 in RANS simulations. The flow inside the HE core ( Fig. 3c) was modelled using the two porous media approaches for the HE core assuming each corrugation layer is a single porous region and employing the Darcy-Forchheimer model where d and f are the Darcy and Forchheimer constants respectfully. The Darcy coefficient d in these simulations was left as 0 eliminating the term as the flows are of low viscosity and forced convection dominant and allowed calculating the f coefficient for a selection of flow rates from the detailed simulation data. An alternative power law model was also implemented which is defined as where C 0 and C 1 are the curve fitted constants which were obtained from the detailed simulation data. The implementation of the two models is quite different within OpenFOAM. The Darcy-Forschheimer model allows manually specifying constant flow resistance values in three directions for a specific case. The power law model, however, considers the mass flow reaching each porous layer and then selects a uniform resistance level in all three directions, which could be beneficial in the cases of HE cores where maldistribution is likely. The heat transfer inside the HE core was modelled using the effectiveness model which for a single fluid is implemented as where _ m 1 ; _ m 2 -primary and secondary mass flows, c p -specific heat, T 1 ; T 2 -primary and secondary fluid temperatures, ð _ m; _ m 2 Þ -effectiveness lookup table. It should be noted that although the HE core is simplified, the flow inside HE headers is resolved using fully detailed CFD.
The HE unit flow volumes (headers and core included) together with the solid region of the HE core that separates the two flows are shown in Fig. 5 were extracted from the HE unit shown in the Fig. 4a. CFD modelling of only the cold flow side of the HE unit was undertaken (blue flow volume in Fig. 5) using the same meshing methodology with the RANS k À xSST [16] model. This was done to reduce the overall computational cost and was appropriate as the cross-flow heat transfer effects were taken into account by the HE section model. Constant fluid and solid properties throughout the simulations were taken based on the inlet temperatures of both fluids during the experiments and are given in Tables 2 and 4. This is acceptable due to the small experimental temperature   Table 3. It should be noted, that prior to using these schemes the solutions for each case were initialised using first order accurate schemes where possible (bounded Gauss upwind in OpenFOAM) to provide a close initial solution (further detail is outlined in [11]).

Experimental setup
The experiments were carried out using an assembled prototype HE unit shown in Fig. 4 at an industrial facility and were split into two parts. Firstly the pressure drop was measured through the HE cold flow side under isothermal conditions. Then the thermal performance of the HE was established by varying the cold flow HE stream rates whilst maintaining a constant low flow rate through the hot side of the HE unit. A Reynolds numbers range of 88 6 Re corrug 6 2957 (assuming uniform flow expansion in the inlet HE header) was covered allowing for an opportunity to observe the unsteady flow occurrence in CFD simulations previously resolved numerically by the author [11] and shown experimentally by Rush et al. [19]. Repeatability of the thermal performance experiments was ensured by running the same test points twice and measuring the heat balance in both cases. It was calculated using the relationship: The heat balance in Fig. 6 shows that the results were repeatable throughout the Reynolds number range tested and was also found above 100% (on average % 105%) which can be explained by some inevitable heat loss to the environment despite insulating the HE unit. It should be noted that the thermocouples used in the experiments had a measurement tolerance of AE0:3 K and if it was assumed that for the both fluids the highest measurement deviation occurred (DT hot þ 0:6 K & DT cold À 0:6 K) it could lead to an % 10% error in Q balance . However, an error this size is unlikely as the experiments were verified by repeating the experimental     points twice (Fig. 6). Figs. 4b and c also indicate that the unit had some imperfections from manufacturing processes employed such as welding and vacuum brazing. Whilst it is difficult to quantify their exact effect, the reduction in cross-section could be assumed to increase the pressure drop across the unit by up to 5-10% and would reduce heat transfer.

Mesh independence of the computational domains
A number of meshes were generated using inbuilt OpenFOAM blockmesh and snappyHexMesh tools for the mesh independence study as summarised in the Tables 5 and 8. Overall, meshing a single channel domain was found less challenging computationally and allowed higher resolution of the mesh boundary layer compared to the HE section model due to the smaller size of the domain. The predictions of the detailed domains were completed at an experimental point of the thermal performance test with the mean Re corrug ¼ 146 for the cold and Re corrug ¼ 35 for the hot side corrugations respectively (assuming even expansion at the HE headers). Inlet temperatures of T inÀcold ¼ 302 K and T inÀhot ¼ 373 K were used for the cold and hot streams respectively for the HE section model calculations and they corresponded to the set experimental inlet temperatures at the HE header inlets with the example outputs of the model given in Fig. 7. As the nature of the single channel modelling approach limits the heat transfer methods the cold flow side with a constant temperature boundary with a mean temperature of T ¼ 0:5ðT inÀcold þ T inÀhot Þ ¼ 338 K used. The summarised mesh independence data for both computational domains can be found in Tables 6 and 7. Both tables compare the overall performance obtained by taking plane averages at the inlets and outlets of the two domains. In both cases the medium meshes were considered sufficiently accurate. However, it should be noted that whilst the pressure drop values in both cases stabilised, the overall temperature data was found to change by a small extent between the medium and fine meshes (% 6%). The medium accuracy in both domains was selected to balance the computational resources and power requirements. The simulations for this article were implemented at the local high performance facility using 48 cores of Intel Broadwell E5-2560 v4. Typically an individual simulation of a HE section model at the medium mesh resolution (as shown in Fig. 7)) took two days to compute under steady-state assumption.
Transient solutions of the corrugation were run and showed the unsteady flow behaviour for both computational domains at the medium mesh resolution and Re corrug ¼ 1213. The steady-state assumption was judged acceptable in this case as the transient    fluctuations were small when compared with [11]. It should be noted though that a previous study [11] identified that the steady-state assumption would result in some under-prediction of the pressure drop. For this particular application it is estimated to be up to %> 5% inside the HE core section. The summarised mesh independence data for the cold side of HE unit model is shown in Tables 8 and 9. Darcy-Forchheimer and effectiveness models were used to simplify the flow and heat transfer inside the HE core for the mesh independence calculations. Uniform friction factor of f ¼ 714 was used which was extracted from one of the HE section mesh independence results together with an effectiveness ¼ 0:61 prescribed across the cold HE core side. Meshes were run using the thermal performance test fluid properties at Re unit ¼ 11120 prescribed at the cold flow inlet (Re corrug ¼ 146 at the corrugation level assuming even expansion of the flow). The overall values between the inlet and the outlet of the domain (Table 9) produced a similar trend to the results of the detailed domains. It was decided to use the Fine2 resolution mesh for the HE unit model predictions -the overall DP and DT showed a trend towards stabilising which supported that the finest mesh is of the sufficient accuracy.

Extending the detailed prediction domains
After the mesh independence simulations were obtained for both single column and HE section corrugation domains, both models were increased twice in the length in the main flow-wise directions. The same mesh density was maintained by keeping identical mesh refinement levels through the meshing settings in OpenFOAM. It resulted in modelling the full cold layer length single column domain (Fig. 2a) whilst the HE section model (Figs. 2b, 7) was doubled in length along the cold stream direction resulting in the 1=8 th layer area model. This step was undertaken to minimise the flow exit effects and also to obtain better quality data across the Reynolds number range. The selective expansion of the HE section model was completed to reduce the computational expense as only the cold stream flow rate was varied during the experiments whilst the hot stream was kept at a constant Re corrug ¼ 35. To compare the different corrugation domains the friction factor f was used derived from the Eq. (4) as where u f was taken as a velocity set at the inlet. The friction factor f data was calculated at every period of the cold side corrugation (Fig. 8) with the extended domain simulations were run at the same flow rate and fluid properties. Full length single column model data was treated as a benchmark (since it is a full length model of the cold stream layer) and showed that the difference between the two single column domains was minimal. The plot also highlights that the shorter of the two HE section model domains was not long enough and was significantly influenced by the flow exit effects. The extended HE section domain, however, agreed relatively well with the full length channel data, although, was still influenced by the flow exit effects (drop in flow resistance) at the last corrugation period. The slightly lower overall friction factors of the HE section model were predicted due to an overall lower resolution of the mesh due to a very large modelling domain. Thus, to compare the HE corrugation predictions further, both double length simulation domains were used to generate the flow resistance data across the experimental Reynolds number range. Results of the last period of the HE section model, however, are purposely omitted from the input data to the HE unit models during validation due to the unwanted exit effects discussed above.

Corrugation predictions
Results from the HE corrugation domains were used to generate the data for the pressure drop test validation, whilst the HE section model was employed to generate the data for the thermal performance test modelling. The overall flow resistance characteristics between the single channel and HE section models are compared in Fig. 9. The single column domain simulations were completed with the laminar flow assumption, whilst the HE section model was implemented using both flow momentum assumptions. A change in solution residual behaviour was observed for the single   column domain predictions at Re corrug % 300 (for the flows below the Re corrug ¼ 300 the residuals converged to 10 À12 whilst the residuals for the predictions above remained steady and at the 10 À3 level). This is believed to coincide with the onset of the unsteady flow which was not resolved directly here due to the computational cost as identified earlier. Instead the flow instability levels of the steady-state single channel predictions were considered by calculating standard deviations of overall pressure drop of iterations 2000 to 8000, sampling the results every 100 iterations. They are shown as error bars in Fig. 9, revealing the low flow unsteadiness levels to Re corrug 6 2010. This onset of the transitional flow agrees well with results from previous work [11] where for a 2D sinusoidal geometry the onset of the transitional flow was observed at Re corrug % 200. Interestingly, the same switch in residuals was not observed for the conjugate heat transfer domain (HE section model) and is thought to be associated with a coarser boundary layer of the domain not being sufficient to capture the phenomena. It highlights the necessity of a sufficient boundary layer resolution in order to predict transient flow at the higher Reynolds numbers. It may also suggest why the transitional Reynolds number regime has not been studied extensively in the literature relating to the compact HE. Despite the overall lower mesh resolution of the HE section domain the results of the two solution domains (Fig. 9) agree well in terms of the flow resistance predictions (around 5% on average). Furthermore, the lower resolution of the HE section domain resulted in a better convergence of the simulations at the higher Reynolds numbers whereas the single column model starts to diverge at Re corrug ¼ 2484 and onwards as shown in Fig. 9. Mainly the differences between the two corrugation domains were found around the flow instability at Re corrug % 600 À 800 where the HE section domain showed a minimal fluctuation in terms of the flow resistance compared to the single column domain. This common fluctuation between the different models also suggests a potentially highly transitional flow in this HE unit around these Reynolds numbers.
The detailed flow and heat transfer predictions for the thermal performance test are presented in Fig. 10. Both were obtained using the HE section model and with the two flow momentum modelling assumptions (laminar and k À xSST) and presented for the f factor and the overall effectiveness data of the cold side of the HE section model which is calculated as [3]: where the C cold ¼ _ m cold C p;cold ; C min ¼ MINðC cold ; C hot Þ and C hot ¼ _ m hot C p;hot . Interestingly, both flow momentum assumptions predicted nearly identical results across the Reynolds number range which is contrary to the previous result by the authors [11] for sinusoidal corrugation where the pressure drop was found to be overpredicted by up to 30% (increasing with the Reynolds number) by the laminar momentum model. There are a number of potential reasons for this: firstly, the serrated corrugation is much less disruptive and leads to less mixing compared to the sinusoidal geometry in [11] which enables validity of the laminar flow assumption at the higher Reynolds numbers. Another reason again is the coarser mesh resolution in 3D (compared to [11] where a 2D geometry was used) which could have caused both the suppression of flow instability and the lack of contrast between the laminar and k À xSST models.
In addition, the start of the difference between the two models in [11] was also observed at the onset of the flow transition which in this study was not resolved due to the excessive computational cost.
The lowest effectiveness values presented in Fig. 10b were also produced around the Re corrug % 300 as suggesting the suspected onset of the flow instability and increasing the role of flow mixing in vicinity of this point. It also shows that with the increasing flow rate the effectiveness of the HE unit increases which supports the importance of mixing in HE. The thermal HE section data was extrapolated for the complete layer prior to being used as the input for the HE unit modelling. This was completed in two directions: parallel and normal to the cold side flow direction. Firstly, it was found that the vast majority of heat transfer in this application was predicted to occur during the flow development phase  (highest pressure drop section in Fig. 8). The extrapolation of the data results also indicates that the flow does not undertake a significant temperature change further down the layer. This is thought to have been due to the highest temperature gradient occurring between the HE streams, which then reduces significantly further downstream. The absence of transient flow mixing in this application, which in [11] was shown to enhance the heat transfer down the channel is identified as another potential reason for the relatively constant downstream temperature. In contrast, the thermally negative cross-flow effects on the cold side HE performance were found to be much more significant. They were measured by sampling the multiple lines across the cold flow HE section model domain, normal to the main flow-wise direction at the end of the first corrugation period (after the highest heat transfer section). The sampled data was then curve fitted across the HE core layer width to produce an estimate for a decrease in DT. The cross-flow heat transfer effects measured were proved significant and increasingly dominant with Reynolds number (% 18 À 30% in terms of DT) and required adjusting the effectiveness values for the HE unit model input.

Comparing HE unit predictions to the experimental data
The two flow resistance models were implemented together with the effectiveness model using the detailed CFD data into the HE unit model (Fig. 5). Outputs of both porous media models were very similar in terms of the overall results despite the different implementation of them (Figs. 12 and 13). Only small visual differences can be observed between the two models as shown in the Fig. 11 where the temperature contours through the cold flow of the HE unit are given. The pressure drop test results in Fig. 12 generally agree well, but with increasing under prediction at the higher Reynolds numbers reaching an % 10% difference at the highest of the Reynolds numbers tested. Some of the discrepancy could be explained by the experimental error which can be assumed to be up to at least 5% (increasing with the flow rate, as estimated by the error bars in Fig. 12) which arises both from experimental equipment and the afore mentioned manufacturing defects. Additionally, the diverging agreement between the numerical predictions and the pressure drop experiments could also be partly contributed to steady-state assumption used throughout the detailed corrugation simulations. A previous study by the authors [11] identified that the pressure drop can be underpredicted in the transitional Reynolds number regime if the steady state solution is maintained. In this case, the effect of a steady state assumption was lower as a less disruptive corrugation was used (sinusoidal versus serrated) and the underprediction could be estimated to be % 5 À 10%.
Overall, a good agreement was found in terms of the thermal quantities presented in Fig. 13. On average, the difference between the DT measured was 0:94 K whilst the average mismatch in terms Fig. 11. Temperature (K) contours with 302-308 K temperature range through the cold HE side inlet/outlet ports at Re = 66680 (HE inlet level) Re = 880 (corrugation level) of (a) Darcy-Forschheimer and (b) Power law models.  of effectiveness was at % 4:8%. The small differences in predictions can be partially explained by the two experimental factors. Firstly, the thermocouples used in experiments had a tolerance of AE0:3 K. Fig. 13b shows the error bars of the experimentally measured effectiveness which increase with the Reynolds number, arising from the smaller DT cold measured, stressing the necessity of both accurate thermocouples and numerics in order to calculate the derived heat transfer quantities accurately. Secondly, during the thermal performance experiments the heat balance measured throughout ( _ Q hot = _ Q cold ) was at % 105% explaining the overprediction of heat transfer at the higher end of the Reynolds range. In terms of the numerical results, in the first third of the Reynolds number range tested the underprediction of heat transfer in terms of both DT and is observed after the first two analysis points and is thought to be associated with the onset of the transitional flow. However, from the Re inlet % 40000 onwards the numerics agree very well in both quantities (especially when the error bars of the experiments are considered). Furthermore, it is interesting that a very similar initial drop in the heat transfer effectiveness was observed experimentally as well as numerically (although slightly delayed in terms of Re), suggesting the onset of significant mixing onwards from low Reynolds numbers. Fig. 11 indicates levels of maldistribution occurring within the HE core. This was further analysed by taking plane averages of flow and heat transfer quantities entering and exiting each layer of the HE core. This output is provided in Fig. 14 where the HE core layers were numbered from 1 to 10 (bottom to the top). A clear preference for the flow entering the channels 2 and 3 is seen (15% and 12% of all the mass flow is entering the layers respectively) as could have been expected due to the shape of the HE headers. It should be noted that the layer 1 was half blocked due to the welding of the headers to the core limiting the mass flow through it as shown in Figs. 4 and 11. A similar behaviour was also observed throughout the simulations for the Reynolds number range tested. Both core flow simplifications (Power law -blue and Forschheimerred lines) also produced a close agreement in terms of the maldistribution of the flow between the HE layers. However, the slightly reduced maldistributed prediction of the Forschheimer model led to higher temperature change prediction across the individual HE core layers. In this case it did not make a significant effect for the overall HE performance, but it could become more important in other applications or designs.

Conclusions
A novel HE corrugation modelling approach has been proposed and evaluated at the Reynolds number range of 25 6 Re corrug 6 2957. It takes a section of an HE core and models both cold and hot HE streams separated by a conducting solid, implemented using a conjugate heat transfer methodology. The method eliminates the need for analytical/empirical relations such as in [3] and provides with potential for a more efficient HE design. The novel HE section domain data was implemented in the macro HE unit models where the HE core was simplified using porous media and effectiveness models for the flow and heat transfer. Main findings of the study are: HE section domain was verified using the flow resistance data of a single channel corrugation model. Strong agreement was observed across a wide Reynolds number range between the predictions with the HE section domain producing marginally lower flow resistance. The single channel domain also detected a suspected onset of transitional flow at Re corrug % 300, levels of which were found significantly smaller than in [11,19]. The two HE unit models with Darcy-Forschheimer and the powerLaw porous media HE core flow simplifications predicted only a marginally different maldistribution level inside the HE core. This leads to only small overall differences being observed, despite the different implementation of them in OpenFOAM for both heat transfer and pressure drop. The HE unit model was compared to the experiments. Overall, a strong agreement was found between the experimental and numerical HE unit results when comparing the overall pressure drop (maximum % 10% disagreement at the highest flow rate) and the temperature change (on average DT % 0:94 K and % 4:8% disagreement) of the HE unit. This agreement validates the novel HE corrugation evaluation methodology and makes this approach an important achievement for the design of such HE. It potentially enables the design of an improved performance HE compared to the existing approaches, being achieved through obtaining detailed information on the crossflow heat transfer effects.