大变形管道流动系统振荡空间的确定方法及系统与流程

本发明涉及人体管壁仿真领域，特别是涉及一种大变形管道流动系统振荡空间的确定方法及系统。

背景技术：

大变形可坍塌管道流动出现在人体的许多生理现象中。例如血液流动引起的血管坍塌在自动把血液提供给各内脏的血液循环中发挥了重要的作用；积极地压缩下肢静脉血管使其坍塌可以作为一种有效的治疗方法来阻止深度静脉血栓的形成；血压计利用血管坍塌而产生的科罗特科夫氏音来测量人的血压；在强迫呼吸中，呼气肌肉的收缩会增大胸膜的压力，当压力超过一定值的时候，也会使近端的气管坍塌。管道坍塌后往往会出现所谓的流动极限现象，伴随着流动极限现象的发生往往出现自激振荡，生理上表现为打呼噜、哮喘等病理现象，因此，这些现象的研究对人类的生命健康有重要意义。

由于大变形管道流动系统是一个强非线性系统，理论求解基本不可能实现，现阶段主要采用数值方法进行求解。在现阶段的求解模型中，对管壁描述都采用线弹性模型，而且材料都是可以压缩的，这并不能很好的描述实际的生物材料。

在采用数值方法求解大变形管道流动系统的现有技术中，有人采用二维模型模拟人体管壁，利用特征值方法来研究大变形管道流动系统的振动特性虽然得到了一些不错的结果，但是二维模型跟实际的三维实体还有很大的差距；而采用三维模型模拟人体管壁时，则通过计算系统瞬态解来研究系统的振动特性，然而逐一求解系统瞬态解会耗费大量的时间，更难绘制出系统的振荡空间。由此可知，二维模型与实际不符，而三维模型中采用的求瞬态解的方法由于计算量太大又不能应用在参数空间中的每个点上，现有的方法都不能准确确定大变形管道流动系统的振荡空间。

技术实现要素：

本发明的目的是提供一种大变形管道流动系统振荡空间的确定方法及系统，以解决三维大变形管道流动系统振荡状态的确定效率低的问题。

为实现上述目的，本发明提供了如下方案：

一种大变形管道流动系统振荡空间的确定方法，采用不可压缩的超弹性材料以模拟大变形管道流动系统；所述大变形管道流动系统振荡空间的确定方法包括：

获取大变形管道流动系统参数；所述大变形管道流动系统参数包括管道几何参数、流体参数以及管壁材料参数；所述管道几何参数包括管道总长度、管道内半径、管道直径、管道厚度、管道上游长度、管道中游长度以及管道下游长度；所述流体参数包括密度以及粘性系数；所述管壁材料参数包括剪切模量以及密度；

根据所述大变形管道流动系统参数确定所述大变形管道流动系统的控制方程并获取边界条件；所述边界条件包括入口处边界条件、刚性壁边界条件、弹性段边界条件、弹性端边界条件以及出口处边界条件；

利用有限元法，根据所述控制方程以及所述边界条件确定所述大变形管道流动系统的稳态解；

根据所述稳态解确定特征值；

根据所述特征值确定所述大变形管道流动系统的振荡空间。

可选的，根据所述大变形管道流动系统参数确定所述大变形管道流动系统的控制方程以及边界条件，具体包括：

根据公式以及detf＝1确定所述大变形管道流动系统的控制方程；所述控制方程为无量纲化处理后的控制方程；其中，re是雷诺数，ui，uj为方向流体速度，其中，i＝1时，ui为x轴方向上的流体速度，i＝2时，ui为y轴方向上的流体速度，i＝3时，ui为z轴方向上的流体速度，j＝1时，uj为x轴方向上的流体速度，j＝2时，uj为y轴方向上的流体速度，j＝3时，uj为z轴方向上的流体速度，t为时间，xj为当前构型坐标，p为流体压强，为剪切模量，fia,a为变形梯度张量对初始构型坐标的导数，为固体等效压强，为变形梯度张量的逆对初始构型坐标的导数，为管壁密度，为管壁加速度，detf为变形梯度张量的行列式；

根据公式u＝v＝0以及确定入口处边界条件；其中，u为流体x方向速度，v为流体y方向速度，w为流体z方向速度，x1为x坐标，x2为y坐标；

根据公式u＝v＝w＝0确定刚性壁边界条件；

根据公式以及确定弹性段边界条件；其中，u为流体x方向速度，为管壁x方向速度，v为流体y方向速度，为管壁y方向速度，w为流体z方向速度，为管壁z方向速度；

根据公式以及确定弹性端边界条件；其中，为管壁x方向位移，为管壁y方向位移，为管壁z方向位移；

根据公式σt＝0以及σn＝-pd确定出口处边界条件，其中，σt为出口界面切应力，σn为出口界面正应力，pd为出口外压。

可选的，所述根据所述稳态解确定特征值，具体包括：

根据所述稳态解确定所述大变形管道流动系统的质量矩阵以及非线性矩阵；

根据所述质量矩阵以及所述非线性矩阵确定特征值。

可选的，所述根据所述特征值确定所述大变形管道流动系统的振荡空间，具体包括：

判断所述特征值的实部是否等于0，得到第一判断结果；

若所述第一判断结果表示为所述特征值的实部等于0，确定振荡空间曲线；

根据所述振荡空间曲线确定所述大变形管道流动系统的振荡空间。

一种大变形管道流动系统振荡空间的确定系统，采用不可压缩的超弹性材料以模拟大变形管道流动系统；所述大变形管道流动系统振荡空间的确定系统包括：

参数获取模块，用于获取大变形管道流动系统参数；所述大变形管道流动系统参数包括管道几何参数、流体参数以及管壁材料参数；所述管道几何参数包括管道总长度、管道内半径、管道直径、管道厚度、管道上游长度、管道中游长度以及管道下游长度；所述流体参数包括密度以及粘性系数；所述管壁材料参数包括剪切模量以及密度；

控制方程以及边界条件确定模块，用于根据所述大变形管道流动系统参数确定所述大变形管道流动系统的控制方程以及边界条件；所述边界条件包括入口处边界条件、刚性壁边界条件、弹性段边界条件、弹性端边界条件以及出口处边界条件；

稳态解确定模块，利用有限元法，根据所述控制方程以及所述边界条件确定所述大变形管道流动系统的稳态解；

特征值确定模块，根据所述稳态解确定特征值；

振荡空间确定模块，用于根据所述特征值确定所述大变形管道流动系统的振荡空间。

可选的，所述控制方程以及边界条件确定模块具体包括：

控制方程确定单元，用于根据公式以及detf＝1确定所述大变形管道流动系统的控制方程；所述控制方程为无量纲化处理后的控制方程；其中，re是雷诺数，ui，uj为方向流体速度，其中，i＝1时，ui为x轴方向上的流体速度，i＝2时，ui为y轴方向上的流体速度，i＝3时，ui为z轴方向上的流体速度，j＝1时，uj为x轴方向上的流体速度，j＝2时，uj为y轴方向上的流体速度，j＝3时，uj为z轴方向上的流体速度，t为时间，xj为当前构型坐标，p为流体压强，为剪切模量，fia,a为变形梯度张量对初始构型坐标的导数，为固体等效压强，为变形梯度张量的逆对初始构型坐标的导数，为管壁密度，为管壁加速度，detf为变形梯度张量的行列式；

入口处边界条件确定单元，用于根据公式u＝v＝0以及确定入口处边界条件；其中，u为流体x方向速度，v为流体y方向速度，w为流体z方向速度，x1为x坐标，x2为y坐标；

刚性壁边界条件确定单元，用于根据公式u＝v＝w＝0确定刚性壁边界条件；

弹性段边界条件确定单元，用于根据公式以及确定弹性段边界条件；其中，u为流体x方向速度，为管壁x方向速度，v为流体y方向速度，为管壁y方向速度，w为流体z方向速度，为管壁z方向速度；

弹性端边界条件确定单元，用于根据公式以及确定弹性端边界条件；其中，为管壁x方向位移，为管壁y方向位移，为管壁z方向位移；

出口处边界条件确定单元，用于根据公式σt＝0以及σn＝-pd确定出口处边界条件，其中，σt为出口界面切应力，σn为出口界面正应力，pd为出口外压。

可选的，所述特征值确定模块具体包括：

矩阵确定单元，用于根据所述稳态解确定所述大变形管道流动系统的质量矩阵以及非线性矩阵；

特征值确定单元，用于根据所述质量矩阵以及所述非线性矩阵确定特征值。

可选的，所述振荡空间确定模块具体包括：

第一判断单元，用于判断所述特征值的实部是否等于0，得到第一判断结果；

振荡空间曲线确定单元，用于若所述第一判断结果表示为所述特征值的实部等于0，确定振荡空间曲线；

振荡空间确定单元，用于根据所述振荡空间曲线确定所述大变形管道流动系统的振荡空间。

根据本发明提供的具体实施例，本发明公开了以下技术效果：本发明提供了一种大变形管道流动系统振荡空间的确定方法及系统，通过求解所述大变形管道流动系统的特征值，并根据该特征值确定大变形管道流动系统的动力学状态，避免现有三维模型求解瞬态解计算量太大的问题，准确确定出大变形管道流动系统的振荡空间。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明所提供的大变形管道流动系统振荡空间的确定方法流程图；

图2为本发明所提供的几何模型示意图；

图3为本发明所提供的旋转线示意图；

图4为本发明所提供的网格化的几何模型示意图；

图5为本发明所提供的大变形管道流动系统振荡空间的确定系统结构图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明的目的是提供一种大变形管道流动系统振荡空间的确定方法及系统，能够准确确定出大变形管道流动系统的振荡空间。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1为本发明所提供的大变形管道流动系统振荡空间确定方法流程图，如图1所示，一种大变形管道流动系统振荡空间的确定方法，采用不可压缩的超弹性材料以模拟大变形管道流动系统；所述大变形管道流动系统振荡空间的确定方法包括：

步骤101：获取大变形管道流动系统参数；所述大变形管道流动系统参数包括管道几何参数、流体参数以及管壁材料参数；所述管道参数包括管道总长度、管道内半径、管道直径、管道厚度、管道上游长度、管道中游长度以及管道下游长度；所述流体参数包括密度以及粘性系数；所述管壁材料参数包括剪切模量以及密度。

如图2所示，考虑流体流经一长为l，内半径为r，直径为d圆柱形管道。管道被分为3段，上游和下游为不变形段，中间为变形段，它的厚度为d，并承受大小为pw的外压。上、中、下游三段的长度分别为lu、lm、ld。流体是不可压缩牛顿流体，它的密度是ρ，粘性系数是μ。流动假设为层流。管壁采用不可压的neo-hookean超弹性材料，剪切模量为密度为在流体和管壁相互作用过程中，管壁是流场的边界，同时受到流场对它的作用力。

步骤102：根据所述大变形管道流动系统参数确定所述大变形管道流动系统的控制方程以及边界条件；所述边界条件包括入口处边界条件、刚性壁边界条件、弹性段边界条件、弹性端边界条件以及出口处边界条件。

首先考虑管壁的控制方程。在初始构形下，不考虑体积力情况下的运动方程为：

其中，t表示时间，表示管壁的位移，表示管壁的加速度，πia,a是第一皮奥拉基尔霍夫应力对初始构型坐标的导数，第一皮奥拉基尔霍夫应力可以被表示为：

其中，w是应变能密度函数，是左(右)柯西格林应变张量的第一，二，三不变量的函数，是由于不可压缩约束引入的一个拉格朗日乘子，可以理解为固体等效压力。fia是变形梯度张量，是变形梯度张量的逆，可以表示为

xi表示在当前构形中的位置矢量，xa表示在初始构形中的位置矢量。为了保证材料的不可压缩性，变形梯度张量的行列式需要保持为1，即：

detf＝1(4)

管壁采用的是neo-hookean材料模型，应变能密度函数是：

其中是剪切模量，i是左(右)柯西格林应变张量的第一不变量。它和变形梯度张量fia间的关系是：

利用式(2)、(3)和式(5)，式(1)变为：

其中，是管壁剪切模量；fia是变形梯度张量；fia,a是变形梯度张量对初始构型坐标的导数；是变形梯度张量的逆；是变形梯度张量的逆对初始构型坐标的导数；是管壁密度；是管壁在三个方向的位移，i＝1时，为x轴方向上的管壁位移，i＝2时，为y轴方向上的管壁位移，i＝3时，为z轴方向上的管壁位移；是管壁在三个方向的加速度，i＝1时，为x轴方向上的管壁加速度，i＝2时，为y轴方向上的管壁加速度，i＝3时，为z轴方向上的管壁加速度。

流体运动采用纳维斯托克斯方程和连续性方程来描述，为了方便，引入无量纲参数，带星号的量为无量纲的量。(即：带星号和不带星号表示的意义是一样的，只是无星号是有量纲的，有星号是无量纲的)。

其中：re是雷诺数，u0是不变形段管道流动的平均速度，xi是三个方向坐标，其中，i＝1时，xi为x轴方向上坐标，i＝2时，xi为y轴方向上的坐标，i＝3时，xi为z轴方向上的坐标；ui是三个方向流体速度，其中，i＝1时，ui为x轴方向上的流体速度，i＝2时，ui为y轴方向上的流体速度，i＝3时，ui为z轴方向上的流体速度；p是流体压强；是三个方向管壁位移，其中，i＝1时，为x轴方向上的管壁位移，i＝2时，为y轴方向上的管壁位移，i＝3时，为z轴方向上的管壁位移；t是时间。

按照上面的方面进行无量纲化后，为了简便去掉表示无量纲的星号，可以得到无量纲化后耦合系统的控制方程为：

detf＝1(12)

其中，uj是三个方向流体速度，其中，i＝1时，uj为x轴方向上的流体速度，i＝2时，uj为y轴方向上的流体速度，i＝3时，uj为z轴方向上的流体速度。

各边界条件如下：

入口处(z＝0)：u＝v＝0,z是沿着管道的坐标轴；

刚性壁(0≤z≤lu,r＝r；lu+lm≤z≤l,r＝r)：u＝v＝w＝0；

弹性段(lu≤z≤lu+lm,r≤r≤r+d)：

弹性端(z＝lu,lu+lm,r≤r≤r+d)：

出口处(z＝l)：σt＝0；σn＝-pd；其中，pd是下游的压强，是一个给定的常数。

步骤103：利用有限元法，根据所述控制方程以及所述边界条件确定所述大变形管道流动系统的稳态解。

采用有限元求解过程中，弹性段流体区域的边界会发生改变，为了避免流体网格在边界移动时出现畸形，采用旋转线方法构造流体的自适应网格，如图3所示。在弹性管道内，同心圆区域内网格保持不变，同心圆和管壁之间有很多相连的旋转线。有一点在内同心圆点与管壁点所连成的旋转线k上。内同心圆上的点是固定不动的。管壁运动时，旋转线k保持为直线段绕固定点旋转，线上的节点根据管壁的移动按比例伸缩。变形过后，管壁点运动到点线上节点运动到点点和点之间的关系是：

其中，

不在旋转线上的点坐标由旋转线上的点通过线性差值得到。

如图4所示，在一个管道截面上使用旋转线方法构造网格后，再沿着管道的方向划分，就可以得到整个模型的网格。

由于网格的移动，纳维斯托克斯方程中，速度对时间的偏导数是参照于固定空间的。我们定义参照于网格移动空间的时间导数δ/δt。利用链规则，δ/δt和之间的关系得到移动网格下的纳维斯托克斯方程为：

其中，w为管壁变形时网格的移动速度，ui是流体速度，p是流体压强，re是雷诺数。

有限元分析中，流体采用15节点三棱柱等参单元，固体采用20节点六面体等参单元。为保证解的收敛性，速度和位移分别采用二次形函数ψ^(f)、ψ^(s)进行插值，流体压强和固体压强分别采用线性形函数φ^(p)，进行插值:

式中l1,l2,l3表示三棱柱等参单元中的面积坐标，ξ,η,ζ表示等参单元中的x,y,z三个方向局部坐标。

因此流体单元的15个节点上都有流体速度自由度u,v,w，只有6个顶点有流体压强自由度p。固体单元的20个节点上都有位移自由度只有8个顶点有压强自由度流体单元和固体单元的公共节点上，既有流体自由度也有固体自由度。

采用伽辽金方法对系统的控制方程进行离散，得到有限元方程。权函数和插值函数是相同的函数。在流体控制方程中，对粘性梯度项和压强梯度项进行分部积分得到方程：

其中n是边界面的外法线单位矢量。

首先得到当前构形中的固体控制方程的有限元形式：

其中σij是柯西应力，σij,j是柯西应力对当前构型坐标的导数，mj是当前构形下边界面的外法线单位矢量，δij是delta算子，是流体对管壁的作用力张量，可以表示为：

初始构形与当前构形下的变量有以下关系：

na是初始构形下边界面的外法线单位矢量。结合式(7)、(18)、(19)和(21)得到初始构形下固体控制方程的有限元形式：

由于流体方程中只有时间的一阶导数，而在固体方程中存在对时间的二阶导数，为了使阶数一致，在固体节点上引入固体速度自由度εi(i＝1,2,3)，表示x，y，z三个方向的固体速度。由于引入了新的未知量εi，需要加入新的节点方程：

把所有方程表示为矩阵的形式：

其中u是整体未知量矢量u,v,w是x，y，z三个方向的流体速度，p是流体压强，是x，y，z三个方向的固体管壁位移，固体等效压强，ε1,ε2,ε3是x，y，z三个方向固体管壁速度，m是质量矩阵，k(u)表示非线性刚度矩阵。f是力向量，r是残差向量，当计算收敛时应该为零。

步骤104：根据所述稳态解确定特征值。

步骤105：根据所述特征值确定所述大变形管道流动系统的振荡空间。

所述步骤105具体包括：判断所述特征值的实部是否等于0，若是，确定振荡空间曲线；根据所述振荡空间曲线确定所述大变形管道流动系统的振荡空间。

其中，当特征值的实部等于0时，表示该特征值为所述振荡空间曲线上的一点；调整管壁参数和流体参数，找出振荡空间曲线上的其他点，并连接为振荡空间曲线。

系统特征值是建立在系统稳态解的基础之上的，因此首先需要求解系统的稳态解，在求解稳态解时，去掉式(24)中的时间项，所需求解的稳态方程为：

k(u)u-f＝r＝0

这是一个非线性方程组，采用牛顿拉夫逊法迭代求解。

得到稳态解后求解修正的orr-sommerfeld特征值问题的具体方法大致如下。设对稳态解作微扰动δu,于是为系统的解。令其中ω和分别为相应的复特征值与特征矢量，把带入方程(24)则得到特征方程为:

确定式中矩阵需要用到稳态解的值。

系统的稳定性由特征值ω的值确定。由(25)式能求得的有效特征值的个数与矩阵m秩相同。当有效特征值ω中有任何一个的实部大于零，则系统是非稳定的；当有效特征值ω中所有的实部都小于零，则系统是稳定的；系统稳定与不稳定的临界中性点则对应于有效特征值ω中实部最大者为零的情况。

arpack中求解的特征方程的形式为ax＝λx，它和方程(25)在形式上不一样。因此对方程(25)进行变形为：

得到需要求解的矩阵所求的特征值

图5为本发明所提供的大变形管道流动系统振荡空间的当前状态确定系统结构图，如图5所示，一种大变形管道流动系统振荡空间的当前状态确定系统，采用不可压缩的超弹性材料以模拟大变形管道流动系统；所述大变形管道流动系统振荡空间的当前状态确定系统包括：

参数获取模块501，用于获取大变形管道流动系统参数；所述大变形管道流动系统参数包括管道参数、流体参数以及流体流经管道时的管壁参数；所述管道参数包括管道总长度、管道内半径、管道直径、管道厚度、管道承受外压、管道上游长度、管道中游长度以及管道下游长度；所述流体参数包括密度以及粘性系数；所述管壁参数包括剪切模量以及密度。

控制方程以及边界条件确定模块502，用于根据所述大变形管道流动系统参数确定所述大变形管道流动系统的控制方程以及边界条件；所述边界条件包括入口处边界条件、刚性壁边界条件、弹性段边界条件、弹性端边界条件以及出口处边界条件。

所述控制方程以及边界条件确定模块502具体包括：

控制方程确定单元，用于根据公式以及detf＝1确定所述大变形管道流动系统的控制方程；所述控制方程为无量纲化处理后的控制方程；其中，re是雷诺数，ui，uj为方向流体速度，其中，i＝1时，ui为x轴方向上的流体速度，i＝2时，ui为y轴方向上的流体速度，i＝3时，ui为z轴方向上的流体速度，j＝1时，uj为x轴方向上的流体速度，j＝2时，uj为y轴方向上的流体速度，j＝3时，uj为z轴方向上的流体速度，t为时间，xj为当前构型坐标，p为流体压强，为剪切模量，fia,a为变形梯度张量对初始构型坐标的导数，为固体等效压强，为变形梯度张量的逆对初始构型坐标的导数，为管壁密度，为管壁加速度，detf为变形梯度张量的行列式；

入口处边界条件确定单元，用于根据公式u＝v＝0以及确定入口处边界条件；其中，u为流体x方向速度，v为流体y方向速度，w为流体z方向速度，x1为x坐标，x2为y坐标；

刚性壁边界条件确定单元，用于根据公式u＝v＝w＝0确定刚性壁边界条件；

弹性段边界条件确定单元，用于根据公式以及确定弹性段边界条件；其中，u(t)为流体x方向速度，为管壁x方向速度，v(t)为流体y方向速度，为管壁y方向速度，w(t)为流体z方向速度，为管壁z方向速度；

弹性端边界条件确定单元，用于根据公式以及确定弹性端边界条件；其中，为管壁x方向位移，为管壁y方向位移，为管壁z方向位移；

出口处边界条件确定单元，用于根据公式σt＝0以及σn＝-pd确定出口处边界条件，其中，σt为出口界面切应力，σn为出口界面正应力，pd为出口外压。

稳态解确定模块503，用于利用有限元法，根据所述控制方程以及所述边界条件确定所述大变形管道流动系统的稳态解。

特征值确定模块504，用于根据所述稳态解确定特征值。

所述特征值确定模块504具体包括：

矩阵确定单元，用于根据所述稳态解确定所述大变形管道流动系统的质量矩阵以及非线性矩阵；

特征值确定单元，用于根据所述质量矩阵以及所述非线性矩阵确定特征值。

振荡空间确定模块505，用于根据所述特征值确定所述大变形管道流动系统的振荡空间的当前状态。

所述振荡空间确定模块505具体包括：

第一判断单元，用于判断所述特征值是否大于0，得到第一判断结果；

发散空间确定单元，用于若所述第一判断结果表示为所述特征值大于0，所述大变形管道流动系统的振荡空间的当前状态为发散空间；

稳定空间确定单元，用于若所述第一判断结果表示为所述特征值不大于0，所述大变形管道流动系统的振荡空间的当前状态为稳定空间。

本发明模拟管壁的变形部分的材料模型为不可压缩的超弹性模型，相比现有模型所采用的线弹性模型来说，本发明所提供的三维模型更贴近人体管壁，振荡空间当前状态的描述更为准确。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。