CN111914448A - 基于控制体有限元方法的流固耦合数值模拟方法 - Google Patents

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

Info

Publication number
CN111914448A
CN111914448A CN202010685058.XA CN202010685058A CN111914448A CN 111914448 A CN111914448 A CN 111914448A CN 202010685058 A CN202010685058 A CN 202010685058A CN 111914448 A CN111914448 A CN 111914448A
Authority
CN
China
Prior art keywords
unit
grid system
fluid
solid
field
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202010685058.XA
Other languages
English (en)
Other versions
CN111914448B (zh
Inventor
万义钊
陈强
欧阳伟平
董林
胡高伟
刘昌岭
吴能友
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Qingdao Institute of Marine Geology
Original Assignee
Qingdao Institute of Marine Geology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Qingdao Institute of Marine Geology filed Critical Qingdao Institute of Marine Geology
Priority to CN202010685058.XA priority Critical patent/CN111914448B/zh
Publication of CN111914448A publication Critical patent/CN111914448A/zh
Application granted granted Critical
Publication of CN111914448B publication Critical patent/CN111914448B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • General Engineering & Computer Science (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Fluid Mechanics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computing Systems (AREA)
  • Algebra (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开一种基于控制体有限元方法的流固耦合数值模拟方法：首先建立计算区域，并对计算区域进行网格离散，得到一套由一定数量的单元组成的网格系统，一个单元由若干节点组成；将上述网格系统中的每个单元以该单元的重心为基点为基础进行平均剖分，剖分的份数为单元的节点数量；一个节点周围的所有上述子单元组成一个控制体，所有的控制体组成对偶网格系统；接下来用网格系统的单元节点上的流体场的物理量计算固体场的物理量，基于有限元方法在网格系统上对固体变形方程进行离散求解；然后用网格系统的单元节点上的固体场的物理量计算流体场的物理量，基于有限体积法在对偶网格系统上对流体方程进行离散求解；最后，增加时间进行重复计算，进而得到流体场和固体场。

Description

基于控制体有限元方法的流固耦合数值模拟方法
技术领域
本发明涉及流固耦合数值模拟技术领域，具体涉及一种基于控制体有限元方法的流固耦 合数值模拟方法。
背景技术
流固耦合研究变形固体在流场作用下的各种行为及固体变形对流场影响的流固相互作用 问题。流固耦合在飞机叶片颤振、涡轮机械设计、海洋工程、动脉血热流动和变形多孔介质 流动广泛存在。总体上，流固耦合问题可以分为两大类：第一类是流体和固体的耦合作用仅 仅发生在流体和固体相的界面上，如飞机叶片颤振等问题。这类问题的方程的耦合是由流体 相和固体相界面上的平衡方程来实现；第二类是固体相和流体相的计算域部分或全部重叠在 一起，如变形多孔介质流动问题。
针对第二类流固耦合问题的数值模拟，由于流体域和固体域存在重合，在耦合计算时， 流体场和固体场的物理量在同一个网格节点上计算，通过相应的流固耦合关系进行耦合，即 固体场的物理量和流体场的物理量需要在同一位置耦合计算。然而，流体方程的模拟由于守 恒性的要求，通常采取有限体积等以单元为中心方法来计算；固体方程的模拟则采取以结点 为中心的有限元方法来实现。在两相耦合过程中，两个场的模拟采取同一套网格体系，则模 拟时流体场的节点在单元中心，固体场在结点上，造成了两个场之间需要通过插值的形式来 实现值的相互传递和耦合。
如在水合物开采的力学响应问题的流固耦合模拟中，水合物开采的流体场采用TOUGH+HYDRTAE软件实现，固体力学场的模拟使用Flac3D软件及TOUGH2Biot等模块， TOUGH+HYDRTAE软件使用有限体积法对含水合物沉积物中的多孔介质渗流、水合物分解 及传热场进行模拟，采取的是单元中心型的离散格式；Flac3D软件及HydrateBiot使用有限元 方法对固体力学场进行计算，采用的是节点中心型的离散格式。在固体变形的力场与多孔介质渗流、水合物分解及传热过程耦合时，使用插值的方法，利用单元中心的渗流场的孔隙压力、饱和度及温度值计算单元节点上的值，再与节点上的力学场耦合(Lei,H.；Xu,T.；Jin,G. TOUGH2Biot–A Simulator for coupled Thermal–Hydrodynamic–MechanicalProcesses in subsurface Flow systems:Application to CO2 Geological storageand Geothermal development. Computers&geosciences 2015.)。
然而，在同一套网格中使用插值引入了额外的插值计算，降低了效率，且插值是近似计 算，存在一定的误差，为了解决第二类流固耦合问题中存在的插值问题，亟待提出一种新的 流固耦合模拟方法。
发明内容
本发明针对多孔介质流固耦合问题，提出一种基于控制体有限元方法的流固耦合数值模 拟方法，该方法可以直接将流体场的计算值和固体场的计算值放在同一节点上，避免了耦合 中的插值计算，有效提高了效率。
本发明是采用以下的技术方案实现的：一种基于控制体有限元方法的流固耦合模拟方法， 包括以下步骤：
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.采用有限元方法在单元上对固体场方程进行离散，形成单元刚度方程；所述的固体 场方程以弹性本构假设下的位移形式为例：
Figure BDA0002587239240000051
式中，
Figure BDA0002587239240000052
A3＝G，G为剪切模量，ν为泊松比；wx、wy分别为x方向和y方向的固体变形位移。
以公式(10)中的第一个方程为例，说明其在单元上的有限元离散。
将方程两端同乘以位移的变分δwx,然后在单元上积分并降阶可得：
Figure RE-GDA0002665828090000061
式中，S为三角形单元的面积，s为三角形单元的三条边；
Figure BDA0002587239240000062
lx，ly为三角形单元边的方向余弦。
将单元插值形函数代入后可得单元刚度方程：
[k11]{wx}+[k12]{wy}+[k13]{wz}+[k14]{p}＝{f1} (12)
其中，
Figure BDA0002587239240000063
Figure BDA0002587239240000064
Figure BDA0002587239240000065
Figure BDA0002587239240000066
{f1}＝-∫(FxNj)ds
S34.在所有的网格中重复上述步骤S31～S33,并将所有单元的单元刚度方程集合形成总 体刚度方程；
S35.求解总体刚度方程得到网格系统节点上的固体场的位移，利用位移计算应力、应变 等其他固体场的物理量；
S4.用网格系统的单元节点上的固体场的物理量计算流体场的物理量，采用有限体积法 在对偶网格系统上对流体方程进行离散求解；
具体处理步骤为：
S41.利用耦合关系，在节点上用步骤S3求得的固体场的物理量计算流场的物理量；所述 的固体场的物理量以有效应力σ′为例，所述流体场的物理量以渗透率K为例，则使用如下关 系计算固体场对流体场流动的影响：
Kσ/K＝cσ′2+dσ′+e (13)
式中，Kσ为应力状态下的渗透率，c,d,e为实验测试的拟合参数；
S42.将上述流场的物理量代入流场的方程，所述的流场的方程以多孔介质单相渗流的压 力方程为例：
Figure BDA0002587239240000071
式中，φ为孔隙度，μ为流体黏度，ct为综合压缩系数，t为时间；
S43.采用有限体积法在控制体上对流场方程进行离散，得到控制体离散方程；在图3所 示的控制体上对方程(14)积分，
Figure BDA0002587239240000072
式中，Sc为图3所示的控制体9面积；
上式的左端使用高斯散度定理可得：
Figure BDA0002587239240000073
式中，
Figure BDA0002587239240000074
为图4控制体的边，即图4中点5、7、6、10、11、12、13、14、15、16组成的线段。
上式右端项表示从边界上流入控制体的净流量，可以使用子单元边界上的流量相加，即
Figure BDA0002587239240000075
L5-6-7表示点5、6、7组成的线段，以此类推。
上式中的压力梯度则使用子单元所在的单元的形函数来计算，以线段5-6-7为例，压力 梯度为：
Figure BDA0002587239240000076
上式中N为单元的插值形函数，Ni＝ai+bix+ciy，
Figure BDA0002587239240000077
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)积分，
Figure BDA0002587239240000091
式中，Vc为图6所示的控制体9体积；
上式的左端使用高斯散度定理可得：
Figure BDA0002587239240000092
式中，
Figure BDA0002587239240000093
为图6控制体的边界面。
上式右端项表示从边界上流入控制体的净流量，可以使用子单元边界上的流量相加，即 将上述积分分解为构成控制体的子单元的面上的积分，其中图6所示的子单元se的积分又可 以写为：
Figure BDA0002587239240000094
S44.重复上述步骤S41～S44,计算所有控制体的离散方程，将所有控制体的方程集成线性 方程组；
S45.求解线性方程组，得到节点上的流场的压力值等物理量，如图7所示，为计算得到 的井筒附近的流畅压力等值线图。
S5.增加时间，重复S3和S4，计算每个时刻的流体场合固体场的物理量，如图8所示为 计算得到井筒固体场的应力等值线图。
通过上述分析，本发明通过在原网格系统中构建了对偶网格系统，原网格系统的节点为 对偶网格系统的中心点；固体场的计算在原网格系统上实现，而流场的计算在对偶网格系统 上开展，可以保证流体场合固体场的计算在同一个点上，因此，流体场和固体场的耦合计算 不需要通过节点插值，避免了插值计算带来的计算效率的损失和误差。
以上所述，仅是本发明的较佳实施例而已，并非是对本发明作其它形式的限制，任何熟 悉本专业的技术人员可能利用上述揭示的技术内容加以变更或改型为等同变化的等效实施例 应用于其它领域，但是凡是未脱离本发明技术方案内容，依据本发明的技术实质对以上实施 例所作的任何简单修改、等同变化与改型，仍属于本发明技术方案的保护范围。

Claims (5)

1.基于控制体有限元方法的流固耦合模拟方法，其特征在于，包括以下步骤：
S1.建立计算区域，对计算区域进行网格离散，得到一套由一定数量的单元组成的网格系统，一个单元由若干节点组成；
S2.将上述网格系统中的每个单元以该单元的重心为基点进行平均剖分，剖分的份数为单元的节点数量，单元剖分的每一份称为一个子单元；一个节点周围的所有上述子单元组成一个控制体，所有的控制体组合在一起组成对偶网格系统；
S3.用网格系统的单元节点上的流体场的物理量计算固体场的物理量，以节点中心型的有限元方法在网格系统上对固体变形方程进行离散求解；
S4.用网格系统的单元节点上的固体场的物理量计算流体场的物理量，采用单元中心型的有限体积法在对偶网格系统上对流体方程进行离散求解；
S5.增加时间，重复S3和S4，进而得到流体场和固体场。
2.根据权利要求1所述的基于控制体有限元方法的流固耦合模拟方法，其特征在于：所述步骤S1具体采用以下方式实现：
S11.根据计算的模型的几何尺寸建立几何模型；
S12.将几何模型划分为成若干单元，一个单元由若干节点和若干条边组成，所有的单元组成网格系统。
3.根据权利要求1所述的基于控制体有限元方法的流固耦合模拟方法，其特征在于：所述步骤S2包括：
S21.计算所述网格系统中每个单元的重心和单元边中点；
S22.连接网格系统中每个单元的重心和单元边中点，将所有单元分别平均剖分；单元的一个节点、与这个节点相邻的所有的边中点和重心组成一个子单元；
S23.将一个节点周围所有的子单元合并在一起组成一个控制体，所有的控制体组合在一起形成对偶网格系统。
4.根据权利要求1所述的基于控制体有限元方法的流固耦合模拟方法，其特征在于：所述步骤S3包括：
S31.基于流固耦合关系，由节点上的流体场的物理量计算固体场的物理量；
S32.采用有限元方法在网格系统上对固体场方程进行离散，形成单元刚度方程；
S33.将所有单元的单元刚度方程集合形成总体刚度方程；
S34.求解总体刚度方程得到节点上的固体场的物理量。
5.根据权利要求1所述的基于控制体有限元方法的流固耦合模拟方法，其特征在于：所述步骤4包括：
S41.利用耦合关系，在节点上用步骤S3求得的固体场的物理量计算流场的物理量；
S42.采用有限体积法在由控制体形成的对偶网格系统上对流场方程进行离散，得到控制体离散方程；
S43.将所有控制体的方程集成线性方程组；
S44.求解线性方程组，得到节点上的流场的物理量。
CN202010685058.XA 2020-07-16 2020-07-16 基于控制体有限元方法的流固耦合数值模拟方法 Active CN111914448B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010685058.XA CN111914448B (zh) 2020-07-16 2020-07-16 基于控制体有限元方法的流固耦合数值模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010685058.XA CN111914448B (zh) 2020-07-16 2020-07-16 基于控制体有限元方法的流固耦合数值模拟方法

Publications (2)

Publication Number Publication Date
CN111914448A true CN111914448A (zh) 2020-11-10
CN111914448B CN111914448B (zh) 2023-06-09

Family

ID=73281801

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010685058.XA Active CN111914448B (zh) 2020-07-16 2020-07-16 基于控制体有限元方法的流固耦合数值模拟方法

Country Status (1)

Country Link
CN (1) CN111914448B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112800700A (zh) * 2021-04-13 2021-05-14 中国空气动力研究与发展中心计算空气动力研究所 低温表面干模态结霜模拟方法、装置、电子设备和介质
CN113642217A (zh) * 2021-08-17 2021-11-12 王永亮 一种多孔弹性岩体热-流-固耦合压裂裂缝扩展模拟方法
CN115688212A (zh) * 2022-11-15 2023-02-03 中国科学院深圳先进技术研究院 一种基于物质点法的软体机器人仿真方法
CN115995277A (zh) * 2023-03-22 2023-04-21 中国空气动力研究与发展中心计算空气动力研究所 一种材料动力学特性评估方法、装置、设备及介质

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20060235667A1 (en) * 2005-04-14 2006-10-19 Fung Larry S Solution method and apparatus for large-scale simulation of layered formations
CN102224502A (zh) * 2008-10-09 2011-10-19 雪佛龙美国公司 用于多孔介质中的流动的迭代多尺度方法
CN108241777A (zh) * 2017-12-27 2018-07-03 青岛海洋地质研究所 基于非结构网格有限元法计算水合物沉积物中渗流速度场的方法
CN108920768A (zh) * 2018-06-07 2018-11-30 天津大学 一种针对弹性薄壁结构的流固耦合算法
CN111339702A (zh) * 2020-02-26 2020-06-26 青岛海洋地质研究所 油藏数值模拟等效井筒半径计算方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20060235667A1 (en) * 2005-04-14 2006-10-19 Fung Larry S Solution method and apparatus for large-scale simulation of layered formations
CN102224502A (zh) * 2008-10-09 2011-10-19 雪佛龙美国公司 用于多孔介质中的流动的迭代多尺度方法
CN108241777A (zh) * 2017-12-27 2018-07-03 青岛海洋地质研究所 基于非结构网格有限元法计算水合物沉积物中渗流速度场的方法
CN108920768A (zh) * 2018-06-07 2018-11-30 天津大学 一种针对弹性薄壁结构的流固耦合算法
CN111339702A (zh) * 2020-02-26 2020-06-26 青岛海洋地质研究所 油藏数值模拟等效井筒半径计算方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
PEIXIAO MAO ET AL.: ""Effect of permeability anisotropy on depressurization‐induced gas production from hydrate reservoirs in the South China Sea"", 《ENERGY SCIENCE AND ENGINEERING》 *
万义钊等: ""南海神狐海域天然气水合物降压开采过程中储层的稳定性"", 《开发工程》 *
杨军征: ""有限体积—有限元方法在油藏数值模拟中的原理和应用"", 《中国博士学位论文全文数据库基础科学辑》 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112800700A (zh) * 2021-04-13 2021-05-14 中国空气动力研究与发展中心计算空气动力研究所 低温表面干模态结霜模拟方法、装置、电子设备和介质
CN112800700B (zh) * 2021-04-13 2021-06-25 中国空气动力研究与发展中心计算空气动力研究所 低温表面干模态结霜模拟方法、装置、电子设备和介质
CN113642217A (zh) * 2021-08-17 2021-11-12 王永亮 一种多孔弹性岩体热-流-固耦合压裂裂缝扩展模拟方法
CN115688212A (zh) * 2022-11-15 2023-02-03 中国科学院深圳先进技术研究院 一种基于物质点法的软体机器人仿真方法
CN115995277A (zh) * 2023-03-22 2023-04-21 中国空气动力研究与发展中心计算空气动力研究所 一种材料动力学特性评估方法、装置、设备及介质

Also Published As

Publication number Publication date
CN111914448B (zh) 2023-06-09

Similar Documents

Publication Publication Date Title
CN111914448A (zh) 基于控制体有限元方法的流固耦合数值模拟方法
Martini et al. Reduced basis approximation and a-posteriori error estimation for the coupled Stokes-Darcy system
CN111859766A (zh) 可变计算域的拉格朗日积分点有限元数值仿真系统及方法
Spekreijse et al. A simple, robust and fast algorithm to compute deformations of multi-block structured grids
CN111125963B (zh) 基于拉格朗日积分点有限元的数值仿真系统及方法
CN109446649B (zh) 致密油藏体积压裂水平井三维渗流模型的建立方法
CN115659908B (zh) 一种印刷电路板换热器的三单元非平衡多孔介质方法
CN115510717A (zh) 一种基于有限差分方法的高精度激波捕捉方法
CN112001109A (zh) 再生核粒子算法实现结构冲击动力学仿真方法
CN109657408A (zh) 一种再生核粒子算法实现结构线性静力学仿真方法
CN109063239B (zh) 一种水热耦合的三维数值模拟方法
CN111241728A (zh) 一种欧拉方程的间断伽辽金有限元数值求解方法
CN117786966A (zh) 一种基于lbfs的裂隙岩体多场耦合数值计算方法
Karahan et al. Time‐dependent groundwater modeling using spreadsheet
Peraire et al. An adaptive finite element method for high speed flows
Cheng et al. A hybrid reconstructed discontinuous Galerkin method for compressible flows on unstructured grids
Niewiarowski et al. Pneumatic storm surge barriers subject to hydrodynamic loading
DOLEJŠÍ et al. Anisotropic mesh adaptation for transonic and supersonic flow simulation
CN116956670B (zh) 一种基于tpfa与mfd混合方法的投影嵌入式离散裂缝模型
Fernandez Implicit conservative upwind schemes for strongly transient flows
Yagawa Free Mesh Method: fundamental conception, algorithms and accuracy study
Bruchon et al. Direct 2D simulation of small gas bubble clusters: From the expansion step to the equilibrium state
Cavendish et al. SOLUTION OF INCOMPRESSIBLE NAVIER‐STOKES EQUATIONS ON UNSTRUCTURED GRIDS USING DUAL TESSELLATIONS
CN116205187B (zh) 多对称声表面波器件的快速仿真方法、系统及相关设备
Johnston et al. An adaptive indirect boundary element method with applications

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant