一种界面双向数据交换的流固耦合计算方法与流程

本发明属于流固耦合计算方法技术领域，具体涉及一种界面双向数据交换的流固耦合计算方法。

背景技术：

流固耦合作用是自然界客观存在的一种特殊现象，是指流体与固体之间的相互作用，在自然界随处可见，但是流体域与固体域均不可单独求解，因此要解决这个问题，需要将流固耦合然后进行求解。流固耦合力学是流体力学与固体力学及交叉形成的一门力学分支，是研究流体固体两种介质相互作用的一种具体手段。目前国内外针对流体域-固体域耦合在一起进行计算的研究均有不同程度的优势和缺陷，传统计算的做法是分别计算水体与固体的作用，在水体域计算中固体域的边界被设置为刚体，在固体域计算中流体域的作用通过简化的荷载形式进行施加，由于水体域与固体域在计算中会相互影响，故单单考虑水体或固体域结构并不能满足实际需求，为了解决这种矛盾，就必须想办法将流体域和固体域相结合，考虑它们之间相互作用的影响，因此对于流固耦合问题的研究势在必行。

技术实现要素：

针对目前流固耦合方法存在的缺陷和问题，本发明提供一种界面双向数据交换的流固耦合计算方法，该方法通过双向耦合，将流体力场与固定基础结构承载能力紧密联系起来，为研究结构损伤发展及变化的规律，以及耦联系统结构的流固耦合作用机理与损伤破坏模式提供技术支持。

本发明解决其技术问题所采用的方案是：一种界面双向数据交换的流固耦合计算方法，包括以下步骤：

step1、建立固体域模型：固体域计算结构动力学采用商用软件abaqus的隐式算法，建立csd固体域有限元模型。

step2、对固体域模型进行有限元网格划分，再设置固体域单元类型和材料属性，进一步对固体域模型分别施加约束条件和施加载荷条件进行结构分析，使其与实际条件相符。

step3、建立流体域模型：计算流体动力学则采用fluent平台进行计算，形成cfd流体分体模型。

step4、对流体域模型进行流体域网格划分，并对网格质量检查，网格质量合格时，对边界条件和初始条件进行设置；网格质量不合格时继续对网格质量进行检查，然后通过动网格技术进行网格重构并对流场进行模拟，使其与实际条件相符。

step5、生成动网格：分别将step1中的csd固体域有限元模型和step2中的cfd流体分体模型通过mpcci联合起来进行计算生成动网格。

step6、建立流固耦合的数据交换模型：通过流固耦合数据交换类型将流场模拟、动网格技术折合成固体域的施加荷载，并施加在固体域上；数据交换是通过设置流固交界面来编辑水体与渠道面板之间的相互作用，动网格间数据交换是在数据传递过程中采用双向耦合的方法，所述动网格间流固耦合过程分为三步：

1、预接触搜索：通过对流体域与固体域交界面处的单元进行分类，分别记录源网格和目标网格上的单元信息，为下一步的网格联系做准备；

2、网格联系：流体域与固体域交界面处的网格是不匹配的，流体域网格要想与固体域网格进行数据传递，每一组源网格均需与对应的一组目标网格建立联系从而为插值计算做准备；

3、通过插值函数进行插值计算：通过插值算法进行耦合界面间的数据交换，插值算法需要同时满足位移连续性条件和动力连续性条件；

a位移连续性条件：uf＝us

b动力连续性条件：σsns＝σfnf

式中：uf、us分别表示流体和固体位移矢量；us分别表示结构与流体的柯西应力张量；σs、σf分别表述结构与流体的外法线方向矢量。

step7、通过以上各步骤将双向流固耦合视为整体的耦联体系结构，对该耦联体系结构进行分析得出耦连体系损伤的体制的结论。

上述的一种界面双向数据交换的流固耦合计算方法，step2中对固体域模型进行有限元网格划分还包括对固体域模型进行cdp本构模型的设置。

上述的一种界面双向数据交换的流固耦合计算方法，所述cdp本构模型指的是结构损伤机制。

上述的一种界面双向数据交换的流固耦合计算方法，step3中还包括在建立流体域模型时，选择合适的湍流模型，所述湍流模型为零方程模型、单方程模型、两方程模型、雷诺应力模型和大涡模拟模型的一种或几种。

上述的一种界面双向数据交换的流固耦合计算方法，step6中所述预接触搜索过程包括：

(1)将流体域与固体域交界面处的每个单元分解为不同的次级单元；

(2)在分解得到的次级单元的外围分别划定1倍与1.5倍的边界；

(3)将(2)中划定得到的边界进行划分为9个小正方形。

上述的一种界面双向数据交换的流固耦合计算方法，step6中所述预接触搜索过程包括：

(1)将流体域与固体域交界面处的每个四边形单元分解为两个不同的三角形单元；

(2)在(1)分解得到的两个三角形单元的外围分别划定1倍与1.5倍的正方形边界得到两个正方形边界；

(3)将(2)中划定得到的两个正方形边界划分为9个小正方形从而将四边形单元分解为9个小单元。

上述的一种界面双向数据交换的流固耦合计算方法，step6中所述网格联系是在预接触搜索后，根据记录的节点信息搜索出与源网格单元节点邻近的目标网格单元节点，再通过节点信息找出源网格与目标网格相对应的最邻近单元。

上述的一种界面双向数据交换的流固耦合计算方法，step6中所述网格联系可以通过单元来寻找最邻近单元，步骤如下：

(1)将流体域与固体域交界面处的所有单元分解为三角形；

(2)建立边界方块，记录建立的边界方块与三角形目标单元相交产生的三角形小方块；

(3)将(2)记录的三角形小方块投射到平面上从而寻找到最邻近的三角形单元。

本发明的有益效果：

1、本方法充分考虑流固耦合的特点，采用多物理场耦合的cfd/csd双向流固耦合求解方法，将固体域与流体域视为整体的耦联体系结构，把流体域力场与固体域基础结构承载能力紧密联系起来，研究其结构损伤发展及变化的规律，从而为更深入地研究该耦联系统结构的流固耦合作用机理与损伤破坏模式提供支撑。

2、由于流体域在一些荷载下会产生波动且每时每刻都在变化，因此流体域的变化导致施加于固体域的力的属性也不一样，本发明根据流体域的变化通过mpcci数据交换，能够将实时的数据进行转化，得出的结果更真实可靠，流场细节展现更为丰富。

3、本方法计算精度较高、流场细节展现丰富，有效保证了耦合界面信息交换时能量守恒和边界层网格质量有效控制等问题，大大提高了模拟精度。

4、充分考虑流体域和固体域之间的相互影响，使流体力场与固体域基础结构承载能力紧密联系。

5、本发明使用的计算流程物理概念清晰，计算过程简单，易于程序设计，便于处理复杂的流固耦合问题。

附图说明

图1为本发明流固耦合计算流程图。

图2为湍流模型示意图；

图3为流固界面数据传递原理图；

图4为弹簧光顺和网格重构结合示意图。

图5为预接触搜索过程示意图。

图6为网格最小距离示意图。

图7为网格单元映射示意图。

具体实施方式

下面结合附图和实施例对本发明进一步说明。

实施例1：流体域与固体域相互接触时，固体域会受到应力和应变作用，这些力会使其结构产生变形，而结构的变形相应的也会对流体域产生影响，进而相互作用，针对该现象，本实施例提供了一种界面双向数据交换的流固耦合计算方法，该方法充分考虑流体域和固体域之间的相互影响，使流体力场与固体域基础结构承载能力紧密联系，研究其结构损伤发展及变化的规律。

在对流固界面进行计算之前需要研究清楚水体与渠道面板之间的相互作用力，这就需要根据不同的研究对象建立流固控制方程，流体域与固体域之间的联系是一种比较复杂的机制，建立流固控制方程以便于了解研究对象流体域与固体域之间的关系。流固控制方程包括流体控制方程和固体控制方程，其中流体控制方程主要包括质量守恒方程、动量守恒方程和能量守恒方程。

其中质量守恒表达式为：

式中：ρf表示为流体密度，t表示为时间，u表示为流体速度矢量；

动量守恒的表达式可根据牛顿第二定律推导得出：

式中：p为压力，fx、fy和fz为体力，u表示为流体速度矢量，τxx、τxy和τxz为表面粘性应力τ的分量，u、v和w分别表述为流体速度矢量u在x、y和z方向的分量；

流体能量守恒方程的表达式为：

式中：t表示为温度，ρf表示为流体密度，st表示为粘性耗散项，cp表示为比热容，k表示为热传导系数。

固体控制方程包括以下方程：

固体域动力学方程：

式中：ρs表示固体结构密度，表示固体结构加速度矢量，σs表示柯西应力张量，fs表示固体结构体积力矢量。

对于固体结构，由温差引起的热变形：

ft＝αt·▽t

式中：αt表示热膨胀系数。

在了解清楚研究对象流体域与固体域之间存在的作用力关系后，本实施例结合附图提供一种界面双向数据交换的流固耦合计算方法，如图1所示，该方法包括以下步骤：

step1、建立固体域模型：固体域计算结构动力学采用商用软件abaqus的隐式算法，建立csd固体域有限元模型，从而得到固体域的位移、应力和损伤等结果。

step2、对固体域模型进行有限元网格划分，再设置固体域单元类型、材料属性和结构损伤机制，进一步对固体域模型分别施加约束条件和施加载荷条件进行结构分析，使其与实际条件相符。

step3、根据研究对象不同的流体运动选择合适的湍流模型，不同湍流模型见图2，然后建立流体域模型，具体为：计算流体动力学采用fluent平台进行计算，形成cfd水体分体模型，从而计算得到流体域的流速、流量、水位、水压力等结果。

由图2可以看出，湍流模型主要分为涡黏模型、雷诺应力模型和大涡模拟模型三大类。

涡黏模型中可以根据建立湍流模型所需微分方程数目的不同，将涡黏模型分为三类，分别为零方程模型(代数模型)、单方程模型(apalart-allmaras模型等)和两方程模型(k-ω、k-ε模型等)。

零方程模型：通过代数关系式求解，往往只适用于简单流动。

单方程模型：在湍流时均连续方程和雷诺方程的基础上建立湍流能k的输运方程。

两方程模型：在湍动能k的基础上再引入ε耗散率方程，主要包括k-模型、剪切力输运模型和k-w模型。

(1)k-模型是计算流体动力学(cfd)中最常用的模型，它包含标准模型、rngk-ε模型、realizablk-ε模型，不同湍流模型的主要适用范围也有所差异，下面先对这三种湍流模型进行介绍。

标准k-ε模型：标准k-ε模型是一种半经验模型，其增加了关于湍动能k和消耗率ε两输运方程，适用于雷诺数大的区域，尤其适用于压力梯度较小的自由剪切层流，是目前工程应用最普遍的湍流模型，但该模型对于强旋流、弯曲边界层、无约束流等模拟上表现不佳。

rngk-ε模型：旋流在湍流中占有重要地位，标准k-ε模型在处理弯曲程度较大的流体中表现不佳，为了弥补这一缺陷，rngk-ε模型在此基础上通过对湍动黏度进行修正，从而可以更好地处理瞬变流与弯曲程度较大的流动状态。

realizablek-ε模型realizablek-ε模型也被称为带旋流修正k-ε模型，它是对标准k-ε模型和rngk-ε模型的补充，realizablek-ε模型的提出，是为了解决标准k-ε对于时均应变率较大的情形有可能导致负的正应力的问题，为符合合理的流体运动规律，需要对这种正应力进行约束。

标准k-ω模型包括标准k-ω模型和sstk-ω模型，其中标准k-ω是一种常见的两方程湍流模型，也被称为低雷诺数模型，是以wilcox模型为基础建立的。

sstk-ω模型：sstk-ω模型是通过将标准k-ω模型与k-ω模型结合一个混合函数相加所构建，该模型综合了标准k-ω模型在近壁低雷诺数区域计算精度高的优点与标准k-ω模型在湍流核心区计算精度高的优点，可以有效地模拟近壁区域的剪切效应，与k-ω模型在远场中保持独立。

雷诺应力模型(rsm)：在求解各分量的输运方程时，若采用雷诺应力模型，则方程数目较多，有较大的精确性，可以准确考虑各向异性效应，单对计算机内存要求较高，计算所需时间较长。

大涡模拟模型：对漩涡进行统计平均时需充分考虑湍流的统计特性，大涡模拟模型根据大涡模拟可以分为大尺度涡及小尺度涡，若想了解大涡的运动特性，则必须求解修正过的navier-stokes方程，采用第一、二类湍流模型计算小尺度涡。

step4、对流体域模型进行流体域网格划分，并对网格质量检查，网格质量合格时，设置边界条件和初始条件的参数；网格质量不合格时继续对网格质量进行检查，然后通过动网格技术进行网格重构并对流畅进行模拟，使其与实际条件相符。

step5、生成动网格：分别将step1中的csd固体域有限元模型和step2中的cfd水体分体模型通过mpcci联合起来进行计算生成动网格。本实施例在流固耦合过程中，为了结合弹簧光顺技术适应变形能力强、计算速度快和网格重构技术网格质量好的特点，如图4所示，将网格重构与弹簧光顺两种动网格方法联合使用来模拟流固耦合过程中的流体域边界变形与固体域变形协调一致。

若每一条边都假定为一个弹簧，且其具有一定刚度系数，则这样的网格系可以应用弹簧光顺方法。弹簧刚度系数如下述公式所示：

式中：rij表示相邻两节点间的距离。

网格运动后的网格节点位移可通过下式求得：

式中：ni表示相邻节点i的节点数，δrj表示节点j的运动位移。网格运动导致网格节点移动，位置产生了变化，具体位置如下式表示：

当弹簧光顺不能保障网格质量与精度时，使用网格重构法动网格技术，通过对局部以及整体动边界面附近网格进行重新定义划分。

在运用时可以根据实际情况选取不同的动网格技术，如网格重构法、弹簧光顺法或动态铺层法，其中弹簧光顺法适用于大多数的网格类型，主要是针对在流体边界处网格数量不变以及保障网格基本质量的情况下通过移动单元节点来改变单元的形状，但是网格的总体数量不会发生改变，这种方法能够较好的节约计算时间与得到精确的结果，但容易引起大部分网格结构产生畸变；网格重构法在流体分析中需要花费较多的时间与内存，同时在网格重新划分得时候重化网格本身的质量也不能够的得到保障；动态铺层法使用条件较为苛刻，只有在保证动边界相邻网格是楔形或六面体网格前提下才能使用，在计算过程中，当边界处网格向一边压缩或者拉伸的时候，动态铺层法通过对网格合并或者切割的方式来保障网格与边界移动的协调变形。

这三种动网格技术的优缺点如表1所示。

表1不同动网格技术的优缺点

step6、建立流固耦合的数据交换模型，通过流固耦合数据交换类型将流场模拟、动网格技术折合成固体域的施加荷载，并施加在固体域上；其中，数据交换是通过设置流固交界面来编辑水体与渠道面板之间的相互作用，动网格间数据交换是在数据传递过程中采用双向耦合的方法，流固耦合界面双向数据传递原理见图3，图3中固体域与流体域的网格单元一一对应，动网格间流固耦合过程分为三步：

1.预接触搜索：采用分块算法对流体域与固体域交界面处的单元进行分类，分别记录源网格和目标网格上的单元节点信息，为下一步网格联系做准备，具体过程包括：

针对二维单元，

(1)将流体域与固体域交界面处的每个四边形单元分解为两个不同的三角形单元；

(2)在(1)分解得到的两个三角形单元的外围分别划定1倍与1.5倍的正方形边界得到两个正方形边界；

(3)将(2)中划定得到的两个正方形边界划分为9个小正方形。

如图5所示，每一个三角形单元都能被一大一小两个边界正方形包围，小正方形可完全包含三角形，大边界包含小正方形。

同理，针对三维单元，每个三维单元均可被分解为四面体单元，包含四面体单元的空间可被分为9个小立方体，记录源网格与目标网格中相邻四面体单元的信息。

在搜索过程中，所有的二维单元均可被分解为三角形单元，所有的三维单元均可被分解为四面体单元，若所用单元类型达到统一，便都可采用同样的搜索算法。

2.网格联系：一般情况下，流体域与固体域交界面处的网格是不匹配的，源网格要想与目标网格进行数据传递，每一组源网格均需与对应的一组目标网格建立联系，为插值计算做准备；

预接触搜索后，在源网格上的单元节点信息与目标网格上的单元节点信息都被记录的前提下，具体步骤如下：

(1)根据记录的节点信息搜索出与源网格单元节点邻近的目标网格单元节点；

(2)通过节点信息找出源网格与目标网格相对应的最邻近单元；

(3)距离搜索中三角形与节点的相对位置及方法，先计算出p′点投影到各个三角单元的距离，通过三个距离可以找出最佳单元，首先是重心距离点p′的距离，由重心坐标u，v，w确定，如下式：

db＝|u|+|v|+|w|-1

当点p′的位置在三角形内时，那么db＝0；如若其位置落于三角形单元外，那么db的距离如图6所示；第二个距离值用dt表示，为离三角形单元最近的距离；第三个距离值表示点p到p′的距离dn。根据τ的取值来选取最佳单元，公式如下所示：

其屮，a表示三角形单元的面积。

3.通过插值函数进行插值计算：通过特定的插值算法进行耦合界面间的数据交换，插值算法需要同时满足位移连续性条件和动力连续性条件；

a位移连续性条件：

uf＝us

b动力连续性条件：

σsns＝σfnf

式中：uf、us分别表示流体和固体位移矢量；us分别表示结构与流体的柯西应力张量；σs、σf分别表述结构与流体的外法线方向矢量。

step7、通过以上各步骤将双向流固耦合视为整体的耦联体系结构，对该耦联体系结构进行分析得出耦连体系损伤的体制的结论，例如根据已知条件，得到固体域的损伤情况(位移、应力等)、流体域的流速、流量、水位、水压力等，为实际的建造、保护以及运行提供理论指导。

实施例2：本实施例与实施例1的相同之处不再赘述，不同之处在于：在step6中通过单元来寻找最邻近单元，具体步骤如下：

1、将所有单元分解为三角形；

2、建立边界方块，记录建立的边界方块与三角形目标单元相交产生的三角形小方块；

3、将2中记录的三角形小方块投射到平面上从而寻找到最邻近的三角形单元，如图7所示，为计算做准备。对于每个投射到平面的三角形，可分别计算出其相对所有交叉区域的面积与权重，其计算公式如下所示：

其中：ai表示交叉面积，表示交叉区域的总面积，权重大小不超过1。

以上所述仅为本发明的较佳实施例，并不限制本发明，凡在本发明的精神和原则范围内所做的任何修改、等同替换和改进，均应包含在本发明的保护范围之内。