基于非平衡壁面模型的心脏瓣膜血液流动LES/IB计算方法

本发明属于计算流体力学领域，具体涉及一种基于非平衡壁面模型的心脏瓣膜血液流动les/ib计算方法。

背景技术：

1、心血管疾病是当今几种主要的致死性疾病之一。由于疾病或先天问题，自然的心脏瓣膜会变得有缺陷，通过外科手术用人工瓣膜代替人体瓣膜已有五十多年的历史。其中，双叶机械心脏瓣膜(bmhv：bileaflet mechanicalheart valve)是最典型的机械瓣膜，如图1所示。然而，由于会产生非生理性的流动结构，mhv容易导致血栓栓塞等并发症。减弱非生理性流动，降低血栓栓塞风险，是mhv流动研究的重点。bhv的研究注重于瓣膜的有限元动力学分析以及血液与瓣膜之间的相互耦合作用。计算流体力学在心脏瓣膜流动研究中的显著优势是能够高分辨率地模拟三维流场，是推进相关血液动力学特性理解的有效选择。

2、mhv的存在使得流场结构变得更加复杂，一个脉动周期内，会出现流动向紊流转捩以及后来的再层流化等复杂现象。ib(immersed boundary)方法使得运动界面流动的数值模拟可以在固定网格上进行，避免了网格的重新生成以及不同网格之间的数值传递，是目前心脏血液流动模拟中的主流方法。les(large-eddy simulation)模拟中，网格计算节点正比于re13/7，壁面模化技术能够降低心脏瓣膜血液流动les模拟时近壁面附近的网格尺度要求，使得基于ib的湍流模拟数值方法可行。壁面模化技术一般都是针对贴体网格方法提出，ib湍流模拟中的壁面模化已引起了广泛的关注。然而，目前针对心脏瓣膜血液流动的数值模拟，均采用基于平衡流假设的、简化边界层方程的平衡性壁面模型。这类壁面模型当中一些与当地边界层状态相关的前提假设(如速度型形状和平衡流动)使得其在心脏瓣膜血液流动的某些区域(如转捩和分离点附近)变得不再准确。

技术实现思路

1、本发明针对现有技术中的不足，提供一种基于非平衡壁面模型的心脏瓣膜血液流动les/ib计算方法，实现心脏瓣膜血液流动的复杂流场的准确预测。

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

3、一种基于非平衡壁面模型的心脏瓣膜血液流动les/ib计算方法，包括：

4、步骤1、根据给定的血管壁面的形状，生成流场域的背景贴体曲面euler网格，建立全局直角坐标系；

5、步骤2、将机械瓣膜表面进行三角离散，生成瓣膜lagrange网格，将瓣膜lagrange网格各点的初始位置设置为瓣膜瓣叶关闭状态的位置，初始速度设置为零，并浸没在背景贴体曲面euler网格中；

6、步骤3、设置背景贴体曲面euler网格的出口和入口边界条件以及其余边界的边界条件，设置初始血液流场解；

7、步骤4、应用流固强耦合方法，以初始血液流场解、瓣膜lagrange网格各点的初始位置和初始速度为基础，考虑边界条件，迭代求解每一个时间步上的血液流场解、瓣膜位置解和瓣膜速度解。

8、为优化上述技术方案，采取的具体措施还包括：

9、进一步地，步骤3具体为：

10、背景贴体曲面euler网格的右边界为出口，采用对流边界条件，uout是出口的流向平均速度，i取值为1、2、3，u1、u2、u3分别为全局直角坐标下的流向速度、法向速度及展向速度；t和x分别表示时间和全局直角坐标系的流向；背景贴体曲面euler网格左边界为入口，入口血液流量为给定的随时间变化的时间函数，初始血液流场解w0＝(p0，u0)为静止流场的解，p0表示初始时刻血液的压强，u0表示初始时刻血液的流动速度矢量；背景贴体曲面euler网格的其余边界设为无滑移边界条件，背景贴体曲面euler网格节点上u1＝0，u2＝0，u3＝0。

11、进一步地，步骤4具体为：基于第n时间步上血液流场解wn＝(pn，un)、瓣膜lagrange网格各点的位置解xn和速度解vn，n＝0～nmax-1，应用强耦合方法求n+1时间步上血液流场解wn+1、瓣膜lagrange网格各点的位置解xn+1和瓣膜速度解vn+1，nmax为最大时间步，pn表示第n时间步上血液的压强，un表示第n时间步上血液的速度矢量，包括以下步骤：

12、步骤4.1、设置时间步n内的子迭代步序号k＝1，令wk＝wn，xk＝xn，vk＝vn，n＝0～nmax-1；wk为第k子迭代步上血液流场解，xk和vk分别为第k子迭代步瓣膜lagrange网格上各点的位置解和速度解；wn＝(pn，un)，wn表示第n时间步的血液流场解，pn表示第n时间步上血液的压强，un表示n时间步上血液的速度矢量，xn表示第n时间步的瓣膜lagrange网格上各点的位置解，vn表示第n时间步的瓣膜lagrange网格上各点的速度解；nmax表示最大时间步；

13、步骤4.2、根据第k子迭代步瓣膜lagrange网格上各点的位置解xk，将背景贴体曲面euler网格点分为解域内部计算点、外部相关点和外部不参与计算的无关点；

14、步骤4.3、根据第k子迭代步瓣膜lagrange网格上各点的位置解xk及对应的速度解vk，利用边界条件及法向积分的湍流边界层方程，求解非平衡壁面模型，获得当地摩擦速度；

15、步骤4.4、利用当地摩擦速度计算瓣膜lagrange网格上各点的壁面剪切应力；利用壁面剪切应力计算外部相关点上的当地流向速度根据第k子迭代步瓣膜lagrange网格上各点的速度解vk、无穿透边界条件和沿瓣膜壁面法向的线性插值求解当地法向动量方程，获得第k子迭代步外部相关点di上的当地法向速度和压强通过当地坐标系与全局直角坐标系之间的转换关系，将当地流向速度和当地法向速度转换到全局直角坐标系上，得到第k子迭代步外部相关点di在全局直角坐标下的速度矢量获得外部相关点di上的第k子迭代步血液流场解

16、步骤4.5、利用三角线性插值确定当前第k子迭代步参考点fi上的湍流粘性系数然后将第k子迭代步外部相关点di上的湍流粘性系数设为与第k子迭代步参考点fi上的湍流粘性系数相等；

17、步骤4.6、基于外部相关点di上的第k子迭代步血液流场解和湍流粘性系数通过虚拟时间推进获得第k+1迭代步上内部计算点上流动方程的解wk+1；

18、步骤4.7、确定瓣膜lagrange网格各点所在的euler网格单元，利用三角线性插值获得瓣膜lagrange网格点上的流体力；

19、步骤4.8、根据瓣膜lagrange网格点上的流体力求解瓣膜的结构运动方程，得到瓣膜打开角度，所述瓣膜的结构运动方程如下：

20、

21、式中，ired为瓣膜约化的转动惯量，θk+1表示第k+1子迭代步上瓣膜打开角度，t为时间，mi为瓣膜表面所受转动矩的积分；

22、采用亚松弛方法优化当前第k+1子迭代步的瓣膜打开角度θk+1，得到第k+1子迭代步上松弛后的瓣膜打开角度

23、通过松弛后的瓣膜打开角度确定第k+1迭代步上瓣膜lagrange网格各点的位置解xk+1，并计算得到第k+1子迭代步上瓣膜lagrange网格各点的运动速度vk+1＝(xk+1-xk)/δτ，δτ为虚拟时间步长；

24、步骤4.9、检验瓣膜位置解的收敛性，判断是否满足||xk+1-xk||<eps或||θk+1-θk||<eps1；eps和eps1分别为瓣膜位置解和瓣膜转角的收敛指标；||·||表示二范数，若瓣膜位置解未收敛，则返回步骤4.2，若瓣膜位置解收敛，令wn+1＝wk+1，xn+1＝xk+1，vn+1＝vk+1，获得n+1时间步上的血液流场解wn+1、瓣膜lagrange网格各点的位置解xn+1和速度解vn+1。

25、进一步地，步骤4.2具体为：

26、根据第k子迭代步瓣膜lagrange网格上各点的位置解xk，将背景贴体曲面euler网格点分为解域内部计算点、外部相关点和外部不参与计算的无关点；处于血液计算域的所有网格节点均为内部计算点，处于瓣膜结构内部，并且通过背景贴体曲面euler网格线与内部计算点相连的网格点为外部相关点；剩下的背景贴体曲面euler网格点为外部无关点。

27、进一步地，步骤4.3具体为：

28、根据第k子迭代步瓣膜lagrange网格上各点的位置解xk及对应的速度解vk，利用边界条件及法向积分的湍流边界层方程，求解非平衡壁面模型，获得当地摩擦速度uτ；所述非平衡壁面模型如下：

29、

30、式中，u表示当地流向速度，uτ为当地摩擦速度，y表示到瓣膜壁面的法向距离，εi为粘性子层高度，δv＝μ/ρuτ为与粘性子层相关的长度尺度，μ为层流粘性系数，ρ表示血液流体的密度，未知参数c保证非平衡壁面模型中的线性层与对数层匹配，k为von karman常数，δy为湍流边界层方程做法向积分的高度，未知参数a表示湍流边界层方程中非平衡项的作用。

31、进一步地，步骤4.4具体为：

32、利用当地摩擦速度uτ计算瓣膜lagrange网格上各点的壁面剪切应力τw，公式表达如下：

33、

34、ρ表示血液流体的密度；

35、利用壁面剪切应力τw计算第k子迭代步外部相关点di上的当地流向速度公式如下：

36、

37、式中，表示第k子迭代步参考点fi的当地流向速度，|fidi|表示参考点fi和外部相关点di之间的距离，μ为层流粘性系数，表示第k子迭代步参考点fi上的湍流粘性系数，ρ表示血液的密度，的值利用三角线性插值确定；

38、根据第k子迭代步瓣膜lagrange网格上各点的速度解vk、无穿透边界条件和沿瓣膜壁面法向的线性插值求解当地法向动量方程，获得第k子迭代步外部相关点di上的当地法向速度和压强所述当地法向动量方程如下：

39、

40、式中，表示压强的法向梯度，γ表示瓣膜壁面的当地法向，表示瓣膜壁面单位法向矢量；

41、通过当地坐标系与全局直角坐标系之间的转换关系，将当地流向速度和当地法向速度转换到全局直角坐标系上，得到第k子迭代步外部相关点di在全局直角坐标下的速度矢量获得外部相关点di上的第k子迭代步血液流场解

42、本发明的有益效果是：

43、提供了一种准确的心脏瓣膜流动的流场预测和流动结构分析的计算方法。该方法具有以下优点：(1)将曲线贴体网格和浸没边界相结合，可以提高血管壁面附近流场的精度，避免了运动瓣膜边界的网格实时更新，与传统的贴体网格方法相比，在模拟兼有运动和静止固壁的流动问题时，具有显著优势；(2)针对心脏瓣膜血管流动问题的ib模拟，有别于之前一直采用的基于平衡性壁面模型，本发明将一种非平衡壁面模型与dfd锐化界面ib方法相结合，提高了复杂瓣膜流场计算的准确性。