Methods of computational modeling of coronary heart vessels for its digital twin

In this work, methods of numerical modelling of the coronary vessels system of the human heart have been studied. This investigation includes transient flow of the liquid – blood and dynamics of zones of shear stress at vessels. The main goal of the research is obtaining of hemodynamic and shear stress for creating the digital twin of coronary heart vessels. The results were obtained for low Reynolds numbers about 20 of three-dimensional laminar flow. With this Reynolds number the turbulent flow of the blood is modelled by Realizable k-ε model, and SST models to the narrowing, expansions, and blocks inside the vessels. Loads caused by the additional energy consumption because of the turbulent flow of the blood (increase in arterial blood pressure) have been analyzed. A twodimensional model of a separated vessel with fixed blood back-flow prevention is developed. Presence of a turbulent flow core is discovered. By the means of stress-strain properties of the model, visual representation of the wearing process of the blood back-flow preventer, and heart diseases progression


Introduction
The concept of the digital twin, suggested in 2002, carries the idea of fusion between the virtual and the physical worlds. Digital twin is a model of a physical product in the field of production and design. «The vision of the digital twin itself includes more or less all information which could be useful in allthe current and subsequentlifecycle phases» [1], such vision is represented in the work of S. Boschert and R. Rosen. Moreover, E. J. Tuegel and co-authors confirmed it. [2] .In their article, a model represented the idea of how a digital twin can be used to predict the lifetime of an aircraft structure and to ensure its structural integrity. In other words, the value of a digital twin is expressed in predicting local damage through the calculations of temperatures and structural deflections in conditions met during a flight. Nowadays, offers of methods for integrating digital twins are being put forward from all over the world. For example, Ramtilak A. [3] in his work describes the Digital Twin Spark Ignition (DTS-i) model successfully launched on two products -150 DTS-i and 180 DTS-i. The main purpose of these models is to improve the characteristics of four-stroke cycle engines in a wide area of usage and to comply with the 2005 emission standard in India.
Anyways, the implementation of ideas embedded in the concept of a digital twin, requires the extensive use of numerical methods of computer modeling and CFDmodeling (Computational Fluid Dynamics) [4]. This method provides the possibility of obtaining a visual representation of the nature of the processes in the object under study and the possibility of investigating its various designs without creating expensive experimental installations. In 2017, PTC and ANSYS introduced a solution that will quickly add the results of ANSYS engineering calculations to applications created on the PTC platform for ThingWorx's Industrial Internet of Things (IIoT) .The development of a connector between these two technical engineering platforms will allow customers to transfer and use the results of calculations when operating physical products in the real world of things. Also, it will allow them to match intelligent finite element models with physical products throughout the life cycle and manipulate the results in the real world. In practice, the software product ANSYS has been frequently used for the development of digital twins. For example, it can be used to create a model of an actually existing pump. That is what Chris MacDonald with his colleagues from ANSYS and PTC exactly did [5] .They believe: "Companies can use digital twins to detect and isolate faults, optimize asset operation, and generate insights to improve the next generation of the product. In 2012, NASA published the work "The digital twin paradigm for the future NASA and US Air Force vehicles" [6], which assessed : "Future generations of NASA and US Air Force vehicles will require lighter mass while being subjected to greater loads and more extreme days. Current approaches for certification will likely be unable to address these extreme requirements. To address the shortcomings of conventional approaches, a fundamental paradigm shift is needed. This paradigm shift, the Digital Twin ..." As mentioned above, the potential for application field of digital twins is quite wide. The digital twin finds its widest application in different technical systems. However, recently, digital twins appear in bioengineering and biomechanical systems [7,8]. In this article, some methods of numerical simulation of the human coronary vessel system for creating its digital twin are considered. For effective application of CFD-modeling one can use number of commercial programs, such as ANSYS, FlowVision, ComsolMultiphysics, etc., and open source programs such as OpenFOAM, Salome, Code Saturn, and others.
How can a digital twin of human vessels be useful? Taking into account the pulsation nature of the velocity inside the vessels, non-Newtonian properties of the blood, and the presence of turbulent flow, it becomes possible to predict the behavior of the object of investigation with some accuracy. Such data as the flow in the areas of narrowing, expansion, and obstacles inside the vessel, the dynamics of the zones of shear stress, etc. is extremely valuable in the study of the development of diseases of the human cardiovascular system. Damage to endothelial cells, given as a wall material, is caused, in particular by turbulent blood flow in bifurcations (Virchow's triad -can be the cause of thrombosis, an extremely dangerous disease). The identification of the turbulent flow core allows observing a direct correlation with the wear of valve flaps that prevent the reverse flow of blood (phlebeurysm, an extremely common disease). The data required for the analysis are an array of hemodynamic and mechanical characteristics. Obtaining this data is the purpose of the work presented. Whether it is necessary to develop a realistic virtual model of biological objects to reflect the real and digital worlds, overcoming the gap between design and production (for example, production of stent for removal of thrombus) or not, will depend on comprehensive research and improvement of various methods of numerical modeling [9].

Methods
One of the most important tasks in the development of a numerical model of the heart-vessel system is the maximally accurate adherence to the physical nature of the vessel. The validation criterion of a numerical model in relation to a physical process or object is the higher the more perfect the specific algorithm and the methodology for performing calculations. This was demonstrated in a study by Xinrong and co-authors, "Influence of the boundary conditions on aorta blood flow simulations" [10], as well as in the study "Oscillatory flow in human arteries", by Giannoglou and his colleagues [11]. In their works, the term "Computational Cardiovascular Dynamics" is used to show the possibility of obtaining high dependability of mathematical methods. This term combines Bioengineering and Computational Fluid Dynamics, as part of Computational Mechanics, as well as phenomena of hemodynamics and the theory of mechanics for the numerical analysis of cardiovascular systems. Thus, in order to apply the digital twin to diagnose and treat cardiovascular diseases, it is required to use a set of specific techniques from the cross disciplines. In the present paper, for these purposes, the various computational domains and numerical setups were used.
The objects of research are part of the system of coronary heart vessels, a two-dimensional vessel with fixed leaves in two positions, and with the presence of obstacle on the way of the blood flow.
For each computational domain are marked following elements: the blood flow inlet and outlet, as well as vessel walls. Creating of geometry is an extremely important stage of the work, which determines the level of correspondence of a particular algorithm with the actual physical process. The geometry is created according to the peculiarities of the human vessel system. Increasing the accuracy of the solution is accomplished by adapting the mesh in the Meshing module of ANSYS. For numerical simulation, a constant increase in the number of elementary cells is often an unjustified solution. Saving the processing power while keeping the quality of the solution is achieved by adapting the mesh in the context of the simulated process. For example, for the region of the computational domain in which there is a transition of energy from the ordered state (the absence of vortices) to the disordered, and then to the heat (dissipation near the solid wall), the Inflation function is used. The Inflation function allows scaling the grid near the wall in a logarithmic order. Such a grid structure allows turbulence models to perform calculations in the boundary layer. The presence or absence of vortices, their dimensions, the degree of dissipation, and the degree of influence of the geometry of the vessel on blood distribution are controlled by the solver settings and the grid structure.
The general view of the geometry and the computational grid is shown in Fig. 1  In the Numerical Setup of the Fluent, the mathematical description of the flow process is expressed by a system of differential equations consisting of continuity equations, the law of conservation of momentum, and the law of conservation of energy, generally described in Ansys Fluent 14.0: Theory Guide [6].
The task is solved in a non-stationary form (Time-Transient) in two-dimensional and three-dimensional computational domains. At the Inlet area, the main calculation parameter for the blood flow (Solver Typepressure-based) is defined. The velocity has a pulsating character [13].  Viscosity. Calculations of blood viscosity are conducted according to Hemorheology. Hemorheology studies the properties of the blood flow and its elements (plasma and cells) [14]. Laminar flow prevails in the vessels. The Reynolds number takes on a value in the region of 50 [15]. With this type of flow, the blood moves with cylindrical layers whose axis coincides with the axis of the vessel. In this case, the bloodstream is quasi-static. Therefore the balance between the inflow of energy from the mean motion and the transfer of energy to heat is maintained, the profile is uniform. However, with a decrease in blood viscosity, constrictions, dilations, bends, obstructions in the flowline, increased blood pressure, and when the inner surface of the vessel becomes coarse and rough, the laminar flow becomes turbulent [16][17]. Turbulent flow is characterized by the presence of vortices. These vortices significantly increase the internal friction of the blood. The volume velocity of the blood flow is no longer proportional to the pressure gradient, but to the square root of the last mentioned. To increase the volume flow velocity, for example, twice, it is necessary to increase the pressure gradient by about 4 times. This means that with a turbulent type of blood flow, the load on the heart increases significantly [18].

Results and discussion
In the calculations, both laminar and turbulent components are included when consider turbulence model equations. One of these components is the turbulent viscosity. It is the component in the equation, which ceases to be constant. Throughout its flow, its variation does not depend on the position in space, except for the wall layer [19][20]. For numerical simulation of processes in a turbulent core, wall layer and turbulence, turbulence models are used. The main ones are related to three types: DNS, RANS, LES.
In the presented study, there is insufficient density of the computational mesh for the use of DNS and LES. In addition, there is no effect of energy exchange between vortices and their pulsations. Of all the models of the RANS type, the universal group is k. In practice, k-ε models are the best for the widest range of Reynolds number.
However, case records of k-ε models have shown that in the calculation of boundary layers with a positive (negative) pressure gradient, all models of this type tend to overestimate the generation of the kinetic energy of turbulence. This fact leads to a fundamentally incorrect description of such flows [21][22]. This disadvantage is quite serious since it makes it practically impossible to correctly calculate the flows with a detachment from a smooth surface. For this reason, the calculation was carried out for two models: k-ε Realizable, and k-ω SST model that combining the main advantages of the groups k-ε Standart and k-ω Standart.
The Realizable model is based on k and ε, where k is responsible for the pulsation of the kinetic energy, and ε is responsible for how fast this energy dissipates into heat. Compared to Standard, it contains an alternative formulation for turbulent viscosity and modified transport equation for the dissipation rate obtained from the exact equation for the transfer of the mean-square vorticity fluctuation. The term "Realizable" means that the model satisfies certain mathematical limitations of the Reynolds stress, consistent with the physics of turbulent flows, in contrast to the k-ε Standard and RNG.
Model SST (Shear Stress Transport). This is a model of the group k-ω. Unlike Standard, it is stable in pure form and, unlike k-ε, underestimates the dissipation at the walls. This model copes well with the calculation of near-wall parameters. The determination of turbulent viscosity is modeled taking into account the transfer of turbulent shear stress [23]. The SST model is considered more accurate and reliable for a wider class of flows (for example, negative pressure gradient flows, aerodynamic profiles, transonic shock waves) than the Standard. The group of models k (both ε and ω) incorrectly reflect heat and mass transfer due to the fact that they are isotropic (properties do not change in the core). Anisotropic models are DNS and LES. The reasons for their incompatibility with the chosen methodology are mentioned above.
The functional of the used models used includes the connection of a variety of improved near-wall functions, which is useful in the numerical modeling of turbulent processes [24]. These models are a two-layer model with the so-called expanded functions of the wall. Their work is possible if the wall mesh is thin enough to be able to solve a viscous sublayer (typically with the first nearwall node placed at y + ≈1) [25]. In the chosen technique, this is meaningless, because the Near-Wall Treatment parameter is set to the Standard Wall Function.
As a fluid material, blood is used, the walls are epidermis. For the epidermis, the following parameters are entered: Density -1300 kg/m 3 , Specific Heat -3650 J/(kg·K), Thermal Conductivity -0.2 W/(m·K). For blood: Density -1060 kg/m 3 , Specific Heat -3750 J/(kg·K), Thermal Conductivity -0.478 W/(m·K) [27]. Blood considered as a non-Newtonian fluid. Hence, the calculation of its Molecular Viscosity cannot be performed in the classical way [28]. There are many models of non-Newtonian fluids, but as in the case of models of turbulence, the universal one has not been found. These are some of them: Carreau (1991), Power Law (1991), Non-Newtonian Power Law (1992), Casson (1993) [29]. They are based on the calculation of the viscosity through the Shear Rate, the temperature (in the study the temperature is constant), and the constants k and n. To specify non-Newtonian blood parameters for both flow modes, one need to use a special command: /define/models/viscous/turbulence-expert > turb nonnewtonian [30]. Then set the constants for viscosity calculation: Consistency Index, k -0.01669 (kg·s n-2 )/m, Power-Law Index, n -0.52815, Viscosity Limits -2.8894e 05 …100 kg/(m·s) [31,32]. In this case, the viscosity is calculated independently of the temperature, according to the Non-Newtonian Power Law (1992). Blood pressure is 13332 Pa. Number of iterations of hybrid initialization is set to 10. The size of the time step is 0.01 s, the number of steps is 50, the number of iterations per time step is 40. The maximum strain rate coincides with the peak of the velocity-time diagram. The presence of vortices in the regions of expanding, bifurcations and after bending is demonstrated. The models of turbulence k-ε Realizable and k-ω SST are equally good at modeling vortices. Fig.  6 shows the velocity profile in the cross-section of the vessel for the models k-ε Realizable and k-ω SST. It can be seen how the SST mode, unlike the Realizable model, slightly underestimates the velocity decrease near the wall (the transition is smoother). The Realizable model increases the velocity of the flow inside the flow core in comparison with the SST model. For both models, there is a general overestimation of the blood flow velocity. This phenomenon is caused by their isotropy. Both models neglect the flow properties in the core, applying averaging, instead of calculating.        7 shows Velocity Contours and shear stress plotted on the model of a two-dimensional vessel with fixed flaps, which prevent the reverse flow of blood. With a high-velocity gradient, the flap walls are subject to strong stress. Fig. 8 shows the formation of a vortex flow inside the vessel. Fig. 9 shows the distribution of shear stress on the flap of the vessel under the blood flow pressure.
The development of a realistic virtual model of the coronary vessel system is based on numerical modeling. During computational modeling, it is important to consider the model not only as an object but also as a process. The methods of mapping the physical process in the digital model are considered above. Turbulent blood flow in bifurcations refers to the causes of damage to the walls of the vessel. As a consequence, thrombus appears. The thrombus reduces the pore of the vessel, acts as an obstacle in the path of blood flow. This phenomenon also leads to the formation of a vortex flow inside the vessel. With a vortex type of blood flow, the load on the heart increases significantly, the blood pressure increases. Under the pressure of blood flow, the walls of the valves, which do not allow blood to flow in the opposite direction, wear down and eventually weaken. If a certain amount of blood is able to flow in the opposite direction, the veins undergo pathological changes. Stresses that cause physical pathologies are shown on a digital model. According to the mechanical characteristics of the coronary vessel model, a visual representation of the progression of heart disease has been obtained. Production, prevention, or timely treatment of diseases of the human cardiovascular system depends on the reliability of a digital twin (for example, production of a stent for the removal of a blood clot).