The Blasius equation dates back to the early twentieth century and describes the boundary layer caused by a uniform flow over a flat plate. Due to the asymptotic condition on the first derivative of the solution, the problem is not an initial value one. The non-linear nature of the equation has led many authors to exploit approximation strategies to obtain an estimate of the solution. Several papers deal with series expansions, but the solutions are partly numerical and do not converge rapidly everywhere. In the present paper, a strategy is proposed to obtain an approximate solution of the Blasius problem in the closed form, i.e., without the use of series expansions. Motivated by the asymptotic behavior of the second derivative, a parametric definition of the third derivative of the Blasius solution is introduced and successive integrations are performed. By imposing the expected asymptotic conditions, a final approximate solution is obtained. This latter has the notable advantage of possessing the same behavior at 0 and at $\u221e$ as the analytical solution up to the third derivative. It will also be shown how this result can be extended to the Falkner–Skan problem with non-zero pressure gradient and provide an accurate estimation of the thermal boundary layer over a flat plate under isothermal and adiabatic conditions.

## I. INTRODUCTION

The modern definition of boundary layer goes back to Prandtl^{1} in his first pioneering work at the beginning of the 20th century. The transition from an asymptotic external (non-viscous) to internal (viscous) solution was formulated by Friedrichs,^{2} who presented a singular perturbation technique to obtain a uniformly approximate solution of a simplified problem. Using Prandtl's approach to deriving the internal solution, Blasius^{3} found (in his Ph.D. thesis) an ordinary differential equation through the definition of a self-similar variable. In this way, the mathematical description of the boundary layer over a flat plate due to a uniform free-stream became possible. Some decades later, Falkner and Skan^{4} achieved a significant generalization of the Blasius equation for non-uniform stream, where the sign and the magnitude of the pressure gradient of the external flow were taken into account.

More than the internal structure of the boundary layer, global information may be more useful in many applications. For these reasons, Kármán^{5} and Pohlhausen^{6} introduced the integral equation describing the spatial variation of the displacement thickness and momentum thickness, used some years later by Gruschwitz.^{7} A further generalization of the Blasius equation, for which similar solutions exist, was later discussed by Goldstein^{8} and Mangler.^{9} According to this formulation, three constants have been taken into account and various solutions have been obtained, depending on particular combinations of them, also in the closed form for some arrangements of the parameters.

In the last century, many authors have dealt with the problem of finding the solution to the Blasius problem. Numerically, this task can be accomplished in terms of a stretched self-similar coordinate. In a given 2D domain, the computational cost depends on the grid refinement and the Reynolds number. Due to the nonlinear term, the transformation of the equation into a differential system is generally preferred (see Ref. 10).

From a theoretical point of view, the problem is quite challenging and a closed-form solution of the original equation has never been found. Rather than looking for a complete solution, many authors have preferred to deal with an approximation. Blasius^{3} himself presented a general approximate method for boundary layers with arbitrary velocity distributions of the outer flow. A power series expansion of the flow velocity with respect to a curvilinear coordinate on the contour of a body was considered. Subsequently, a large literature has emerged for the evaluation of the terms of the series. Many of these results are collected in Schlichting^{11} and recalled in Schlichting and Gersten,^{12} where an approximate evaluation of the wall shear stress is found.

More recently, Parlange *et al.*,^{13} starting from the Blasius approximation, have found a series expansion for the numerical computation of the original solution. Conversely, two different iterative procedures, one based on a variational method and one based on a numerical method, have been suggested by He^{14} and Lin,^{15} respectively. Later, He^{16} developed a simple perturbation approach coupled with an iterative technique to obtain accurate solutions.

In Datta,^{17} another perturbation technique has been applied to the Blasius problem, while in Allan and Syam,^{18} the non-homogeneous problem has been tackled through the homotopy analysis method (HAM). The latter exploits a partly numerical and partly analytical procedure, but, for some values of the initial condition, the method provides nonunique solutions. In Yun,^{19} a simple approach was taken to approximate Blasius' solution, based on a sum of matching functions. An iterative algorithm, based on the fixed point method (FPM), was discussed in Xu and Guo,^{20} and a semi-analytic solution was obtained. A contractive map has replaced the original problem with a set of simpler linear differential problems. In Zheng *et al.,*^{21} a closed analytical solution was presented. The authors found a series expansion of the solution with respect to its first derivative, so that the closed attribute refers to the fact that no external calculations are needed. The series expansion is then evaluated through a numerical procedure.

Among the most recent studies, Jafarimoghaddam *et al.*^{22} investigated a generalization of the Blasius problem using a contraction mapping theorem. Subsequently, Jafarimoghaddam *et al.*^{23} used a decomposition method for solving of the Falkner–Skan problem.

A summary of the works dealing with the solution of the Blasius problem over the past 20 years is given in Table I.

Author . | Year . | Approach . |
---|---|---|

Lin^{15} | 1999 | Iterative method |

He^{16} | 2003 | Perturbative |

Yun^{19} | 2010 | Approximate analytic |

Xu and Guo^{20} | 2013 | Fixed point method (FPM) |

Zheng et al.^{21} | 2017 | Power series |

Jafarimoghaddam et al.^{22} | 2021 | Adomian decomposition method (ADM) |

In the present paper, a different approach based on the asymptotic behavior of the second derivative is illustrated. The final solution is obtained in parametric form, where the parameter values are deduced, once and for all, through the asymptotic trend of the solution and its derivatives. The proposed approach provides an accurate solution up to the third derivative and has the same asymptotic behavior of the analytic solution and the same behavior close to the wall. The advantage of the present approach mainly lies in its possible exploitation in all theoretical applications where the structure of the boundary layer is needed, such as, for example, wall models. Moreover, since many recent research papers deal with the power law shear-driven boundary layer,^{24–28} the same approach can be exploited for the approximation of the Falkner–Skan solution, as it will be shown later.

The paper is organized as follows. In Sec. II, the Blasius problem is briefly recalled and some properties of the solution and of its derivatives are discussed; in Sec. III, the solution strategy is described; in Sec. IV, the velocity and vorticity fields are explicitly expressed and the vorticity is compared with an accurate numerical solution; in Sec. V, the same procedure used for obtaining the Blasius solution approximation is applied to the Falkner–Skan problem; finally, in Sec. VI, the application to the study of the thermal boundary layer over a flat plate is discussed. Conclusions are addressed in Sec. VII.

## II. BRIEF RECALL OF THE BLASIUS PROBLEM

*x*axis, $ u = ( u ( x , y ) , v ( x , y ) )$ is the velocity field (assumed a $ C 2$ function at least), $ p \xaf = p / \rho $ is the pressure field divided by the density

*ρ*of the fluid, and

*ν*is its kinematic viscosity.

^{2}the solution is split into an external and an internal one, where the internal one tends toward the external for $ \nu \u2192 0$. The quantities

*u*and $ p \xaf e$ are the velocity and the pressure of the free stream and represent the external solution of the problem. In order to find the internal solution, which is actually the boundary layer, a coordinate transformation is introduced along the normal direction

_{e}*x*.

*b*is a constant. It is worth noting that the Blasius problem is not an initial value problem, due to the asymptotic condition on $ f \u2032$, so its numerical integration is not trivial. The limiting value

*b*is found in the literature with iterative procedures, but can easily be found in many fluid dynamics textbooks, such as, for example, in Ref. 11 or 10.

*v*is not 0 [see relations (4)], as one would expect, but it is related to the constant

*b*as follows:

*b*, also the value of $ f \u2033 ( 0 ) \u2261 a$ may be found in the literature.

*a*and

*b*are here evaluated. Many strategies are available but, in the present work, the procedure described in Ref. 10 is adopted, so as to obtain the values

^{20}and Zheng

*et al.*

^{21}

## III. SOLUTION STRATEGY

*α*,

*β*, and

*γ*are constants to be determined. It is worth noting that these constants must be all positive in order to not alter the correct behavior of the solution. The form (8) of the solution complies with the correct value at wall of $ f \u2034 ( 0 ) \u2261 0$, but does not respect the right behavior close to the wall. In fact, as it has been previously stressed, it should go as $ \eta 2$, although it actually behaves as

*η*, because of the form 8 of the present strategy. Despite it, the approximate solution

*f*and its first derivative $ f \u2032$ possess the correct trends of

*η*and $ \eta 2$, respectively, as shown below.

*η*= 0 [i.e., $ f \u2033 ( 0 ) = a$] is satisfied once the calculated constants

*α*,

*β*, and

*γ*comply with the relation (10). It must be underlined that, given the asymptotic trend of the error function [ $ erfc ( x ) \u223c \u2009 exp \u2009 ( \u2212 x 2 ) / x$, Ref. 29], $ f \u2033$ behaves as $ exp \u2009 ( \u2212 \eta 2 / \beta + \gamma \eta )$ for $ \eta \u2192 \u221e$ in agreement with the expected decay of the solution, discussed above.

*γ*can be easily expressed as

*η*, Eq. (12) leads to the following estimation:

*a*and

*b*are known

*a priori*(see Sec. II), the coefficients

*α*,

*β*, and

*γ*are found from the relations (10), (13), and (16) by means of an iterative algorithm. The calculated constants are

It is important to underline that in the present approach, the solution and its derivatives are given in the closed form, thus meaning (see Ref. 30): a solution expressed using a finite number of standard operations, which may contain constants, variables, and classical functions, without limit, differentiation, or integration.

The asymptotic behavior of $ f \u2033$, given in Eq. (7), indicates that the value of *β* is 14.17% different from the estimated value of 4, whereas the value of *γ* is 55.84% greater than $ b / 2$. These differences are obviously related to the modeled definition of $ f \u2034$.

In order to quantify the errors on *f* and its derivatives, accurate numerical solutions, obtained with a Runge–Kutta fourth-order method, are used for comparison. The integration step adopted for the numerical solution guarantees a truncation error at machine precision (i.e., $ 10 \u2212 16$). Despite the differences highlighted before, the approximate solution *f* and its derivatives show a remarkable superposition with the numerical output, as sketched in Figs. 1–4.

In terms of the global behavior, the calculated quantities are in good agreement with their numerical counterparts, with a maximum deviation (calculated as the modulus of the percentage error) of 0.017% for *f* (see Fig. 1) and of 0.011% for $ f \u2032$ (see Fig. 2). The average errors in the range $ \eta \u2208 [ 0 \u2212 10 ]$ are 0.005% for *f* and 0.003% for $ f \u2032$. For the derivatives, the deviation is referred to as the largest value (in modulus) of the functions, because their decays for $ \eta \u2192 0$ or for $ \eta \u2192 \u221e$ would overestimate the percentage error in some ranges, leading to the incorrect conclusion of a poorly accurate solution.

Figure 3 shows a satisfying agreement between the second derivative of the approximate solution, defined in Eq. (11), and the numerical output, with the largest error lower than $ a \xb7 0.0004 \u2248 1.32 \xd7 10 \u2212 4$ [ $ f \u2033 ( 0 ) = a$ is the maximum of $ f \u2033$].

The last comparison is performed between the approximate and numerical third derivatives and is sketched in Fig. 4. The agreement is still satisfactory, although for $ \eta \u2192 0$ the numerical solution behaves as $ \u2212 0.25 a 2 \eta 2$ (see the expansion of *f* around 0 discussed in Sec. II), whereas the corresponding approximate function decays as $ \u2212 \eta / \alpha $ [see Eq. (8)]. The maximum deviation observed is lower than the 0.2% with four evident maxima, so it may be concluded that the found solution is a very accurate approximation of the theoretical one.

In terms of the fulfillment of the Blasius equation, Fig. 5 shows the value of quantity $ f \u2034 + 0.5 \u2009 f f \u2033 = \epsilon $ calculated with the approximate solution. The value of *ε* is zero for $ \eta > 8$, whereas for $ 0 < \eta < 8$, the error oscillates between ±0.0002. The largest $ | \epsilon |$ is attained for $ \eta \u2248 1.8$, where a value around 0.0002 is found. Modeling the third derivative is the main source of error, with respect to errors on *f* and $ f \u2033$, so it should also be the dominant term. In order to clarify this point, the errors on $ f \u2034$ and on $ f f \u2033$ are also drawn separately in the same figure. As expected, the error on the term $ 0.5 \u2009 f f \u2033$ is significantly lower: in particular for $ \eta \u2208 [ 0 \u2212 2 ]$, the global error *ε* is almost superimposed on the $ f \u2034$ error.

For $ 2 < \eta < 6$, the error on the product $ f f \u2033$, lower than 7 × 10^{−5}, is mainly related to the deviation of the approximate $ f \u2033$ [given in Eq. (11)] from its numerical counterpart, as shown in Fig. 3.

It is worth to underline that the present procedure is not intended for the estimation of limiting values ( $ \eta \u2192 0$ and $ \eta \u2192 \u221e$) of *f* and of its derivatives, thus implying that the values of $ f \u2033 ( 0 )$ and of $ \eta \u2212 f ( \eta )$ for $ \eta \u2192 \u221e$ have to be assigned in advance. The values of *a* and *b* are indeed given with high precision to obtain a good agreement between numerical and approximate solutions.

Furthermore, it is interesting to make some comparisons between the present approach and the recent literature. In a recent work of Zheng *et al.*,^{21} the solution is obtained with a series expansion of the function *f* and of the self-similar variable *η* with respect to the auxiliary variable $ T = f \u2032$. This approach depends on a single parameter previously calculated (instead of two parameters *a* and *b*, as in the present approach) and the expression of $ f ( \eta )$ is not explicit, as the authors underlined. Respect to the solution of Zheng *et al.*,^{21} the present approach describes an accurate approximation given in a single compact expression, so that it is more practical and versatile to handle. Moreover, the solution of Zheng *et al.*^{21} requires about 200.000 terms to give reliable values of *f* and $ f \u2032$.

In order to make a quantitative comparison with other approaches of the recent literature, the iterative algorithm of Xu and Guo^{20} is taken into account. The convergence of the algorithm was reported in a table where 50, 200, and 700 iterations were used to collect some values of the solution and of its first and second derivative. The values of *f*, $ f \u2032$, and $ f \u2033$ reported by the authors for *n* = 700 (in Table II)^{20} are used for comparison on top frames of Fig. 6: the differences with the present method are lower of $ 1.6 \xd7 10 \u2212 3$ for *f*, lower than $ 1.05 \xd7 10 \u2212 3$ for $ f \u2032$, and lower than $ 1.2 \xd7 10 \u2212 3$ for $ f \u2033$ in the range $ \eta \u2208 ] 0 \u2212 6 ]$. For *η* = 0 the differences are 0, while for *η* higher than 6, they become negligible.

Re . | L_{2} norm
. |
---|---|

1 | $ 3.68 \xd7 10 \u2212 7$ |

10 | $ 5.06 \xd7 10 \u2212 7$ |

100 | $ 1.95 \xd7 10 \u2212 6$ |

1000 | $ 4.89 \xd7 10 \u2212 6$ |

Re . | L_{2} norm
. |
---|---|

1 | $ 3.68 \xd7 10 \u2212 7$ |

10 | $ 5.06 \xd7 10 \u2212 7$ |

100 | $ 1.95 \xd7 10 \u2212 6$ |

1000 | $ 4.89 \xd7 10 \u2212 6$ |

By considering that the accuracy of the Xu and Guo^{20} approach is very similar to the present one, it can be concluded that the present method is comparable in precision with others in the literature, but its compactness makes it easier to handle and to implement. Moreover, as it will be discussed later, the present approach is generalizable to more complex problems as the Falkner–Skan one in a rather straightforward way.

Finally, in order to complete the comparison with other recent solutions found in the literature, the solution (15) is drawn together with the approximations of Lin,^{15} He,^{16} and Yun^{19} on the bottom frames of Fig. 6. The present approach shows the best agreement with the reference numerical solution discussed above.

## IV. VELOCITY AND VORTICITY EXPRESSIONS

In Fig. 7, the comparisons between accurate numerical solutions at different Reynolds numbers and their analytical approximation are portrayed in terms of vorticity contour lines. By considering that a reference length is missing and the Reynolds number depends on the local abscissa, the ratio $ u e / \nu $ has been set (corresponding to *Re _{x}* at

*x*= 1). As clear from Fig. 7, the level of agreement is remarkable and the differences become negligible for higher Reynolds numbers as graphically visible by comparing the vorticity isolines. The

*L*

_{2}norm of the error at different Reynolds numbers is reported in Table II. The error is calculated for every grid point as the difference between numerical and approximate solution. Because the absolute average value of the vorticity increases with the Reynolds number, the error norm increases as well, passing from order $ 10 \u2212 7$ to $ 10 \u2212 6$. However, it remains very low, so that it can be concluded that the approximation found here is very reliable in terms of accuracy.

## V. FALKNER–SKAN GENERALIZATION OF THE BLASIUS PROBLEM

*x*

_{0}the distance from the leading corner where the free stream velocity is known, the boundary layer is described by the Falkner–Skan problem

^{4}

*m*is a constant. By differentiation of $ u e ( x )$, the link between the exponent

*m*and the pressure gradient is easily deduced

*m*= 0 for the Blasius problem. Using the same strategy described in Sec. II, an estimation of the asymptotic behavior of $ F \u2034$ can be easily deduced as

*b*(

*m*) is $ lim \eta \u2192 \u221e ( \eta \u2212 F )$ and, similarly to

*a*and

*b*for the Blasius problem, $ F \u2033 ( 0 ) \u2261 a ( m )$. In this case, the value of $ F \u2034 ( 0 )$ is equal to −

*m*from Eq. (20) and a parametric expression is considered

*β*and

*γ*. Unlike the Blasius problem, the parameters

*α*,

*β*, and

*γ*are here functions of the pressure gradient and must be evaluated for varying

*m*. By recursively integrating the parametric function (22), the following quantities are obtained:

*χ*.

_{m}*α*,

*β*, and

*γ,*which must be evaluated for varying

*m*. In fact, together with the relations (23)–(25), the following constraints must be enforced:

*a*,

*b*,

*α*,

*β*, and

*γ*by

*m*is not made explicit. The quantity $ q m \u2261 m \u2009 ( \beta / 2 \u2212 a \alpha )$ with $ q 0 = 0$ for the Blasius problem.

The asymptotic matching condition $ ( \eta \u2212 F ) \u2192 b ( m )$ leads to the condition (28) that is more complicated than the corresponding formula (17) for *m* = 0. Indeed, the presence of a pressure gradient leads to a third-order algebraic equation in *β* rather than a simple fractional expression. However, the above conditions are not sufficient to obtain the entire set of parameters *α*, *β*, and *γ* because they still depend on *a*(*m*) and *b*(*m*), which, similarly to the Blasius problem, must be known *a priori*.

*a*(

*m*) and

*b*(

*m*) are numerically found by iterative procedures (see also Ref. 10). From Schlichting

^{11}and Riccardi and Durante,

^{10}it is evident that the dependence of $ F \u2033 ( 0 )$ from

*m*is monotonically increasing, as depicted in the left frame of Fig. 8. It should be emphasized that the minimum admissible value of $ m \u223c \u2212 0.090 \u2009 43$ corresponds to $ F \u2033 ( 0 ) = \u2202 \eta u | \eta = 0 = 0$, so the

*τ*is zero and the separation occurs (for more details, the interested reader is addressed to Ref. 12). The flow condition at separation cannot be modeled by the Falkner–Skan problem because the solution is not more a similar one. By considering the monotony of the function

_{w}*a*(

*m*), its inverse

*m*(

*a*) is taken into account and expanded with

*a*up to the second order

*a*= 0.332 057 3 for

*m*= 0 from the Blasius problem (see Sec. III) and

*a*= 0 for $ m = \u2212 0.090 \u2009 43$, as explained above. Moreover, exploiting the Hartree profiles,

^{35}the following expression holds:

*a*(

*m*). By comparing the obtained data with accurate numerical outcomes, it turns out that the quadratic approximation (29) is affected by very small errors (lower than $ 10 \u2212 5$) at least in the range $ ( m ( 0 ) , 1 )$. A comparison between the numerical solution and the analytic expression (29) is depicted in the right frame of Fig. 8.

*b*(

*m*) is still needed. By considering that for $ m \u2260 0$, the following identity holds:

*b*(

*m*) when combined with the expressions (23) and (25).

Once calculated all the parameters, the Falkner–Skan approximate solution is evaluated for varying *m* and compared with accurate numerical solutions. The comparison is shown in Fig. 9.

The superposition of the approximate solution (25) and of its derivatives with the corresponding numerical ones is very encouraging, thus demonstrating the accuracy of the approximation found with the present strategy.

## VI. THERMAL BOUNDARY LAYER IN A UNIFORM STREAM

The approximate solution of the Blasius equation (15) allows to obtain accurate estimations of a wide class of problems. One is the solution of the thermal boundary layer generated by a uniform stream over a flat plate. The thickness of this region, where the advective thermal fluxes are comparable with the diffusive ones, is indicated with *δ _{T}* in analogy with the thickness of the kinetic boundary layer.

*c*is the specific heat of the fluid,

_{p}*T*is the temperature field,

*T*is the plate temperature, and

_{w}*T*,

_{e}*p*, and

_{e}*ρ*are the temperature, the pressure, and the density of the free stream, respectively.

_{e}*f*represents the state equation of the fluid, while $ K / ( \rho c p )$ is the thermal diffusivity and

_{s}*μ*the dynamic viscosity. Dirichlet or homogeneous Neumann boundary conditions are related to isothermal or adiabatic plate, respectively. By neglecting the density variations due to the non-uniformity of the temperature field, the velocity field is given by the similar solution of the Blasius equation. This suggests to look for the temperature distribution as a function of the similar variable

*η*. By defining the non-dimensional temperature,

*T*is the plate temperature in the isothermal case or a reference temperature otherwise. Once repeated the same considerations of Sec. II, the following equation is obtained:

_{w}*f*in Eq. (33) is the Blasius solution for which $ f = \u2212 2 f \u2034 / f \u2033$ holds (see Sec. II). Making this substitution and dividing by $ f \u2033 ( Pr )$ [consider that: $ f \u2033 ( \eta ) > 0$ for $ \eta \u2265 0$], the following expression is obtained:

*C*

_{1}the first integration constant. With a second integration, the final solution is found

*T*<

_{w}*T*then Br is always lower than 0 and

_{e}*h*> 0, so that the plate is always heated by the flow. If

*T*>

_{w}*T*, the plate is heated for Br > 2 and cooled for $ 0 <$ Br < 2.

_{e}## VII. CONCLUSIONS

In the present work, an analytical approach to the solution of the Blasius equation is shown. Starting from the asymptotic behavior of the second derivative of this solution, the third derivative is modeled with three parameters, the determination of which is found, in the closed form, by enforcing the expected asymptotic conditions on the solution and on its first and second derivatives. The approximate solution found has a high accuracy when compared to the numerical results, thus representing a good approximation of the analytical one. Compared with other analytical approaches, the current one offers comparable or superior precision, negligible computation times and no convergence issues.

Two constants *a* and *b*, representing the values of $ f \u2033 ( 0 )$ and $ lim \eta \u2192 \u221e ( \eta \u2212 f ( \eta ) )$, respectively, have to be known *a priori*. The knowledge of these constants does not represent an hindrance, because they can be easily found in many papers and textbooks and are here given once and for all with high accuracy.

In order to show the reliability of such a solution when derived quantities are calculated, the vorticity field at different Reynolds numbers has been compared with an accurate numerical solution, showing a remarkable superposition and negligible errors in *L*_{2} norm.

Furthermore, to underline the direct applicability of the present approach to more complex problems, an application to the Falkner–Skan equation is also addressed, showing again an encouraging overlap to the numerical approximation of the solution for varying pressure gradient. Finally, the solution of the thermal boundary layer over a flat plate was also discussed and the case Pr = 1 was investigated.

In order to assess the reliability of the present approach, a comparison with other analytical and approximate solutions of the recent literature has been provided. The present solution, although approximate, reveals greater accuracy than other approximate approaches.

Numerical applications where knowledge of the velocity field near a wall is important, such as wall models used for computational fluid dynamics (CFD), can benefit from the compactness of the present solution, which offers a remarkable precision in a single formula.

## ACKNOWLEDGMENTS

The research activity has been developed within the Project Area Applied Mathematics of the Department of Engineering, ICT and Technology for Energy and Transport (DIITET) of the Italian National Research Council (CNR).

The present research was not supported by any funding.

## AUTHOR DECLARATIONS

### Conflict of Interest

The author has no conflicts to disclose.

### Author Contributions

**Danilo Durante:** Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Validation (lead); Writing – original draft (lead); Writing – review & editing (lead).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

*Verhandlungen des III Internationalen Mathematiker Kongresses*

*Die turbulente Reibungsschicht in Ebener Strömung Bei Druckabfall Und Druckanstieg*

*Elementi di Fluidodinamica: Un'introduzione per L'Ingegneria*

*Handbook of Mathematical Functions With Formulas, Graphs, and Mathematical Tables*

*Mathematical Proceedings of the Cambridge Philosophical Society*

*Fluid Mechanics*