2. 中国空气动力研究与发展中心空气动力学国家重点实验室, 四川绵阳 621000
2. State Key Laboratory of Aerodynamics, China Aerodynamics Research and Development Center, Mianyang 621000, China
在自然界中存在很多的流/固耦合现象, 如飞行器的机动飞行气动弹性问题复杂多体分离水流导致的桥墩振动台风导致的大跨度桥梁振动和变形鱼类的自主游动鸟类和昆虫的自主飞行等.这些问题的一个突出特点是流体和固体的非线性非定常耦合[1].对于气动弹性水流导致的桥墩振动台风导致的大跨度桥梁振动和变形等问题, 涉及流体动力学方程和结构动力学方程的耦合求解, 而对于飞行器的机动飞行复杂多体分离鱼类的自主游动鸟类和昆虫的自主飞行等问题则涉及到流体动力学方程和运动/动力学方程的耦合求解问题.当流场变化引起的流体动力变化的时间尺度与耦合运动的时间尺度大致接近时, 传统的线性叠加分析方法难以描述流体动力的非线性变化规律, 因此往往不能得到令人满意的结果[2], 此时人们往往采用实验手段[3]进行研究，若要进行数值模拟则需要采用非定常的耦合计算方法.本文重点关注流动/运动耦合问题的求解.
根据耦合策略的不同, 耦合计算方法一般可分为全耦合松耦合紧耦合3种[4-6].
全耦合即将不同子系统的控制方程联立并按统一的形式进行求解, 该方法在气动弹性领域已经得到了较为深入的研究和应用[7-10].但对于流动/运动耦合问题, 该方法比较难以实现, 因为流动控制方程与运动/动力学方程在性质上有较大的差异, 其求解方式有很大不同, 相关的研究工作也未见有报道.
目前, 对于流动/运动耦合问题大多采用松耦合方法, 该方法将流体动力学方程和运动/动力学方程解耦求解, 分别采用各自的时间积分, 然后交错推进时间获得系统的耦合解.该方法最大的优点是利于程序的模块化, 只需要对独立的流场求解模块和六自由度运动求解模块进行简单改造即可, 当某个模块需要替换为更先进的算法时, 对程序的修改量也较小[7].但该方法也有非常显著的缺点:由于积分时间的不同步, 两组方程之间存在一个时间步长的延时, 无法实现实时动态平衡, 需要选择较小的时间步长才能保证整个系统的稳定性和时间方向的计算精度.
2001年Söding [11]通过理论分析证明, 当流体运动产生的力严重依赖于物体的运动加速度时, 松耦合将导致计算不稳定, 尤其是在流体密度与运动物体的密度相当的情况, 如鱼体的自主游动问题.一些学者在随后的数值模拟中也验证了该理论分析, 他们均认为在研究鱼体的自主游动等问题时必须采用紧耦合的计算方法[11-13].对于飞行器在空中高速机动的问题, 虽然飞行器的质量密度较空气密度大很多, 但是流动的非定常效应有可能使得非定常气动力严重依赖于飞行器的运动加速度, 此时采用松耦合算法也有可能导致计算不稳定的现象.
为了解决松耦合方法存在的问题, 可以将流体动力学方程和运动/动力学方程都构造为含子迭代的隐式计算格式以进行耦合计算, 这种方法称为紧耦合方法, 其概念由Strganac等[14]在求解气动弹性问题时提出, 后被Alonso等[4]推广至流动/运动耦合问题, 并称之为全隐式方法.该方法在每个时间步内对流场控制方程和运动/动力学方程进行多次数据交换, 可以视为通过预估和校正迭代在一定程度上消除了时间步延时导致的误差, 从而提高了整个耦合系统的稳定性, 因此能够在不影响数值模拟精度的情况下取得比松耦合更长的时间步长.对于非线性耦合特性很强的流动/运动耦合问题, 紧耦合方法无疑能提高时间方向的数值模拟精度和计算的稳定性.但是, 由于涉及到隐式处理和内迭代, 每个真实时间步内的紧耦合计算的计算量无疑会大幅增加.因此, 针对不同的流动问题, 如何在稳定性计算精度和计算效率等方面需求平衡, 就是一个值得深入研究的问题.
为研究上述问题, 本文以课题组自主研发的CFD软件HyperFLOW [15-16]为基础, 建立了耦合动态混合网格生成非定常流场计算和飞行力学方程求解的一体化计算方法, 在统一框架内同时实现了松耦合与紧耦合方法, 并采用不同耦合方法, 对二维圆柱涡致自激振荡进行了数值研究.通过选取不同的运动参数, 探讨了松耦合与紧耦合方法的适应性, 并对其稳定性进行了理论分析.最后, 将耦合算法应用于二维鱼体自主游动和钝锥三自由度自由飞过程的数值模拟, 进一步证实了理论分析的结论.
1 一体化耦合计算方法流动/运动一体化计算需要耦合动态网格生成非定常流场计算飞行力学方程求解, 甚至还可能耦合飞行控制律.为此, 建立了如图 1所示的耦合计算框架, 以下就其中的关键环节进行简要介绍.
|
| 图 1 一体化耦合计算示意图 Fig.1 Procedure of coupling simulation |
因为混合网格较好地兼顾了几何适应性以及数值模拟精度, 代表了未来网格技术的发展方向, 所以本文工作基于混合网格开展.首先根据求解问题的特点, 对整体计算网格进行分区, 对于每个区定义不同的运动属性, 区分出远场固定网格随机体刚性或柔性运动的贴体网格随局部控制部件运动的贴体网格以及过渡区的变形网格等; 然后对不同的网格区域根据运动属性的不同选用不同的动网格策略, 如弹簧松弛法[17] Delaunay背景网格映射法[18]径向基函数方法(radial basis functions, RBF) [19]等.对于由于物体相对运动剧烈导致动网格的质量逐渐下降甚至出现交错的问题, 则采用局部网格重构的策略来弥补动网格技术的不足.关于动态混合网格生成的具体细节参见文献[20-21].
1.2 非定常N-S方程计算方法本文采用基于任意Lagrange-Euler (arbitrary Lagrange-Euler, ALE)体系的非定常计算方法.在ALE坐标系中, 非定常N-S方程可以写成如下的积分形式:
| $ \frac{\partial }{{\partial t}}\int\limits_{ V} \mathit{\boldsymbol{Q}} {\rm{d}}V + \oint\limits_S {({\mathit{\boldsymbol{F}}^{\rm{i}}}\left( \mathit{\boldsymbol{Q}} \right) - \mathit{\boldsymbol{Q}}{\mathit{\boldsymbol{v}}_{\rm{g}}}\cdot\mathit{\boldsymbol{n}}){\rm{d}}S} = \oint\limits_S {{\mathit{\boldsymbol{F}}^{\rm{v}}}} \left( \mathit{\boldsymbol{Q}} \right){\rm{d}}S. $ | (1) |
其中, S为控制体边界, V为控制体体积, Q为守恒变量, Fi, Fv分别为对流通量和黏性通量, vg为控制体边界运动速度, n为控制体边界外法向矢量.
本文采用格心型2阶精度有限体积算法求解公式(1).无黏项离散采用原始变量型的NND格式, 单元交界面上的通量计算采用近似Riemann解(如Roe格式[22]); 通过线性重构得到2阶精度的物理量分布, 而单元物理量的梯度由节点型的Gauss-Green公式求得.为了抑制激波附近的振荡保持单元内物理量分布的单调性, 一般采用Vankata Krishnan限制器[23]或Barth限制器[24].而在动网格上进行非定常计算需要满足几何守恒律, 通过推导得到了针对不同时间离散精度的网格交界面运动速度计算方法, 代入原离散方程即可自动满足几何守恒律[25].黏性项采用中心格式计算, 而湍流效应通过松耦合湍流模型方程的形式进行求解.可以选用SA一方程湍流模型[26]或SST两方程湍流模型[27], 湍流模型方程本身的计算与RAN-S方程的求解类似.对于不可压缩非定常流, 采用虚拟压缩法进行计算.时间推进采用了双时间步[28]和BLU-SGS隐式计算方法[29]; 为了提高计算效率, 发展了基于网格分区的大规模并行计算技术和多重网格加速收敛技术.更多具体细节参考文献[30-31].
1.3 松耦合与紧耦合计算方法如前所述, 对于非定常N-S方程和飞行力学方程，全耦合方法难以实现, 因此本文采用解耦方法进行计算.同时, 为了满足不同耦合程度的流动/运动耦合问题的计算需求, 本文建立了统一的框架, 可通过参数选取实现不同时间精度的松耦合和紧耦合计算.
非定常N-S方程动力学方程和运动学方程如式(2) 所示.因为非定常N-S方程的求解需要知道n+1时刻的网格, 所以需要先求解动力学方程和运动学方程.
| $ \left\{ \begin{array}{l} \frac{{\partial \mathit{\boldsymbol{Q}}}}{{\partial t}} = \mathit{\boldsymbol{f}}(\mathit{\boldsymbol{Q}},\mathit{\boldsymbol{U}},\mathit{\boldsymbol{X}})\\ \frac{{d\mathit{\boldsymbol{U}}}}{{dt}} = \mathit{\boldsymbol{g}}(\mathit{\boldsymbol{Q}},\mathit{\boldsymbol{U}},\mathit{\boldsymbol{X}})\\ \frac{{d\mathit{\boldsymbol{X}}}}{{dt}} = \mathit{\boldsymbol{h}}(\mathit{\boldsymbol{U}},\mathit{\boldsymbol{X}}) \end{array} \right.. $ | (2) |
其中，U为运动速度(角速度)矢量，X为质心的位置(姿态角)矢量.
若动力学/运动学方程采用显式推进, 在时间推进的每一步内, 动力学/运动学方程与N-S方程之间只进行一次数据交互, 即为前面提到的松耦合方式.以时间1阶精度的Euler显式格式为例, 松耦合时间推进交替迭代策略如式(3) 和图 2所示:
|
| 图 2 松耦合数据交互示意图 Fig.2 Sketch of data exchanges between CFD and RBD in loose coupling computation |
| $ \left\{ \begin{array}{l} {\mathit{\boldsymbol{U}}^{n + 1}} = {\mathit{\boldsymbol{U}}^n} + {\mathit{\boldsymbol{g}}^n}\Delta t\\ {\mathit{\boldsymbol{X}}^{n + 1}} = {\mathit{\boldsymbol{X}}^n} + {\mathit{\boldsymbol{h}}^n}\Delta t \end{array} \right.. $ | (3) |
RBD(rigid body dynamics)表示动力学和运动学方程的求解, CFD表示N-S方程的求解.
从图 2中可以看出, N-S方程和飞行力学控制方程可以各自做到高阶时间精度, 但由于积分时间的不同步, 使得两组方程之间始终存在一个时间步长的延时, 因此需要选择较小的时间步长才能保证整个系统的稳定性.所以, 对于耦合程度较高的问题, 需要发展紧耦合的计算方法.
若动力学/运动学方程采用双时间步隐式方法推进, 在物理时间推进的每一步内, 动力学/运动学方程与N-S方程之间进行多次数据交互, 即为前面提到的紧耦合方式.以时间1阶精度的Euler后差隐式离散格式为例, 采用双时间步方法, 紧耦合时间推进交替迭代策略如式(4) 和图 3所示, 而子迭代中各模块间的流程如图 4所示.
|
| 图 3 紧耦合数据交互示意图 Fig.3 Sketch of data exchanges between CFD and RBD in fully implicit coupling computation |
|
| 图 4 紧耦合计算中子迭代流程 Fig.4 Sub-iteration procedure of implicit coupling solver |
| $ \left\{ \begin{array}{l} \frac{{{\mathit{\boldsymbol{U}}^{m + 1}} - {\mathit{\boldsymbol{U}}^m}}}{{\Delta \tau }} + \Delta {\mathit{\boldsymbol{g}}^m} =- \frac{{ {\mathit{\boldsymbol{U}}^{m + 1}} - {\mathit{\boldsymbol{U}}^n}}}{{\Delta t}} + {\mathit{\boldsymbol{g}}^m}\\ \frac{{{\mathit{\boldsymbol{X}}^{m + 1}} - {\mathit{\boldsymbol{X}}^m}}}{{\Delta \tau }} + \Delta {\mathit{\boldsymbol{h}}^m} = - \frac{{{\mathit{\boldsymbol{X}}^{m + 1}} - {\mathit{\boldsymbol{X}}^n}}}{{\Delta t}} + {\mathit{\boldsymbol{h}}^m} \end{array} \right.. $ | (4) |
事实上, 物体运动的运动学和动力学控制方程为常微分方程组, 可以采用线性多步法求解, 由此可写成统一形式:
| $ \begin{array}{l} \mathit{\boldsymbol{Q}}_{{\rm{flight}}}^{^{m + 1}} = \left( {1 - {\boldsymbol{\omega}} } \right)\mathit{\boldsymbol{Q}}_{{\rm{flight}}}^{^m} + \omega [a\cdot\mathit{\boldsymbol{Q}}_{{\rm{flight}}}^{^n} + b\cdot\mathit{\boldsymbol{Q}}_{{\rm{flight}}}^{^{n - 1}} + \\ (c\cdot\mathit{\boldsymbol{F}}_{{\rm{flight}}}^{^m} + d\cdot\mathit{\boldsymbol{F}}_{{\rm{flight}}}^{^n} + e\cdot\mathit{\boldsymbol{F}}_{{\rm{flight}}}^{^{n - 1}})\Delta t]. \end{array} $ |
其中, Qflight为运动学参数U和X，Fflight为动力学控制方程的通量g和h，ω为松弛因子.根据上述统一形式, 选取不同的系数, 即可得到不同时间精度和耦合方式的计算方法，见表 1.这样即在统一的框架下实现了不同形式的松耦合和紧耦合算法, 而且非常容易实现编程.
|
下载CSV 表 1 统一形式飞行力学离散控制方程中各种格式的系数 Tab.1 Parameters in different coupling algorithms for the solution of flight dynamics equations |
对流动/运动耦合问题进行数值模拟需要将动态网格生成非定常流场计算六自由度运动方程求解等模块有机结合, 为实现各个功能模块协调统一的运作及相互调度, HyperFLOW软件采用了C++面向对象的程序设计思想, 将各个模块有机整合并实现了各模块之间的灵活调度.
2 一体化耦合计算方法的验证 2.1 物理模型本节通过两自由度圆柱涡致自激振荡算例对该一体化算法进行考核与分析, 其物理模型如图 5所示.由于该考核算例中仅涉及圆柱的刚体运动, 因此采用了刚性运动网格和弹簧松弛法等动网格技术进行耦合计算.
|
| 图 5 圆柱涡致自激振荡系统示意图 Fig.5 Sketch of 2DOFs VIV model of oscillating cylinder |
圆柱仅沿水平和垂直方向发生平动而没有转动, 且在这两个方向受到具有结构阻尼的线性弹簧约束.在初始时刻, 圆柱处于静止状态; 随后在水平来流作用下, 圆柱尾迹中会出现Karman涡街, 进而诱导出圆柱的自激振荡.其运动学和动力学控制方程可以简化为
| $ m\ddot x + c\dot x + kx = D,m{\ddot y} + c{\dot y} + ky = F. $ | (5) |
其中, m, D, F分别为圆柱单位展长的质量阻力和升力, 方程采用与N-S方程相同方式无量纲化后为
| $ \begin{array}{l} {{\ddot x}_{{\rm{ea}}}} + 2\zeta \left( {\frac{2}{{\bar u}}} \right){{\dot x}_{{\rm{ea}}}} + {\left( {\frac{2}{{\bar u}}} \right)^2}{x_{{\rm{ea}}}} = \frac{2}{{\mu _{\rm{s}}{\rm{\pi }}}}{C_{\rm{d}}},\\ {{\ddot y}_{{\rm{ea}}}} + 2\zeta \left( {\frac{2}{{\bar u}}} \right){{\dot y}_{{\rm{ea}}}} + {\left( {\frac{2}{{\bar u}}} \right)^2}{y_{{\rm{ea}}}} = \frac{2}{{\mu _{\rm{s}}{\rm{\pi }}}}{C_{\rm{l}}}. \end{array} $ | (6) |
式(6) 中各变量定义如下,
| $ \omega = \sqrt {\frac{k}{m}} ,\bar u = \frac{{{u_\infty }}}{{r\omega }},{\mu _s} = \frac{m}{{{\rm{\pi }}{\rho _\infty }{r^2}}},\zeta = \frac{c}{{2\sqrt {mk} }}. $ |
其中，r为圆柱半径.
来流条件取与文献[32-33]中相同的条件, 即M=0.1, Re=200, 且圆柱脱涡频率与弹簧阻尼系统固有频率相同, 此时整个系统发生一种频率耦合现象.为保证与文献计算状态相同, 先通过非定常计算获得定态圆柱脱涡频率St=0.1940, 此时的减缩速度ū应满足2/ū=2πSt, 于是得到ū=1.6408.物体与流体的密度比μs与文献中保持一致, 即μs=12.7324 (或40/π), 可以认为物体密度远远大于流体密度.
2.2 计算结果分析首先对阻尼比ζ=0.1的情况进行不同耦合方式的无量纲时间步长(Δt)收敛性研究, 表 2所示为不同Δt下各种耦合方式圆心y方向无量纲位移yea的最大最小值之差取常用对数(lg(2ymax)).在Δt较大时, 各种耦合方式得到的垂直方向位移峰值差有一定差异, 计算精度主要取决于N-S方程与RBD方程的时间精度; 若Δt过大, 某些松耦合方法甚至不能得到收敛解.为获得更可靠的结果, 应尽量选取CFD/RBD高阶时间离散的耦合方法, 本文后续计算也主要采用4种CFD/RBD均为时间2阶以上精度的耦合方法.随着Δt的减小, 4种方法逐渐趋于同一值(见图 6), 说明计算方法具有良好的时间步长收敛性.在这4种格式中, 有3种是紧耦合方法: RBD系统采用BDF2方法的时间精度最好; C-N方法次之; 而A-M3方法虽为3阶精度, 但因为受到N-S方程时间精度的限制, 整个系统只能达到2阶, 其精度反而低于前面两种2阶格式.相比之下, 采用松耦合策略的FDF2方法虽同为2阶格式, 但其精度与紧耦合方法相比较差.换言之, 要得到同样精度的解, BDF2紧耦合方法可以取较大的时间步长, 而FDF2松耦合方法却需要取较小的时间步长.从计算效率的角度来看, 虽然紧耦合方法需要采用隐式内迭代, 单个真实时间步的耗时较多; 但是由于其可以取较大的真实推进时间步长, 因此, 整体的计算效率可能并不比松耦合方法低.
|
|
下载CSV 表 2 不同Δt下各种耦合方式计算得到的lg(2ymax) Tab.2 lg(2ymax) computed by different coupling algorithms with different time steps |
|
| 图 6 4种CFD/RBD方法lg(2ymax)随Δt的变化 Fig.6 lg(2ymax) vs. time steps (computed by four types of second-order CFD/RBD coupling algorithms) |
图 7为该状态下典型x和y方向的位移变化规律曲线, 各种能够得到收敛解的方法获得的曲线都与该图类似, 曲线近似为正弦函数, 与文献[32-33]中描述相符.在物体与流体的密度比较大时, 不同耦合方式的收敛过程, 包括位移和相位基本一致, 紧耦合方法并未体现出明显的优势, 松耦合方法亦能得到较好的计算结果.这与Söding[11]的理论分析结论一致.
|
| 图 7 典型位移变化规律 Fig.7 Displacements in the x and y directions |
通过上述计算, 我们认为Δt取0.02已足以满足精度需求, 后续计算若未特殊说明均采用此值.将阻尼比ζ取0.001~0.2共15种状态, 分别采用前面提及的4种耦合方法对这些状态进行模拟.本文计算结果与参考文献[32-33]计算结果对比见图 8, 其中横坐标为减缩阻尼(2π2St2ζm/ρ∞r2)取常用对数, 各种耦合方法计算结果基本一致, 且与文献吻合较好, 进一步证明了方法的正确性.
|
| 图 8 不同耦合方式计算的lg(2ymax) vs. lg(2π2St2ζm/ρ∞r2) Fig.8 Displacements in the y-direction vs. the reduced damping coefficient (ζ) with different coupling algorithms |
在固定Reynolds数下, 当ζμs为常值时, 水平方向位移会随μs变化, 而垂直方向位移在进入稳定的周期运动后变化规律应基本相同, 经过计算验证确实如此.但当μs减小到1.273(即4/π)时, 松耦合计算的相位与紧耦合计算结果出现差异, 且水平方向位移出现不合理的振荡(如图 9); 当物体与流体密度相同时(μs=1), 松耦合方法计算发散, 减小时间步长和CFL数等条件也不能得到收敛的结果(如图 10); 而紧耦合方法即使在μs=0.1的情况下也能得到较好的结果.
|
| 图 9 μs=4/π时松耦合方法计算得到的位移变化曲线 Fig.9 Displacements in the x and y directions by loose coupling algorithm when μs=4/π |
|
| 图 10 μs=1时松耦合方法计算的位移变化曲线 Fig.10 Displacements in the x and y directions by loose coupling algorithm when μs=1 |
以下就松耦合计算发散而紧耦合方法稳定性更好的现象进行简化的理论分析.为便于分析, 将物体运动控制方程式(5) 进行简化, 仅考虑一维情形.与Söding [11]文献中所述方法类似, 将流体对加速运动物体的作用力近似视为附加质量力与附加阻尼力:
| $ \ddot x = - \frac{A}{m}\ddot x - \frac{{c + B}}{m}\dot x - \frac{k}{m}x + \frac{B}{m}{v_{{\rm{fluid}}}}. $ |
采用松耦合方法, 如FDF1对运动方程离散:
| $ \frac{{{{\dot x}^{n + 1}} - {{\dot x}^n}}}{{\Delta t}} = - \frac{A}{m}\frac{{{{\dot x}^n} - {{\dot x}^{n - 1}}}}{{\Delta t}} - \frac{{c + B}}{m}{{\dot x}^n} - \frac{k}{m}{x^n} + \frac{B}{m}v_{{\rm{fluid}}}^{^n}. $ |
不妨假设位移x和流体速度vfluid不存在误差, 而速度存在误差ε, 可得到误差关系方程:
| $ \frac{{{\varepsilon ^{n + 1}} - {\varepsilon ^n}}}{{\Delta t}} = - \frac{A}{m}\frac{{{\varepsilon ^n} - {\varepsilon ^{n - 1}}}}{{\Delta t}} - \frac{{c + B}}{m}{\varepsilon ^n}. $ | (7) |
再引入误差增长率G=εn+1/εn, 式(7) 化为
| $ \left( {G - 1} \right)\left( {G + \frac{A}{m}} \right) = - \frac{{\left( {c + B} \right)\Delta t}}{m}G. $ |
若算法是稳定的, 则要求|G| < 1, 于是得到稳定性条件:
| $ m >A + \frac{{ \left( {c + B} \right)\Delta t}}{2}. $ |
在本算例中, 二维圆柱的附加质量等于与物体所占体积相同的流体之质量, 换言之, 当采用该松耦合方法时, 物体密度必须大于流体密度，计算才可能稳定.
类似地, 若采用紧耦合方式, 如BDF1, 误差增长率为
| $ G = \frac{{1 + \frac{A}{m}}}{{1 + \frac{A}{m} + \frac{{\left( {c + B} \right)\Delta t}}{m}}}. $ |
可见, G < 1, 因此该紧耦合方法是恒稳的.
以上分析对于其他时间离散格式的耦合方法也可得到类似的结论, 并且对于更为普通的(不包括弹簧阻尼系统)问题也适用.这与Söding [11]及其他学者给出的结论一致:松耦合方法在物体与流体的密度接近或物体密度小于流体密度时会出现计算不稳定现象, 此时必须采用紧耦合方法.对于多维问题, 则需要考虑附加质量矩阵对稳定性的影响.
3 耦合计算方法的应用 3.1 二维鱼体简化模型S型启动过程数值模拟鱼类的自主游动是生物流体力学领域中典型的流动/运动耦合问题, 在启动过程中, 鱼体仅通过自身的摆动就可以在很短的时间内获得很大的加速度.本节利用前述的耦合计算方法, 针对二维情况下简化鱼体外形的S型启动过程进行了数值模拟.
鱼体外形近似为NACA0012翼型, 加速过程中的摆动规律如下:
| $ y(x,t) = A(x,t){\rm{sin}}\left[ {\frac{{2{\rm{\pi }}}}{\lambda }\left( {x - ct} \right)} \right]. $ |
其中，y为鱼体中轴线的法向位移, x为距离头部顶点的距离, λ为摆动波长, c为波速, 摆动振幅A为x和时间t的函数如下:
| $ A(x,t) = {\rm{min}}(1,\frac{t}{T})(0.02 - 0.0825x + 0.1625{x^2}). $ |
T为加速过程的时间, t=0~T的过程中振幅线性增加, t > T则进入正常巡游状态下的摆动.鱼体的运动及质量特性参数如表 3所示, 取T=2 s.由此可得, Reynolds数Re=cl/υ=10000(水的动力黏性系数υ近似为10-6 m2/s).
|
|
下载CSV 表 3 鱼体运动及质量参数 Tab.3 Parameters of 2D fish swimming |
分别对鱼体的单自由度运动及三自由度运动过程进行了数值模拟.其中单自由度运动中对法向以及俯仰方向均进行了限制, 鱼体只做水平方向的单自由度加速运动; 而三自由度情况下鱼体除了存在水平方向运动外, 还会在铅垂以及俯仰方向运动.计算过程中所采用的网格如图 11所示, 采用了四边形/三角形混合网格, 并在鱼体尾部进行了局部加密处理.从图中可以看出, 在整个启动过程中, 网格质量保持良好.具体的动网格生成方法见文献[31].
|
| 图 11 S型启动过程中的动态混合网格 Fig.11 Hybrid dynamic meshes during the S-type starting process |
在计算过程中, N-S方程采用2阶时间精度求解, 而RBD方程分别采用松耦合(FDF2) 和紧耦合(BDF2) 求解.正如2.3节中所述, 鱼体的自主游动问题中, 鱼体的密度与流体介质(水)的密度非常接近, 因此松耦合算法将导致计算发散, 必须采用紧耦合算法.计算过程也证实了这一理论分析结论.对于单自由度运动, 由于在铅垂和俯仰方向进行了限制, 相当于施加了外部作用力, 所以松耦合算法在时间步长极小(无量纲时间0.002) 的情况下还可以得到收敛解.但是利用紧耦合算法, 即使是在时间步长较大(无量纲时间0.1) 的情况下, 也能得到收敛解.而对于三自由度运动, 松耦合算法会很快发散, 只能采用紧耦合算法.因此, 后续计算中仅给出紧耦合的计算结果.
图 12为鱼体x方向质心位置及运动速度随时间的变化曲线.鱼的摆动幅度在2 s后达到最大值, 但水平方向的速度在大约5 s后才趋于稳定并进入巡游状态, 且单自由度运动的速度略大于三自由度运动.图 13为鱼体铅垂和俯仰方向的位移和(角)速度, 俯仰角呈现出的轻微不对称性是由前2 s不对称的启动动作导致的.
|
| 图 12 鱼体质心在水平方向的位移及速度变化情况 Fig.12 Translations and velocities in horizontal direction |
|
| 图 13 鱼体铅垂和俯仰方向的位移和(角)速度 Fig.13 Kinematics variables in the case of 3DOFs |
为了验证计算结果, 另外设计了一组算例:假设鱼体的质心固定, 在均匀来流中以不同频率(不同的摆动波速)摆动, 来流Reynolds数根据单自由度情况巡游速度确定为6600, 其他计算参数保持不变.前面提到的自主游动算例相当于在无来流的弹道靶中进行自由飞实验; 而这种固定质心算例则相当于在有来流的风洞中对不同摆动规律的鱼进行风洞实验[34-35].当阻力为0时, 鱼的摆动波速与巡游速度(来流速度)之比应与单自由度自主游动情况下的比值大致相同.图 14为固定质心时阻力随摆动波速的变化情况, 表 4给出了固定质心和自主游动的巡游速度与摆动波速之比.可以看出:固定质心得到的结果和单自由度自主游动基本一致; 而在三自由度情况下, 由于鱼体存在横向运动, 鱼体自身摆动输入能量一部分耗散至侧向, 因而导致摆动的效果减弱.图 15为单自由度和三自由度计算得到的流场涡量云图, 在一个摆动周期内, 有方向相反的一对旋涡脱落到尾迹中(系反向Karman涡街, 这是鱼体推力产生的主要原因), 且单自由度情况下尾涡强度稍大.
|
| 图 14 固定质心方式阻力随摆动波速的变化情况 Fig.14 Drag vs. wave speed when the center of gravity is fixed |
|
|
下载CSV 表 4 3种方法得到的鱼体巡游速度 Tab.4 Cruising speeds calculated by three different methods |
|
| 图 15 鱼体启动过程中若干典型位置流场的涡量云图 Fig.15 Vorticity contours during the S-type starting process |
本小节模拟了一个更为典型的气动/运动耦合问题——三维钝锥超声速自由飞, 来进一步考核前述的耦合方法.
钝锥尺寸外形和质量特性见图 16和表 5.飞行Mach数为5, 单位长度Reynolds数约为8.45×106, 其他计算参数见表 6.在实验中, 钝锥由气枪水平地逆风发射到超声速流中, 因为阻力的影响, 其水平方向速度渐渐降低, 最终反向至与流动方向相同.在气动力和重力的作用下, 钝锥运动轨迹大致呈抛物线, 并通过垂直于气流方向水平安装的相机将轨迹记录下来.
|
| 图 16 钝锥外形尺寸示意图 Fig.16 Configuration of blunted cone model |
|
|
下载CSV 表 5 钝锥质量特性 Tab.5 Mass characteristics of the blunt cone |
|
|
下载CSV 表 6 钝锥自由飞过程计算条件 Tab.6 Flow parameters of the blunt cone free-flight |
钝锥在刚脱离气枪时还带着一个用于保持气密性和发射稳定的夹块, 这个夹块会在阻力作用下与钝锥分离.为了消除夹块的影响, 根据0.3115 s时的状态来确定数值模拟时钝锥的初始六自由度参数, 即水平方向速度为-1.83 m/s，铅锤方向速度为-1.12 m/s，俯仰角为4.4°.因为实验没有提供壁面温度, 采用绝热壁条件进行模拟(事实上, 整个实验仅持续了大约0.5 s, 气动加热效应并不会对结果产生显著的影响).
与3.1节类似, 求解N-S方程采用2阶时间精度, 而求解RBD方程分别采用松耦合(FDF2) 和紧耦合(BDF2) 计算, 所采用的无量纲时间步长分别为dt=0.20, 0.10, 0.05.钝锥在相机平面内的三自由度运动情况和轨迹如图 17所示.可以看到, 不管是松耦合还是紧耦合, 计算结果都呈现出较好的时间步长收敛性, 不同耦合方式的计算结果基本相同, 并与实验基本吻合.这与2.3节理论分析的结论一致, 在物体密度与流体密度相比较大时, 松耦合和紧耦合都能较好地模拟问题, 计算精度仅取决于求解N-S和RBD控制方程所采用的时间离散格式.垂直方向位移和钝锥运动轨迹在后半段与实验略有不同, 这可能是初始状态的误差时间积累造成的.不同时间步长的俯仰角变化情况与其他变量相比有较大的差异, 这可能是多方面原因造成的:一是因为时间步长过大而不足以捕捉俯仰方向的高频振荡; 二是因为球锥底部网格密度不够而不足以捕捉流场中不停变化的涡系结构; 而且, 由于不在相机平面内的另外3个自由度难以通过相片辨识, 数值模拟中直接将这3个自由度的位移(角度)和初始(角)速度设为0, 这也会对俯仰力矩产生一定的影响.但根据现有的结果我们认为, 在合适的计算条件下, 能够得到随时间步长加密而收敛的俯仰角变化规律.图 18为典型时刻的实验纹影图和计算所得的压力云图, 可见, 计算与实验得到的激波位置基本吻合, 且流场结构合理.
|
| 图 17 位移(角度)变化情况及轨迹 Fig.17 3DOFs parameters vs. time and pathway of the cone |
|
| 图 18 典型状态时的纹影图和压力云图 Fig.18 Schlierens and pressure contours of typical states |
本文在HyperFLOW软件基础上建立了耦合动态混合网格生成非定常流场计算六自由度运动方程求解的一体化计算方法, 并在统一框架内同时实现了松耦合与紧耦合方法, 利用圆柱涡致自激振荡的模拟验证了该方法的适应性, 并对松耦合和紧耦合的稳定性进行了理论分析; 最后用于二维鱼体的自主游动和三维钝锥超声速自由飞过程的数值模拟.通过对上述算例的数值模拟和理论分析, 可以看出:
(1) N-S方程与RBD方程的时间精度都会对流动/运动耦合问题的计算精度产生影响, 当时间步长过大时某些松耦合方法甚至不能得到收敛解.在相同时间步长情况下, 紧耦合方法在时间精度方面较松耦合有一定优势.
(2) 对于本文研究的问题, 当物体密度与流体密度相比较大时, 不同耦合方式的计算结果基本一致, 紧耦合并未体现出明显的优势, 且由于紧耦合每一子迭代步都需要重新获得计算网格信息, 计算效率反而会降低.
(3) 当物体密度与流体密度较为接近或小于流体密度时, 流体运动产生的力严重依赖于物体的运动加速度, 此时需要采用紧耦合方法; 如采用松耦合方法, 可能会出现计算不稳定现象.
下一步我们将对上述现象开展更为深入的研究和分析, 开展多自由度耦合计算稳定性分析, 并将进一步耦合飞行控制律计算结构力学动力系统等计算方法和计算模块, 最终实现复杂外形飞行器的数值虚拟飞行模拟.
致谢 本文工作得到国家自然科学基金(11532016，11202227) 的资助，钝锥自由飞实验数据由中国空气动力研究与发展中心高超声速研究所马晓宇钟俊郭雷涛朱涛等提供, 在此表示衷心的感谢.| [1] |
Dean J P, Clifton J D, Bodkin D J, et al.High resolution CFD simulations of maneuvering aircraft using the CREATE/AV-Kestrel solver[R].AIAA 2011-1109, 2011.
|
| [2] |
Allan M R, Badcock K J, Richards B E.CFD based simulation of longitudinal flight mechanics with control[R].AIAA 2005-0046, 2005.
|
| [3] |
赵忠良, 吴军强, 李浩, 等. 高机动导弹气动/运动/控制耦合的风洞虚拟飞行试验技术[J]. 空气动力学学报, 2016, 34(1): 14-19. Zhao Z L, Wu J Q, Li H, et al. Wind tunnel based virtual flight testing of aerodyanmics, flight dynamics and flight control for high maneuver missle[J]. ACTA Aerodynamica Sinica, 2016, 34(1): 14-19. (in Chinese) |
| [4] |
Alonso J J, Jameson A.Fully-implicit time-marching aeroelastic solutions[R].AIAA 1994-0056, 1994.
|
| [5] |
Kamakoti R, Shyy W. Fluid-structure interaction for aeroelastic applications[J]. Progress in Aerospace Sciences, 2004, 40(8): 535-558. DOI:10.1016/j.paerosci.2005.01.001 |
| [6] |
陈龙, 伍贻兆, 夏健. CFD/CSD紧耦合及新型动网格方法在气动弹性模拟中的应用[J]. 计算力学学报, 2011, 28(5): 807-812. Chen L, Wu Y Z, Xia J. CFD/CSD closely coupled and new dynamic grid method in application of aeroelastic simulation[J]. Chinese Journal of Computational mechanics, 2011, 28(5): 807-812. DOI:10.7511/jslx201105027 (in Chinese) |
| [7] |
Farhat C, Van Der Zee K G, Geuzaine P. Probably second-order time-accurate loosely-coupled solution algorithms for transient nonlinear computational aeroelasticity[J]. Computer Methods in Applied Mechanics and Engineering, 2006, 195(17/18): 1973-2001. |
| [8] |
Etienne S, Pelletier D, Garon A.A monolithic formulation for steady-state fluid-structure interaction problems[R].AIAA 2004-239, 2004.
|
| [9] |
Etienne S, Pelletier D, Garon A.A monolithic formulation for unsteady fluid-structure interactions[R].AIAA 2006-694, 2006.
|
| [10] |
Tan B H, Lucey A D, Howell R M. Aero-/hydro-elastic stability of flexible panels:prediction and control using localized spring support[J]. Journal of Sound and Vibration, 2013, 332(26): 7033-7054. DOI:10.1016/j.jsv.2013.08.012 |
| [11] |
S9ding H.How to integrate free motions of solids in fluids[A].//Fourth Numerical Towing Tank Symposium[C].Hamburg, 2001.
|
| [12] |
Leroyer A, Visonneau M. Numerical methods for RANS simulations of a self-propelled fish-like body[J]. Journal of Fluids and Structures, 2005, 20(7): 975-991. DOI:10.1016/j.jfluidstructs.2005.05.007 |
| [13] |
Liu G, Yu Y L, Tong B G. Flow control by means of a traveling curvature wave in fishlike escape responses[J]. Physical Review E:Statistical Nonlinear and Soft Matter Physics, 2011, 84(2): 2269-2278. |
| [14] |
Strganac T W, Mook D T. Numerical model of unsteady subsonic aeroelastic behavior[J]. AIAA Journal, 1990, 28(5): 903-909. DOI:10.2514/3.25137 |
| [15] |
He X, He X Y, He L, et al. HyperF LOW:a structured/unstructured hybrid integrated computational environment for multi-purpose fluid simulation[J]. Procedia Engineering, 2015, 126: 645-649. DOI:10.1016/j.proeng.2015.11.254 |
| [16] |
赫新, 赵钟, 马戎, 等. HyperF LOW亚跨声速流动验证[J]. 空气动力学学报, 2016, 34(2): 267-275. He X, Zhao Z, Ma R, et al. Validation of HyperF LOW in subsonic and transonic flow[J]. ACTA Aerodynamica Sinica, 2016, 34(2): 267-275. (in Chinese) |
| [17] |
Batina J T. Unsteady Euler airfoil solutions using unstructured dynamic meshes[J]. AIAA Journal, 1990, 28(8): 1381-1388. DOI:10.2514/3.25229 |
| [18] |
Liu X Q, Qin N, Xia H. Fast dynamic grid deformation based on Delaunay graph mapping[J]. Journal of Computational Physics, 2006, 211(2): 405-423. DOI:10.1016/j.jcp.2005.05.025 |
| [19] |
Rendall T C, Allen C B. Efficient mesh motion using radial basis functions with data reduction algorithms[J]. Journal of Computational Physics, 2009, 228(17): 6231-6249. DOI:10.1016/j.jcp.2009.05.013 |
| [20] |
张来平, 邓小刚, 张涵信. 动网格生成技术及非定常计算方法进展综述[J]. 力学进展, 2010, 40(4): 424-447. Zhang L P, Deng X G, Zhang H X. Reviews of moving grid generation techniques and numerical methods for unsteady flow[J]. Advances in Mechanics, 2010, 40(4): 424-447. DOI:10.6052/1000-0992-2010-4-J2009-123 (in Chinese) |
| [21] |
张来平, 赫新, 常兴华, 等. 复杂外形静动态混合网格生成技术研究新进展[J]. 气体物理, 2016, 1(1): 42-61. Zhang L P, He X, Chang X H, et al. Recent progress of static and dynamic hybrid grid generation techniques over complex geometries[J]. Physics of Gases, 2016, 1(1): 42-61. (in Chinese) |
| [22] |
Roe P L. Approximate Riemann solvers, parameter vectors, and difference schemes[J]. Journal of Computational Physics, 1981, 43(2): 357-372. DOI:10.1016/0021-9991(81)90128-5 |
| [23] |
Venkatakrishnan V.On the accuracy of limiters and convergence to steady state solutions[R].AIAA 1993-0880, 1993.
|
| [24] |
Barth T J, Jespersen D C.The design and application of upwind schemes on unstructured meshes[R].AIAA1989-0366, 1989.
|
| [25] |
Chang X H, Ma R, Zhang L P, et al. Further study on the geometric conservation law for finite volume method on dynamic unstructured mesh[J]. Computers & Fluids, 2015, 120: 98-110. |
| [26] |
Spalart P R, Allmaras S R.A one equation turbulence model for aerodynamic flows[R].AIAA 1992-0439, 1992.
|
| [27] |
Menter F R.Zonal two equation turbulence models for aerodynamic flows[R].AIAA 1993-2906, 1993.
|
| [28] |
Jameson A.Time dependent calculations using multigrid, with applications to unsteady flows past airfoils and wings[R].AIAA 1991-1596, 1991.
|
| [29] |
Chen R F, Wang Z J. Fast block lower-upper symmetric Gauss-Seidel scheme for arbitrary grids[J]. AIAA Journal, 2000, 38(12): 2238-2245. DOI:10.2514/2.914 |
| [30] |
Zhang L P, Wang Z J. A block LU-SGS implicit dual time-stepping algorithm for hybrid dynamic meshes[J]. Computers & Fluids, 2004, 33(7): 891-916. |
| [31] |
Zhang L P, Chang X H, Duan X P, et al. Applications of dynamic hybrid grid method for three-dimensional moving/deforming boundary problems[J]. Computers & Fluids, 2012, 62: 45-63. |
| [32] |
Blackburn H M, Karniadakis G E.Two and three-dimensional simulations of vortex-induced vibration of a circular cylinder[A].//Proceedings of the International Society of Offshore and Polar Engineers(ISOPE)93 Conference[C].Singapore, 1993, 3:715-720.
|
| [33] |
Morton S A, Melville R B, Visbal M R.Accuracy and coupling issues of aeroelastic Navier-Stokes solutions on deforming meshes[R].AIAA 1997-1085, 1997.
|
| [34] |
Zhang L P, Chang X H, Duan X P, et al. A block LU-SGS implicit unsteady incompressible flow solver on hybrid dynamic grids for 2D external bio-fluid simulations[J]. Computers & fluids, 2009, 38(2): 290-308. |
| [35] |
Chang X H, Zhang L P, He X. Numerical study of the thunniform mode of fish swimming with different Reynolds number and caudal fin shape[J]. Computers & fluids, 2012, 68: 54-70. |