On the  Existence of Logarithmic Terms in the Drag Coefficient and   Nusselt Number of a Single Sphere at High Reynolds Numbers
1,2]Yousef M. F. El Hasadi[Email address for correspondence: G.DamianidisAlChasanti<EMAIL_ADDRESS>1] Johan T. Padding[Email address for<EMAIL_ADDRESS>[1]Process and Energy department, Delft University of Technology,Leeghwaterstraat 39,2628 CB Delft, the Netherlands
[2] Civil Engineering and Geosciences Department,Delft University of Technology, Stevinweg 1, 2628 CN Delft, the Netherlands

At the beginning of the second half of the twentieth century, Proudman and Pearson (J. Fluid. Mech.,2(3), 1956, pp.237-262) suggested that the functional form of the drag coefficient ($C_D$) of a single sphere subjected to uniform fluid flow consists of a series of logarithmic and power terms of the Reynolds number ($Re$). In this paper, we will explore the validity of the above statement for Reynolds numbers up to $10^{6}$ by using a symbolic regression machine learning method. The algorithm is trained by available experimental data and data from well-known correlations from the literature for $Re$ ranging from $0.1$ to $2\times 10^5$. Our results show that the functional form of the $C_D$ contains powers of $\log(Re)$, plus the Stokes term, fulfilling partially the statement made above. The logarithmic $C_D$ expressions can generalize (extrapolate) beyond the training data and are the first in the literature to predict with acceptable accuracy the rapid decrease (drag crisis) of the $C_D$ at high $Re$. We also find a connection between the root of the $Re$-dependent terms in the $C_D$ expression and the first point of laminar separation. We did the same analysis for the problem of heat transfer under forced convection around a sphere and found that the logarithmic terms of $Re$ and Peclect number $Pe$ play an essential role in the variation of the  Nusselt number $Nu$. The machine learning algorithm independently found the asymptotic solution of Acrivos and Goddard (J. Fluid. Mech., 23(2),1965, pp.273-291).

§ INTRODUCTION

Predicting the drag force on an object fixed in a planar flow has been the  subject of extensive investigation from the early days  of fluid mechanics when it emerged as an independent discipline. The analytical solution for the drag force experienced by a  rigid sphere for  creeping flow conditions, found by Stokes [1] in 1851, is one of  the first known analytical expression in the fluid mechanic's community. Stokes assumed in his solution that  inertial effects of the fluid could be neglected throughout the solution domain. However, Oseen
[2] found an inconsistency in the Stokes solution. Specifically, he found that inertial fluid effects cannot be neglected far away from the sphere.  He derived a new form of  equations, known as Oseen equations [2],  that can handle this inconsistency, and he came up with an improved approximation  for the drag coefficient, defined as  $C_D = F_D/(\frac12 \rho v_{\infty}^2 \frac{\pi}{4} d^2)$,
where $F_D$ is the drag force, $\rho$ the fluid density, $v_{\infty}$ the fluid flow velocity far away from the sphere, and $d$ the sphere diameter
[3]. There are additional solutions to the  Oseen equations,  such as those of Goldstein
[4] and Faxén

Proudman and Pearson[6] and Kaplun and Lagerstrom [7] used the matched asymptotic method to solve the Navier-Stokes equations to resolve the fluid flow around different blunt bodies.
Proudman and Pearson [6] divided the flow field around the sphere into two stream function expansions. The first one, which they called the  Stokes expansion,  controls the flow near the surface of the sphere. The  second expansion, which they called  the  Oseen expansion,  controls the flow far from the surface of the sphere. Both expansions are based on the Navier-Stokes equations, and the two expansions are matched at a certain distance from the sphere using the method of matched asymptotics. Evaluating stresses from the Stokes expansion they arrived at  the following expression for the $C_D$ of a sphere:

\begin{equation} \label{eq1}
C_D= \frac{24}{Re}(1+\frac{3}{16}Re+\frac{9}{160}Re^2\log(\frac {Re}{2}))
\end{equation}

Here $Re = \rho v_{\infty} d/\mu$  is the Reynolds number. They  made the following statement (conjecture) about the expansions that govern the flow field
[6]: “The non-linearity of the Navier-Stokes equation then shows that both expansions must involve powers of $\log(Re)$, and it seems reasonable to suppose that both expansions are in powers of $Re$, each term of which is multiplied by polynomial in $\log(Re)$”. This statement also reflects on the functional form of the drag coefficient. However, the authors did not mention the $Re$ range for which the statement is valid. From now on, we  will call this conjecture   P&P.
Graebel  [8] supported the P&P statement by mentioning that the $C_D$ functional form that will result from  asymptotic expansions of the Navier-Stokes equations will always be a function of $\log(Re)$. A few years later, Chester et al.[9]  added an extra term to Eq.(<ref>), which was the last addition that came from the expansion of the Navier-Stokes equations.

The appearance of  logarithmic terms (alternatively  known as logarithmic switchback
terms [10]) in the asymptotic expansions have intrigued the scientific community, because in some instances  they were not   forced by the  governing equations [11] .
Van Dyke [12] dedicated  a section in his book describing the proliferation of logarithmic terms in different fluid mechanics problems, and he made the following comment: “one can philosophize that description by fractional powers fails to exhaust the myriad phenomena in the universe, and logarithms are the next simplest function". Initially,   the logarithms were tied with paradoxes in  fluid mechanics, or to  the singular perturbation techniques themselves. However, Lagerstrom and Reinelt [10]  showed that  logarithmic terms are part of the  solution of the governing equations, and the asymptotic expansion  method is just one way to reach to the solution. This view is supported by other investigations using different mathematical methods [13, 14].

There are  analytical solutions for the Stokes and Oseen regimes for some non-spherical particles such as oblate or prolate  spheroids, circular cylinders  and few other particle geometries
[15, 16, 17, 18]. Eq.(<ref>) and all  other analytical solutions, regardless of the  shape of the particles,  are valid up to  $Re\approx$ 1.0.

For higher $Re$, analytical solutions for the Navier-Stokes equations cease to exist due to its non-linearity. The flow around a sphere at high $Re$ consists of a mosaic of different flow morphologies, depending on  $Re$ as described by Achenbach [19], and  Kambel and Girimaji [20]. High Reynolds number flows ($Re\geq 10^4$) are usually classified into four flow regimes. In the subcritical flow regime, the $C_D$ value is independent of $Re$. In contrast, in the critical flow regime,  $C_D$ starts to decrease rapidly as  $Re$ increases until a minimum is reached at a critical Reynolds number. For a smooth sphere  $Re_{cr} \approx 3.7\times 10^5$. This critical flow regime some times is referred to as the drag crisis. Beyond the critical $Re$, in the so-called supercritical regime, the drag coefficient slowly increases with increasing $Re$ until it reaches a maximum value. Further increasing  $Re$, the drag coefficient stays constant and this regime is called transcritical. For  the prediction of $C_D$ at high $Re$ one usually resorts to numerical simulations [21, 22, 23, 24] or experiments [19, 25, 26].
The results of these numerical simulations and experiments are translated into  fitting correlations,  with a range of applicability limited to the range of the data that is used in the fitting process. This has resulted in a zoo of  correlations that take different mathematical  forms  [27, 28, 29, 30, 31, 32, 33, 34], as shown  in the  extensive list published in the recent review by Goossens <cit.>.  The majority of  correlations focus on the subcritical regime and  take the following functional form:

\begin{equation} \label{eq2}
C_D=\underbrace{ \underbrace{\frac{24}{Re}(C_1+C_2Re^a)}_{\text{Schiller and Naumman}}+\frac{C_3}{1+\frac{C_4}{Re}}}_{\text{Brown and Lawler}}
\end{equation}

The second term of  Eq.(<ref>) arises from  boundary layer theory [36], which accounts for the inertial effects of the fluid. The value of the  exponent $a$  ranges from  0.5 to 0.68. These type of correlations are suitable for $Re$  up to  $2\times 10^5$, right before  drag-crisis.

Concerning  the heat transfer rate from a particle fixed in a fluid, most  investigations available in the literature are related to the case of  forced convection. In this type of flow, the velocity profile is decoupled from that of the temperature.  For further simplification, there is also no variation in the transport properties of the fluid with temperature. These simplifications pave the way of obtaining several analytical solutions for a single sphere [37] for   limited cases of low $Re$ and Peclet number $Pe = v_{\infty} d / \alpha$, where $\alpha$ is the thermal diffusivity of the fluid. Acrivos and Taylor [37]  used  asymptotic expansions  and the  velocity profile of the Stokes solution to find  the following relation for the  Nusselt number $Nu = hd/k$, where $h$ is the (convective and surface mean) heat transfer coefficient and $k$ is the thermal conductivity of the fluid (linked to the thermal diffusivity through $k = \alpha \rho c_p$, with $c_p$ the specific heat capacity of the fluid),
for the case of $Pe\to$ 0 and $Re\to$ 0 :

\begin{equation}\label{eq3}
Nu = 2+\frac{1}{2}Pe+\frac{1}{4}Pe^2\log(Pe)+0.034 Pe^2 +\frac{1}{16}Pe^2\log(Pe)
\end{equation}

In practice, this solution is limited to $Re \lesssim$ 0.03. Rimmer [38] added an extra term to Eq.(<ref>) from asymptotic expansions, and as far as we know this is the last term that evolved from the matched asymptomatic expansions in the low $Pe$ and  $Re \to 0$ regime.  Conversely, for
$Pe\to \infty$ and  $Re \to 0$,  Acrivos and Goddard [39] used the matched asymptotic expansions  to  arrive at the following relation for  $Nu$:

\begin{equation}\label{eq4}
Nu =0.922+ 1.249Pe^{\frac{1}{3}}
\end{equation}

As for the case of drag, for higher $Re$ we need to rely on semi-empirical relations to express the variation of  $Nu$ with the flow field parameters. Whitaker[40] provided a correlation, which is still considered one of the most accurate available in  literature [41]:

\begin{equation} \label{eq5}
Nu =  2+ (C_4Re^{a_1}+C_5Re^{a_2})Pr^{a_3}
\end{equation}

where $Pr = c_p \mu / k$ is the Prandtl number (note that $Pe = Re Pr$). The values of $a_1$, $a_2$, and $a_3$ are $\frac{1}{2}$, $\frac{2}{3}$, and 0.4, respectively. The Whitaker correlation is valid for 1 $\leqslant Re \leqslant 10^5$ and a wide range of $Pr$. The second, and  third terms represent  inertial fluid effects, and their functional form is inspired by  boundary layer theory. Although the first term comes from the analytical solution for  pure conduction from a sphere, all exponents in Eq.(<ref>) are obtained from empirical fitting.

In summary, almost all correlations for drag and heat transfer found in literature are expressed as power law expansions, similar to Eqs. (2), (4) and (5). Correlations with logarithmic terms, such as Eqs. (1) and (3), are extremely rare and seem to have been largely overlooked.

The improvement of high-performance computer architectures, plus the availability of data from numerical simulations and experiments, sparked an increase  in interest to use machine learning methods to solve  problems in many scientific disciplines. This has led to label machine learning as the fourth paradigm in science, next to experimentation, theory and simulation [42]. When it comes to fluid mechanics, applying machine learning methods constitutes a challenge for several reasons, such as the transient nature of most fluid mechanics problems, the heterogeneity of most available data, the extensive non-linearities that govern fluid mechanics, and the multi-scale nature  of most problems in hand [43]. To deal with these challenges, an  ideal machine learning algorithm for fluid mechanics, should possess features such as interpretability, explainability, generalisability, and convergence [43]. One of the most popular machine learning frameworks that are used extensively in different fluid mechanics problems, from solving partial differential equations [44, 45], discovering physics [46], learning the active-nematic hydrodynamics [47], and predicting  physical properties [48] are the artificial neural networks (ANN). Other machine learning methods that are used for scientific discovery are  sparse identification of nonlinear dynamical systems for discovering differential equations from sparse data [49], and symbolic regression that is used for discovering laws of nature [50], discovering new materials [51], and solving
fluid flow problems [52].

In this paper we will use symbolic regression, which is a modern tool for unbiased determination of correlations, to re-investigate known data on drag and heat transfer. We will show that symbolic regression actually rediscovers the logarithmic terms, suggesting that logarithmic expansions may represent the physics better than power law expansions. As a side result, we will show that there is an intriguing connection between the found logarithmic terms and the point of first boundary layer separation.

§ METHODOLOGY

In this paper, we will use the symbolic regression machine learning method proposed by Koza [53]. Symbolic regression is a powerful tool for searching the mathematical space for an approximate functional relation between a certain number of input and output variables, and it is based on genetic programming proposed by
Holland [54]. The framework of genetic programming is probabilistic, and is not based on mathematical principles, such as correctness, consistency, justifiability, certainty, orderliness, and decisiveness as outlined by Koza [53], but solely on the principles of Darwinian evolution [55]. The idea of the genetic programming is simple, and it is based on transforming an initial population (in our case a population of  mathematical functions) to a new population that survived a particular fitness constraint. The main operators that are used to create the new population are similar to those found in nature, namely that of reproduction and crossover
[53] .

The algorithm first generates a random pool of functions, that undergo genetic operations such as crossover, which corresponds to the combination of two functions  to give a new offspring function. Another operation is a mutation in which a certain part of the mathematical function is changed randomly.  Two indexes measure the fitness of the newly obtained functions.    The first index  is minimizing the mean square difference between the training and predicted dependent values. The  second index is to check the mathematical complexity of functions, and select the simplest ones, to prevent over-fitting. We used the Eureqa software [50]  as symbolic regression    platform. A  rigorous description of the symbolic regression algorithm  in use in the current investigation  is given in [52].

In  Appendix A we illustrate that the machine learning algorithm we use can capture a known function's series expansion. It shows the ability of symbolic regression to find expansions of functions, that are valid beyond the training data used to obtain them, which gives symbolic regression an advantage compared to, for instance, artificial neural networks.

§ RESULTS

In the first subsection, we will explore the  dependence of the drag coefficient  $C_D$ on $Re$ for a fixed sphere. We will devote the second subsection to explore the  dependence of the Nusselt number   $Nu$ of a sphere on $Re$ and $Pe$ (or $Pr$) for the case of forced convection with constant transport properties.

§.§ Drag coefficient $C_D$

We will start by exploring the $C_D$ dependency on $Re$ for the case of a sphere. We will create three data sets for the regression process. The first one will be generated from the correlation of Brown and Lawler [56] which has the functional  shape of Eq.(<ref>). This data set contains about 8500 points in the range $ 0.1 < Re < 1.9\times 10^5$, which is enough to capture the smallest details in the $C_D$ variation. The second data set that we will use is the exact experimental data  that Brown and Lawler [56] used themselves  to derive their correlation. It contains about 450 points in the range $0.1 < Re < 1.975\times 10^5$ . The final data set is based on the Schilller and Naumann [57] correlation, and contains  of 5020 points in the range $0.1 < Re < 700$.

We will start by examining the first data set, and  we will  let the symbolic regression  algorithm guess about the functional form of the $C_D$ dependence on $Re$. We can do this by specifying the most general initial functional form:

\begin{equation} \label{eq6}
C_D =f(Re)
\end{equation}

The algorithm derived several regression equations, but here we will show  two, one because it accurately fits the results, and the other  because it is simple. The equations are the following:

\begin{equation} \label{eq7}
C_D = a_1 + \frac{a_2}{Re} + a_3\sqrt{Re} + \frac{a_4}{\sqrt{Re}}+ \frac{a_5}{(a_6+Re)} +a_7Re
\end{equation}

\begin{equation}\label{eq8}
C_D=a_1 + \frac{a_2}{Re} + \frac{a_3}{\sqrt{Re}}
\end{equation}

The coefficients of Eq.(<ref>), and Eq.(<ref>) are listed in Table <ref>. Eq.(<ref>) contains  the Stokes  $\frac{1}{Re}$ term, and the first-order term from  boundary layer theory $\frac{1}{\sqrt{Re}}$. The first known dependency of  $C_D$ on  $\frac{1}{\sqrt{Re}}$  came from the Blasius solution [58]  of the boundary layer equations proposed by Prandtl [59] for the case of a flat plate. The $C_D$ for blunt bodies, like a sphere, has a similar dependency on $Re$  [60, 61].  A similar form as Eq.(<ref>) was  obtained previously by fitting experimental data  [62, 63], and also by using concepts of boundary layer theory [61].
Refs [62, 63]  used
non-linear fitting tools to obtain their correlations,  which require  a priori knowledge of the functional structure. A comparison between the the coefficients of Eq.(<ref>), and those of  Refs     [62, 63, 61] is given in Table <ref>. The coefficients of Eq.(<ref>) have similar values to those of [62]. Compared to those of [63] there is only significant difference in the value of $a_3$. There is  also a significant difference between the coefficients of Eq.(<ref>) and those of Abraham [61]. This may be due to the pure theoretical nature of the equation proposed by Abraham.

It is important to note that both the Stokes term and the boundary layer term have been found without using any sophisticated mathematical approach. On the contrary, they have been found by a probabilistic genetic  algorithm. The emergence of the boundary layer term in Eqs. (<ref>) and (<ref>) without human intervention can be added to the experimental and numerical results that support boundary layer theory, even though there is no general mathematical proof of its existence, as mentioned by Batchelor [64].

We will now try to explore the existence of  logarithmic switchback terms for the drag on a sphere for the  higher $Re$ regime. We will use for this the first data-set  (i.e. data from the Brown and Lawler [56] correlation). We will start by imposing the following initial functional form:

\begin{equation}\label{eq9}
C_D = f(\frac{24}{Re}, \log(Re), Re\log(Re), \log^2(Re))
\end{equation}

We choose this form of the initial function because we want to ensure that logarithmic switchback terms similar to Eq.(<ref>) will be part of the initial soup of functions that the symbolic algorithm will further evolve. The symbolic regression algorithm converged to the
following equation :

\begin{equation}\label{eq10}
C_D = a_1 +{\frac{a_2}{Re}} +{a_3\log(Re)+ a_4\log^2(Re) + a_5\log^4(Re)}
\end{equation}

The values of the coefficients of Eq.(<ref>) are listed in Table <ref>. Eq.(<ref>) depends on  powers of $\log(Re)$ and also contains the  Stokes law term. The form of Eq.(<ref>) is partially fulfilling the  P&P conjecture [6] for $Re$ as high as $2\times 10^{5}$. Overall, Proudman and Pearson [6] made a profound statement more than 64 years ago, using only  mathematical intuition, and they may have been right when they  suspected that  logarithmic switchback terms are part of the solution.
It may be difficult for the current form of the genetic algorithm to spot the entire logarithmic switchback  series,  because  reducing the complexity of the equations is part of its optimization process. Therefore,   terms that do not play a  significant role in the variation of the dependent variable ($C_D$) will die out during the evolution process.  The failure of detection of  $Re^n\log^n(Re)$ terms,  where $n$ is an integer, after a  significant number of mathematical formula evaluations  exceeding $10^{11}$,  suggests that their signal is weak (a metaphor for their insignificant role in the dependence of $C_D$ on $Re$). If we read more carefully the conjecture, we  find that Proudman and Pearson [6] used the following wording: “ It seems reasonable to suppose that both expansions are in powers of Re".  They used the word `reasonable to suppose', expressing doubt, while for the   $\log(Re)$ terms they used the word `must' which reflects that the authors were sure about their appearance in the two expansions. Adding to that, Proudman [9] was frustrated about the poor convergence  of his equation, mainly because it is only valid for extremely low values of $Re$. He suggested that the expansion in  powers of $Re$ may be a poor idea [9, 65].

In order to further validate the ecosystem of the equations that we obtained, we will compare their predictions with various sources in the literature, as shown in Figure <ref>. The first insight from  Figure <ref> is that  Eq.(<ref>) is valid only at low $Re$, and this was one of the main reasons we believe that the scientific community did not  further explore the use of  logarithmic terms, even as fitting functions. Eq.(<ref>) and Eq.(<ref>) follow closely the correlation of Brown and Lawler [56], and also the experimental data used to obtain  their correlation. The average relative errors  between the predictions of Eq.(<ref>) and Eq.(<ref>)  with respect to the experimental results of [56] are 3.87% and 3.39%, respectively. We see that Eq.(<ref>) follows closely the results of [63, 62], while it deviates from the predictions  of Abraham [61] especially for values of $Re$  above $10^3$. This is expected because the equation provided by Abraham [61] is valid for $Re$ up to $10^3$. Also, Eq.(<ref>) and those of
references [62, 63, 61] cannot capture the local minimum for Re between $10^3$ and $10^4$ that the experimental results of [56] show.

Comparing   Eq.(<ref>) and Eq.(<ref>), we find that  their complexity  index is  34 and 19, respectively. The complexity index shows that the logarithmic series  representation of  $C_D$ is mathematically simpler compared to the power series representation, making Eq.(<ref>) more favourite to represent the physical phenomena of the $C_D$ variation  according to  Occam's razor statements [66]. One of these statements is: “Given two models with the same generalization error, the simpler one should
be preferred because simplicity is desirable in itself."

Now we will use the second (experimental) data set, to explore the feasibility of getting predictive equations for $C_D$ from a limited amount of noisy experimental data. We will start by letting the algorithm guess the $C_D$ dependence:

\begin{equation} \label{eq11}
C_D =f(Re)
\end{equation}

The symbolic regression algorithm found the following equation:

\begin{equation} \label{eq12}
C_D = a_1 + \frac{a_2}{Re} + \frac{a_3}{\sqrt{Re}}
\end{equation}

The coefficients of Eq.(<ref>) are listed in Table <ref>. Using the second data set we next  explore if the data show any  logarithmic dependence  by  imposing the following initial set of functions:

\begin{equation} \label{eq13}
C_D = f(\frac{24}{Re}, \log(Re), Re\log(Re), \log^2(Re))
\end{equation}

We got the following equation for $C_D$:

\begin{equation} \label{eq14}
C_D =  a_1 +\frac{a_2}{Re} +\frac{ a_3\log^2(Re)}{Re}+   a_4\log(Re)+ a_5\log^2(Re)
\end{equation}

The values of the coefficients are listed in Table <ref>. Eq.(<ref>) is of a similar form as Eq.(<ref>), but the coefficients are not identical, because the second data set contains far less data, and also  contains some noise.  The derivation of Eq.(<ref>) from pure experimental data, without imposing  knowledge of any physics, except the definition of $Re$, shows that the symbolic regression algorithm discovered the Stokes limit and the term attributed to  boundary layer theory  without any external help. The algorithm needed less than an hour to discover what took human intellect  hundreds of years to achieve. However, the human factor is still required  since we have to select the equations that we think  represent physical reality from the population of equations  that the algorithm suggests. Eq.(<ref>) shows that we can get the logarithmic dependence from a pure experimental data set, and  it partially fulfils the P&P conjecture. Eq.(<ref>) and Eq.(<ref>) are quite similar. We believe that Eq.(<ref>) failed to capture the $\log^4(Re)$ term because this term influences  $C_D$  in the high $Re$ regime where  there are significant fluctuations in the experimental  data set.  Probably if there were a higher volume of data, especially at higher $Re$, the $\log^4(Re)$ term could also be captured from pure experimental results. A comparison of the  performance of the power expansion  Eq.(<ref>) and the logarithmic expansion Eq.(<ref>) against existing data in the literature is shown in Figure. <ref>. The average relative error for Eq.(<ref>) and Eq.(<ref>) is 13.7% and 12.0%, respectively, against the experimental results of [56].  Eq.(<ref>) shows a local minimum in the range of the $Re$ close to that of the experimental results of  [56], while Eq.(<ref>) fails to show any local minimum.

We will use the third and final data set from the Schiller and Naumann [57] correlation  which contains information about the variation of $C_D$  for $Re$ ranging from 0.1 to 700.  We will use the following general initial functional form:

\begin{equation}\label{15extra}
C_D = f(Re)
\end{equation}

The symbolic regression algorithm found the following equation for  $C_D$:

\begin{equation}\label{16extra}
C_D =  a_1 + \frac{a_2}{Re} +a_3\log(Re)+ a_4\log^2(Re)
\end{equation}

The   coefficients of Eq.(<ref>) are listed in Table <ref>. The genetic algorithm came up with the logarithmic dependence of $C_D$ on $Re$  without any external help, and it discovered the P& P conjecture partially. The value of $a_1$ = 3.1406, differs from the value of $\pi$  by only about 0.03%. It will be very interesting in the future to investigate the value of $a_1$ by fitting to very accurate numerical or experimental data.  Eq.(<ref>) follows the Brown and Lawler correlation [56] up to  $Re$ of $10^3$, as shown in Figure <ref>. This behaviour is expected because  higher power  logarithmic terms  are missing from Eq.(<ref>), since the training data was limited to $Re$ up to 700.

Up to this point we have discussed the drag without referring to the flow around the sphere. The flow around a sphere is a rich mosaic of phenomena, and usually drag correlations, fail to predict them. Among these  phenomena is the  emergence of a laminar separation point, which  is well known to occur for sufficiently blunt  objects, including a sphere. The point of  laminar separation is identified by the formation of a closed recirculating ring eddy at the rear of the sphere, as indicated in Figure <ref>. The first emergence of separation  is difficult to detect either experimentally or theoretically. For this reason, there is some discrepancy  in the literature on the value of the reported critical $Re_s$, and  corresponding drag $C_{Ds}$, at first separation. The first experimental observations by Nisi and Porter [67] suggested  that $Re_s$ = 10. This was  confirmed by numerical simulations of Rimon and Cheng [68].  On the other hand, Proudman and Pearson [6], and  Van Dyke [12], by using the Stokes second  expansion,   estimated that $Re_s$ = 16, close to the numerical results of Bourot [69] and Jenson [22] of  15.2 and 17, respectively, and the experiments of Payard and Countanceau [70] indicating  $Re_s$ = 17. Other  simulation results    [21, 71]     show that $Re_s$ is equal to approximately 20, and the experiments of Taneda [72] predict that $Re_s$ = 24.

If we inspect $a_1$ of the logarithmic expansion Eq.(<ref>) in Table <ref> we  see that its value is  3.286, which is quite similar to the value of the drag coefficient $C_{Ds}$ at the initial laminar separation reported by [70], which  is 3.306. If the constant $a_1$  is the drag coefficient at initial laminar separation, then the following transcendental equation must have a positive root at the corresponding Reynolds number $Re_s$:

\begin{equation}\label{eq15}
\frac{a_2}{Re} + a_3\log(Re)+ a_4\log^2(Re) + a_5\log^4(Re)=0
\end{equation}

By solving Eq.(<ref>) we find that  $Re_{rt}$ = 14.06 is its only root. That makes  $Re_{rt}$  the only $Re$ value that zeroes off all terms beyond the constant $a_1$. This  $Re_{rt}$  is close to values of  $Re_s$ reported in literature. For example, the relative error with respect to the results of Bourot[69] and Chang and Maxey [71] is 8% and 30%, respectively. We conjecture that $Re_{rt}$ is representing $Re_s$, even though   we do not have any  proof for this. We believe we are witnessing an instance where the machine learning algorithm found a mathematical description of a physical phenomenon, which needs human  abilities to be interpreted  in terms  of  physical laws. Otherwise, it will be a good approximation, that can describe some of the physics involved in the process of flow separation. As far as the authors are aware, there is  only one analytical prediction for the  point of first flow separation, from  slow motion viscous theory [6, 73]. However, that result was disputed  by the authors of   [6, 73], as we will show later.   In practice, we depend on  numerical simulations to find the point of zero local shear stress, as described by  boundary layer theory [36]. However, Batchelor [64] raised serious doubts about estimating the onset of  separation by this method.

Beyond this point, we will assume that  (the smallest, real) root $Re_{rt}$ is equal to $Re_s$. Using the same procedure to calculate  $Re_s$,  from the logarithmic  Eq.(<ref>) by solving the following  transcendental equation:

\begin{equation}\label{eq16}
\frac{a_2}{Re} +\frac{ a_3\log^2(Re)}{Re}+   a_4\log(Re)+ a_5\log^2(Re)=0
\end{equation}

we found the two following roots: $Re_s$ =  15.76, and $9.52\times 10^7$. The large root value of $ 9.52\times 10^7$,  is a non-physical result, which we believe is caused by the missing higher power $\log(Re)$ term from Eq.(<ref>). However, $Re_s$ = 15.76 compares very well with the results  of Bourot[69] and Chang and Maxey [71], with a relative difference of 3.68% and 21.2%, respectively. If we do the same analysis for the logarithmic  Eq.(<ref>), we will find that $Re_s$ =  15.19, and $3.518\times 10^6$. For the smallest root,   the relative difference with the results of Bourot[69] and  Chang and Maxey [71] is 0.13%, and 24.0%, respectively.

We will  next calculate $Re_{rt}$ from the more popular power-law expressions  Eq.(<ref>) and Eq.(<ref>) in  the same way. For Eq.(<ref>) we find the following roots $Re_{rt1} = -2461-767i$, $ Re_{rt2} =-2461+767i$, and $Re_{rt3} = 3\times 10^5$. The first two roots are non-physical, while the third root, closely approximates the critical Reynolds number ($Re_{cr}\thickapprox  3.7\times 10^5$)  for the critical flow regime (drag crisis) as reported by Achenbach [19]. We will further discuss the physical significance of   $Re_{rt3}$ in the generalization subsection since the value of$Re_{rt3}$  is outside the training data range.  As for power-law Eq.(<ref>), it does not have any roots, neither in the real nor in the complex domain.

Returning to the logarithmic ecosystem of equations,  in their seminal works, Proudman and
Pearson [6] and  Van Dyke [12]  calculated the $Re_s$ value to be  16 analytically from the first and second terms in the Stokes expansion. Proudman
and Pearson[6] made the following comment:  “This Reynolds number is far too large to make estimates based on only two terms of the Stokes expansion at all reliable. In fact, it cannot seriously be claimed that slow-motion theory gives even a qualitative expansion of the phenomena.” However, Van Dyke [12]  and Ranger[74] tried to confirm the result of Proudman and Pearson[6], by using  extra terms in the Stokes expansion that contain the logarithmic terms from the results of Proudman and Pearson [6] and those of Chester et al. [9]. They failed because the Stokes expansion equation that includes the logarithmic terms has only complex roots. Van Dyke [12] commented on this issue saying that “ the logarithm needs reinterpretation." In our work we now see that  the values of $Re_s$ from Eq.(<ref>), Eq.(<ref>), and Eq.(<ref>) are converging with different degree of accuracy toward a value of approximately 16.

§.§ Generalization beyond the training data

In this subsection, we will  test our newly derived equations generalisation behaviour, for flow regimes that were not  included in the training data. Specifically, we will test their behaviour for the low Reynolds number regime for $Re$ down to $10^{-4}$, and for the critical flow regime for $Re$ up to $10^6$.

§.§.§ Low $Re$ flow regime

In the low $Re$ regime, $\dfrac{24}{Re}$ is the dominant term for the drag coefficient, which will make it difficult to assess the performance of our equations, against the existing correlations, analytical solutions, experimental and numerical results. For this reason, we will use the way Maxworthy [75] plotted his drag coefficient data. He plotted the quantity $\dfrac{C_D}{C_{Ds}}-1$ against  $Re$,where $C_{Ds}$ is  the Stokes drag ($\dfrac{24}{Re}$). This way, we eliminate the divergence of the Stokes term, which makes the comparison with different sources from the literature more precise. From  low Reynolds number theory we know that $\dfrac{C_D}{C_{Ds}}-1$ converges  to $\dfrac{3}{16}Re$  (Oseen term) for extremely  low $Re$.

The predictions for the variation of $\dfrac{C_D}{C_{Ds}}-1$  against $Re$ from our models and numerous sources from  literature are shown in Figure <ref>. In the range of $Re$ $10^{-1}$ to $10$, which is within the range of the  training data, all our derived equations, plus the Brown and Lawler [56] correlation, follow with reasonable accuracy the experimental results of Maxworthy [75], and Veysey and Goldenfeld [76], in addition to the numerical results of Jenson [22], and Dennis and Walker [21]. In the same $Re$ range, the analytically derived equations of Proudman and Pearson [6], Goldstein [77], and Oseen [2] deviate from experimental, and numerical results, because of their limited applicability range.

Next we turn to the $Re$ range between $10^{-4}$ to $10^{-1}$, which is beyond the training data range. In this flow regime, the logarithm-based equations(<ref>) and
(<ref>) follow closely the analytical results of <cit.>, and the semi-empirical and empirical correlations of [78, 79], and the numerical simulations of [80].  On the contrary, the power-based equations Eq.(<ref>), and Eq.(<ref>), as well the Brown Lawler [56] correlation, divert significantly from the analytical, experimental, and numerical data. For example  the relative difference for the prediction of $\dfrac{C_D}{C_{Ds}}-1$ between Eq.(<ref>) and the analytical solution of Proudman and Pearson [6] is 240% at $Re$ = $10^{-4}$.  At the same conditions, the relative difference between the Brown and Lawler [56] correlation and Proudman and Pearson [6] is 1410%, which is significantly higher than the error generated by both logarithmic equations. The five times increase in the accuracy of the logarithmic based equations Eqs(<ref>,<ref>) compared to the power-based equations Eqs(<ref>,<ref>)  suggests that  the  logarithmic equations contain terms that describe the physical reality better. Another interesting aspect of the results of Figure <ref> is that it shows that we can improve the accuracy of  machine learning models for the same training data set, by using previous physical knowledge about the problem at hand.
The observation from Figure <ref> is similar to our observations for the Maclaurin expansion of the $\sin(x)$ function in the Appendix A. In both cases,   only  equations that have similar terms to the actual representation of a  function, or the  physical law  that they are approximating, generalize well beyond their training data. The results from Figure <ref>, show that the popular power-based representation of $C_D$ fails to extrapolate beyond the range of $Re$ that is used for its training, which indicates that the power-based  representation may have only  been a convenient  mathematical fit, rather than having  physical significance. Finally, we want to explain why Eq.(<ref>) diverges even though it consists of logarithmic terms similar to the previous two. The reason for the divergence is the  $\dfrac{a_3\log^2(Re)}{Re}$ term which increases its value as the value of $Re$ decreases. This term can be considered an overfitting parameter, which it is easy to spot , due to the interpretable nature of the results of  symbolic regression.

§.§.§ Critical flow regime

The critical flow regime is less well investigated, neither experimentally or numerically, compared to the subcritical or lower $Re$ regimes. There are not  any analytical approximations for  $C_D$ in the critical flow regime. Even  direct numerical simulations (DNS) are limited to the onset of the subcritical flow regime at $Re$ = $10^4$ [81]. Current computational fluid dynamics (CFD) simulations  that deal with the critical flow regime use different approximations to deal  with turbulence. Constantinescu et al. [24]  use  Detached-Eddy-Simulations (DES), which is a hybrid method that combines Reynolds-Averaged Navier-Stokes (RANS), and Large Eddy Simulations (LES).  Nakhostin and Gilijahus [23]  used RANS turbulence models for their simulations, and  Muto et al.[82]  used  Large Eddy Simulations coupled with the a subgrid-scale turbulence model. The most extensive numerical simulations in the critical and supercritical regime have been  conducted by Geier et al.
[83]  using a Cumulant Lattice Boltzmann method, and they do not use any turbulence models. Their high fidelity model uses a fourth-order accurate diffusion approach, suitable for low viscosity  of high $Re$ flows. The accuracy of the Cumulant Lattice Boltzmann depends on the optimization of  its parameters. The authors used a spectrum of three different mesh grid schemes, namely a course one with $40\times 10^6$ nodes, a medium one with $75\times 10^6$ nodes, and a fine grid mesh with $133\times 10^6$ nodes.

Figure <ref> explores the performance of the power-based Eq.(<ref>) and logarithm-based Eq.(<ref>) in  the subcritical,   critical,  and supercritical flow regimes, and compares their performance against experimental and numerical results. The training data for Eq.(<ref>) and Eq.(<ref>) was limited to $Re$ up to $2\times10^5$. There is a significant discrepancy between the different experimental results,  for different reasons, such as the turbulence intensity the positions of the sensors around the sphere [64].  Eq.(<ref>) follows the anticipated trend in the critical flow regime in which the $C_D$ is decreasing with increasing  $Re$. Note that on the contrary, the value of     $C_D$ from the correlation of Brown and
Lawler [56]  stays constant for $Re$ values higher than $10^4$. The onset of the critical flow regime for the power-based equation Eq.(<ref>) starts at approximately   $Re\approx 10^5$, earlier than most  experimental and numerical results, except the experimental data of Maxworthy  [26], in which the critical flow regime starts at much lower $ Re$. At approximately  $Re = 3\times 10^5$ Eq.(<ref>) drops to zero, and its values resembles  the experimental values of Achenbach [19]. The drop of Eq.(<ref>) to zero at $Re$ = $3\times 10^5$ was already predicted algebraically in the previous section, and Eq.(<ref>) is the first in literature that predicts  with good accuracy  the value of $Re_{cr}$  reported by the  experiments of  Achenbach [19]. From Figure <ref>,  we can see that even the high fidelity simulations of Gerier [83] with fine grid failed to predict $Re_{cr}$ since they failed to resolve the Kolmogorov length scale at such high $Re$. The numerical results for the medium grid scheme of Gerier et al. [83] are close to the predictions of the Eq.(<ref>) for $Re$ until the critical Reynolds number. The logarithmic based equation (<ref>)  predicts the onset of the critical flow regime  with great accuracy since it follows the $C_D$ values from the experiments of Suryanarayana et al.
[84], and Achenbach [19] from $Re$ = $5\times10^4$ to about $3\times 10^5$. Eq.(<ref>) does not does not drop to zero at $Re_{cr}$ as Eq.(<ref>), however it follows very closely the high fidelity numerical results of Gerier et al. [83], for the coarse grid case for  $Re$ up to $10^6$. This shows that Eq.(<ref>) follows an approximately physical reality for $Re$ up to  $10^6$, since the results of Gerier et al.
[83] are generated by solving an approximate form of the Navier-Stokes equations.  Both Eq.(<ref>) and Eq.(<ref>) fail to predict the increase of  $C_D$  after the end of the critical flow regime, and the start of the supercritical flow regime at which the boundary layer attached at the surface of the sphere changes from being partly laminar  to being fully turbulent. This failure is attributed to the fact the training data used to obtain Eq.(<ref>), and Eq.(<ref>) are far from the critical flow regime. Predicting  $C_D$ for the critical flow regime is  difficult even for high fidelity solvers. For example, the non-optimized (Nonopt) solver of
Gerier et al. [83] failed to predict the drag crisis. Instead, it predicts that $C_D$ does not change with $Re$, similar to what the correlation of Brown and Lawler [56] predicts. Eq.(<ref>), and Eq.(<ref>) perform better in the critical regime than the fitting correlation of Morrison [85] which is a result of fitting experimental data from the literature.  Another interesting observation is that the rate of change of  $C_D$ with  $Re$ in the critical flow regime, for both Eq.(<ref>) and Eq.(<ref>), follows the smooth trend similar to the experiments of Maxworthy [26] and the high fidelity simulations of
Geier et al. [83], rather than the sharp nearly discontinuous  change of  $C_D$ observed in the experiments of [19, 84, 86].

Both Eq.(<ref>), and Eq.(<ref>) predict different stages of the critical flow regime with surprising accuracy. They are the first in  literature to make such predictions without being exposed to the critical flow regime, but only by using a limited amount of physics stored in the training data and the imposed functional forms. The question  may arise whether  these predictions are just a product of chance? Our short answer is no, for several reasons. The first reason is that the $Re$ number changes by orders of magnitude in the critical flow regime, which gives many  possibilities for the output of the predictive function, but Eq.(<ref>)  predicts with nearly zero error the experimental results of
Suryanarayana et al. [84] concerning  the onset of the critical flow regime. The same applies to the $Re_{cr}$ predicted by Eq.(<ref>) compared to the experimental results of Achenbach [19]. The second and  more supportive reason is that  symbolic regression  can generalize and predict the approximated function's unexpected behaviour, similar to the example  shown in  Appendix A about the $\sin(x)$ approximation. The algorithm was trained to predict the peaks; however, it also accurately predicts the existence of valleys.  We strongly believe that  Eq.(<ref>) contains terms that approximate the fundamental physical law that  $C_D$ is following, which is why it managed to generalize both  the Stokes and  critical flow regimes. This makes the logarithmic representation of  $C_D$ a serious candidate an the analytical  mathematical formulation that governs the  variation of $C_D$  with the $Re$.

In summary, in this section we showed that the functional form of $C_D$ could be represented  by both powers and logarithmic functions of $Re$. However, the logarithmic representation conveys  the physics in a different way than the power representation, and illuminates new physical phenomena, which are beyond the reach of  current analytical, or empirical  $C_D$ formulas. Because of  the logarithmic equations' good generalization behaviour, especially Eq.(<ref>), such equations should not be considered as merely  fitting equations, but rather as semi-analytical equations. When appealing to mathematical aesthetics, our results suggest that the drag coefficient of a sphere might be well described by the form $C_D = \pi + 24/Re + f(\log Re)$, with $C_D = \pi$ at the first point of separation, occurring at a Reynolds number $Re_s$ given by the transcendental equation  $24/Re_s + f(\log Re_s) = 0$. Van Dyke [12] described the appearance of logarithms in the asymptotic expansions as obscure, but it appears that these obscure entities can speak the language of fluid dynamics much better than  powers. A similar situation exists in the field of turbulence, especially regarding  channel flow, where there is an open debate in the scientific community whether   power or logarithmic expansions  bests describe the velocity at the wall in certain flow regimes [87]. Note that the logarithmic dependence of the drag coefficient $C_D$ also exists for geometries different than a sphere such as spherocylinders, and prolate spheroids, as shown by our previous work [52].

§.§ Nusselt number $Nu$

In this section, we will explore the possibility of a logarithmic dependence of the Nusselt number  $Nu$ on the Peclet number  $Pe$ and  Reynolds number $Re$. For this purpose we will create a data set of 26,796 points from the Whitaker [40] correlation Eq.(<ref>) for   $Pr$ in the ranging   from  0.74  to 7.0, and  $Re$ in the range of 10$^{-1}$ to $10^4$. We will start with the simplest assumption by allowing the  symbolic regression algorithm to guess about the dependency of $Nu$ on $Re$, $Pr$ and/or $Pe$, through the following initial function:

\begin{equation} \label{eq20}
Nu = f(Re, Pr, Pe)
\end{equation}

The resulting $Nu$ correlation is the following:

\begin{equation} \label{eq21}
Nu = a_1 + a_2\sqrt{Pe} + a_3\sqrt{Re}\sqrt{a_4 + a_5\sqrt{Pe}} + a_6Pe + a_7Re
\end{equation}

The coefficients are listed in Table <ref>. Most equations that the algorithm produces show that  $Nu$ is a function of $Re$ and $Pe$, and excludes the explicit dependence  on $Pr$. This is  different from the source of our data (the Whitaker correlation Eq.(<ref>)), which explicitly depends on $Pr$ and $Re$. Even when  we used a  substantial amount of data, the algorithm failed to predict the exact structure of the Whitaker correlation [40].  The recent investigation of Udrescu and Tegmark [88] showed, consistent with our results, that Eureqa failed to predict the exact functional structure of many  functions included in the Feynman lectures [89]. They attributed this failure due to the complexity of those functions, and the number of variables that they contain.

Examining the properties of Eq.(<ref>), we find  that as $Re \to 0$, Eq.(<ref>) reduces to $a_1 + a_2\sqrt{Pe}$,  which bears similarities with Eq.(<ref>) for the $Pe$ dependency, because for both cases the power of $Pe$ is less than one, and both equations show that even at very low $Re$  convection affects the heat transfer rate. This type of dependency did not exist in the Whitaker correlation Eq.(<ref>), where for  $Re \to 0$ (outside the range of validity of the Whitaker correlation)  $Nu$   converges to a value of 2.0,   corresponding to pure conduction from a single sphere.

We will now  examine the full dependence of  $Nu$ on  logarithms of $Pe$, $Re$, and $Pr$. This structure of dependency is based on our previous knowledge of the physics of the problem of forced convection over a sphere. We know that for $Re \to 0$  and $Pe <1$,  $Nu$ depends on  $\log(Pe)$ [37] (Eq.<ref>), so there may exist an intermediate $Pe$ regime where logarithms  will play a role as well, until we reach a high $Pe$ regime where Eq.(<ref>) is dominant. For the high $Re$ regime we already showed that the drag  coefficient $C_D $ is a function of  logarithms of $Re$, so because of the tight relation between flow and heat transfer [90] we expect that logarithms of $Re$ will play a role in the convective heat transfer process as well. The initial  function has  the following form:

\begin{multline} \label{eq24}
Nu = f(\log(Pe), Pe \log(Pe), \log^2(Pe), \log(Re), Re\log(Re),\\
\log^2(Re), \log(Pr), Pr\log(Pr), \log^2(Pr))
\end{multline}

As initial guess  we gave equal weight to all  functional forms, to avoid any bias, toward any of the independent variables. The symbolic regression algorithm  found the following two correlations:

\begin{equation} \label{eq25}
Nu = a_1 + a_2\log^2(Re)\log(Pe)Pe^{a_3} + a_4Pe^{a_5}
\end{equation}

\begin{equation} \label{eq26}
Nu = a_1+ a_2\log^2(Re) + a_3Pe^{a_4} + a_5\log^2(Re)\log(Pe)Pe^{a_6} +a_7\log(Pe)
\end{equation}

The second equation is more complex than the first. The coefficients of both Eq.(<ref>), and Eq.(<ref>) are  listed in Table <ref>. Both equations posess  very interesting features. We will  start with Eq.(<ref>), where  the term $a_1+a_3Pe^{a_4}$ resembles closely the  approximation of Eq.(<ref>). The relative difference of  the $a_1$, $a_3$ coefficients and those of Eq.(<ref>) is 15%, and 8%, respectively. The relative error is remarkably small, if we take into  account that the source of the data set is coming from an empirical correlation that has an average predictive error of 30%.

We believe that the combination of the logarithmic dependence of $Pe$ and $Re$ plays an essential role in the emergence of an asymptotic solution. It seems there  are very few possible ways to fit the data of [40] using logarithms of $Pe$ and $Re$ and one of those few is using terms similar to   Eq.(<ref>). Our findings show the essential role played by  previous physical knowledge of the problem in specific regimes, to help the machine learning algorithm to reach a physically meaningful result.

The genetic algorithm predicted the asymptotic solution for the high $Pe$ (Eq.<ref>) case, rather than for low $Pe$ (Eq.<ref>), probably  because our training data  is more biased toward the high $Pe$ regime. Since the lowest $Re$ and $Pr$  used are 0.1 and 0.7 respectively,  the lowest  $Pe$ we used is 0.07, which lies at the boundary of the high $Pe$ regime. We could not use lower $Pe$ because  the Whitaker correlation [40] is based on $Re$ ranging between 3.5
and  $7.6\times10^4$, and $Pr$ ranging between 0.7 and 380. Note that we did use  the Whitaker correlation [40] also for lower Re, $0.1 < Re < 3.5$, to generate our training data. We test its validity against the experimental data of
Will et al.  [91] for the lowest Prandtl number  that we used, $Pr$ = 0.7, and for  $Re$ as low as 0.1, and we found that the Whitaker correlation [40] follows closely the results of [91], as shown in Figure <ref>. An indication that the hydrodynamics in the highly inertial regime may be governed by logarithmic terms of $Re$, is the the appearance of $\log^2(Re)$ terms both in Eq.(<ref>) and Eq.(<ref>), similar to the case of $C_D$ (see Eqs.(<ref>), (<ref>) and (<ref>)). Also, the $\log^2(Re)$ terms for both  $Nu$ and $C_D$ share the same sign, and their pre-factors are of the same order of magnitude.

We  compare the performance of  our predictor equations for different $Pr$, and $Re$ numbers,  in Figure <ref>.   We select four cases,  two of them lie within the training data set ($Pr$ = 0.7 and 7.0) that we supplied to the algorithm. The other two test cases ($Pr$ =50 and 300)  lie outside the training data set to test the extrapolation capabilities of our predictor equations. For $Pr$ = 0.7, Eqs (<ref>), (<ref>) and (<ref>) perfectly follow the Whitaker[40] correlation and the experimental results of Will et al.[91]. At high $Re$ they  also follow  the numerical results of Feng and Michaelides [92]. As expected, our ecosystem of equations do not follow the asymptotic solution of Acrivos and Goddard [39] since their solution is only valid in the  low $Re$ and high $Pe$ regime. For the case of $Pr$= 7.0, our ecosystem of equations predicts the evolution of $Nu$ with great accuracy. For the cases of $Pr$ = 50, and 300, Eqs.(<ref>) and (<ref>) predict with  great accuracy the results of the Whitaker[40] correlation, except in a very narrow region at  low $Re$. The conditions in this low $Re$ - high $Pr$ regime are applicable to  the asymptotic   solution of Acrivos and Goddard [39]. This is why  the whole ecosystem of our equations deviate from the results of the Whitaker[40] correlation, and follow by different degrees of accuracy  the asymptotic solution of Acrivos and Goddard [39], Eq.(<ref>). All of our equations are functions of $Pe$ and $Re$. However, for  low $Re$ the $Nu$ correlations switch to a   dependency on $Pe$ only, which is consistent with the physics of Eq.(<ref>) and (<ref>).

The above  shows that  symbolic regression can find an asymptotic solution by using previous physical knowledge, rather than depending completely on the training data set.  Feeding machine learning algorithms  previous physical knowledge for the problem that they try to optimize, increases substantially  the probability of better extrapolation predictions. For further discussion on how to implement previous knowledge into symbolic regression, the readers is referred to our recent publication [52].

§ CONCLUSIONS

In this investigation, we  explored the possibility of a logarithmic dependence of the drag coefficient $C_D$ on the Reynolds number $Re$, and the  Nusselt number $Nu$ on  $Re$ and  Peclet number $Pe$, inspired by  asymptotic solutions for  creeping flow conditions. We used a symbolic regression machine learning algorithm, and our training data are based on experiments, and data from well-known empirical correlations available in the literature.   We can make  the following conclusions :

* The drag coefficient $C_D$ can be expressed as a function of powers in  $\log(Re)$, partially fulfilling the Proudman and Pearson [6] conjecture P&P.

* If an expansion in terms of $\log(Re)$ is made for the drag coefficient $C_D$, the value of the $Re$ at which all the $Re$ dependent terms go to zero is closely resembling the $Re$ at the first emergence  of  laminar separation, as  predicted analytically by Proudman and Pearson

* The logarithmic dependence of $C_D$ on $Re$ is found independently, without any prior knowledge, by the symbolic regression algorithm.

* The logarithmic based Eq.(<ref>) can generalize in both low, and high $Re$ regimes. In the high $Re$ regime Eq.(<ref>) can predict the drag crisis, its results  closely following  experimental, and numerical predictions from  literature.

* Since Eq.(<ref>)is derived from the experimental data of Brown and Lawler [56], the appearance of the logarithmic terms in  $C_D$ equations is independent of the  correlation that is  used as a source of the training data.

* The Nusselt number of a single sphere depends on  logarithms of $Re$,  $Pe$, as well as  powers of $Pe$.

* If logarithmic functions of  $Re$ and $Pe$ are used as initial functions for the symbolic regression algorithm, the  algorithm produces with high accuracy the asymptotic solution derived by Acrivos and Goddard [39] from the matched asymptotic method, in the low $Re$ and high $Pe$ regime. Interestingly, the training data that we used  does not follow the asymptotic solution of  Acrivos and Goddard [39].

* There is a connection between the appearance of the logarithmic terms in both $C_D$, and $Nu$ expressions, and the ability of those expressions to generalize outside the training data range. This connection makes the logarithmic representation a strong candidate for the functional form of $C_D$ and $Nu$ that could result from solving the Navier-Stokes equations analytically for the problem of flow over a single sphere at high $Re$, and be a result of a generalized fluid mechanics theory that applies to both low and high $Re$ regimes.

The bigger picture of our results is that, although our method cannot give answers as rigid mathematical proofs, it is highly probably that if one day we manage to solve in a closed form the Navier-Stokes equations, combined with a heat equation around a sphere, this solution will be expressed in terms of logarithms rather than powers. The logarithmic terms that symbolic regression found are related to the velocity and pressure fields around the sphere. Symbolic regression is an excellent candidate to further investigate the functional form of these fields, and we intend to conduct a future study toward this goal. Finally, we note that the machine learning framework that we developed is general and can be used in different scientific disciplines with the condition that experimental and numerical data exists, plus the availability of some limited analytical solutions.

§ ACKNOWLEDGEMENTS

The first author  thanks Dimitra Damianidou for the enlightening discussions about the subject. The authors thank Lorenzo Botto for the discussion about  matched asymptotic methods. Finally, the authors thank the European Research Council for its financial support under its consolidator grant scheme, contract no. 615096 (NonSphereFlow).

§ APPENDIX A

§ MACLAURIN EXPANSION OF SIN FUNCTION

A well-known result of applied mathematics is the representation of  continuous functions by the  Taylor expansion[93]:

\begin{equation}\label{eq.1s}
f(x)= \sum_{n=0}^\infty \frac{f^n(a)(x-a)^n}{n!}
\end{equation}

When $a$ = 0, the Taylor series reduces to the Maclaurin series. The following expansion gives the Maclaurin series for $\sin(x)$:

\begin{equation}\label{eq.2s}
\sin(x) = \sum_{n=0}^\infty\frac{(-1)^n}{(2n+1)!}x^{2n+1} = x-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}+...
\end{equation}

One of the reasons we choose the $\sin(x)$ function as our test case for the symbolic regression algorithm is its non-monotonic nature, specifically its  transition from an increasing to a decreasing function. This feature will help us assess the generalization behaviour of the algorithm. We generated 5000 uniform training points in the range  [0,$\dfrac{\pi}{2}$]. We selected this specific range because we wanted to feed the algorithm  only the monotonically increasing part of the $\sin(x)$ function, and see if it can generalize, and predict the decreasing part of the function between [$\dfrac{\pi}{2},\pi$]. The algorithm does not possess any prior knowledge of the $\sin(x)$ function and starts by assuming the most primitive initial function for the symbolic regression algorithm:

\begin{equation}\label{eq.3s}
y = f(x)
\end{equation}

The symbolic regression algorithm suggested many equations, including the following two:

\begin{equation}\label{eq.4s}
y(x) =  a_1x-a_2x^3+a_3x^5-a_4x^7
\end{equation}

\begin{equation}\label{eq.5s}
y(x) = a_1x-a_2x^3+a_4x^4+a_5x^5
\end{equation}

The values of the coefficients of Eq.(<ref>), and Eq.(<ref>) are listed in Table <ref>. Eq.(<ref>) contains the first four terms of the Maclaurin series for the $\sin(x)$ function. Although this may seem to be trivial, to the best of our knowledge this is the first time that a machine-learning algorithm managed to derive a Taylor or a Maclaurin series out of pure data. For the derivation of any Taylor series of a function we need to use the  calculus  invented simultaneously by Newton [94] and  Leibniz [95].

First, we want to illustrate the effect of the different terms of Eq.(<ref>) on its accuracy and generalization, as shown in Figure <ref>. For the [$0,\dfrac{\pi}{2}$] domain, except for the first linear term, regardless of the number of terms we add, the  decreasing nature of  $\sin(x)$  for $x>\dfrac{\pi}{2}$ is predicted.  Adding more terms increases the accuracy. While the first three terms are enough to predict with great accuracy the training data, the fourth term plays a significant role for values of $x > \dfrac{\pi}{2}$ which is beyond the range of the training data. We chose Eq.(<ref>) not only because of its accuracy but due to its resemblance of the Maclaurin series, thus our selection is based on our own previous knowledge. What is missing is a generalization theorem which can tell us about the generalization behaviour of a specific machine learning algorithm, trained at a specific range of data.  Without this theorem, we will always be hesitant to use machine learning predictions beyond their training range, specifically when dealing with problems for which we have minimal knowledge about the behaviour outside the training range.

Finally, we want to compare the performance of the symbolic regression algorithm with other popular machine algorithms in  literature, such as polynomial regression and artificial neural networks (ANN) for the same $\sin(x)$ case.  Polynomial regression may be considered as one of the oldest machine learning algorithms  [96], inspired by Legendre and Gauss's works, and  implemented in a robust algorithm by Gregonne in 1815 [97]. Polynomial regression  is the most appropriate “traditional” regression method to arrive at polynomials such as the Maclaurin series. In polynomial regression, the structure of the fitting equation  and the degree of the polynomial are predefined. For our case we will use two different polynomials one with a degree of $n$ = 3, and other one with $n$ = 7. We use the same training data set that we used for the symbolic regression, and for implementation, we will use the  Polyfit function from the open-source  Numpy library written in python [98]. The main output of the algorithm is the coefficients  of the following equation:

\begin{equation}\label{eq.6s}
y(x) = a_0+a_1x+... a_nx^n
\end{equation}

The coefficients for the two polynomials that we used are listed in Table <ref>.

We selected the artificial neural network because it is considered as a universal function approximators [99, 100], but also because it does not need any prior knowledge about the structure of the equation to best fit the training data, similar to the symbolic regression algorithm. Contrary to symbolic regression, the product of a neural network approach is not a function but the trained neural network itself. We will use a  feed-forward deep neural network, with eight hidden layers. The first hidden layer consists of 64 neurons, while, the remaining hidden layers contain 32 neurons, and finally an output layer containing a single neuron [96]. In each hidden layer we use the Relu activation function, and also we apply L2 regularization to avoid overfitting. The algorithm minimizes the mean square  difference between the predicted and training data,  using a  gradient descent algorithm.  We use the open-source library TensorFlow [101] to implement the artificial neural network framework. For training, we use 40,000 training points, which is a  much higher volume compared to the other two algorithms, because deep neural networks require a large amount of data to be trained appropriately[43].

A comparison between the performance of the  three  algorithms is shown in Figure <ref>. Symbolic regression and polynomial regression were the only algorithms that predict the peaks and valleys of the $\sin(x)$ function within the range of [-$\pi$,$\pi$]. This success can be attributed to the fact that both algorithms represent the $\sin(x)$ function as a  polynomial. For the case of the symbolic regression, it discovered the polynomial representation by itself. On the contrary, the ANN failed to generalize beyond the training data. We hoped that by making the network deeper, we could help the network extract sufficient features from the training data, and generalize. However, what we observe is that the ANN memorizes the training data instead of generalizing it. For example for $x$ $>$ $\dfrac{\pi}{2}$ the output of the ANN is always a constant value of one, which is the  value of $\sin(\dfrac{\pi}{2})$, and  for
$x<0$ the output of the ANN is always a constant value of  zero, which is the value of $\sin(0)$. This type of  memorization by an  ANN is also observed in several other studies such as [102]. Also, the work of Kim et al. [103] showed that if  feed-forward ANN is integrated  with symbolic regression, one obtains a better generalization behaviour compared to pure ANN. Another interesting observation  is that despite the fact that both symbolic regression and ANN optimize the mean  square difference, they come up with totally different generalization behaviour.

This Appendix A showed that symbolic regression can generalize beyond the training data, and can predict a  change in the original function occurring beyond the training range. This shows the usefulness of using interpretable machine learning results, as recommended by [104], and it helps us understand the output function behaviour within and beyond the training range.

[1]
Stokes, G. G.
On the effect of the internal friction of fluids on the motion
of pendulums, volume 9.
Pitt Press Cambridge, Cambridge, UK,  (1851).

[2]
Oseen, C. W.
Arkiv Mat., Astron. och Fysik 6, 1 (1910).

[3]
Oseen, C. W.
Ueber den gueltigkeitsbereich der stokesschen
Friedländer,  (1913).

[4]
Goldstein, S.
Proceedings of the Royal Society of London. Series A, Containing
Papers of a Mathematical and Physical Character 123(791), 225–235

[5]
Faxen, H.
Arkiv for Matemetik Astronomi och Fysik 17, 1–28 (1923).

[6]
Proudman, I. and Pearson, J.
Journal of Fluid Mechanics 2(3), 237–262 (1957).

[7]
Kaplun, S. and Lagerstrom, P.
Journal of Mathematics and Mechanics 6(5), 585–593

[8]
Graebel, W.
Advanced fluid mechanics.
Academic Press, Burlington, MA, USA,  (2007).

[9]
Chester, W., Breach, D., and Proudman, I.
Journal of Fluid Mechanics 37(4), 751–760 (1969).

[10]
Lagerstrom, P. and Reinelt, D.
SIAM Journal on Applied Mathematics 44(3), 451–462

[11]
Popović, N.
In Journal of Physics: Conference Series, volume 22,  164. IOP
Publishing,  (2005).

[12]
Van Dyke, M.
NASA STI/Recon Technical Report A 75 (1975).

[13]
Holzer, M., Kaper, T. J., et al.
Advances in Differential Equations 19(3/4), 245–282

[14]
Popović, N. and Szmolyan, P.
Nonlinear Analysis: Theory, Methods & Applications 59(4), 531–565 (2004).

[15]
Happel, J. and Brenner, H.
Low Reynolds number hydrodynamics: with special applications to
particulate media, volume 1.
Springer Science & Business Media,  (2012).

[16]
Cox, R.
Journal of Fluid Mechanics 23(4), 625–643 (1965).

[17]
Breach, D.
Journal of Fluid Mechanics 10(2), 306–314 (1961).

[18]
Aoi, T.
Journal of the Physical Society of Japan 10(2), 119–129

[19]
Achenbach, E.
Journal of Fluid Mechanics 54(3), 565–575 (1972).

[20]
Kamble, C. and Girimaji, S.
Physics of Fluids 32(10), 105110 (2020).

[21]
Dennis, S. and Walker, J.
Journal of Fluid Mechanics 48(4), 771–789 (1971).

[22]
Jenson, V.
Proceedings of the Royal Society of London. Series A.
Mathematical and Physical Sciences 249(1258), 346–366 (1959).

[23]
Nakhostin, S. and Giljarhus, K.
In IOP Conference Series: Materials Science and Engineering,
volume 700,  012007. IOP Publishing,  (2019).

[24]
Constantinescu, G., Pacheco, R., and Squires, K.
In 40th AIAA Aerospace Sciences Meeting & Exhibit,  425,

[25]
Deshpande, R., Kanti, V., Desai, A., and Mittal, S.
Journal of Fluid Mechanics 812, 815–840 (2017).

[26]
Maxworthy, T.
Journal of Applied Mechanics 36, 598–607 (1969).

[27]
Rouse, H.
Technical report, Dover Publications Inc, New York, USA,  (1961).

[28]
Engelund, F. and Hansen, E.
A monograph on sediment transport in alluvial streams.
TEKNISKFORLAG Skelbrekgade 4 Copenhagen V, Denmark.,  (1967).

[29]
Clift, R.
Proc. Chemeca'70 1, 14 (1970).

[30]
Morsi, S. and Alexander, A.
Journal of Fluid mechanics 55(2), 193–208 (1972).

[31]
Graf, W. H.
Hydraulics of sediment transport.
Water Resources Publication,Highlands Ranch, CO USA,  (1984).

[32]
Flemmer, R. L. and Banks, C.
Powder Technology 48(3), 217–221 (1986).

[33]
Khan, A. and Richardson, J.
Chemical Engineering Communications 62(1-6), 135–150

[34]
Swamee, P. K. and Ojha, C. S. P.
Journal of Hydraulic Engineering 117(5), 660–667 (1991).

[35]
Goossens, W. R.
Powder Technology 352, 350–359 (2019).

[36]
Schlichting, H. and Gersten, K.
Boundary-layer theory.
Springer, Berlin Heidelberg, Germany,  (2016).

[37]
Acrivos, A. and Taylor, T. D.
The Physics of Fluids 5(4), 387–394 (1962).

[38]
Rimmer, P. L.
Journal of Fluid Mechanics 32(1), 1–7 (1968).

[39]
Acrivos, A. and Goddard, J.
Journal of Fluid Mechanics 23(2), 273–291 (1965).

[40]
Whitaker, S.
AIChE Journal 18(2), 361–371 (1972).

[41]
Sparrow, E. M., Abraham, J. P., and Tong, J. C.
International Journal of Heat and Mass Transfer 47(24),
5285–5296 (2004).

[42]
Butler, K. T., Davies, D. W., Cartwright, H., Isayev, O., and Walsh, A.
Nature 559(7715), 547–555 (2018).

[43]
Brunton, S. L., Noack, B. R., and Koumoutsakos, P.
Annual Review of Fluid Mechanics 52, 477–508 (2020).

[44]
Dissanayake, M. and Phan-Thien, N.
communications in Numerical Methods in Engineering 10(3),
195–201 (1994).

[45]
Raissi, M. and Karniadakis, G. E.
Journal of Computational Physics 357, 125–141 (2018).

[46]
Iten, R., Metger, T., Wilming, H., Del Rio, L., and Renner, R.
Physical Review Letters 124(1), 010508 (2020).

[47]
Colen, J., Han, M., Zhang, R., Redford, S. A., Lemma, L. M., Morgan, L.,
Ruijgrok, P. V., Adkins, R., Bryant, Z., Dogic, Z., et al.
Proceedings of the national academy of sciences 118(10),
3932–3937 (2016).

[48]
Kushvaha, V., Kumar, S. A., Madhushri, P., and Sharma, A.
Journal of Composite Materials 54(22), 3099–3108 (2020).

[49]
Brunton, S. L., Proctor, J. L., and Kutz, J. N.
Proceedings of the national academy of sciences 113(15),
3932–3937 (2016).

[50]
Schmidt, M. and Lipson, H.
science 324(5923), 81–85 (2009).

[51]
Weng, B., Song, Z., Zhu, R., Yan, Q., Sun, Q., Grice, C. G., Yan, Y., and Yin,
Nature communications 11(1), 1–8 (2020).

[52]
El Hasadi, Y. M. and Padding, J. T.
AIP Advances 9(11), 115218 (2019).

[53]
Koza, J. R.
Genetic programming: on the programming of computers by means of
natural selection, volume 1.
MIT press, Cambridge, MA, USA,  (1992).

[54]
Holland, J. H. et al.
Adaptation in natural and artificial systems: an introductory
analysis with applications to biology, control, and artificial intelligence.
MIT press,Cambridge, MA, USA,  (1992).

[55]
Darwin, C.
On the origin of species.
Johon Murray, London, UK,  (1859).

[56]
Brown, P. P. and Lawler, D. F.
Journal of environmental engineering 129(3), 222–231

[57]
Schiller, L. and Naumman, A.
Z. Vereines Deutscher Inge. 77, 318–321 (1933).

[58]
Blasius, H.
Grenzschichten in Flüssigkeiten mit kleiner Reibung.
University of Gottingen, Germany,  (1908).

[59]
Prandtl, L.
In Verhandlg. III. Intern. Math. Kongr. Heidelberg,  484–491.
Math Congress,  (1904).

[60]
Leal, L. G.
Advanced transport phenomena: fluid mechanics and convective
transport processes, volume 7.
Cambridge University Press, Cambridge, UK,  (2007).

[61]
Abraham, F. F.
The Physics of Fluids 13(8), 2194–2195 (1970).

[62]
Brauer, H. and Mewes, D.
Chemie Ingenieur Technik 44(13), 865–868 (1972).

[63]
Hölzer, A. and Sommerfeld, M.
Powder Technology 184(3), 361–365 (2008).

[64]
Batchelor, G.
An introduction to fluid dynamics.
Cambridge university press, Cambridge, UK,  (2000).

[65]
Hunter, C., Tajdari, M., and Boyer, S.
SIAM Journal on Applied Mathematics 50(1), 48–63 (1990).

[66]
Domingos, P.
Data mining and knowledge discovery 3(4), 409–425

[67]
Nisi, H. and Porter, A. W.
The London, Edinburgh, and Dublin Philosophical Magazine and
Journal of Science 46(275), 754–768 (1923).

[68]
Rimon, Y. and Cheng, S.
The Physics of Fluids 12(5), 949–959 (1969).

[69]
Bourot, J.
A 269, 1017–1020 (1969).

[70]
Payard, M. and Coutanceau, M.
CR Academie des Sciences, Paris B 278, 369–372 (1974).

[71]
Chang, E. J. and Maxey, M. R.
Journal of Fluid Mechanics 277, 347–379 (1994).

[72]
Taneda, S.
Journal of the Physical Society of Japan 11(10),
1104–1108 (1956).

[73]
Van Dyke, M.
Journal of Fluid Mechanics 44(2), 365–372 (1970).

[74]
Ranger, K.
SIAM Journal on Applied Mathematics 23(3), 325–333

[75]
Maxworthy, T.
Journal of Fluid Mechanics 23(2), 369–372 (1965).

[76]
John Veysey, I. and Goldenfeld, N.
Reviews of Modern Physics 79(3), 883 (2007).

[77]
Goldstein, S.
Aeronautical Research Council (Great Britain)  (1965).

[78]
Lewis, J. and Carrier, G.
Quarterly of Applied Mathematics 7(2), 228–234 (1949).

[79]
Beard, K. and Pruppacher, H.
Journal of the Atmospheric Sciences 26(5), 1066–1072

[80]
Le Clair, B., Hamielec, A., and Pruppacher, H.
Journal of the Atmospheric Sciences 27(2), 308–315

[81]
Beratlis, N., Balaras, E., and Squires, K.
Journal of Fluid Mechanics 879, 147–167 (2019).

[82]
Muto, M., Tsubokura, M., and Oshima, N.
Physics of Fluids 24(1), 014102 (2012).

[83]
Geier, M., Pasquali, A., and Schönherr, M.
Journal of Computational Physics 348, 889–898 (2017).

[84]
Suryanarayana, G., Pauer, H., and Meier, G.
Experiments in Fluids 16(2), 73–81 (1993).

[85]
Morrison, F. A.
An introduction to fluid mechanics.
Cambridge University Press,  (2013).

[86]
Wieselsberger, C. v.
physics journal 23, 219–224 (1922).

[87]
Schultz, M. P. and Flack, K. A.
Physics of Fluids 25(2), 025104 (2013).

[88]
Udrescu, S.-M. and Tegmark, M.
Science Advances 6(16), eaay2631 (2020).

[89]
Feynman, R. P., Leighton, R. B., and Sands, M.
American Journal of Physics 33(9), 750–752 (1965).

[90]
Duan, Z., He, B., and Duan, Y.
Scientific reports 5, 12304 (2015).

[91]
Will, J. B., Kruyt, N. P., and Venner, C. H.
International journal of Heat and Mass transfer 109,
1059–1067 (2017).

[92]
Feng, Z.-G. and Michaelides, E. E.
International Journal of Heat and Mass Transfer 43(2),
219–229 (2000).

[93]
Taylor, B.
Methodus incrementorum directa & inversa.
Inny,  (1717).

[94]
Newton, I.
Philosophiae naturalis principia mathematica, volume 1.
G. Brookman,  (1833).

[95]
Leibniz, G. W.
Acta eruditorum , 467–473.

[96]
Brunton, S. L. and Kutz, J. N.
Data-driven science and engineering: Machine learning, dynamical
systems, and control.
Cambridge University Press,  (2019).

[97]
Stigler, S. M.
Historia Mathematica 1(4), 431–439 (1974).

[98]
www.numpy.org/doc/stable/user/ .

[99]
Cybenko, G.
Mathematics of control, signals and systems 2(4),
303–314 (1989).

[100]
Hornik, K.
Neural networks 4(2), 251–257 (1991).

[101]
Abadi, M., Agarwal, A., Barham, P., Brevdo, E., Chen, Z., Citro, C., Corrado,
G. S., Davis, A., Dean, J., Devin, M., et al.
arXiv preprint arXiv:1603.04467  (2016).

[102]
Zhang, C., Bengio, S., Hardt, M., Recht, B., and Vinyals, O.
arXiv preprint arXiv:1611.03530  (2016).

[103]
Kim, S., Lu, P. Y., Mukherjee, S., Gilbert, M., Jing, L., Čeperić,
V., and Soljačić, M.
IEEE Transactions on Neural Networks and Learning Systems  (2020).

[104]
Rudin, C.
Nature Machine Intelligence 1(5), 206–215 (2019).

to 0.8 | X[c] | X[c] | X[c] | X[c] |
Coefficients      Eq.(<ref>)     Eq.(<ref>)    Eq.(<ref>)

$a_1$      0.251      0.412    0.505

$a_2$      23.620      23.311   23.224

$a_3$      0.001      4.119   2.762

$a_4$      3.255      -     -

$a_5$      49.291      -     -

$a_6$      97.537     -      -

$a_7$      -$2.709\times 10^{-6}$      -     -

Coefficients for Eq.(<ref>)  Eq.(<ref>), and Eq.(<ref>)

to 0.8 | X[c] | X[c] | X[c]  | X[c] |
Coefficients      Ref[62]     Ref[63]     Ref [61]

$a_1$      2.9%      -1.94%     29.01%%

$a_2$      -2.95%      -2.95%     -2.87%

$a_3$      2.88%      27.16%      -28.40%

Relative difference in the values of coefficients of Eq.(<ref>) to that of Brauer and Mewes [62],  Holzer and Sommerfeld [63], and Abraham[61].

to 0.8 | X[c] | X[c] | X[c] |
Coefficients      Eq.(<ref>)     Eq.(<ref>)

$a_1$      3.286     3.272

$a_2$      24.205      23.26

$a_3$      -0.818      0.112

$a_4$       0.064     -0.652

$a_5$      -0.000107      0.035

Coefficients for Eq.(<ref>) and Eq.(<ref>)

to 0.8 | X[c] | X[c] |
Coefficients      Eq.(<ref>)

$a_1$      3.140

$a_2$      24.270

$a_3$      -0.716

$a_4$       0.047

Coefficients for Eq.(<ref>)

to 0.8 | X[c] | X[c] | X[c] | X[c] |
Coefficients      Eq.(<ref>)     Eq.(<ref>)    Eq.(<ref>)

$a_1$      2.0     1.582      1.063

$a_2$       0.343    0.003     0.0067

$a_3$      0.0454     0.326      1.351

$a_4$       9.341     1.0     0.299

$a_5$      1.0      0.322      0.0028

$a_6$      $-7.0\times10^{-5}$     -    0.332

$a_7$      -0.00131      -     -0.128

Coefficients for Eq.(<ref>), Eq.(<ref>), and Eq.(<ref>)

to 0.8 | X[c] | X[c] | X[c] |
Coefficients      Eq.(<ref>)     Eq.(<ref>)

$a_1$      0.9999     1.0001

$a_2$      0.1665      0.1682

$a_3$      0.00826      0.0031

$a_4$      0.000173     0.0065

Coefficients for Eq.(<ref>) and Eq.(<ref>)

to 0.8 | X[c] | X[c] | X[c] |
Coefficients      $n$ = 3     $n$ = 7

$a_0$      -0.002     -4.70$\times 10^{-8}$

$a_1$      1.027      1.0

$a_2$      -0.069      -2.339$\times 10^{-5}$

$a_3$      -0.138    -0.166

$a_4$	   -    -2.45$\times 10^{-4}$

$a_5$	   -    0.008

$a_6$ 	   -    -2.046$\times 10^{-4}$

$a_7$ 	   -    -1.377$\times 10^{-4}$

Coefficients of polynomials of degree $n$ = 3,and 7

Comparison  between the drag coefficient $C_D$ predicted by  Eq.(<ref>), Eq.(<ref>),  Eq.(<ref>), Eq.(<ref>)and, different sources from the literature. Dashed lines indicate literature correlations. Symbols indicate experimental values.

Comparison  between drag coefficient $C_D$ predicted by  Eq.(<ref>)  Eq.(<ref>), and different sources from the literature. Dashed lines indicate literature correlations. Symbols indicate experimental values.

Schematic of separated flow around a sphere

Comparison between the $C_D$ predictions in the low $Re$ limit by Eq.(<ref>), Eq.(<ref>), Eq.(<ref>),Eq.(<ref>), and Eq.(<ref>), and  different sources from the literature for  low $Re$ regime.  Circles represents experiments, and squares represents numerical simulations.

Comparison between the $C_D$ predictions by Eq.(<ref>),Eq.(<ref>), and different sources in  the high $Re$ regime where the drag crisis occurs.  Circles represents experiments, and squares represents numerical simulations.

Comparison  between the results of different predictor equations  for the Nusselt number $Nu$   with those  from  literature for four different Prandtl numbers $Pr$.

The influence of different terms of Eq.(<ref>) on its variation with $x$. Blue bars indicate the range of training data.

Comparison between different machine learning methods for the $\sin(x)$ example. Blue bars indicate the range of training data.