Physics-to-Geometry Transformation to Construct Identities between Reynolds Stresses

: Modeling has become ﬁrmly established as a methodology to close the Reynolds-averaged Navier–Stokes (RANS) equations, owing to theoretical and empirical efforts towards a complete formulation of the Reynolds stress tensor and, recently, breakthroughs in data-processing technology. However, mathematical exactness is not generally ensured by modeling, which is an intrinsic reason why the reliability of RANS closure models is not supposed to be consistent for all kinds of turbulent ﬂow. Rather than straightforwardly overcoming this inherent limitation, most of the studies to date were reasonably directed towards broadening the range of turbulent ﬂows, where reliable prediction accuracy can be obtained via modeling. In this paper, we present three identities between components of the Reynolds stress tensor, constructed via spatial mapping on the basis of the differential version of the Gauss–Bonnet formula. Further, we present a constraint condition that gives a set of equations as numerous as the parameters within a RANS model.


Introduction
In deriving the Reynolds-averaged Navier-Stokes (RANS) equations by the application of the Reynolds decomposition to the Navier-Stokes equations, six unclosed terms, known as Reynolds stresses, appear in the derivation process [1]. As a methodology to construct a closure for these unclosed terms, modeling has been firmly established so far. In a typical process to model the Reynolds stress tensor, the construction of a mathematical framework for model representation relies on physical intuition and empiricism [2]. Furthermore, the calibration of model parameters is indeed limited in practice to a few benchmark turbulent flows. Therefore, it should be underlined that the reliability of RANS closure models is not supposed to be consistent for all kinds of turbulent flow.
Mathematical exactness is not generally ensured by modeling, which is an intrinsic reason why RANS closure models inevitably have uncertainties in predictive capability. Examples include one-and two-equation eddy-viscosity models [3][4][5][6], and the isotropy eddy-viscosity assumption equipped in these models was empirically found to not accurately predict Reynolds stress anisotropy for flow over curved surfaces [7]. Turbulent flows with Reynolds stress anisotropy are common in engineering applications, and RANS closure models were correspondingly developed by use of more adaptive functional forms or mathematical frameworks; examples include quadratic and cubic k − models [8,9], V2F model [10], algebraic stress models [11][12][13][14][15][16], and algebraic structure-based models [17][18][19].
Recently, data-driven modeling has been attracted as a promising methodology owing to its high prediction accuracy for turbulent flows where the linear eddy-viscosity models have significant shortcomings. For example, secondary flows in a square duct and/or separated flows over a smooth surface were successfully predicted via machine [20][21][22][23] or deep [24] learning technologies. For existing data-driven RANS closures, they are also tied to the inherent limitations of modeling in two respects: data dependence still seems to be unavoidable (see, e.g., Figure 2 in Ref. [25]) and validation assessment is inevitably limited to a few benchmark cases. As can be seen above, existing RANS closure models suffer in common from the uncertainty in predictive capability. This long-standing problem may be resolved by constructing mathematically exact relations between Reynolds stresses and then incorporating those relations into the RANS models. To the best of our knowledge, the construction of Reynolds stress identities has not been reported to date. In the present paper, we construct three identities between components of Reynolds stress tensor via spatial mapping. Further, we give a set of equations that is supposed to be applied to dynamically determine the RANS model parameters.

Conceptual Basis
We describe a conceptual basis for the construction of a generic mathematical framework that can provide the identities of physical quantities in a geometrical sense. Let us consider two fundamental physical problems: one is Hagen-Poiseuille flow at a sufficiently low Reynolds number under a favorable pressure gradient, and the other is two-dimensional steady heat conduction in a sufficiently long solid cylinder subjected to uniform internal heat generation and isothermal boundary. The distribution of axial velocity for the former and that of temperature for the latter are visualizable as a three-dimensional (3D) surface at an axial station. These two surfaces are geometrically analyzable by means of a geometric object; both surfaces are parabolic in the centerline, and this is characterized as a positive value of Gaussian curvature. From this estimation two facts are noticed: the geometric approach transcends the boundary between the completely different realms of physics and the geometric object is able to play as a common indicator of estimation without being constrained by that boundary. Extending this geometry point of view to all physical quantities, particularly in R 3 , their spatial distributions may be viewed as surfaces. Further, these surfaces may be characterizable in terms of their geometric objects. From the geometric universality reflected on this shift in perspective, we see a possibility that a mathematical relation on surface geometry may provide a generic mathematical framework to describe spatial distributions of those physical quantities.

Mathematical Framework
To embody the above-described possibility, it is the key to set a geometrical circumstance from the field of a physical quantity. Let T(t, x) be a time-dependent, physical scalar field in the (x, y, z) Cartesian coordinate system and time t and D ⊂ R 3 be the domain for T. To convert the unit of T into that of length makes geometric analysis possible on the converted field. The fundamental form of this unit conversion is written as follows: where σ is the unit conversion function (note that construction of this function depends on the unit of T and the associated physics). Hereafter, let the overtilde sign (.) denote such a converted field. Then, the next step is to extract surfaces from T, and these will be subject to application of the differential version of the Gauss-Bonnet formula below. For this regard, the field of T is subdivided in a stereological manner; three orthogonal coordinate planes are defined at a single point of D at an instant in time, and the respective three graphs on these planes are then extracted from T. We refer to such a graph as a coordinate plane subfield. By extending the subdivision to the entire region of D, three graphs are extracted at each point of D. For example, three coordinate plane subfields at x p (x p , y p , z p ) ∈ D and t = t o are written as T(t o , x p , y, z), T(t o , x, y p , z), and T(t o , x, y, z p ). The differential version of the Gauss-Bonnet formula, which has been proven in Ref. [26], is taken as a mathematical tool to describe coordinate plane subfields of T, since it allows for differential analysis on surface geometry. Theorem 1 therein states that if an orientable smooth surface S is parameterized such that r : where K is the Gaussian curvature over S, N is the normal to S, F a and F b are the products of the geodesic curvatures of the coordinate curves v = const and u = const and those speeds, respectively, and φ is the angle of intersection from v = const to u = const on S.
Hereafter, let Γ, Λ, and Φ denote the three budgets of Equation (2) in order, respectively: In the Cartesian frame, the respective graphs of two-variable function on the three orthogonal coordinate planes are considered: (x, y, f (x, y)), (x, g(z, x), z), and (h(y, z), y, z). In particular for a surface given by parametrization r(u, v) = (u, v, f (u, v)), the geometric entities of Equation (2) are expanded in terms of the derivatives of f (u, v) [27], as follows: where n is the unit normal to S, and the subscripts u, v, uu, vv, and uv denote the first-and second-order derivatives with respect to u and v. For the other two graphs, Equation (2) is similarly expanded. In applying the geometrical relation of Equation (2) to each of the above-considered three graphs and then expanding the geometric entities of each resultant equation in terms of the derivatives of the two-variable function as introduced above, the differential Gauss-Bonnet formula ramifies into three equations: where the subscripts f , g, and h denote that the differential operators correspond to the above-considered graphs, respectively. Adding up the above equations yields a combined one by numerical equality: where the subscript s denotes the sum of the respective operators. Back to the three coordinate plane subfields of T defined at a point of D, the differential Gauss-Bonnet formula is concerned to those fields at their intersection point, respectively. If the differential Gauss-Bonnet formula holds for all coordinate plane subfields of T at each point of D and every instant in time, the combined differential relation of Equation (10) naturally holds for T. Accordingly, it can play a role as an integrated, generic mathematical framework for pointwise analysis on the whole field of T.

Examples
As an example of the above-introduced physics-to-geometry transformation, first in the regime of laminar flow, we consider a fully developed laminar square duct flow. The analytical solution of its streamwise velocity is given [28]: where µ is the molecular viscosity, h is half of the duct side length, dp/dx is pressure gradient in the streamwise direction. The unit of u(y, z) is converted by means of σ = µl c /h 2 (−dp/dx) where l c is a length scale. We choose l c = 1 for simplicity, and consider a subset of R 2 , , as a domain for u(y, z). This is uniformly discretized with (N y , N z ) = (101, 101), where N y and N z are the numbers of grid point in yand z-directions, respectively. Two approximate solutions of u(y, z) are obtained by computing 150 and 200 terms of its infinite series and then adding up those terms, respectively. The approximate relative error of the former in bulk velocity is small enough as 8.684 × 10 −7 that the validity of truncation with 150 terms is demonstrated. From the shift in perspective from physics to geometry, u(y, z) is analyzed on the surface S L given by parametrization ( u(y, z), y, z), (y, z) ∈ D L , in a geometrical sense. For this laminar flow case, only the streamwise velocity component is non-trivial and only its coordinate plane subfield on the yz-plane is not homogeneous. This suggests that the differential relation of Equation (10) is reduced to Equation (9) for S L . Thus, all terms of Equation (9) were computed numerically with second-order finite-difference method. Figure 1a,b display S L and the distributions of the three budgets of Equation (9) at y = 0, respectively. The deviation of the sum of Γ, Λ, and Φ from zero can play a role as an indicator of accuracy. More precisely, it can be shown by this indicator that depending on whether the differential relation of Equation (9) is numerically well resolved, the absolute values of the residual of Equation (9) are less than 9.278 × 10 −5 , as observed in Figure 1c. Since as can be seen in Figure 1a S L is parabolic as a whole, its Gaussian curvature is supposed to be positive in the whole range of D L . Taking into account that sign of Γ entirely depends on that of K, this is identified from the profile of Γ h . On the other hand, the angle of intersection has a non-uniform distribution on S L since the associated parametrization is non-orthogonal; Φ h reaches the maximum at z = 0, then monotonically decreases from the peak location towards both ends of the considered interval, and changes sign from positive to negative at z = ±0.576. To give another example of the physics-to-geometry transformation, but now in the regime of turbulent flow, we deal with the mean axial velocity for fully developed turbulent pipe flow at Re D = 5300 in the (r, θ, z) cylindrical coordinates. For this task, we use the direct numerical simulation (DNS) data of Wu and Moin (2008) [29]. As a preliminary setup, the dimensional mean axial velocity U z is converted into a quantity of length: U z = U z R/u τ , where R is pipe radius and u τ is the friction velocity. As the distribution of U z is axisymmetric about the pipe centerline, it is convenient to deal with this task in the cylindrical coordinate system (note that there is no constraint for a choice of coordinate system in describing the differential version of the Gauss-Bonnet formula).
The spatial distribution of U z is represented as a surface of revolution in the cylindrical coordinate system. Setting the coordinate x of the Cartesian frame to coincide with the axial one, the surface is given by parametrization ( U z , r cos θ, r sin θ) with parameters (r, θ) on a finite domain D T = [0, 1] × [0, 2π). Let S T denote this parameterized surface. Rather than using Equation (10), Equation (2) is directly applied to S T , which makes it possible to obtain single profiles of Γ, Λ, and Φ as a function of r. Figure 2a displays the distributions of the three budgets of Equation (2) as a function of (1 − r) + , which is the wall-normal distance normalized by the viscous wall unit. As can be seen in Figure 2b, the absolute values of residual of Equation (2) are less than 8 × 10 −4 at 0 < (1 − r) < 0.9, which suggests that the differential relation is numerically well resolved at least up to the log law region. Otherwise, the residual sharply increases and then drops steeply near the pipe centerline, which qualitatively agrees with the deviation of the DNS statistics of Wu and Moin (2008) at Re D = 44,000 from the Reynolds-averaged mean axial momentum transport equation (see Figure 6 in Ref. [29]). They attribute the relatively large deviation to the imposed centerline boundary condition. Even taking into account the difference in Re D , the sharp deviation of the residual in Figure 2b might be due to the centerline boundary condition. Accordingly, the following analyzes are focused on the viscous sublayer, buffer layer, and log law region classified for this turbulent pipe flow [29,30], except for the region with relatively high residual values. In the viscous sublayer, the mean axial velocity scaled by u τ is linearly proportional to (1 − r) + , so that S T has a cone-like shape at (1 − r) + ≤ 5. Geometrically, this implies that the Gaussian curvature of S T in this layer must be zero. Taking into account that magnitude of the normal to an orientable surface cannot be zero by definition [31], this geometrical feature is identified in the profile of Γ. In the buffer layer, 5 ≤ (1 − r) + ≤ 30, the parabolicity of S T is identified from the positive sign of Γ, and Γ has a local maximum value of 0.5151 at (1 − r) + = 32.4. In the log law region, 30 ≤ (1 − r) + and (1 − r) ≤ 0.9, and Γ monotonically decreases from the first peak location toward (1 − r) + = 100, and then steeply increases till the end of this layer. Taking into account that U z is axisymmetric about the pipe centerline, the parametrization ( U z , r cos θ, r sin θ) is orthogonal. This means that φ is uniform as π/2 on S T and actually, the maximum absolute value of Φ is negligible as 1.5 × 10 −14 . The distributions of Γ and Λ are hence symmetric to each other about the horizontal line of value zero in the considered whole range of (1 − r) + .

Mathematical Construction
On the basis of the geometrical relation of Equation (10), we now construct identities between components of the Reynolds stress tensor R ij . In this construction process, the key is how to set a quantity of length from R ij . Taking the unit of R ij into account, the unit conversion has a fundamental form: where t c , l c , and u c are time, length, and velocity scales, respectively. Then, the next step is to derive scalar quantities from this converted tensor. The principal invariants of R ij afford three scalar quantities, which are written as follows: Note that these have different units due to the associated mathematical operations. Thus, we return their units back to that of length by scaling: By applying the differential Gauss-Bonnet formula to the coordinate plane subfields of I * i , where R s is the sum of Γ s , Λ s , and Φ s . Depending on a choice of the unit conversion function, these three equations can have different forms; they may consist only of components of R ij or contain other averaged flow quantities. For the latter, the identities may involve other turbulence quantities (e.g., the dissipation rate of turbulent kinetic energy) as well as Reynolds stresses. The applicability of these identities depends on whether coordinate plane subfields of I * i , i = 1, 2, 3, satisfy all the conditions stated in the differential Gauss-Bonnet theorem. If applicable, their mathematical exactness is not subjected to verification because the mathematically proven, numerical equality of the differential Gauss-Bonnet formula ensures it deterministically.

Example
We show an example of the proposed identities through the second-order statistics of velocity fluctuations for the turbulent pipe flow considered earlier. Preliminarily, all components of the Reynolds stress tensor are converted into length quantities with u c = 1 and t c = 1. The variations of the scaled principal invariants of R ij are displayed as a function of (1 − r) + in Figure 3a-c. The distribution of I * i on the yz plane is viewed as a surface of revolution since I * i is axisymmetric. This surface is parameterized as (I * i , r cos θ, r sin θ) on the same domain D T . The three budgets of Equation (2) were computed numerically for the three parameterized surfaces, and their distributions are displayed in Figure 3d-f. Prior to detailed analysis, it has been checked that the residuals of Equation (2) were not significantly deviated from zero at 0 ≤ (1 − r) + ≤ 180, which are plotted in the subfigures of Figure  The geometric properties of the above-described surfaces are analyzed by means of Γ. Three geometrical similarities are found from these surfaces: local concavity in the immediate vicinity of the wall, convexity, and monotonous concavity in order from the wall toward the pipe centerline. Accordingly, one can know from this sequence of appearances that Γ undergoes two sign changes within the considered interval of (1 − r) + . The first concavities very close to the wall are too subtle to be clearly identified, so these are highlighted in Figure 3a-c. However, these concavities are confirmed apparently from the profiles of Γ in Figure 3d-f; the signs of Γ are clearly negative from the wall to the end of the viscous sublayer. While all of I * i increases further toward each maximum after the appearance of the first concavities, their variations are now convex within the buffer layer. Accordingly, the profiles of I * i have inflection points close to the beginning of the buffer layer and, correspondingly, Γ undergoes the first sign changes. After the maximum, all of I * i monotonously decreases till the pipe centerline in a concave tendency. The profiles of I * i thus have other inflection points within the buffer layer and, correspondingly, Γ undergoes the second sign changes. Common to all of I * i , these are subtle since the degrees of those concavities are monotonously weak in a wide range between the beginning of the buffer layer and the end of the log law region.

Relations between the Model Parameters and Reynolds Stresses
In order to dynamically determine the parameters within a RANS closure model, as many equations as those parameters are necessary. We suggest a constraint condition to construct those equations from the three identities constructed earlier. Suppose that the Reynolds stress tensor is modeled with a total of n model parameters, which is written as a schematic form, as follows: where v i is an array of independent variables and c i is a set of model parameters. In this form, all of c i is assumed to be variable in space and time. By repeating the same process to obtain I * i from R ij , the three scaled, principal invariants of τ ij , denoted by J * i , are also obtained. The differential Gauss-Bonnet formula yields three differential relations for coordinate plane subfields of J * i , i = 1, 2, 3: In general, the three budgets of Equations (18)- (20) for I * i are not equal to those of Equations (22)-(24) for J * i , respectively. The errors between them, which arise from modeling R ij , are written as follows: The sum of the squares of these nine errors is written as follows: The constraint condition to minimize S r with respect to each of c i yields a set of equations, whose number is the same as that of the model parameters, as follows: ∂S r ∂c i = 0, i = 1, 2, . . . , n.
If the unit conversion function for R ij and τ ij is set as a constant or a function of known variables, additional unknowns are not associated in R ij and τ ij . Then, the equations given by Equation (29) represent an implicit relation between the undetermined model parameters and Reynolds stresses. Since Equation (29) is highly nonlinear and implicit, it is computationally burdensome to numerically solve the simultaneous equations of Equations (21) and (29) for a benchmark turbulent flow problem. Accordingly, it is clarified that numerical validation for the constraint condition is beyond the scope of the present study.

Summary and Concluding Remarks
In summary, the present study has established a generic methodology to describe the spatial distributions of physical quantities in a geometrical sense. Further, the methodology has been embodied as a mathematical framework on the basis of the differential version of the Gauss-Bonnet formula. This framework yields three identities between components of the Reynolds stress tensor. In the aim to relieve the uncertainties arising from the parameters within a RANS closure model, a constraint condition that makes it possible to construct equations as numerous as the model parameters has been suggested.
Taking into account that subgrid-scale stresses can be converted into a quantity of length, the entire process to obtain Equation (29) could be similarly applied to dynamically determine the model parameters within a subgrid-scale (SGS) stress model for large-eddy simulation. It is open to an empirical study whether the parameters within a RANS or SGS model can be dynamically determined in practice by numerically solving Equations (21) and (29) and further, they are able to enhance the predictive capability of the model.