Variable-Order Fractional Models for Wall-Bounded Turbulent Flows

Modeling of wall-bounded turbulent flows is still an open problem in classical physics, with relatively slow progress in the last few decades beyond the log law, which only describes the intermediate region in wall-bounded turbulence, i.e., 30–50 y+ to 0.1–0.2 R+ in a pipe of radius R. Here, we propose a fundamentally new approach based on fractional calculus to model the entire mean velocity profile from the wall to the centerline of the pipe. Specifically, we represent the Reynolds stresses with a non-local fractional derivative of variable-order that decays with the distance from the wall. Surprisingly, we find that this variable fractional order has a universal form for all Reynolds numbers and for three different flow types, i.e., channel flow, Couette flow, and pipe flow. We first use existing databases from direct numerical simulations (DNSs) to lean the variable-order function and subsequently we test it against other DNS data and experimental measurements, including the Princeton superpipe experiments. Taken together, our findings reveal the continuous change in rate of turbulent diffusion from the wall as well as the strong nonlocality of turbulent interactions that intensify away from the wall. Moreover, we propose alternative formulations, including a divergence variable fractional (two-sided) model for turbulent flows. The total shear stress is represented by a two-sided symmetric variable fractional derivative. The numerical results show that this formulation can lead to smooth fractional-order profiles in the whole domain. This new model improves the one-sided model, which is considered in the half domain (wall to centerline) only. We use a finite difference method for solving the inverse problem, but we also introduce the fractional physics-informed neural network (fPINN) for solving the inverse and forward problems much more efficiently. In addition to the aforementioned fully-developed flows, we model turbulent boundary layers and discuss how the streamwise variation affects the universal curve.


Introduction
Reynolds [1] was the first to statistically describe turbulence by decomposing the instantaneous velocity vector into an average field and its fluctuation. Upon substitution into the Navier-Stokes equations and averaging, assuming quasi-stationarity, a new modified equation emerged for the average velocity that includes an additional term, namely, the averaged dissipation tensor leading to the turbulence-closure problem [2]. Addressing the closure complexity has been a century-long pursuit, starting with the seminal work of Prandtl [3], who proposed a simplified mixing length model analogous with Fick's law of local diffusion. Interestingly, at about the same time, Richardson [4], in an attempt to unify turbulent diffusion with molecular diffusion, combined geophysical measurements with Brownian motion to produce the famous scaling law on turbulent pair diffusivity. While ingenious, both approaches assume implicitly locality in turbulent interactions, which limits the universality of the derived correlations-an open standing question for over a century. As stated by Kraichnan [5], Prandtl's approach is valid only when the

Variable-Order Fractional Models for Turbulent Flows
The first fractional model for the Reynolds averaged Navier-Stokes equations was developed by Chen [11], who proposed a fractional Laplacian to model the Reynolds stresses and to account for intermittency [18,19] as follows: where U is the average velocity and γ is the turbulent diffusion coefficient. Hence, the effective fractional order in this model is fixed at α = 2/3. This value is consistent with the Richardson superdiffusion scaling for homogeneous turbulence that leads to a t 3 scaling for the mean square displacement, but it is not valid for wall-bounded turbulence where anisotropy and the distance from the wall determine the effective rate of turbulent diffusion. Defining a fractional Laplacian in multiple dimensions and in bounded domains is still an open issue in fractional calculus and extending it to variable orders is challenging [15]. However, other somewhat equivalent definitions based on tempered fractional calculus [20] may lead to satisfactory nonlocal representations as well; specifically, in a Boltzmannian framework, Samiee et al. [13] developed a tempered fractional subgrid-scale model to capture high-order structures at the inertial and dissipative ranges. As Richardson first noted, the velocity field in the atmosphere shares a number of properties with the Weierstrass function, i.e., it appears to be continuous but non-differentiable, and this provides a strong case for fractional modeling of turbulence in the atmosphere but also in wall-bounded flows in engineering applications.
In this section, we present a variable-order fractional model for turbulent flows. We firstly consider a one-sided model for channel and pipe flows. Furthermore, we formulate an inverse problem for the fractional order α(y). We present a finite difference method and design a physics-informed neural network (PINN) to obtain the fractional order. Finally, we propose a divergence variable fractional (two-sided) model for turbulent flows.

One-Sided Fractional Derivative Modeling
For wall-bounded turbulence, the effective rate of diffusion varies with distance from the wall. Hence, we exploit the power of fractional calculus that allows variable fractional order, and we propose a variable-order fractional differential equation for modeling the Reynolds stresses, i.e., α(y), where y is the distance from the wall. In particular, we consider fully developed turbulent flows with one-dimensional (dimensionless) averaged velocity U(y) = u/V (where V is the characteristic velocity), including channel flows and pipe flows for which we apply a unified fractional modeling approach. Specifically, assuming that the flow direction is along x and y is the wall-normal direction (distance from the wall), we consider the variable fractional model (VFM-I) in the normalized interval [0, 1]: with α(0) = 1, 0 ≤ α(y) ≤ 1, D α y is the (Caputo) fractional derivative, f = − 1 ρ ∂P/∂x is a constant pressure gradient, U(y) is the mean velocity we want to model, and ν 0 is the kinematic viscosity. The Caputo derivative is defined as: and it is identical to the Riemann-Liouville left-sided derivative because U(0) = 0. Interestingly, we can obtain the scalar coefficient ν(y) (we refer to it as turbulent diffusivity, although it does not have the correct units) explicitly in terms of the fractional order α(y) from: where Re τ = u τ R/ν 0 is the friction Reynolds number, R is the radius of the pipe (or the half channel width), and u τ is the wall friction velocity u τ = τ w /ρ, where τ w = µ∂U/∂y| y=0 is the wall shear stress with µ being the dynamic viscosity.
We discuss an alternative model, where the variable fractional order α(y) is between one and two instead of the VFM-I we presented, where 0 < α(y) ≤ 1; this model is analogous to VFM-I and is defined by: with α(0) = 2, and the variable-order 1 ≤ α(y) ≤ 2 is an unknown function to be determined by the data. The scalar coefficient ν(y) can also be computed from a similar formula as before, i.e.,

Numerical Method
We assume that we know the mean velocity U(y) (also U + (y + )) from the DNS data or experimental results. The VFM-I can be written in the form: where f = − 1 ρ ∂P/∂x. Since the fractional order α(y) is unknown in Equation (6), we need to solve a nonlinear problem to obtain α(y). Alternatively, we consider the following optimization problem: given U and f , find the α(y) that satisfies where, S(Λ) := {0 ≤ a(y) ≤ 1, a(y) ∈ C 0 (Λ)}. If α * (y) satisfies Equation (6), then we obtain J(α * (y)) ≡ 0. Next, we present a numerical method for solving the optimization problem (7). The fractional derivative is discretized with the finite difference method. Then, the fractional order α(y) can be solved point-by-point; for each point y n = n∆y, ∆y = 1/N, n = 1, 2, · · · , N, we calculate the fractional derivative D α(y n ) y U n with the DNS data using the finite difference method [21] D α(y n ) y where b n j := (j + 1) α(y n ) − j α(y n ) and U n = U(y n ). The discrete optimization problem can now be written as Finally, we formulate the fractional physics-informed neural network (fPINN) for the inverse problems of the proposed turbulence model; see Figure 1.
The aim of the inverse problem is to estimate the fractional order α(y) given the mean velocity profile U in the DNS data. We approximate the variable fractional order α(y) by a multi-layer feedforward neural network α NN (y; θ = {W j , b j } l j=1 ), where θ are a collection of parameters of the NN. The locations y are the input of the NN, and the output U is computed by a recursive formula The weight matrix between the (j − 1)th and jth layers has the dimension W j ∈ R n j ×n j−1 , and the bias vector b j in the jth layer. The column vectors Y j−1 ∈ R n j−1 ×1 and Y j ∈ R n j ×1 denote the input and output of the jth layer, respectively. The input vector Y j−1 is first subject to a linear transformation and then an element-wise nonlinear function σ(·), which is called the activation function. The NN consists of one input layer (j = 0), l − 1 hidden layers (j = 1, 2, · · · , l − 1), and one output layer (j = l). The depth of the NN is l, and the width of the jth layer is n j . To determine the parameters θ, we minimize the following loss function with respect to θ The first term on the right-hand side is the equation residual, and the second term is the constraint on the fractional order at the wall, i.e., α(0) = 1. We select N t training points, {y i } N t i=1 , to enforce the equation residual on them to be zero. The fractional derivative is evaluated using the finite difference method (8). We optimize the loss function with respect to θ, employing a stochastic gradient descent, Adam, written in TensorFlow. Finally, we estimate the variable fractional order using α NN (y; θ). Figure 1. Basic structure of fPINN in 1D for the inverse fractional-order problem. The left uninformed DNN processes data to predict the fractional order, which also has to satisfy the correct physics of turbulence for the channel fully developed flow, represented by the right informed DNN induced by the fractional governing equation.

Fractional Modeling in Divergence Form
We consider the Reynolds averaged momentum equation for incompressible fully developed channel flow; the governing equation is as follows where ρ is the density; and P and U are the mean pressure and velocity, respectively. The process of Reynolds averaging introduces the unclosed Reynolds stress, τ ij = −ρu v . The total shear stress on the wall is τ w . Integrating the above equation from wall to an arbitrary position in wall-wise y, we obtain a new formula as follows We assume the dimensionless wall shear τ w and pressure gradient ∂P ∂x = C are constants. Additionally, we introduce a symmetric divergence variable fractional model for approximating the total shear stress, with the boundary conditions α(0) = α(2) = 1, where the fractional derivative is defined as follows D α(y) and D α(y) y and y D α(y) U are left and right Caputo derivatives, respectively. The definitions are given as follows and it is identical to the Riemann-Liouville derivatives because U(0) = 0 and U(2) = 0. We also propose the eddy viscosity in the fractional momentum equation, and the explicit formula is as follows where Re τ = u τ R/ν 0 is the friction Reynolds number, R is the radius of the pipe (or the half channel width), and u τ is the wall friction velocity, u τ = τ w /ρ, where τ w = µ∂U/∂y| y=0 is the wall shear stress with µ being the dynamic viscosity.

Numerical Method
We assume that we know the mean velocity U(y) (also U + (y + )) from the DNS data or experimental results. Since the fractional order α(y) is unknown in Equation (13), we need to solve a nonlinear problem to obtain α(y). Alternatively, we consider the following optimization problem: given U and f , find α(y) that satisfies where f = 1 − y and S(Λ) := {0 ≤ a(y) ≤ 1, a(y) ∈ C 0 (Λ)}. If α * (y) satisfies Equation (13), then we obtain J(α * (y)) ≡ 0. Next, we present a numerical method for solving the optimization problem (16). The fractional derivative is discretized with the finite difference (FD) method. Then, the fractional order α(y) can be solved point-by-point; for each point y n = n∆y, ∆y = 1/N, n = 1, 2, · · · , N, we calculate the fractional derivatives D α(y n ) |y| U n with the DNS data using the finite difference method [21] Left: D and where b n j := (j + 1) α(y n ) − j α(y n ) , c n j = b n j and U n = U(y n ). The discretized optimization problem can be now written as Here, we use N ≈ Re τ points to solve the above optimization for the channel flow at a given Reynolds number Re τ .
Alternatively, we propose the fractional fPINN for solving the inverse DVFM with the loss function

Turbulent Boundary Layer and Couette Flow
For a boundary layer and Couette flow with zero pressure gradient, the mean twodimensional continuity and stream-wise momentum reduce to ∂UU ∂x If we assume that the convective effects are small near the wall for the boundary layer problem, then the above equation reduces to Here, U is viewed as a function of y due to ∂U ∂x = 0. Since the two plates are infinitely long for the Couette flow, the flow properties cannot change with x and all partial derivatives with respect to x vanish. Flow motion only occurs in the x direction, and thus, V = 0. After simplifying the RANS equations, the turbulent Couette flow is governed by Equation (22) too.
Further integrating the above equation provides where C is a constant and uv = 0 at the wall, while ν ∂U ∂y is simply the wall shear stress τ w /ρ. Then, we have the following equation with α(0) = 1, 0 < α ≤ 1, D α y is the (Caputo) fractional derivative, and ν(y) is the eddy viscosity defined as

Numerical Method
We solve the fractional order α(y) for the turbulent boundary layer problem and Couette flowy using fPINN ( Figure 1) with the loss function where U is the DNS data. It changes with Re θ for the boundary layer problem, so there is (implicit) x dependence as well.

Numerical Results
In this section, we present the results for the turbulent channel, pipe, Couette, and boundary layer flows.

Numerical Results of the One-Sided Models
We first consider turbulent channel flow for which DNS data are available up to Re τ = 5200 [22]. Here, we use the FD scheme with N ≈ Re τ points to solve the aforementioned inverse problem for the channel flow at a given Reynolds number Re τ . Solving for α(y), which uniquely determines the Reynolds stresses, Figure 2a depicts the profiles of the fractional order α(y) for different Re τ as a function of the non-dimensional distance from the wall y ∈ [0, 1]. We see a strong dependence of α(y) on Re τ ; however, if we re-plot all data in terms of the viscous wall units, i.e., y + = yu τ /ν 0 we see a collapse of all results into a single universal curve, as shown in Figure 2b. Moreover, we employ the empirical Spalding formula [23] for U + = u/u τ in order to extend the results up to high Re τ = 10 6 , and again we obtain a similar universal scaling with the exception of low Re τ for which the Spalding formula is known to be somewhat inaccurate. We fit the fractional order using these numerical results to obtain the fractional order α(y + ) in wall units as follows where φ(y + ) = tanh(ln(y + /9.5)/1.049) and a(y + ) = 1/(b + κ| ln(y + )| 0.9 ) with b = 0.855, κ = 0.301 are constants. This is a remarkable result as it goes beyond the logarithmic profile and seamlessly connects the viscous sublayer with the buffer zone, the logarithmic profile, and the wake region. Although at first it appears to be a perfect fitting exercise, it has important consequences due to the nonlocal interpretation of the fractional derivative involved, i.e., it shows that nonlocality is stronger away from the wall and at high Reynolds numbers. Using the same data for U(y), we show that the alternative model VFM-II with 1 ≤ α(y) ≤ 2 also leads to the same type of universality ( Figure 3). However, unlike the aforementioned VFM-I, we are unable to obtain an explicit formula for ν(y), relating it to the Reynolds number as in the first model (i.e., α(y) ∈ (0, 1]); instead, we can compute it numerically from the DNS data of turbulent channel flow. As shown in Section 3, this alternative fractional model also exhibits a universal scaling if plotted in terms of wall units, with the lowest value of α(10 5 +) ≈ 1.3.  To evaluate the predictability of the universal scaling, we now solve the forward Equation (2) to obtain U(y) at Re τ = [4200, 6000, 8600], which are cases not used in the training of the model for α(y + ). The results presented in Figures 4 and 5 are in good agreement with DNS and experimental data. We also include the turbulent channel flow results obtained by nested LES [24]. Figures 4 and 5 show that the mean velocity profiles predicted by VFM-I exhibit the correct behavior throughout the channel for Reynolds numbers up to Re τ = 8600, including the correct slope in the logarithmic layer, and agree with DNS and experimental data in the wake region for all Re τ = [4200, 6000, 8600].   We used fPINN to investigate the turbulent channel flows. We used different training points for investigating the convergence using DNS data at Re τ = 2000. Figure 6 shows the training results with uniform training points in the interval for N t = 500, 1000, 2000. Figure 7 shows the training results with log-uniform training points in wall units scaling for N t = 10, 20, 40, 80. The corresponding loss histories are listed in Table 1. Figure 7 presents the comparison profiles between the training sets. We can observe that the results trained by the log-uniform are smoother than the uniform training points near the wall. (a) α(y) for N t = 500.      (c) α(y) for N t = 40.  Next, we test the accuracy of the forward problem and the loss function error with the training fractional order predicted by log-uniform training points N t = 20. We solve the fractional equation as follows: with U(0) = 0, and the fractional orderwasis obtained by training fPINN with N t = 20 and Equation (24). The corresponding loss functional error is defined as follows Finally, we use the simplified one-dimensional equation where the R uv denotes the Reynolds stress R uv = u v , τ uv denotes the viscous shear stress τ uv = ν 0 ∂U/∂y, and U is the mean velocity, which is the solution to the above fractional Equation (26). Then, we obtain the Reynolds stresses by integration, We can compare the predicted Reynolds stresses to their counterparts, R D from DNS data for turbulent channel flow, and using the corresponding viscous shear stress denoted by τ D = µ∂U D /∂y, where U D denotes the mean velocity from the DNS database. In Figure 9, we plot the predicted and DNS profiles for Reynolds numbers Re τ = 4000, 5200 and the corresponding pointwise error. We can observe that they are all in very good agreement. The numerical results of the mean velocities and shear stresses for all Reynolds number Re τ match very well with the DNS data; here, we only show the high Reynolds number cases due to space limitations.   Here, τ uv denotes the wall shear stress for the fractional order predicted by the finite difference (FD) method, τ NN uv denotes the wall shear stress predicted by the NN, and τ D is the corresponding profile from DNS data. −R uv denotes the Reynolds shear stress predicted by Equation (24), −R uv denotes the wall shear stress predicted by the NN, and −R D is the corresponding profile from DNS data.

Numerical Results of the Two-Sided Models
In this subsection, we focus on the two-sided models. Solving for α(y), which uniquely determines the total shear stresses, Figure 10 plots the profiles of the fractional order α(y) for different Re τ as a function of the non-dimensional distance between the two walls y ∈ [0, 2]. We see a strong dependence of α(y) of Re τ , which is the same conclusion as for the previous variable fractional model. Furthermore, we re-plot all data in terms of the viscous wall units, i.e., y + = yRe τ , and we see an approximate collapse of all results into a single universal curve in the half-plane excluding the wake region (i.e., near the centerline), as shown in Figure 10. Next, we test the accuracy of the forward problem with the fractional order provided by the inverse optimization problem (19). We solve the divergence variable fractional equation as follows − D ν(y)D α(y) |y| U = 1, ∀y ∈ (0, 2), (28) with U(0) = U(2) = 0. Figure 11 plots the solutions (left) of the above equation and the pointwise error (right) of the mean velocity in each subfigure for several Re τ . We can observe that this model predicts the mean velocity well. Moreover, it can obtain a smooth mean velocity profile in the whole domain along the wall-wise direction.
We also use fPINN (20) to solve the inverse problem to obtain the variable order α(y). The two results from the two different methods (i.e., FD and fPINN) agree well for all Reynolds numbers.

Turbulent Pipe Flow
In this subsection, we consider turbulent pipe flow and again test the universal variable fraction order α(y + ) against DNS and experimental data. First, we examine the highest Reynolds number available from the superpipe experiment [27,28] at Re τ = 5 × 10 5 , estimated at Re R ≈ 3.525 × 10 7 based on the pipe radius R. As the experimental data were only available for y + > 10,000, we synthesized an entire profile from the pipe wall to centerline using multifidelity Gaussian process regression (M-GPR) [29] as follows: we considered as high fidelity data the superpipe data in the outer region together with the highest DNS data for channel flow at Re τ = 5200. We then employed the Spalding curve to provide the low-fidelity data and, using M-GPR, we constructed the final profile as shown in Figure 12a. Having this profile and the VFM-I model transformed in polar coordinates, we can then solve the inverse problem and obtain a new variable fractional order α(y + ). Figure 13a shows that the variable fractional order we obtain for this problem is identical to the function defined by Equation (24). This finding further confirms the universality of the variable fractional order even at very high Reynolds numbers. Having validated the accuracy of the variable fractional order, we can now solve the forward fractional differential problem to obtain predictions of the entire velocity profiles from Re τ = 10 5 to Re τ = 5 × 10 5 . Figure 12b plots the results, showing that there is excellent agreement with all available data from the superpipe experiment. Figure 13b plots the mean velocity profiles from the DNS data base [30] at low Reynolds numbers, the corresponding VFM predictions, and the Spalding profile. The universal defect law for pipe flows is not valid for the low Reynolds number range, and this is also in agreement with [27], who argued that the lowest Re τ for universality is approximately 5000.

Turbulent Couette Flow
In reference [12], the authors proposed the double-log profile to predict the mean velocity for the Couette flow as follows where d is a small number (d 1) that represents a viscous sublayer or roughness height. The non-dimensional boundary conditions are U(0) = 0 and U(1) = 1.
Here, we consider the predictions from the universal scaling fractional order α * (y + ), and we also compare it against the double-log profile. The variable fractional order α * (y + ) is between zero and one in our turbulence model. So, we work in the half-plane y ∈ [0, 0.5] (see the dashed square in Figure 14a). We then obtain the results in the other half of the domain with U(y) = 1 − U(1 − y), y ∈ (0.5, 1]. Figure 14 shows the mean velocity profiles predicted using (29) and the mean velocity, which is predicted by the variable fractional order α * (y + ). We can observe that the variable fractional model is in agreement with the experiment data as well as the double-log profile. However, the double-log profile is unable to capture the correct mean velocity near the wall. We also tested the profiles for low Reynolds number Re τ = 52, where the numerical data were obtained from reference [31]. For the double-log profile, we could not find a suitable parameter d to obtain a good fit for the low Re τ = 52. Finally, we show the comparisons between the TCM predicted mean velocities and DNS data at Re τ = 250 obtained from reference [32]. Figure 15 shows that the fractional predictions are correct almost everywhere, especially near the wall regions for high Reynolds numbers.

Turbulent Boundary Layer Flow
In this subsection. we focus on the boundary layer problem. We use data available from the KTH turbulence group from the turbulent boundary layer DNS [34,35]. We first investigate the correlations between Re θ (x-variable) and Re τ (y-variable); Figure 16 shows the downstream variations in the friction Reynolds number Re τ , and unlike the channel flow, here, Re τ is a function of the streamwise distance x.
Then, we test if the mean velocity of the boundary layer problem exhibits any universality as the channel and pipe flow. We solve the forward boundary layer problem with the fractional order predicted by Equation (24) (i.e., the formula is the same as the channel flow case) including the wake region. Figure 17 presents the mean velocity profiles from the DNS [34] and fractional modeling near the wall for several Re θ from 670 to 4060, with the corresponding Re τ varying from 252 to 1200. We observe that the mean velocities are different in the wake region for different Re τ . Figure 18 plots the wake region, which is between δ + 99 and the error E = 1%. We define this error as the difference in the mean velocity between the DNS data and the fractional model as follows: where U is the DNS data and U f presents the numerical results from the fractional model.           Since the mean velocity does not exhibit universality in the wake region, we solve the fPINNs to investigate the variations in the fractional order in the wake region. In Figure 19, we plot the fractional order inferred by fPINN based on the DNS data for Re θ = 670 to 4060. We can observe that the fractional order varies for different Re θ in the wake region. Then, we train the fractional order in the wake region selecting the data set Re θ = 670 to 4060 but excluding Re θ = 2000. In Figure 20, we present the factional order in the 2D plane for the x-axis and y + -axis. Finally, we solve the fractional turbulent boundary layer model with the fractional orders presented in Figure 20. The comparison between the numerical results and the DNS data set is presented in Figure 21.

Summary
We proposed multiple fractional models for wall-bounded turbulent flows in benchmark cases where the mean flow is either one-dimensional (channel, pipe, and Couette flows) or two-dimensional (boundary layer). The main idea is to employ a variable-order fractional gradient that depends on the distance from the wall, starting with an integer order at the wall. The computational problem we addressed is the discovery of the fractional variable-order profile given DNS or experimental data for the mean velocity profile. To this end, we formulated an inverse problem for the fractional order as a function of the distance from the wall, and we solved it using a finite difference method point-by-point and through a new fractional physics-informed neural network (fPINN) that encodes the physics of turbulence expressed via the fractional derivative of variable order. The fractional order is a function of the distance from the wall, a unique capability enabled by fractional calculus. We discovered that this fractional order function is universal for all Reynolds numbers and for different geometries.
The main contributions of this work are: (1) new fractional turbulent models with variable order are presented to model the total shear stress of RANS; (2) two solution methods for the non-trivial inverse problem, a FD method, and a fPINN for obtaining the fractional order function; (3) a universal fractional order profile was discovered for the channel and pipe flows that allowed us to accurately predict the fractional order for the boundary layer flows.