Machine learning for rapid discovery of laminar flow channel wall modifications that enhance heat transfer

Heat transfer in fluid flow is an important physical phenomenon, with relevance across many science and engineering disciplines ranging from microfluidic devices in chemical engineering and biomedical implants, all the way to high-temperature physics and cosmology. In this proof-of-concept study, we explore an interesting engineering question, which can be posed as follows: “Is it possible to introduce structural changes to the wall of a channel that increases heat transfer, without having a corresponding increase in the pressure drop?” This fundamental question linked to the ultimate goal of dissimilar flow control or dissimilar heat transfer enhancement that has been asked for decades by various research groups in different application fields.

Dissimilar heat transfer enhancement is proven to be extremely challenging due to similarity in the mechanisms of momentum and heat transfer.¹ Investigations of various surfaces including specially designed fins,^2,3 dimples,⁴ or vortex generators⁵ report that an increase in heat transfer (where heat transfer is described by the Stanton number St) is always accompanied by inevitable manifold increase in pressure drop (described by the drag coefficient C_f). This eventually results in a decrease in the Reynolds analogy factor RA = 2St/C_f⁶ in comparison to a flat channel configuration. It is, however, known that a dissimilar modification of momentum and heat transfer is possible when sophisticated flow control methods are applied. Those control methods, for instance, can be based on the introduction of flow perturbations or optimally distributed blowing/suction profiles from the wall surface.^7–11 These studies confirm that a significant enhancement of the Reynolds analogy factor (tripling RA in comparison to the uncontrolled channel flow) is possible when an appropriate flow manipulation is created. It is found that an introduction of large-scale spanwise rolls significantly promotes heat transfer while the drag coefficient remains less affected. This concept has also been successfully tested in the framework of turbulent channel flows, where RA > 2 can be achieved instead of RA = 1 in an uncontrolled flow configuration.^12,13 Recent studies in turbulent flows also report a possibility of RA modification using streamwise elongated structures, leading to the formation of turbulence-driven secondary motions.¹⁴ The modification of RA is, however, limited in this case to several percent due to increase in the wetted area and the corresponding increase in C_f.

Most of the above studies were focused on turbulent flow regimes. To simplify the scenario, in the present study, we consider the laminar flow regime in a two-dimensional channel with heat transfer. The immersed boundary method is used for introducing the wall structuring. This approach allows for rapid execution of direct numerical simulations (DNSs), making it possible to investigate a wide range of arbitrarily generated surfaces.

Due to rapid growth in ML method development, there is also a growing interest in applying ML methods to challenges in the field of fluid dynamics. Seminal work on ML in fluid mechanics was focused on image processing based on convolutional neural networks (CNNs). CNNs have been successfully utilized as a straightforward method for learning from flow field data, reducing the computational cost.^15,16 For instance, deep CNNs have been used to construct velocity fields from image pairs of particle image velocimetry.¹⁷ Furthermore, the CNN-based superresolution algorithm has contributed to turbulent flow reconstruction.¹⁸ A recent study established two different models based on CNNs and encoder–decoder (ED) architectures to predict characteristics of the flow and heat transfer around the NACA.¹⁹ The established CNN predicts the aerodynamic coefficients and the Nusselt number. However, there are still open challenges approaching dissimilar heat transfer enhancement and applying ML algorithms to this challenge. Furthermore, interpreting ML models can offer additional insights to complement human understanding and engineering intuition.

In this study, we present a workflow consisting of numerical simulations (Sec. II A), the key parameters for quantifying the performance of a thermo-fluid system (Sec. II B), and dataset generation along with the ML model (i.e., the CNN model) training (Sec. II C). In Sec. III, We show that the ML model can predict fluid flow and heat transfer characteristics with a large speedup compared to numerical simulations and with a high enough accuracy to prescreen a large database of possible channel geometries. Following an analysis of the velocity and temperature fields of some of the relevant structures in the training dataset (Sec. IV A), we interpret the trained ML model using SHapely Additive exPlanations (SHAP), which is a widely used explainable artificial intelligence method (Sec. IV B). Finally, we summarize the conclusions of the current study in Sec. V.

Due to the CFR approach, the bulk Reynolds number is fixed to Re_b = 2U_bδ/ν = 200 for all the considered simulations, where $Ub=12δ∫02δudy$ is the bulk mean velocity. This means that any modification of the flow is translated into an alteration of the resulting mean streamwise pressure gradient P_x required to maintain the chosen flow rate. Periodic boundary conditions are applied in the streamwise directions while the wall-normal extension of the flow domain is bounded by no-slip boundary conditions at the lower and upper domain walls (y = 0, 2δ). Subscript l and u are used throughout the manuscript to denote quantities on the lower and upper walls, respectively.

The solver implementation is based on a spectral solver for incompressible boundary layer flows.²⁴ The Navier–Stokes equations are numerically integrated using the velocity-vorticity formulation by a spectral method with Fourier decomposition in the horizontal directions and Chebyshev discretization in the wall-normal direction. For temporal advancement, the convection and viscous terms are discretized using the third-order Runge–Kutta and Crank–Nicolson methods, respectively. The flow domain is discretized with N_x × N_y = 256 × 129 grid nodes, while the immersed boundary method is applied on the dealiased grid (3/2 rule) with M_x × M_y = 384 × 129 grid nodes.

The obtained supporting points are then interpolated with cubic Bézier curves. The distance of the control points from the supporting points is determined by the parameter r. For small values of r, the radius of the curves can get very small, resulting in sharp features. The parameter a controls the smoothness of the curve. The values of both r and a are varied in the range [0, 0.2] as mentioned in Table I. For a = 0, the angle through one supporting point is determined by the mean of directions to both neighboring points. At higher a, the direction to one neighboring point is weighted higher and hence the curve features stronger edge.²⁶ 

The obtained bottom and top wall curves were then projected onto the dealiased grid. Using 20 random initializations per parameter variation, 108 000 random wall structures were generated, hereafter called the repository set.

We calculated C_f and St using the simulation method described in Sec. II A for a randomly sampled subset of 10 800 channel geometries (10% of the repository set). The 10 800 calculations each requiring 20–30 min single core execution on an Intel Xeon Gold 6230 CPU result in a total cost of 3600–5400 CPU h. Out of these, 9185 met specific criteria related to temperature convergence and geometric validity. We set aside 5% (459) of these data to be included in the test set.

We used varying fractions of the remaining 8726 channel geometries (hereafter called the labeled set) for hyperparameter optimization of CNNs. Hyperparameters were determined without flat channel geometries. We used a total of 8776 channel geometries as the training set and 464 channel geometries as the test set. We added 50 flat channels with various channel heights to the training dataset and 5 to the test dataset. The inputs for the CNN are binary images where 1 and 0, respectively, represent solid and fluid domains. Each binary image with 384 × 129 × 1 pixels represents a cross section of the channel geometry, i.e., exactly the same input that is also used in the numerical simulations. The input is passed through a varying number of convolution steps, each consisting of a convolution with padding, a varying kernel size and a varying number of filters, followed by a ReLU-activation and a 2x2x max.-pooling with stride 2. The output of these convolutions is then flattened and passed through one hidden dense layer with ReLU-activation and a varying number of neurons. This hidden layer is additionally regularized by a varying dropout. C_f and St are then predicted with an output layer with two neurons and linear activation. An overview of the optimized CNN model is shown in Fig. 2(a).

The model was implemented using TensorFlow²⁷ and Keras,²⁸ and training was done using the Adam optimizer.²⁹ The model was trained for 100 epochs with a batch size of 256 using the mean squared error (MSE) loss while logarithmically reducing the learning rate from 1 × 10⁻³ at the tenth epoch to a varying final learning rate at the last epoch.

The hyperparameter optimization was done with Bayesian optimization to find the most optimal set of hyperparameters, where optimal was defined as the lowest validation set error (see Table II). The reported errors of the model were then calculated on a separate test set. It ensures that overfitting to the training set is minimized and potential overfitting to the validation set during hyperparameter optimization is not included in the reported test set errors. For this, SigOpt was employed^30,31 with 140 experiments on eight asynchronous parallel channels. Each channel has access to a single Tesla A100 GPU. To ensure stable predictions, for each evaluation the mean MSEs for the predicted C_f and St of a fivefold cross-validation were used as metrics.

Under the present flow and temperature conditions, where no volume forces are considered and temperature is treated as a passive scalar, C_f and St of a channel do not change when the channel geometry is flipped along the channel centerline or shifted along the streamwise direction. We used this information to augment the data during the ML model training. We divided input images into batches. For instance, 35 batches containing 256 images were created. Each half of the 256 images in each batch was flipped and shifted randomly during the training. The channel geometry was flipped about the y-axis and also randomly shifted in the x-direction within the range of $Mxs$ and $Mxe$ [see Fig. 2(b)]. Input images are shuffled every epoch so that the same images are not always flipped and shifted in the same order. Hence, multiple shifted and/or flipped copies of channel geometries from a single geometry were generated without the need for any additional numerical simulation. The effect of shifting and flipping on the test data is discussed in detail in the  Appendix. In the end, the computational cost of training the model is negligible (around 20 min on a NVIDIA Tesla V100 GPU).

We trained the CNN model and obtained MAEs of 1.90 × 10⁻³ (3.16% of C_f,ref) and 1.29 × 10⁻⁴ (0.81% of St_ref) for C_f and St predictions on the test dataset. This corresponds to an expected MAE of 1.25 × 10⁻² in RA (2.35% of RA_ref). The respective r²-scores for C_f and St are 0.951 and 0.925. Figure 3 shows training and testing loss as a function of epochs. The curves converge from around epoch 20 toward 100. A comparison of CNN predictions with simulated ground truth values for the test set is shown in Fig. 4. Most values of C_f and St are well correlated. For a perfect model, the predictions and the ground truth values are the same, so all points would lie on the unity line. Deviations for large values can be seen (C_f = 0.12∼, St = 0.020∼), but at the same time, this only affects very few channels that are not of practical relevance.

To demonstrate the expected exponential behavior of the learning curve and to estimate necessary training set sizes for given target errors, we generated learning curves (see Fig. 5), where we observe the MAE in C_f and St as a function of the training set size. The hyperparameters were kept constant, and the amount of training data was varied from 5% to 90% of a total of training and test set, 9240 channel geometries. With each diverse dataset, every model is trained independently, and fluctuations are to be expected due to randomly sampled initial weights. We observe an exponential decrease in MAE with increase in the training set size. It can be seen from Fig. 5 that larger datasets increase the model accuracy. In the case of C_f, the MAE score for the test dataset is approximately twice smaller than that of 5%. In the case of St, the factor is slightly smaller. As the plots continue to decrease exponentially, further improvement in the model accuracy can only be expected by substantially increasing the number of data. However, we do not observe strong saturation, indicating that the model is not limited by capacity or incomplete input information.

The scatter plot of 8726 labeled data points from the repository set used for training the CNN is shown in Fig. 6(a). We would like to point out that the training set contains geometries that are well representative of the flow configurations usually encountered in fluid dynamics. Three of such structures that are highlighted in Fig. 6(a) are discussed in detail in Sec. IV. Almost all the points fall below the RA/RA_ref = 1 line in Fig. 6(a), indicating that deviations from a flat channel lead to an increase in C_f that outweighs the simultaneous increase in St. We then exploited the speedup of the surrogate machine learning model (⁠$<100ms$ per channel) compared to the numerical simulation (⁠$≈20−30min$ per channel) to explore the flow and heat transfer characteristics (C_f and St) of a much larger set of unlabeled channel geometries. The term unlabeled set is part of the repository set for which the ground truth information is absent, signifying that no simulations have been conducted for these geometries. The calculation of C_f and St for the unlabeled channel geometries took only 2.5 GPU hours. Thus, the benefit of training a machine learning model becomes evident when the model is used for inference on new structured channel geometries. The scatter plot of C_f vs St for the repository set, split up into labeled and unlabeled data points, is shown in Fig. 6(b). In general, it is very difficult to find geometries that can result in dissimilar heat transfer enhancement, i.e., RA/RA_ref > 1. However, as can be seen from Figs. 6(a) and 6(b), there are some points very close to (and above) the RA/RA_ref = 1 line. However, these points fall within the combined uncertainty level from the simulation results plus the ML prediction. The structures corresponding to these points closely resemble either a flat channel, the trivial solution, or the nontrivial geometry of structure C (see Sec. IV A). To come up with more innovative structures with RA/RA_ref > 1, we can combine the speedup of the surrogate model with evolutionary optimization techniques.³² It is essential to note that evolutionary optimization algorithms in general require the optimization of their own set of hyperparameters. In this context, having a surrogate ML model that can provide values for the quantities of interest (RA for example) in a fraction of a second is highly beneficial. Figures 6(a) and 6(b) also show that the training data adequately represent the repository set. The histogram representation of the repository set shown in Figs. 6(c) and 6(d) confirms the same.

In this section, we analyze in detail the pressure loss and heat transfer characteristics of the three labeled structures highlighted in Fig. 6(a) (namely, structures A, B, and C). The mean flow and temperature fields of these structures are shown in Fig. 7, and their C_f, St, and RA values are listed in Table III. Interestingly, these three structures resemble three canonical flow configurations encountered in fluid dynamics: vortex generator, converging–diverging nozzle, and backward-facing step. For the three considered structures, the percentage increase in surface area with respect to a flat channel is 25.1%, 5.4%, and 6.3%, respectively.

Among the three structures, the heat transfer achieved is maximum for structure A with the vortex-generating wing-like protrusion. Vortex generators of varying shapes, similar to the wing-like protrusion of structure A, are commonly used in heat exchanger devices to introduce unsteady swirling motions that can increase heat transfer.³³ Though there are two vortices inside the flow field as can be seen in Fig. 7(a), these vortices are part of the recirculation regions. Such regions in fact isolate the wall from the bulk of the fluid and are detrimental to effective heat transfer. Nonetheless, the developing thermal boundary layers generated on both the top and bottom walls [refer to the temperature field in Fig. 7(b)] result in an 80% increase in heat transfer. At the same time, the shape of structure A together with its increased surface area of 25% introduces significant pressure loss leading to a 60% reduction in RA when compared to the flat channel.

For structure B, the flow inside gets accelerated inside the converging section followed by fluid deceleration inside the diverging section. Unlike the case of structure A, a strong variation in the thermal boundary layer thickness is present only along the bottom wall. The presence of a recirculation region, covering the entire bottom wall except for the converging–diverging section, further limits heat transfer from the bottom wall. As a result, the increase in heat transfer is only 16%. The considerable pressure loss due to the absence of a well-streamlined converging–diverging section leads to a 60% reduction in RA with respect to the flat channel.

Structure C is unique owing to the fact that it contains the backward-facing step together with wall meandering. The structure results in a RA/RA_ref value approximately equal to unity but with a 30% increase in heat transfer with respect to the flat channel. This means that the structure results in a proportional increase in pressure loss as well. With the exception of the immediate area downstream of the backward step, the streamlines in the wall-normal direction exhibit behavior akin to those found within a rectilinear channel, maintaining a uniform and evenly spaced distribution. The absence of a developing thermal boundary layer indicates that wall meandering with a 6% increase in the surface area should be the primary reason for increased heat transfer. This suggests the possibility of achieving a value greater than unity for the RA/RA_ref value with carefully selected parameters for wall meandering.^34–36 Some of the structures in Figs. 6(a) and 6(b) also indicate this trend. Nevertheless, experimental validation is required to assess the practical efficacy of these geometries, and this constitutes one of our future research objectives.

Given a specific structured channel geometry, it is worth knowing which local aspects of this given geometry influence the CNN prediction. In this regard, we used Shapely Additive exPlanations (SHAP) developed by Lundberg and Lee.³⁷ The resulting Shapley value maps for structure A are shown in Fig. 8. Since our CNN model predicts two outputs, C_f and St, we have two plots that provide local approximations for these two outputs. The outline of structure A is superimposed on top of these plots to obtain a better understanding of the local distribution of the Shapley values.

The Shapley values generated using the DeepExplainer module within the SHAP Python package are computed using the CNN model alongside a distribution of background samples. Here, we are interested in understanding how C_f and St change with respect to that of a flat channel configuration. Hence, we used around 90 flat channels that are located symmetrically or asymmetrically about y = δ as background samples. It is important to note that the summation of the Shapley value is equal to the difference between the model prediction for the given structure and the mean prediction derived from the background samples.

The Shapley value plots depicted in Fig. 8 suggest that the CNN model predominantly emphasizes the detection of geometric characteristics, such as well-defined edges, corners, and recessed regions when making predictions—especially those geometric characteristics that are responsible for the appearance of flow features, such as boundary layer growth, flow separation, and recirculation regions, that influence the local drag and heat transfer. For structure A, the wing-like protrusion on the bottom wall and the sharp corner on the top wall are the predominant features responsible for the developing boundary layer and also for flow separation leading to the generation of recirculation regions. Since these features contribute to an increase in both C_f and St, the Shapely value plots for them also look similar. Thus, even though the CNN model has been trained on binary representations of structured channel geometries, the ML model has the ability to capture geometrical features of wall structuring related to the flow physics as indicated by the Shapely value distributions. The Shapley value distribution not only helps to understand what the ML model is focusing on, but it can also be used to optimize the shape of the structured channel to achieve a desirable outcome, for example, to increase heat transfer, which we will try to exploit in future work.

We have trained a convolutional neural network (CNN) based on fluid flow and heat transfer simulations in two-dimensional laminar structured channel geometries. The trained CNN model accurately predicts the target properties, i.e., drag coefficient C_f and Stanton number St, for a structure geometry based on its binary representation. As a surrogate model, the CNN can predict the target properties at a fraction of the time (⁠$<100ms$ per channel) required by numerical simulations (⁠$≈20−30min$ per channel). The significant speedup provided by the surrogate model was then used for predicting C_f and St values of around 80 000 randomly generated geometries. This was done to find out geometries with normalized Reynolds analogy factor RA greater than unity, where RA = 2St/C_f is an indicator of dissimilar heat transfer enhancement. However, our study suggests that to achieve such a goal, we have to exploit the speedup of the surrogate model by combining it with powerful optimization techniques.³² 

One of the highlights of the present CNN model is how we used domain-specific knowledge to augment the original training dataset. Recognizing that C_f and St remain invariant when a structure geometry is shifted along the flow direction or mirrored along the channel centerline, we populated the training dataset accordingly.

We showed that the training data generated using Bézier splines adequately represent structures that are commonly used in thermo-fluid applications. We then used Shapely Additive exPlanations as an explainable artificial intelligence tool to understand the CNN model. When applied to a structure, the Shapley values provide insights into the contributions of specific geometry features to the model’s output in comparison to the background reference. In this regard, the Shapley values are not only valuable in understanding the ML models but can also be used to improve their performance.

The current CNN model is trained for a specific set of flow conditions, in terms of the Reynolds number and boundary conditions. As part of future work, we are interested in extending the current ML model to cover a wide range of these conditions. Additionally, in order to further exploit ML models in general and CNNs in particular for the design of chemical engineering unit operations, we plan to implement active learning approaches and generative models to reliably explore the possible design space of structured channel geometries and directly solve the inverse problem, i.e., the suggestion of channel structures given desired target properties. Moreover, trying to predict not just the values of C_f and St but also the velocity and temperature fields inside these arbitrary geometries, for example, using a physics-informed neural network³⁸ or graph neural network,³⁹ can be thought of as a natural extension of the present work.

The general approach presented here can be extended to more complex tasks, such as three-dimensional multiphase or even reactive unit operations in chemical engineering. The limitation will be the availability of data or the associated computational cost of the underlying simulations.

This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—FR 4072/3—within the Priority Programme “SPP 2331: Machine Learning in Chemical Engineering.” The authors acknowledge the support received from the state of Baden-Württemberg through bwHPC. Parts of this work were performed on the HoreKa supercomputer funded by the Ministry of Science, Research and the Arts Baden-Württemberg and by the Federal Ministry of Education and Research.

The authors have no conflicts to disclose.

Y.K. and A.J.K. contributed equally to this paper.

A.S., P.F., and B.P.L. contributed to conceptualization, funding acquisition, and supervision. Y.K., A.J.K., and M.S. contributed to the investigation, software, and visualization. All authors contributed to the methodology, validation, and writing of the original as well as the final draft. Y.K. and A.J.K. contributed equally.

Yuri Koide: Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). Arjun J. Kaithakkal: Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). Matthias Schniewind: Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). Bradley P. Ladewig: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Supervision (equal); Validation (equal); Writing – original draft (equal). Alexander Stroh: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Supervision (equal); Validation (equal); Writing – original draft (equal). Pascal Friederich: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Supervision (equal); Validation (equal); Writing – original draft (equal).

The data that support the findings of this study are openly available in GitHub at https://github.com/aimat-lab/ChemEngML.⁴⁰ 

As mentioned in Sec. II B, we investigated the effect of shifting and flipping of structures on the CNN prediction of C_f and St. Each of the test data points, excluding the flat channels (i.e., 459 data points), were examined as follows: For each test data point, 32 shifted and flipped copies were generated. Within this set, 16 structures were systematically produced by shifting the original structure along the x-direction in a total of 16 distinct increments within $Mxs$ and $Mxe$⁠. The remaining 16 structures were obtained by creating flipped images of the previously shifted structures. Subsequently, the (C_f, St) values of these 32 structures were averaged and compared with their respective ground truth values as shown in Fig. 9. The error in mean prediction is within ±5% for the majority of structures in the case of C_f and for nearly all structures in the case of St. Similarly, the percentage error in mean prediction for (C_f, St) for structures A, B, and C is within a reasonable range and is estimated to be (−13.87, 4.04), (−15.17, −0.53), and (12.91, 2.85), respectively.