The problem of turbulent fluid motion description has not been solved yet. It creates difficulties when oil, gas pipelines design calculation is performed. Besides, there are no theoretical methods for description of bodies motion in turbulent environment. These methods would be necessary for description of motion of aircrafts, submarines or above-water ships in the turbulent mode. Without simulation of bodies motion in wind tunnels or water basins, design of the bodies moving in the viscous environment is impossible.
There are approximate formulas for pipeline resistance coefficient at some ranges of Reynolds numbers, see [1, 2]. But they are empirical approximate formulas and they are applicable only for particular Reynolds number ranges.
Classical experimental Nikuradze curves of round pipeline resistance coefficient versus Reynolds number and roughness degree are well known. Approximation of convective term reducing Navier – Stokes problem to linear one with effective turbulent viscosity is applied. But such transformation distorts solution of Navier – Stokes equations and for matching to experiment the turbulent viscosity coefficient can have any value, up to negative. Galerkin method which brings hydrodynamic problem solution to system of non-linear ordinary differential equations is applied. But in the case of turbulent mode, this non-linear equation system has complex balance positions, i.e. the solution is complex. Indeed, hydrodynamics equation system, in turbulent mode, in real plane does not have any solutions, the equation solution tends for infinity, see [3] and main part I of the paper. But complex solution is finite. About physical meaning of the complex solution and oscillatory behavior of its imaginary part see [4, 5] or article III of this paper. Thus, it is necessary to solve hydrodynamic problems for turbulent mode in complex plane. At that, the turbulent solution is not single-valued, there are finite number of the solution branches.
Calculation of Round Section Pipeline Resistance Coefficient for Incompressible Fluid
This algorithm has been used for calculation of resistance coefficient for pipeline with round cross section. The algorithm is described in article [6] in English. We will seek the solution of problem for cross section round pipeline in form 
, in cylindrical system of coordinates. As external factor acts only along longitudinal axis
where P2, P1 – pressure in initial and final part of the pipeline; L – length of the pipeline, radial and angular velocities components are neglected. External action is equal to 
. According to formula (6), the pressure gradient is equal to 
. So we have the equation
Substituting velocity value we obtain
(1a)
Multiplying this equation by radius and performing integration over radius, as we use cylindrical coordinate system, we have
To obtain finite number of solutions we will multiply equation (1a) by function 
and integrate over volume. Then we will obtain finite number of turbulent solutions both for smooth and rough surfaces. Stationary laminar solution satisfying condition 
is single-valued as in the equation (1a) the laminar solution is identical for different values of 
. At the same time, likewise Schrodinger equation, finite number of turbulent solutions is found, each has its own energy. At transition from one state to another, discrete energy is radiated. The own energy minimum value defines the solution choice.
After calculating a module of the right part and averaging module of deviation angle tangent, we obtain
(1b)
It will be seen that when minus sign is chosen for value of average module of deviation angle tangent 
, roughness presence increases flow velocity as the full derivative 
increases and this is not correct, flow velocity has to decrease due to roughness presence.
When turbulent viscosity is taken into account, negative value of average velocity associated with process velocity correlation function 
, see [1], is used and this leads to plus sign for average module of roughness inclination tangent. The movement equation taking into account disturbances is
That is, convection term should be taken with minus, at right part of (1b) should be taken with plus.
Besides, it is necessary to choose plus for average module of roughness inclination tangent to obtain complex turbulent solution. Otherwise, solution describing pulse turbulent mode will not be steady.
Changing pipeline radius to diameter and dividing by value 
, we obtain
(2)
If you use another branch of root mean square unsteady solution and following equation will be obtained
(3)
Thus, steady solution for large difference in pressure is
Laminar solutions of these two equations at small pressure difference are the same. For turbulent mode with big pressure the solution has linear dependence of Reynolds number versus pressure square root. At small pressure increase, Reynolds number also grows and, as it follows from (3), pressure is increased. So the solution is not steady. In case of the complex solution it is equal to
At that, when pressure increases, imaginary part of velocity increases too and this does not lead to increase of real pressure, the real pressure keeps the value unchanged.
If micro roughness 
is distributed all over the pipeline surface, it is also present on macro roughness and defines critical Reynolds number and resistance coefficient at Reynolds number 2300. Micro roughness has the molecular nature, it is defined by average atom size equal to average geometrical difference between the nuclear size rA and size of Bohr orbit 
when the distance between atoms a = 3,043A is equal to some value determined by properties of pipeline boundary, iron, titanium and carbon. Distance between iron atoms is aFe = 2,87A, between titanium atoms – aTi = 3,46A, between carbon atoms – aC = 3,567A, see [7]. At the same time, the absolute value of tangent of micro roughness height inclination for metal surface of the pipeline is determined by formula
The average tangent of inclination is equal to
In this paper, critical Reynolds number was calculated with respect to radius. Critical Reynolds number with respect to diameter is equal to Rcr = 2300. But why critical Reynolds number for the sphere is equal to 3•105? This is due to different definition of critical Reynolds number. This value is equal to
where leff – effective hydrodynamic size of the body, including medium, a – true geometrical body size, and 
– molecular tangent of roughness inclination. And the ratio 
can be equal to 
.
Critical Reynolds number is equal to Rcr = 2300. Macro-roughness elements 
are rarer and this causes increase of resistance coefficient at Reynolds numbers which is 12 or more times more.
So we obtained a stationary criterion for Navier – Stokes equations taking into account one term of the solution series for one-dimensional case:
For one-dimensional case, on condition of pipeline cross section area constancy, the continuity equation is the same. Laminar solution of this equation is
For external pressure equal to 
, a complex solution and turbulent mode take place as Reynolds number from this point is equal to critical value. From experiment and calculation, we have critical Reynolds number for round pipeline
The pipeline resistance coefficient for round cross section pipe is determined by formula (we substituted to the formula the pressure difference expressed through dimensionless pressure)
The average velocity used for Reynolds number is equal to
The pipeline resistance coefficient λlam asymptotic for laminar mode in round cross section pipeline is calculated truly.
Asymptotic behavior of the pipeline resistance coefficient is obtained for small Reynolds numbers when the convective term is small.
In case of large pressure difference, we have a complex turbulent solution
Computing more precisely, contribution of rotary imaginary part to forward velocity of flow movement corresponds to square root of imaginary part according to formula (4)
(4)
and it is necessary to use value of ratio of Reynolds number to square root of dimensionless pressure as value of order 1 in the turbulent mode. At infinite pressure, Reynolds number for the flow is proportional 
. At that, the smoothest surface is the surface with average module of inclination tangent equal to inverse value of critical Reynolds number. For solution in the form of series, another value of α will be calculated. This value is defined from identical values of resistance coefficients at large Reynolds numbers and molecular roughness. The smoothest surface corresponds to average module of tangent of inclination equal to the inverse value of critical Reynolds number as the smallest modules of tangent of inclination correspond to molecular type of roughness. At that, effective diameter is less than true diameter. The average module of tangent of inclination angle can not be less than molecular roughness and its minimum value is equal to 
. That is, 1 is the maximum value of ratio of effective diameter to true diameter because α = 2. For external problem, effective diameter will increase, and the coefficient will be determined by formula
Coefficient β is proportional to
At zero macro roughness, effective diameter is equal to 1, that is, when roughness is increased, effective diameter decreases. Value 
was obtained from numerical experiment that corresponds to fourth root of mean square deviation. At zero macro roughness, micro roughness presents. And ratio of tangent of macro inclination roughness to micro roughness is more than 
At l/k = 30, we have value of effective pipeline diameter
At the same time, diameter is changed only for coefficient of pulsing part of the solution, i.e. for imaginary part from where the multiplier 
originates as the imaginary term is proportional to 
which is averaged. At that, value 
corresponds to fourth root of mean square deviation.
Here, influence of walls roughness in turbulent flow on imaginary part of Reynolds number of the flow is taken into account. To obtain curves with constant roughness height, it is necessary to enter effective average module of tangent of roughness inclination angle. The effective average module of tangent of roughness inclination angle has to depend on external pressure 
.
And in points of infinite Reynolds numbers or dimensionless pressure, we have the roughness corresponding to constant roughness height
where k – mean square root of the roughness height; r0 – radius of round cross section of the pipeline.
The formula is chosen in such a way that it defines correctly dependence of Reynolds number versus external pressure and pipeline resistance coefficient at infinite Reynolds numbers and external pressure
at resistance coefficient equal to
When average module of tangent of roughness inclination angle 
is constant but roughness height k is varying we obtain a curve which differs from Nikuradze curve.
Fig. 1. Curve of round pipeline resistance coefficient versus Reynolds number for different mean square root tangent of roughness inclination angle
But Nikuradze formula is obtained for constant ratio of pipeline radius r0 to average roughness height k. The formula (4) contains effective average module of tangent of roughness inclination angle expressed through ratio of pipeline radius to average roughness height using dimensionless pressure
Value 
. Influence of effective average module of tangent of roughness inclination on flow property depends on Reynolds number or pressure difference.
Empirical formula for finding of coefficients α(ξ0), β(ξ0), γ(ξ0) is following
At the same time, at the beginning of formation of the complex solution imaginary part 
, or at the beginning of turbulent solution, roughness inclination tangent is equal to approximately 1, and curves for different roughness inclination tangents coincide.
At that, flow resistance coefficient for round pipeline is determined by formula 
, Reynolds number calculated based on the average velocity of flow movement is equal to 
. Resistance coefficient at infinite pressure is proportional to 
. Here we demonstrate curves for solution obtained using one term of the series.
To compare theoretical and experimental curves of resistance coefficient dependence versus flow Reynolds number, experimental curve by Nikuradze is given in Fig. 2, on the right. Error of the theoretical curve relative to experimental one is about 10 %. But for laminar mode two solutions (2) and (3) are possible. Averaged solution will yield zero convection term and dependence 
that is not taken into account at the computation. In the theoretical curve convective term is taken into account which became equal to zero after averaging in laminar mode.
This curve was calculated for constant flow temperature over the flow cross section therefore in case of weak dependence of kinematic viscosity on temperature the formula will not change. For turbulent mode, it is necessary to substitute into the formula normalized pressure and ratio of pipeline radius to roughness height
Fig. 2. Calculated and measured dependence of round pipeline resistance coefficient versus Reynolds number for different roughness
And the formula is constructed so that 
. In case of the laminar mode there is a simple formula for Reynolds number:
Algorithm for Solution of Internal Hydrodynamic Problem for Arbitrary Flow Geometry
Navier – Stokes equations in Cartesian coordinates is
(5)
We will solve a three-dimensional laminar stationary problem without convective term for defined external action gl
Let us transform this equations to dimensionless form by dividing it by v2/d3, as a result we obtain dimensionless equation
Following function is a solution of this problem
We will seek solution of continuity equation for external action, where ri – response to external action
(6)
From this we obtain equation for finding of flow pressure
We will seek the pressure value in the form 
. Then we will substitute pressure into expression under the integral sign, multiply by φm(y1, y2, y3), and perform integration over the volume, then we obtain a system of linear equation
bm = Amnan.
Expressions for coefficients are
where hl(y1, y2, y3) is defined by external action. Let us transform Navier – Stokes equations to dimensionless form by dividing it by v2/d3, and we have dimensionless equation
Then we multiply Navier – Stokes equations by area of flow tube cross section, write the equations along laminar solution and enter flow tube with constant flow, see [8]. 
. In the convection term and in pressure gradient, we enter the derivative in the direction corresponding to the direction of flow lines in laminar solution. When substituting of the solution into equation
(7)
where Ss – flow tube cross section in laminar mode, expression Rs[y1(α, β), y2(α, β), y3(α, β)] is a stationary solution of Navier – Stokes equations without convection term which is equal to zero for flow tube as it does not depend on longitudinal coordinate.
We built these flow tubes for any external action which affects pressure difference. Further we consider roughness and under certain conditions obtain complex turbulent solution which is associated with influence of quadratic convection term with small multiplier, taking into account roughness, which yields complex solution at large pressure difference. At the same time, we reject real solution which was obtained for another sign of the module of average deviation, as it does not define fluctuating, turbulent solution. And imaginary part of the solution defines the solution pulsations.
If another sign of square root is chosen and correlation function of the process 
, where 
is a velocity deviation from its average value, is taken into account, turbulent viscosity becomes negative.
Let us substitute the solution (7) into Navie – Stokes equation, integrate it over flow pipe and divide by pipeline cross-sectional area. Then the convective term will be equal to
Taking roughness into account results in dependence of the pipeline radius a(s) on macro-roughness. Further we will extract the term da/ds associated with roughness and will find average value of its module. At the same time, we will make averaging of the equation with respect to s. It can be found out that convection term in laminar mode for smooth surface is equal to zero, and roughness has to be taken into account for non-zero value. So, we have the equation
To take into account roughness of pipeline surface and obtain turbulent solution, it is necessary to consider the average module of tangent of roughness inclination angle. Then this convective term will have a small multiplier, and the convection term will be non-zero. This term is proportional to average value of tangent of module inclination at roughness 
. At the same time, there is a term depending on variable pipeline cross section area 
. And flow lines of complex turbulent solution corresponding to flow lines of laminar solution will remain the same but there will be a solution pulsing around laminar flow lines. At that, the pulsations are defined by imaginary part of velocity, and the imaginary part of the solution, equal to a constant, means pulsations with amplitude equal to imaginary part of velocity.
Now, we will substitute the solution (7) into Navier – Stokes equations and will integrate along flow tubes, will multiply by Rcr in domain where this value meets a condition 
and where 
– average module of inclination tangent for not removable micro roughness with envelope forming macro-roughness, and we will obtain the following equation
where Rs(y1, y2, y3), p(y1, y2, y3) are determined from laminar solution and continuity equation, and function of external action hl(y1, y2, y3) is defined. So it was found out that micro roughness located along all length of the pipeline defines critical Reynolds number. This micro roughness is less than macro-roughness which affects resistance coefficient at large Reynolds numbers. But as Reynolds number depends on pipeline geometry through its diameter, then critical Reynolds number is inversely proportional to the average module of tangent of micro roughness inclination and depends on pipeline geometry. At the same time, reduction of pipeline radius results in negative da/ds value and, therefore, absence of complex turbulent solution in the narrow place, i.e. the critical Reynolds number raises. On the contrary, the pipeline widening causes increase of da/ds and, therefore, reduction of critical Reynolds number and can result in earlier occurrence of complex solution, i.e. the turbulent mode.
And, as Reynolds number depends on temperature through dependence of kinematic viscosity on temperature, it is obvious that occurrence of critical Reynolds number depends on environment temperature.
Coordinates of balance position are defined from a quadratic equation
At the same time, the laminar solution 
, as the convective term has different signs in laminar flow. In turbulent conditions the solution the convective term of one sign is not sustainable.
In this case, turbulent formula for roughness calculation is applicable due to identical averaging method in turbulent mode
where 
– effective average tangent of roughness inclination, ξ0 – ratio of roughness height to pipeline radius and critical Reynolds number 
is value of Reynolds number corresponding to the beginning of the complex solution. At the same time, for small Reynolds number we obtain a laminar solution. But difficulties in obtaining of turbulent solution do not come to an end. It is necessary to define effect of the surface roughness and for this use of experimental data is still inevitable. In principle, exact dependence of Reynolds number for smooth surface on macro-roughness is necessary to be learned. But external problem has some features associated with existence of resistance crisis which is caused by presence of a trace behind the body placed into the flow. This trace does not present in internal problems such as flow in pipeline.
Specificity of Flow Velocity Calculation for Sphere
Let us find out solution of Navier – Stokes equations for external problem. We have laminar solution for sphere motion in fluid for small Reynolds number. It yields the following velocity distribution, see [8]:
At that, pressure dependence on flow parameters is
Motion equations in spherical coordinate system for solutions which do not depend on angle φ can be written as
Let us change coordinate system to ξ, τ, θ with unknown Rr, Rθ, P, the coordinate system is defined by formula 
, after division of the equation system by 
we will have equation system
At that, in dimensionless constants, solution can be expressed as
But if you consider solution for one domain θ ∈ [0, π], zero value will be obtained for coefficient Rx. So, the domain should be divided into two parts θ ∈ [0, θ0], θ ∈ [θ0, π] and value θ0 should be found out of equality of Rx coefficients computed for different domains. At that Rx – common for either of Reynolds number components as laminar solution.
Here we will show how to find solution for the first equation, solution of the second equation can be found similarly. For this, for internal problem, we will multiply equation by r2sinθdrdθ. For external problem, we will enter variable 
for ξ ∈ [1, 0] and the multiplier will be following 
. Let us write down the equation with all multiplies:
Integration over the angle [0, π] yields zero right part of the equation. So, it is necessary to divide this solution into two domains and match solutions at the boundary. At low velocity, this solution will be real but it is possible that the angle is complex.
Let us integrate this equation over two domains θ ∈ [0, θ0], 
and θ ∈ [ θ0, π], 
, then
Equation for another domain is
For laminar mode and very small Reynolds number R0 << Rcr, we have following expression for Reynolds number
Solution obtained is symmetrical: θ0 = π/2, Rx = 0,8214. If non-linearity is taken into account:
Where parameter 
is defined for area of Reynolds number increase. Another solution is:
And complex Reynolds number Rx corresponds to beginning of turbulent mode.
If you take into account all coefficients, solution θ0 = π/2 will not be obtained but you will have two values for coefficient θ0. It will be found that two angles θ1, θ2 exist for each Navier – Stokes equation which correspond to two different variants of domain division. In case R0 → 0 angles θl = π/2 are equal, we have Rl(θl) = 1/2. Coefficients R1(θ1) = R2(θ2), R3(θ3) = R4(θ4) will be found from two Navier – Stokes equations which will be integrated separately over domains [0, θ0], [θ0, π]. At that, the two first of the angles will be found from the first Navier – Stokes equation, and the third and the fourth – from the second one.
Final solution will be found in the form
We substitute the decision in two equations of Navier – Stokes and in the continuity equation, we average on space and we find the stationary solution.
And Cartesian components of velocity are equal to
For the following examples initial data were taken which do not match the solution. Fig. 1 shows a plot for real angles versus two angles and on condition
Rl(θl) = 1; R0 = 1,5.
And for all plots ℜr = ℜθ = 1.
Fig. 3
Fig. 4
Fig. 5
Fig. 6
Next Fig. 4 shows results for Reynolds number R0 = 150, angles
Rl(θl) = 0,1.
The more is Reynolds number, the more is deviation of angles θl from π/2. Fig. 5 shows flow with Reynolds number R0 = 5000 and complex angles
Rl(θl) = 0,1.
Two singular domains are seen in front of the sphere and behind the sphere. In these areas, velocity corresponds to tangent line. Plot in Fig. 6 was calculated for the same parameters as plot in Fig. 5 but Reynolds number of the body is equal to R0 = 50000. The flow parameters are maximal, pattern remains the same as for parameter R0 = 5000.
Fig. 7 was plotted for parameters
Rl(θl) = 0,1; R0 = 500.
Fig. 7
Velocity distribution is so that there is a singular domain in front of the sphere – incompressible fluid can not penetrate into this area. And this area, inaccessible for fluid flow, has large length that provides conditions for origin of long vortex path.
To change pattern, it is necessary to change angular boundaries and relation between coefficients Rl(θl). Besides, at large Reynolds number, imaginary part increases and, hence roughness effect is rather large.
For incompressible liquid the equation of continuity along a current tube with longitudinal coordinate s has an equation 
. As the normal derivative from a normal component of speed is equal to zero for border of a special zone, we have constant longitudinal speed on border of a special zone. The convective term on border of a special zone is equal to zero. At that critical Reynolds number for external region off the body is equal to 
, where a – specific body size, lcr – length of the smooth body envelope when condition of complex coefficients Rl(θl) beginning is satisfied. Ratio 
is found from non-linear equation for lcr finding, which corresponds to beginning of complex solution.
For the plots computing, following equation system was resolved in dimensionless coordinate system
x0 = –2; y0 ∈ [–4, 4].
For this, we write down new formula which is necessary to substitute to Navier – Stokes and in the continuity equation, to average the solution and to define new multiplies ℜr, ℜθ, P by which the solution will be multiplied
Let us draw the curves for real boundaries of the area definition.
Vertical axis characterizes module of difference between coefficients calculated for two different areas. On horizontal axis the real angle θ0 is shown. In Fig. 8, the only root for small Reynolds number is shown. In Fig. 9, there are two real roots corresponding to the laminar mode with Reynolds number equal to R0 = 100.
Fig. 8
Fig. 9
In Fig. 10, 11 complex roots existence is shown, the roots are equal to
The imaginary axis values change in interval [0, 2], real axis values – in interval [0, π].
Description of Singular Domain
At that, solution for fluid flow has discontinuous zones, velocity perpendicular to boundaries of these zones is zero. Therefore fluid in these zones is independent of main flow. But tangential velocity components on boundary have to coincide. Now we will find the solution in these zones. Real part of the solution R = R1 + iR2 corresponds to component z, the imaginary part – to component x, and the x axis rotates around the axis 0z with change of angle φ. But the solution is to be found for fixed angle and should not be dependent of this angle. Then the solution of Navier – Stokes equation will be
(8)
Where new scaled angular variable 
is entered, where θmax, θmin – extreme values of turbulent zone boundaries. Besides, we will enter the scaled radius
where a max(θ) , a min(θ) – maximum and minimum value of radius of the turbulent zone boundary. In case if denominator is zero, value 
should be used for r. Then lnρ will be continuous and equal to π in this point. Coefficients bnm will be defined from values of the laminar solution within turbulent zone boundaries r = a min(θ); s = a max(θ), where θ∈ [θmin, θmax].
Fig. 10
Fig. 11
Coefficients bnm will be determined by formula
As boundary values at the beginning and the end of the period differ and area boundaries expressed in coordinates r, θ are not rectangular (in coordinates Φ, lnρ velocity on the boundary is variable), a series will be discontinuous, that is, the coefficient bnm decreases as 
when n, m → ∞, i.e. this solution is discrete. In singular domain, in coordinates lnρ, Φ, the solution is discrete due to discretization of functions R(lnρ0, Φ) in the form of discrete series. But as the description of singular domain is performed relative to coordinates lnρ0, Φ, the singular domain is discrete. Vortex path or pulsing turbulent mode with variable boundary is formed in this area at laminar mode.
The formula (8) can be rewritten in the form
(9)
where in this case we have
and then step with amplitude Anm and phase 
will be found from equations
where indexes n, m = –N, ..., –1, 1, ..., N.
It should be noted that A00 = b00. If the series in the left part of (9) is not summarized directly as this requires too large number of terms, then the right part of (9) will determine its discrete sum for finite number of terms. It should be noted that
is almost periodic coordinate of the step.
Why the turbulent solution in singular domain has the pulsing character with variable boundaries? The turbulent area boundary is not smooth function due to discreteness of the turbulent solution, unlike the laminar solution. This results in non-equality of tangential component of the solution and boundary pulsation in case of turbulent mode.
For description of laminar flow, it is necessary to enter dependence of specified radius on time
where Sh is a Strouchal number. At that, the pattern will fluctuate with Strouchal frequency according to value of ln ρ0 and this will lead to vortexes rotation in opposite directions as the frequencies under condition 
have different signs. At the same time, on the area boundary, frequency is zero, i.e. the solution on boundary is continuous in laminar mode.
Solution of the Flow Problem for Arbitrary Smooth Body in Spherical Coordinate System
Laminar solution of the flow problem for arbitrary body in spherical coordinate system we regard resolved in the form of final formula. That is, value of Reynolds number and pressure for laminar mode is found:
We resolve each Navier – Stokes equation by multiplying by 
, integration over inverse radius and angle θ, over two areas, which have one of the boundaries θl, l = 1, 2. We defined this boundary from equation 
. As the equation for these angles finding is the second degree one, two angles, θk1, θk2, are found. We define value θ0r(φ) for laminar solution and consider this in formula for Reynolds number taking area boundaries into account.
We do the same operation with other components of Reynolds numbers. Further we find out the solution by entering four unknown constants
(10)
We substitute these functions into Navier – Stokes equations and continuity equation, we integrate over the volume and then we obtain 4 constants ℜr, ℜθ, ℜφ, P. These coefficients can be complex describing the complex turbulent solution. Real part of the solution will be an average solution, and imaginary part – mean square deviation. At that, as the angle enters into solution function in non-linear way, it is possible to integrate on periodic angle φ without obtaining of zero integral. When solving non-linear equation, there can occur complex function θrl(φ), θθl(φ), θφl(φ), l = 1,2. Similarly, it is possible to find the problem solution for sphere, determining not laminar pressure, but such solution will be complicated. It is possible to add angle dependence of the sphere solution versus angle φ in Cartesian coordinate system and to solve a problem defining θ1(φ), θ2(φ), then dependence of the solution on angle φ will be found. At the same time, it is necessary to keep dependence on spherical coordinate system at Cartesian components versus velocity and pressure. In curvilinear coordinate system, the derivative is determined by formula
Where
From this we define 
through dependence 
. The second derivatives with respect to xl can be found similarly but in this case dependence on mixed derivatives with respect to r, θ, φ will occur.
At that, as 
velocity component Rφ will occur. As θxl = θ0x, θyl = θ0y, θzl = θ0z, this dependence vanishes at small Reynolds number.