基于控制体有限元方法的流固耦合数值模拟方法与流程

本发明涉及流固耦合数值模拟技术领域，具体涉及一种基于控制体有限元方法的流固耦合数值模拟方法。

背景技术：

流固耦合研究变形固体在流场作用下的各种行为及固体变形对流场影响的流固相互作用问题。流固耦合在飞机叶片颤振、涡轮机械设计、海洋工程、动脉血热流动和变形多孔介质流动广泛存在。总体上，流固耦合问题可以分为两大类：第一类是流体和固体的耦合作用仅仅发生在流体和固体相的界面上，如飞机叶片颤振等问题。这类问题的方程的耦合是由流体相和固体相界面上的平衡方程来实现；第二类是固体相和流体相的计算域部分或全部重叠在一起，如变形多孔介质流动问题。

针对第二类流固耦合问题的数值模拟，由于流体域和固体域存在重合，在耦合计算时，流体场和固体场的物理量在同一个网格节点上计算，通过相应的流固耦合关系进行耦合，即固体场的物理量和流体场的物理量需要在同一位置耦合计算。然而，流体方程的模拟由于守恒性的要求，通常采取有限体积等以单元为中心方法来计算；固体方程的模拟则采取以结点为中心的有限元方法来实现。在两相耦合过程中，两个场的模拟采取同一套网格体系，则模拟时流体场的节点在单元中心，固体场在结点上，造成了两个场之间需要通过插值的形式来实现值的相互传递和耦合。

如在水合物开采的力学响应问题的流固耦合模拟中，水合物开采的流体场采用tough+hydrtae软件实现，固体力学场的模拟使用flac3d软件及tough2biot等模块，tough+hydrtae软件使用有限体积法对含水合物沉积物中的多孔介质渗流、水合物分解及传热场进行模拟，采取的是单元中心型的离散格式；flac3d软件及hydratebiot使用有限元方法对固体力学场进行计算，采用的是节点中心型的离散格式。在固体变形的力场与多孔介质渗流、水合物分解及传热过程耦合时，使用插值的方法，利用单元中心的渗流场的孔隙压力、饱和度及温度值计算单元节点上的值，再与节点上的力学场耦合(lei,h.；xu,t.；jin,g.tough2biot–asimulatorforcoupledthermal–hydrodynamic–mechanicalprocessesinsubsurfaceflowsystems:applicationtoco2geologicalstorageandgeothermaldevelopment.computers&geosciences2015.)。

然而，在同一套网格中使用插值引入了额外的插值计算，降低了效率，且插值是近似计算，存在一定的误差，为了解决第二类流固耦合问题中存在的插值问题，亟待提出一种新的流固耦合模拟方法。

技术实现要素：

本发明针对多孔介质流固耦合问题，提出一种基于控制体有限元方法的流固耦合数值模拟方法，该方法可以直接将流体场的计算值和固体场的计算值放在同一节点上，避免了耦合中的插值计算，有效提高了效率。

本发明是采用以下的技术方案实现的：一种基于控制体有限元方法的流固耦合模拟方法，包括以下步骤：

s1.建立计算区域，对计算区域进行网格离散，得到一套由一定数量的单元组成的网格系统，一个单元由若干节点组成；

s2.将上述网格系统中的每个单元以该单元的重心为基础进行平均剖分，剖分的份数为单元的节点数量，单元剖分的每一份称为一个子单元；一个节点周围的所有上述子单元组成一个控制体，所有的控制体组合在一起组成对偶网格系统；

s3.用网格系统的单元节点上的流体场的物理量计算固体场的物理量，基于以节点中心型的有限元方法在网格系统上对固体变形方程进行离散求解；

s4.用网格系统的单元节点上的固体场的物理量计算流体场的物理量，采用单元中心型的有限体积法在对偶网格系统上对流体方程进行离散求解；

s5.增加时间，重复s3和s4，进而得到流体场和固体场。

进一步的，所述步骤s1具体采用以下方式实现：

s11.根据计算的模型的几何尺寸建立几何模型；

s12.将几何模型划分为成若干单元，一个单元由若干节点和若干条边组成，所有的单元组成网格系统。

进一步的，所述步骤s2包括：

s21.计算上述网格系统中每个单元的重心和单元边中点；

s22.连接网格系统中每个单元的重心和边中点，将单元平均剖分；由单元的一个节点、与这个节点相邻的所有的边中点和重心组成一个子单元，将网格中所有的单元都进行上述剖分；

s23.将一个节点周围所有的子单元合并在一起组成一个控制体，所有的上述控制体组合在一起形成对偶网格系统。

进一步的，所述步骤s3包括：

s31.利用流固耦合关系，由节点上的流体场的物理量计算固体场的物理量；

s32.采用有限元方法在单元上对固体场方程进行离散，形成单元刚度方程；

s33.将所有单元的单元刚度方程集合形成总体刚度方程；

s34.求解总体刚度方程得到节点上的固体场的物理量。

进一步的，所述步骤4包括：

s41.利用耦合关系，在节点上用步骤s3求得的固体场的物理量计算流场的物理量；

s42.采用有限体积法在控制体上对流场方程进行离散，得到控制体离散方程；

s43.将所有控制体的方程集成线性方程组；

s44.求解线性方程组，得到节点上的流场的物理量。

与现有技术相比，本发明的优点和积极效果在于：

本方案所提出的流固耦合数值模拟方法，通过对网格系统中的每个单元以该单元的重心为基础进行平均剖分，剖分的每一份称为一个子单元；并将一个节点周围的所有上述子单元组成的区间称为一个控制体，基于控制体有限元方法进行具体的分析；实现在原网格系统中构建对偶网格系统，原网格系统的节点为对偶网格系统的中心点；固体场的计算在原网格系统上实现，而流场的计算在对偶网格系统上开展，这样可以保证流体场合固体场的计算在同一个点上，流体场和固体场的耦合计算不需要通过节点插值，避免了插值计算带来的计算效率的损失和误差，有效提高计算效率和精度，具有更好的实际应用价值。

附图说明

图1为本发明实施例1所述模型和网格系统示意图；

图2为本发明实施例1所述子单元构建示意图；

图3为本发明实施例1所述控制体构建示意图；

图4为本发明实施例1所构建的网格系统和对偶网格系统示意图；

图5为本发明实施例2所述三维模型的四面体网格系统示意图；

图6为本发明实施例2所述的四面体单元及子单元构建示意图；

图7为本发明实施例2计算得到的流场压力场等值线示意图；

图8为本发明实施例2计算得到的固体场的应力等值线示意图；

图9为本发明所述流固耦合数值模拟方法的流程示意图；

其中，1、2、3为单元的三个结点；4、5、6为单元的边中点；7为单元的重心；8为子单元；9为控制体。

具体实施方式

为了能够更加清楚地理解本发明的上述目的、特征和优点，下面结合附图及实施例对本发明做进一步说明。在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是，本发明还可以采用不同于在此描述的其他方式来实施，因此，本发明并不限于下面公开的具体实施例。

本发明公开一种基于控制体有限元方法的流固耦合模拟方法，如图9所示，包括以下步骤：

s1.建立计算区域，对计算区域进行网格离散，得到一套由一定数量的单元组成的网格系统，一个单元包括若干节点；

s2.将所述网格系统中的每个单元以该单元的重心为基础进行平均剖分，剖分的份数为该单元的节点数量，单元剖分的每一份称为一个子单元；一个节点周围的所有上述子单元组成一个控制体，所述控制体组成对偶网格系统；

s3.用网格系统的单元节点上的流体场的物理量计算固体场的物理量，采用有限元方法在网格系统上对固体变形方程进行离散求解；

s4.用网格系统的单元节点上的固体场的物理量计算流体场的物理量，采用有限体积法在对偶网格系统上对流体方程进行离散求解；

s5.增加时间，重复步骤s3和s4，得到流体场和固体场。

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合二维三角形单元和三维四面体单元两个实施实例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

实例1：以二维三角形为例，具体的：

s1.建立计算区域，对计算区域进行网格离散，得到一套由一定数量的单元组成的网格系统，所述的一个单元由若干节点组成；

具体处理的步骤为：

s11.根据计算的模型的几何尺寸建立几何模型，设模型为长300米、宽50米的二维模型；

s12.将上述二维的几何模型划分为网格系统，所述网格系统由1339个单元和740个节点组成，所述单元为三角形单元，所述的每个三角形单元由三个节点和三条边组成，如图1所示；

s2.将上述网格系统中的每个单元用该单元的重心进行平均剖分，剖分的份数为单元的节点数量，单元剖分的每一份称为一个子单元；一个节点周围的所有上述子单元组成一个控制体，上述控制体组成对偶网格系统；

具体处理步骤为：

s21.取一个三角形单元，所述三角形的三个节点的编号分别为1、2、3，三角形的三条边的中点分别为4、5、6，三角形单元的重心为7。如图2所示。所述三个结点的坐标分别为(x1,y1),(x2,y2),(x3,y3)。所述的三角形单元三条边的三个中点4、5、6采用如下公式计算：

x4＝(x1+x2)/2(1)

y4＝(y1+y2)/2(2)

x5＝(x2+x3)/2(3)

y5＝(y2+y3)/2(4)

x6＝(x1+x3)/2(5)

y6＝(y1+y3)/2(6)

所述三角形单元的重心点7的坐标采用如下公式计算：

x4＝(x1+x2+x3)/3(7)

y4＝(y1+y2+y3)/3(8)

s22.将单元的重心7与三角形单元三条边的中点4、5、6分别连接，如图2所示；

s23.在一个三角形单元内，由节点3、节点3相邻的两个边中点5、6和重心7相连围成子单元，如图2所示；同理，节点2与边中点4、5和重心7围成一个子单元；节点1与边中点4、6和重心7围成一个子单元；一个单元分为三个子单元；

s24.对网格系统中的1339个单元重复上述步骤s21～s23，可以得到1339×3＝4017个子单元；

s25.在每个单元都进行上述剖分后，一个节点周围存在若干个子单元，一个节点周围的所有子单元合并在一起组成一个控制体，如图3所示，节点3周围有5个子单元，分别为子单元3-6-7-5、子单元3-5-15-15、子单元3-15-14-13、子单元3-13-12-11、子单元3-11-10-6，这5个子单元合并形成一个由5、7、6、10、11、12、13、14、15、16围成的控制体9，在上述网格系统中可以形成740个控制体，这740个控制体形成对偶网格系统，如图4所示；

s3.用网格系统的单元节点上的流体场的物理量计算固体场的物理量，采用有限元方法在网格系统上对固体变形方程进行离散求解；

具体处理步骤为：

s31.利用流固耦合关系，由节点1、2、3上的流体场的物理量计算固体场的物理量；所述的流体场的物理量以压力p为例，所述固体场的物理量以有效应力σ′为例，则计算关系式为：

σ′i＝σi-αpi(9)

式中，i为三角形单元的三个节点,取值分别为1、2、3；σ为总应力；α为biot系数。

s32.将计算得到的固体场物理量代入固体场的方程；

s33.采用有限元方法在单元上对固体场方程进行离散，形成单元刚度方程；所述的固体场方程以弹性本构假设下的位移形式为例：

式中，a3＝g，g为剪切模量，ν为泊松比；wx、wy分别为x方向和y方向的固体变形位移。

以公式(10)中的第一个方程为例，说明其在单元上的有限元离散。

将方程两端同乘以位移的变分δwx,然后在单元上积分并降阶可得：

式中，s为三角形单元的面积，s为三角形单元的三条边；

lx，ly为三角形单元边的方向余弦。

将单元插值形函数代入后可得单元刚度方程：

[k11]{wx}+[k12]{wy}+[k13]{wz}+[k14]{p}＝{f1}(12)

其中，

{f1}＝-∫(fxnj)ds

s34.在所有的网格中重复上述步骤s31～s33,并将所有单元的单元刚度方程集合形成总体刚度方程；

s35.求解总体刚度方程得到网格系统节点上的固体场的位移，利用位移计算应力、应变等其他固体场的物理量；

s4.用网格系统的单元节点上的固体场的物理量计算流体场的物理量，采用有限体积法在对偶网格系统上对流体方程进行离散求解；

具体处理步骤为：

s41.利用耦合关系，在节点上用步骤s3求得的固体场的物理量计算流场的物理量；所述的固体场的物理量以有效应力σ′为例，所述流体场的物理量以渗透率k为例，则使用如下关系计算固体场对流体场流动的影响：

kσ/k＝cσ′²+dσ′+e(13)

式中，kσ为应力状态下的渗透率，c,d,e为实验测试的拟合参数；

s42.将上述流场的物理量代入流场的方程，所述的流场的方程以多孔介质单相渗流的压力方程为例：

式中，φ为孔隙度，μ为流体黏度，ct为综合压缩系数，t为时间；

s43.采用有限体积法在控制体上对流场方程进行离散，得到控制体离散方程；在图3所示的控制体上对方程(14)积分，

式中，sc为图3所示的控制体9面积；

上式的左端使用高斯散度定理可得：

式中，为图4控制体的边，即图4中点5、7、6、10、11、12、13、14、15、16组成的线段。

上式右端项表示从边界上流入控制体的净流量，可以使用子单元边界上的流量相加，即

l5-6-7表示点5、6、7组成的线段，以此类推。

上式中的压力梯度则使用子单元所在的单元的形函数来计算，以线段5-6-7为例，压力梯度为：

上式中n为单元的插值形函数，ni＝ai+bix+ciy，i,j,k为单元的三个节点，即图3所示的1，2，3。

s44.重复上述步骤s41～s44，计算所有控制体的离散方程，将所有控制体的方程集成线性方程组；

s45.求解线性方程组，得到节点上的流场的压力值等物理量；

s5.增加时间，重复s3和s4，计算每个时刻的流体场合固体场的物理量。

本方案在一般的网格系统中，通过连接单元中心的方式构建了针对流场方程计算的控制体，这些控制体组成的对偶网格系统所表示的计算区域与原网格系统表示的计算区域完全相同，而控制体的中心点则为原网格的节点。构建这样的对偶网格系统后，固体场的计算仍然采用以节点中心型的有限元方法在原网格系统(剖分之前)的单元上进行，求得的固体场的物理量位于原网格的节点上；而流体场的计算则采用单元中心型的有限体积法计算，计算得到的流体场的物理量位于构建的控制体的中心上，也即原网格系统的节点上。构建的控制体确保了流体场的计算结果和固体场的计算结果位于同一个点上，因此，流体场和固体场的物理量之间的耦合计算不需要通过节点插值，避免了插值计算带来的计算效率的损失和误差。

需要说明的是，上述详细说明中，使用了二维三角形单元进行计算，但本发明方法不限定在二维三角形单元上，对于二维和三维任意网格系统都适用；同时，上述方法也适用于多相渗流的情况。

实例2：以三维问题的四面体网格为例说明本发明在三维情况下的具体实施过程：

s1.建立计算区域，对计算区域进行网格离散，得到一套由一定数量的单元组成的网格系统，所述的一个单元由若干节点组成；

具体处理的步骤为：

s11.根据计算的模型的几何尺寸建立几何模型；设模型半径为2米、厚度为0.1米的三维圆柱模型，模型中间被半径为0.2米、厚度为0.1米的圆柱掏空，如图5所示；该模型模拟油气开采时井壁附近的流固耦合问题；

s12.将上述三维的几何模型划分为网格系统，所述网格系统由114949个单元和27317个节点组成，如图5所示；所述单元为四面体单元，所述的四面体单元由四个节点、四个面和六条边组成，如图6所示，构成四面体单元的四个节点为p1、p2、p3、p4，四面体单元四个面分别由p1-p2-p3、p1-p2-p4、p1-p3-p4和p2-p3-p4，四面体单元的六条边分别为p1-p2、p1-p3、p1-p4、p2-p3、p2-p4、p3-p4；

s2.将上述网格系统中的每个单元用该单元的重心进行平均剖分，剖分的份数为单元的节点数量，单元剖分的每一份称为一个子单元；一个节点周围的所有上述子单元组成一个控制体，上述控制体组成对偶网格系统；

具体处理步骤为：

s21.取一个四面体单元，所述四面体单元的四个节点的编号分别为p1、p2、p3、p4，所述四面体单元的面的中心点用fijk表示，其中ijk为构成该面的三个节点的编号，如图6所示，面p1-p2-p3的中心点为f123，面p1-p2-p4的中心点为f124，依次类推。所述四面体的边的中点用eij表示，其中ij为构成该条边的两个节点的编号，如图6所示，边p1-p2的中点为e12,边p1-p3的中点为e13,依次类推；所述四面体单元的重心为c；

s22.将四面体单元的重心与单元的四个面的中心点相连，将面的中心点与该面的三条边的中点相连，如图6所示；

s23.在一个四面体单元内，在节点p4周围构建子单元sub4，所述子单元sub4由四个面组成：sf24、sf34、sc34、sc24；所述的面sf24由p4、e14、f134、e34构成；所述的面sf34由节点p4、e14、f124、e24构成；所述的面sc34由e34、f234、c、f134构成；所述的面sc24由e24、f124、c、f234构成；一个四面体单元可以分为四个子单元；

s24.对网格系统中的114949个单元重复上述步骤s21～s23，可以得到114949×4＝459796个子单元；

s24.在每个单元都进行上述剖分后，一个节点周围存在若干个子单元，一个节点周围的所有子单元合并在一起组成一个控制体，如图6所示，节点p4周围有8个子单元，这8个子单元合并形成一个控制体9。在上述网格系统中可以形成27317个控制体，这27317个控制体形成对偶网格系统。

s3.用网格系统的单元节点上的流体场的物理量计算固体场的物理量，采用有限元方法在网格系统上对固体变形方程进行离散求解；

具体处理步骤与实例1中的s31～s35原理相似，在此不做赘述；

s4.用网格系统的单元节点上的固体场的物理量计算流体场的物理量，采用有限体积法在对偶网格系统上对流体方程进行离散求解；

具体处理步骤中的前三步与实例1中的s41～s43原理相似，在此不做赘述；

s43.采用有限体积法在控制体上对流场方程进行离散，得到控制体离散方程；在图3所示的控制体上对方程(14)积分，

式中，vc为图6所示的控制体9体积；

上式的左端使用高斯散度定理可得：

式中，为图6控制体的边界面。

上式右端项表示从边界上流入控制体的净流量，可以使用子单元边界上的流量相加，即将上述积分分解为构成控制体的子单元的面上的积分，其中图6所示的子单元se的积分又可以写为：

s44.重复上述步骤s41～s44,计算所有控制体的离散方程，将所有控制体的方程集成线性方程组；

s45.求解线性方程组，得到节点上的流场的压力值等物理量，如图7所示，为计算得到的井筒附近的流畅压力等值线图。

s5.增加时间，重复s3和s4，计算每个时刻的流体场合固体场的物理量，如图8所示为计算得到井筒固体场的应力等值线图。

通过上述分析，本发明通过在原网格系统中构建了对偶网格系统，原网格系统的节点为对偶网格系统的中心点；固体场的计算在原网格系统上实现，而流场的计算在对偶网格系统上开展，可以保证流体场合固体场的计算在同一个点上，因此，流体场和固体场的耦合计算不需要通过节点插值，避免了插值计算带来的计算效率的损失和误差。

以上所述，仅是本发明的较佳实施例而已，并非是对本发明作其它形式的限制，任何熟悉本专业的技术人员可能利用上述揭示的技术内容加以变更或改型为等同变化的等效实施例应用于其它领域，但是凡是未脱离本发明技术方案内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与改型，仍属于本发明技术方案的保护范围。