流固耦合算法研究报告

流固耦合算法研究报告
1 流固耦合的概念
流固耦合力学是流体力学与固体力学交叉而生成的一门力学分支，它是研究变形固体在流场作用下的各种行为以及固体位形对流场影响这二者相互作用的一门科学。

流固耦合力学的重要特征是两相介质之间的相互作用，变形固体在流体载荷作用下会产生变形或运动。

变形或运动又反过来影响流，从而改变流体载荷的分布和大小，正是这种相互作用将在不同条件下产生形形色色的流固耦合现象。

流固耦合问题可由其耦合方程定义，这组方程的定义域同时有流体域与固体域。

而未知变量含有描述流体现象的变量和含有描述固体现象的变量，一般而言具有以下两点特征：
1）流体域与固体域均不可单独地求解
2）无法显式地削去描述流体运动的独立变量及描述固体现象的独立变量从总体上来看，流固耦合问题按其耦合机理可分为两大类：
第一类问题的特征是耦合作用仅仅发生在两相交界面上，在方程上的耦合是由两相耦合面上的平衡及协调来引入的如气动弹性、水动弹性等。

第二类问题的特征是两域部分或全部重叠在一起，难以明显地分开，使描述物理现象的方程，特别是本构方程需要针对具体的物理现象来建立，其耦合效应通过描述问题的微分方程来体现。

实际上流固耦合问题是场（流场与固体变形场）间的相互作用：场间不相互重叠与渗透其耦合作用通过界面力（包括多相流的相间作用力等）起作用，若场间相互重叠与渗透其耦合作用通过建立不同与单相介质的本构方程等微分方程来实现。

流固耦合的数值计算问题，早期是从航空领域的气动弹性问题开始的，这也就是通过界面耦合的情况，只要满足耦合界面力平衡，界面相容就可以。

气动弹性开始主要是考虑机翼的颤振边界问题，计算采用简化的气动方程和结构动力学方程，从理论推导入手，建立耦合方程，这种方法求解相对容易，适应性也较窄。

现在由于数值计算方法，计算机技术的发展，整个的求解趋向于N S
方程与非线性结构动力学。

一般使用迭代求解，也就是在流场，结构上分别求解，在
各个时间步之间耦合迭代，收敛后再向前推进。

好处就是各自领域内成熟的代码稍作修改就可以应用。

其中可能还要涉及一个动网格的问题，由于结构的变形，使得流场的计算域发生变化，要考虑流场网格随时间变形以适应耦合界面的变形。

流固耦合求解时有两种方式：两场交叉迭代和直接全部同时求解。

也就是我们平时说的单向耦合和双向耦合。

单向流固耦合是在不同的耦合场进行交叉迭代，通过场间耦合媒介交换耦合信息，这种方式应用的比较多，但是操作起来比较麻烦，而且由于不是同步求解，所以不能十分精确的表达流固之间的相互影响；双向耦合是把流场与固体场耦合到同一个有限元方程中进行计算，这种方法反映的情况较单向耦合准确，缺点是数值处理的难度较大，在ANSYS 11.0中可以很好的解决这个问题。

故在本文中采用的是双向耦合的求解方法。

在CFX11.0和ANSYS Workbench11.0中进行双向流固耦合计算的计算流程如图1所示：
图1 双向流固耦合计算流程 2 网格变形理论
在双向流固耦合分析中，ANSYS 中进行的结构分析的结果作为载荷传递到CFX 中进行求解。

同样的，CFX 中的分析结果又返回到ANSYS 当中。

例如，流体的压力可以被当作边界的载荷加载到结构计算中，从结构计算中得到的变
是
是
形、速度、加速度也可以作为载荷返回到CFX 的流体计算中。

这样的过程将会一直持续到两种计算都达到平衡为止。

在双向流固耦合过程中要实现固体的变形，其先决条件是周围的网格发生移动、变形。

我们可以把每一个网格看成是一个微小的控制体，以控制体的观点来分析网格的变形。

这种方式是利用网格把整个空间域离散成有限的控制体。

控制方程对每一个控制体积分，使其有关的量以不连续的形式存在。

图2.1是典型的网格单元。

很明显，每一个节点都被一组控制体的表面包围着，所有的求解量和流体特性都集中在节点上。

在直角坐标系下，质量守恒、动量守恒和守恒标量方程表述如下： ()0j j U t x ρρ∂
∂+=∂∂ (2-38) ()()j i i j i eff j i j
j i U U P U U U t x x x x x ρρμ⎡⎤⎛⎫∂∂∂∂∂∂+=-++⎢⎥ ⎪ ⎪∂∂∂∂∂∂⎢⎥⎝⎭⎣⎦ (2-39) ()()j eff j j j U S t x x x φφρφρφ⎡⎤⎛⎫∂∂∂∂+=Γ+⎢⎥ ⎪ ⎪∂∂∂∂⎢⎥⎝⎭⎣⎦
(2-40) 将上述方程在控制体上积分，并应用高斯定理将体积分转换成面积分。

当控制体没有及时变形的情况下，时间导数可以挪到积分外，得到：
图2 网格单元
0j j V s
d dV U dn dt ρρ+=⎰⎰ (2-41) i
j i i j i j j eff j U j i V s s s V U U d U dV U U dn Pdn dn S dV dt x x ρρμ⎛⎫∂∂+=-+++ ⎪ ⎪∂∂⎝⎭
⎰⎰⎰⎰⎰ (2-42)
j j eff j j V s s V d dV U dn dn S dV dt x φ
φρφρφ⎛⎫∂+=Γ+ ⎪ ⎪∂⎝⎭
⎰⎰⎰⎰ (2-43) 这里V 和S 分别表示积分的体和面，j dn 是正交面外法线方向矢量的微分。

体积分代表控制体内原有的项，面积积分代表的是总的流量。

当控制体发生及时的变形时，积分形式的守恒方程将要做适当变换，遵从莱布尼兹法则得到下列变换：
()j j V t V s
d dV dV W dn dt t φφφ∂=+∂⎰⎰⎰ (2-44) 其中，j W 是控制体边界的速度
和前面一样，将微分形式的守恒方程在指定的控制体上积分。

此时莱布尼兹法则适用，积分形式的守恒方程变换为：
()()0j j j V t s
d dV U W dn dt ρρ+-=⎰⎰ (2-45) ()()i j i i j j i j j eff j U j
i V t s s s V U U d U dV U W U dn Pdn dn S dV dt x x ρρμ⎛⎫∂∂+-=-+++ ⎪ ⎪∂∂⎝⎭⎰⎰⎰⎰⎰(2-46) ()()j j j eff j j V t s s V d dV U W dn dn S dV dt x φ
φρφρφ⎛⎫∂+-=Γ+ ⎪ ⎪∂⎝⎭
⎰⎰⎰⎰ (2-47) 衰减项说明变形的控制体的内部储存发生了变化，对流项说明网格越过了控制体的移动边界。

3 ANSYS 多物理场求解器中流固耦合的原理
ANSYS 多物理场求解器可以将耦合载荷在不同的网格中自动传递。

有两种插值方法可以应用到载荷的传递，它们分别是截面守恒法和完全守恒法。

在截面守恒插值时，接受方的每一个节点映射到发送方()i α的一个单元上。

传递变量内差值替换i α，传递值()i i T ϕα=。

因而，接受方的所有节点询问发送方。

下图反映了截面守恒插值的过程：
图3 截面守恒插值
在完全守恒插值时，发送方的每一个节点X 映射到接受方的一个单元上。

因而发送方的传递变量分裂成两份加到接受方上。

如下图所示，节点4上的力分
裂到节点3'、4'上：
图4 完全守恒插值
如下图所示，在完全守恒插值时，最好是在发送方有比较精细的网格，而在接收方有较粗糙的网格。

当发送方的网格精细时，接收方可以很充分的捕捉到其传递过来的应力。

反之如果接收方的网格比发送方的网格精细则不容易捕获。

完全守恒插值时发送方为精细网格
图5 完全守恒插值时接收方为精细网格
上面两点适用与发送方或者接收方是由高阶单元组成的情况。

如果想要制造高阶单元和低阶单元指尖的点对点映射，则需要谨慎。

例如下图所示，一个完全守恒的载荷通过有相同单元数的双方交界面传递的情况。

当接收方是高阶时，是不能产生正确的截面的。

图6 3个低阶单元
如果想得到正确的截面，需要把发送方的低阶单元数加倍。

如下图所示。

同样需要注意的是不能落下接触面的中部节点。

图7 6个低阶单元
对于力、热流、和热生成，可以指定完全守恒或者截面守恒的插值方法；对于变形和温度传递问题只能用截面守恒。

为了在不同的网格交界面指尖传递载荷，一个网格的节点必须映射到相应的另一个网格的单元上。

多物理场求解的求解法则必须在每个交界面的两个映射中履行。

例如流固耦合问题中为了传递变形量，流体节点必须映射到固体单元上。

同样的道理，在传递应力时，固体节点也必须映射到流体单元上。

图8 流固耦合载荷传递。