一种基于物质点法的流固耦合高效插值计算方法与流程

1.本发明涉及力学领域，更具体的说是涉及一种基于物质点法的流 固耦合高效插值计算方法。


背景技术：

2.求解流固耦合问题时，数据在流体网格和结构网格的插值传递至 关重要，因为流体和固体所用的离散方法一般不同(流体用有限体积 法，固体用有限元法)，加上流体一般要考虑边界层，在变量梯度较 大的区域需要局部加密，导致流体网格比固体网格要密得多，在网格 插值后网格重构的过程，如果变形过大会出现单元面积、体积等出现 负数和不光滑、失真的现象，特便是大变形(断裂、穿透等)很容易 出现，计算中如果出现这种现象，计算程序就终止，一方面影响流固 耦合研究进展，另一方面影响流固耦合的工程化应用。因此，研究一 种流固耦合计算过程的高效插值方法意义重大。


技术实现要素：

3.有鉴于此，本发明提供了一种基于物质点法的流固耦合高效插值 计算方法，建立了基于物质点法的流固耦合数据交换机制，成功将物 质点法应用于流固耦合的分析中，其特点是不用划分网格，在精度和 效率不损失的前提下，可以适应大变形的计算中。
4.为了实现上述目的，本发明采用如下技术方案：
5.一种基于物质点法的流固耦合高效插值计算方法，包括以下步 骤：
6.导入网格；
7.进行流体求解；
8.将流体求解结果插值到网格上，获得结构载荷；
9.求解结构力学方程，获得结构有限元网格节点位移以及位移插 值；
10.将结构壁面网格节点位移插值至流体网格边界，得到流体网格边 界节点的位移，重构网格后再进行流场计算；
11.输出计算结果。
12.其中，重构网格后需重新进行流体方程的求解，输出的计算结果 包括：流场边界处的速度矢量、压力及温度等。
13.可选的，进行流体求解前需要进行流体计算设置，包括设置物性 参数、计算模型、离散格式的选择、控制参数。
14.可选的，还包括获得每个网格的位移后，重构网格，得到新的网 格。
15.其中，网格重构的目的在于获取高质量网格，避免流体网格和固 体网格反复插值出现的单元面积、体积等出现负数和不光滑、失真的 现象，以及求解过程中的数值发散。
16.可选的，所述结构力学方程
[0017][0018]
其中，m为质量矩阵，c为阻尼矩阵，k为刚度矩阵，q为插之 后有限元网格上的气动
力，x为节点位移。
[0019]
可选的，将节点位移插值至网格，单元插值模式包括线单元、三 角形单元和四边形单元。
[0020]
可选的，所述线单元的形函数如下：
[0021]
n1(ξ1,ξ2)＝ξ1＝x/l；
[0022]
n2(ξ1,ξ2)＝ξ2＝1-ξ1；
[0023]
其中，ni(ξ1,ξ2)表示参数坐标表示的形函数，ξ1,ξ2表示参数坐标。
[0024]
可选的，所述三角单元的形函数如下
[0025][0026][0027][0028]
其中，a
t
＝a1+a2+a3，ξ1,ξ2表示参数坐标；其中，n1(ξ1,ξ2)、 n2(ξ1,ξ2)、n3(ξ1,ξ2)表示三角形单元的3个形函数；
[0029]at
表示三角形单元的面积；
[0030][0031][0032][0033]
a1、b1、c1分别为待定常数，(x,y)为节点坐标；
[0034]
其中，各形函数可以理解为对应节点发生单位位移而其它节点不 动时，单元内的位移分布情况。
[0035]
可选的，所述四边形单元的形函数如下：
[0036]
n1(ξ1,ξ2)＝(1-ξ1)(1-ξ2)；
[0037]
n2(ξ1,ξ2)＝ξ1(1-ξ2)；
[0038]
n3(ξ1,ξ2)＝ξ1ξ2；
[0039]
n4(ξ1,ξ2)＝(1-ξ1)ξ2；
[0040]
ξ1,ξ2表示参数坐标。
[0041]
其中，n1(ξ1,ξ2)、n2(ξ1,ξ2)、n3(ξ1,ξ2)、n4(ξ1,ξ2)表示四边形单元的 4个形函数；ξ1，ξ2表示参数坐标。
[0042]
各形函数可以理解为对应节点发生单位位移而其它节点不动时， 单元内的位移分布情况。
[0043]
经由上述的技术方案可知，与现有技术相比，本发明公开提供了 一种基于物质点法的流固耦合高效插值计算方法，无需划分网格，并 且在精度不损失的前提下计算效率比
常规的流固耦合计算要高，并且 能够避免常规计算方式中流体网格和固体网格反复插值出现的单元 面积、体积等出现负数和不光滑、失真的现象，因此可以有效解决工 程和科研中的流固耦合计算过程的数据传递问题，可以适应大变形的 计算中。
附图说明
[0044]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面 将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而 易见地，下面描述中的附图仅仅是本发明的实施例，对于本领域普通 技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附 图获得其他的附图。
[0045]
图1为本发明流固耦合的计算方法的详细步骤；
[0046]
图2为本发明网格数据/物理数据插值流程示意图；
[0047]
图3为本发明插值单边单元基函数参数示意图；
[0048]
图4为本发明插值三角形单元基函数参数示意图；
[0049]
图5为本发明插值四边形单元基函数参数示意图。
具体实施方式
[0050]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方 案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部 分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普 通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例， 都属于本发明保护的范围。
[0051]
本发明实施例公开了一种基于物质点法的流固耦合高效插值计 算方法，具体实施方式如下：
[0052]
如图1所示，首先导入网格数据，求解流体方程，获得流体网格 壁面压力信息，将此压力信息插值到结构网格壁面上，求解结构方程、 获得结构网格上的位移、将此位移插值到流体网格上，进行网格重构， 计算下一时间点或输出计算结果，返回求解流体方程直到计算终止。 计算步骤如下：
[0053]
步骤1导入网格数据，此步骤由第三方软件(如： ansys-icem,gridgen)完成并存成plot3d格式，结构模型网格由 ansa软件计算生成。
[0054]
步骤2进行流体计算设置，包括设置物性参数、计算模型、离散 格式的选择、控制参数等，此过程由公司自行开发的adi.simwork4.0 中的计算流体力学模块的参数设置界面上设置完成。
[0055]
步骤3进行流体求解，包括求解过程中的残差曲线观测，此过 程由公司自行开发的adi.simwork4.0中的计算流体力学模块完成。
[0056]
步骤4气动力插值，此过程将步骤3计算出的壁面气动力从流体 壁面插值到结构表面网格上，获得结构载荷。此过程自编程序完成。
[0057]
步骤5求解结构力学方程，获得结构有限元网格节点位移，此过 程编程完成，结构力学方程
[0058][0059]
其中(1)式中m为质量矩阵，c为阻尼矩阵，k为刚度矩阵， q为插之后有限元网格上
的气动力，x为节点位移。
[0060]
步骤6获得节点位移，位移插值，将步骤5求解出来的节点位 移再插值回流体网格上去，得到流体网格每个节点的位移。三种单元 插值的示意图可以如图3-图5所示，其中ξ1，ξ2为当地坐标。此 过程由自编程序完成。
[0061]
步骤7网格变形(重构)，完成步骤6后，获得了每个网格气动 单元的位移，然后进行气动网格重构，并得到新的气动网格。此过程 由自编程序完成。
[0062]
步骤8输出计算结果或返回步骤3进行下一时间点计算，进行 后处理，此过程由公司自行开发的simwork4.0后处理功能或第三方 软件完成，如tecplot等。
[0063]
如图2-5所示，步骤6中插值时由于气动网格和结构网格一般由 不同的前处理器生成，且网格和数据格式类型很不一致，由于本插值 程序不兼容其他数据格式，需要将插值前数据转化为ddf格式，插 值后再转化为其他格式，数据接口由自编程完成，具体的过程如下：
[0064]
ⅰ
)将目标网格数据及源网格数据、物理数据等转化为ddf格 式。
[0065]
ⅱ
)由目标网格数据和源网格数据生成插值转化矩阵。
[0066]
ⅲ
)插值程序进行插值。
[0067]
ⅳ
)生成目标网格ddf格式数据。
[0068]
ⅴ
)将目标网格上的ddf格式数据进行转换
[0069]
ⅵ
)得到目标网格上的数据。
[0070]
插值过程简述如下：
[0071]
插值时从源网格上找出离每个目标网格节点距离最近的单元，然 后找出目标节点到该单元的最短距离点(投影点)，用该投影点的物 理量来表示目标点上的物理量,投影点上的物理量用所在单元的形函 数求得(如要将气动网格节点上的气动荷载传到结构节点上，结构节 点便是目标节点，气动网格为源网格)。由于网格只是对曲面几何的 近似逼近，结构网格与气动网格疏密程度一般不同，所以目标节点一 般不在某个源网格单元内。然后用形函数计算投影点上的物理量，将 其近似作为目标点上的物理量。
[0072]
对于源网格单元，可以用参数坐标和节点物理坐标表示其它点的 物理坐标，其表达式如(2)式
[0073][0074]
其中，(2)中的表示单元内任意点j的物理坐标， 表示单元节点物理坐标，i表示单元节点编号，ni(ξ1,ξ2)表示参数 坐标表示的形函数，ξ1,ξ2表示参数坐标。对于结构几何表面插值 时线单元、三角形单元和四边形单元足以满足要求，其三者的形函数 分别为：
[0075]
线单元的形函数
[0076]
n1(ξ1,ξ2)＝ξ1＝x/l
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0077]
n2(ξ1,ξ2)＝ξ2＝1-ξ1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0078]
三角形单元的形函数
[0079][0080][0081][0082]
其中：a
t
＝a1+a2+a3[0083]
四边形单元的形函数
[0084]
n1(ξ1,ξ2)＝(1-ξ1)(1-ξ2)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0085]
n2(ξ1,ξ2)＝ξ1(1-ξ2)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)
[0086]
n3(ξ1,ξ2)＝ξ1ξ2ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(10)
[0087]
n4(ξ1,ξ2)＝(1-ξ1)ξ2ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(11)
[0088]
目标点的投影坐标应该满足式(12)，其中目标点的物理坐标， 为已知量。
[0089][0090]
由极值点出的导数为0可得投影点满足的等式约束方程，如式 (13)
[0091][0092]
对于线单元和三角形单元方程(13)为一个线性方程组，直接用求 解线性方程组的方法可得到投影点的参数坐标。对于四边形单元，方 程为双线性方程组，需要用牛顿-拉夫逊方法求解。对于三角形单元 将式(5-7)代入方程(13)可得线性方程组如下方程(14)。
[0093][0094]
一旦求得了目标点在源单元中的投影点，便可用形函数求得投影 点的相应物理量，如式(15)
[0095][0096]
其中为投影点上的物理量，可将其作为目标点上的物理 量；为源单元节点上的物理量。比如：可为结构气动单元节点 上的位移量。将方程(16)写成矩阵形式，如式(16)
[0097]
[fa]＝[t
as
][fs]
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(16)
[0098]
其中，
[0099]
对于流固耦合问题，主要涉及到将结构表面节点上的位移插值到 气动表面节点上，然后将气动表面节点上的气动荷载插值到结构表面 节点上，整个过程如(17-20)式所示。方程(17)表示结构变形后的气动 力计算，(18)为气动网格节点上的气动力向结构网格传递，(19)为通 过结构方程求解结构位移，(20)表示将结构网格节点上的位移传递到 气动网格节点上。
[0100][0101][0102][0103][0104]
为了保证数据传递过程的精度，首先得保证传递过程中功守恒。 由虚功原理可得：
[0105][0106]
将(18)代入(21)式得：
[0107][0108]
比对方程(22)和方程(18)可得：
[0109]
[t
sa
]＝[t
as
]
t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(23)
[0110]
方程(23)是一个重要的结论，它表明从结构网格到气动网格的传 递矩阵与从气动网格到结构网格的传递矩阵互为转置矩阵。
[0111]
上述插值格式也能保证插值前气动网格上的总合力等于插值后 结构网格上的总合力，如方程(24)：
[0112][0113]
由于气动网格和结构网格不能完全切合，所以力矩传递时会有相 应的精度损失。
下面对此做简要说明：
[0114][0115]
方程(25)中的为气动网格节点在结构网格上的投影点到矩心 的距离，而非气动网格节点到矩心的距离，两者的差别是由于气动表 面与结构网格表面不完全切合造成的，如式(26)。
[0116][0117]
注：如果结构太简化，否则该方法插值的效果会不光滑。在做模 型时尽量保证气动外形和结构外形一致，有必要时在结构模型上做一 些虚拟单元(不增加结构的刚度和质量，只是保证气动外形和结构外 形一致，方便插值)来弥补其不足。
[0118]
上面全部说明了一种流固耦合计算的全部过程以及在计算中新 的插值技术为了验证这种插值的有效性，我们选取了f-16战斗机作 为算例，其几何模型使用公开的数据。结构模型由机身、机翼、挂架 和挂载四部分组成，机身、挂架和挂载都采用四节点壳单元，材料属 性见表1；采用集中质量元来模拟挂架和挂载的质量特性(其属性如 表2)；机翼分为三部分，根部、挂架部分和尖部，分别采用四节点 壳单元，材料属性如表3。气动表面网格采用结构化的多块网格，总 数为344块，网格节点数为106万。飞行条件为ma数0.85、攻角2.12 度、侧滑角0度、来流速度为279.1m/s、高度为2436m。
[0119]
s1：利用商业软件分别生成结构外形网格和气动外形网格，此过 程由ansa和ansys-icem软件完成。
[0120]
s2进行气动力计算，进行插值计算出结构模态，此过程由公司 研制的simwork4.0完成气动力计算，进行插值，导入nastran计 算结构模态。
[0121]
s3将结构模态(取前八阶模态)插值到气动网格上。
[0122]
s4输出结果文件，利用simwork4.0进行后处理。
[0123]
通过计算验证，并与实验结果进行对标，证实了我们建立的流固 耦合分析的一般方法和新的插值方法的正确性和可执行性。
[0124]
对所公开的实施例的上述说明，使本领域专业技术人员能够实现 或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来 说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的 精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被 限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新 颖特点相一致的最宽的范围。