Method and system for simulating hemodynamics of moyamoya disease with non-newtonian fluid

Abstract

The present application relates to the field of biomedical engineering, and more particularly to a method and system for simulating hemodynamics of moyamoya disease by fusing non-Newtonian fluid. The method comprises the following steps: reconstructing a three-dimensional geometric model of a target cranial and cerebral vascular network from medical image data of a moyamoya disease patient, dividing the model into a large vessel calculation domain and a small vessel calculation domain covering moyamoya vessels according to the diameter characteristic parameters, respectively configuring a Newtonian fluid constitutive equation with constant viscosity and a non-Newtonian fluid constitutive equation with nonlinear viscosity function, setting an outlet impedance boundary condition, fusing the fluid and the solid dynamics equation of the vessel wall to construct a fluid-structure interaction mathematical physics model, and iteratively outputting three-dimensional dynamic blood flow parameters by using a parallel solving algorithm. The present application improves the distortion problem of the traditional simulation which excludes small vessels and assumes a single fluid, realizes high-fidelity joint simulation of the complex flow behavior of macro-arteries and micro-moyamoya networks, and restores the dynamic interaction physical process of blood flow and the vessel wall.

Description

Technical Field

[0001] This invention relates to the field of biomedical engineering, and in particular to a method and system for simulating hemodynamics in Moyamoya disease by incorporating non-Newtonian fluids. Background Technology

[0002] Moyamoya disease is a serious cerebrovascular disease characterized by progressive stenosis or occlusion of large intracranial arteries and the development of an abnormal network of microcollateral vessels. Superficial temporal artery-middle cerebral artery bypass surgery is currently the main intervention to improve cerebral blood flow perfusion in patients. Postoperatively, there is intense competition between intracranial and extracranial blood flow, and the bypass vessels undergo morphological remodeling and evolution under complex stress. Since the dynamic interaction between postoperative global cerebral blood flow distribution and vascular morphology directly affects whether patients develop hyperperfusion syndrome or long-term ischemic recurrence, establishing a computational model that can accurately quantify and dynamically simulate the complex hemodynamic behavior of the entire brain is crucial for personalized clinical treatment and long-term prognostic assessment.

[0003] Existing traditional medical imaging techniques are unable to quantitatively capture complex transient hemodynamic changes. Conventional numerical simulation methods generally use simplified models, directly eliminating small, fumed blood vessels and applying only a single Newtonian fluid assumption, ignoring the non-Newtonian nonlinear rheological properties of blood in microvessels. At the same time, they lack strong coupling considerations for the elastic deformation of the vessel wall and the physiological impedance at the end of the vessel, resulting in distorted simulation data that is difficult to meet the needs of accurate prediction. Summary of the Invention

[0004] To overcome the above shortcomings, this invention provides a method and system for simulating hemodynamics in Moyamoya disease by incorporating non-Newtonian fluids. This aims to improve the problem of distorted blood flow assessment caused by traditional simulations that remove tiny moyamoya blood vessels and ignore non-Newtonian characteristics.

[0005] In a first aspect, the present invention provides the following technical solution, which integrates a non-Newtonian fluid hemodynamic simulation method for Moyamoya disease, including: S1. Obtain brain medical imaging data of patients with moyamoya disease, extract vascular branch structures including intracranial and extracranial major arteries, bypass vessels and moyamoya vessels, construct the target brain vascular network, and reconstruct the three-dimensional geometric model of the target brain vascular network. S2. Extract the diameter feature parameters of each blood vessel branch structure in the three-dimensional geometric model, and divide the three-dimensional geometric model into a large blood vessel calculation domain and a small blood vessel calculation domain covering the smoke blood vessels based on the diameter feature parameters. S3. Configure a Newtonian fluid constitutive equation with a constant viscosity parameter in the large blood vessel computational domain, and configure a non-Newtonian fluid constitutive equation with a nonlinear function of local flow field variables as the viscosity parameter in the small blood vessel computational domain. S4. Set the outlet impedance boundary conditions based on the blood flow physiological mechanism matching the three-dimensional geometric model, and integrate the Newtonian fluid constitutive equation, the non-Newtonian fluid constitutive equation and the blood vessel wall solid dynamics equation to construct a fluid-structure coupling mathematical and physical model for the target cranial vascular network. S5. Using the outlet impedance boundary condition as an input constraint, the fluid-structure interaction mathematical-physical model is coupled and iteratively solved using a computational fluid dynamics parallel solution algorithm, and the three-dimensional dynamic hemodynamic parameters of the target cranial vascular network are output.

[0006] Preferably, in step S1, the step of reconstructing the three-dimensional geometric model of the target cranial vascular network includes: Magnetic resonance angiography image data is extracted from the acquired cranial medical imaging data, and image segmentation processing is performed on it to obtain the pixel matrix of the target region; Spatial topology reconstruction is performed on the pixel matrix to generate a three-dimensional point cloud contour containing the vascular branch structure of the intracranial and extracranial aorta, bypass vessels, and smoke vessels. Perform surface smoothing and mesh generation on the three-dimensional point cloud contour to output the three-dimensional geometric model.

[0007] Preferably, in step S2, the step of dividing the three-dimensional geometric model into a large blood vessel computational domain and a small blood vessel computational domain covering the smoke-filled blood vessels based on the pipe diameter characteristic parameters includes: Extract the centerline of each vascular branch structure in the three-dimensional geometric model, and calculate the equivalent diameter of its cross section as the tube diameter characteristic parameter; The tube diameter feature parameters are compared with a preset blood vessel scale threshold, and spatial grids that are greater than or equal to the scale threshold and spatial grids that are less than the scale threshold are marked respectively. Spatial grids larger than or equal to the scale threshold are divided into the large blood vessel computation domain, and spatial grids smaller than the scale threshold are divided into the small blood vessel computation domain.

[0008] Preferably, in step S3, the step of configuring a Newtonian fluid constitutive equation with a constant viscosity parameter into the large blood vessel computational domain includes: The fluid control equations are defined within the computational domain of the large blood vessel, and the strain rate tensor of the flow field is extracted. The stress tensor in the computational domain of the large blood vessel is set to have a linear mapping relationship with the strain rate tensor; A fixed dynamic viscosity value is assigned to the linear mapping relationship to generate a Newtonian fluid constitutive equation with a constant viscosity parameter.

[0009] Preferably, in step S3, the step of configuring the non-Newtonian fluid constitutive equation with viscosity parameters as nonlinear functions of local flow field variables in the computational domain of the small blood vessels includes: Extract the flow field deformation parameters within the computational domain of the tiny blood vessels, and calculate the second invariant of the flow field strain rate tensor; Extract the preset blood rheology consistency index and shear cutoff index, and construct a piecewise power-law nonlinear viscosity function with the second invariant as the independent variable; By embedding the piecewise power-law nonlinear viscosity function into the stress tensor calculation, a non-Newtonian fluid constitutive equation is constructed in which the viscosity parameter is a nonlinear function of the local flow field variable.

[0010] Preferably, in step S4, the step of setting the outlet impedance boundary conditions based on the blood flow physiological mechanism matching the three-dimensional geometric model includes: Identify the geometric outlet of the corresponding end microcirculation bed in the three-dimensional geometric model and configure the lumped parameter model; Based on the aforementioned blood flow physiological mechanism, the dynamic resistance value and compliance parameter at each of the geometric outlets are calculated; The dynamic resistance value and compliance parameter are mapped to a time-varying pressure constraint term to generate the outlet impedance boundary condition.

[0011] Preferably, in step S4, the step of fusing the Newtonian fluid constitutive equation, the non-Newtonian fluid constitutive equation, and the blood vessel wall solid dynamics equation includes: A Saint-Venant-Kirchhoff hyperelastic material model incorporating the Lamé coefficient is introduced to construct the solid dynamic equations of the blood vessel wall that characterize the deformation features of the blood vessel wall under stress. A continuity constraint that satisfies motion coordination and interfacial stress balance is established at the interface between the internal fluid and the solid blood vessel wall, and a dynamic mesh equation is introduced to track the boundary deformation of the computational domain. By using the continuity constraint, the state variables of the Newtonian fluid constitutive equation and the non-Newtonian fluid constitutive equation, as well as the displacement variables of the blood vessel wall solid dynamics equation, are integrated into a fully implicit system matrix, thus completing the construction of the fluid-structure interaction mathematical physics model.

[0012] Preferably, in step S5, the step of using a computational fluid dynamics parallel solution algorithm to perform coupled iterative solution of the fluid-structure interaction mathematical-physical model includes: The fluid-structure interaction mathematical physics model is mapped to a parallel computing architecture, and the computational domain is spatially partitioned using a domain decomposition strategy. A nonlinear solver is introduced to perform strongly coupled iterative calculations on the partitioned fluid-structure interaction mathematical-physical model. The iterative residuals during the solution process are monitored, and the three-dimensional dynamic hemodynamic parameters are extracted and output when the convergence condition is met.

[0013] Secondly, the present invention provides the following technical solution: a hemodynamic simulation system for Moyamoya disease that integrates non-Newtonian fluids, comprising: The model reconstruction module is used to acquire brain medical imaging data of patients with moyamoya disease, extract vascular branch structures including intracranial and extracranial major arteries, bypass vessels and moyamoya vessels, construct the target brain vascular network, and reconstruct the three-dimensional geometric model of the target brain vascular network. The computational domain partitioning module is used to extract the diameter feature parameters of each blood vessel branch structure in the three-dimensional geometric model, and divide the three-dimensional geometric model into a large blood vessel computational domain and a small blood vessel computational domain covering the smoke blood vessels based on the diameter feature parameters. The constitutive equation configuration module is used to configure a Newtonian fluid constitutive equation with a constant viscosity parameter for the large blood vessel computational domain, and at the same time configure a non-Newtonian fluid constitutive equation with a viscosity parameter that is a nonlinear function of the local flow field variable for the small blood vessel computational domain. The fluid-structure interaction modeling module is used to set the outlet impedance boundary conditions based on the blood flow physiological mechanism matching the three-dimensional geometric model, and to integrate the Newtonian fluid constitutive equation, the non-Newtonian fluid constitutive equation and the vascular wall solid dynamics equation to construct a fluid-structure interaction mathematical and physical model for the target cranial vascular network. The simulation solution output module is used to take the outlet impedance boundary condition as input constraint, and use a computational fluid dynamics parallel solution algorithm to perform coupled iterative solution on the fluid-structure interaction mathematical and physical model, and output the three-dimensional dynamic hemodynamic parameters of the target cranial vascular network.

[0014] The present invention has the following beneficial effects: 1. In this invention, the whole cranial vascular model is divided into computational domains of different scales by extracting the diameter features. The constitutive equations of Newtonian fluid with constant viscosity and non-Newtonian fluid with piecewise power-law nonlinearity are integrated in the same system. This improves the distortion problem caused by the traditional single fluid assumption or the exclusion of microvessels, and realizes high-fidelity joint simulation of the complex rheological behavior of macroscopic large artery blood flow and microscopic abnormal capillary network in Moyamoya disease.

[0015] 2. In this invention, the lumped parameter outlet impedance boundary condition is configured according to the blood flow regulation mechanism, and the solid dynamic equation of the blood vessel wall based on the hyperelastic material model and the dynamic mesh constraint are deeply coupled. The complex blood flow and the stress deformation of the vessel wall are unified and integrated into the fully implicit system algebraic matrix, which restores the physiological buffering effect of the terminal microcirculation and the dynamic nonlinear fluid-structure interaction physical process of blood vessels under the impact of pulsating blood flow.

[0016] 3. In this invention, for the extremely complex large-scale fluid-structure interaction system of the whole brain, a domain decomposition strategy is adopted to divide the space and map it to a parallel architecture. A Newton-Krylov-Schwarz nonlinear solver with an additive Schwarz preprocessing mechanism is introduced to perform strongly coupled iteration, which effectively optimizes the condition number of the large sparse equation system and achieves efficient and stable solution of the large-scale nodal numerical model while ensuring the convergence accuracy of the multiphysics field. Attached Figure Description

[0017] Figure 1 The flowchart is a simulation method for hemodynamics of moyamoya disease that incorporates non-Newtonian fluids, as proposed in this invention. Figure 2 This is a flowchart illustrating the basic steps of the hemodynamic numerical simulation proposed in this invention. Figure 3 This is a schematic diagram of the hemodynamic mathematical modeling proposed in this invention; Figure 4 This is an architecture diagram of the Moyamoya disease hemodynamic simulation system that integrates non-Newtonian fluids proposed in this invention. Detailed Implementation

[0018] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0019] Example 1 In a first embodiment of the present invention, the present invention provides a method for simulating the hemodynamics of moyamoya disease by incorporating non-Newtonian fluids, such as... Figure 1 As shown, it includes: S1. Obtain brain medical imaging data of patients with moyamoya disease, extract vascular branch structures including intracranial and extracranial major arteries, bypass vessels and moyamoya vessels, construct the target brain vascular network, and reconstruct the three-dimensional geometric model of the target brain vascular network. Further, step S1, the step of reconstructing the three-dimensional geometric model of the target cranial vascular network, includes: Magnetic resonance angiography image data is extracted from the acquired cranial medical imaging data, and image segmentation processing is performed on it to obtain the pixel matrix of the target region; Spatial topology reconstruction is performed on the pixel matrix to generate a three-dimensional point cloud contour containing the vascular branching structures of intracranial and extracranial aortas, bypass vessels, and moss vessels. Perform surface smoothing and mesh generation on the 3D point cloud contour to output a 3D geometric model.

[0020] Specifically, high-resolution time-of-flight magnetic resonance angiography (TFAA) images of patients with Moyamoya disease were acquired. This data includes the gray-scale spatial distribution characterizing blood flow signal intensity. To extract vascular branch structures including intracranial and extracranial aortas, bypass vessels, and Moyamoya vessels, a threshold segmentation algorithm was used to perform image segmentation processing on the acquired image data. The input three-dimensional image data spatial domain is assumed to be a three-dimensional matrix discretized from a continuous function. Indicates spatial coordinates as , , The grayscale intensity value at that location. Construct a binary pixel matrix for the target region. The discriminant formula is as follows: ; in A preset signal segmentation threshold is used to incorporate blood imaging characteristics. In this scheme, the signal segmentation threshold... The preferred value is 0.35 of the normalized grayscale intensity value. When the grayscale intensity value is greater than or equal to the set threshold, the pixel matrix value of the corresponding coordinate is recorded as 1, indicating that the location belongs to the vascular blood flow region; otherwise, it is recorded as 0, belonging to non-vascular background brain tissue. The obtained target region pixel matrix is ​​composed of the set of coordinate points with a value of 1.

[0021] Subsequently, spatial topological reconstruction is performed on the pixel matrix to generate a 3D point cloud contour of the vascular branching structure. In this step, the moving cube algorithm is used to extract isosurfaces from the discrete voxel field. Treating the pixel matrix as a scalar field, the algorithm traverses a voxel grid composed of eight adjacent pixels, calculating the boundary intersections where the scalar field equals a specific isosurface through linear interpolation. Connecting these intersections forms multiple triangular patches, which are then combined and stitched together to form a spatial surface that approximates the boundaries of the actual vascular anatomy, generating a 3D point cloud contour containing intracranial and extracranial aortas, bypass vessels, and moyamoya vessels. ,in Indicates the first The three-dimensional spatial coordinate vectors of the vertices of the contour. This represents the total number of contour vertices.

[0022] Because the directly extracted 3D point cloud contours exhibit step discontinuities at pixel boundaries due to image discretization, surface smoothing and mesh generation are necessary. The Laplacian smoothing algorithm is used to adjust the positions of the vertices on the point cloud contour surface, applying this smoothing to any vertex on the contour. The update position after a single smooth iteration The calculation formula is as follows: ; in Represents the vertex The set of directly connected topologically adjacent vertices. A weighted coefficient is used to measure the influence of adjacent vertices, and its value is the reciprocal of the total number of adjacent vertices. The relaxation factor is used to control the amount of smooth displacement in a single step. In this scheme, the relaxation factor... The optimal value is 0.5. After multiple iterative calculations, the vascular contour transition is smoothed. Finally, a volume mesh generation algorithm based on Delaunay triangulation is used to fill the internal space of the smoothed spatial boundary, discretizing the connected complex vascular lumen into mutually adjacent finite tetrahedral mesh elements, and outputting a three-dimensional geometric model of the target cranial vascular network.

[0023] The above steps achieve a precise conversion from the original grayscale medical image to the finite element computational mesh, eliminating image staircase errors while completely preserving the geometric topological connectivity of the tiny smoke-like blood vessels.

[0024] S2. Extract the diameter feature parameters of each blood vessel branch structure in the three-dimensional geometric model, and divide the three-dimensional geometric model into a large blood vessel computation domain and a small blood vessel computation domain covering the smoke blood vessel based on the diameter feature parameters. Further, in step S2, the step of dividing the three-dimensional geometric model into a large vessel computational domain and a small vessel computational domain covering the smoke-filled vessels based on the pipe diameter characteristic parameters includes: Extract the centerline of each vascular branch structure in the three-dimensional geometric model and calculate the equivalent diameter of its cross section as the tube diameter characteristic parameter; The tube diameter feature parameters are compared with the preset blood vessel scale threshold, and the spatial grids that are greater than or equal to the scale threshold and the spatial grids that are less than the scale threshold are marked respectively. Spatial grids with a scale threshold or greater than the scale threshold are divided into large vessel computational domains, and spatial grids with a scale threshold less than the scale threshold are divided into small vessel computational domains.

[0025] Specifically, the boundary surface mesh data of the 3D geometric model is acquired, and a spatial Euclidean distance field is constructed from the boundary surface into the vascular network. The topological skeleton is extracted along the ridges of the local maxima of the Euclidean distance field to obtain the set of centerline nodes for each vascular branch structure. For any node position on the centerline, an orthogonal normal plane perpendicular to the tangent vector of the current local centerline is constructed. The closed two-dimensional contour formed by the intersection of this orthogonal normal plane and the boundary surface of the 3D geometric model is calculated, and the true vascular cross-sectional area inside this closed two-dimensional contour is obtained by integrating on the plane using Green's theorem. Let the vascular cross-sectional area at any node on the centerline be... The equivalent diameter of the corresponding cross section is calculated based on the principle of equivalent circular area. The calculation formula is as follows: ; in Pi is a constant. The equivalent diameter is obtained by calculating the centerline node set point by point, and this equivalent diameter is used as a characteristic parameter describing the local morphological geometry of the blood vessel.

[0026] Next, domain attribute determination is performed on the tetrahedral spatial mesh elements generated by discretization within the 3D geometric model. The centroid coordinates of any spatial mesh element are extracted, and the nearest node to that centroid is searched in the set of central axis nodes. The equivalent diameter carried by this nearest node is assigned to the current spatial mesh element, ensuring that each computational mesh in the whole-brain 3D vascular network possesses a local scale attribute representing the thickness of its location. Let the diameter characteristic parameter inherited by a certain spatial mesh element be... A preset vascular scale threshold is introduced based on the structural differences in moyamoya disease imaging. In this scheme, based on the critical change scale of blood rheological properties in hemodynamics, the vascular scale threshold is used. The preferred size is 0.5 mm. A computational domain partitioning determination function is constructed for all tetrahedral spatial mesh elements. The logical classification formula for this function is as follows: ; When the output of the region partitioning decision function for a certain spatial grid is 1, the grid cell is determined to belong to the macroscopic flow category of Newtonian fluid, and the spatial grid is marked and partitioned into the large blood vessel computational domain; when the output of the region partitioning decision function is 0, the grid cell is determined to belong to the category of smoke-like capillary groups with significant microscopic non-Newtonian characteristics, and the spatial grid is marked and partitioned into the small blood vessel computational domain covering the smoke-like blood vessels.

[0027] This step enables precise binary classification decoupling of the three-dimensional whole-skull complex network based on real physical scale.

[0028] S3. Configure a Newtonian fluid constitutive equation with constant viscosity parameters in the computational domain of large blood vessels, and configure a non-Newtonian fluid constitutive equation with nonlinear functions of local flow field variables as viscosity parameters in the computational domain of small blood vessels. Further, step S3, configuring the Newtonian fluid constitutive equation with a constant viscosity parameter into the large blood vessel computational domain, includes: The fluid control equations are defined within the computational domain of the large blood vessel, and the strain rate tensor of the flow field is extracted. The stress tensor in the computational domain of large blood vessels is set to have a linear mapping relationship with the strain rate tensor; By assigning a fixed dynamic viscosity value to the linear mapping relationship, a Newtonian fluid constitutive equation with a constant viscosity parameter is generated.

[0029] Further, step S3, which involves configuring the non-Newtonian fluid constitutive equation with viscosity parameters as nonlinear functions of local flow field variables in the computational domain of the small blood vessels, includes: Extract the flow field deformation parameters within the computational domain of small blood vessels, and calculate the second invariant of the flow field strain rate tensor. Pre-defined blood rheology consistency index and shear cutoff index are extracted, and a piecewise power-law nonlinear viscosity function with the second invariant as the independent variable is constructed. By embedding a piecewise power-law nonlinear viscosity function into the stress tensor calculation, a non-Newtonian fluid constitutive equation is constructed in which the viscosity parameter is a nonlinear function of the local flow field variables.

[0030] Specifically, the Navier-Stokes fluid control equations for incompressible blood are defined within the computational domain of large blood vessels. The strain rate tensor, characterizing the degree of deformation of fluid micro-elements at the spatial grid nodes within the computational domain of large blood vessels, is extracted. Let the blood flow velocity field vector in three-dimensional space be... The strain rate tensor of the flow field is calculated by obtaining the symmetric part of the velocity field vector gradient. The Cauchy stress tensor of the flow field within the large vessel computational domain is set as a linear mapping relationship directly proportional to the strain rate tensor, and a fixed dynamic viscosity value is assigned to the proportionality coefficient in this linear mapping relationship. In this scheme, a fixed dynamic viscosity value is used. The preferred value is 0.0035 Pa·s. This generates a Newtonian fluid constitutive equation with a constant viscosity parameter, calculated as follows: ; in The Cauchy stress tensor of the flow field within the computational domain of the large blood vessel. For fluid pressure, For unit tensors, The spatial gradient matrix of the velocity field vector. Let be the transpose of the velocity gradient matrix. In this computational domain, the blood viscosity remains constant and does not change with the local velocity gradient, characterizing the macroscopic Newtonian fluid rheological properties of high-velocity blood in a large artery.

[0031] Within the computational domain of small blood vessels, flow field deformation parameters characterizing the microscopic aggregation and shearing state of erythrocytes are extracted. Based on the extracted velocity field vector, the second invariant of the flow field strain rate tensor is calculated. To address the abnormally proliferating microvascular network in patients with Moyamoya disease, characteristic parameters of blood materials were introduced. Extract the preset blood rheological consistency index With shear cutoff index In this scheme, the characteristic parameters of blood materials... The preferred value is 0.78, a blood rheology consistency index. Preferably 0.015 Pa·s n , cutoff index Preferred 100s -1 Construct a second invariant using the flow field strain rate tensor. Piecewise power-law nonlinear viscosity function with independent variables The formula for calculating this function is as follows: ; Where parameters Characterizing the degree of shear-thinning flow properties of blood materials. When the second invariant of the local flow field... Strictly greater than the set shear cutoff index At that time, the fluid viscosity decreases exponentially with the deformation parameter; when the second invariant... Less than or equal to the shear cutoff index At that time, the fluid viscosity is kept constant to avoid singular divergence in numerical solutions at extremely low shear rates.

[0032] The piecewise power-law nonlinear viscosity function obtained from the calculation As a dynamically evolving rheological coefficient, it is embedded in the stress tensor calculation rule of the local fluid, thereby constructing a non-Newtonian fluid constitutive equation whose viscosity parameter is a nonlinear function of the local flow field variables, and its expression is as follows: ; in The updated Cauchy stress tensor of the flow field within the computational domain for small blood vessels is calculated. Because the viscosity parameter is dynamically updated with real-time feedback of the velocity gradient, this equation can accurately capture the rheological differences of blood in different velocity regions within the smoke-like collateral vascular network.

[0033] This step involves matching stress calculation rules for different regions as needed within the same whole-skull model, enabling precise separation and characterization of macroscopic high-speed blood flow and microscopic complex capillary rheological behavior.

[0034] S4. Based on the blood flow physiological mechanism of the matching three-dimensional geometric model, set the outlet impedance boundary conditions, and integrate the Newtonian fluid constitutive equation, the non-Newtonian fluid constitutive equation and the vascular wall solid dynamics equation to construct a fluid-structure coupling mathematical and physical model for the target cranial vascular network. Further, step S4, which involves setting the outlet impedance boundary conditions based on the blood flow physiological mechanism of the matching three-dimensional geometric model, includes: Identify the geometric outlet of the corresponding end microcirculation bed in the 3D geometric model and configure the lumped parameter model; Based on the physiological mechanisms of blood flow, the dynamic resistance values ​​and compliance parameters at each geometric outlet were calculated; The dynamic resistance value and compliance parameter are mapped to time-varying pressure constraint terms to generate the outlet impedance boundary condition.

[0035] Furthermore, step S4, which involves fusing the Newtonian fluid constitutive equation, the non-Newtonian fluid constitutive equation, and the blood vessel wall solid dynamics equation, includes: A Saint-Venant-Kirchhoff hyperelastic material model incorporating the Lamé coefficient was introduced to construct a solid dynamic equation for the blood vessel wall that characterizes the deformation features of the blood vessel wall under stress. A continuity constraint that satisfies motion coordination and interfacial stress balance is established at the interface between the internal fluid and the solid blood vessel wall, and a dynamic mesh equation is introduced to track the boundary deformation of the computational domain. By using continuity constraints, the state variables of the Newtonian fluid constitutive equation and the non-Newtonian fluid constitutive equation, as well as the displacement variables of the blood vessel wall solid dynamics equation, are integrated into a fully implicit system matrix, thus completing the construction of a fluid-structure interaction mathematical physics model.

[0036] Specifically, the truncated arterial terminal sections in the 3D geometric model are identified as the geometric exits of the corresponding cerebral microcirculation bed. Lumped parameter models representing the downstream unmodeled microvascular network are configured at each geometric exit. Based on the blood flow physiology of collateral circulation in patients with Moyamoya disease, dynamic microvascular resistance values ​​allocated to each geometric exit are calculated. With vascular compliance parameters Extract the cross-sectional volumetric blood flow through the geometric exit at any given time. Using the input source variable, the time-varying pressure constraint term at the outlet is solved by combining the impedance equivalence principle. Its first-order differential equation is as follows: ; The time-varying pressure constraint terms obtained from the solution The outlet normal stress boundary condition of the fluid domain is set as follows, applied to the outlet boundary of the fluid control equation. ,in The unit external normal vector at the fluid outlet is used to generate the outlet impedance boundary condition to simulate the physiological regulatory resistance and elastic buffering effect of the terminal microvascular bed.

[0037] For the blood vessel wall region with elastic deformation capability, a Saint-Venant-Kirchhoff hyperelastic material model is introduced to characterize the stress deformation features of the blood vessel wall. Let the spatial displacement vector of the solid blood vessel wall be... Its deformation gradient tensor is defined as Based on Lamé's first coefficient Lamé's second coefficient, which characterizes shear resistance. Construct a solid stress tensor to describe the stress state of blood vessel wall deformation. In this scheme, the Lamé first coefficient of the blood vessel wall The preferred pressure is 2 MPa, and the Lamé second coefficient is... The preferred pressure is 0.33 MPa. The formula for calculating its nonlinear constitutive equation is as follows: ; Constructing the solid stress tensor Substituting into the elastic dynamics conservation system, we can construct the complete solid dynamic equations of the blood vessel wall: ; in The set baseline density of blood vessel wall tissue, This represents the vector distribution of volume forces acting within a solid structure. In this scheme, the basic density of the blood vessel wall tissue is... The preferred value is 1050 kg / m 3 .

[0038] A continuity constraint is established at the boundary between the internal fluid domain and the solid domain of the blood vessel wall. To ensure the realism of the physical interaction boundary, the normal stress applied by the fluid and the normal stress fed back by the solid must be equal in magnitude and opposite in direction to achieve interfacial stress balance. Simultaneously, the velocities of the fluid boundary particles and the solid wall must be completely consistent to achieve motion coordination. The interfacial continuity equations are expressed as follows: ; ; in Let be the unit outward normal vector of the solid wall. Let be the fluid velocity vector at the interface. To track the changes in the boundary space of the fluid computational domain caused by pipe wall stress deformation, a dynamic mesh displacement variable representing the spatial node displacement is introduced. Set dynamic mesh displacement matching constraints at the fluid-structure interface Simultaneously solve the Laplace smoothing equations of the moving mesh within the fluid domain. ,in The virtual mesh stress tensor is used to control the mesh deformation stiffness.

[0039] Through the aforementioned joint constraint mechanism of interface stress balance and displacement velocity, the state variables of the Newtonian fluid constitutive equation describing the flow field properties and the state variables of the non-Newtonian fluid constitutive equation, as well as the displacement variables and dynamic mesh displacement variables of the vascular wall solid dynamic equation describing the wall deformation, are strongly coupled and correlated at the boundary condition level. Then, they are uniformly mapped and integrated into the fully implicit system algebraic matrix, thus completing the overall construction of the fluid-structure interaction mathematical and physical model for the target cranial vascular network.

[0040] This step enables the integration of non-Newtonian blood flow with the physical phenomenon of hyperelastic blood vessel wall deformation into a tightly coupled overall system architecture.

[0041] S5. Using the outlet impedance boundary condition as an input constraint, a computational fluid dynamics parallel solution algorithm is used to perform coupled iterative solution of the fluid-structure interaction mathematical physics model, and output the three-dimensional dynamic hemodynamic parameters of the target cranial vascular network.

[0042] Further, step S5, which involves using a computational fluid dynamics parallel solution algorithm to perform coupled iterative solutions on the fluid-structure interaction mathematical-physical model, includes: The fluid-structure interaction mathematical physics model is mapped to a parallel computing architecture, and the computational domain is spatially partitioned using a domain decomposition strategy. A nonlinear solver is introduced to perform strongly coupled iterative calculations on the partitioned fluid-structure interaction mathematical-physical model. Monitor the iterative residuals during the solution process, and extract and output three-dimensional dynamic hemodynamic parameters when the convergence condition is met.

[0043] Specifically, the integrated fluid-structure interaction mathematical physics model and the set outlet impedance boundary conditions are discretized spatially and temporally to construct a global, fully implicit, nonlinear algebraic equation system describing the physical state of the entire system within a single time step. This global nonlinear equation system is then expressed as a system function. , where the state variable vector It includes fluid velocity, fluid pressure, and solid displacement variables of the blood vessel wall at all nodes in a three-dimensional mesh system. This system of nonlinear algebraic equations is mapped to a high-performance parallel multi-core processor architecture, and the computational domain is spatially partitioned using a domain decomposition strategy. The overall massive target cranial vascular computational domain is cut and allocated into multiple local sub-computational domains with some spatial geometric overlap, allowing each parallel processor core to undertake the numerical evolution task within its respective local sub-computational domain.

[0044] A Newton-Krylov-Schwarz nonlinear solver with preprocessing mechanisms is introduced to perform globally strongly coupled iterative calculations on the partitioned fluid-structure interaction mathematical physics model. In the outer nonlinear iterative loop, the Newton-Krylov iterative algorithm is used for approximate linearization expansion. For the ... The next nonlinear iteration requires solving the following large sparse linear correction equation: ; in The global Jacobian matrix generated by taking the partial derivatives of a system of nonlinear algebraic equations with respect to the state variable vector. This represents the increment vector of state variables obtained in the current iteration step. To improve the solution efficiency and numerical stability of this system of equations under a parallel architecture, a generalized minimum residual algorithm is used for iterative approximation in the inner layer, and an additive Schwarz preprocessor is introduced simultaneously. The condition number of the linear correction equation is optimized. The construction formula of this preprocessing operator is as follows: ; in The total number of local sub-computational domains formed by spatial partitioning. To reduce the global state variable space mapping to the nth Boolean constraint operator matrices for local sub-computational domains For the corresponding extension operator matrix, For the first The inverse matrix of the local Jacobian submatrix extracted over each sub-computational domain. This preprocessing operator matrix is ​​left-multiplied and applied to both sides of the linear correction equation to accelerate the calculation of variable increments, and based on... The state variable vector of the entire system is updated synchronously.

[0045] During the nonlinear solution process at each time step, the iterative residuals of the algebraic equation system are continuously monitored. A dimensionless convergence criterion based on the relative decrease ratio of the initial residuals is established, and its criterion formula is as follows: ; in The rules for vector 2 norm operations, This is the initial state residual vector calculated before the first Newton iteration at the current time step. A preset minimum convergence tolerance threshold is used to control computational accuracy. In this scheme, the convergence tolerance threshold is... Preferred When the system detects the current number When the iterative residual of each step strictly satisfies the above criteria, the fully coupled system is considered to have achieved physical equilibrium convergence at the current time. At this point, iteration is stopped, and the state variable vector at this time is extracted. The physical quantities hidden within it output three-dimensional dynamic hemodynamic parameters, including the three-dimensional spatial velocity distribution field, pressure drop across vascular stenosis lesions, and shear stress of the vessel wall, which are dynamically evolved over time within the whole brain grid system.

[0046] This step, through a preprocessing acceleration mechanism under a distributed architecture, enables efficient and stable numerical computation of a large-scale blood-flow-structure coupling system in the whole brain containing tens of millions of nodes and exhibiting extremely strong nonlinearity.

[0047] Example 2 Existing traditional medical imaging techniques struggle to quantitatively capture complex transient hemodynamic changes. Conventional numerical simulation methods typically employ simplified models, directly excluding tiny moyamoya blood vessels and applying only a single Newtonian fluid assumption, neglecting the non-Newtonian nonlinear rheological properties of blood within microvessels. Furthermore, they lack strong consideration of the elastic deformation of the vessel wall and the physiological impedance at the distal end, leading to distorted simulation data and failing to meet the requirements for accurate prediction. To address these issues, this invention provides a hemodynamic simulation system for moyamoya disease that incorporates non-Newtonian fluid dynamics, the structure of which is as follows: Figure 4 As shown. The specific implementation process of this system is as follows: The model reconstruction module is used to acquire brain medical imaging data of patients with moyamoya disease, extract vascular branch structures including intracranial and extracranial major arteries, bypass vessels and moyamoya vessels, construct the target brain vascular network, and reconstruct the three-dimensional geometric model of the target brain vascular network. The computational domain partitioning module is used to extract the diameter feature parameters of each blood vessel branch structure in the three-dimensional geometric model, and divide the three-dimensional geometric model into a large blood vessel computational domain and a small blood vessel computational domain covering the smoke blood vessels based on the diameter feature parameters. The constitutive equation configuration module is used to configure a Newtonian fluid constitutive equation with a constant viscosity parameter for the large blood vessel computational domain, and at the same time configure a non-Newtonian fluid constitutive equation with a viscosity parameter that is a nonlinear function of the local flow field variable for the small blood vessel computational domain. The fluid-structure interaction modeling module is used to set the outlet impedance boundary conditions based on the blood flow physiological mechanism matching the three-dimensional geometric model, and to integrate the Newtonian fluid constitutive equation, the non-Newtonian fluid constitutive equation and the vascular wall solid dynamics equation to construct a fluid-structure interaction mathematical and physical model for the target cranial vascular network. The simulation solution output module is used to take the outlet impedance boundary condition as input constraint, and use a computational fluid dynamics parallel solution algorithm to perform coupled iterative solution on the fluid-structure interaction mathematical and physical model, and output the three-dimensional dynamic hemodynamic parameters of the target cranial vascular network.

[0048] Specifically, the system first transforms the input magnetic resonance angiography image into a three-dimensional geometric model containing an internal finite element tetrahedral mesh using image thresholding and topology reconstruction algorithms. Next, it calculates the equivalent diameter of the cross-section of the vascular branch centerline. Using a preset critical physical scale threshold as a boundary, the whole-brain spatial mesh is precisely decoupled into two computational domains: macroscopic large blood vessels and microscopic fine, foggy blood vessels. Subsequently, the underlying physics engine calls a constant dynamic viscosity algorithm and a piecewise power-law nonlinear algorithm based on the second invariant of the flow field strain rate to dynamically update the fluid stress tensor under different computational domain scales. Furthermore, by combining the lumped-parameter impedance model for simulating terminal microcirculation, the Saint-Venant-Kirchhoff hyperelastic material model for characterizing vessel wall elasticity, and dynamic mesh interface constraint technology, the complex physical quantities of fluid-solid interaction are integrated into a unified global nonlinear algebraic equation system. Finally, through a spatial domain decomposition strategy, this system is mapped to a multi-core parallel computing architecture. A nonlinear solver with an additive Schwarz preprocessing mechanism is used to perform strongly coupled iterative evolution, thereby efficiently and faithfully calculating and outputting three-dimensional dynamic hemodynamic data that reflects the patient's real blood rheological differences and vessel wall deformation.

[0049] Finally, it should be noted that the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for simulating hemodynamics in Moyamoya disease using non-Newtonian fluids, characterized in that, include: S1. Obtain brain medical imaging data of patients with moyamoya disease, extract vascular branch structures including intracranial and extracranial major arteries, bypass vessels and moyamoya vessels, construct the target brain vascular network, and reconstruct the three-dimensional geometric model of the target brain vascular network. S2. Extract the diameter feature parameters of each blood vessel branch structure in the three-dimensional geometric model, and divide the three-dimensional geometric model into a large blood vessel calculation domain and a small blood vessel calculation domain covering the smoke blood vessels based on the diameter feature parameters. S3. Configure a Newtonian fluid constitutive equation with a constant viscosity parameter in the large blood vessel computational domain, and configure a non-Newtonian fluid constitutive equation with a nonlinear function of local flow field variables as the viscosity parameter in the small blood vessel computational domain. S4. Set the outlet impedance boundary conditions based on the blood flow physiological mechanism matching the three-dimensional geometric model, and integrate the Newtonian fluid constitutive equation, the non-Newtonian fluid constitutive equation and the blood vessel wall solid dynamics equation to construct a fluid-structure coupling mathematical and physical model for the target cranial vascular network. S5. Using the outlet impedance boundary condition as an input constraint, the fluid-structure interaction mathematical-physical model is coupled and iteratively solved using a computational fluid dynamics parallel solution algorithm, and the three-dimensional dynamic hemodynamic parameters of the target cranial vascular network are output.

2. The method for simulating hemodynamics of Moyamoya disease using non-Newtonian fluids according to claim 1, characterized in that, Step S1, the step of reconstructing the three-dimensional geometric model of the target cranial vascular network, includes: Magnetic resonance angiography image data is extracted from the acquired cranial medical imaging data, and image segmentation processing is performed on it to obtain the pixel matrix of the target region; Spatial topology reconstruction is performed on the pixel matrix to generate a three-dimensional point cloud contour containing the vascular branch structure of the intracranial and extracranial aorta, bypass vessels, and smoke vessels. Perform surface smoothing and mesh generation on the three-dimensional point cloud contour to output the three-dimensional geometric model.

3. The method for simulating hemodynamics of Moyamoya disease using non-Newtonian fluids according to claim 1, characterized in that, In step S2, the step of dividing the three-dimensional geometric model into a large blood vessel computational domain and a small blood vessel computational domain covering the smoke-filled blood vessels based on the pipe diameter characteristic parameters includes: Extract the centerline of each vascular branch structure in the three-dimensional geometric model, and calculate the equivalent diameter of its cross section as the tube diameter characteristic parameter; The tube diameter feature parameters are compared with a preset blood vessel scale threshold, and spatial grids that are greater than or equal to the scale threshold and spatial grids that are less than the scale threshold are marked respectively. Spatial grids larger than or equal to the scale threshold are divided into the large blood vessel computation domain, and spatial grids smaller than the scale threshold are divided into the small blood vessel computation domain.

4. The method for simulating hemodynamics of Moyamoya disease using non-Newtonian fluids according to claim 1, characterized in that, Step S3, the step of configuring a Newtonian fluid constitutive equation with a constant viscosity parameter in the large blood vessel computational domain, includes: The fluid control equations are defined within the computational domain of the large blood vessel, and the strain rate tensor of the flow field is extracted. The stress tensor in the computational domain of the large blood vessel is set to have a linear mapping relationship with the strain rate tensor; A fixed dynamic viscosity value is assigned to the linear mapping relationship to generate a Newtonian fluid constitutive equation with a constant viscosity parameter.

5. The method for simulating hemodynamics of Moyamoya disease using non-Newtonian fluids according to claim 1, characterized in that, Step S3, the step of configuring the non-Newtonian fluid constitutive equation with viscosity parameters as nonlinear functions of local flow field variables in the computational domain of the small blood vessels, includes: Extract the flow field deformation parameters within the computational domain of the tiny blood vessels, and calculate the second invariant of the flow field strain rate tensor; Extract the preset blood rheology consistency index and shear cutoff index, and construct a piecewise power-law nonlinear viscosity function with the second invariant as the independent variable; By embedding the piecewise power-law nonlinear viscosity function into the stress tensor calculation, a non-Newtonian fluid constitutive equation is constructed in which the viscosity parameter is a nonlinear function of the local flow field variable.

6. The method for simulating hemodynamics of Moyamoya disease using non-Newtonian fluids according to claim 1, characterized in that, In step S4, the step of setting the outlet impedance boundary conditions based on the blood flow physiological mechanism matching the three-dimensional geometric model includes: Identify the geometric outlet of the corresponding end microcirculation bed in the three-dimensional geometric model and configure the lumped parameter model; Based on the aforementioned blood flow physiological mechanism, the dynamic resistance value and compliance parameter at each of the geometric outlets are calculated; The dynamic resistance value and compliance parameter are mapped to a time-varying pressure constraint term to generate the outlet impedance boundary condition.

7. The method for simulating hemodynamics of Moyamoya disease using non-Newtonian fluids according to claim 1, characterized in that, Step S4, the step of fusing the Newtonian fluid constitutive equation, the non-Newtonian fluid constitutive equation, and the blood vessel wall solid dynamics equation, includes: A Saint-Venant-Kirchhoff hyperelastic material model incorporating the Lamé coefficient is introduced to construct the solid dynamic equations of the blood vessel wall that characterize the deformation features of the blood vessel wall under stress. At the interface between the internal fluid and the solid blood vessel wall, a continuity constraint that satisfies both motion coordination and interfacial stress balance is established, and a dynamic mesh equation is introduced to track the boundary deformation of the computational domain. By using the continuity constraint, the state variables of the Newtonian fluid constitutive equation and the non-Newtonian fluid constitutive equation, as well as the displacement variables of the blood vessel wall solid dynamics equation, are integrated into a fully implicit system matrix, thus completing the construction of the fluid-structure interaction mathematical physics model.

8. The method for simulating hemodynamics of Moyamoya disease using non-Newtonian fluids according to claim 1, characterized in that, Step S5, the step of using a computational fluid dynamics parallel solution algorithm to perform coupled iterative solution of the fluid-structure interaction mathematical-physical model, includes: The fluid-structure interaction mathematical physics model is mapped to a parallel computing architecture, and the computational domain is spatially partitioned using a domain decomposition strategy. A nonlinear solver is introduced to perform strongly coupled iterative calculations on the partitioned fluid-structure interaction mathematical-physical model. The iterative residuals during the solution process are monitored, and the three-dimensional dynamic hemodynamic parameters are extracted and output when the convergence condition is met.

9. A simulation system for hemodynamics of moyamoya disease incorporating non-Newtonian fluids, characterized in that, The system for the hemodynamic simulation method of Moyamoya disease incorporating non-Newtonian fluids as described in any one of claims 1-8 comprises: The model reconstruction module is used to acquire brain medical imaging data of patients with moyamoya disease, extract vascular branch structures including intracranial and extracranial major arteries, bypass vessels and moyamoya vessels, construct the target brain vascular network, and reconstruct the three-dimensional geometric model of the target brain vascular network. The computational domain partitioning module is used to extract the diameter feature parameters of each blood vessel branch structure in the three-dimensional geometric model, and divide the three-dimensional geometric model into a large blood vessel computational domain and a small blood vessel computational domain covering the smoke blood vessels based on the diameter feature parameters. The constitutive equation configuration module is used to configure a Newtonian fluid constitutive equation with a constant viscosity parameter for the large blood vessel computational domain, and at the same time configure a non-Newtonian fluid constitutive equation with a viscosity parameter that is a nonlinear function of the local flow field variable for the small blood vessel computational domain. The fluid-structure interaction modeling module is used to set the outlet impedance boundary conditions based on the blood flow physiological mechanism matching the three-dimensional geometric model, and to integrate the Newtonian fluid constitutive equation, the non-Newtonian fluid constitutive equation and the vascular wall solid dynamics equation to construct a fluid-structure interaction mathematical and physical model for the target cranial vascular network. The simulation solution output module is used to take the outlet impedance boundary condition as input constraint, and use a computational fluid dynamics parallel solution algorithm to perform coupled iterative solution on the fluid-structure interaction mathematical and physical model, and output the three-dimensional dynamic hemodynamic parameters of the target cranial vascular network.

Images