CN102436514B

CN102436514B - 一种流固耦合网格更新的方法

Info

Publication number: CN102436514B
Application number: CN 201110247524
Authority: CN
Inventors: 卢旦; 杨联萍; 李承铭
Original assignee: Shanghai Xiandai Archtectural Design (group) Co Ltd
Current assignee: East China Architectural Design and Research Institute Co Ltd
Priority date: 2011-08-25
Filing date: 2011-08-25
Publication date: 2013-05-08
Anticipated expiration: 2031-08-25
Also published as: CN102436514A

Abstract

本发明提出一种流固耦合网格更新的方法，应用于土木工程领域，该方法包括下列步骤：分别记录流体网格和固体网格两套网格的节点与单元信息；以每个固体网格三角形单元的重心坐标系为参考坐标系，计算出每个流体网格中的节点在该坐标系中的局部坐标；根据判定条件确定流体网格节点与固体网格单元间的从属关系；根据实际情况使固体网格的单元节点发生改变；计算流体网格节点的新坐标，并将所述流体网格节点更新。本发明提出的流固耦合网格更新的方法，利用三角形重心坐标进行网点的定位，与传统方法相比具有计算量小，判别条件少，计算过程简洁等特点。

Description

一种流固耦合网格更新的方法

技术领域

本发明涉及土木工程行业的结构分析领域，且特别涉及一种流固耦合网格更新的方法。

背景技术

流固耦合现象是工程中常见的一种自然现象，如土壤与基础的耦合、建筑物振动与风荷载的耦合。虽然，目前的固体软件和流体软件均能在各自的领域做出精确的分析，但当考虑两者的相互耦合计算时却碰到了许多困难（见文献【1】：沈世钊，武岳，大跨度张拉结构风致动力响应研究进展，第十届全国风工程学术会议论文集，2001年，P：10~16）。其中之一便是：网格数量多的一方（流体网格）需要跟踪网格数量少的一方（固体网格）进行变形。以往的方法，如样条插值法、守恒法等，见文献【2】（王彬，流固耦合作用的弱耦合算法及风与薄膜结构的耦合分析[D]，北京交通大学博士学位论文，2008年7月）、文献【3】（孙晓颖，薄膜结构风振响应中的流固耦合效应研究[D]，哈尔滨工业大学博士学位论文，2007年1月），虽然保证了数据传递的准确性，但由于计算方法过于繁琐，导致计算速度成倍下降，使得这种先进的技术无法应用于实际工程，只能停留在小规模的科学研究（如文献【3】采用基于薄板样条法的插值计算程序，由于“计算量巨大，缺乏工程可操作性”，只好采用“简化的数值模拟方法”）。

如上所述，目前大多数流固耦合方法在实际工程尤其是特大型工程的应用中，当进行变形插值计算时需要耗费漫长的机时，因而是无法满足实际工程需要的。为此，需要提出一种既满足计算精度要求，又能极大程度提高计算速度的变形跟踪计算方法。

本发明提出的这种方法既能保证计算数据的准确，同时最明显的特征是：流程清楚、计算公式简单，因此计算速度得到极大提高，使之在实际工程中的应用成为可能。

发明内容

本发明提出一种流固耦合网格更新的方法，应用于土木工程行业的结构分析领域，在土木工程的流固耦合计算分析中，流体网格需要根据固体网格的变动而做出相应的变化，本方法用于快速、准确地实现这一过程。

为了达到上述目的，本发明提出一种流固耦合网格更新的方法，应用于土木工程领域，该方法包括下列步骤：

分别记录流体网格和固体网格两套网格的节点与单元信息；

以每个固体网格三角形单元的重心坐标系为参考坐标系，计算出每个流体网格中的节点在该坐标系中的局部坐标；

根据判定条件确定流体网格节点与固体网格单元间的从属关系；

根据实际情况使固体网格的单元节点发生改变；

计算流体网格节点的新坐标，并将所述流体网格节点更新。

进一步的，所述记录两套网格的节点与单元信息步骤包括：每个固体网格节点与单元的从属关系以及它们的坐标；每个流体网格的节点坐标。

进一步的，所述根据判定条件确定流体网格节点与固体网格单元间的从属关系步骤包括：如果某一流体网格节点对于某一固体网格单元的三角形满足判别条件，则定义其从属于上述固体网格单元。

进一步的，所述计算流体网格节点的新坐标步骤采用以下公式计算：

u = \frac{(v_{1} \cdot v_{1}) (v_{2} \cdot v_{0}) - (v_{1} \cdot v_{0}) (v_{2} \cdot v_{1})}{(v_{0} \cdot v_{0}) (v_{1} \cdot v_{1}) - (v_{0} \cdot v_{1}) (v_{1} \cdot v_{0})};

v = \frac{(v_{0} \cdot v_{0}) (v_{2} \cdot v_{1}) - (v_{0} \cdot v_{1}) (v_{2} \cdot v_{0})}{(v_{0} \cdot v_{0}) (v_{1} \cdot v_{1}) - (v_{0} \cdot v_{1}) (v_{1} \cdot v_{0})},

其中，v₀＝C-A,v₁＝B-A,v₂＝P-A，三角形ABC代表固体网格中的某一三角形单元，A、B、C分别为从属于该单元的节点，点P代表落在三角形ABC内部的流体网格的一个节点。

进一步的，如果系数u或v为负值，那么相当于朝相反方向移动，即BA或CA，因此，如果：u≥0,v≥0，且u+v≤1则可以判断点P位于三角形ABC内部。

本发明提出一种流固耦合网格更新的方法，利用三角形重心坐标进行网点的定位，与传统方法相比具有计算量小，判别条件少，计算过程简洁等特点，因而其直接的效果是计算速度非常快。本方法巧妙地采用三角形重心坐标原理，使得原先位于一个稀疏网格（固体网格）三角形单元内部的稠密网格（流体网格）节点会随着稀疏网格（固体网格）节点的变形依然贴附于原三角形内。三角形内的点会随着三角形在平面内、外发生大变形，但这些点对所从属的整个三角形的拓扑关系不发生变化，三角形重心原理在本领域内的使用尚属首次，且真是由于该原理算法朴素、计算简单，很容易使本领域的普通技术人员能够掌握和实施。

附图说明

图1所示为本发明较佳实施例的流固耦合网格更新的方法流程图。

图2所示为重心坐标计算原理示意图。

图3所示为变形跟踪系统实际使用效果对比图。

图4所示为流固耦合计算流程图。

具体实施方式

下面结合附图对本发明的实施例作详细说明：本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

请参考图1，图1所示为本发明较佳实施例的流固耦合网格更新的方法流程图。本发明提出一种流固耦合网格更新的方法，应用于土木工程领域，该方法包括下列步骤：

步骤S100：分别记录流体网格和固体网格两套网格的节点与单元信息；

步骤S200：以每个固体网格三角形单元的重心坐标系为参考坐标系，计算出每个流体网格中的节点在该坐标系中的局部坐标；

步骤S300：根据判定条件确定流体网格节点与固体网格单元间的从属关系；

步骤S400：根据实际情况使固体网格的单元节点发生改变；

步骤S500：计算流体网格节点的新坐标，并将所述流体网格节点更新。

以上五步即为本方法的主要过程与思路。通常，第三步仅在计算初始化时使用，而第五步计算中仅涉及少量的四则运算，因此整套计算的效率非常高。

根据本发明较佳实施例，所述记录两套网格的节点与单元信息步骤S100包括：每个固体网格节点与单元的从属关系以及它们的坐标；每个流体网格的节点坐标。所述根据判定条件确定流体网格节点与固体网格单元间的从属关系步骤S300包括：如果某一流体网格节点对于某一固体网格单元的三角形满足判别条件，则定义其从属于上述固体网格单元。

进一步的，请参考图2，图2所示为重心坐标计算原理示意图。图2中三角形ABC代表稀疏网格（固体网格）中的某一三角形单元，A、B、C分别为从属于该单元的节点，点P代表落在三角形ABC内部的稠密网格（流体网格）的一个节点。假设以点A作为起点，那么点B相当于在AB方向移动一段距离得到，而点C相当于在AC方向移动一段距离得到。则对于平面内任意一点，可由如下方程来表示：

P＝A+u×(C-A)+v×(B-A) （1）

整理方程（1）得到

P-A＝u×(C-A)+v×(B-A) （2）

令v₀＝C-A,v₁＝B-A,v₂＝P-A则，

v₂＝u×v₀+v×v₁ （3）

等式（3）两边分别点乘v₀和v₁得到两个等式：

(v₂)·v₀＝(u×v₀+v×v₁)·v₀ （4）

(v₂)·v₁＝(u×v₀+v×v₁)·v1 （5）

求解上述方程（4）、（5）得到：

v = \frac{(v_{1} \cdot v_{1}) (v_{2} \cdot v_{0}) - (v_{1} \cdot v_{0}) (v_{2} \cdot v_{1})}{(v_{0} \cdot v_{0}) (v_{1} \cdot v_{1}) - (v_{0} \cdot v_{1}) (v_{1} \cdot v_{0})} - - - (6)

v = \frac{(v_{0} \cdot v_{0}) (v_{2} \cdot v_{1}) - (v_{0} \cdot v_{1}) (v_{2} \cdot v_{0})}{(v_{0} \cdot v_{0}) (v_{1} \cdot v_{1}) - (v_{0} \cdot v_{1}) (v_{1} \cdot v_{0})} - - - (7)

所述计算流体网格节点的新坐标步骤S500根据公式（6）、（7）计算：

u = \frac{(v_{1} \cdot v_{1}) (v_{2} \cdot v_{0}) - (v_{1} \cdot v_{0}) (v_{2} \cdot v_{1})}{(v_{0} \cdot v_{0}) (v_{1} \cdot v_{1}) - (v_{0} \cdot v_{1}) (v_{1} \cdot v_{0})};

v = \frac{(v_{0} \cdot v_{0}) (v_{2} \cdot v_{1}) - (v_{0} \cdot v_{1}) (v_{2} \cdot v_{0})}{(v_{0} \cdot v_{0}) (v_{1} \cdot v_{1}) - (v_{0} \cdot v_{1}) (v_{1} \cdot v_{0})},

再请参考图3，图3所示为变形跟踪系统实际使用效果对比图。可见本方法巧妙地采用三角形重心坐标原理，使得原先位于一个稀疏网格（固体网格）三角形单元内部的稠密网格（流体网格）节点会随着稀疏网格（固体网格）节点的变形依然贴附于原三角形内。三角形内的点会随着三角形在平面内、外发生大变形，但这些点对所从属的整个三角形的拓扑关系不发生变化。

流固耦合计算中有以下两个量需要传递：

（1）流体域需要将荷载传递给固体域；

（2）固体域需要将变形传递给流体域。

这两组信息传递的思路分别是：

（1）对每一个有限元单元，搜索与其单元中心距离最短的流体力学单元，将该流体力学单元上的荷载传递给有限元单元。

（2）对每一个流体力学膜面节点，首先判断其原始坐标位置落入哪个有限元单元范围内；其次，根据对应有限元单元的变形来修改该流体力学膜面节点的变形，即所采用“三角形重心坐标”法。

整个流固耦合计算系统的流程如图5所示：

具体计算过程由以下步骤和程序组成：

1、分别建立流体和固体模型。

2、在有限元中用“creat-femct.inp”程序获取结构单元的形心坐标文件“femct.txt”。

3、在流体力学中用“creat-cfdct.c”程序获取流体单元的形心坐标文件“mo_center.txt”和“mos_center.txt”。

4、在有限元中用“creat-area-index.inp”程序生成由流体单元向固体单元传递荷载所需的对应文件“mo-index.txt”和“mos-index.txt”。

5、在流体力学中用“creat-grid.c”程序生成流体力学中循环膜面节点的节点坐标文件“mo-grid.txt”和“mos-grid.txt”。

6、在有限元中用“creat-grid-index.inp”程序生成膜面节点对应文件“mo-grid-index.txt”和“mos-grid-index”。

7、同时运行流体和固体力学计算程序，进行流固耦合计算。

以下为流固耦合计算所需各程序及相应说明，本例所用有限元软件为ansys，流体力学软件为fluent。

一、creat-femct.inp

二、creat-cfdct.c

三、creat-area-index.inp

四、creat-grid.c

五、creat-grid-index.inp

综上所述，本发明提出一种流固耦合网格更新的方法，利用三角形重心坐标进行网点的定位，与传统方法相比具有计算量小，判别条件少，计算过程简洁等特点，因而其直接的效果是计算速度非常快。本方法巧妙地采用三角形重心坐标原理，使得原先位于一个稀疏网格（固体网格）三角形单元内部的稠密网格（流体网格）节点会随着稀疏网格（固体网格）节点的变形依然贴附于原三角形内。三角形内的点会随着三角形在平面内、外发生大变形，但这些点对所从属的整个三角形的拓扑关系不发生变化，三角形重心原理在本领域内的使用尚属首次，且真是由于该原理算法朴素、计算简单，很容易使本领域的普通技术人员能够掌握和实施。2011年的上海世博轴索膜结构，结构（固体）模型计算单元数量数十万，流体模型计算单元上千万，如果采用过去的方法，需要几个月时间才能完成流固耦合计算，而采用本文方法，可在十几天内完成计算，从而保证了重要工程的如期进行。

虽然本发明已以较佳实施例揭露如上，然其并非用以限定本发明。本发明所属技术领域中具有通常知识者，在不脱离本发明的精神和范围内，当可作各种的更动与润饰。因此，本发明的保护范围当视权利要求书所界定者为准。

Claims

1.一种流固耦合网格更新的方法，应用于土木工程领域，其特征在于，该方法包括下列步骤：

分别记录流体网格和固体网格两套网格的节点与单元信息；

根据实际情况使固体网格的单元节点发生改变；

计算流体网格节点的新坐标，并将所述流体网格节点更新，所述计算流体网格节点的新坐标步骤采用以下公式计算：

u = \frac{(v_{1} \cdot v_{1}) (v_{2} \cdot v_{0}) - (v_{1} \cdot v_{0}) (v_{2} \cdot v_{1})}{(v_{0} \cdot v_{0}) (v_{1} \cdot v_{1}) - (v_{0} \cdot v_{1}) (v_{1} \cdot v_{0})};

v = \frac{(v_{0} \cdot v_{0}) (v_{2} \cdot v_{1}) - (v_{0} \cdot v_{1}) (v_{2} \cdot v_{0})}{(v_{0} \cdot v_{0}) (v_{1} \cdot v_{1}) - (v_{0} \cdot v_{1}) (v_{1} \cdot v_{0})},

其中，v₀＝C-A,v₁＝B-A,v₂＝P-A，三角形ABC代表固体网格中的某一三角形单元，A、B、C分别为从属于该单元的节点，点P代表落在三角形ABC内部的流体网格的一个节点；如果系数u或v为负值，那么相当于朝相反方向移动，即BA或CA，因此，如果：u≥0,v≥0，且u+v≤1则可以判断点P位于三角形ABC内部。

2.根据权利要求1所述的流固耦合网格更新的方法，其特征在于，所述记录两套网格的节点与单元信息步骤包括：每个固体网格节点与单元的从属关系以及它们的坐标；每个流体网格的节点坐标。

3.根据权利要求1所述的流固耦合网格更新的方法，其特征在于，所述根据判定条件确定流体网格节点与固体网格单元间的从属关系步骤包括：如果某一流体网格节点对于某一固体网格单元的三角形满足判定条件，则定义其从属于上述固体网格单元。