一种计算机实现的非牛顿血流仿真方法、系统及介质与流程

本发明涉及计算流体力学、医学仿真、生物力学、图神经网络和计算机辅助诊疗，尤其涉及一种计算机实现的非牛顿血流仿真方法、系统及介质。

背景技术：

1、血流动力学仿真在心血管疾病诊断、手术规划、支架和导管设计、药物输送评估以及医学教育中具有重要应用价值。传统的计算流体动力学（computational fluiddynamics，cfd）方法，如基于navier-stokes方程求解的有限体积法或有限元法，虽然精度较高，但通常需要复杂的前处理（网格划分、边界层设置）和较长的计算时间，不适合术中实时仿真或大规模参数探索。

2、平滑粒子流体动力学（smoothed particle hydrodynamics，sph）是一种无网格拉格朗日方法，天然适合处理自由表面、大变形和复杂边界问题，在血流仿真中已有应用。然而，传统sph方法在处理壁面边界条件时通常依赖罚力模型或虚拟粒子，容易在复杂几何中出现粒子穿透或边界条件松弛现象；此外，血液作为一种非牛顿流体，其黏度随剪切率变化（剪切稀化效应），在微血管、狭窄血管或低剪切率区域，采用牛顿黏度假设会引入显著误差。

3、基于图神经网络（graph neural network，gnn）的学习型求解器已在固体力学和流体力学中展现出潜力，如图网络仿真器、连续卷积网络和消息传递网络可用于加速传统数值仿真。然而，将图神经网络应用于sph血流仿真面临以下挑战：

4、（1）sph粒子相互作用具有天然的成对反对称性，即粒子i对粒子j的作用力与粒子j对粒子i的作用力大小相等、方向相反，而普通图神经网络的消息传递未必满足该条件，训练收敛后可能得到局部能量梯度正确但全局动量不守恒的力场；现有方法通常仅在损失函数中以惩罚项的形式近似施加该约束，在推理阶段不提供任何保证；

5、（2）非牛顿流体的黏度依赖于局部剪切率，需要模型在推理时动态确定有效黏度；现有的学习型sph求解器多将黏度视为固定的全局常数或仅作为节点标量特征附加于网络输入，未使有效黏度直接参与成对相互作用的计算，因而难以在剪切率剧烈变化的狭窄段和近壁面区域正确再现剪切稀化行为；

6、（3）壁面、入口和出口等边界粒子在仿真中的运动受到物理约束，仅仅依靠损失函数进行软约束无法严格保证前向仿真满足边界条件，在长时间步推进中误差会持续累积。

7、因此，需要一种计算机实现的非牛顿血流仿真方法、系统及介质，以实现准确、稳定且实时的非牛顿血流仿真。

技术实现思路

1、本发明提出的一种计算机实现的非牛顿血流仿真方法、系统及介质，以解决上述背景技术中提出的问题。

2、为了实现上述目的，本发明采用了如下技术方案：

3、一种计算机实现的非牛顿血流仿真方法，包括以下步骤：

4、s1：接收表示血管内血液和边界的粒子图，所述粒子图的节点为平滑粒子流体动力学粒子，所述粒子图的边由支撑半径内的粒子邻域关系生成，所述支撑半径为2h，h为粒子光滑长度；

5、s2：基于粒子相对位置、粒子相对速度、平滑粒子流体动力学核函数值和核函数梯度构造边(i，j)的边特征e_ij；

6、s3：根据非牛顿黏度模型，基于粒子i处的局部剪切率γd_i确定该粒子的血液有效黏度μ_eff,i，并将边(i，j)两端粒子的有效黏度μ_eff,i和μ_eff,j作为附加分量并入所述边特征e_ij；

7、s4：利用具有反对称消息传递机制的图神经网络，根据所述粒子图和并入有效黏度后的边特征生成每个流体粒子的加速度；所述反对称消息传递机制对边(i，j)以共享消息函数φ计算第一消息φ(h_i，h_j，e_ij)，并对交换端点后的边特征e_ji计算第二消息φ(h_j，h_i，e_ji)，再以第一消息减去第二消息形成边(i，j)上的定向消息m_ij ＝ φ(h_i，h_j，e_ij) - φ(h_j，h_i，e_ji)，其中e_ji由e_ij按下述规则得到：相对位置、相对速度和核函数梯度分量取相反数，粒子间距离和核函数值保持不变，两端粒子的有效黏度分量互换位置；由此对于边(j，i)恒有m_ji ＝ -m_ij；粒子i的加速度由作用于其上的成对相互作用求和后除以粒子质量得到：a_i ＝ (1/m_i)·σ_(j∈n(i)) m_ij，且σ_i m_i·a_i ＝0；

8、s5：在模型输出阶段对非流体边界粒子施加加速度、速度或位置的硬门控：，其中m_i为粒子i的类型掩码，为该粒子类型所规定的加速度值。

9、优选的，所述s1中，粒子图包括流体、壁面、入口和出口中的至少两类粒子标签，每个粒子具有位置、速度、密度和压力的物理属性，粒子图的节点特征包括位置、速度、密度、压力和粒子类型，粒子图的边特征包括三维相对位置、距离、三维相对速度、核函数值和三维核函数梯度。

10、优选的，所述s2中，平滑粒子流体动力学核函数为wendland-c2核函数、三次样条核函数或高斯核函数；所述wendland-c2核函数为：当0≤q≤2时，；当q＞2时，w(q) ＝ 0；其中q ＝ ||r_ij||/h，α_d为与空间维数相关的归一化常数，三维情形下。

11、优选的，所述s3中，非牛顿黏度模型为carreau-yasuda模型、casson模型、幂律模型或quemada模型；所述carreau-yasuda模型为；局部剪切率γd_i ＝ √(2·d_i:d_i)，。

12、优选的，所述s5中，除对加速度施加所述硬门控外，还对壁面、入口或出口粒子的速度或位置施加投影，使其投影到满足无滑移边界、入口速度边界、出口压力边界或周期边界的子空间。

13、优选的，其还包括通过辛欧拉法、速度verlet法、runge-kutta法或投影积分法更新流体粒子状态，并基于状态方程或连续性方程更新密度和压力；其中状态方程为，γ_t为tait状态方程指数，人工声速满足c_s≥10·u_max。

14、优选的，所述s4中，图神经网络包括编码器和l个反对称消息传递层，所述共享消息函数φ为多层感知器；第l个反对称消息传递层的共享消息函数输出三维的成对相互作用力，聚合之后不再设置解码器。

15、本发明还提出了一种计算机实现的非牛顿血流仿真系统，包括医学图像输入模块、血管分割与粒子化模块、图构建模块、材料与边界条件赋值模块、仿真模块和输出模块。

16、本发明还提出了一种非暂态计算机可读存储介质，其上存储有指令，所述指令在由处理器执行时使计算设备执行上述所述的方法。

17、与现有技术相比，本发明的有益效果是：

18、1、将平滑粒子流体动力学核函数的核函数值和梯度作为图神经网络的边特征，使图神经网络学习到的粒子相互作用天然建立在平滑粒子流体动力学的数学框架之上，网络无需从零学习邻域权重的空间衰减规律；

19、2、采用反对称消息传递机制，使交换端点后的消息满足m_ji＝-m_ij的代数恒等式，从而σ_i m_i·a_i ＝ 0在每一次前向传播中严格成立，全局动量守恒由网络结构保证而非依赖损失函数的软约束；

20、3、将血液有效黏度建模为局部剪切率的函数，并将边两端粒子的有效黏度作为附加分量并入边特征，使有效黏度直接参与成对相互作用的计算，而非仅作为节点标量附加于网络输入，从而使仿真能够在剪切率剧烈变化的狭窄段和近壁面区域正确再现血液的剪切稀化行为；

21、4、在前向传播的输出阶段直接对壁面、入口和出口粒子施加硬门控，使边界粒子的加速度、速度或位置在每一时间步均严格满足规定边界条件，从根本上防止边界条件违反和长时间步推进中的误差累积；

22、5、训练完成后，单个时间步的求解为一次前向图传播，无需组装和求解全局线性系统，在稳态泊肃叶流的验证中相对于解析解的速度剖面误差为0.1%。