A Closed-Form Analytical Solution to the Complete Three-Dimensional Unsteady Compressible Navier-Stokes Equation

The three main approaches to exploring fluid dynamics are actual experiments, numerical simulations, and theoretical solutions. Numerical simulations and theoretical solutions are based on the continuity equation and Navier-Stokes equations (NSE) that govern experimental observations of fluid dynamics. Theoretical solutions can offer huge advantages over numerical solutions and experiments in the understanding of fluid flows and design. These advantages are in terms of cost and time consumption. However, theoretical solutions have been limited by the prized NSE problem that seeks a physically consistent solution than what classical potential theory (CPT) offers. Therefore, the current author embarked on a doctoral research on the refinement of CPT. He introduced the Refined Potential Theory (RPT) that provides the Kwasu function as a physically consistent solution to the NSE problem. The Kwasu function is a viscous scalar potential function that captures known and observable unsteady features of experimentally observed wall bounded flows including flow separation, wake formation, vortex shedding, compressibility effects, turbulence and Reynolds-number-dependence. It is appropriately defined to combine the properties of a three-dimensional potential function to satisfy the inertia terms of the NSE and the features of a stream function to satisfy the continuity equation, the viscous vorticity equation and the viscous terms of the NSE. RPT has been verified and validated against experimental and numerical results of incompressible unsteady sub-critical


Introduction
At the 2018 American Institute of Aeronautics and Astronautics (AIAA) Aerospace Sciences Meeting in Florida, United States (US), the current author presented a work in progress on the analytical solution to the complete Navier-Stokes equations (NSE) [1]. This formed the basis for the current author's doctoral research at the Daniel Guggenheim School of Aerospace Engineering, Georgia Institute of Technology, Atlanta, US that was completed in 2020 [2]. The present publication is a first in a series to communicate the research findings to the community. It describes an engineer's viewpoint on the existence and smoothness of the NSE solution(s) [3]. This viewpoint is expressed by Anderson[4,p. 185] when he writes that "we have to obtain a solution of the governing equations...-a solution that will give us the results for the dependent flow field variables p [pressure], ρ [density], V [velocity], etc., as a function of the independent variables of spatial location and time [t]" of the flow from which practical quantities of engineering interest can be obtained [4].
Traditionally, the steady inviscid incompressible flow over a body can be obtained by solving the steady Euler equations with the appropriate boundary conditions or by the superposition of elementary flows with known solutions using potential theory [2,5]. This is in recognition of the fact that the velocity field of such a flow obeys the Laplace equation which is linear and lends itself to analytical solutions by the method of separation of variables. The applicability of the Laplace equation on incompressible flows has experimental evidence, and its solution for some bodies and configurations share some qualitative characteristics with experimentally observed flows as exemplified by the streamlines of the unseparated flow over a circular cylinder at the freestream Reynolds number, Re ∞ = 1.54 (based on the freestream velocity, V ∞ , the cylinder diameter, D, and the fluid kinematic viscosity, ν, as Re ∞ = V ∞ D/ν) presented in Fig. 1. In this figure, the flows are from left to right. This was, hitherto, the state of the art in the closed-form mathematical solution of the cylinder flow. To make a distinction between this and the theory presented herein, it will be referred to with the adjective 'classical'. At first, the qualitative comparison between Figs. 1a and 1b seems accurate. However, a quantitative comparison reveals some differences. Due to the assumption of an ideal fluid in the classical theory, it would be expected that its prediction should only approach reality in the inviscid limit as the Re becomes so high that the effects of the viscous terms vanish from the equations of motion of an experimentally observed flow contrary to the figures. In addition to this mismatch in Re limits, the classical theory always realizes a slip velocity at the wall and assumes the flow to be irrotational in the entire domain, in contrast to the mechanism which gives rise to the observed presence of the boundary layer and vorticity in actual viscous flows. Finally, the gap between the prediction of zero drag for a steady flow in the classical theory and actual experiments is dubbed d'Alembert's Paradox [2,6].
Classical potential theory (CPT) falls short of reconciling the actions of viscosity and frictional forces in an experimentally observed flow with the theoretical analyses of such a flow. As such, it is unable to resolve the boundary layer and predict the especially important flow separation phenomenon that results in the pressure drag experienced by a body in the flow [2]. Actual flows (a) The streamlines from CPT at the inviscid limit, Re∞ → ∞ for a cylinder flow [2].
(b) The streamlines for the experimentally observed creeping flow over a cylinder [4] (Re∞ = 1.54) Fig. 1: A comparison of streamlines from theory and experiment are viscous flows whose physics are governed by the NSE [2][3][4]7]. An attempt to refine CPT of a finite circular cylinder flow to capture important features of a bluff body flow including flow separation and reattachment, wake formation, vortex shedding, turbulence, compressibilty effects and Re-dependence is synonymous with solving the NSE problem [2].
Adequate solutions of the NSE should appropriately predict flow separation and result in a finite drag on a body within the flow. However, the NSE are highly complex and non-linear. Up until the current author's thesis, there was no known physically acceptable analytical solution of the full NSE for arbitrary bodies [2]. At the turn of the 21 st century, the Clay Mathematics Institute of Cambridge Massachussets (CMI) included finding a solution or the proof of a solution or the proof of no solutions to these equations as one of its seven millennium prize problems [3]. CMI has been at the forefront of the coordination of the worldwide search for the general solution(s) (or proof of no solution(s)) to these relations. A check on the status of the problem on its website (accessed on November 28, 2021) still reveals it to be unsolved despite attempts by various researchers [2,8,[8][9][10][11][12][13][14][15][16][17][18][19][20][21]. However, with limiting but suitable assumptions and approximations, exact and approximate analytical solutions of these relations have been obtained for some special flows. Additionally, direct numerical solution (DNS) of these flow field equations within a discretized domain can be obtained with the use of high-speed computers. Depending on the complexity of the flow and the computational speed, such simulations can take weeks of wall time [2].
Efforts at solving the NSE eventually led to the formulation in the early 1900s of Boundary Layer Theory (BLT) which simplifies the NSE into the boundary layer equations (BLE) [2,4,[22][23][24]. In his proposition of BLT, Prandtl recognised that the viscous effects in a wall-bounded flow are largely confined to the thin boundary layer just adjacent to the wall. Prandtl's premise on BLT is the experimentally backed assumption of the no-slip condition at the wall. He surmised that the no-slip condition, imposed on the flow by friction and viscosity, ensures that the fluid adheres to and also assumes the velocity of the wall. Since the flow velocity field is continuous, the fluid rapidly accelerates to the local freestream velocity at the edge of the boundary layer where inviscid theory reasonably holds [2,22]. The boundary layer hypothesis has proven crucial in understanding and predicting viscous phenomena including flow separation, skin friction drag and aerodynamic heating [2]. With Prandtl's submission, the flow field around bodies have been obtained with inviscid theory applied outside the boundary layer coupled with the solution of the BLE. This has proven successful in over a century of applications, and it has been generally considered as the resolution of the d'Alembert's Paradox. Nonetheless, there are conflicting discussions on the validity of the no-slip boundary condition in the literature [25,26] as exemplified by the work of Hoffman et al. [12,27]. They computed drag with a slip boundary condition from the computational solution of the unsteady Euler equations and obtained results comparable to actual experiments [2].
Various attempts have been made at applying the proper orthogonal decomposition (POD) to the circular cylinder flow. Particular interests have been taken in understanding and characterizing the wake dynamics as well as 'training' the resulting model in a Galerkin projection [2]. This is to capture the flow dynamics at different Re than the one used for the POD modes extraction [2,28]. This seems promising; however, Deane et al. [28] observed that the 'trained' model is highly sensitive to the Re used to extract the POD modes for the cylinder wake flow. This is because the mean flow used for the Galerkin projection is itself Re-dependent. Therefore, the Re-dependence of the mean flow must be extracted from a series of experimental and/or numerical data which will serve as inputs to the POD process. This further amplifies the volume of data, cost and time needed to obtain a sufficiently accurate reduced order model [2].
The gap in resolving these cylinder flow dynamics is necessary to improve the advances in the theoretical prediction of fluid dynamic phenomena. Many approaches have been tried, but they have been limited in the dimensions of the flow and the regime of the Re. Therefore, the current author proposes the refined potential theory (RPT) that bridges the gap between CPT and experimentally observed canonical bluff body flows [2]. RPT introduces a viscous potential function, the Kwasu function, that simultaneously satisfies the continuity equation and the NSE [2]. Section 2 of this publication reviews the literature. Section 3 discusses how the Kwasu function conserves both mass and momentum and satisfies the flow boundary conditions. Section 4 discusses some results. Section 5 concludes the article.

Literature Review on Extant Solutions and Models
The set of the NSE for the compressible flow of a Newtonian fluid in an arbitrary control volume, V, are expressed in single vector equations, using the indicial notation as [2] V ∂ρ ∂t where g is a body force vector, µ is the dynamic viscosity, and λ is the bulk viscosity. A Newtonian fluid is one for which there exists a linear relationship between stress and strain [2,24]. Equation 1 is a statement of the conservation of mass, and Eqn. 2 embodies the conservation of momentum for fluid flows [2][3][4]24]. The indicial notation in Eqn. 2 is necessary to express the stress tensor [24]. The general form of this tensor including a normal stress component, thermodynamic pressure (p), and the tangential stress (shear stress,τ ) is given as for a Newtonian viscous fluid [24]. When an analogy between Eqn. 3 and the strain relation in elasticity is drawn, the summation of Eqn. 3 is employed to define the mechanical pressure,p, as identifying the pressure as a tensor invariant [24, p. 67]. Thus, the term λ + 2 3 µ ∇ · V in Eqn. 4 accounts for the difference between the thermodynamic pressure in an undisturbed freestream and the mechanical pressure due to a flow in question [23,24,29,30]. It is important for compressible flows and flows in which the viscous normal stresses are not negligible. However, when the flow is incompressible (∇ · V = 0) as per Fefferman [3], there is no distintion betweenp and p, and λ vanishes from the NSE [2,24].
It is usually taken that the integrands in Eqns. 1 and 2 are identically equal to zero at every point in the flow since the choice of the control volume must be arbitrary [4]. Then, the issue with the NSE, formulated by CMI [3], is to obtain a realistic closed-form solution to simultaneously satisfying the continuity equation subject to the appropriate boundary conditions with an initial divergence-free velocity field and an assumption of incompressibility. A solution is accepted as physically reasonable if it defines the velocity, V, and the pressure, p, as smooth functions of the spatio-temporal dimensions and has bounded energy for t ≥ 0 [2,3]. Equation 5 is rewritten as to highlight the role of the vorticity vector, ω, in the convective acceleration term through the application of some vector identity [2,24]. g, is also assumed to be the gradient of a scalar force potential, g [2].
Fefferman reports that some of the proposed solutions to these relations in two-dimensions are unstable in the time domain and offer no clear path to solutions in three-dimensions. He makes no mention of whether an accepted solution would be that of a free flow or a wall-bounded flow, and if wallbounded, no indications as to the acceptable boundary condition at the wall. This provides some leeway to researchers to work with the most appropriate formulations. He concludes that the problem persists when using known standard mathematical methods for solving partial differential equations and calls for new ideas [2,3].
For two-dimensional "incompressible" flows, − → ψ , on the one hand, satisfies the continuity equation identically. Its physical realization has been the basis for its well-documented usage in the literature in the vorticity-stream function formulation for the solution of flow fields around objects (Eqn. 12) [2,24,[34][35][36][37]. On the other hand, − → ψ cannot satisfy the NSE by identity in due to the cross-derivatives and its rotationality condition (ω = 0) [2].
A number of different computational approaches to solving Eqn. 11 exists based on the variables used. Either the solution is sought in terms of the primitive variables of V and p; or in terms of the ψ and ω; or in terms of ψ alone [2,38].
In terms ψ and planar vorticity, ω z , the governing NSE for two-dimensional constant property flow is presented as which is known as the vorticity equation. The cylindrical polar coordinate provides a convenient coordinate system to examine the circular cylinder flow [2,24,[34][35][36][37]. Numerical solutions of Eqn. 12 can be obtained using computational fluid dynamics (CFD) methods. ψ and ω z can also be expanded as truncated Fourier series, inserted in Eqn. 12, and numerically solved in a semi-analytical method [2,35]. An accepted general analytical solution of the vorticity equation does not exist even though it is one dimension less complicated than the full NSE. However, an exact analytical solution of this equation exists for a viscous line vortex. CPT also provides an analytical solution for the trivial case when ω z equals zero [2]. Additionally, Talaei and Garrett [16] propose a general analytical solution to this equation for an axially symmetric two-dimensional flow around a sphere assuming that ψ is subject to a separation of variables [2].
The disagreement of CPT with actual experiments is generally attributed to the loss of the no-slip condition at the wall [2,7]. When the fluid is assumed to have no viscosity, the only wall condition imposed on Eqn. 12 is the impermeability of the normal velocity [2,7]. However, Hoffman et al. [12] employed the slip condition at the wall to avoid the difficulty of numerically resolving boundary layers in high Re flows. They report values for aircraft lift and drag comparable to actual experiments. They suggest that the integrated properties of the flow are pressure-based, rather than viscous driven [2,12].
Vorticity is introduced into the flow and convected downstream due to the very high velocity gradient between the wall and the local freestream [2,4,24,27,33]. The vorticity effect of the boundary layer is idealized in theoretical solutions with inviscid vortex sheets [2,4,33] as exemplified by thin airfoil theory [2,4]. This is also the approach proposed by Helmholtz, and later, Levi-Civita, in his contributions to the theoretical prediction of finite drag on a non-accelerating moving body in an inviscid fluid -the wake hypothesis [2,39]. Whilst being successful in the prediction of lift on the airfoil, it is less so with drag [4]. Traditionally, the non-lifting flow over bodies are simulated with source/sink sheets, and the lifting flow with vortex or doublet sheets through panel methods. The combination of sources, sinks and vortices are sometimes used in these methods [2,4]. In RPT, the boundary layer flow is idealized as one consisting of contributions from local vortices and sources/sinks that are mutually concentric at every location on the wall [2]. Each of these elementary flows is an exact and smooth solution of the NSE [2].
Hoffmann et al. weigh in on Fefferman's analyses as they conclude against the existence of a "well-posed smooth solution to the [NSE] with smooth data." [12,27] Buckmaster and Vicol [13] are in agreement with this conclusion. They buttress this point of view with their findings and proofs that weak solutions to the NSE are non-unique [13,21]. Ramm [18][19][20] also maintains that the NSE problem is not "physically acceptable" [18]. He writes that "the NSP [Navier-Stokes problem] is not a correct description of the fluid mechanics problem and the NSP does not have a solution." [19] However, Chen and Kratka "establish the existence of global solutions to the initial boundary value problem for the [multidimensional] NSE for compressible heat-conducting flow with large spherically symmetric initial data between the solid core and the free boundary connected to the surrounding vacuum state." [40] They also show that, no cavities develop between the solid core and the free boundary that expands with some finite speed [2,40].
Muriel [9,10] used the time evolution of a one-particle distribution function initialized with spatially uniform data along with a choice of a Gaussian pair potential between particles to arrive at a divergence-free velocity field, but he encountered problems in the integration of the NSE to obtain a scalar pressure field. He concludes with a suggestion that the pressure field is a higher-order tensor field [2].
Scholle et al. [11] also encountered a difficulty in obtaining a viscous velocity potential that is self-consistent with the pressure field. Thus, they reformulate the two-dimensional NSE in terms of complex variables. This is to derive a first integral in the form of a real-valued scalar potential and a complexvalued velocity field from which the pressure can be obtained. Furthermore, they obtain a real-valued tensorial representation of the complex-valued first integral of the reformulated NSE. A scalar pressure field is then computed after solving a set of second-order partial differential equations subject to the flow-specific boundary condition(s) [2,11]. They offer an extension of the two-dimensional derivation to the general case of an unsteady threedimensional viscous flow using tensor calculus [2,14]. Their approach results in a tensor-valued and a vector-valued field equations that are constrained by a vector potential introduced for the velocity and involves five unknown fields. These were solved analytically (for a translating disc in a viscous fluid and a non-axisymmetric stagnation flow) and numerically (for a viscous flow in a cubic domain) in a successful verification of their methodology. To achieve these results, some of the equations involving the unknown potentials were eliminated using the guage invariance or symmetry of the NSE [2,14].
Amoloye [1,2] explores a viscous-potential function,κ, that satisfies the physical relationship between the scalar pressure field and the velocity field in as a realistic closed-form solution of the incompressible NSE [1,16]. The existence of such a function is well-known and has been researched by many [1,2,11,14,16,31,32,41]. Although, widely considered an oxymoron, the idea of a viscous potential flow or function is not new [2,32]. Joseph [32] advocates this in his historical review of the potential flow of viscous fluids [2]. The idea can be traced back to both George Stokes [2,32] and Horace Lamb [2,32], and later to Joseph and Wang [2,32]. In addition, a potential-based Lagrangian solution to the NSE has been obtained for some compressible flows for which µ is neglected in favor of λ [2,14]. However, generalizing these formulations to obtain the full solutions of unsteady viscous three-dimensional compressible flows could not be guaranteed to be analytically tractable [2,14]. Stokes' earlier attempt at the theoretical study of the cylinder flow had resulted in the paradox named after him which was resolved by Oseen [2,24]. Neglecting the inertia terms in the NSE, Stokes sought and derived a stream function that solved the steady motion of a sphere in a flow, but not that of a cylinder. Oseen partially incorporated the inertia terms in a proposed set of linearized equations. Lamb solved these linearized equations and found his solutions valid only for small Re [2,31]. The scope of the asymptotic theoretical analyses on the NSE, building on the works of Oseen and Lamb, have been generally limited to the low-Re flows [2,16,42,43]. However, Talaei's and Garrett's [16] proposed analytical solution to the two-dimensional vorticity equation for a moving sphere in a quiescent and viscous fluid extends these analyses to high Re [2].
White [24] writes that an assumption of an ideal zero viscosity fluid is not always necessary to define a potential flow, but when ω is identically zero regardless of the viscosity, V must be the gradient of a potential function [2]. This is the case everywhere in the high Re flow past a solid body except in the boundary layer [2,24]. Therefore, he concludes that CPT describes exact solutions of the full NSE [24, p. 85]. Schlichting [23, pp. 72-73] and Anderson [4, p. 209] agree with White's analyses. Anderson also employs the continuity equation for an incompressible potential flow to establish the applicability of the Laplace equation to such a flow [4, p. 236]. The summation of particular solutions of the Laplace equation is also one of its solutions because the Laplace equation is linear [2,4]. Thus, complex flows of practical importance are derived from the analytical superposition of individual elementary flows that are solutions of the Laplace equation with appropriate boundary conditions [4, pp. 238-239]. Lamb [31] also discusses the Clebsch transformation that extends the definition of a velocity potential to inviscid flows with nonzero ω. In that formulation, ω is identified as an invariant for a barotropic flow [2,11,14,41].
Closed-form analytical solutions have the advantages that their development illuminates the physics behind a problem. They give information on the relationships between the pertinent variables, and they are very useful for rapid and less expensive preliminary design purposes [2,4]. In the absence of comprehensive theoretical solutions, the need to obtain accurate qualitative and quantitative characterizations of the dynamics of fluid flows precipitated the development of reduced order models of high fidelity numerical solutions and/or experiments [2,44]. Some of the methodologies in use to obtain these models include POD and the Dynamic Mode Decomposition (DMD) [2]. The idea behind POD in fluid dynamics is that spatial velocity correlations can be orthogonally decomposed to identify the dominant eigen modes that characterize the coherent structures in the flow [45,46]. Developing a dynamic model from these spatial modes, i.e., computing the temporal coefficients, requires projecting these modes onto the NSE. Alternatively, these coefficients are obtained directly from the input data set by computing the ratio of a single snapshot matrix with the generated modes [2]. In a two-pronged approach, DMD decomposes the flow field into its spatial modes and the corresponding eigenvalues which signifies the dynamic evolution of these modes with time without the need to project the spatial modes onto the governing equation of the system [2,47]. Attempts are usually made to train the reduced order models derived from these data analysis methodologies using the NSE. These are geared towards making them useful for other geometries and Re than the ones from which they were obtained. However, their success is not always guaranteed [2,28,48]. This process is also reminiscent of the search for an analytical solution(s) to the complete NSE [2].
Raissi et al. [49] present the hidden fluid mechanics algorithm which is claimed to be able to encode the NSE from an input data set and provide a trained model that is agnostic to geometry, initial and boundary conditions. However, the success of this and other input-data-dependent methods depends highly on the fidelity of the input data that in some cases are very expensive to collect [2].
These literature findings motivated Amoloye's doctoral research on the Kwasu function, κ, in a refinement of CPT [2]. The Kwasu function is a viscous scalar potential function that captures known and observable unsteady features of experimentally observed wall bounded flows including flow separation, wake formation, vortex shedding, turbulence, compressibility effects and Re-dependence [2]. It is appropriately defined to combine the properties of a three-dimensional potential function to satisfy the inertia terms of the NSE and the features of a stream function to satisfy the continuity equation, the viscous vorticity equation and the viscous terms of the NSE [2].

The Eulerian Kwasu Function: A Viscous
Velocity Potential
This means the velocity field is perpetually divergence-free [2]. The flow is governed by the Laplace equation because φ can be defined such that ψ also identically satisfies the divergence-free condition since However, it only satisfies the Laplace equation when the flow is considered irrotational [2]. An irrotational flow is one in which ω is a null vector [2,4,6]. For a two-dimensional flow, the condition for irrotation is For ψ, this means That is the Laplacian of the stream function and it satisfies the Laplace equation [2]. φ identically satisfies the irrotational condition [2,4,6] since Experimentally observed flows are rotational although some of their regions can be assumed irrotational [2,4,6,24,50]. It can be seen from Eqns. 15-18 that when the flow is rotational, φ is not defined. That is to say it is not physical and cannot be used to adequately describe experimentally observed flows. However, ψ is defined. It is observable with flow visualization techniques that reveal streamlines [2,51]. Streamlines are level sets of ψ. ψ identically satisfies the continuity equation. For an experimentally observed rotational flow, it satisfies the Poisson's equation (∇ 2 ψ = −ω z ) but not the Laplace equation [6]. Both of these differential equations are linear and obey the superposition principle [2,52]. Therefore, the use of a superposed streamfunction to study experimentally observed flows is justified [2].
Separated-variables solutions to the Laplace equation abound in mathematics literature [2,52]. In a cylindrical polar coordinate system, the general solution is obtained as an infinite Fourier series for ψ and φ as where m is a positive integer and θ m is an arbitrary constant angle [2]. ψ and φ are harmonic functions because they satisfy the Laplace equation [2,52]. Depending on the values for m and the arbitrary constants (A 0 , A m , B 0 and B m ) in these equations, these solutions correspond to the inviscid irrotational incompressible lifting flow of a uniform freestream over a number of two-dimensional geometries ranging from the circular cylinder (m = 1) [2,4] to semi-infinite wedges (m > 1) [2,53]. These general solutions are identifiable as superposition of elementary flows [2,4,53]. This was the state of the art in the closed-form mathematical solution of the steady incompressible flow field around infinitely spanned geometries -potential flow theory [2,6]. Any arbitrary definition of the constants in these solutions do not preclude the functional relationship they describe between the dependent and independent variables that makes them viable solutions of the fluid dynamics equations (refer to Eqns. 10-11) [2]. Schlichting writes that "frictionless flows [potential flows] may also be regarded as exact solutions of the Navier-Stokes equations, because in such cases the frictional terms vanish identically." [23, p. 72] However, the gap between the prediction of zero drag therefrom for a steady flow and results from experimentally observed flows is dubbed d'Alembert's Paradox [2,4,6].
Other elementary flows such as sinks, sources and vortices are also present in complex flows like the experimentally observed cylinder flow [2]. These flows are illustrated in Fig.2. The potential solution for an inviscid vortex can be obtained as where Γ is the circulation or strength of the vortex [2,4]. However, experimentally observed flows are viscous [2,24]. For a line viscous vortex, the incompressible NSE reduce to the heat equation for which various analytical solutions exist in the literature [2,37,52]. Amongst such solutions is the Lamb-Oseen vortex model [37, p. 262]. The circulation for this model is given as where r 2 c = 4νt, ω z is the planar vorticity, Γ 0 is the initial circulation, ν is the kinematic viscosity of the fluid and r c is the vortex core radius [2]. This is a time-dependent model of the viscous decay of the vortex [24, p. 207]. When the time component in the equation is fixed, this vortex solution becomes the Burgers' vortex model [2,37]. Burgers' vortices have been used to model eddies in turbulent flows successfully [37, p. 265]. Therefore, a viscous stream function can be defined for a vortex as using Eqn.23 [2]. Similarly, the stream function for a potential sink/source and the radial velocity can be expressed as V r = Λ 2πr (26) where Λ is the strength of the sink/source [2,4].
To highlight the contributions from the elementary uniform, doublet and vortex flows, Eqn. 19 is rewritten as in which R is the body radius, A 0 = −B 0 ln R and R 2m = −B m /A m . It can be seen that B 0 is associated with the strength, Γ, of the vortex flow [2]. If it is a constant value and m = 1, Eqn. 27 corresponds to an inviscid flow over a rotating cylinder that produces lift depicted in Fig. 3   To obtain the components of velocity, ψ is differentiated in the relevant directions to give For the non-lifting flow (NLF) over a circular cylinder in which the uniform flow is parallel to the body x-axis as illustrated in Fig. 4, B 0 = θ m = 0, m = 1, A 1 = V ∞ , B 1 = k/2π and R 2 ≡ k/2πV ∞ in Eqn. 27 and the corresponding expressions for its velocity field [2,4]. That is, the stream function is the velocity components are and the speed is In Fig. 4, κ stands for the doublet strength [4]. This is different from the predominant use of κ for the Kwasu function in this publication. The freestream is assumed to be steady, uniform and unperturbed by the cylinder. That is, at r = ∞, This is consistent with experimentally observed flows [2,4]. At the cylinder surface where r = R, Eqns. 29 and 30 yield ψ N LF = 0 identifies the cylinder surface as a streamline of the flow. V r N LF = 0 predicts that there is no flow through the surface since it is impermeable [2,4,6]. Both of these are in agreement with an experimentally observed flow [2,4]. However, the prediction of a slip tangential velocity, V θ N LF = 0, on this non-rotating wall violates the no-slip condition that has been observed in experiments [2,4,51]. The no-slip condition gives rise to the boundary layer [2]. A boundary layer is that region of an experimentally observed flow adjacent to the surface where viscous effects including flow rotation are dominant [2,4,24,50]. ψ N LF is inviscid or non-viscous and features no boundary layer [2,4].
Due to wall friction, an experimentally observed local freestream velocity is retarded to the surface velocity. This results in a vertical velocity profile illustrated in Fig. 5 and a crossflow velocity gradient [2,4,24,50]. Boundary layers are usually small in comparison to the body, and the crossflow velocity gradient occurs over a relatively short distance [2,22,23,31]. There is also a streamwise gradient of velocity [2].

Boundary Conditions
It is assumed that the cylinder is a sole body of gravity subjected to a flow in an infinite domain and referenced to an inertial spatio-temporal frame [2]. With a further assumption that an experimentally observed flow over a nonrotating cylinder is a subset of the general solution of the Laplace equation (Eqn. 27), the viscous stream function should have the form [2] ψ viscous (r, θ, V ∞ , ν, t, Re) = V ∞ r sin θ 1 − R 2 r 2 +B 0 (r, θ, V ∞ , ν, t, Re) ln and its velocity should be It should also satisfy the freestream boundary condition and the wall no-slip condition [2]. B 0 ln (R) enforces the no-slip condition at the wall where r = R as and (38) Equation 38 shows that B 0 is at least f (r, θ, V ∞ ), and its form should be based on the tangential velocity of the NLF to satisfy the no-slip condition. B 0 ln (R) can be understood as the reaction (based on Newton's third law) at the wall to the introduction of a vortex into the flow [2].
The limit of Eqn. 34 as r → ∞ is The freestream uniform flow should be recovered for Eqn. 39 to be bounded. However, the natural logarithm term precludes this except if B 0 → 0 at a faster rate in the limit as ln (r/R) → ∞. This suggests that B 0 should consist of an exponent of the negative of r that tends to zero as r → ∞ in a viscous and time-dependent way [2]. When this is the case, the vorticity, is also bounded and confined to the boundary layer/wake. It is observable in Eqn. 40 that the inviscid term in ψ viscous does not contribute to ω z because it is irrotational [2]. Thus, based on the preceding discussion and the viscous stream function for a vortex described with the Lamb-Oseen model in Eqn.24, a basic functional (a) A sink-source-vortex sheet on the cylinder surface [2] (b) The streamlines of the superposition of a vortex plus a sink (with arrowheads), and vortex plus a source (without arrowheads).
However, when an experimentally observed flow impinges a surface, the local change in momentum of the fluid elements is propagated radially outward from the surface [4, p. 604]. This is assumed to be similar to local source/sink flows situated at the surface. Therefore, to theoretically simulate the boundary layer effect in the NLF, the cylinder surface is replaced with a viscous sink-sourcevortex sheet as depicted in Fig. 6a. The sheet has the sinks/sources and the vortices positioned at the same locations. The resulting flow pattern of this superposition at a specific location is illustrated in Fig. 6b. The strengths of the sinks/sources and vortices are dictated by the velocity field of the NLF. Additionally, their locations are defined by the dividing streamline equation for the NLF. These locations are continuous and not discrete. The dividing streamline, ψ = 0 in Fig. 4 delineates the external flow from the internal flow and ensures the impermeability of the surface. The strengths of a sink /source and vortex at a specific location on the viscous sheet are not the same but vary in a manner that reflects the physics of the experimentally observed flow [2]. These considerations are employed to modify B 0 further from Eqn. 41

Conservation of Momentum
The integral compressible NSE, Eqn. 2, are expressed in an equivalent differential form in A solution to the incompressible NSE for the pressure field is sought in terms of a known stream function/velocity components in However, a number of issues immediately arise. One is that the stream function ψ viscous [2, p. 90] is two-dimensional whereas experimentally observed cylinder flows are three-dimensional [43, p. 246]. Stream functions are usually difficult to define for general three-dimensional space [2,4,37]. Additionally, integrating the pressure out of Eqn. 42 by vector identity manipulations is impossible due to the terms involving the ω and the local temporal derivative [2,11,14]. Although, ψ viscous [2, p. 90] is continuously differentiable in space and time, present mathematical methods with the use of a symbolic mathematical software like Maple™ do not offer anti-derivatives for such an involved equation [2]. Numerical integration would not be consistent with solving the problem analytically. However, when that is done, different results are obtained from integrating in either of the two spatial directions. Muriel [9] encountered this problem, and he suggests that the pressure is a higher-order tensor which converges to a scalar value in a long-time limit. However, the current literature on pressure appears to suggest otherwise [2,4,24]. Sheng proposes a review of the fundamental physical bases of the NSE with which he envisions easier solutions of NSE [1,15]. However, his proposition paves no path towards resolving ω [2,15]. A possible option is to carry out the foregoing theoretical development in terms of φ which is well defined for three-dimensional flows and would satisfy the unsteady Bernoulli equation (Eqn. 10). It cannot be overemphasized, however, that φ is mathematically undefined in its present form for a rotational flow [2]. For viscous fluid flows, the condition of irrotationality must be dropped [2,33]. This essentially negates the physical existence of φ. However, the presence of pressure as a scalar quantity in the NSE suggests a potential function for a rotational velocity field exists [2,14,31]. That function is a stream function and a velocity potential simultaneously. This is the central idea to the Kwasu function [1,2]. Obtaining the incompressible Kwasu function,κ from requires a direct integration process that is not supported by available modern mathematical methods [1,2]. However, this is accomplished by resolving ψ viscous [2, p. 90] in the Cartesian coordinate system as follow [1,2]. The cylindrical polar coordinate variables are related to the Cartesian x-y axis as This is illustrated in Fig. 7. There are at least 360 otherx-ỹ axes that can be defined in relation to the same cylindrical polar coordinates in Fig. 7. These correspond to successive units of θ. So, the angle in ψ viscous [2, p. 90] is redefined as θ = arctan (y, x) ⇒ arctan (x, y) (44) in which the negative windx-axis is coincident with the positive body y-axis. Thus, the functionκ [2, p. 96] that satisfies the continuity equation and the incompressible NSE is summarized in Table 1. Table 1 also comparesκ with ψ and φ. This is illustrated in Figs. 8 -11 that present plots of the level sets of the three functions and their velocity fields for the inviscid non-lifting cylinder flow [2]. The freestream flow parallel to the negative windx-axis that is coincident with the negative body x-axis in Figs. 8a and 8b now approaches from the positive body y-axis in Fig. 8c. The streamlines and the potential lines in Figs. 8a and 8b respectively are orthogonal. However, they describe the same velocity field [2,4]. Whereas ψ is differentiated in a perpendicular direction to a flow to obtain the velocity components in Fig. 9, φ is differentiated in the flow direction to obtain the same velocity components in Fig. 10. This is the reason φ results in an irrotational flow, and ψ does not [2,4].
In the polar coordinate system, the velocity components are obtained from κ in a similar way to ψ as depicted in Figs.11a and 11b. However in the (c)κ = r sin (arctan (r cos θ, r sin θ)) 1 − R 2 r 2 Fig. 8: A visual comparison of the classical stream function, ψ, Kwasu function,κ and the classical velocity potential, φ for the cylinder NLF [2].
Cartesian coordinate system,κ is defined on a principal axis for a viscous rotational flow about which ω is identically zero [23, pp. 57-58]. This axis is almost aligned with the wind axis with its negative x-axis pointing into the incident flow [2].κ is a potential function in this axis, and its velocity components are obtained as    These are illustrated in Figs.11c and 11d for the z = 0 plane. Positive z-axis extends out of the page. Figure 12 shows that the radial velocity vanishes at the surface because the surface is impermeable [2].
The ω components in the principal axis are computed as and ω z = 0 is a solution to the two-dimensional vorticity equation However, the flow remains rotational especially towards the surface and in the wake [2,23].
In the polar coordinate system,κ identically satisfies the two-dimensional instantaneous continuity equation, and it yields the component of ω for the flow as [2] and The pressure, p, which is a scalar quantity (zeroth-order tensor) is invariant to the foregoing rotation of the Cartesian axes [2,23]. The Laplacian ofκ required to obtain the viscous term of the pressure field is computed as [2] and (61) The foregoingκ is a two-dimensional potential stream function for a crossflow in the center plane of a finite cylinder. To incorporate the cylinder span,κ axis is transformed to ensure the flow's rotational symmetry in the wind axis. This is accomplished in Eqn. 44 with the substitution [2] Fig. 13: Present theoretical streamlines of an unseparated crossflow over a finite cylinder [2]. .
where θ 3D = arctan (y 3D + y m , x 3D ) In Eqn. 63, y m is a geometric constant (Fig. 14), and B is the cylinder span. The positive directions of the geometric variables are indicated with the arrowheads in Fig. 14. All the angular coordinates are in radians. R 2D and R 3D are the body radial coordinate, R, in two-and three-dimensions respectively [2]. Equation 63 defines the two-dimensionalκ in all ϕ = constant planes of a three-dimensional domain x 3D = r cos θ sin ϕ y 3D = r sin θ sin ϕ z 3D = r cos ϕ .
For all of these planes,κ is simultaneously a viscous potential function and a stream function following the properties highlighted in Table 1 and Eqn. 69. This is also illustrated in Fig. 13 in which the flow is from left to right [2]. y m controls the flow axis about the cylinder as illustrated by the red arrowheads in Fig. 15. In these figures, the flow is also from left to right. When y m is zero in Fig.15a, the direction of the maximum gradient of the streamlines highlighted by the arrowhead is vertically inwards and directly above the cylinder crest. This is contrary to an experimentally observed flow depicted in the lower image that has this direction inclined to the cylinder crest. When  [4,54]) at Re = 1.54 [2]. Fig.15b, there is a consistence with the experimentally observed flow pattern [2]. In summary, and by inspection of the rotated circumferential coordinate, θ W in Eqn. 62, the wind axis is defined as [2] ϕ = π 2x = r cosθ sinφỹ = r sinθ sinφz = r cosφ = 0 Fefferman allows an assumption of no body forces (g = 0) or body forces defined by potential functions (g = ∇g) [2,3]. With this assumption, the pressure field that identically satisfies the incompressible NSE can be obtained from The relationship between ψ, φ andκ is established as follows. For a twodimensional flow, lines of constant ψ are streamlines [2]. Along a streamline, Equivalently, lines of constantκ are equipotential streamlines, and surfaces of constantκ are equipotential streamsurfaces (or equipotential streamtubes) [2]. On such a surface, dκ = ∂κ ∂r dr + ∂κ r sin ϕ∂θ r sin ϕdθ + ∂κ ∂r ∂r r∂ϕ rdϕ (69) To ensure that the equipotential streamsurfaces remain physical, no flow must be allowed to cross them [2]. This means that the cross-product of the surface vector and the velocity vector must be equal to zero as Additionally, the continuity equation must be satisfied at all times as κ is not an axisymmetric solution in the body axis since [2] V ϕ = ∂κ ∂r ∂r ∂z The derivation of its functional form from CPT for a viscous flow is elaborated on in Amoloye's thesis [2, pp. 46-98].

Results and Discussion
RPT and the Kwasu function are discussed extensively in Amoloye's thesis [2]. RPT revisits Stokes' hypothesis and provides a theoretical determination of the bulk viscosity for high Re bluff body flows [2, pp. 121-124]. It resolves d'Alembert's Paradox [2, pp. 140-150] . It is used to explain the sign changes in the shear stress distribution at the rear of a circular cylinder in a crossflow [2, pp. 128-131]. It is verified against experimental and numerical data on a cylinder in an incompressible crossflow at Re ∞ = 3, 900 [2, pp. 151-170]. Its drag prediction is within the error bound of measured data and CFD prediction [2, p. 158]. The predictions of the pressure distribution, separation point and Strouhal number are also within acceptable ranges [2, pp. 157-158]. Its prediction of the force coefficients over the range 25 ≤ Re ∞ < 300, 000 is validated against experimental and theoretical data on a cylinder in crossflow [2, pp. 158-160]. There is a good agreement in the magnitude and trend for Re ∞ > 100 [2]. For Re ∞ < 100, there is a disparity in magnitude that is unsafe for design purposes [2]. Similarly, it under-predicts the coefficient of drag in some of the explored cylinder axial flow configurations [2, pp. 170-172]. However, at Re ∞ = 96, 000 and an aspect ratio of 2 for a cylinder axial flow, the RPT drag prediction falls within 1.2% of validated computational result [2, p. 171].
The drag coefficient of a sphere for 25 ≤ Re ∞ < 300, 000 is explored to demonstrate the application of refined potential theory [2, p. 176]. Additionally, the flow over a stationary sphere at Re ∞ = 100, 000 is explored in detail [2, pp. 174-175]. A generally good agreement is observed in the prediction of the experimental trend for Re ∞ ≥ 2, 000 [2, p. 176]. The transitional incompressible flows over a 6 : 1 prolate spheroid at 45 • angle of attack for Re ∞ = 3, 000 and Re ∞ = 4, 000 are also explored [2, pp. 177-180]. The RPT pressure distribution has a close agreement with the DNS result in the starboard rear of the spheroid [2, p. 178]. However, the magnitude of the predicted force coefficients are generally less than five times the corresponding DNS results [2, p. 180]. The asymmetry of the DNS pressure distribution in the meridian plane is not captured [2, p. 178]. Therefore, more in-depth analyses of the spheroid flow including the separation locations are subjects of further studies [2, p. 186].

Conclusion
The three main approaches to exploring fluid dynamics are actual experiments, numerical simulations, and theoretical solutions [2]. Numerical simulations and theoretical solutions are based on the continuity equation and NSE that govern experimental observations of fluid dynamics. Theoretical solutions can offer huge advantages over numerical solutions and experiments in the understanding of fluid flows and design. These advantages are in terms of cost and time consumption [2]. However, theoretical solutions have been limited by the prized NSE problem [3] that seeks a physically consistent solution than what CPT offers.
In CPT, the steady inviscid incompressible flow over a body can be obtained by the superposition of elementary flows that are exact solutions to the NSE. However, CPT falls short of reconciling the actions of viscosity in an experimentally observed flow with the theoretical analyses of such a flow. As such, it is unable to resolve the boundary layer and predict the especially important flow separation phenomenon that results in the pressure drag experienced by a body in the flow. This relegated the use of potential theory to idealized flows of limited practical importance [2]. Therefore, the current author embarked on a doctoral research on the refinement of CPT [2]. He introduced RPT that provides the Kwasu function as a physically consistent solution to the NSE problem. The Kwasu function is a viscous scalar potential function that captures known and observable unsteady features of experimentally observed wall bounded flows including flow separation, wake formation, vortex shedding, compressibility effects, turbulence and Re-dependence [2]. It is appropriately defined to combine the properties of a three-dimensional potential function to satisfy the inertia terms of the NSE and the features of a stream function to satisfy the continuity equation, the viscous vorticity equation and the viscous terms of the NSE [2].
RPT has been verified and validated against experimental and numerical results of incompressible unsteady sub-critical Re flows on stationary finite circular cylinder, sphere and spheroid. It is concluded that the Kwasu function is a physically consistent and closed-form analytical solution to the NSE problem.