Heat pipes (HP’s) are very interesting heat transfer devices, currently widely used in almost all cases of heat-stressed cooling of spacecraft and satellites, military facilities and household devices and computers. Together with the wide practical application, the HP’s are of great interest for the scientific study of moist vapour currents, the occurrence of pulsations in the vapour channel during the intermittent vaporization operation of the capillary-porous evaporator, the formation of toroidal vortices near the cooled condensation surface and many other interesting effects   . The perform of HP’s vapour channel profiled in the Laval-liked nozzle form allows to increase the heat transfer coefficient at high heat power loads, when boiling begins and a large amount of vapour is formed in the capillary-porous evaporator   . The numerical analysis of the vapour flows in the profiled vapour channel shows the occurrence of toroidal vortex formations near the cooled HP’s surface  . These vapour vortices are of great interest because their toroidal rotation direction depends on the thermal power entering the HP’s evaporator and affects the thickness of the condensate film of the working fluid. The earlier measurements of the condensate film thickness  were compared with the calculated values. The original HP’s developed for all condensate film thickness measurements are presented below in Figure 1.
Figure 1. HP’s diagram: 1―top cover; 2―cylinder body of the HP’s; 3―cone-shaped turbulator; 4―capillary-porous insert; 5―bottom cover; 6―capillary injector channels, 7―bottom flat capillary-porous insert-evaporator. There are capacitive sensors 8, 9, installed inside the top cover     .
2. Experimental Results
In this study was applied stainless steel HP’s with shaped vapour channel in the converging-diverging nozzle form of nonlinear varying diameter and with flat top and bottom covers   . In the upper cover capacitive sensors are installed to measure the working fluid condensate film thickness, as which the diethyl ether was used with a boiling point at atmospheric pressure 308.65 K (35.5˚C).
All details of the HP’s manufacture, profiled vapour channel formation along the entire HP’s length and hydraulic diethyl ether delivery system to the evaporator, design and calibration of capacitance measuring sensors, development of high-frequency electronic generators and measuring devices are given in the .
The dimensions are chosen in such a way that it is possible to establish capacitive sensors for measuring film thickness. Experimental results on the diethyl ether condensate film thickness measurements depending on the heat power value, entering to the HP’s evaporator are shown in Figure 7   .
3. Numerical Modeling Results
The present paper represents the numerical modeling of the vortex flows inside a profiled vapour channel of short HP’s, finite volume and pressure based CFD 10.0 code, for predicting compressible vapour flow with friction, heat transfer and converging-diverging nozzle area change. Navier-Stokes and heat-transfer equations with measured boundary conditions were solved, i.e. using fixed temperature values of heat source and heat outlet. A flat resistive heater is used as a located at the HP’s bottom heat source, and a original design flow calorimeter was used as a heat drain   , the measurement error of calorimeter thermal power does not exceed 1.7%
The model was studied as a longitudinal section along the axes of the two injector channels, which helps to preserve all the specific features of whirling instability under the conditions of continuous circulation motion of the working fluid during liquid and vapour phases. CFD 10.0 code employs a finite volume formulation of mass, momentum and energy conservation equations in a network consisting of nodes and branches. The mesh discretizes the nozzle space into small tetrahedral solid elements second order of accuracy and high quality, with ten nodes and enough accurate in the problem modeling. The forming vortex flow in the HP’s upper part, cooled by the flow calorimeter, are shown on Figure 2 and Figure 3.
The cross-section profile of the vapour toroidal vortex ring near the cooled by the flow calorimeter upper HP’s cover 1, Figure 1 is asymmetric and non-circular, what is primarily due to friction on the condensate film surface of the diethyl ether.
On the moving in a radial direction on the cover to the capillary-porous insert 4, Figure 1 due to capillary forces, the diethyl ether film is additionally affected by the vapour vortex ring, rotating in the toroidal direction and shown in Figure 2 and Figure 3.
Figure 2. The occurrence of a toroidal vortex ring near the condensation surface inside of Laval-liked HP’s at low heat loads. Moving vapour jets due to the impact of the Coanda effect sticks to the HP’s walls and vortex ring formed from the periphery to the center of the vapour channel.
Figure 3. The occurrence of a toroidal vortex ring near the condensation surface inside of Laval-liked HP’s at higher values heat loads entering to the HP’s evaporator. In this case, the vortex ring formed from the center to the periphery of the vapour channel.
If the toroidal rotation direction of the vapour vortex due to the Coanda effect occurs from the periphery to the longitudinal axis of the HP’s vapour channel as shown in Figure 2, the vortex due to the surface friction slow down the diethyl ether film moving to the periphery, and its thickness δf in the near-axis area on the HP’s upper cover will have inflated values.
With increasing the heat power value, the direction of toroidal rotation of the vapour vortex acquires the opposite direction, from the longitudinal axis to the periphery of the vapour channel, as it shown in Figure 3, therefore, due to the surface friction сf of the film diethyl ether on the cover will accelerate and its thickness δf will have to decrease.
Figure 4 shows the cross section of the toroidal vapour vortex ring at the overheat of the HP’s evaporator relative to the boiling point of the diethyl ether at atmospheric pressure δT = T − TB = 15 K, and the temperature difference between the evaporator and the condensation surface Tev − Tcond = 25 K.
Giving HP’s vapour channel profiled Laval-liked nozzle form at a large heat power value and the boil beginning in the capillary-porous evaporator and the formation of large vapour amount leads to the toroidal vortices formation near the cooled condensation surface.
The vortices occurrence is a decisive factor that changes the longitudinal components of vapour velocity and pressure distribution in the Laval-liked nozzle vapour channel, and explains the increase in the heat transfer coefficient KHP, in comparison with standard HP’s with a cylindrical vapour channel and a simple vapour return currents in them at equal overall dimensions.
The calculated confirmation of the heat transfer coefficient KHP increase for the HP’s with a Laval nozzle-liked vapour channel is shown in Figure 5. The
Figure 4. The computational cross section profile of the toroidal vapour vortex near the cooled in the flow calorimeter HP’s upper cover, the velocity is in m/s. Toroidal vortex has a large contact surface with the underlying liquid condensate film of diethyl ether and due to the surface friction cf affects on the radial velocity of the film flow.
Figure 5. The results of numerical analysis of the distribution of the vapour flow longitudinal velocity component inside the Laval-like HP’s vapour channel of compressible moist vapour at high heat loads.
upper part of the Figure 5 represents the calculated vapour velocity component along the longitudinal axis z and shows the occurrence of two peaks―one about 85 m/s in the critical section of the vapour nozzle and the second one peak with a velocity value of 30 - 35 m/s in the right corner near the HP’s condensation surface.
In the center of the toroidal vortex ring there is a rarefaction (−δPv), the value of which is about |−δPv| ~ 2 × 103 Pa, or 0.5% of the condensation vapour pressure 4.1 × 105 Pa in the cooled part of the HP’s vapour channel.
The lower part of the Figure 5 shows the velocity distribution in the vapour channel also along the longitudinal axis z in the critical section of the Laval nozzle and the peripheral part of the toroidal vapour ring with a negative current velocity value (counter flow) near the channel walls.
The counter current (negative value) velocity on the toroidal vortex outer surface forms a negative dynamic pressure about―2 × 103 Pa, and the total pressure drop on the vortex ring near the cooled HP’s condensation surface reaches 4 × 103 Pa, or about 1% of the condensing vapour pressure. The vapour pressure gradient inside the toroidal vortex with the vapour channel radius 1×10−2 m is no less than 4 × 105 Pa/m. This value is more than enough to increase the condensation intensity of the diethyl ether along the vortex constant pressure gradient line (surface).
4. Analytical Results
Let’s consider the interaction of a vapour jet with a HP’s flat cover. The origin of the coordinate system is placed at the junction central axis of the vapour channel with a HP’s flat cover, the z coordinate is measured in the opposite direction to the vapour jet direction.
The gas dynamics equations in the Cartesian coordinate system represent a description of a unsteady viscous compressible vapour flow interaction using the following equations:
Equation (1) is supplemented by the ideal gas equation of state:
Vector of conservative variables Q and vectors of vapour flows F, G, H have the following form:
The components of the viscous stress tensor and the components of the heat flux vector are written in a standard way:
In the above equations, t―time, s; ρvp―the moist vapour density, kg/m3; u, v, w, m/s, components of the vapour flow velocity in the coordinate directions x, y, z in the HP’s channel; p―pressure, Pa; e―total energy per unit mass of moist vapour, kJ/kg; Tvp―vapour temperature, K; γ―the specific heats ratio; λeff―effective thermal conductivity, W/m・K; μeff―effective viscosity is calculated in the standard way as the sum of molecular μ and turbulent μt viscosity:
where ―the increase in the mass of the condensed film in a unit of time, kg/s; Jvp―the condensing vapour flow, kg/s.
In the vapour flow braking area on the HP’s upper cover near the central axis (r = 0), the speed of the vapour jets significantly reduced, in the turning and vortex formation section are strongly curved. The depth (thickness) parameter of vortex is b.
In the near-wall section (vortex area formation), the shear stresses τs become dominant compared to normal ones τn. Spreading in the radial (and tangential) directions vapour jet loses momentum (amount of motion), the thickness of the boundary layer and the condensate film thickness δf increase and the static pressure increases also. The calculation results show the occurrence of a large vortex formation formed in the mixing layer of the vapour jet in confined conditions of the vapour channel with HP’s flat cover surface. The vapour pressure pulsations Reynolds number  in the vapour channel nozzle throat reaches Revp ~ (1.5 - 1.65) × 105, the Prandtl number Pr = 0.77.
The structure of the vortex formation on the HP’ upper cover have a defining effect on the heat exchange in this region of the vapour jets interaction with a flat end cap (flat wall). The heat transfer characteristics are non-uniform even at small values of the Reynolds number Revp ~ 5 × 102. The Nusselt number Nu has a characteristic maximum at the breaking point (r = 0), and the minimum value in the upper cover area near the walls where the inverse (counter) vortex flow is formed .
The justification for the heat exchange characteristics change in the vapour jets interaction point (r = 0) with the HP’s flat top cover and the occurrence of a Nusselt number local maximum are explained by the formation of a large-scale toroidal vortex  and a laminar-turbulent transition in the boundary layer and an turbulence kinetic energy increase in the wall jets  .
To describe the flows arising from the interaction of a vapour jets with HP’s flat top cover, the Navier-Stokes equations (averaged over the Reynolds number) are usually used, in solving which the shortcomings of various turbulence models are apparent manifested. In the event of vapour toroidal vortices, an additional source of turbulence energy generation of the condensed vapour occurs due to the large curvature of the current lines near the HP’s condensation surface.
To assess the effect of the moving vapour jets curvature in the toroidal vapour vortices occurrence and a more complete assessment of the turbulent energy generation, the so-called Kato-Loner correction  is used, which is an additional source of turbulent energy in the toroidal vortex vapour flow forming, that can be written in both integral and differential form . This increases the stability of the iterative process, although the time-averaged characteristics of the flows were studied.
Figure 3 and Figure 4 shows a simplified diagram of a condensable vapour flow, is included in the inner hole of a vapour toroidal ring vortex of radius rv, with a temperature Tv and axial component velocity uv and tangential component velocity vv. We consider the circumferential (tangential) velocity component in the Laval-liked nozzle vapour channel is not equal to be zero. Static vapour pressure at the inlet of toroidal vapour ring vortex is equal to Pv, the cylindrical coordinate system is used.
As a main model variant, to simplify the analysis we consider the stationary state flow of a diethyl ether film condensate between the stationary state vapour toroidal ring vortex and the flat surface of the HP’s top cover. The vapour toroidal vortex ring and the condensate film move in the radial and tangential directions. The task is to determine the thickness δf of the diethyl ether condensate film under the toroidal ring vortex.
The proposed mathematical model uses the principle of the Karman and Pohlhausen   , the basic idea of which is the solvability of the differential equations system on average thickness of the boundary layer δ, as which the thickness of the condensate film δf is considered. The peculiarity of the problem under consideration is that the longitudinal and tangential components of the velocity in vapour channel are commensurate.
Therefore, the momentum theorem is applied twice―for radial and circumferential directions of the vapour flow:
where uf, vf ―radial and tangential condensate film flow velocity under the vapour vortex, m/s; uv, vv―radial and tangential vapour velocity in the vapour vortex, m/s; δf ―the hydrodynamic thickness of the diethyl ether condensate film under the vapour vortex ring, m, in general depends on the vapour flow rate; ―the condensate film thickness of flow momentum loss due to the friction in the radial direction, m; ―the condensate film displacement thickness, taking into account the displacement of the current lines in the condensate film due to the viscosity of working fluid, m; ―the condensate film thickness of flow momentum loss due to the friction in the tangential direction, m; τrs = μ(duv/dz)z = δf―radial shear friction stress on the outer surface of condensate film under the vapour vortex ring, Pa.
Momentum conservation equation for tangential (circumferential) direction of the condensate film flow look like this:
where τφs = μ(dvv/dz)z = δf―tangential shear friction stress on the outer surface of the condensate film under the vortex ring, Pa; ―the condensate film thickness of flow momentum loss due to the friction in the radial and tangential directions, m.
The thicknesses of the liquid condensate film δ* and δ** are determined in the standard way   :
where ρl―the density of liquid condensate, kg/m3.
Equations (7) and (8) contains more than two unknowns, therefore additional physically based linking equations are needed. The velocity vector deviates from the geometric radial direction, and the degree of deviation in the flux core and near the surface of the film are different.
Tangential (circumferential) vapour velocity v in the flux core in the in the shaped vapour channel we can define in the simplest way:
where m is a numerical parameter of unity order; rv―the inner radius of the toroidal vortex, m.
The relationship between the pressure, density and temperature of the condensing vapour in the first approximation can be given by the ideal gas equation of state in the following form:
where R = 8.3144 J/mol・K―universal gas constant; k = 1.31―value of the adiabatic index of diethyl ether vapour; Tv―the temperature inside the toroidal vapour vortex, K; ρv―the vapour density inside the toroidal vortex, kg/m3; Pv―the vapour pressure inside the toroidal vortex, Pa.
After integrating the differential Equation (7) for the vapour core flow in the Laval nozzle shaped vapour channel along the radius r and taking into account the equation of state (10), we obtain an equation to estimate the pressure distribution over the radius r in the core flow:
From the continuity equation of the vapour flow in the integral form and subject to the vortex size we obtain the next correlation:
where b―the thickness of the vapour toroidal ring vortex, m. Figure 6 shows the vapor density values of the of diethyl ether in the HP’s with a profiled vapour channel and a cylindrical vapour channel.
Figure 6. The calculated diethyl ether vapour data density ρvp in the HP’s vapour channels along the central axis z at the temperature of the evaporator Tev = 323.65 ± 0.01 K (50.5˚C ± 0.01˚C). The temperature of the HP’s condensing surface in the flow vortex calorimeter is 295 ± 0.03 K (22˚C ± 0.03˚C). 1―solid line, the vapour density in the Laval-liked vapour channel at a pressure P = 6.13 × 105 Pa; 2―dotted line, the vapour density in the HP’s with cylindrical vapour channel.
From this relation, taking into account the density ρvp from (7) and (8), an equation can be obtained to estimate the radial velocity in the core of the vapor flow at the inlet to the vapour toroidal ring vortex:
In general, the vapour flow in the Laval nozzle shaped vapour channel is swirling flow, which is due to both the uneven jet character of operation of capillary-porous evaporator with injection channels and the variable area of the HP’s passage section. For swirled vapour flows, the spin parameter is the determining in the HP’s vapour channel parameter, characterizing the ratio between the circumferential (tangential) and longitudinal (axial) velocity components inside the channel.
To determine the local twist of the passing into the liquid phase (condensing) vapour stream through the thickness of the diethyl ether film condensate (boundary layer), the parameter of the local spin angle is used  :
where tgφ―parameter of local spin angle.
Within the variable thickness of the of liquid condensate film under vapour toroidal ring vortex parameter of local twist tgφ depends on the radius on the flat surface of the HP’s condensation top cover and the linear parameter of the film thickness z. At the vapour stream core in the diffuser part of the vapour channel tgφ0 parameter characterizes the degree of deviation of the vapour flow from the radial direction:
where tgφ0―parameter, defining the degree of deviation in the diffuser part of the HP’s vapour channel from the radial direction.
The liquid condensate film, formed under the vapour toroidal ring vortex, due to surface interfacial friction is fond of it and moves in the radial and tangential directions on the HP’s smooth condensation flat (wall) surface, moreover, the wall shear stress τw ~102 Pa  on the top cover surface (condensation surface) is considered the same for both directions and is taken into account in all calculations.
The parameter tgφmax represents the limiting (surface) tangent of the twist angle of the formed liquid condensate film and is the ratio of surface shear stresses of friction in tangential (circumferential) and radial directions on the film surface of the film of condensate of variable thickness and can be recorded as follows:
For evaluation the circumferential velocity component in the turbulent diethyl ether condensate film under the vapor toroidal ring vortex in accordance with the Nikuradze’s experiments results  , can be represented as the following parabolas family:
where the exponent of the parabolas’ n depends on the Reynolds number Revp of the vapour flow and changes in the interval 1/6 … 1/10 (0.17 - 0.1).
The radial component of the liquid condensate film velocity uf under the vapor toroidal vortex ring can be defined as follows:
where tgφ―twist parameter of the surface layer of the liquid condensate film, formed on the HP’s condensation surface under the vapour toroidal vortex; χ― friction parameter of the toroidal vapour vortex ring.
The twist parameter tgφ decreases within the thickness of the condensate z, varies according to the quadratic parabola law   :
where ε―the tangent of the limiting (surface) angle of twist of the vapour flow is the ratio of surface shear stresses of friction in the circumferential (tangential) and radial directions on the condensate film surface.
In the equations for the momentum conservation of the vapour moving (7) and (8), the projection of the radial and tangential friction stress τs on the condensate film surface must be determined. To do this, we use the Prandtl hypothesis  of a two-layer flow model, including condensing two-phase vapour, introduce a criterion for the stability β of liquid condensate film flow under the vapour toroidal vortex ring in the manner of a two-dimensional vector with magnitude of the shear stress τs (with components τφs, τrs) on the film surface as follows:
where β―the criterion of diethyl ether film condensate flow stability under the vapour toroidal ring vortex; νl―kinematic coefficient of diethyl ether viscosity, m2/s.
The limit value of the working fluid film thickness δf under the toroidal vortex ring can be determined from the condition of velocities equality of the outside film surface and the toroidal vapour ring vortex. Given that the tangential constituent of the condensate film flow velocity provided the sticking on the inner surface of the top cover (v = 0) and can be represented as:
therefore, the calculation of the diethyl ether film flow rate can be carried out by the formula:
where μl―dynamic coefficient of diethyl ether viscosity, Pa・s; δ―condensate film thickness at immovable (stationary) vapour, m.
Measurements of the film flow rate were carried out earlier   , the solution of this equation is as follows:
After replacing the limit δf in Equation (22), we can obtain an expression for the radial component of the vapour―condensate film shear stress τrs as the following:
where Reδ―is the Reynolds number, calculated on the radial velocity uf and thickness δf of the diethyl ether condensate film flow.
The friction coefficient at the interfacial surface is determined in a standard way using the Blasius ratio:
When the film condensation occurs on the oriented in the horizontal direction HP’s top cover, the projection of gravity on the film flow direction is zero, and the provision of shear motion of the film occurs only due to the dynamic effects of vapour flow moving along the film outer surface. Motion in a thin film is considered to be inertialess.
Measurement average values results of the of the diethyl ether condensate film thickness on the cooled surface of the HP’s upper cover depending on the thermal power on the HP’s evaporator are shown in Figure 7. The absolute error of the thickness measuring does not exceed 1.5 × 10−3 mm   .
A sharply non-linear and decreasing values of the film thickness of the diethyl ether on the cooled smooth surface of top cover in dependence of the HP’s evaporator’s overheating compared to the boiling point of working fluid δT = Tev − TB. At low thermal power loads and the vapour vortex flow from the walls to the longitudinal axis of the vapour channel, condensate film has a large thickness. As the thermal power on the HP’s evaporator increases, the direction of the vortex flow changes to the opposite, the direction of vapour and condensate film flow become wake (in one direction), and the film thickness begins to decrease sharply.
Figure 7. Experimental value of the diethyl ether film thickness on the cooled in the calorimeter surface of the HP’s upper cover in depending on the amount of evaporator overheating relatively to the boiling point of diethyl ether as the working fluid, δT = Tev − TB, K at a semi-logarithmic scale.
The standard mean-square time average error of the film thickness measurement is ~ 1 × 10−3 mm, but begins to increase substantially starting from the evaporator superheat value δT = T - TB ~ (11 - 12) К, at which begins the boiling and a two-phase moist vapour form in the capillary porous evaporator. The condensing of the two-phase moist vapour leads to a noticeable increase in the experimental spread in the capacitive measurements of the diethyl ether film thickness on the cooled condensation surface.
The Reynolds number of the moving film of diethyl ether with thickness δ = 10−2 mm, νl = 0.032 × 10−4 m2/s, Prandtl number Pr = 0.77, at equal velocities of the vortex (30 m/s) and the film is equal to:
The calculated values of the film thickness are in the range of 10−2 - 10−4 mm when the value of friction coefficient is cf = 10−2 - 10−3, condensing vapor Reynolds number Revp = 103 - 104 and the vortex velocity up to 30 m/s. That is close to the experimentally obtained values of the film thickness of the diethyl ether condensate   .
1) For the first time by the calculation method, it was noted the possibility of changing the rotation direction of toroidal vortices of condensed vapour near the cooled HP’s condensation surface with Laval-liked nozzle vapour channel. At low values of the thermal power entering to the evaporator, the toroidal (axial) rotation of the vapour vortex due to the Coanda effect occurs from the periphery to the longitudinal axis of the vapour channel. A rotating toroidal vapour vortex acts on the moving condensate film and changes its thickness. With low thermal power and vortex rotation from the periphery to the center of the vapour channel, the condensate film movement due to the vortex friction about the surface is braked and slowing down, whereby the film thickness in the radial area of the upper cover effectively increases, has inflated values.
2) With increasing thermal power, the toroidal rotation direction of the vapour vortex changes to the opposite, from the z axis to the vapour channel walls and due to the friction on the outer surface of moving to the capillary-porous insert condensate film, the film radial flow velocity increases and its thickness decreases. The toroidal vortex additionally ‘blows off’ condensate film to the channel periphery along the HP’s upper flat cover surface.
3) The experimental results of the diethyl ether condensate film thickness depending on the HP’s thermal power (the evaporator overheat relative to the boiling point of the diethyl ether at atmospheric pressure) show a sharp decrease in the film thickness with an increase in thermal power. Such a sharp dependence of the film thickness on the thermal power indirectly confirms the change in the toroidal rotation direction of the HP’s vapour vortex. Increased film thickness (and HP’s high thermal resistance) at low loads and a sharp decrease in film thickness (and significantly reduced HP’s thermal resistance) with increasing thermal power can be associated with a change in the direction of the toroidal rotation of vapour vortex.
 Seryakov, A.V., Konkin, A.V. and Belousov, V.K. (2011) The Intensification of Heat-Transfer Characteristic of Heat Pipes. Proceedings of the VIII Minsk International Seminar of Heat Pipes, Heat Pumps, Refrigerators, Power Sources, Minsk, 12-15 September 2011, 59-65.
 Seryakov, A.V. (2016) Characteristics of Low Temperature Short Heat Pipes with a Nozzle-Shaped Vapour Channel. Russian Journal of Applied Mechanics and Technical Physics, 57, 69-74.
 Seryakov, A.V., Shakshin, S.L. and Alekseev, A.P. (2017) The Droplets Condensate Centering in the Vapour Channel of Short Low Temperature Heat Pipes at High Heat Loads. Journal of Physics: Conference Series, 891, 012123.
 Seryakov, A.V. (2017) The Study of Condensation Processes in the Low-Temperature Short Heat Pipes with a Nozzle-Shaped Vapour Channel. Engineering, 9, 190-240.
 Seryakov, A.V. (2016) The Application of Capacitance Transducer for Measuring Local Thickness of Condensate Film in Low-Temperature Range Heat Pipes. International Journal on Heat and Mass Transfer Theory and Application, 4, 1-13.
 Chung, Y.M., Luo, K.H. and Sandham, N.D. (2002) Numerical Study of Momentum and Heat Transfer in Unsteady Impinging Jets. International Journal of Heat Fluid Flow, 23, 592-600.
 Behnia, M., Parneix, S. and Durbin, P.A. (1998) Prediction of Heat Transfer in an Axisymmetric Turbulent Jet Impinging on a Flat Plate. International Journal of Heat and Mass Transfer, 41, 1845-1855.
 Craft, T.J., Iacovides, H. and Yoon, J.H. (2000) Progress in the Use of Non-Linear Two-Equation Models in the Computation of Convective Heat-Transfer in Impinging and Separated Flows. Flow, Turbulence and Combustion, 63, 59-80.
 Pohlhausen, K. (1921) Zur naherungweisen integration der differential-oleichung der laminaren reibungschicht. Zeitschrift fur Angewandte Mathematik und Mechanik, 1, 252-268.