A calculation method for a coronary artery geometric multi-scale model based on fluid-structure interaction

Abstract

A calculation method for a fluid-structure interaction-based coronary artery geometric multi-scale model belongs to the hemodynamics domain. When performing hemodynamic simulations, there is a mutual influence between the vessel wall and the blood flow. To fully consider the influence of vessel wall elasticity on blood flow, a fluid-structure interaction model is generally used. However, when the fluid-structure interaction model simulates the non-invasive fractional flow reserve, it cannot simulate the change in microvascular resistance under the hyperemic state. To simultaneously consider the vessel wall and microvascular resistance, a geometric multi-scale calculation model based on fluid-structure interaction is established. The resistance R of the electronic component m and the inductor L are connected in series and coupled to the three-dimensional fluid-structure interaction model, where the resistance R m simulates the microcirculation resistance of the vessel branch, and the inductor L simulates the inertia of the microcirculation blood flow. This method can not only simulate the influence of the vessel wall on blood flow, but also provide reliable microvascular resistance boundary conditions for the three-dimensional fluid-structure interaction model, making it closer to the real physiological environment and improving the calculation accuracy of hemodynamics.

Description

Technical Field

[0001] The present invention provides a calculation method for a coronary artery geometric multi-scale model based on fluid-structure interaction, belonging to the field of numerical simulation of hemodynamics. Background Art

[0002] Fractional flow reserve. Fractional flow reserve (FFR) is the gold standard for clinically judging whether functional myocardial ischemia occurs, and is defined as: the ratio of the maximum blood flow in the myocardium supplied by the stenotic coronary artery to the maximum blood flow that the myocardium can provide assuming no stenosis. Based on the fact that myocardial blood flow is proportional to perfusion pressure, the calculation formula can be simplified to where Pa is the mean pressure at the aortic root when the coronary artery is in the maximum hyperemic state, and Pd is the mean pressure at 30 mm distal to the stenotic lesion, as Figure 1 shown.

[0003] With the rapid development of computational fluid dynamics and medical images, FFR can be calculated by numerical simulation methods. Based on coronary computed tomography angiography (CTA) images, computational fluid dynamics methods are used to simulate the hemodynamics in the coronary artery, obtain the pressure distribution under coronary hyperemia, and then calculate the ratio of the pressure at the distal end of the coronary stenosis to the aortic pressure, that is, the non-invasive numerical calculation of FFR (FFR CT ). Currently, a large number of researchers perform numerical calculations on FFR. For example, HeartFlowNXT showed at the Transcatheter Cardiovascular Therapeutics (TCT) conference that FFR CT is more accurate than CTA (79% vs. 34%). When calculating FFR using current existing scientific means and research methods, a fluid-structure interaction model or a geometric multi-scale model is usually adopted. When only using the fluid-structure interaction model for calculation, the change in microcirculation resistance after hyperemia cannot be fully considered. When using the traditional geometric multi-scale model for simulation calculation, the blood vessel wall is assumed to be rigid, but the real human blood vessels are elastic, and the blood vessel wall and blood interact in real time during blood flow. Therefore, only using the traditional geometric multi-scale model for numerical simulation will ignore the influence of the blood vessel wall on blood flow. Whether using the fluid-structure interaction model or the geometric multi-scale model for calculating FFR, the accuracy needs to be improved. Therefore, the present invention combines the two to establish a calculation method for a geometric multi-scale model based on fluid-structure interaction. This method can consider both the influence of the blood vessel wall on blood flow and the change in microvascular resistance under hyperemic conditions, is closer to the real blood vessel physiological environment, and provides a more accurate hemodynamic calculation method for fractional flow reserve. The calculation model is as Figure 2 shown.

[0004] Fluid-structure interaction model. When numerically simulating FFR, since the coronary artery wall is elastic and the blood flow pulsates with the heartbeat, the numerical calculation of hemodynamics requires solving the governing equations of three-dimensional unsteady flow in a deformable blood vessel. During the blood circulation process in the cardiac cycle, the arterial wall will deform and expand, and the blood flow region will also change accordingly. Therefore, the arterial blood flow and the blood vessel wall form a transient fluid-structure interaction mechanical system, and the Arbitrary Lagrange-Eulerian (ALE) method in continuum mechanics is needed to describe the motion and dynamic characteristics of the system. However, the general fluid-structure interaction simulation calculation method cannot provide the true physiological boundary conditions after coronary artery branching when calculating the fractional flow reserve, so a geometric multi-scale model needs to be established on the basis of fluid-structure interaction.

[0005] Geometric multi-scale model. Since the blood flow in the coronary microcirculation structure cannot be measured clinically in real time, continuously, and non-invasively, the true boundary conditions of coronary artery branches cannot be given, and a special strategy for simulating the blood circulation system is needed. The blood vessel distribution in the human body presents a complex network structure. Based on the similarity between the laws of the circuit system and the laws of blood flow, Frank established the Windkessel theory, and a lumped parameter model (zero-dimensional model) containing electrical components can be used to simulate the cardiovascular system. Equivalent the blood pressure and blood flow to voltage and current, which has the characteristics of simplicity and high efficiency. During the model establishment process, the flow rate is considered an important factor connecting the three-dimensional physiological parameters of the patient and the parameter values of each component of the model, and it is determined by the resistance value of the zero-dimensional model corresponding to the patient's personalized three-dimensional geometric model. The equivalent relationship between hemodynamic parameters and physical electrical parameters is described in Table 1.

[0006] Table 1 Equivalent relationship between hemodynamic parameters and physical electrical parameters

[0007] Tab.1 Equivalent relationship between hemodynamic parameters and physical electrical parameters

[0008]

[0009] The geometric multi-scale model uses the characteristics of models with different dimensions to simulate different parts of the circulatory system respectively. The three-dimensional model is used to simulate the hemodynamic environment of local details, while the blood flow and blood pressure in the peripheral circulatory system are simulated by the reduced-dimensional one-dimensional or zero-dimensional model, providing reliable boundary conditions for the calculation of the three-dimensional model.

[0010] When focusing on the details of 3D fluid-solid interactions without using artificially determined fixed boundary conditions, a 0D / 3D coupling model is often used. This model typically uses a 3D model to simulate the local flow field and solid deformation of interest, while using a 0D model to simulate the peripheral circulatory system. In this way, when the structure of the 3D model changes, the boundary conditions provided by the peripheral 0D model for the 3D model will also adaptively change accordingly, thus avoiding the adverse effects brought by fixed boundary conditions. The specific structural schematic diagram is as shown in Figure 3 shown. Summary of the Invention

[0011] When calculating the fractional flow reserve (FFR), most existing studies usually use a fluid-structure interaction model or a geometric multi-scale model for simulation calculations. When using a fluid-structure interaction model for calculation, the change in microcirculation resistance after hyperemia cannot be fully considered, and a geometric multi-scale model cannot provide reliable boundary conditions for it. When using a traditional geometric multi-scale model for simulation calculations, the blood vessel wall is defaulted to be rigid, but real human blood vessels are elastic, and the blood vessel wall and blood interact in real time during blood flow. Therefore, only using a traditional geometric multi-scale model for numerical simulation will ignore the influence of the blood vessel wall on blood flow. Therefore, in order to more realistically and accurately simulate the calculation of FFR, the present invention combines the two to establish a geometric multi-scale model calculation method based on fluid-structure interaction. This method can not only consider the influence of the blood vessel wall on blood flow but also take into account the change in microvascular resistance under hyperemic conditions, being closer to the real vascular physiological environment and providing more accurate hemodynamic calculations for the fractional flow reserve.

[0012] In order to implement the geometric multi-scale model calculation method based on fluid-structure interaction, the present invention is achieved through the following techniques.

[0013] (1) 0D / 3D coupling algorithm. The pressure obtained by the 0D model calculation provides the resistance boundary condition for the 3D model calculation, while the outlet flow obtained by the 3D model calculation serves as the inlet flow for the 0D model calculation. After flow distribution calculation, the microcirculation resistance value R m after coronary artery stenosis branches at rest is obtained. Since the fractional flow reserve is calculated under hyperemic conditions, the microcirculation resistance will become 0.24 times that at rest.

[0014] According to the theoretical formula:

[0015] Pa = Δp + Pd (1)

[0016] Pd = Q outlet *R m *0.24 (2)

[0017] Using quadratic interpolation, calculate the pressure drop of the stenosis simulation calculation model at the current flow rate. When the pressure drop of the stenosis simulation calculation model makes Equation 1 hold, it can be considered that the calculation balance is achieved. Where Pa is the pressure inlet boundary condition, giving the true physiological pressure waveform, R m is the microvascular resistance, Δp is the pressure loss generated by flowing through the stenosis, and Q outlet is the flow rate at the outlet of the stenosis model.

[0018] In the coupling algorithm, it is a transient calculation. A 0D-3D data exchange and a two-way fluid-solid data exchange are performed at each time step, and at the same time, residual detection is carried out. Define the error of the outlet pressure of the 3D model fluid domain, the inlet flow rate, the displacement of the solid domain, and the data transmission at the fluid-solid interface during different cardiac cycles as the residual detection item. When the residual is less than the preset value of 1e-6, it is considered that the calculation result converges and the simulation ends. The specific process is as Figure 4 shown.

[0019] (2) Geometric multi-scale fluid-structure interaction. The entire calculation platform is based on ANSYS-Workbench. ANSYS-Fluent is used for the fluid, and Transient Structural is used for the solid part. Two-way fluid-structure interaction simulation calculations are performed in System Coupling. ANSYS-Fluent is a relatively complete computational fluid dynamics simulation software owned by ANSYS. Based on ANSYS-Fluent, secondary development of ANSYS-Fluent is required to implement the 0D / 3D coupling algorithm for the simulation calculation of the 0D / 3D coupling model.

[0020] User Defined Functions (UDF) is a secondary development interface provided by Fluent. It can be dynamically connected to the Fluent solver to improve the solver performance and embed user-defined calculation methods. Standard C language library functions can be used in UDF, and predefined macros provided by Fluent Inc. can also be used.

[0021] The present invention uses a user-defined program based on C language to enhance the boundary conditions, uses a large number of macro definitions to achieve the interaction of data between the user and the solver, and enhances the functions of Fluent and the application scope of the model. When the program is used, it is regarded as a compiled function, embedded in the shared library and connected to Fluent, and added to the three-dimensional model outlet boundary condition setting by loading the compiled subroutine.

[0022] A calculation method for a geometric multi-scale model based on fluid-structure interaction, the technical solution includes the following steps:

[0023] (1) Establish a fluid domain of an ideal coronary artery and a simulation model of the solid structure of the blood vessel wall. After assembling the two, a three-dimensional fluid-structure interaction model is obtained;

[0024] (2) Connect an inductor L and a resistor R in series at the outlet of the established three-dimensional fluid-structure interaction model m , R m The resistance value is the microcirculation resistance value determined by the process of quantifying the relationship between the flow rate, pressure and flow resistance of a section of blood vessel by combining Kirchhoff's law and the Windkessel theory established by Frank. The inductor is used to adjust the data transfer process between the coupling models, control the convergence speed and improve the convergence condition. The inductance value of the present invention is 0.1H.

[0025] (3) Set the boundary conditions and the properties of the solid material. The inlet boundary condition is given the real physiological pressure waveform of the patient. The outlet boundary uses UDF to implement the 0D / 3D coupling algorithm, and extracts the outlet flow rate Q calculated by the 3D model outlet as the inlet flow rate for the 0D model calculation. After the flow distribution calculation, the microcirculation resistance value R m is obtained. Since the fractional flow reserve is calculated under the hyperemic state, the microcirculation resistance will become 0.24 times that under the resting state. After writing the entire microcirculation resistance with UDF, it is compiled and loaded into the Fluent outlet boundary condition as the outlet resistance boundary of the three-dimensional model. Set the material properties of the blood vessel wall in the solid domain, and determine the Young's modulus and Poisson's ratio of the blood vessel wall by referring to the literature.

[0026] (4) Use the finite element simulation calculation method to calculate the above model. According to the theoretical formula:

[0027] Pa = Δp + Pd (1)

[0028] Pd = Q outlet * R m * 0.24 (2)

[0029] Use quadratic interpolation to calculate the pressure drop of the stenosis calculation simulation model at the current blood flow rate. When the pressure drop of the stenosis calculation simulation model makes formula (1) hold, it can be considered that the fluid-structure interaction calculation for one time step reaches equilibrium. Where Pa is the pressure inlet boundary condition, giving the real physiological pressure waveform, R m is the microvascular resistance, Δp is the pressure loss generated by flowing through the stenosis, and Q outlet is the flow rate at the outlet of the stenosis calculation simulation model.

[0030] In the coupling algorithm, transient calculations are performed. During each iterative calculation, a 0D-3D data exchange and a two-way fluid-solid data exchange are carried out, and residual detection is also performed simultaneously. Define the outlet pressure of the fluid domain of the 3D model, the inlet flow rate, the displacement of the solid domain, and the error in data transmission at the fluid-solid interface during different cardiac cycles as the residual detection term. When the residual is less than the preset value of 1e-6, the calculation result is considered to converge, and the simulation ends. Then, take the hemodynamic parameters of the blood vessel, including the pressures at the inlet, outlet, and wall, the wall shear stress (WSS), the deformation displacement of the blood vessel wall, and the blood flow velocity, and calculate the non-invasive fractional flow reserve.

[0031] This method not only considers the influence of the blood vessel wall on blood flow but also takes into account the changes in microvascular resistance under hyperemic conditions, making it closer to the real blood vessel physiological environment and providing more accurate hemodynamic calculations for fractional flow reserve. Brief Description of the Drawings

[0032] Figure 1 : Principle of fractional flow reserve calculation.

[0033] Figure 2 : Geometric multi-scale model based on fluid-structure interaction

[0034] Figure 3 : Geometric multi-scale model

[0035] Figure 4 : Algorithm flow chart Detailed Implementation Modes

[0036] The present invention will be explained below in conjunction with the detailed implementation modes, but the present invention is not limited to the following embodiments.

[0037] Embodiment 1

[0038] A calculation method for a geometric multi-scale model based on fluid-structure interaction includes the following steps:

[0039] (1) Based on the real coronary artery anatomical structure parameters, construct an ideal coronary artery blood vessel calculation simulation model. According to clinical statistics, it is found that coronary artery stenosis is prone to occur in the left anterior descending branch. The present invention creates an ideal left anterior descending coronary artery blood vessel simulation calculation model with a diameter of 4 mm, a stenosis rate of 40%, a blood vessel wall thickness of 0.5 mm, and a blood vessel length of 100 mm;

[0040] (2) Determine the values of the inductance L and resistance R connected in series at the outlet of the three-dimensional fluid-structure interaction model m of R m The resistance value is the microcirculation resistance value determined by quantifying the relationship between the flow rate, pressure, and flow resistance of the blood vessel by combining Kirchhoff's law and the Windkessel theory established by Frank. The inductance is used to regulate the data transfer process between the coupling models, control the convergence speed, and improve the convergence condition. The inductance value in the present invention is 0.1 H.

[0041] (3) Build a fluid-structure interaction calculation platform and perform two-way fluid-structure interaction simulation calculations using Fluent + Transient Structural + System Coupling in ANSYS-Workbench. Import the geometric model, suppress the blood vessel wall part in Fluent, mesh the fluid model, and define the fluid density as 1050 kg / m 3 , viscosity 0.0035 Pa·s, set the boundary conditions, assign the patient's real physiological pressure waveform to the inlet boundary condition of the three-dimensional model, and the outlet boundary condition is provided by the program written in UDF. Extract the outlet flow rate Q obtained from the 3D model calculation outlet as the inlet flow rate for the 0D model calculation. The microcirculation resistance value R after coronary artery stenosis bifurcation under resting state m is the microcirculation resistance value determined by the process of combining Kirchhoff's law and the Windkessel theory established by Frank to quantify the relationship between the flow rate, pressure, and flow resistance of a section of blood vessel. Since the fractional flow reserve is calculated under the hyperemic state, the microcirculation resistance will become 0.24 times that under the resting state. In order to regulate the data transfer process between the coupling models, control the convergence speed, and improve the convergence condition, an inductor L = 0.1 H is added as a section of damping and connected in series to the microcirculation resistance. After the entire microcirculation resistance is written in UDF, it is loaded into the Fluent outlet boundary condition in a compiled manner as the resistance boundary at the outlet of the three-dimensional model. Define the outer surface of the fluid as the data transfer surface with the solid, set the time step to 0.005 s, and the calculation duration to 1.6 s;

[0042] (4) Suppress the fluid part under the Transient Structural module, mesh the blood vessel wall, define the blood vessel wall material as linear elastic with Young's modulus 1 Mpa and Poisson's ratio 0.45, set the inner side of the blood vessel wall as the surface for transferring with the fluid part, set the time step to 0.005 s, and the calculation duration to 1.6 s;

[0043] (5) Couple the fluid and solid data transfer surfaces in the System Coupling module, and set the total duration of the entire fluid-structure interaction calculation to 1.6 s and the time step to 0.005 s;

[0044] (6) Extract the calculation results including the pressures at the inlet, outlet, and wall, as well as the wall shear stress (WSS), blood vessel wall deformation displacement, and blood flow velocity, and calculate the non-invasive fractional flow reserve. The FFR value obtained in this embodiment is 0.94.

Claims

1. A calculation method for a geometric multi-scale model based on fluid-structure interaction, characterized in that, It includes the following steps: (1) Based on the real coronary artery anatomical structure parameters, construct a fluid-structure interaction three-dimensional simulation calculation model of the blood fluid domain of the ideal coronary artery and the vascular wall solid structure; (2) Establish a multi-scale simulation model based on fluid-structure coupling to achieve 0D / 3D coupling. An inductor L and a resistor R are connected in series at the outlet of the established three-dimensional fluid-structure coupling model. m , determine the inductance L and resistance R connected in series at the outlet of the three-dimensional fluid-structure interaction model m The value of (3) Build a fluid-structure interaction calculation platform, and set the boundary conditions and the material properties of the three-dimensional model on the calculation platform; Impose the real physiological pressure waveform of the patient at the inlet of the three-dimensional model, and load the microcirculation resistance boundary condition at the outlet; (4) Use the finite element simulation calculation method to calculate the above model formula, and extract the calculation results including the pressures at the inlet, outlet and wall surface, the wall shear stress, the deformation displacement of the vascular wall, the blood flow velocity, and calculate the non-invasive fractional flow reserve; Step (2) includes: R m The resistance value is the microcirculation resistance value determined by the process of quantifying the relationship among the flow rate, pressure and flow resistance of blood vessels by combining Kirchhoff's law and the elastic cavity theory established by Frank; the inductance is used to regulate the data transfer process between coupling models, control the convergence speed and improve the convergence condition; the inductance value adopted is 0.1H; Step (3) includes: building a calculation platform using Fluent + Transient Structural + System Coupling in ANSYS-Workbench to perform two-way fluid-structure interaction simulation calculations; importing the geometric model, suppressing the blood vessel wall part in Fluent, meshing the fluid model, defining the fluid density as 1050 kg / m 3 , viscosity as 0.0035 Pa·s, setting boundary conditions, assigning the real physiological pressure waveform to the inlet boundary condition of the three-dimensional model, and providing the outlet boundary condition by the program written with UDF, and extracting the outlet flow rate Q obtained from the 3D model calculation outlet as the inlet flow rate for the 0D model calculation; the microcirculation resistance value R after coronary artery stenosis branching under the resting state m is the microcirculation resistance value determined by the process of combining Kirchhoff's law and the Windkessel theory established by Frank to quantify the relationship between the flow rate, pressure, and flow resistance of a section of blood vessel; the microcirculation resistance will become 0.24 times that under the resting state, adding an inductor L = 0.1 H as a section of damping, connecting it in series to the microcirculation resistance, after writing the entire microcirculation resistance with UDF, loading it into the Fluent outlet boundary condition by compilation as the resistance boundary at the outlet of the three-dimensional model, defining the outer surface of the fluid as the data transfer surface with the solid, setting the time step as 0.005 s, and the calculation duration as 1.6 s; determining the blood density as 1050 kg / m 3 and viscosity as 0.0035 Pa·s; suppressing the fluid part under the Transient Structural module, meshing the blood vessel wall, defining the blood vessel wall material as linearly elastic and isotropic, with Young's modulus of 1 Mpa and Poisson's ratio of 0.45, setting the inner side of the blood vessel wall as the surface for fluid part transmission, setting the time step as 0.005 s, and the calculation duration as 1.6 s; coupling the fluid and solid data transfer surfaces in the System Coupling module, setting the total duration of the entire fluid-structure interaction calculation as 1.6 s and the time step as 0.005 s.

2. The calculation method of a geometric multi-scale model based on fluid-structure interaction according to claim 1, characterized in that In step (1), based on the real coronary artery anatomical structure parameters, an ideal coronary artery three-dimensional vascular calculation simulation model is constructed, and an ideal left anterior descending artery simulation calculation model with a diameter of 4 mm, a moderate stenosis rate of 40%-70%, a vascular wall thickness of 0.5 mm, a total vascular length of 100 mm, a stenosis inlet length of 30 mm, a stenosis length of 20 mm, and a stenosis outlet length of 50 mm is created.

3. A calculation method of a geometric multi-scale model based on fluid-structure interaction according to claim 1, characterized in that, (4) includes: using the finite element simulation method to calculate the model in step (3), and after calculation, extract the hemodynamic parameters of the blood vessel. The hemodynamic parameters include the pressures at the inlet, outlet and fluid wall surface, the wall shear stress, the deformation displacement of the vascular wall, the blood flow velocity, and calculate the non-invasive fractional flow reserve.

Images