基于投影浸入边界法的流固耦合计算模型

基于投影浸入边界法的流固耦合计算模型
王文全，梁林，闫妍
【摘要】摘要:针对建立的刚体和流体相互作用的浸入边界法数学模型，借助求解不可压缩N-S方程组的分步投影方法的思想来求解基于浸入边界法的耦合系统方程.在空间离散上，对流项采用Quick格式，扩散项采用中心差分格式.时间推进采用显式欧拉法，固体对流体的作用力源项引入流体体积(VOF)方法中的体积比函数，通过龙贝格算法高效求解每个网格单元中浸入刚体所占的体积比.利用VC++编写投影浸入边界法的数值计算程序，并以单圆柱绕流为基准数值算例，通过与其他文献和实验结果的对比，验证了数值计算结果的准确性和可靠性，并进一步分析了不同雷诺数下圆柱绕流场的涡结构分布特征.
【期刊名称】北京工业大学学报
【年(卷),期】2014(040)006
【总页数】6
【关键词】关键词:流固耦合;投影浸入边界法;分步法;圆柱绕流
非线性动态流固耦合问题是学术界和工程界最具挑战性的研究内容之一.借助现代数值手段，深入研究结构振动和非定常流场之间的相互作用，对于避免流体诱发的结构共振从而提高结构的可靠性和安全性具有重要意义［1-2］.在流固耦合中，有相当一部分涉及刚体运动与流体的相互作用，当刚体发生大位移时，对常规的基于动网格技术的流固耦合求解方法提出了潜在的挑战，如网格动态更新技术、边界层网格控制等［3-4］.而基于浸入边界法的流固耦合方法采用固定网格技术，无须在计算过程中更新网格，在具有大位移的流固耦合计算中显示出巨大优越性［5-7］.在各类浸入边界法发展过程中，结合传统分步投影
法的思想［8-10］，投影浸入边界法也引起了研究者的广泛兴趣.如Minion ［11］在考虑不可压缩流动时采用Godunov投影法求解Navier-Stokes方程，并采用光滑过的交错投影计算力源项，此方法比传统浸入边界法有更高阶的收敛精度.Griffth等［12］利用有限差分的cell-centered投影浸入边界法研究生物流体力学中的流固耦合问题.Le等［13］利用两阶投影法满足无滑移边界条件，结合笛卡儿网格，采用有限差分离散控制方程，模拟不规则区域二维圆柱绕流和三维圆球绕流问题.Taira等［14］发展了一种等同于传统分步法的浸入边界法，将浸入边界外加的边界力作为拉格朗日乘子去满足无滑移条件，通过投影分离速度场，其数值方法可得到二阶精度.Lai等［15］求解Navier-Stokes方程采用半隐式压力投影，离散过程可保持对称性，并模拟二维难溶物质与流体的相互作用.文献［16］研究了连续显力法和离散隐力法的改进方法，并比较了2类方法的求解精度.综上可见，对投影浸入边界法的研究已取得了一些丰硕的成果，但总体来说，其研究起步较晚，尤其在国内研究报道非常少见.为了深入开展这方面的基础性研究工作，本文以传统投影法为基础，建立了刚体和流体相互作用的浸入边界法数学模型和数值求解策略，并利用VC++编写投影浸入边界法的数值计算程序，以单圆柱绕流为基准数值算例，验证数值计算结果的准确性和可靠性.
1 投影浸入边界法的控制方程
1.1 数学模型
包括刚体和流体在内的整个物理区域，都假定为不可压缩牛顿流体的黏性流动，其连续性和动量方程可分别表示为
式中:u代表速度矢量;p为流体压力;Re为流体雷诺数;f表示虚拟的力源项，这
里指刚体存在对流体产生的作用力密度.
1.2 时间推进方法
借鉴Chorin分步投影求解不可压缩Navier-Stokes方程的思路［17］，求解刚体与流体相互作用的浸入边界法控制方程式(1)、(2)，具体步骤如下(以二维流动为例):
1)忽略压力项和虚拟的力源项，求得中间速度u*和v*，即
2)忽略虚拟的力源项，求解中间速度 u'和v'，即
式中 pn+1是下一时间步的压力值，通过解耦的Poisson方程求出，即
式(5)要求压力泊松方程在区域边界上严格满足Neumann边界条件，即
3)加上力源项，更新速度，得到下一时间步速度 un+1和 vn+1，即
式中η为流体网格单元中浸没刚体的体积，定义η=S R/S F.其中:S R是在当前网格单元中刚体所占的体积(面积);S F是当前网格单元的体积(面积).如图1所示，在完全被刚体浸没的单元，η=0;在纯流体单元，η=1;而在刚体与流体耦合界面穿过单元，0＜η＜1.相应可求得力源项的大小，如果在纯流场区域;在刚体存在区域其中，us、vs分别为刚体在x、y方向的速度.通过式(3)～(11)，将第n时间步的流场推进到第n+1时间步.
1.3 空间离散
设φ为流场任意变量，在内部节点，对流项采用如下格式离散，即
式中:β=0.125.在边界节点，采用一阶迎风格式.
扩散项采用中心差分格式，即
在边界上，使用一阶迎风格式，如在左边界(i=0)，其离散形式可表示为
压力泊松方程空间离散采用标准的五点差分格式，即
对于Neumann边界条件，在边界点处，采用增设虚点的方法，具体参考文献［18］.对于式(15)右边导数项，当离散节点为内部节点，采用中心差分，即当离散节点为边界节点，采用迎风格式，如在左边界(i=0)，有
2 数值算例
2.1 问题描述
如图2所示，求解区域为一矩形域，L×H=30×15，流场中浸没一刚性圆柱，圆柱圆心位置在(7.5，7.5)，圆柱直径为D=2.左边为入流边界条件，为均匀来流U∞=2，v=0，上下为壁面边界条件，即u=v=0;右边为出流边界条件网格为均匀四边形网格，网格间距为0.05.时间步长Δt=0.001.
2.2 计算结果分析
一般在Re小于40时，圆柱后会形成对称的脱落涡对，在圆柱尾部的分离气泡长度(一般定义圆柱后侧点到逆时针旋转涡后驻点之间的距离［19］)有丰富的变化.表1为本文计算出的气泡长度与其他文献实验和数值结果的比较，吻合较好，证实本文采用的数值计算方法是可靠的.
图3为计算稳定后不同雷诺数下流场的流线图.从图可见，在Re小于40时，在圆柱后侧形成上下对称的回流漩涡，随着Re的增加，分离气泡长度越来越大，但Re接近42时，圆柱后侧回流漩涡的对称性遭到破坏.随着Re的进一步增加，接近80时，圆柱后侧形成卡门涡街.
图4为Re=15、t=90 000Δt时速度、压力和涡量的等值线图.可以看出，流场计算结果在y方向保持了较好的对称性，x方向速度在圆柱断面达到最大值，在圆柱后侧涡对处速度较小.压力在圆柱前端达到最大值，并沿圆柱表面向下游逐渐减小，在逆向压力梯度作用下，接近圆柱上、下顶点时，流体脱落，形成
圆柱后侧的脱落涡对.从图4(d)可见，上下涡对保持较好的对称性，并沿流向逐渐减小.计算结果较为真实地再现了圆柱涡脱落这一物理现象和涡结构的分布特征.
图5显示了不同雷诺数下沿水平方向6个截面的x方向速度分布情况.从图可见，x方向速度分布对y=7.5断面保持具有较好的对称性，但随雷诺数的不同(尽管雷诺数差异不大)，圆柱后各断面速度大小差异较大，如在x=9断面，Re=10时，x方向速度的数值出现从正到负的转变，而对Re=3，x方向速度一直为正值，说明不同雷诺数下圆柱后侧涡对前后驻点的变化，反映涡脱落位置以及形成的涡结构强度的差异，同时也证实圆柱涡脱落特征对雷诺数极其敏感.
3 结论
1)提出了一种基于投影法的浸入边界法来处理流体－刚体耦合的计算模型.当考虑刚体对流体的作用时，通过VOF方法的体积比思想，在刚体附近的网格中，通过龙贝格数值积分来求解刚体浸没在当前网格中的体积.在网格密度比较适宜的情况下，添加的力源项可以很好地体现刚体对流体的耦合作用.
2)通过分步投影，有效地处理了带有力源项的浸入边界法的控制方程.不仅实现了速度与压力的解耦，也很好地实现了对力源项的处理，更新了压力场和速度场.
3)利用VC++程序设计语言编写数值计算程序，并以单圆柱绕流为基准数值算例，通过与其他文献和实验结果的对比，验证数值计算结果的准确性和可靠性，并进一步分析不同雷诺数下圆柱绕流的涡结构分布特征.
参考文献:
［1］张立翔，杨柯.流体结构互动理论及其应用［M］.北京:科学出版社，
2004:1-25.
［2］张阿漫，戴绍仕.流固耦合动力学［M］.北京:国防工业出版社，2011:1-20.
［3］LIEW K M，WANG W Q，ZHANG L X，et al.A computational approach for predicting hydroelasticity of flexible structures based on the pressure Poisson equation［J］.International Journal for Numerical Methods in Engineering，2007，72(13):1560-1583.
［4］WANGW Q，HE X Q，ZHANG L X，et al.Strongly coupled simulation of fluid-structure interaction in a Francis hydroturbine［J］. International Journal for Numerical Methods in Fluids，2009，60(5):515-538.
［5］MITTAL R，IACCARINO G.Immersed boundary methods ［J］.Annual Review of Fluid Mechanics，2005，37:239-261.
［6］JIC，MUNJIZA A，WILLIAMS J JR.A novel iterative direct-forcing immersed boundary method and its finite volume applications［J］. Journal of Computational Physics，2012，231(4):1797-1821.
［7］YANG JM，STEM F.A simple and efficient direct forcing immersed boundary framework for fluid-structure interactions［J］.Journal of Computational Physics，2012，231(15):5029-5061.
［8］KIRISC，KWAK D.Numerical solution of incompressible Navier-Stokes equations using a fractional-step approach［J］.Computer＆Fluids，2001，30:829-851.
［9］ GALBRAITH M，ABDALLAH S.Implicit solutions of incompressible Navier-Stokes equations using the pressure gradient method［J］.AIAA Journal，2011，49:2491-2501.
［10］SAN O，STAPLESA E.A coarse-grid projection method for accelerating incompressible flow computations［J］.Journal of Computational Physics，2013，233:480-508.
［11］ MINION M L.On the stability of Godunov-projection methods for incompressible flow［J］. Journal of Computational Physics，1996，123(2):435-449.
［12］GRIFFITH B E，PESKIN C S.On the order of accuracy of the immersed boundary method:higher order convergence rates for sufficiently smooth problems［J］.Journal of Computational Physics，2005，208(1):75-105.
［13］LE D V，KHOO B C，LIM K M.An implicit-forcing immersed boundary method for simulating viscous flows in irregular domains ［J］.Computer Methods in Applied Mechanics and Engineering，2008，197(25):2119-2130.
［14］TAIRA K，COLONIUS T. The immersed boundary method:a projection approach［J］. Journal of Computational Physics，2007，225(2):2118-2137.
［15］LAIM C，TSENG Y H，HUANG H.An immersed boundary method for interfacial flows with insoluble surfactant［J］.Journal of
Computational Physics，2008，227(15):7279-7293.
［16］谢胜百，单鹏.两种浸入式边界方法的比较［J］.力学学报，2009，41(5):618-627.XIE Sheng-bai，SHAN parison of two types of immersed boundary method in uniform incompressible flows ［J］.Chinese Journal of Theoretical and Applied mechanics，2009 ，41(5):618-627.(in Chinese)
［17］LEONARD B P. Bounded higher-order upwind multidimensional finite-volume convection-diffusion algorithms［M］∥MINKOWYCZW J，SPARROW E M.Advances in Numerical Heat Transfer.London:Taylor＆Francis，1997:432.
［18］胡建伟，汤怀民.微分方程数值方法［M］.2版.北京:科学出版社，2007:85-96.
［19］VAKIL A，GREEN S I.Two-dimensional side-by-side circular cylinders at moderate Reynolds numbers［J］.Computers＆Fluids，2011，51:136-144.
［20］TENEDA S.Experimental investigation of the wakes behind cylinders and plates at low reynolds numbers［J］.Physics Social of Japan，1956，11:302-307.
［21］ TAKAMI H，KELLER H B.Steady two-dimensional viscous flow of an imcompressible fluid past a circular cylinder［J］.Physics of Fluid Supply，1969，12:51-56.
［22］姜知亨.浸入边界法及在大涡模拟中的应用［D］.北京:清华大学航空航
天学院，2011.JIANG Zhi-heng.Immersed boundary method and it’s application in large eddy simulation［D］.Beijing:School of Aerospace of Tsinghua University，2011.(in Chinese)
(IBM)was validated by the basic numerical example of flow over a circular cylinder.Finally，the distributions of the velocity field and the pressure field as well as the vortex structures were analyzed in different Reynolds numbers.
基金项目:国家自然科学基金资助项目(11262008)
(责任编辑张蕾)。