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

Numerical simulation of fluids plays an essential role in modeling many physical phenomena, which enables technological advancements, contributes to sustainable practices, and expands our understanding of various natural and engineered systems. The calculation of heat transfer in fluid flow in simple flat channels is a relatively easy task for various simulation methods. However, once the channel geometry becomes more complex, numerical simulations become a bottleneck in optimizing wall geometries. We present a combination of accurate numerical simulations of arbitrary, flat, and non-flat channels and machine learning models predicting drag coefficient and Stanton number. We show that convolutional neural networks (CNN) can accurately predict the target properties at a fraction of the time of numerical simulations. We use the CNN models in a virtual high-throughput screening approach to explore a large number of possible, randomly generated wall architectures. Data Augmentation was applied to existing geometries data to add generated new training data which have the same number of parameters of heat transfer to improve the model's generalization. The general approach is not only applicable to simple flow setups as presented here but can be extended to more complex tasks, such as multiphase or even reactive unit operations in chemical engineering.


I. INTRODUCTION
Heat transfer in fluid flow is an important physical phenomenon, with relevance across all areas of science and engineering ranging from microfluidic devices in chemical engineering and biomedical implants, all the way to high-temperature physics and cosmology.In this proofof-concept study, we explore an interesting engineering question, which could be posed as "Is it possible to introduce structural changes to the wall of a channel that increases heat transfer, without a corresponding increase in the pressure drop?".This fundamental question linked to the ultimate goal of dissimilar flow control or dissimilar heat transfer enhancement 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 2 .Investigations of various surfaces including specially designed fins 18,23 , dimples 7 or vortex generators 9 report that an increase in heat transfer (described by Stanton number St) is always accompanied by inevitable manifold increase in the drag coefficient C f , which eventually results in a decrease of the Reynolds analogy factor RA = 2St/C f 29 in comparison to a flat channel configuration.It is, however, known that a dissimilar modification of momentum and heat transfer is possible when more 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 13,15,16,26 .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,31 .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 30 .The modification of RA is however limited in this case to several percent due to the increase of the wetted area and the corresponding increase in C f .
To simplify the scenario, we consider a twodimensional channel with laminar flow, heat transfer, and immersed boundary method for the introduction of surface structuring.This allows to quickly execute direct numerical simulations (DNS), where a large set of arbitrarily generated surfaces can be investigated.In this proof-of-principle study, we present a workflow consisting of numerical simulations (Section II A and Section II B), the generation of a dataset from numerical simulations, and the training of machine learning (ML) models (Section II C).The "hybrid" approaches can be used to explore heat transfer enhancement efficiently compared to only doing numerical simulations.The idea of using "pure" ML is to replace the entire Navier-Stokes simulation with approximations based on deep neural networks. 3,19,22To utilize this cost-efficient approach, precise calculation of geometries generation for ML training is needed.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 screen a large database of possible channel geometries (Section III).We believe that fluid dynamics is central to transportation, health, and defense systems, and it is, therefore, essential that ML solutions are interpretable, explainable, and generalizable.

A. Numerical procedure
For the problem setup, we consider a laminar channel flow with arbitrary wall structuring.The coordinate system of the numerical domain and its geometry (L x × L y = 10δ × 2δ with δ being the half channel height) are illustrated in Figure 1, where (x, y) = (x 1 , x 2 ) correspond to the streamwise and wall-normal directions respectively.The velocity components in the two directions are denoted by (u, v) = (u 1 , u 2 ).The analysis is carried out using flow and temperature fields produced by DNS in a channel flow driven at a constant flow rate (CFR).Assuming an incompressible flow, the velocity field is required to satisfy continuity: and the Navier-Stokes equations for a constant property Newtonian fluid: Here p is the fluctuating pressure part, ρ is density, ν is the kinematic viscosity and F IBM ,i represents the external volume force per unit mass required for the immersed boundary method (IBM) with which the structured surface is introduced into the flow domain 11 .In the present configuration F IBM,i corresponds to the frictional drag between the flow and the part of the surface reproduced by the immersed boundary method, i.e. the structured wall surface.P x is the absolute value of the mean streamwise pressure gradient added to the equation in order to drive the flow through the channel.
Due to the CFR approach the bulk Reynolds number is fixed to Re b = 2U b δ/ν = 200 for all considered simulations, where U b 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 wall (y = 0, 2δ).Subscript l and u are used throughout the manuscript to denote quantities on the lower and upper walls, respectively.
Temperature T is treated as a passive scalar and has to satisfy the scalar transport equation: where α denotes the thermal diffusivity.Periodic boundary conditions are applied for the thermal field in xdirection, while constant temperature is prescribed on both the lower and upper walls of the flow domain.
The non-dimensionalized temperature is defined as θ = (T − T l )/∆T w with ∆T w = T u − T l , such that θ l = 0 and θ u = 1.The Prandtl number is chosen to be Pr = ν/α = 1.Q IBM is proportional to the heat transfer rate between the flow and the structured wall and can be considered as a counterpart to the volume force F IBM,i in the momentum equation.This term is adjusted to fulfill the temperature boundary condition on the structured wall.Due to the use of periodic boundary condition for temperature, the absolute value of the heat transfer rate on the two walls should be identical once the solution reaches the thermal equilibrium.For the same reason, the mean heat flux in the wall-normal direction is constant in the channel.The present thermal boundary condition is chosen following other studies of heat transfer above structured walls 21,25,27 .The solver implementation is based on a spectral solver for incompressible boundary layer flows 5 .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 384 × 129 grid nodes.

B. Performance indices
Contrary to the laminar flow in a flat channel no universal analytical solution can be derived for a channel with arbitrary structuring at both channel walls.Utilizing the melt-down heights of the imposed structure for both walls (h u , h l ) and splitting the flow into two halves based on the position of the maximal spatially averaged velocity denoted with y c (Fig. 1), the balance between pressure drop P x and the average effective wall shear stress τ ef f is given by where δ u and δ l define the upper and lower effective channel half heights with respect to y c .Based on the wall θ shear stress the mean drag coefficient is given as where the effective bulk mean velocity The brackets ⟨⟩ denote a quantity averaged in x-direction so a split-up into the mean part ⟨ϕ⟩ (y) and spatial fluctuation part ϕ ′ (x, y) can be performed for any quantity ϕ(x, y): Due to the asymmetry in the temperature boundary condition, the heat transfer properties have to be separately evaluated for each wall.Hence, the hydraulic diameter is defined for the upper and lower walls as The Nusselt number for both walls can be estimated with where q tot denotes the total heat flux and the bulk mean temperature differences are defined as ∆θ b,l = 1 and ∆θ b,u = 1 The average of Nu l and Nu u is computed to determine the resultant Nusselt number of a particular case.The effective bulk mean velocity for each channel half is given by The total heat flux q tot , which is a constant as mentioned previously, can be estimated as the sum of where the three terms are respectively named the laminar, total fluctuation, and IBM contributions 30 .Here c p denotes the specific heat capacity and Reynolds analogy factor RA relates Stanton number to the drag coefficient and is used to evaluate the similarity between drag coefficient and heat transfer 4 .An increase in RA highlights a stronger enhancement in heat transfer compared to that in the drag coefficient and hence is desirable in the design of an energy-efficient fluidic system.It has to be noted that for the chosen boundary conditions RA = 0.533 with St = 0.016 and C f = 0.06 in the flat channel configuration.

C. Dataset and machine learning model
To generate a diverse dataset of wall structuring, we used a random walk algorithm combined with spline interpolation and discretization on the simulation grid.Each wall structure consists of n supporting points between a start and an end point at x = 0 and x = 384 with the same y-position for periodic boundary conditions.The x-coordinates of the supporting points are sampled from: With y ′ being the previous y position, σ y varied from table I and ∆ the available build space in y-direction.A minimum of 50% of the channel height is kept empty for the flow.To allow for larger meanders, the generation algorithm of the first wall surface can use the full 50% of build space, so here ∆ remains at constant 64.5 at all x positions.For the second channel surface, ∆ is adjusted according to the first channel surface.The first y-coordinate at x = 0 to initialize y ′ is drawn from a uniform distribution in the interval [0, ∆].
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 r's also the radius of the curves can get very small, resulting in sharp features.The parameter a controls the smoothness of the curve.For a = 0 the angle through one supporting point is determined by the mean of the directions to both neighbouring points.At higher a the direction to one neighbouring point is weighted higher and hence the curve features stronger edges 28 .
Those obtained curves were then projected on the dealiased grid.While the x-coordinates where linearly spaced, the y-coordinates of the grid nodes were obtained using: With the integer i in the interval [0, 128].Using 20 random initializations per parameter variation 108, 000 random wall structures were generated, hereafter called the repository set.
We calculate the drag coefficient C f and Stanton number St using the simulation method described above for a subset of 10, 800 randomly sampled channel geometries.From that, 9, 185 passed a set of filters regarding temperature convergence and geometric validity.5 % (459) of this data was put aside as the test set.
We used varying fractions of the remaining 8, 726 channel geometries (hereafter called the labeled set) for hyperparameter optimization of CNN.Hyperparameters were determined without flat channel geometries.The training and test sets include flat channels to boost geometry To account for periodical boundary condition and invariance of vertical orientation the channel geometry was randomly flipped about the y-axis, and also randomly shifted in the xdimension within the range of x = 0 to 384.The model was implemented using TensorFlow 1 and Keras 6 , and training was done using the Adam optimizer 20 .The model is trained for 100 epochs with a batch size of 256 using the mean squared error (MSE) as loss while logarithmically reducing the learning rate from 1e −3 at the 10 th epoch to a varying final learning rate at the last epoch.
To find out whether the prediction of C f and St can be improved jointly this search space (table II) is optimized using multi-metric Bayesian optimization.For this SigOpt is employed 14,24 with 140 experiments on 8 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 5-fold cross-validation are used as metrics.

III. RESULTS
We trained the CNN model described above and obtained MAE (mean absolute errors) (and r 2 -scores) of M AE = 1.90 • 10 −3 (r 2 = 0.951) and M AE = 1.29 • 10 −4 (r 2 = 0.925) for predictions of C f and St.A comparison of CNN predictions and simulated ground truth on the test set is shown in Figure 2. Most values of C f and St are well correlated.However, deviation for large values can be seen(C f = 0.12 ∼, St = 0.020 ∼).
In order to evaluate how well the CNN performs on smaller datasets, we generated learning curves (see Figure 3), 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, 9,240 channel geometries.We trained the CNN model with each data set and evaluated the model performance by plotting mean absolute errors.We observe an exponential decrease in the mean absolute error with the increase in training set size.It can be observed from Figure 3 that larger datasets increase the model accuracy.In C f case, the MAE score for all data is approximately twice smaller than that of 5 %, meaning a large number of geometry channels is necessary to predict C f and St with high accuracy.In St case, a linear decrease as the size of the dataset is larger can be seen from the plot.As the plots continue to decrease exponentially, further improvement in the accuracy of the model can be expected by increasing the number of data.
The scatter plot of 8, 726 labeled data points from the repository set used for training the CNN is shown in Fig- ure 4a.Deviations from a flat channel (C f,ref = 0.06 and St ref = 0.016) necessarily lead to an increase in C f that outweighs the simultaneous increase in St.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 these structures that are highlighted in Figure 4a are discussed in detail in Section IV.We then exploited the speed up of the surrogate machine learning model (< 100 ms per channel) compared to the numerical simulation (≈ 20−30 min per channel).This allowed us to explore the flow and heat transfer characteristics (C f and St) of a much larger set of unlabeled channel geometries from the repository set.A scatter plot of C f vs St for the repository set, split up into labeled and unlabeled data points, is shown in Fig- ure 4b.In general, it is very difficult to find geometries with RA value exceeding that of a flat channel as can be seen in Figures 4a and 4b.The histogram representation of the repository set shown in Figures 4c and 4d shows that the training data adequately represents the repository set.

IV. DISCUSSIONS
In this section, we analyze in detail the pressure loss and heat transfer characteristics of the three labeled structures highlighted in Figure 4a (namely structures  A, B, and C).The mean flow and temperature fields of these structures are shown in Figure 5, 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 is 25.1%, 5.4%, and 6.3% respectively.
Among the three structures, 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 8 .Though there are two vortices inside the flow field as can be seen in Figure 5a, 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 Figure 5b) result in an 80% increase in heat transfer.At the same time, the shape of the structure-A together with its increased surface area introduces significant pressure loss leading to a 60% reduction in RA when compared to the flat channel.
For the structure B, flow inside gets accelerated inside the converging section followed by fluid deceleration inside the diverging section.Unlike the case of structure A, a developing thermal boundary layer is present only along the bottom wall.The presence of the recirculation region, covering the entire bottom wall except for the converging-diverging section, further limits the 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 convergingdiverging 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 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 wallnormal 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 the increased heat transfer.This suggests the possibility to achieve a value greater than unity for the RA value with carefully selected parameters for wall meandering 10 .

V. CONCLUSIONS AND OUTLOOK
We presented a combination of accurate numerical simulations of fluid flow and heat transfer in arbitrary, nonflat channels and machine learning models predicting drag coefficient C f and Stanton number St. We found that C f and St are well predicted from channel geometries by the CNN model with data augmentation.However, prediction is limited for the complex geometries.The higher numbers of C f and St are difficult to predict the ground truth numbers.We show that once trained the CNNs can predict the target properties at a fraction of the time (< 100 ms per channel) required by numerical simulations (≈ 20 − 30 min per channel).This can be exploited for exploration and optimization tasks 17 .
The general approach is not only applicable to simple flow setups as presented here but 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.Since the current CNN is trained for a spe- cific set of flow conditions, in terms of Reynolds number and boundary conditions, it will be interesting to know how to modify the current model to cover a wide range of these conditions.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 channel structures and directly solve the inverse problem, i.e. the suggestion of channel architectures given desired target properties.Also, trying to predict not only 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, can be thought of as a natural extension of the present work.

FIG. 2 :FIG. 3 :
FIG. 2: Predictions of the CNN model of a) drag coefficient C f and b) Stanton number St compared to the ground truth on the validation set.
FIG. 4: Scatter plot of C f vs St for the a) training (labeled) and b) repository set.Histogram representation of the training and repository set for c) C f and d) St.

FIG. 5 :
FIG. 5: Mean velocity and temperature fields for the structures A, B, and C. Streamlines are added on top of the velocity fields for a better understanding of the flow field, especially the recirculation regions.The white dashed lines in the temperature fields indicate θ = 0.10 and 0.90 and are representative of the thermal boundary layer.

TABLE I :
Parameters for the structure generation algorithm.istheith supporting point in the interval[1, n].σ x is varied according to table I.The y-coordinate at each respective x is sampled from: i

TABLE II :
Parameters and final optimum for the hyperparameter optimization.
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.We trained the CNN model with augmented channel geometries to improve the model's interpretability for different shapes of geometries.The final hyperparameter optimum was used to train the model.Data augmentation was implemented on each batch and epoch during training.

TABLE III :
C f , St, and RA values for the structures A, B, and C. Also shown are the values normalized with the corresponding values for the flat channel, i.e.C f ref = 0.06 and St ref = 0.016.