搜索

文章查询

x

留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

基于直角网格法的单个和阵列布置下柔性水翼绕流数值模拟

辛建建 陈振雷 石凡 石伏龙

基于直角网格法的单个和阵列布置下柔性水翼绕流数值模拟

辛建建, 陈振雷, 石凡, 石伏龙
PDF
HTML
导出引用
导出核心图
  • 研究柔性水翼在不可压缩流体中的水动力特性, 对于船舵和减摇鳍等海洋结构物的设计和性能优化具有重要意义. 本文将自主开发的径向基函数虚拟网格法求解器扩展到模拟绕单个或多个柔性水翼的不可压缩流动问题. 数值模型基于虚拟网格有限差分法考虑浸入边界对流场的影响, 引入紧支径向基函数(compact supported radial basis function, CSRBF)以物面Lagrangian质点追踪复杂的柔性动边界. 基于该方法, 首先模拟了均匀流中主动拍动的柔性水翼, 升阻力系数良好的网格收敛性结果验证了本文方法的精度和可靠性. 并研究了柔性水翼在不同振荡频率下的水动力特性, 阐述了柔性水翼的推力生成机制. 然后模拟了绕阵列布置柔性水翼的流动现象, 研究了不同间距和不同振荡频率下水翼表面的升阻力系数变化规律和尾涡特性, 观察到紧密布置的柔性水翼在高频振荡下推力系数存在显著的放大效应, 同时推力为零时的临界频率提前.
      通信作者: 陈振雷, chenzhenlei@nbu.edu.cn
    • 基金项目: 国家级-国家青年科学基金项目(51909124)
    [1]

    Zabihi M, Lari K, Amiri H 2017 J. Mech. Sci. Technol. 31 3539

    [2]

    Yang J M 2016 J. Hydrodyn. Ser. B 28 713

    [3]

    王力, 田方宝 2018 中国科学: 物理学 力学 天文学 48 094703

    Wang L, Tian F B 2018 Sci. China Phys. Mech. Astron. 48 094703

    [4]

    Al-Marouf M, Samtaney R 2017 J. Comput. Phys. 337 339

    [5]

    Huang W X, Chang C B, Sung H J 2011 J. Comput. Phys. 230 5061

    [6]

    Tullio M D D, Pascazio G 2016 J. Comput. Phys. 325 201

    [7]

    吴晓笛, 刘华坪, 陈浮 2017 物理学报 22 224702

    Wu X D, Liu H P, Chen F 2017 Acta Phys. Sin. 22 224702

    [8]

    Yeo K S, Ang S J, Shu C 2010 Comput. Fluids 39 403

    [9]

    Tian F B, Wang W, Wu J, Sui Y 2016 Comput. Fluids 124 1

    [10]

    Bergmann M, Iollo A, Mittal R 2014 Bioinspiration Biomimetics 9 046001

    [11]

    Khalid M S U, Akhtar I, Dong H 2016 J. Fluids Struct. 66 19

    [12]

    Khalid M S U, Akhtar I, Imtiaz H, Dong H, Wu B 2018 Ocean Eng. 157 108

    [13]

    Mittal R, Dong H, Bozkurttas M, Najjar F M, Vargas A, Loebbecke A V 2008 J. Comput. Phys. 227 4825

    [14]

    王亮, 吴锤结 2011 力学学报 43 18

    Wang L, Wu C J 2011 Chin. J. Theor. Appl. Mech. 43 18

    [15]

    王文全 2014 计算力学学报 31 646

    Wang W Q 2014 Chin. J. Comput. Mech. 31 646

    [16]

    辛建建, 石伏龙, 金秋 2017 物理学报 66 186

    Xin J J, Shi F L, Jin Q 2017 Acta Phys. Sin. 66 186

    [17]

    Xin J J, Shi F L, Jin Q, Lin C 2018 Comput. Fluids 176 210

    [18]

    Leer B V 1979 J. Comput. Phys. 32 101

    [19]

    Van d V H A, Kohn R V, Iserles A, Ciarlet P G, Wright M H 2006 Iterative Krylov Methods for Large Linear Systems (Cambridge: Cambridge University Press) pp133–135

    [20]

    Liu G R, Gu Y T 2005 An Introduction to Meshfree Methods and Their Programming (Netherlands: Springer) p73

    [21]

    Sui Y, Chew Y T, Roy P, Low H T 2007 Int. J. Numer. Methods Fluids 53 1727

  • 图 1  虚拟网格法的二维插值格式 (a)无虚拟网格; (b)一个虚拟网格; (c)两个虚拟网格

    Fig. 1.  Two-dimensional interpolation stencil of the ghost cell method: (a) No ghost cells; (b) one ghost cell; (c) two ghost cells.

    图 2  静止或运动边界流动求解算法的流程图

    Fig. 2.  Flowchart of numerical algorithm for stationary or moving boundary flows.

    图 3  五个频率下的力系数比较 (a) 平均阻力系数; (b) 均方根升力系数

    Fig. 3.  Comparison of force coefficients under five oscillation frequencies: (a) Average drag coefficient; (b) root mean square lift coefficient.

    图 4  三个频率下升阻力系数随无因次时间的变化趋势 (a) 阻力系数; (b) 升力系数

    Fig. 4.  Variation of drag and lift coefficients with the dimensionless time under three oscillation frequencies: (a) Drag coefficient; (b) lift coefficient.

    图 5  $\omega $ = 0.15时在两个时刻的局部瞬时涡等值线 (a) $\theta = $$ {\text{π}}/2$; (b) $\theta = 3{\text{π}}/2$

    Fig. 5.  Local instantaneous vortex contours for $\omega $ = 0.15 at two moments: (a) $\theta = {\text{π}}/2$; (b) $\theta = 3{\text{π}}/2$.

    图 7  $\omega $ = 0.75时在两个时刻的局部瞬时涡等值线 (a) $\theta = $${\text{π}}/2$; (b) $\theta = 3{\text{π}}/2$

    Fig. 7.  Local instantaneous vortex contours for $\omega $= 0.75 at two moments: (a) $\theta = {\text{π}}/2$; (b) $\theta = 3{\text{π}}/2$.

    图 6  $\omega $ = 0.45时在两个时刻的局部瞬时涡等值线 (a) $\theta = $$ {\text{π}}/2$; (b) $\theta = 3{\text{π}}/2$

    Fig. 6.  Local instantaneous vortex contours for$\omega $ = 0.45 at two moments: (a) $\theta = {\text{π}}/2$; (b) $\theta = 3{\text{π}}/2$.

    图 8  阵列布置水翼游动的几何模型

    Fig. 8.  Geometrical model of the undulating hydrofoil in array arrangement.

    图 9  三个间距下群游中间位置水翼和单水翼在五个频率下的力系数比较 (a)平均阻力系数; (b)均方根升力系数

    Fig. 9.  Comparison of force coefficients between the central hydrofoil in array arrangement for three gap distances and the single hydrofoil under five frequencies: (a) Average drag coefficient; (b) root mean square lift coefficient.

    图 10  群游中间位置水翼在三个频率下升阻力系数随无因次时间的变化趋势 (a)阻力系数; (b)升力系数

    Fig. 10.  Variation of drag and lift coefficients on the central hydrofoil in array arrangement with the dimensionless time under three oscillation frequencies: (a) Drag coefficient; (b) lift coefficient.

    图 11  $\omega $ = 0.15时在两时刻的水翼群游局部涡场图  (a) $\theta = {\text{π}}/2$; (b) $\theta = 3{\text{π}}/2$

    Fig. 11.  Local instantaneous vortex contours of undulating hydrofoils in array arrangement for $\omega $ = 0.15 at two moments: (a) $\theta = {\text{π}}/2$; (b) $\theta = 3{\text{π}}/2$.

    图 12  $\omega $ = 0.45时在两时刻的水翼群游局部涡场图 (a) $\theta = $$ {\text{π}}/2$; (b) $\theta = 3{\text{π}}/2$

    Fig. 12.  Local instantaneous vortex contours of undulating hydrofoils in array arrangement for $\omega $ = 0.45 at two moments: (a) $\theta = {\text{π}}/2$; (b) $\theta = 3{\text{π}}/2$.

    图 13  $\omega $ = 0.75时在两时刻的水翼群游局部涡场图 (a) $\theta = $$ {\text{π}}/2$; (b) $\theta = 3{\text{π}}/2$

    Fig. 13.  Local instantaneous vortex contours of undulating hydrofoils in array arrangement for $\omega $ = 0.75 at two moments: (a) $\theta = {\text{π}}/2$; (b) $\theta = 3{\text{π}}/2$.

    表 1  本文力系数结果与Sui等[21]的参考结果比较

    Table 1.  Comparison of the force coefficients between the present method and Sui’s method[21].

    方法平均阻力系数均方根升力系数
    150 × 800.2000.91
    300 × 1500.1841.15
    600 × 3000.1821.16
    Sui等[21]0.1801.20
    下载: 导出CSV
  • [1]

    Zabihi M, Lari K, Amiri H 2017 J. Mech. Sci. Technol. 31 3539

    [2]

    Yang J M 2016 J. Hydrodyn. Ser. B 28 713

    [3]

    王力, 田方宝 2018 中国科学: 物理学 力学 天文学 48 094703

    Wang L, Tian F B 2018 Sci. China Phys. Mech. Astron. 48 094703

    [4]

    Al-Marouf M, Samtaney R 2017 J. Comput. Phys. 337 339

    [5]

    Huang W X, Chang C B, Sung H J 2011 J. Comput. Phys. 230 5061

    [6]

    Tullio M D D, Pascazio G 2016 J. Comput. Phys. 325 201

    [7]

    吴晓笛, 刘华坪, 陈浮 2017 物理学报 22 224702

    Wu X D, Liu H P, Chen F 2017 Acta Phys. Sin. 22 224702

    [8]

    Yeo K S, Ang S J, Shu C 2010 Comput. Fluids 39 403

    [9]

    Tian F B, Wang W, Wu J, Sui Y 2016 Comput. Fluids 124 1

    [10]

    Bergmann M, Iollo A, Mittal R 2014 Bioinspiration Biomimetics 9 046001

    [11]

    Khalid M S U, Akhtar I, Dong H 2016 J. Fluids Struct. 66 19

    [12]

    Khalid M S U, Akhtar I, Imtiaz H, Dong H, Wu B 2018 Ocean Eng. 157 108

    [13]

    Mittal R, Dong H, Bozkurttas M, Najjar F M, Vargas A, Loebbecke A V 2008 J. Comput. Phys. 227 4825

    [14]

    王亮, 吴锤结 2011 力学学报 43 18

    Wang L, Wu C J 2011 Chin. J. Theor. Appl. Mech. 43 18

    [15]

    王文全 2014 计算力学学报 31 646

    Wang W Q 2014 Chin. J. Comput. Mech. 31 646

    [16]

    辛建建, 石伏龙, 金秋 2017 物理学报 66 186

    Xin J J, Shi F L, Jin Q 2017 Acta Phys. Sin. 66 186

    [17]

    Xin J J, Shi F L, Jin Q, Lin C 2018 Comput. Fluids 176 210

    [18]

    Leer B V 1979 J. Comput. Phys. 32 101

    [19]

    Van d V H A, Kohn R V, Iserles A, Ciarlet P G, Wright M H 2006 Iterative Krylov Methods for Large Linear Systems (Cambridge: Cambridge University Press) pp133–135

    [20]

    Liu G R, Gu Y T 2005 An Introduction to Meshfree Methods and Their Programming (Netherlands: Springer) p73

    [21]

    Sui Y, Chew Y T, Roy P, Low H T 2007 Int. J. Numer. Methods Fluids 53 1727

  • [1] 辛建建, 石伏龙, 金秋. 一种径向基函数虚拟网格法数值模拟复杂边界流动. 物理学报, 2017, 66(4): 044704. doi: 10.7498/aps.66.044704
    [2] 程玉民, 戴保东. 势问题的径向基函数——局部边界积分方程方法. 物理学报, 2007, 56(2): 597-603. doi: 10.7498/aps.56.597
    [3] 司马文霞, 刘 凡, 孙才新, 廖瑞金, 杨 庆. 基于改进的径向基函数神经网络的铁磁谐振系统混沌控制. 物理学报, 2006, 55(11): 5714-5720. doi: 10.7498/aps.55.5714
    [4] 刘飞飞, 魏守水, 魏长智, 任晓飞. 基于速度源修正的浸入边界-晶格玻尔兹曼法研究仿生微流体驱动模型. 物理学报, 2014, 63(19): 194704. doi: 10.7498/aps.63.194704
    [5] 吴晓笛, 刘华坪, 陈浮. 基于浸入边界-多松弛时间格子玻尔兹曼通量求解法的流固耦合算法研究. 物理学报, 2017, 66(22): 224702. doi: 10.7498/aps.66.224702
    [6] 李强, 李五明. 带嵌件型腔内熔接过程的数值模拟研究. 物理学报, 2016, 65(6): 064601. doi: 10.7498/aps.65.064601
    [7] 张妮, 刘丁, 冯雪亮. 直拉硅单晶生长过程中工艺参数对相变界面形态的影响. 物理学报, 2018, 67(21): 218701. doi: 10.7498/aps.67.20180305
    [8] 巫梦丹, 周胜林, 叶安娜, 王敏, 张晓华, 杨朝晖. 基于中性水凝胶/取向碳纳米管阵列高电压柔性固态超级电容器. 物理学报, 2019, 68(10): 108201. doi: 10.7498/aps.68.20182288
    [9] 刘佳, 徐玲玲, 张海霖, 吕威, 朱琳, 高红, 张喜田. 一步水热法在Al掺杂ZnO纳米盘上可控自组装合成ZnO纳米棒阵列. 物理学报, 2012, 61(2): 027802. doi: 10.7498/aps.61.027802
    [10] 于涛, 尹成友, 刘汉. 基于特征基函数的球面共形微带天线阵列分析. 物理学报, 2014, 63(23): 230701. doi: 10.7498/aps.63.230701
    [11] 白安琪, 胡迪, 丁武昌, 苏少坚, 胡炜玄, 薛春来, 樊中朝, 成步文, 俞育德, 王启明. 用PAA模板法实现硅基纳米孔阵列结构. 物理学报, 2009, 58(7): 4997-5001. doi: 10.7498/aps.58.4997
    [12] 李宏斌, 刘文清, 张玉钧, 丁志群, 赵南京, 魏庆农, 王玉平, 杨立书. 基于径向基函数网络的激光诱导荧光特征光谱分离算法. 物理学报, 2005, 54(9): 4451-4457. doi: 10.7498/aps.54.4451
    [13] 武浩, 朱拓, 孔艳, 陈卫, 杨建磊. 基于径向基函数神经网络的荧光光谱技术在菌种识别中的应用. 物理学报, 2010, 59(4): 2396-2400. doi: 10.7498/aps.59.2396
    [14] 曾喆昭. 不确定混沌系统的径向基函数神经网络反馈补偿控制. 物理学报, 2013, 62(3): 030504. doi: 10.7498/aps.62.030504
    [15] 申金媛, 刘, 常胜江, 贾 佳, 张文伟, 张延, 母国光. 基于径向基函数的多目标旋转不变分类及其光电实现. 物理学报, 1998, 47(12): 1968-1975. doi: 10.7498/aps.47.1968
    [16] 魏德志, 陈福集, 郑小雪. 基于混沌理论和改进径向基函数神经网络的网络舆情预测方法. 物理学报, 2015, 64(11): 110503. doi: 10.7498/aps.64.110503
    [17] 刘 丁, 任海鹏, 孔志强. 基于径向基函数神经网络的未知模型混沌系统控制. 物理学报, 2003, 52(3): 531-535. doi: 10.7498/aps.52.531
    [18] 郭会军, 刘君华. 基于径向基函数神经网络的Lorenz混沌系统滑模控制. 物理学报, 2004, 53(12): 4080-4086. doi: 10.7498/aps.53.4080
    [19] 陈帝伊, 柳烨, 马孝义. 基于径向基函数神经网络的混沌时间序列相空间重构双参数联合估计. 物理学报, 2012, 61(10): 100501. doi: 10.7498/aps.61.100501
    [20] 杨青, 曹曙阳, 刘十一. 基于浸入式边界方法的串联双矩形柱绕流数值模拟. 物理学报, 2014, 63(21): 214702. doi: 10.7498/aps.63.214702
  • 引用本文:
    Citation:
计量
  • 文章访问数:  319
  • PDF下载量:  7
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-11-08
  • 修回日期:  2019-12-10
  • 刊出日期:  2020-02-01

基于直角网格法的单个和阵列布置下柔性水翼绕流数值模拟

  • 1. 宁波大学海运学院, 宁波 315211
  • 2. 长沙理工大学水利工程学院, 长沙 410114
  • 通信作者: 陈振雷, chenzhenlei@nbu.edu.cn
    基金项目: 国家级-国家青年科学基金项目(51909124)

摘要: 研究柔性水翼在不可压缩流体中的水动力特性, 对于船舵和减摇鳍等海洋结构物的设计和性能优化具有重要意义. 本文将自主开发的径向基函数虚拟网格法求解器扩展到模拟绕单个或多个柔性水翼的不可压缩流动问题. 数值模型基于虚拟网格有限差分法考虑浸入边界对流场的影响, 引入紧支径向基函数(compact supported radial basis function, CSRBF)以物面Lagrangian质点追踪复杂的柔性动边界. 基于该方法, 首先模拟了均匀流中主动拍动的柔性水翼, 升阻力系数良好的网格收敛性结果验证了本文方法的精度和可靠性. 并研究了柔性水翼在不同振荡频率下的水动力特性, 阐述了柔性水翼的推力生成机制. 然后模拟了绕阵列布置柔性水翼的流动现象, 研究了不同间距和不同振荡频率下水翼表面的升阻力系数变化规律和尾涡特性, 观察到紧密布置的柔性水翼在高频振荡下推力系数存在显著的放大效应, 同时推力为零时的临界频率提前.

English Abstract

    • 鱼类游动是学者和工程师开发先进高效推进设备的重要灵感来源. 例如, 仿生机器鱼以类似游鱼摆动尾鳍的方式获得推力, 相比传统螺旋桨推进器, 具有机动性强、噪声低及隐蔽性好等优点, 在海洋勘测和水下侦查等民用和军用等领域具有重要的应用价值. 另外, 许多鱼类趋向于成群游动, 通过从邻近同伴产生的涡流中获得能量, 以提高游动效率和节省体力. 因此, 研究游鱼或鱼群在水中的游动特性, 阐明鱼尾鳍横向摆动的推力产生机制和水动力相互作用机制, 有助于船舵、减摇鳍和仿生推进器等海洋结构物的设计和性能优化, 并为小型自主潜航器或无人机群的最优排布提供技术指导. 然而, 游鱼或鱼群游动是一个复杂的水动力现象, 涉及周期性的旋涡脱落, 涡流与鱼体之间的能量交换以及游鱼之间的水动力相互干涉作用.

      随着计算机硬件和数值算法的迅速发展, 计算流体力学(computational fluid dynamics, CFD)技术为鱼类等柔性翼型结构物的水动力特性研究提供了强有力的手段. 然而, 数值模拟此流固耦合现象对于CFD研究是极具挑战性的, 由于涉及复杂几何形状、界面大变形等难点. 典型的求解技术是贴体网格方法[1], 界面附近网格随物体一起运动以精确地描述运动界面, 适合于模拟高雷诺数流动. 然而, 贴体网格生成对经验要求高, 拟合柔性多体运动边界涉及复杂的网格生成和重构过程, 也会极大增加计算量和引入额外数值误差. 另一方面, 浸入边界法[2]在近些年受到极大的关注. 浸入边界法在固定背景(一般为直角)网格上求解Navier-Stokes方程, 通过在流体动量方程中附加力源项[3]或重构浸入边界附近插值模板[4]以考虑浸入边界对流场的影响. 其网格生成简单, 无需网格重构, 适于模拟复杂几何形状、多体动边界流动问题. 例如, Huang 等[5]采用惩罚浸入边界法, Tullio 和Pascazio [6]采用直接力浸入边界法研究柔性丝带的流致振荡特性, 吴晓笛等[7]采用基于浸入边界的格子玻尔兹曼求解器模拟了双圆柱两自由度涡激振动现象.

      针对鱼类游动的水动力机理研究, 许多学者基于直角网格方法对绕单个或多个鱼形水翼的流动问题进行了数值模拟, 并取得卓有成效的进展. Yeo等[8]基于混合直角网格方法研究了单个鱼体的自主推进性能和鱼体大变形时的操纵性, 他们展示了仿生鲤鱼在体长范围内执行70°急转弯的能力. Tian等[9]基于贴体动网格技术的空间时间有限元法模拟了一对母子鱼在均匀流中的振荡特性, 研究了鱼体间相对距离、雷诺数、相位差和相速度对力系数、能效和涡结构的影响. Bergmann等[10]基于虚拟区域法模拟了单个和并列布置两鱼体的游动特性. Khalid等[11,12]采用Mittal等[13]的多维虚拟网格法研究了串列布置下非同步振荡的两游鱼的水动力性能, 他们发现不同振荡频率的两游鱼产生非线性相互作用, 生成了新的水动力分量, 而且, 在下游鱼体的导边处尾涡会发生分裂, 导致上游鱼体的阻力降低. 国内方面, 王亮和吴锤结[14]研究了三鱼群游动以揭示“槽道效应”在鱼群游动中的节能机制, 王文全[15]研究了鱼体自主游动的尾涡模式和鱼体运动特性.

      目前, 学者对单鱼体、串列或并列布置的双鱼体的流体特性进行了大量研究. 然而, 尚未知串列或并列布置时形成的力模式和尾涡结构是否适用于鱼群游动, 当并列和串列布置同时存在时是否有新的力模式和尾涡结构出现. 因此, 研究流体中以阵列布置的鱼群游动能更好地再现鱼群游动的水动力特性. 鉴于此, 本文拟数值研究单个或阵列布置下鱼形水翼在流场中的水动力特性, 以揭示鱼体游动的推力产生机制和鱼群游动的水动力优势. 数值模型基于自主开发的径向基函数虚拟网格法[16,17], 在固定直角交错网格上求解不可压缩Navier-Stokes方程, 引入紧支径向基函数到虚拟网格法以处理任意复杂的柔性边界, 一个面积分数表示方法通过提高浸入界面的局部质量守恒性以抑制动边界的压力振荡问题. 采用该方法模拟了不同激励频率下以鱼形游动的水翼在均匀流中的水动力特性, 并进行了网格收敛性的验证分析. 然后, 研究了9个以阵列布置的鱼形水翼在不同间距和不同激励频率下的力系数变化规律和尾涡形态.

    • 非定常黏性不可压缩流体的控制方程如下:

      $\nabla \cdot {{u}} = 0,$

      $\frac{{\partial {{u}}}}{{\partial t}} + {{u}}\nabla \cdot {{u}} = - \frac{1}{\rho }\nabla p + \nabla \cdot \left( {\nu \nabla {{u}}} \right),$

      ${{u}}\left| {_\varGamma } \right. = {{V}},~~ \frac{{\partial p}}{{\partial {{n}}}}\left| {_\varGamma } \right. = - {{n}} \cdot {{a}},$

      其中${{u}}$是速度向量, 具有直角网格分量$\left( {u, v} \right)$, $\rho $是密度, $t$是时间, $p$是压力, $\nu $是流体的运动黏性系数, Va是固体的速度和加速度; n是固体表面的外法向向量; $\varGamma $表示固体表面.

      流体控制方程(1)和方程(2)采用有限差分法在交错直角网格上离散, 采用虚拟网格法施加浸入动边界的无滑移边界条件(3)式. 在本文求解器中, 以分步法解耦速度和压力, 以二阶Total Variation Diminishing Runge-Kutta(TVD-RK2) 格式进行时间积分, 采用时间显式格式处理对流项, 黏性项以Crank–Nicolson格式半隐式处理. 对于空间离散格式, 采用TVD-MUSCL(Monotonic Upstream-centered Scheme for Conservation Laws)[18]格式离散对流项, 黏性项以中心差分格式离散. 采用不完全Cholesksy预处理的双共轭梯度法(Incomplete Cholesksy Biconjugate gradient stabilized method, BiCGSTAB)[19]求解压力泊松方程离散形成的大型稀疏线性方程组.

    • 为了考虑静止或运动物体对流场的影响, 本文采用虚拟网格法. 该方法的基本理念是, 通过在浸入物体内部布置有限数量的虚拟网格以表示静止或运动边界的存在, 以适当的插值模块及技术(如双线性插值方法), 施加浸入动边界条件.

      虚拟网格的插值重构过程可细分为三步: 标记镜像点、插值重构和满足镜像条件. 首先, 对于每个虚拟网格(GC), 沿着物体表面的外法向从虚拟网格单元中心延伸到流体域以确定唯一的一个镜像点(MP), 如图1所示. 然后是采用双线性插值格式计算每个镜像点的流体变量, 为了避免内插模板点是虚拟网格自身的情况, 镜像点的计算需要分如图1中的三种情况进行讨论. 最后是根据镜像点变量通过满足镜像条件得到虚拟网格变量值, 更多细节参见文献[17]. 需要注意的是由于交错网格布置, 两个方向速度分量或压力的离散单元位置并不重合, 因此, 需要对每个方向的速度分量或压力分别应用以上的重构过程. 在整个计算域(液体和固体)求解Navier-Stokes方程, 边界条件便得到隐式满足.

      图  1  虚拟网格法的二维插值格式 (a)无虚拟网格; (b)一个虚拟网格; (c)两个虚拟网格

      Figure 1.  Two-dimensional interpolation stencil of the ghost cell method: (a) No ghost cells; (b) one ghost cell; (c) two ghost cells.

    • 浸入边界法的网格生成过程简单, 但由于浸入边界本身并未包含在计算网格中, 精确表示背景网格中浸入边界的存在是个难点. 本文采用基于CSRBF的隐式等值面方法重构任意界面, 以光滑的基函数拟合离散的物面控制点构造任意复杂物体的等值曲面, 根据等值面的特性识别场点属性. 本文CSRBF界面描述方法简单、以有限的物面控制点获得较高的精度, 适于描述不规则、大变形甚至凹形界面.

      针对任意浸入界面, 基于有限数量的界面位置信息的离散控制点, 引入径向基函数构造如下场函数:

      $\psi \left( {{{{x}}_i}} \right) = - \left[ {\sum\nolimits_{j = 1}^M {{\beta _j} \cdot R\left( {\left\| {{r_{ij}}} \right\|} \right)} + P\left( {{{{x}}_i}} \right)} \right],$

      其中M为物面控制点的数量, ${\beta _j}$为权重系数, $P\left( {{{{x}}_i}} \right)$为线性多项式. 在文献[16]中, 我们采用MQ (Multi-quadrics)核函数拟合复杂刚性界面. 然而, 当涉及复杂薄壁边界时, 由于矩阵病态, 拟合过程可能失败. 为了解决该问题, 本研究采用紧支Wendland’s C2函数[20], 得到的矩阵为正定稀疏矩阵, 矩阵的病态性得到缓解, 收敛性得到提高. Wendland’s C2函数表达如下:

      $R\left( {\left\| r \right\|} \right) = \left\{ \begin{aligned} & {\left(1 - \dfrac{r}{\delta }\right)^4}\left(1 + 4\frac{r}{\delta }\right),\;\;\;\;r \leqslant \delta \\ & 0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\quad\;\;\quad\;r > \delta, \end{aligned} \right.$

      其中搜索距离$r = \sqrt {{{(x - {x_j})}^2} + {{(y - {y_j})}^2}} $, $\delta $是基函数的支持半径, 当前研究中令$\delta = 2\max (\Delta x, \Delta y)$, $\Delta x$$\Delta y$分别是xy方向的网格间距. 该当描述物体表面的等值面函数得到后, 背景网格相位状态如流体网格、固体网格和虚拟网格可根据等值面函数特性直接识别.

    • 本文径向基函数虚拟网格法的计算程序采用Fortran 90语言编写, 数值算法的流程图见图2. 以圆柱绕流、翼型绕流和圆柱振荡等算例验证了本文方法的精度、可靠性和模拟刚性动边界的能力[16,17].

      图  2  静止或运动边界流动求解算法的流程图

      Figure 2.  Flowchart of numerical algorithm for stationary or moving boundary flows.

    • 这里模拟了均匀流中柔性水翼的横向振荡运动以验证本文方法模拟柔性动边界流动的能力. 水翼剖面为NACA0012机翼, 采用径向基函数拟合141个样品点进行追踪. 为了模拟水中游鱼的运动, 柔性水翼的横向振荡采用一个波浪函数描述:

      $y(x,t) = A(x)\sin (kx - \omega t), \;x \in \left[ {0,1} \right],$

      $y(x, t)$是水翼表面质点在y方向的位移, $\omega $是振荡角频率, 波数$k = 2{\text{π}}$. 振荡幅值$A(x)$表示为$A(x) = $$ {A_1} + {A_2}x + {A_3}{x^2}$, 其中系数分别为${A_1} = 0.04$, ${A_2} = $$ - 0.16$${A_3} = 0.32$. 因此, 水翼运动在导边和随边的最大幅值分别为0.04和0.2 m. 雷诺数 (Re)和Strouhal数(St)是两个重要的流动参数, Strouhal数定义为${S_t} = 2{A_{\max }}f/{U_\infty }$, 其中${A_{\max }}$是随边的最大幅值, $f$是振荡幅值($f = \omega /2{\text{π}}$). 根据Sui等[21]的工作, 雷诺数设置为Re = 400, 选择了五个振荡角频率($\omega $ = 0.15, 0.30, 0.45, 0.60和0.75)以研究振荡频率对力系数的影响, 分别对应的Strouhal数为St = 0.191, 0.382, 0.573, 0.764和0.955. 计算区域设置为[–3b, 9b] × [–3b, 3b], 翼弦长度b = 1 m. 入口边界为均匀流, 均匀流速度${U_\infty }$ = 0.05 m/s, 在出口采用出流边界条件, 上下边界采用辐射边界条件. 对于压力, 在出口指定压力为零, 其余边界采用Neumann边界条件, 选择固定的时间步长$\Delta t$ = 0.002 s.

      首先对激励频率$\omega $= 0.45时鱼体振荡运动进行了网格收敛性研究, 采用非均匀网格对局部区域[–0.2b, 1.3b] × [–0.3b, 0.3b]进行网格细化, 细化区域分别采用150 × 80, 300 × 150, 600 × 300, 对应在x-方向的最小网格间距$\Delta x$分别为0.01b, 0.005b和0.0025b, 对应在y-方向的最小网格间距$\Delta y$分别为0.0075b, 0.004b和0.002b. 表1比较了本文三套网格下计算得到的力系数与Sui等[21]的参考结果. 从表1知, 随着网格细化, 平均阻力系数和均方根升力系数逐渐趋近Sui等[21]的Lattice Boltzmann方法结果, 而且, 两套细网格上的结果基本重合, 说明本文方法获得收敛解, 其中中间网格(300 × 150)上计算的平均阻力系数相比参考结果略大2.22%, 均方根升力系数略小4.17%.

      方法平均阻力系数均方根升力系数
      150 × 800.2000.91
      300 × 1500.1841.15
      600 × 3000.1821.16
      Sui等[21]0.1801.20

      表 1  本文力系数结果与Sui等[21]的参考结果比较

      Table 1.  Comparison of the force coefficients between the present method and Sui’s method[21].

      鉴于300 × 150的中间网格上已计算得到相对合理的收敛解, 将300 × 150的中间网格用于其他工况的计算. 图3将五个频率下的力系数与Sui等[21]的Lattice Boltzmann方法的计算结果进行了比较. 平均阻力系数和均方根升力系数与参考结果合理地符合. 对于角频率$\omega \leqslant 0.45$, 本文力系数与参考结果较好地符合, 当$\omega > 0.45$, 可观察到力系数与参考结果的微小偏差. 以角频率$\omega = 0.75$为例, 本文平均阻力系数略大于参考结果7.1%, 均方根升力系数略小于参考结果3.3%. 从图3可观察到, 随着振荡频率的增加, 平均阻力系数迅速降低. 临界频率为$\omega = 0.6$, 此时平均阻力系数从正值变为负值, 即水翼在流体中运动所遭受的阻力变为推力, 这是由于高频振荡产生的推力即为水中鱼类向前运动的动力来源. 当角频率从0.6增加到0.75时, 推力系数从大约为0显著增加到0.21. 另外, 当角频率从0.15增加到0.75时, 均方根升力系数呈现近乎线性的增加趋势(从0.21 到 3.1).

      图  3  五个频率下的力系数比较 (a) 平均阻力系数; (b) 均方根升力系数

      Figure 3.  Comparison of force coefficients under five oscillation frequencies: (a) Average drag coefficient; (b) root mean square lift coefficient.

      图4展示了三个角频率下($\omega = 0.15$, 0.45, 0.75)升阻力系数随无因次时间的变化趋势. 初始阶段经历一些不稳定的波动后, 三个频率下的升力和阻力系数都呈现出周期性的变化趋势. 升力系数的振荡周期是阻力系数的两倍, 该特性与圆柱绕流的卡门涡街现象类似. 随着振荡频率从$\omega = 0.15$增加到$\omega = 0.75$, 平均阻力系数的平衡位置从0.36降低到–0.3左右, 但是升力系数的振荡幅值从0.4显著增加到4.3左右, 与图3的平均阻力系数和均方根升力系数的变化趋势一致. 对于$\omega = 0.15$, 阻力和升力系数呈现小幅值和规则的正弦曲线变化. 当频率增加到$\omega = 0.45$, 阻力系数曲线的波峰变得略微平坦, 波谷变得较为尖锐, 同时, 升力系数曲线的正弦曲线特性发生轻微偏斜. 当频率进一步增加到$\omega = 0.75$, 阻力系数观察到更加平坦的波峰和更加尖锐的波谷, 类似于脉冲形状, 而且, 升力系数曲线观察到更明显的波峰波谷偏斜现象.

      图  4  三个频率下升阻力系数随无因次时间的变化趋势 (a) 阻力系数; (b) 升力系数

      Figure 4.  Variation of drag and lift coefficients with the dimensionless time under three oscillation frequencies: (a) Drag coefficient; (b) lift coefficient.

      进一步分析了柔性水翼振荡的水动力特性. 图5图7分别展示了在三个振荡频率下相位角$\theta = {\text{π}}/2$$3{\text{π}}/2$时的局部瞬时涡等值线, 根据(5)式, 相位角$\theta = {\text{π}}/2$$3{\text{π}}/2$分别对应随边位移处于正的和负的最大位移时. 由于柔性水翼周期性的摆动, 沿着水翼表面的剪切流变得不稳定, 一对涡从水翼表面交替脱落并沿着流动方向流向尾流远方. 当水翼尾部在前半个周期摆动到正的最大位移时, 一个顺时针旋转的涡从水翼上表面生成. 当水翼尾部在后半个周期摆动到负的最大位移时, 一个逆时针旋转的涡从水翼下表面脱落. 这对从水翼尾部交替脱落的反向旋转的涡列就是卡门涡街, 因而产生周期性的升力系数变化. 在这个涡列的尾流形成一个能量缺损区, 由此产生一个正向阻力. 随着振荡频率的增加, 观察到$\omega = 0.45$时尾涡变换为逆von-Karman涡街形态, 此时, 在水翼尾流区产生一个能量增强区, 阻力降低, 当振荡频率增加到$\omega = 0.6$时, 阻力甚至降低到负值, 即为推力, 如图3所示. 在更高的振荡频率如$\omega = 0.75$时, 水翼快速的拍动运动导致逆卡门涡街过渡为偏斜尾涡, 如图5所示. 另外, 随着振荡频率增加, 附着在水翼尾部的涡列随着水翼的加速运动以更快的速度脱落, 如图7所示. 在$\omega = 0.15$, 0.45, 0.75时分别可观察到两对、四对和七对涡列, 当水翼以高频振荡时($\omega = $ 0.75), 迅速脱落的涡列在尾流区紧致排列, 而且, 涡之间发生相互排挤, 导致涡形状的变形, 在尾流远方, 涡列甚至发生偏斜.

      图  5  $\omega $ = 0.15时在两个时刻的局部瞬时涡等值线 (a) $\theta = $$ {\text{π}}/2$; (b) $\theta = 3{\text{π}}/2$

      Figure 5.  Local instantaneous vortex contours for $\omega $ = 0.15 at two moments: (a) $\theta = {\text{π}}/2$; (b) $\theta = 3{\text{π}}/2$.

      图  7  $\omega $ = 0.75时在两个时刻的局部瞬时涡等值线 (a) $\theta = $${\text{π}}/2$; (b) $\theta = 3{\text{π}}/2$

      Figure 7.  Local instantaneous vortex contours for $\omega $= 0.75 at two moments: (a) $\theta = {\text{π}}/2$; (b) $\theta = 3{\text{π}}/2$.

      图  6  $\omega $ = 0.45时在两个时刻的局部瞬时涡等值线 (a) $\theta = $$ {\text{π}}/2$; (b) $\theta = 3{\text{π}}/2$

      Figure 6.  Local instantaneous vortex contours for$\omega $ = 0.45 at two moments: (a) $\theta = {\text{π}}/2$; (b) $\theta = 3{\text{π}}/2$.

    • 本节模拟了9个以阵列布置水翼在均匀流中的振荡运动以研究鱼群游动中振荡水翼表面的力模式和尾涡形态. 为了简化计算模型, 研究聚焦于中心位置的水翼, 该水翼以(6)式指定的模式进行振荡运动, 其余水翼保持静止. 几何模型如图8所示, 计算区域设置为[–4b, 12b] × [–4b, 4b]. 水翼之间的横向间距和纵向间距相同, 研究了三个间距下鱼群的力系数特性, 间距G分别为G = 0.25, 0.5和1 m. 根据3.1节的网格收敛性结果, 采用均匀网格对水翼附近[–2.5b, 2.5b] × [–1.5b, 1.5b]的区域进行细化, 在xy方向的最小网格间距分别为0.005 和0.004 m, 对应于3.1节的收敛网格. 边界条件和计算参数的设置类似于3.1节, 本节研究了五个振荡角频率($\omega $ = 0.15, 0.30, 0.45, 0.60和0.75)下中间位置水翼的力模式和尾涡形态.

      图  8  阵列布置水翼游动的几何模型

      Figure 8.  Geometrical model of the undulating hydrofoil in array arrangement.

      首先对中间位置水翼表面的升阻力系数进行了定量分析. 图9比较了三个位置下群游中间位置水翼和单水翼游动中振荡频率和水翼表面力系数之间的关系. 在水翼间距G = 1 m时, 中间位置水翼的平均阻力系数和升力均方根系数与3.1节单水翼绕流算例近似, 由于距离远, 邻近水翼对中间位置水翼的影响较为微弱, 中间水翼保持相对独立的游动模式. 当间距缩小为G = 0.5 m时, 水翼之间存在显著的水动力相互作用, 力系数对振荡频率的敏感性增加, 当振荡频率从ω = 0.15增加到ω = 0.75时, 平均阻力系数从 0.32降低到–0.38, 升力均方根系数从0.41增加到4.11. 当缩小到最窄间距G = 0.25 m时, 此时水动力相互作用最为显著, 在振荡频率从$\omega = 0.15$增加到$\omega = 0.75$时, 平均阻力系数从0.32降低到–0.67, 升力均方根系数从0.75增加到10.10. 值得注意的是, 相比单水翼绕流, 水翼表面阻力变为推力的临界频率从大约0.60减小到0.45左右, 说明水翼群游时, 水翼可以以更少的能量消耗获得同样推力. 在低振荡频率ω = 0.15时, 三个间距下的水动力相互作用力系数变化不大, 因为此频率下脱落的尾涡结构相对稳定独立, 受到周围水翼的影响小. 在振荡频率ω = 0.75时, G = 0.25 m时群游中间水翼相比单水翼绕流可获得2.5倍的平均推力增加, 升力均方根系数可获得3倍增加, 说明增大激励频率可放大水翼之间的水动力相互作用, 也说明群游水翼在相同的运动形式下能获得更高的游动效率, 这反映出鱼类群游的水动力学优势.

      图  9  三个间距下群游中间位置水翼和单水翼在五个频率下的力系数比较 (a)平均阻力系数; (b)均方根升力系数

      Figure 9.  Comparison of force coefficients between the central hydrofoil in array arrangement for three gap distances and the single hydrofoil under five frequencies: (a) Average drag coefficient; (b) root mean square lift coefficient.

      图10展示了G = 0.5 m时群游中间位置水翼游动中三个频率下升力和阻力系数随无因次时间的变化趋势. 群游水翼的升阻力系数变化趋势与单水翼类似, 然而相比单水翼游动, 随着激励频率的增加, 升力系数以更大的幅值和更大的偏斜呈现周期性的振荡特性, 阻力系数呈现更加明显的波谷平坦和波峰尖锐的现象.

      图  10  群游中间位置水翼在三个频率下升阻力系数随无因次时间的变化趋势 (a)阻力系数; (b)升力系数

      Figure 10.  Variation of drag and lift coefficients on the central hydrofoil in array arrangement with the dimensionless time under three oscillation frequencies: (a) Drag coefficient; (b) lift coefficient.

      图11图13展示了群游水翼在间距G = 0.5 m时三个频率下两个时刻的局部涡场. 中间位置水翼的游动特性和力系数受两个因素影响. 首先, 中间水翼的头部位于上游水翼尾部的低压区, 其尾部位于下游水翼头部的高压区, 此压差增加了中间水翼的推力; 其次, 中间水翼与其并列的上下两水翼形成一个半封闭区域, 两水翼之间的流体受到限制而沿着水翼游动方向流动, 减少了能量向两侧的耗散, 提高了推进效率. 当中间水翼尾部向上摆动到最大幅值时(如图11(a)), 中间水翼上表面的流体由于伯努利效应而得到加速, 导致压力降低, 下表面的流体减速, 导致压力升高, 因此, 上下表面之间压差形成向上的升力, 相比单水翼游动时增大了升力幅值, 如图4所示; 当中间水翼尾部向下摆动到最大幅值时(如图11(b)), 中间鱼体形成向下的升力. 随着Sr数增加, 鱼体之间的剧烈水动力相互作用, 扭曲了周期性的对称交替涡泄, 导致升力系数的正弦振荡特性发生显著的偏离, 如图10中的升力系数曲线. 而且, 由于激励频率增加, 相速度增加, 涡列以更快的速度脱落, 在下游鱼体的阻碍作用下, 涡列沿着下游鱼体上下表面发生分叉, 形成两列紧密布置的涡列向远处, 如图12图13所示.

      图  11  $\omega $ = 0.15时在两时刻的水翼群游局部涡场图  (a) $\theta = {\text{π}}/2$; (b) $\theta = 3{\text{π}}/2$

      Figure 11.  Local instantaneous vortex contours of undulating hydrofoils in array arrangement for $\omega $ = 0.15 at two moments: (a) $\theta = {\text{π}}/2$; (b) $\theta = 3{\text{π}}/2$.

      图  12  $\omega $ = 0.45时在两时刻的水翼群游局部涡场图 (a) $\theta = $$ {\text{π}}/2$; (b) $\theta = 3{\text{π}}/2$

      Figure 12.  Local instantaneous vortex contours of undulating hydrofoils in array arrangement for $\omega $ = 0.45 at two moments: (a) $\theta = {\text{π}}/2$; (b) $\theta = 3{\text{π}}/2$.

      图  13  $\omega $ = 0.75时在两时刻的水翼群游局部涡场图 (a) $\theta = $$ {\text{π}}/2$; (b) $\theta = 3{\text{π}}/2$

      Figure 13.  Local instantaneous vortex contours of undulating hydrofoils in array arrangement for $\omega $ = 0.75 at two moments: (a) $\theta = {\text{π}}/2$; (b) $\theta = 3{\text{π}}/2$.

    • 本文采用自主开发的径向基函数虚拟网格法研究了绕单个或阵列布置柔性水翼的不可压缩流动问题. 采用虚拟网格法在浸入边界内部布置有限数量的虚拟网格以施加物面边界条件, 引入紧支径向基函数通过物面Lagrangian质点拟合物面等值面函数以追踪任意复杂的物面边界. 在绕单柔性水翼流动计算中, 网格收敛性研究表明当网格细化时升阻力系数迅速趋近于参考结果, 同时, 不同激励频率下的升阻力系数与参考结果符合较好, 验证了本文方法的精度和可靠性. 而且, 平均阻力系数随振荡频率的增加迅速降低, 临界频率为$\omega = 0.6$, 此时水翼在流体中运动所遭受的阻力变为推力, 即水翼高频横向振荡产生了自身向前运动的推力. 同时, 在流场中观察到尾涡交替脱落的卡门涡街现象, 随着振荡频率增加, 涡列以更快的速度脱落, 涡之间甚至发生相互排挤和变形. 然后, 分析了不同间距和不同激励频率下阵列布置水翼的水动力特性, 随着间距缩小, 力系数呈现显著的放大效果, 在间距G = 0.25 m时, 推力系数的临界频率减小为$\omega = 0.45$, 说明紧密布置的水翼在相同的运动形式下能获得更高的游动效率.

参考文献 (21)

目录

    /

    返回文章
    返回