12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

MACHINE LEARNING FOR THE PREDICTION OF THE LOCAL DRAGFORCES AND HEAT TRANSFER RATES IN TURBULENT FLOWS PASTROUGH SURFACES

Rafael Diez Sanhueza
Faculty of Mechanical, Materials and Maritime Engineering
Delft University of Technology
Leeghwaterstraat 39, 2628 CB, Delft, the Netherlands
R.G.DiezSanhueza-1@tudelft.nl

Ido Akkerman
Faculty of Mechanical, Materials and Maritime Engineering
Delft University of Technology
Mekelweg 2, 2628 CD, Delft, the Netherlands
I.Akkerman@tudelft.nl

Jurriaan W. R. Peeters
Faculty of Mechanical, Materials and Maritime Engineering
Delft University of Technology
Leeghwaterstraat 39, 2628 CB, Delft, the Netherlands
J.W.R.Peeters@tudelft.nl

ABSTRACTof rough surfaces can substantially increase the drag resistanceof transportation systems, and lower the efficiency of thermodynamic cycles. One of the main challenges while working
with rough surfaces is to predict the impact of a given surface
topography on the drag resistance and heat transfer rates ofa system. Most traditional correlations are based on standard
surface metrics, such as the root-mean-squared height variations, skewness, kurtosis, effective slope, forward-facing angle, or different auto-correlation functions.
In this work, we develop a machine learning system basedon convolutional neural networks to predict the local behaviourof turbulent flows past irregular surfaces. These systems areable to describe the complex physical effects observed whenturbulent flows interact with rough walls. The complex nonlinear functions required to approximate the behavior of turbulence are automatically reconstructed using a database containing rough surfaces with unique groups of topological features producing different flow effects.

Turbulent flows past rough surfaces can create substantial energy losses in engineering equipment. During the last
decades, developing accurate correlations to predict the thermal and hydrodynamic behavior of rough surfaces has provento be a difficult challenge. In this work, we develop a convolutional neural network architecture to perform a direct imageto-image translation between the height map of a rough surface
and its detailed local drag resistance and heat transfer rates.
Various techniques are discussed to improve the computational
efficiency of the machine learning architecture proposed, andeven to reduce its time and space complexity. The main study
is based on a new DNS database formed by 24 flow cases at afriction Reynolds number of Reτ = 180 obtained by applying arandom shift to the Fourier spectrum of the grit-blasted surface
scanned by Busse et al. (2015,). The results show that machinelearning can accurately predict the global values of the dragresistance and heat fluxes across a rough surface. The local
predictions for both momentum and heat transfer also show aconsiderable improvement upon increasing the dataset size. A
detailed analysis of the global skin friction values and Stanton numbers predicted by deep learning further reveals thatthe results surpass the accuracy of traditional correlations by a
substantial margin in the dataset analyzed.

DNS DATABASE
The study is based on a new DNS database formed by
24 flow cases at Reτ = 180. Numerical simulations are usedinstead of experimental measurements, since the local dragforces and heat fluxes correspond to quantities of interest. The
flow cases are generated by applying a random phase shift
to the Fourier spectrum obtained from the grit-blasted surface
scanned by Busse et al. (2015,). The random shift is given by
the variable φi in the following equation for the height function
H(x, z):INTRODUCTION
Turbulent flows past rough surfaces can be found in alarge variety of engineering applications. Irregular surfaces
are often caused by external degradation processes, such as
bio-fouling, abrasion, machining, or corrosion. The presence1

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

-15.7%

-7.8%

0%

7.8%

15.7%

Figure 1. Examples of height maps generated by applying a random shift to the Fourier spectrum of the grit-blasted surface scanned
by Busse et al. (2015,). All percentages are with respect to the half-channel height δ = Ly /2.

 


x
z
H(x, z) = ∑ Ri cos 2π Mi
+ Ni
− φi
Lx
Lz
i

Tb+ = 11.21 respectively. The maximum differences observedwith respect to these quantities are 3.84% for Reb and 2.82%
for Tb+ . Therefore, all DNS cases employed similar mass
fluxes and other bulk flow parameters. The velocity shifts observed in the log-layer (∆U + ) for all DNS cases ranged from4.09 to 5.07. The latter is equivalent to variations in the dimensionless Nikuradse sand-grain roughness height (ks+ ) from24 to 33 according to the equations provided in the work of
Thakkar et al. (2018,). Therefore, all DNS simulations contained similar levels of turbulence activity.

(1)

The terms Ri , Mi , Ni correspond to constants extracted
from the original Fourier spectrum for the grit-blasted surface. Figure 1 presents a collection of 8 examples with roughsurfaces generated according to this technique. Despite thesimplicity of the methodology used, every DNS case contains
unique groups of topological features, which in turn create
flow fields with different characteristics. Therefore, predicting the local drag forces and heat fluxes present in each DNS
case corresponds to a challenging machine learning task.
The dimensionless Navier-Stokes equations and the energy equation solved for each DNS case are the following:∇ · u = 0,(2)

2
∂t u + u · ∇u = −∇p + Re−1
τ ∇ u+Sf ,(3)

2∂t T + u · ∇T = Pe−1τ ∇ T + Sq .

(4)

WALL FORCE INTERPOLATION
The local drag forces (F) and heat transfer rates (Q̇) acting over irregular surfaces are obtained by using a finite element integration scheme with face elements located at therough walls. The equations integrated are the following:

F=

¨ 
 
1 
−P n +∇u + ∇uT · n dA,
Reτ
A

Q̇ =
In eqs. (2-4), the variables u, p, T correspond to the dimensionless velocity components, pressure and temperature
respectively. The variable Peτ = Reτ Pr is the friction Peclet
number, which is defined using the molecular Prandtl number
Pr. The source terms S f and Sq for the momentum and energy equations are equal to unity in all DNS simulations, aswell as the Prandtl number Pr = 1. The velocity and ptemperature fields are scaled using the friction velocity uτ = τw /ρ
and the friction temperature Tτ = q̇w /(ρc p uτ ). Since the reference density ρ and specific heat capacity c p are equal toone, it can be further proven that uτ = Tτ = q̇w = 1 for all
the DNS cases considered. The numerical simulations are carried out using the in-house DNS code described by Peeters& Sandham (2019,). Identical discretization schemes are employed as in the original study. The friction Reynolds number (Reτ ) is set to 180. The domain size and the grid resolution have dimensions (Lx × Ly × Lz ) = (5.63 × 2 × 2.815) and(Nx × Ny × Nz ) = (280 × 280 × 140) in the streamwise, wallnormal and spanwise directions respectively. The scaling ratio
for the height of the rough surfaces (k/δ ) is kept at a constantvalue of 1/6, according to the definition of the original study.
In all the DNS cases generated, the bulk Reynolds numbers Reb = Ub Ly /ν and the dimensionless bulk temperatures(Tb+ ) presented nearly constant values of Reb = 4002 and

1
Peτ

(5)

¨
∇T · n dA.

(6)A

In eqs. (5-6), the variable A refers to the area of the roughsurfaces in contact with the fluid. The vector n corresponds tothe normal of the rough surfaces:−1∂ H 2 ∂ H 2
∂H T
∂H
1+
, 1, −
+
−
. (7)
∂x
∂z
∂x
∂z

s
n=

The orientation of eq. (7) is intended to work with roughsurfaces located at the bottom side of channel flows, where
the wall-normal direction corresponds to the y-axis. For roughsurfaces present at the top side of channel flows, the orientationof the vector n must be adjusted according to the transformations considered for the respective height function H(x, z). The
differential area dA for the rough surface elements is given by:
s
dA =

1+

∂H 2 ∂H 2+
dx dz.
∂x∂z(8)Afterqcombining eqs. (5-8), it can be shown that the expression 1 + ∂ H/∂ x2 + ∂ H/∂ z2 cancels out. Therefore, if

212th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
Depthwise
Convolution
3×3

BN + PReLU

Pointwise
Convolution
1×1

BN + PReLU

BN + PReLU

Figure 2. Schematic representation of the depthwise separable convolution (DSC) modules employed in the machine learning study.
The abbreviation BN refers to the 2-D batch normalization operations applied before each PReLU activation function.
Table 1. Order of the polynomial terms considered by theshape functions in every direction.

polynomial shape functions are used to represent the variables(u, P, T, H), then the final integrals given by eqs. (5-6) containpolynomial terms exclusively. Thus, a Gauss-Legendre integration scheme with a sufficient number of intermediate points
can be used to find the analytical solution of the resulting integrals for the drag forces and heat fluxes.
In order to perform finite element integration, the variables (u, P, T, H) were reconstructed using shape functions
based on a symmetric stencil of data points surrounding eachfinite element. The degree of the polynomial terms chosen
in each direction can be found in Table 1. Here, it can benoted that all interpolation schemes considered linear terms in
the wall-normal direction (y), whereas mixed-order linear or
quadratic terms were considered in the horizontal directions(x − z). The latter was necessary because the DNS solver employed a staggered discretization scheme to handle the velocity components u = (ux , uy , uz ). Therefore, the number of data
points required to form a symmetric stencil was different ineach direction. The usage of linear interpolation schemes in
the wall-normal direction was validated empirically, since it
was noted that quadratic terms yielded slightly higher errors
with respect to both the global force and heat transfer balances.
All FEM integration areas considered for eqs. (5-6) were centered around the original x − z coordinates for the pressure andtemperature fields (P, T ). The pressure field was extrapolated
in the wall-normal direction (y) by considering the first two
data points available above the rough surfaces. However, thevelocity and temperature fields (u, T ) considered that the firstlayer of data points was located exactly over the rough surfaces
at the corresponding H(x, z) locations, with a value equal to thehomogeneous boundary conditions (u = T = 0).
The final numerical implementation was written in PyTorch, since it can be proven that the entire numerical procedure can be expressed as a sequence of parallelizable array
operations. The process of gathering data points located inneighboring x − z locations can be expressed using standard
array shift operations due to the presence of periodic boundary
conditions for the rough surfaces. The linear systems of equations required to fit separate shape functions for the flow variables (ux , uy , uz , P, T ) surrounding each FEM integration area
can be solved using 4-D tensors containing batches with small
systems of equations.

Variable

Streamwise(x)

Vertical(y)Spanwise
(z)

P

2

1

2

T

2

1

2

ux

1

1

2

uy

2

1

2

uz

2

1

1

H

2

-

2

the input height maps H(x, z) of the rough surfaces.
The final machine learning study is based on a deep learning architecture formed by depthwise separable convolution
modules (Chollet, 2017). The overall shape of the ML architecture is presented in Table 2, whereas the details of eachconvolution module are described in Figure 2. Every module has an internal dimension of 20 channels, 3×3 depthwise
convolutions, and three PReLU activation functions (He et al.,2015) preceded by batch normalization. The usage of ReLU
operators is avoided during the machine learning study, since
PReLU neurons increase the expressivity of neural networks
with a minimal computational cost, and they can facilitate theconvergence of the training procedure. Max-pooling operators
are also avoided because they lose the exact spatial location of
the features being processed. Every machine learning model
created has a total of 5,780 trainable parameters.
The training procedure for each machine learning configuration is driven by L1 loss functions instead of traditional L2
metrics. This choice helps the optimizer to train models that
focus on predicting the global values for the drag forces andheat transfer rates with greater precision than L2 loss functions, which tend to over-penalize small regions with outlier
points. Physics-informed data augmentation is performed bymirroring all DNS cases with respect to the spanwise axis. The
final study is based on the leave-one-out cross-validation technique, which is beneficial for small datasets (Wong, 2015).
The overall time and space complexity of the machinelearning system can be substantially improved by creating a
direct image-to-image translation system between the input
height map H(x, z) for a DNS case and the output predictions
for the local drag forces fx (x, z) or heat transfer rates q̇(x, z).
A schematic representation of this approach can be found in
Figure 3. Here, it can be noted that the standard convolutional
neural network shown in Figure (3.a) systematically reduces
an input image until a single output prediction is obtained forfx (x, z) or q̇(x, z), discarding all intermediate information. Thisapproach is highly inefficient while generating predictions for

MACHINE LEARNING
In order to predict the local drag forces ( fx ) and heat
fluxes (q̇) acting over irregular surfaces, it is necessary to build
advanced non-linear models based on techniques such as machine learning. A preliminary assessment revealed that linear
regression models presented a low degree of correlation between the input height map H(x, z) for a DNS case and theobserved ground-truth data for fx (x, z) or q̇(x, z). The average
L1-errors with respect to the total drag forces and heat transfer
rates reached magnitudes up to 114% and 42.2% using linear
models respectively. Therefore, it was decided to build predictive systems based on convolutional neural networks to scan

3

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
Standard System
Convolutional Stridea)

Efficient System
Convolutional Dilation

b)

1
2
1
1
2
1
1
Input Image

Figure 3.

4
2
2
2
1
1
1
Input Image

Comparison in 1-D between a standard convolutional neural network and the efficient computer vision system developed.
Table 2. Layers present in the deep learning architecture.
The columns Cin , Cout , D and AF&B refer to the numberof input channels, output channels, dilation and the presenceof activation functions with bias respectively. The column
Input Image refers to the equivalent input image size considered by each convolutional layer, with respect to a total image
size of 280×140 for the entire DNS domain. The final layer
Conv. 2D corresponds to a simple convolutional filter with akernel size of 3×3. The architecture of the DSC modules ispresented in Figure 2.

neighboring (x, z) spatial locations sharing a large number of
input pixels in their local height maps, since the same convolutional operations must be repeated again. In contrast, theadvanced computer vision system shown in Figure (3.b) avoids
redundant operations, since the results of intermediate convolutional layers are computed only once, and the upcoming layers share information by utilizing a modified dilation systemthat produces the same effect as applying reduction strides ina traditional convolutional neural network. As a result, an array with all output predictions for fx (x, z) or q̇(x, z) for a DNS
case is computed after only one pass through the neural network. The empty rectangles shown in Figure (3.b) correspondto information that is not connected to the output prediction
highlighted at the top, but is used to process other neighboringdata points. In terms of computational efficiency, the approachshown in Figure (3.b) is able to reduce the baseline time complexity required to process an input height map H(x, z) discretized as an image with N pixels from quadratic time O(N 2 )to linear time O(N). Therefore, the running times using theapproach described in Figure (3.b) are orders of magnitude
lower than traditional computer vision architectures, which
process one image for every data point predicted. Moreover,the methodology described in Figure (3.b) can also reduce thespace complexity and RAM memory requirements of the machine learning system during both the training and evaluation
stages, since only one image is processed per DNS case. Finally, circular padding is applied to the convolutional architecture presented in Figure (3.b) to account for the periodic
boundary conditions considered for the rough surfaces.

Layer

Cin

Cout

D

AF&B

Input
Image

DSC

1

20

1

Yes

119×61

DSC

20

20

2×1

Yes

59×59

DSC

20

20

4×2

Yes

29×29

DSC

20

20

4×2

Yes

27×27

DSC

20

20

8×4

Yes

13×13

DSC

20

20

8×4

Yes

11×11

DSC

20

20

16×8

Yes

5×5

Conv. 2D

20

1

16×8

-

3×3

12 DNS cases show mixed trends, since the errors for the local
drag forces are similar to the study using 24 DNS cases, yet theerrors for the local heat fluxes only present moderate improvements. The relatively large errors obtained for the momentum
predictions, compared to the heat fluxes, are partially causedby the fact that the local drag forces reach magnitudes up to 23
times larger than the average values for the DNS data, whereas
the largest heat fluxes only reach values up to 8.1 times higher.
Therefore, minor differences in areas where strong drag forces
are expected will create much greater errors than similar percentual differences in regions where large heat fluxes occur.
In Figure 4, the streamwise patterns observed in the errors for the heat fluxes are likely caused by the reduced scopeof the input images considered to perform predictions. While
it is ideal from a physical perspective to build a predictive systems considering the largest input images available, the usageof images with more pixels can lead to neural networks withan increased number of parameters, and the associated risk ofover-fitting. Therefore, a trade-off must be considered between
the size of the input images, and the increased number of parameters that the neural network will require to process more

RESULTSThe results of the machine learning study can be found in
Figures 4 and 5. In Figure 4, a comparison is presented between the local drag forces and the heat fluxes obtained by the
deep learning system with respect to the original DNS data.
As it can be seen in Figure 4, the machine learning predictions show physically realistic trends, with the majority of theerrors located in areas where large spatial variations in the local drag forces or heat transfer rates are expected. A detailed
analysis of the global errors obtained for the force and heattransfer balances revealed that most of the errors were lower
than 5% during the cross-validation study, whereas the maximum errors recorded were 13.21% for the drag forces and5.99% for the heat fluxes. Both of these results can be regarded as satisfactory. Figure 5 presents histograms comparing the distributions of the local errors obtained using datasetswith 7, 12 and 24 DNS samples. These results reveal that theerrors obtained for the momentum and heat transfer predictions present a substantial improvement upon increasing thedataset size from 7 to 24 DNS samples. The predictions using

4

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
Drag Forces - Worst Case

Heat Transfer - Best Case

Heat Transfer - Worst Case

Difference

MachineM
Learning
L

DNS

Drag Forces - Best Case

-17.5

-8.7

0

8.7

17.5

-12.5

-6.2

0

6.2

12.5

0

1.4

2.8

4.2

5.6

0

1.7

3.4

5.1

Figure 4. Comparison between the local drag forces and heat transfer predictions obtained by the deep learning system with respect
to the DNS data. The cases presented were selected based on the highest and the lowest values for the L1-norm of the local errors. The
values of each colormap indicate the ratio between the local quantities and the average drag forces or heat fluxes present in the DNS
data.
input pixels. The influence of the input images becomes moreevident for the heat transfer predictions than for momentum,
since the higher accuracy of the deep learning system implies
that the assumptions established during the modelling process
start to become noticeable. Another important insight is thatthe results obtained during the worst-case scenario for the momentum predictions originated from a rough surface which had
yielded favorable results while using 12 DNS samples insteadof 24 cases. Therefore, it can be concluded that the momentum predictions present signs of over-fitting, and that further
research is required to determine if the errors can be reduced
by incorporating more training data, or if the architecture of
the deep learning model can be further optimized.
After obtaining the deep learning predictions, a study was
conducted to analyze the trends observed for the skin friction
factors (C f ) and the Stanton numbers (St). Among the traditional correlations proposed in the literature to predict the skinfrictions factors (C f ) or the equivalent Nikuradse sand-grainroughness (ks,eq ), it was determined that the correlation developed by Flack & Schultz (2010,) yielded accurate predictions:ks,eq = 4.43 Sq (1 + Ssk )1.37 .

Dipprey & Sabersky (1963,). As it can be seen in the subfigures, the deep learning system is substantially more accurate
than traditional correlations while predicting both momentum(C f ) and heat transfer (St) parameters. The maximum errors
for the skin friction coefficients C f in the deep learning systemreached values of 13.21%, whereas the system with traditionalcorrelations reached values up to 15.61%. Therefore, it can beconcluded that deep learning constitutes a valid alternative togenerate improved predictions for flow parameters, such as C f
or St, if enough training data is collected.

CONCLUSION
This study presented a deep learning architecture capable of predicting detailed maps for the local drag forces andheat fluxes acting over irregular surfaces. The results showthat machine learning is able to achieve reliable results while
predicting the global force and heat transfer balances. A sensitivity study with respect to the dataset size also revealed asignificant reduction in the errors for the momentum and heattransfer predictions once the dataset size is increased from 7
to 24 DNS cases. The comparisons performed with respect totraditional correlations proved that deep learning is a valid alternative to generate improved predictions for important flow
parameters, such as the skin friction factors (C f ) or the Stanton
number (St). Therefore, it is recommended to perform further
research regarding the creation deep learning models to predict
the behavior of turbulent flows past rough surfaces.

(9)

In eq. (9), the variables Sq and Ssk correspond to the rootmean squared height variations and the surface skewness respectively. A detailed review of the numerical procedure to
calculate these variables can be found in the work of Thakkar
et al. (2017,). The system of equations to calculate the skinfriction factor (C f ) is completed using the formulas describedby Peeters & Sandham (2019,), and the non-linear function∆(U + ) = f (ks+ ) provided by Thakkar et al. (2018,). The values for the Stanton number (St) are obtained by applying thecorrelation proposed by Dipprey & Sabersky (1963,). Theseequations establish a methodology to predict the Stanton number (St) of a rough surface based on a skin friction factor (C f )previously calculated. In Figure 6, the histograms present
the errors obtained for C f and St based on the deep learning
predictions and the traditional formulas previously described,
which are based on the work of Flack & Schultz (2010,) and5

6.9

Frequency

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
20%

20%

15%

15%

10%

10%

5%

5%

0%

N=7
N = 12N = 24

0%
0%

25%
50%
75%
Absolute Errors, Local Forces

100%

0%

7.5%
15%
22.5%
Absolute Errors, Local Heat Transfer

30%

Figure 5. Changes in the distribution of the absolute errors for the momentum and heat transfer predictions with respect to the dataset
size. All errors are scaled using the average DNS values for the drag forces or heat transfer rates.

Freq. Density

15

Machine learning
Traditional correlations

20

10

10

5
0

0
0%

2.5%

5%

7.5% 10%Error C f

12.5%

16%

0%

3%

6%

9%
12%
Error St

18.5%

Figure 6. Histograms comparing the distribution of errors for the skin friction values C f and the Stanton numbers St between the deeplearning model and the traditional correlations of Flack & Schultz (2010,) and Dipprey & Sabersky (1963,).

REFERENCES

agenet classification. In IEEE International Conference on
Computer Vision, pp. 1026–1034.
Peeters, J.W.R. & Sandham, N.D. 2019, Turbulent heat transfer in channels with irregular roughness. International Journal of Heat and Mass Transfer 138, 454–467.
Thakkar, Manan, Busse, Angela & Sandham, Neil 2017, Surface correlations of hydrodynamic drag for transitionallyrough engineering surfaces. Journal of Turbulence 18 (2),
138–169.
Thakkar, M., Busse, Angela & Sandham, Neil 2018, Direct
numerical simulation of turbulent channel flow over a surrogate for nikuradse-type roughness. Journal of Fluid Mechanics 837.
Wong, Tzu-Tsung 2015, Performance evaluation of classification algorithms by k-fold and leave-one-out cross validation. Pattern Recognition 48 (9), 2839–2846.

Busse, A., Lützner, M. & Sandham, N. D. 2015, Direct numerical simulation of turbulent flow over a rough surface basedon a surface scan. Computers & Fluids 116, 129–147.
Chollet, F. 2017, Xception: Deep learning with depthwise separable convolutions. In IEEE Conference on Computer Vision and Pattern Recognition, pp. 1800–1807.
Dipprey, D.F. & Sabersky, R.H. 1963, Heat and momentum
transfer in smooth and rough tubes at various prandtl numbers. International Journal of Heat and Mass Transfer 6 (5),329–353.
Flack, Karen A. & Schultz, Michael P. 2010, Review of Hydraulic Roughness Scales in the Fully Rough Regime. Journal of Fluids Engineering 132 (4), 041203.
He, K., Zhang, X., Ren, S. & Sun, J. 2015, Delving deep
into rectifiers: Surpassing human-level performance on im-

6