CN107133397B

CN107133397B - 一种基于ale法对生物瓣膜进行双向流固耦合分析的方法

Info

Publication number: CN107133397B
Application number: CN201710287293.XA
Authority: CN
Inventors: 袁泉; 朱宏伟; 申炳申; 王志超; 唐丹
Original assignee: Shandong University
Current assignee: Shandong University
Priority date: 2017-04-27
Filing date: 2017-04-27
Publication date: 2018-10-19
Anticipated expiration: 2037-04-27
Also published as: CN107133397A

Abstract

本发明公开了一种基于ALE法对生物瓣膜进行双向流固耦合分析的方法。本发明对血液与瓣叶的耦合分析采用双向耦合分析，在处理问题时，每一步均计算血液对瓣叶的力的作用，同时也计算瓣叶变形也对血液流动造成的影响，采用基于ALE法对生物瓣膜进行双向流固耦合分析，使得力学性能分析同时兼顾流体(血液)和固体(瓣膜及血管壁)，同样使得分析更为精确可靠。方法包括如下步骤：第一步：建立血液及生物瓣膜几何模型；第二步：基于流体域与固体域相互耦合，构建流体控制方程、固体控制方程；第三步：构建任意变量的流固耦合方程，对瓣叶及血管壁的变形、生物瓣膜表面应力变化和瓣叶开口面积变化进行分析。

Description

一种基于ALE法对生物瓣膜进行双向流固耦合分析的方法

技术领域

本发明属于流固耦合数值模拟控制领域，具体涉及一种基于ALE法对生物瓣膜进行双向流固耦合分析的方法。

背景技术

心脏是人体的重要器官，心脏瓣膜一旦发生病变就会危及人的生命。目前还没有方法能对功能失调的心脏瓣膜进行药物治疗，修复或更换瓣膜仍是唯一能降低心脏瓣膜病发病率和死亡率的治疗方法。生物瓣膜具有良好的力学性能逐渐成为瓣膜置换手术的首选，但仍存在着耐久性方面的问题，主要表现在瓣膜材料的疲劳和破坏，以及多发于儿童和青年患者的钙化现象。利用计算流体力学的方法对生物瓣膜与血液的耦合过程进行分析，获得瓣叶的变形特点及应力分布情况，对研究瓣膜损坏及钙化的原因有重要的理论意义和参考价值。

由于生物瓣膜的临床应用中，必须通过缝合环与血管壁连在一起，因此在分析中必须考虑血管壁的形态及受力变化，然而，现有研究中，往往只是单独分析瓣膜应力变化，而忽视了血管壁对瓣叶与血液流固耦合的影响，因此造成仿真结果与实际结果出现偏差；同时，现有技术中对生物瓣膜研究中普遍采用单向耦合分析，这使得流体与固体间的参数传递过于简化，这一分析方法在较小变形时可以采用，但对于生物瓣膜研究中的复杂大变形，结果则会存在一定失真。因此，建立新的分析方法，使得力学性能分析能够同时兼顾流体(血液)和固体(瓣膜及血管壁)，对生物瓣膜的设计制作无疑具有重要的指导意义。

发明内容

针对上述现有技术，本发明提供一种基于ALE法对生物瓣膜进行双向流固耦合分析的方法。本发明对血液与瓣叶的耦合分析采用双向耦合分析，在处理问题时，每一步均计算血液对瓣叶的力的作用，同时也计算瓣叶变形也对血液流动造成的影响，再将计算的结果带入到下一次的计算中，依次迭代，计算整个流固耦合过程，得到瓣叶与血液相互作用的结果，从而使结果更加准确。

具体的，本发明的技术方案如下：

一种基于ALE法对生物瓣膜进行双向流固耦合分析的方法，包括如下步骤：

第一步：建立血液及生物瓣膜几何模型；

第二步：基于流体域(血液)与固体域(生物瓣膜)相互耦合，构建流体(血液)控制方程、固体(生物瓣膜)控制方程；

第三步：基于ALE法对生物瓣膜进行双向流固耦合过程进行分析，构建任意变量的流固耦合方程，使用LS-DYNA对瓣叶及血管壁的变形、生物瓣膜表面应力变化和瓣叶开口面积变化进行分析。

其中，所述第一步中，对生物瓣膜构建几何模型时，采用不可压缩弹性材料对瓣叶进行描述，具体步骤包括：

(1)给出瓣叶的位移方程，X表示空间中单元的初始位置，x表示单元受力移动后的位置向量，y表示位移向量；

x(X,t)＝X+y(X,t)

则单元的变形梯度可以表示为：

(2)由步骤(1)得到的单元的变形梯度，得到单元的柯西应力张量C、应变张量E以及变形张量J，由于瓣叶模型被认为是不可压缩的，所以变形张量J等于1；具体柯西应力张量C、应变张量E以及变形张量J的描述方程如下：

C＝F^TF

J＝detF(X,t)≡1；

所述第一步中，对血液构建几何模型时，用牛顿流体来描述，将血液描述为不可压缩的粘性牛顿流体，具体方程如下：

▽·u＝0

所述第二步中，流体控制方程由流体的连续性方程和流体的动量方程组成，分别为：

所述第二步中，由于生物瓣膜的瓣叶属于超弹性的材料，在分析时对生物瓣膜的本构方程进行简化，采用弹性体进行描述，则固体控制方程即为固体的连续方程，即为：

式中，X表示拉格朗日坐标，ρ_s表示瓣叶结构的密度，f_i为体积力，u为瓣叶的位移；

进一步的，由于几何模型是连续的，为了进行有限元分析，需要对血液模型进行离散处理，采用有限元法对流体控制方程进行离散，具体的，将流体的连续性方程和流体的动量方程转化成时间离散的形式，描述方程如下：

进一步的，由于上述微分方程不能直接进行求解，需要应用有限元中伽辽金法对上式进行转化，伽辽金法可以将微分方程转化为线性方程组，进而方便求解，上述两式转化后的方程写成矩阵的形式如下：

MΔρ_n+1＝Δtf_ρ；

M(ΔU_i)_n+1＝Δtf_u；

其中Δρ，ΔU_i为具体的质量分量和动量分量，M为给定的质量矩阵，表达式为：

所述第三步中，采用罚函数法处理流体(血液)与固体(瓣叶)的耦合，构建任意变量f的流固耦合方程如下：

其中，χ为ALE坐标，w表示瓣叶速度与网格速度之差；

上述方程具体构建方式如下：

对于罚函数方式，在结构积分耦合点找到相应的流体材料点，并跟踪他们的相对位移d(又称穿透距离)然后根据相对位移的大小分别对结构和流体施加节点力，这个耦合过程是同时计算结构体和流体的受力，所以是双向的；

具体的，流体单元和结构单元的相对位移d可表示为：

d＝u_s-u_f；

式中，u_s和u_f分别表示相互接触的结构单元和流体单元的位移；

对于耦合界面的平衡条件，表示为：

F_S＝F_f＝F；

F表示在耦合界面瓣叶与血液的相互作用力，可以通过罚函数计算：

F＝kd+cd；

其中，k表示罚函数的刚度，c表示阻尼系数；

利用罚函数法，即得上述流固耦合方程。

进一步的，第三步中对生物瓣膜表面应力变化分析包括对瓣叶等效应力变化分析、最大主应力变化分析及最大剪切应力变化分析。

本发明的有益效果：本发明中对生物瓣膜与血液的耦合分析模型进行改进，在结构模型中加入了血管壁模型，考虑了血管壁对瓣叶与血液流固耦合的影响，使仿真结果更为准确；同时，本发明采用基于ALE法对生物瓣膜进行双向流固耦合分析，使得力学性能分析同时兼顾流体(血液)和固体(瓣膜及血管壁)，使得分析更为精确可靠，这对于生物瓣膜的设计制作无疑具有重要的指导意义。

附图说明

图1为罚函数耦合法则示意图；

图2为血液几何模型；

图3为椭球面瓣叶几何模型；

图4为分离瓣叶模型和血管壁模型；

图5为包含血管壁的瓣叶分析模型；

图6为血液分析模型，其中(a)为血管壁模型，(b)血液模型；

图7为SOLID164单元示意图；

图8为扫略方式得到的血管壁网格；

图9为扫略方式得到的瓣叶网格；

图10为流域和血液的网格；

图11为网格划分后各个PART的网格数量示意图；

图12为流量表示的主动脉瓣血流速度图；

图13为简化的血液入口速度；

图14为瓣叶的变形过程图；

图15为血管壁的变形过程图；

图16为瓣叶等效应力的变化图；

图17为瓣叶最大主应力的变化图；

图18为瓣叶最大剪切应力的变化图；

图19为开口面积比变化曲线。

具体实施方式

应该指出，以下详细说明都是例示性的，旨在对本申请提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本申请所属技术领域的普通技术人员通常理解的相同含义。

需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本申请的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”或“包括”时，其指明存在特征、步骤、操作、器件、组件或它们的组合。

正如背景技术所介绍的，现有研究中，往往只是单独分析瓣膜应力变化，而忽视了血管壁对瓣叶与血液流固耦合的影响，因此造成仿真结果与实际结果出现偏差；同时，现有技术中对生物瓣膜研究中普遍采用单向耦合分析，这使得流体与固体间的参数传递过于简化，导致结果失真。

有鉴于此，本申请的一种典型的实施方式中，提供一种基于ALE法对生物瓣膜进行双向流固耦合分析的方法，包括如下步骤：

第一步：建立血液及生物瓣膜几何模型；

其中，生物瓣膜应与原生心脏瓣膜有类似结构，在血液作用下发生大变形的同时也控制流经血液的方向。原生心脏瓣膜是生物组织，由弹性纤维和胶原纤维组成。这样的组成结构使得瓣膜具有非均质性、粘弹性和各向异性。瓣叶的不同区域在厚度和成分上存在差异，这也是导致瓣叶各向异性的原因。生物瓣膜取自哺乳动物的心包瓣(通常是指牛或者猪)，经过化学试剂戊二醛处理鞣制，使之具备一定的形态和力学特性。

故在第一步中，对生物瓣膜构建几何模型时，采用不可压缩弹性材料对瓣叶进行描述，具体步骤包括：

x(X,t)＝X+y(X,t)

则单元的变形梯度可以表示为：

(2)由变形梯度，得到单元的柯西应力张量C、应变张量E以及变形张量J，由于瓣叶模型被认为是不可压缩的，所以变形张量等于1；单元的柯西应力张量、应变张量以及变形张量的具体描述方程如下：

C＝F^TF

J＝detF(X,t)≡1；

而流体具有一定程度的可压缩性，当作用在流体上的压力增加时，其体积或密度将会减小。对于血液来说，体积弹性模量和水相似，可以认为是不可压缩流体。血液具有一定的粘性，通常用卡森关系来描述血液流变学。血液呈现非线性的切应力与切应变率关系，在低切应变率时表现的更为明显，而在相对较高的切应变率下，粘度趋于常数，称为牛顿流体。对于大血管中的血液流动，如主动脉的血液流动，切应变率相对较高。

故在第一步中，对血液构建几何模型时，用牛顿流体来描述；将血液描述为不可压缩的粘性牛顿流体，具体描述方程如下：

▽·u＝0。

拉格朗日法和欧拉法常被用来模拟各种工程问题。在拉格朗日方法中，节点和材料点同时移动，因此该方法适合于固体力学模拟。在欧拉法中，网格在空间中固定而材料点在变形期间移动，因此它适合于应用流体力学模拟。Hughes T J R等人开发了一种称为任意拉格朗日-欧拉法的更通用的方法，简称ALE法。ALE法其包含拉格朗日法和欧拉法的优点，同时避免它们的缺点。在ALE分析中，有限元网格既不依附在到材料点中，也不固定在空间中，并且网格和材料点可以单独移动。ALE法中的网格运动可以减少网格失真并将网格传送到需要更高网格密度的区域，而不增加求解参数的数量。

ALE法引入了拉格朗日和欧拉坐标之外的第三个任意参照坐标，根据流体区域的边界构造新网格，并进行每步的计算和迭代，通过网格的不断划分避免计算中产生网格的畸变，本质上是一种坐标描述的方法。以此为基础，在血液的边界条件下先求解血液的相关性能参数；然后将流体计算的参量，如剪切应力等传递到固体网格中，求解固体的位移和压力值；接着进行循环不断迭代更新，流体参量进一步求解，直至计算收敛结束。

故在第二步中，将流体控制方程写成ALE法描述的形式，具体的，流体控制方程由流体的连续性方程和流体的动量方程组成，分别为：

进一步的，由于几何模型是连续的，为了进行有限元分析，需要对血液模型进行离散处理，采用有限元法对流体控制方程进行离散，具体的，将流体的连续性方程和流体的动量方程转化成时间离散的形式，具体方程如下：

MΔρ_n+1＝Δtf_ρ；

M(ΔU_i)_n+1＝Δtf_u；

由于任意的变量f在参照坐标下的全导数可由下式导出，其中X为拉格朗日坐标，x为欧拉坐标，w为相对速度并以此出控制方程的ALE表达。

其中，X，x分别表示拉格朗日坐标和欧拉坐标，w表示拉格朗日坐标下运动速度υ与欧拉坐标下速度u之间的相对速度。

而在ALE求解中，需要定义耦合边界的网格运动，一般使用表示网格运动的速度。在本申请中网格运动的速度用参考点的速度ω_i进行定义：

本申请第三步中采用罚函数法处理血液与瓣叶的耦合。对于罚函数方式，首先在结构积分耦合点找到相应的流体材料点，并跟踪他们的相对位移d(又称穿透距离)，然后根据相对位移的大小分别对结构和流体施加节点力，这个耦合过程是同时计算结构体和流体的受力，所以是双向的。罚函数法耦合过程如附图1所示。

流体单元和结构单元的相对位移d可表示为：

d＝u_s-u_f

式中，u_s和u_f分别表示，相互接触的结构单元和流体单元的位移。

对于耦合界面的平衡条件，可以表示为：

F_S＝F_f＝F

F＝kd+cd

上式中，k表示罚函数的刚度，c表示阻尼系数

利用罚函数法，可以得到求解任意变量f的流固耦合方程：

其中，χ为ALE坐标，w表示瓣叶速度与网格速度之差。

为了使得本领域技术人员能够更加清楚地了解本申请的技术方案，以下将结合具体的实施例详细说明本申请的技术方案。

实施例1

基于ALE的双向耦合基本分析流程

使用ANASYS/LS-DYNA软件进行仿真分析一般可以分为以下步骤：

(1)建立几何模型：模型可以在ultraedit中直接建立，也可以通过三维建模软件建立然后倒入到分析软件中。

(2)划分网格：几何模型是连续的，为了进行有限元分析，需要对模型进行离散处理。根据模型的特点，选择合适的划分方法，网格质量对计算结果有着直接影响。

(3)设定边界条件及材料属性：结合分析模型的实际情况施加载荷或者给定初始速度。分析求解前需要定义单元类型、材料模型，需要定义密度、杨氏模量、泊松比等参数。

(4)输出K文件并进行修改：一些参数的定义无法在图形界面中完成，需要借助Ultraedit或文本编辑器对K文件进行修改

(5)LS-DYNA Solver求解：选择LS-DYNA Solver求解器进行求解前应设置仿真的持续时间，设置与输出有关的信息，这些都要在关键字中进行定义。下面对各步骤作详细分析。

模型的建立

血液几何模型

通过对MRI(磁共振)数据进行处理，得到瓣叶所在动脉窦的三维模型。因扫描所得的模型与血液几何模型的形状相同，所以利用得到的尺寸参数对血管窦部分的血液进行三维重建。图2中旋转对称结构模型即为血液几何模型。

生物瓣膜几何模型

生物瓣膜模型的建立采用几何建模的方法，利用椭球面建立与血管壁相匹配的圆锥面，进而得到椭球面生物瓣膜的曲面模型。在三维建模软件中对曲面进行阵列处理，设置瓣叶的厚度为0.4mm，得到图3所示椭球面瓣叶几何模型。

双向流固耦合分析模型的建立

1 瓣叶分析模型的建立

导入的模型并不能直接进行双向流固耦合分析，需要对模型进行一定的处理，便于后续步骤的进行。

首先将瓣叶模型从装配模型中分离出来，如图4所示，左侧为血管壁模型，右侧为瓣叶模型。再将分离的瓣叶模型略放大，使瓣叶的边缘与血管有相交，以便下一步划分网格。最后将放大后的瓣叶与血管壁模型重新装配，即可得到图5中的包含血管壁的瓣叶分析模型。从模型中可以看到，瓣叶的缝合边缘已经与血管壁的外侧重合。

2 血液分析模型的建立

LS-DYNA中的模型建立流体域需将血管壁模型进行拓展，在血液入口及出口区域增加长度，并将壳结构转化为实体结构，得到的血液流域模型如图6所示，(a)为血管壁模型，(b)为血液模型。

双向耦合分析前处理

1 单元类型及算法选择

(1)单元类型

LS-DYNA提供了多种单元材料可供选择，如：LINK160，BEAM161，SHELL163，SOLID164等。在对瓣膜进行分析时，我们采用SOLID164三维实体单元描述瓣叶，因为使用三维实体单元可以较为准确的分析模型的应力与应变。

在LS-DYNA中使用关键字*ELEMENT_SOLID对三维实体单元进行定义，描述了SOLID164单元的几何图形、节点位置和坐标系。该单元有8个节点，如图7所示。SOLID164输入概括:节点——(I,J,K,L,M,N,O,P)，自由度——(UX,UY,UZ,VX,VY,VZ,AX,AY,AZ)。对于显式动力分析，V(X,Y,Z)提供节点的速度以及A(X,Y,Z)提供节点的加速度。虽然V(X,Y,Z)和A(X,Y,Z)出现在固定约束的地方，但他们不是真实的物理约束。

(2)算法选择

在LS-DYNA中三维单元有三种基本的算法：Lagrangian、Eulerisn和ALE(任意拉格朗日欧拉算法)算法，是由关键字*SECTION_SOLID中的ELFORM控制的。选择ELFORM的数值为1，表示中心单点积分的ALE多物质单元(一个单元内可以包含多种物质)。

Lagrange算法的单元网格附着在材料上，随着材料的流动而产生单元网格的变形。在结构变形过于巨大时，有可能使有限元网格造成严重畸变，引起数值计算的困难，甚至程序终止运算。ALE算法可以克服单元严重畸变引起的数值计算困难，并实现流体-固体耦合的动态分析。ALE算法先执行一个或几个Lagrange步计算，此时单元网格随材料的移动而产生变形，然后执行ALE步计算：(1)保持变形后的物体边界条件，对内部单元进行重分网格，网格的拓扑关系保持不变，称为Smooth Step；(2)将变形网格中的单元变量(密度、能量、应力张量等)和节点速度矢量输运到重分后的新网格中，称为Advection Step。用户可以定义ALE步的开始和终止时间，以及其频率。

LS-DYNA一般有两种耦合方式，一种为约束加速度或速度的方式，一种为罚函数方式。本发明中选择ALE法描述流体部分，固体部分选择用Lagrange法进行描述，通过罚耦合算法，使用关键字为*CONSTRAINED_LAGRANGE_IN_SOLID对流固耦合进行定义，其中MASTER定义为1表示ALE算法标书流体，PFAC定义罚因子，选择默认的0.1，DIREC表示耦合的方向选择，选择2表示在法线方向耦合—仅在压缩方向，这种耦合方向求解更为稳定。

多物质单元的定义：LS-DYNA用*ALE_MULTI-MATERIAL_GROUP关键字把各种材料绑定在一个单元。对于拉格朗日算法而言，知道单元里面只能是一种物质，即单元是附着在物质上的，而对于欧拉算法或ALE算法而言，一个单元中可以包含不同的材料，也只有这样才能在空间网格中完成物质的输送，如流体分析中的多相流等。

2 网格划分

网格划分的总体思路：可将包含瓣叶的血管壁和流域分别划分网格，最后移动到一起。网格划分的目的是得到六面体网格，这样可以得到较好的计算结果。

通常，采用扫略方式形成网格是一种非常好的方式，对于复杂几何实体，经过一些简单的切分处理，就可以自动形成规整的六面体网格，它比映射网格划分方式具有更大的优势和灵活性。对于由面经过拖拉、旋转、偏移等方式生成的复杂三维实体而言，可先在原始面上生成壳单元形式的面网格，然后在生成体的同时自动形成三维实体网格。对于已经形成好了的三维复杂实体，如果其在某个方向上的拓扑形式始终保持一致，则可用扫略方式划分功能来划分网格。这两种方式形成的单元几乎都是六面体单元。对于血管壁和瓣叶均采用此种方式，得到的网格如图8和图9所示。

要注意的是，在初始时刻各PART的网格划分问题，在各PART的界面上，一定要确保网格一致，此时网格上节点是公共的。

对流域和血液进行划分需要遵循同样的原则，在划分的时候需要注意的是流域和血液必须共节点。划分后的网格如下图10所示，图中深色区域为血液的网格，浅色部分为流域的网格。

在统计中查看网格状态，六个part分别包含的单元数目如图11所示，其中瓣叶为PART2，PART3，PART4，血管为part1，流场为part5，血液为part6。

3 材料属性定义

(1)瓣叶及血管材料属性定义

生物瓣膜的瓣叶由生物组织构成，一般为牛或猪的心包经过化学处理得到。这两种材料包含按一定规律排列的胶原纤维构成，在各个方向上的力学性能是不同的。由于测量的差异，在非线性分析中所采用的材料参数也存在一定差异，没有明确的标准。由于生物瓣膜的三片瓣叶是对称分布的，本文中我们将其简化为各向同性线弹性材料，定义瓣叶的泊松比为0.4，弹性模量7×10⁶pa。血管材料的泊松比为0.3，弹性模量为5×10⁷Pa。

(2)血液材料属性

在对流体材料处理的过程中，就需要同时使用两种方式来描述材料，用本构模型和状态方程(EOS)来同时描述一种材料的特性。在LS-DYNA中提供一种空材料模式*MAT_NULL用来描述具有流体行为的材料(如空气，水等)。

结合本发明中血液流体力学性质，流体采用Null空材料模型和Gruneisen状态方程来定义血液流体。在本文中将血液模型简化为不可压缩的粘性流体，血液密度为1.105g/cm³，粘度为0.00466Pa·s，弹性模量为3×10⁸N/m²。

4 边界条件的施加及沙漏控制

在LS-DYNA中对于边界条件的设置较为简单，不需要设置入口出口以及流固耦合面，只需要将初始条件施加在对应的PART即可。在血液与瓣叶的流固耦合问题分析中，需要设置血液的初始速度。经过核磁共振测速得到真实人体主动脉瓣血液流速如图12所示。该曲线是在横截面积约400mm²的条件下得到的，用流量除以横截面积可将上述血液入口速度曲线的流量形式进行转化。不考虑返流的情况下得到简化后的血液速度曲线，如图13所示。

出口处选择用压力来定义血液的边界条件，本文中对主动脉压强进行简化，设置出口的压强为0Pa，血液运动主要取决于血液入口速度。

除了入口速度，其他的参数定义如下：血液在本文中被看作不解压缩的粘性流体，血液密度为1.105g/cm³，血液的粘度为0.00466Pa·s，弹性模量为3×10⁸N/m²，温度按照人体正常温度定义为309.15K。

设置流体速度的关键字为*DEFINE_CURVE。

单元的形状经常为四面体或六面体，出于计算速度与稳定性考虑,常采用单点积分，这样虽然使计算量减少，但是也导致了沙漏(零能模式)。血液冲击瓣叶的过程时间较短，同时血液具有一定的粘性，所以我们在*HOURGLASS关键字中将IHQ定义为1，表示粘性的沙漏控制。沙漏系数采用默认值0.1。

生物瓣膜双向流固耦合过程分析

1 瓣叶及血管壁的变形分析

(1)瓣叶的变形分析

在8核心主频为2.1HZ，16G内存的工作站上运行LS-DYNA Solver进行求解，运算16小时后求解完成。使用后处理软件LS-PrePost打开分析结果，得到瓣叶变形过程如图14所示。

瓣叶最大变形为图中的深色区域，随着时间增大，瓣叶的变形量也逐渐增加。三片瓣叶的变形是对称的，在25ms之前，瓣叶的变形较快，中心开口区域近似三角形。在45ms时，瓣叶完全打开，开口区域近似圆形，此时瓣叶的变形量最大，为11.88mm。在瓣叶开启的过程中，生物瓣膜瓣叶的最大变形区域由瓣叶腹部区域逐渐转移至瓣叶自由边的边缘。在20ms时瓣叶自由边的中间区域出现了一定的卷曲现象，随着瓣叶变形的增大，卷曲的程度逐渐减小。

(2)血管壁变形如图15所示：

血管壁的材料相对于瓣叶材料有较高的弹性模量，当血液与瓣叶进行耦合时，由于瓣叶与血管壁相连，血管壁也会受到力的作用产生一定的变形。血管壁的变形集中在瓣叶开启的0-25ms，在瓣叶完全开启时，血管壁的变形较小。

2 生物瓣膜表面应力变化分析

(1)瓣叶等效应力变化

从图16中可以看出，瓣叶的等效应力随着瓣叶变形而增大。在45ms时，等效应力出现最大值3.349Mpa。最大的等效应力主要分布在瓣叶的缝合边以及瓣叶缝合边与自由边的交界处，而瓣叶自由边的边缘等效应力较小。可以观察到，双向流固耦合分析中瓣叶等效应力的分布呈阶梯状，从瓣叶的缝合边到自由边，等效应力逐渐减小，这一规律特征是在以往单向耦合分析中没有体现的。

(2)最大主应力变化

最大主应力由于分布在瓣叶的结合边边缘，而且由于方向的问题，从俯视图上看不到最大主应力的分布位置，所以选择主视图查看瓣叶的主应力变化，如图17。可以看出，最大主应力随着瓣叶的变形逐渐增大。最大主应力出现在瓣叶缝合边的边缘。整体上看最大主应力的分布较为均匀，在瓣叶的腹部以及自由边的边缘，最大主应力较小。

(3)最大剪切应力变化

由图18可以观察到：瓣叶表面的最大剪切应力主要集中在瓣叶的缝合边及缝合边与自由边的交界处，最大剪切应力随着瓣叶的变形增加而增大，在45ms时达到最大值。

3 开口面积变化分析

本发明使用开口面积比来描述瓣叶的开口面积变化过程。开口面积比是指瓣叶自由边面积即开口面积与瓣叶总面积的比值。不同时刻瓣叶的有效开口面积比如表格1所示。由于双向耦合中，血管壁的变形较大，瓣叶的面积也发生一定的变化，所以为了得到准确的开口面积比，在表格中加入了瓣叶面积。

表1 不同时刻瓣叶开口面积与开口面积比表

时间/ms	0	5	10	15	20
						开口面积(cm²)	0	3.82	11.62	25.89	41.30
瓣叶面积(cm²)	86.00	86.98	88.55	87.83	87.05
						开口面积比	0	0.0439	0.1312	0.2948	0.4744
时间/ms	25	30	35	40	45
						开口面积(cm²)	51.74	55.48	55.56	55.58	55.61
瓣叶面积(cm²)	86.32	85.91	84.56	84.10	84.00
						开口面积比	0.5994	0.6458	0.6570	0.6609	0.6620

根据表格中瓣膜开口面积比的数据，绘制双向流固耦合分析开口面积比变化的折线图如图19所示。

折线图可以清晰的看到开口面积比的变化：0-25ms的过程中，开有面积比急剧增加。最终在45ms时，开口面积比达到最大，此时的开口面积即为有效开口面积。双向耦合分析中，开口面积的变化并不是线性的，在初始开启阶段的增加较快，然后缓慢增加至最大值，这与单向耦合分析中，开口面积比近似线性的均匀增加明显不同。

以上所述仅为本申请的优选实施例而已，并不用于限制本申请，对于本领域的技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本申请的保护范围之内。

Claims

1.一种基于ALE法对生物瓣膜进行双向流固耦合分析的方法，其特征在于，包括如下步骤：

第一步：建立血液及生物瓣膜几何模型；

第二步：基于流体域与固体域相互耦合，构建流体控制方程、固体控制方程；

第三步：基于ALE法对生物瓣膜进行双向流固耦合过程进行分析，构建任意变量的流固耦合方程，使用LS-DYNA对瓣叶及血管壁的变形、生物瓣膜表面应力变化和瓣叶开口面积变化进行分析；

所述第一步中，对生物瓣膜构建几何模型时，采用不可压缩弹性材料对瓣叶进行描述，具体步骤包括：

(1)给出瓣叶的位移方程，X表示空间中单元的初始位置，x表示单元受力移动后的位置向量，y表示位移向量，t表示单元变形时间；

x(X,t)＝X+y(X,t)

则单元的变形梯度可以表示为：

其中I为变形常量，值为1，

(2)由步骤(1)得到的单元的变形梯度，得到单元的柯西应力张量C、应变张量E以及变形张量J，由于瓣叶模型被认为是不可压缩的，所以变形张量等于1；单元的柯西应力张量C、应变张量E以及变形张量J的具体描述方程如下：

C＝F^TF

J＝detF(X,t)≡1；

其中，det为行列式运算符；

所述第三步中，采用罚函数法处理流体与固体的耦合，构建任意变量f的流固耦合方程如下：

其中，χ为ALE坐标，w表示瓣叶速度与网格速度之差；

所述流固耦合方程具体构建方式如下：

对于罚函数方式，在结构积分耦合点找到相应的流体材料点，并跟踪他们的相对位移d然后根据相对位移的大小分别对结构和流体施加节点力，这个耦合过程是同时计算结构体和流体的受力，所以是双向的；

具体的，流体单元和结构单元的相对位移d可表示为：

d＝u_s-u_f；

对于耦合界面的平衡条件，表示为：

F_S＝F_f＝F；

F＝kd+cd；

其中，k表示罚函数的刚度，c表示阻尼系数；

利用罚函数法，即得所述流固耦合方程。

2.如权利要求1所述的分析方法，其特征在于，所述第一步中，对血液构建几何模型时，用牛顿流体来描述；将血液描述为不可压缩的粘性牛顿流体，具体描述方程如下：

其中，u表示占据某空间点的物质点的位移，t是时间，ρ是血液密度，μ是血液粘度，p是流体压力。

3.如权利要求1所述的分析方法中，其特征在于，所述第二步中，流体控制方程由流体的连续性方程和流体的动量方程组成，分别为：

式中，x为欧拉坐标，u为欧拉坐标下速度，w表示拉格朗日坐标下物质点的相对速度，t是时间，ρ是血液密度，p是流体压力，τ为切应力，g为重力加速度。

4.如权利要求1所述的分析方法中，其特征在于，所述第二步中，由于生物瓣膜的瓣叶属于超弹性的材料，在分析时对生物瓣膜的本构方程进行简化，采用弹性体进行描述，则固体控制方程即为固体的连续方程，即为：

式中，X表示拉格朗日坐标，ρ_s表示瓣叶结构的密度，f_i为体积力，u为瓣叶的位移，σ为正应力，t是时间。

5.如权利要求4所述的分析方法，其特征在于，将流体的连续性方程和流体的动量方程转化成时间离散的形式，转化后的微分方程如下：

式中，n为迭代次数，x表示欧拉坐标，u为欧拉坐标下速度，w表示拉格朗日坐标下物质点的相对速度，t是时间，ρ是血液密度，p是流体压力，τ为切应力，g为重力加速度。

6.如权利要求5所述的分析方法，其特征在于，使用伽辽金法将权利要求5所述微分方程转化为线性方程组，进而方便求解，上述两式转化后的方程写成矩阵的形式如下：

MΔρ_n+1＝Δtf_ρ；

M(ΔU_i)_n+1＝Δtf_u

其中Δρ，ΔU_i为具体的质量分量和动量分量，Δt为时间增量，f为体积力，M为给定的质量矩阵，表达式为：

N^TN表示n阶矩阵，表示球坐标系下的参数分量。

7.如权利要求1所述的分析方法，其特征在于，所述第三步中对生物瓣膜表面应力变化分析包括对瓣叶等效应力变化分析、最大主应力变化分析及最大剪切应力变化分析。