*Article* **Obtaining Expressions for Time-Dependent Functions That Describe the Unsteady Properties of Swirling Jets of Viscous Fluid**

**Eugene Talygin \* and Alexander Gorodkov**

Bakulev National Medical Research Center for Cardiovascular Surgery, 121552 Moscow, Russia; euegene.talygin@skmenergy.com

**\*** Correspondence: skalolom@gmail.com

**Abstract:** Previously, it has been shown that the dynamic geometric configuration of the flow channel of the left heart and aorta corresponds to the direction of the streamlines of swirling flow, which can be described using the exact solution of the Navier–Stokes and continuity equations for the class of centripetal swirling viscous fluid flows. In this paper, analytical expressions were obtained. They describe the functions *C*0(*t*) and Γ0(*t*), included in the solutions, for the velocity components of such a flow. These expressions make it possible to relate the values of these functions to dynamic changes in the geometry of the flow channel in which the swirling flow evolves. The obtained expressions allow the reconstruction of the dynamic velocity field of an unsteady potential swirling flow in a flow channel of arbitrary geometry. The proposed approach can be used as a theoretical method for correct numerical modeling of the blood flow in the heart chambers and large arteries, as well as for developing a mathematical model of blood circulation, considering the swirling structure of the blood flow.

**Keywords:** potential swirling flow; Navier–Stokes equations; unsteady swirling flow; tornadolike jets

#### **1. Introduction**

In our previous work, we investigated the dynamic geometry of the left heart and aorta to find consistency between the configuration of the flow channel and the directions of the swirling streamlines of the TLJ (Tornado-Like Jets) class [1–3]. The found correspondence allowed us to assume that the blood flow in the heart and large vessels belong to the TLJ class. These flows were described using the explicit solutions [4], which determine the conditions for a swirling jet of Newtonian fluid in space to appear and evolve.

It has been experimentally proved that swirling flows of this class, under certain conditions, are formed on a concave surface streamlined by the viscous medium. Apparently, an important role in the formation of such a swirling flow is played by the vortex boundary layer arising on the concave surface. This layer should consist of some small-scale vortex structures such as Taylor–Görtler vortices and differ in their properties from the classical shear boundary layer of L. Prandtl. These vortices are cylindrical structures, the axis of rotation of which is parallel to the incoming flow. This shape allows the swirling flow to rely on these vortex structures, conjugating with them only at one point. In this case, the movement of the swirling flow relative to the concave surface causes the appearance of rolling stresses, which are much less than the shear stresses in the boundary layer of L. Prandtl. The geometric shape of the generatrix of the concave surface determines the shape of the streamlines of the swirling flow. In this case, the vortex boundary layer is able to change its thickness along the concave surface dynamically. This allows it to compensate for local inconsistencies between the real surface and the geometric shape built along streamlines. Such compensation is necessary when the swirling blood flow moves into a section of the vascular bed with pathological geometry disorders.

**Citation:** Talygin, E.; Gorodkov, A. Obtaining Expressions for Time-Dependent Functions That Describe the Unsteady Properties of Swirling Jets of Viscous Fluid. *Mathematics* **2021**, *9*, 1860. https:// doi.org/10.3390/math9161860

Academic Editor: Mauro Malvè

Received: 23 June 2021 Accepted: 29 July 2021 Published: 5 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Obviously, blood has complex rheology and, in the general case, cannot be regarded as a Newtonian fluid. However, the available experimental data allow us to assert that in large-caliber vessels at relatively high speeds, the dynamic viscosity of blood hardly changes, and blood can be considered a Newtonian fluid.

Calculating the Reynolds number for the blood flow in the heart and great vessels, we obtain a value that allows us to make an unambiguous conclusion about the turbulent nature of the flow under consideration. However, a turbulent flow is characterized by active mixing of the liquid volume and spontaneous vortex formation. It is obvious that these properties of turbulent flows do not correspond to the physiological characteristics of blood circulation. Considering the blood flow as TLJ, the flow structure can be preserved without noticeable swirls and perturbations at such Reynolds numbers that, for other principles of flow organization, would lead to flow turbulization. However, it was not possible to take into account the nonstationarity blood flow because, in the exact solution, it is determined by arbitrary time-dependent functions that do not have a formal analytical record. At the same time, it is obvious that the real blood flow is strictly unsteady. A formal description of the unsteady properties of the swirling blood flow in terms of analytical functions would make it possible to construct a model of blood circulation with a higher degree of accuracy.

Earlier it was found that the geometric characteristics of the bloodstream (if we assume that they correspond to the geometric characteristics of the flow channel) are described with a high degree of accuracy by quasi-stationary solutions [1–4]. In these solutions, the non-stationary properties of the flow are determined by the behavior of the time-dependent functions *C*0(*t*) and Γ0(*t*). These functions could be determined experimentally if measuring the vector field of flow velocities. However, so far, making such a measurement with sufficient accuracy is rather difficult. The task is complicated by the fact that the bloodstream at all stages of its evolution interacts with the movable walls of the flow channel, taking the direction of movement given by the instantaneous geometric configuration of the channel. At the same time, the bloodstream is a submerged stream, and its structure changes upon interaction and merging with residual blood volume that retains a certain movement after the previous cardiac cycle.

In previous studies, we considered the swirling blood flow in the heart and great vessels as quasi-stationary—a continuous set of values of the functions *C*0(*t*) and Γ0(*t*) on a time interval was replaced by a discrete subset of values, each of which corresponded to a certain configuration of the flow channel. This made it possible to replace the functions *C*0(*t*) and Γ0(*t*) with a set of constants that reflect the nonstationarity behavior of the flow. The analysis of the dynamic geometry of the left heart and aorta confirmed this hypothesis; the values of some parameters of the swirling blood flow in this segment of blood circulation were calculated. However, the accuracy of the obtained results was limited by the resolution of discretization.

Obtaining formal relations for the functions *C*0(*t*) and Γ0(*t*) will make it possible to unambiguously relate the structure of an unsteady swirling flow with the conditions that form it—the dynamics of the inflowing blood flow, the geometry of the flow channel, and the interaction with residual volumes at each stage of the evolution of the bloodstream both in time and along the length of the flow channel.

Therefore, the aim of this work was to obtain formal mathematical expressions for the functions *C*0(*t*) and Γ0(*t*), which describe in general form the unsteadiness of the swirling flow. Due to the calculation difficulties, we could not validate the proposed equations with the experimental results. This work will be conducted in the next stage of our study.


• References.

#### **2. Materials and Methods**

According to [1–4], velocity components of the TLJ in a cylindrical coordinate system are expressed by the following relations:

$$\begin{cases} \begin{aligned} \boldsymbol{\mu}\_{z} &= 2\mathbf{C}\_{0}(t)z \\ \boldsymbol{\mu}\_{r} &= -\mathbf{C}\_{0}(t)\boldsymbol{r} \\ \boldsymbol{\mu}\_{\boldsymbol{\theta}} &= \frac{\Gamma\_{0}(t)}{2\pi r} \* \left(1 - e^{-\frac{\boldsymbol{C}\_{0}(t)r^{2}}{2\overline{v}}}\right) \end{aligned} & \text{(1a)} \end{cases}$$

where *uz*, *ur*, and *u<sup>ϕ</sup>* are the velocity components, ν is the kinematic viscosity, and the unsteadiness and evolution of the flow are determined by the behavior of the time-dependent functions *C*0(*t*) and Γ0(*t*). *C*0(*t*) is an arbitrary function of time, which, by its meaning, represents the gradient of the longitudinal component of the velocity (sec−1); Γ0(*t*) is an arbitrary function of time, corresponding to the physical meaning of the circulation of the medium (m2/sec).

It can be seen from the presented evidence that the unsteady properties of the flow under study are determined by the behavior of the time-dependent functions *C*0(*t*) and Γ0(*t*). While obtaining relations for the swirling flow velocity (1a) from the original system of Navier–Stokes equations and continuity, the functions *C*0(*t*) and Γ0(*t*) were considered as arbitrary functions of time.

These ratios were analyzed using the methods of linear differential equations solving and standard methods of differentiation for chain functions with several arguments.

#### **3. Results**

To obtain an analytical form of the functions *C*0(*t*) and Γ0(*t*), relations (1a) were substituted into the Navier–Stokes equations.

For this, the equations for the three velocity components were considered separately in a cylindrical coordinate system.

Azimuthal velocity component of the Navier–Stokes equation can be expressed as follows [5]:

$$\begin{cases} \frac{u\_r u\_\theta}{r} + \frac{\partial u\_\theta}{\partial r} u\_r + \frac{1}{r} \frac{\partial u\_\theta}{\partial \varphi} u\_\varphi + \frac{\partial u\_\theta}{\partial z} u\_z + \frac{\partial u\_\theta}{\partial t} = F\_\varphi - \frac{1}{\rho r} \frac{\partial p}{\partial \varphi} + \nu \left( \Delta u\_\varphi - \frac{u\_\theta}{r^2} + \frac{2}{r^2} \frac{\partial u\_r}{\partial \varphi} \right) \\ \qquad \qquad \qquad \qquad \Delta = \frac{\partial^2}{\partial r^2} + \frac{1}{r} \frac{\partial}{\partial r} + \frac{1}{r^2} \frac{\partial^2}{\partial \varphi^2} + \frac{\partial^2}{\partial z^2} \end{cases} \tag{1b}$$

In the case of axial symmetry of the considered swirling flow (this condition means that all partial derivatives with respect to *ϕ* must be equal to 0), the Equation (1b) has the following form:

$$\frac{u\_r u\_\phi}{r} + \frac{\partial u\_\phi}{\partial r} u\_r + \frac{\partial u\_\phi}{\partial z} u\_z + \frac{\partial u\_\phi}{\partial t} = \upsilon \left( \frac{\partial^2 u\_\phi}{\partial r^2} + \frac{1}{r} \frac{\partial u\_\phi}{\partial r} + \frac{\partial^2 u\_\phi}{\partial z^2} - \frac{u\_\phi}{r^2} \right) \tag{1c}$$

It is seen from relations (1a) that *uϕ* does not depend on z. This means that the partial derivative *uϕ* with respect to *z* equals 0. Then Equation (1c) can be changed as follows:

$$\frac{u\_r u\_\phi}{r} + \frac{\partial u\_\phi}{\partial r} u\_r + \frac{\partial u\_\phi}{\partial t} = \upsilon \left( \frac{\partial^2 u\_\phi}{\partial r^2} + \frac{1}{r} \frac{\partial u\_\phi}{\partial r} - \frac{u\_\phi}{r^2} \right) \tag{1d}$$

We substitute relations (1a) into Equation (1d), writing out each term separately:

$$\frac{\partial u\_{\varphi}}{\partial t} = \frac{d\Gamma\_0(t)}{dt} \ast \frac{1}{2\pi r} \ast \left(1 - e^{-\frac{\zeta\_0(t)r^2}{2v}}\right) + \frac{\Gamma\_0(t)r}{4\pi v} \ast \frac{d\mathbb{C}\_0(t)}{dt} \ast e^{-\frac{\zeta\_0(t)r^2}{2v}}$$

$$\frac{u\_r u\_\varphi}{r} = -\frac{\mathbb{C}\_0(t)\Gamma\_0(t)}{2\pi r} \ast \left(1 - e^{-\frac{\mathbb{C}\_0(t)r^2}{2r}}\right)$$

$$\frac{\partial u\_\varphi}{\partial r} u\_r = -\mathbb{C}\_0(t)\tau \ast \left(\frac{\mathbb{C}\_0(t)\Gamma\_0(t)}{2\pi v} \ast e^{-\frac{\mathbb{C}\_0(t)r^2}{2r}} - \frac{\Gamma\_0(t)}{2\pi r^2} \ast \left(1 - e^{-\frac{\mathbb{C}\_0(t)r^2}{2r}}\right)\right)$$

$$\frac{\partial^2 u\_\varphi}{\partial r^2} = \frac{\Gamma\_0(t)}{\pi r} \ast \left[\frac{1 - e^{-\frac{\mathbb{C}\_0(t)r^2}{2r}}}{r^2} + \frac{\mathbb{C}\_0(t) \ast e^{-\frac{\mathbb{C}\_0(t)r^2}{2r}} \ast \left(1 - e^{-\frac{\mathbb{C}\_0(t)r^2}{2r}}\right)}{2v} - \frac{\mathbb{C}\_0(t) \ast e^{-\frac{\mathbb{C}\_0(t)r^2}{2r}}}{v}\right]$$

$$\frac{1}{r}\frac{\partial u\_\varphi}{\partial r} = \frac{\mathbb{C}\_0(t)\Gamma\_0(t)}{2\pi rvr} \ast e^{-\frac{\mathbb{C}\_0(t)r^2}{2r}} - \frac{\Gamma\_0(t)}{2\pi r^3} \ast \left(1 - e^{-\frac{\mathbb{C}\_0(t)r^2}{2r}}\right)$$

$$\frac{u\_\varphi}{r^2} = \frac{\Gamma\_0(t)}{2\pi r^3} \ast \left(1 - e^{-\frac{\mathbb{C}\_0(t)r^2}{2r}}\right)$$

Substituting the obtained relations into the original Equation (1d), the left-hand side is transformed as follows:

$$\begin{split} &\frac{d\Gamma\_{0}(t)}{dt} \ast \frac{1}{2\pi r} \ast \left(1 - e^{-\frac{\mathsf{C}\_{0}(t)r^{2}}{2\overline{r}}}\right) + \frac{\Gamma\_{0}(t)r}{4\pi\upsilon} \ast \frac{d\mathbb{C}\_{0}(t)}{dt} \ast e^{-\frac{\mathsf{C}\_{0}(t)r^{2}}{2\overline{r}}} - \frac{\mathsf{C}\_{0}(t)\Gamma\_{0}(t)}{2\pi r} \ast \left(1 - e^{-\frac{\mathsf{C}\_{0}(t)r^{2}}{2\overline{r}}}\right) - \mathsf{C}\_{0}(t)r \\ &\ast \left(\frac{\mathsf{C}\_{0}(t)\Gamma\_{0}(t)}{2\pi r} \ast e^{-\frac{\mathsf{C}\_{0}(t)r^{2}}{2\overline{r}}} - \frac{\Gamma\_{0}(t)}{2\pi r^{2}} \ast \left(1 - e^{-\frac{\mathsf{C}\_{0}(t)r^{2}}{2\overline{r}}}\right)\right) \\ &= \frac{d\Gamma\_{0}(t)}{dt} \ast \frac{1}{2\pi r} \ast \left(1 - e^{-\frac{\mathsf{C}\_{0}(t)r^{2}}{2\overline{r}}}\right) + \frac{\Gamma\_{0}(t)r}{4\pi v} \ast \frac{d\mathbb{C}\_{0}(t)}{dt} \ast e^{-\frac{\mathsf{C}\_{0}(t)r^{2}}{2\overline{r}}} - \frac{\mathsf{C}\_{0}(t)\Gamma\_{0}(t)r}{2\pi v} \ast e^{-\frac{\mathsf{C}\_{0}(t)r^{2}}{2\overline{r}}} \end{split}$$

Similarly, the right-hand side of Equation (1d) can be transformed as follows:

$$\begin{split} \frac{\Gamma\_{0}(t)v}{\pi r} \ast \left[ \frac{1 - e^{-\frac{\mathsf{C}\_{0}(t)r^{2}}{2\nu}}}{r^{2}} + \frac{\mathsf{C}\_{0}(t) \ast e^{-\frac{\mathsf{C}\_{0}(t)r^{2}}{2\nu}} \ast \left(1 - e^{-\frac{\mathsf{C}\_{0}(t)r^{2}}{2\nu}}\right)}{2v} - \frac{\mathsf{C}\_{0}(t) \ast e^{-\frac{\mathsf{C}\_{0}(t)r^{2}}{2\nu}}}{v} \right] + \frac{\mathsf{C}\_{0}(t)\Gamma\_{0}(t)}{2\pi r} \ast e^{-\frac{\mathsf{C}\_{0}(t)r^{2}}{2\nu}} \\ -\frac{\Gamma\_{0}(t)v}{2\pi r^{3}} \ast \left(1 - e^{-\frac{\mathsf{C}\_{0}(t)r^{2}}{2\nu}}\right) - \frac{\Gamma\_{0}(t)v}{2\pi r^{5}} \ast \left(1 - e^{-\frac{\mathsf{C}\_{0}(t)r^{2}}{2\nu}}\right) = -\frac{\mathsf{C}\_{0}(t)\Gamma\_{0}(t)}{2\pi r} e^{-\frac{\mathsf{C}\_{0}(t)r^{2}}{2\nu}} \end{split}$$

The final Equation (1d) will be written as follows:

$$\begin{cases} \frac{d\Gamma\_{0}(t)}{dt} \frac{1}{2\pi r} \left(1 - e^{-\frac{\mathbb{C}\_{0}(t)r^{2}}{2\nu}}\right) + \frac{d\mathbb{C}\_{0}(t)}{dt} \frac{\Gamma\_{0}(t)r}{4\pi\nu} e^{-\frac{\mathbb{C}\_{0}(t)r^{2}}{2\nu}} - \frac{\mathbb{C}\_{0}^{2}(t)\Gamma\_{0}(t)r}{2\pi\nu} e^{-\frac{\mathbb{C}\_{0}(t)r^{2}}{2\nu}} = \\ \quad - \frac{\mathbb{C}\_{0}(t)\Gamma\_{0}(t)}{2\pi r} e^{-\frac{\mathbb{C}\_{0}(t)r^{2}}{v}} \end{cases} \tag{1e}$$

The Navier-Stokes equation for calculating the radial velocity component is written as follows [5]:

$$\frac{\partial \boldsymbol{u}\_{\rm tr}}{\partial t} + \boldsymbol{u}\_{\rm tr} \frac{\partial \boldsymbol{u}\_{\rm tr}}{\partial r} - \frac{\boldsymbol{u}\_{\rm \varphi}}{r} \frac{\partial \boldsymbol{u}\_{\rm r}}{\partial \boldsymbol{\varrho}} - \frac{\boldsymbol{u}\_{\rm \varphi}^{2}}{r} + \boldsymbol{u}\_{\rm \tilde{z}} \frac{\partial \boldsymbol{u}\_{\rm r}}{\partial z} = -\frac{1}{\rho} \frac{\partial p}{\partial r} + \nu \left( \frac{\partial^{2} \boldsymbol{u}\_{\rm r}}{\partial r^{2}} + \frac{1}{r} \frac{\partial \boldsymbol{u}\_{\rm r}}{\partial r} - \frac{\boldsymbol{u}\_{\rm r}}{r^{2}} + \frac{\partial^{2} \boldsymbol{u}\_{\rm r}}{\partial z^{2}} \right) \tag{2a}$$

Considering the axisymmetric flow, we have:

$$\frac{\partial u\_r}{\partial t} + u\_r \frac{\partial u\_r}{\partial r} - \frac{u\_\phi^2}{r} + u\_z \frac{\partial u\_r}{\partial z} = -\frac{1}{\rho} \frac{\partial p}{\partial r} + \nu \left( \frac{\partial^2 u\_r}{\partial r^2} + \frac{1}{r} \frac{\partial u\_r}{\partial r} - \frac{u\_r}{r^2} + \frac{\partial^2 u\_r}{\partial z^2} \right).$$

Substituting the expressions for the velocity Components (1a) into this equation, we get:

$$-r\frac{d\mathbb{C}\_0(t)}{dt} + \mathbb{C}\_0^2(t)r - \frac{\Gamma\_0(t)^2}{4\pi^2 r^3} \left(1 - e^{-\frac{\mathbb{C}\_0(t)^2}{2r}}\right)^2 = -\frac{1}{\rho}\frac{\partial p}{\partial r} \tag{2b}$$

The Navier-Stokes equation for calculating the longitudinal velocity component will then be written as follows [5]:

$$\frac{\partial u\_z}{\partial t} + u\_r \frac{\partial u\_z}{\partial r} + \frac{u\_\theta}{r} \frac{\partial u\_z}{\partial \phi} + u\_z \frac{\partial u\_z}{\partial z} = -\frac{1}{\rho} \frac{\partial p}{\partial z} + \nu \left( \frac{\partial^2 u\_z}{\partial r^2} + \frac{1}{r} \frac{\partial u\_z}{\partial r} + \frac{1}{r^2} \frac{\partial^2 u\_z}{\partial \phi^2} + \frac{\partial^2 u\_z}{\partial z^2} \right) \tag{3a}$$

In the case of axial symmetry of the considered swirling flow, Equation (3a) is written as follows:

$$\frac{\partial \mu\_z}{\partial t} + \mu\_r \frac{\partial \mu\_z}{\partial r} + \mu\_z \frac{\partial \mu\_z}{\partial z} = -\frac{1}{\rho} \frac{\partial p}{\partial z} + \nu \left( \frac{\partial^2 \mu\_z}{\partial r^2} + \frac{1}{r} \frac{\partial \mu\_z}{\partial r} + \frac{\partial^2 \mu\_z}{\partial z^2} \right).$$

Substituting relations (1a) into Equation (3a), writing out each term separately, we get:

$$2z\frac{d\mathbb{C}\_0(t)}{dt} + 4\mathbb{C}\_0^2(t)z = -\frac{1}{\rho}\frac{\partial p}{\partial z} \tag{3b}$$

Thus, finding the functions *C*0(*t*), Γ0(*t*) is reduced to solving the following system of differential equations:

$$\begin{cases} \begin{array}{c} \frac{d\Gamma\_{0}(t)}{dt} \frac{1}{2\pi r} \left(1 - e^{-\frac{\mathsf{C}\_{0}(t)r^{2}}{2r}}\right) + \frac{\Gamma\_{0}(t)r}{4\pi v} e^{-\frac{\mathsf{C}\_{0}(t)r^{2}}{2v}} \left(\frac{d\mathsf{C}\_{0}(t)}{dt} - 2\mathsf{C}\_{0}^{2}(t)\right) = -\frac{\mathsf{C}\_{0}(t)\Gamma\_{0}(t)}{2\pi r} e^{-\frac{\mathsf{C}\_{0}(t)r^{2}}{v}} \\ \qquad\qquad - r\frac{d\mathsf{C}\_{0}(t)}{dt} + \mathsf{C}\_{0}^{2}(t)r - \frac{\Gamma\_{0}(t)^{2}}{4\pi^{2}r^{3}} \left(1 - e^{-\frac{\mathsf{C}\_{0}(t)r^{2}}{2v}}\right)^{2} = -\frac{1}{\rho} \frac{\partial p}{\partial r} \\ \qquad \qquad 2z\frac{d\mathsf{C}\_{0}(t)}{dt} + 4\mathsf{C}\_{0}^{2}(t)z = -\frac{1}{\rho} \frac{\partial p}{\partial z} \end{array} \tag{4}$$

It should be noted that the equations in the written system (4), in addition to the sought-for time-dependent functions, also contain the derivatives of pressure with respect to coordinates and these coordinates. Thus, the conditions for which the exact solutions were derived are satisfied if and only if any function of the partial derivatives of pressure with respect to coordinates and of the coordinates themselves is a function only of time.

The last equation of system (4) contains only one unknown function—*C*0(*t*); therefore, such an equation is most conveniently solved using the method of Lie groups. The solution to this equation gives the following expression for the desired function:

$$C\_0(t) = -\frac{\varkappa\_1}{2\text{tanh}(a\_1\*(C\_1-t))}\tag{5a}$$

where *C*<sup>1</sup> is the constant of integration and can be taken as 0.

Then the resulting expression for *C*0(*t*) was substituted into the second equation of the system (4). After this, the equation is a differential equation from one unknown function Γ0(*t*). Its solution gives the following result:

$$\Gamma\_0(t) = \frac{\pi \sqrt{\frac{\alpha\_3 \ast \tanh^2(a\_1 \ast (C\_1 - t)) + a\_4}{\exp\left(\frac{a\_2}{4\nu \ast \tanh(a\_1 \ast (C\_1 - t))}\right) - 1}}}{\tanh(a\_1 \ast (C\_1 - t))} \tag{5b}$$

The following notations are introduced in Expressions (5a) and (5b):

$$\begin{cases} \begin{array}{c} \alpha\_{1} = \sqrt{-\frac{1}{\rho z} \ast \frac{\partial p}{\partial z}} \left[ \frac{1}{s} \right] \\ \alpha\_{2} = r^{2} \ast \sqrt{-\frac{1}{\rho z} \ast \frac{\partial p}{\partial z}} \left[ \frac{m^{2}}{s} \right] \\ \alpha\_{3} = \frac{4r^{3}}{\rho} \frac{\partial p}{\partial r} - \frac{2r^{4}}{\rho z} \ast \frac{\partial p}{\partial z} \left[ \frac{m^{4}}{s^{2}} \right] \\ \alpha\_{4} = \frac{3r^{4}}{\rho z} \ast \frac{\partial p}{\partial z} \left[ \frac{m^{4}}{s^{2}} \right] \end{array} \tag{5c}$$

If the pressure change along the longitudinal and radial coordinates is constant, the parameters (*α*1, ... , *α*4) show only the change in the geometric configuration of the jet over time. If the jet is enclosed in a flow channel that takes its shape, then the parameters (*α*1, ... , *α*4) can be determined experimentally by measuring the dimensions of the channel.

In [6,7], the empirical dependences on time of Γ0(*t*) and of *C*0(*t*)/Γ0(*t*) were obtained. Comparing these dependencies, the approximate order of values of the constants was determined as (*α*1, ..., *α*4) = (10<sup>−</sup>1, 10−6, 10−4, 10−2).

Using the suggested values of the parameters, the graphs of the dependence of the functions *C*0(*t*) and Γ0(*t*) on time were plotted (Figures 1 and 2).

**Figure 2.** Plot of the Γ0(*t*) function.

#### **4. Discussion**

As a result, we managed to obtain a general analytical form of the functions *C*0(*t*) and Γ0(*t*). The ratios for these functions contain several parameters, the exact value of which should be determined within the framework of a specific problem. As can be seen from the graphs, the studied functions *C*0(*t*) and Γ0(*t*) are strictly positive and decrease monotonically with time. This behavior expresses the fulfillment of the energy conservation law for the considered swirling flow. Indeed, solution (1a) shows that if the values of the functions under study could increase in absolute value over time, the kinetic energy of the swirling flow would have increased in the absence of external action, which is impossible. The dynamics of the functions *C*0(*t*) and Γ0(*t*) approximately corresponds to the experimental curves plotted, based on studies of the blood flow structure in the left ventricle [7].

Analysis of the constants (*α*1, ..., *α*4) allows the following statements:

The dimensions of the constants correspond to the dimensions of the studied functions *C*0(*t*) and Γ0(*t*). In particular, the dimension *α*<sup>1</sup> coincides with the dimension *C*0(*t*), the dimension *α*<sup>2</sup> coincides with the dimension Γ0(*t*) and the dimensions *α*3, *α*<sup>4</sup> coincide with the dimensions Γ0(*t*) 2 .

The swirling flow belongs to the TLJ class and can be described by exact solutions (1a) if and only if the form of the dynamic pressure dependence on the channel coordinates does not change with time.

Two swirling flows belonging to the TLJ class and having different geometric characteristics can be determined by the same functions *C*0(*t*) and Γ0(*t*) and, accordingly, have an identical dynamic structure if the set of constants (*α*1, ..., *α*4) is the same for both.

Therefore, the obtained relations (5a) and (5b) make it possible to describe in a general form the unsteady swirling flow of a viscous fluid, the velocity components of which are expressed by solution (1a). Considering the obtained expressions for the functions *C*0(*t*) and Γ0(*t*), the exact solution can be used for numerical modeling of the swirling blood flow in the heart and large vessels without considering the non-Newtonian properties of blood.

Attempts to formally describe the unsteady blood flow have been undertaken earlier by many researchers [8–11]. Computational fluid dynamics methods were mainly used in these works. With the development of computing power, it has become possible to simulate the geometry of heart cavities using neural networks [12]. However, despite great advances in the field of numerical modeling, these models turn out to be rather complex and lose the necessary clarity. This is even more true for models based on neural networks.

In this work, based on the well-studied non-stationary equations of hydrodynamics, relatively simple relations were obtained that describe the non-stationary properties of the swirling blood flow in the heart and great vessels. The use of the obtained ratios in CFD models and for calculating initial approximations of the geometry of flow channels for training neural networks will enable the development of relatively simple and effective blood circulation models.

#### **5. Conclusions**

Most modern circulatory models are based on numerical modeling of blood flow and based on methods of computational hydrodynamics and deep learning [8–12]. This approach simplifies the calculations but does not allow derivation of the physical and mathematical meaning of the functional dependencies, which are included in the ratios for the components of the blood flow velocity. So, it was possible to obtain analytical relations for functions reflecting the unsteady properties of the swirling blood flow. These relations, on the one hand, make it possible to perform a more accurate analysis of the physical characteristics of an unsteady swirling blood flow. On the other hand, they allow formulating characteristic flow parameters that can be verified by experiment. Thus, the results of this study allow us to deepen the physical understanding of the peculiarities of the evolution of the swirling blood flow, as well as to provide convenient characteristic parameters for validating experimental data.

It should be noted that the obtained ratios do not consider the pulsating nature of the blood flow in the heart and large vessels. The necessary improvements are planned to be added to the model at subsequent stages of the study.

The next steps in the development of this study will be to carry out experimental measurements of the field of velocities and pressures and to validate the obtained results using the functional relationships obtained in this work.

**Author Contributions:** Development of the concept of the article and collection of the material, E.T.; Development of a research algorithm, finding suitable mathematical methods for solving the assigned tasks, E.T. and A.G.; Carrying out mathematical transformations, solving equations, E.T., Description of the physiological side of the study, A.G.; Supervision and Project Administration, A.G.; Writing, E.T. and A.G. Both authors have read and agreed to the published version of the manuscript.

**Funding:** The study was supported by the Russian Scientific Foundation (grant No. 16-15-00109).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

