2020-09ML 目录
-
-
绕流是自然界以及实际工程运用中一个十分普遍的物理现象,也是一个经典的流体力学问题。例如:机翼在流动空气中的振动问题[1-3];导弹无侧滑大攻角飞行时,背风面分离涡的不对称运动导致的附加偏航力问题;水轮机导叶、叶片在水流作用下的摆动和振动问题[4-5];鱼在水中通过鱼鳍的摆尾提供推力[6];昆虫利用翅膀的拍动使其在空中飞行[7]等等。这些运动及工程都有涉及到翼型、椭圆在流场中的绕流问题。当流体绕过不具有流线型的钝体结构时,在钝体尾部会发生周期性的漩涡脱落,产生周期性的上下拖曳力,也就是升、阻力。随着来流情况的变化,当涡脱频率接近结构的自振频率时,漩涡脱落和结构振动相互锁定,形成共振,也就是涡激振动,其在实际生活中产生的危害十分常见并且影响也十分广泛。例如:1940年刚通车数月的塔科马大桥发生颤振风毁事件[8],震惊世界;以及近期的虎门大桥涡激振动事件[9]等。因此,流固耦合问题一直是计算流体力学领域较为关注的一类重要问题。
一般地,研究流固耦合问题时,通常采用贴体网格,这类网格质量良好且生成速度快,但是在实际计算中,为了适应不断变化的边界条件,在每一个时间步中都要重新划分网格,这样就大大增加了计算成本。鉴于这些原因,Peskin[10]提出了浸入边界法,并应用到模拟人类心脏血液流动中的瓣膜运动中。浸入边界法(immersed boundary method,IBM)在整个物理计算区域中采用一套笛卡尔网格,避免了在每一个时间步更新网格的繁琐,使用力场来代替固体边界,使流固耦合问题中的动边界处理起来更加简单和精确。因此,力源项的获得是浸入边界法求解流体问题的关键。Peskin方法,早期是将固体边界离散为一组由弹簧连接的单元,当边界发生变形时,就会产生一个意图使其回到原位置的回复力,通过边界上的回复力结合delta函数求得力源项。但该方法涉及到固体的弹性变形,主要适用于柔性体,而且该变形不易通过流场物理量求得,不具有一般性。Goldstein等[11]在此基础上发展了虚拟边界法来求解刚性边界问题,主要使用反馈力法,配合谱方法来渐进地施加所需的边界条件,该法也称反馈力法。在浸入边界法的发展过程中,Jamaludin等[12]提出了直接力法,即通过直接构造受力区的速度场来推导力源项,从而确保满足边界条件,基本实现了对边界的准确描述,成为学者们的主要研究方向之一。如Fadlum等[13]将其应用于流固耦合研究,获得了与贴体网格下计算结果一致性较好的结果;Le等[14]提出了一种隐式直接力浸入边界法,Wu等[15]、王文全等[16]提出了一种对速度进行二次修正的隐式直接力浸入边界法;Uhlmann[17]等提出了显式直接力法。直接力法相比于Peskin的经典IBM法,无须通过繁琐的delta函数[18,19]完成拉格朗日点与欧拉点上信息交换,通过直接施加浸入边界处的速度边界,由此推导出力源项,不依赖于固体的本构关系,也不进行反馈调整,故名直接力法。
按照是否使用delta函数或能否对界面进行清晰描述,可将IBM法分为扩散界面(diffused-interface)浸入边界法和锐利界面(sharp-interface)浸入边界法。理论上来讲,经典浸入边界法、反馈力法和直接力法均属于扩散界面浸入边界法。锐利界面法[20-23]不直接计算力源项来施加边界条件,主要是在界面处通过局部插值,重构边界网格点附近的流动参数(速度),以此作为边界条件将其直接施加到计算域进行求解。理论上可以具有高阶精度[24-25],在可压缩流体[26]以及激波冲击[27]问题中均得到了利用。本文将亦采用一种改进的锐利界面浸入边界法,将整个物理区域划分成纯流体区域以及包含固体的次流体区域。在流体次区域边界的法向,构造“虚拟点—受力点—垂足点”的计算结构。受力点速度由流体域内某一虚拟点和边界上垂足点之间的速度线性插值确定,实现复杂动边界的流动数值模拟。本文利用C++编程,通过圆柱绕流经典算例与其他方法的文献结果及试验结果进行对比,验证该方法的有效性和准确性。基于该方法,进一步探究不同轴长比AR、不同攻角
θ 下的椭圆柱绕流的涡结构分布特征、流线、压力等流场现象和升阻力系数、涡脱频率等水力不稳定现象。1. 锐利界面浸入边界法求解过程
1.1 控制方程
将整个物理区域(包括流体、固体)看作不可压缩牛顿流体的粘性流动。在整个计算域内,流体采用Euler描述,固体采用Lagrange描述。其连续性方程和动量方程可表示为
∇⋅u=0 (1) ∂u∂t + (u⋅∇)u=−1ρ∇p+μρ∇2u+f (2) 式中:
u 代表速度矢量(u,v),p代表流体的压力,μ 为动力黏度,ρ 为密度,f 表示力源项。1.2 时间推进
采用二阶时间分布投影法对时间进行离散,来求解控制方程(1)、(2),具体步骤如下。
(1)忽略力源项,计算考虑初始压力作用下的中间速度
u* ,即:u*−unΔt=32Hn−12Hn−1−∇pn (3) 式中:
H=u∂u∂t+v∂u∂t−1Re(∂2u∂x2+∂2u∂y2) (4) (2)忽略压力项,考虑计算所得力源项作用下的中间速度
u′ ,即:u′−unΔt=32Hn−12Hn−1+fn+1 (5) 式中:
fn+1 力源项采用如下格式进行离散:fn+1=uf−u*Δt (6) 式中:
uf 为考虑了力源项的Euler网格点速度。首先在次区域内通过算法判断和标识流体节点、固体节点,确定受力点。在任意的流体网格点周围找到相邻的4个网格点,若相邻的这4个网格点中只要有1个为固体点,那么认为该点为考虑固体边界影响的流体网格点,即受力点,如图1(a)所示;然后,构造如图1(b)所示的“虚拟点—受力点—垂足点”计算结构:过受力点F向固体边界做垂线交固体边界于垂足点B,反向延长垂线到虚拟点I,本文使点I到点F的距离等于受力点F到垂足B的距离h。首先,在局部插值得到虚拟点I的速度,再结合垂足点B的速度利用线性插值确定受力点F的速度
uf ,其中垂足点满足固体边界无滑移条件。虚拟点I的速度利用其周围4个Euler网格点速度进行双线性插值得到:
uI=4∑n=1ϕiun (7) 式中:
un 为距离虚拟点最近的4个Euler网格上的点(3点不可在同一直线上),ϕi 为常系数:ϕiT = [xIyIxIyI1]A−1 (8) 式中:
A = [x1y1x1y11x2y2x2y21x3y3x3y31x4y4x4y41] (9) 式中:下标1~4分别指距离虚拟点I最近的四个Euler网格点。
(3)考虑新的压力项,更新流场,求得下一步的速度
u ,即:un+1−u′Δt=−∇pn+1 (10) ∇2p = 1Δt∇⋅u′ (11) 式中:
pn+1 为通过解耦压力Poisson方程(11)得到的下一步压力值。1.3 空间离散
对于对流项,设
ϕ 为流场内的任意变量,对于流体区域内部点,采用如下格式进行离散:ϕ∂ϕ∂x=ϕi,j2Δx(ϕi+1,j−ϕi−1,j)−ϕ+λΔx(ϕi+1,j−3ϕi,j+3ϕi−1,j−ϕi−2,j)−ϕ−λΔx(ϕi+2,j−3ϕi+1,j+3ϕi,j−ϕi−1,j) (12) 式中:
ϕ+=12(ϕi,j+|ϕi,j|) ;ϕ−=12(ϕi,j−|ϕi,j|) ;λ=0.125 。在边界位置,对流项采用一阶迎风格式离散。对于扩散项,在边界上采用一阶迎风格式,中间节点处采用中心差分格式,即:
∂2φ∂x2=φi+1,j−2φi,j+φi−1,jΔx2,∂2φ∂y2=φi,j+1−2φi,j+φi,j−1Δy2 (13) 在对压力的离散过程中,采用五点差分离散压力泊松方程:
pi+1,j−2pi,j+pi−1,jΔx2+pi,j+1−2pi,j+pi,j−1Δy2=∂∂x(u′Δt)+∂∂y(v′Δt) (14) 对于Neumann边界条件,在边界点处采用增设虚点的方法[28]。对式(14)等号右边项,内部节点采用中心差分格式:
∂∂x(u′Δt)+∂∂y(v′Δt)=1Δt(u′i+1,j−u′i−1,j2Δx+v′i,j+1−v′i,j−12Δy) (15) 边界节点采用一阶迎风格式,如在左边界上i=0,有:
∂∂x(u′Δt)+∂∂y(v′Δt)=1Δt(u′1,j−u′0,jΔx+v′0,j+1−v′1,j−12Δy) (16) 2. 验证算例
如图2所示,计算域为一矩形区域,长×宽=15D×10D,流体次区域为边长1.5D的正方形,次区域左侧距流场入口的距离为4.25D,次区域中央浸没一个刚性圆柱,其约束方式为全固定边界,圆柱直径为D,圆心坐标为(0,0)。计算域左侧为均匀来流入口边界,采用Dirichlet边界条件,即
u=U,v=0 ;上、下两侧均为无穿透边界;右侧为自由出流边界,采用Neumenn边界,即∂u/∂x=0,∂v/∂x=0 。整个流场采用一套间距为Δx=Δy=0.025 的均匀四边形网格,固体边界采取与流体网格尺度相等的等间距离散,时间步长为Δt=0.002 。图3为
t=160 s、雷诺数Re=300时圆柱绕流的流场速度分布(范围:−1≤u/U≤1 )。从图中可看出,速度等值线表现光滑,说明计算过程中对速度的更新、修正并没有造成流场的震荡。由此可知,采用设置中间速度并不断更新速度的方法是适用的。同时,固体边界周围的速度分布也符合真实的圆柱绕流流场分布规律。说明利用这种锐利界面法描述固体对流体的耦合作用和运用双线性插值方法计算虚拟点的速度是可行的。图4为一个涡脱落周期T内圆柱尾迹涡随时间的演化过程(范围:
−8≤wzD/U≤8 ,wz 为涡量;实线为正值,虚线为负值)。从中可看出柱体表面的边界层分布、分离剪切层的卷起和圆柱后方的分离区域清晰可见。在一个周期内圆柱尾部出现明显的正、负漩涡交替脱落,向下游发展,呈现出典型的卡门涡街现象,与试验得到的卡门涡街现象一致。说明本文对边界附近的流动参数插值的计算方法及程序是可行的。图5为Re=300时尾涡脱落稳定后的升力系数及阻力系数的时程曲线,可以看出,随着时间的推进,流场逐渐趋于稳定,在流场稳定后升阻力系数都随时间的发展而呈现出稳定的周期摆动。其阻力系数、Strouhal数(St)的全时域统计结果见表1所示。其中升力系数
cL = 2Fy/ρU2D ,阻力系数cD = 2Fx/ρU2D ,Strouhal数St = fsD/U ,fs 为涡的脱落频率。对比表中结果可知,相比于文献[16,29]的数值解,本文的计算结果与文献[30]的实验值更为接近。验证了本方法及程序的准确性。3. 椭圆柱主动运动算例
3.1 计算模型
目前,以圆柱为对象的钝体绕流的研究已较为充分,对椭圆柱的关注度相对要少一些,但椭圆柱作为介于圆柱和平板之间更具有代表性的钝体,其形状更具流线型,流动阻力低,具有较好的对抗流体的作用,而且剪切边界层的卷起、尾迹涡的脱落、转捩和耗散等也很复杂,其研究对工程实际具有重要的现实意义。取计算域尺寸、边界条件、网格尺寸、时间步长等均与验证算例一致,雷诺数Re=800。不同之处在于:流体次区域中央浸没的固体为一系列不同轴长比AR=b/a的椭圆柱(见图6),其绕枢轴点(原点O)旋转摆动时攻角
θ (摆角)的变化规律由下式控制:θ={θ0sin(12βωt)0≤t≤πβωθ0πβω<t≤π ω(2−1β)θ0sin[12βωt+π(1−β)]πω(2−1β)<t≤2πω (17) 式中:
θ0 为最大攻角值,取80°;ω=2πf 为振荡角频率,f 为运动频率,β 为曲线形状参数。该椭圆柱的摆动形式由水翼在正弦振荡的基础上演变而来(当β=1 时为正弦波的半个波长)。β=2.5 代表非正弦振荡,前后0.2 s为匀速开、关时间段,中间持续0.6 s,其1个开、关周期T内的方波曲线如图7所示。3.2 不同轴长比对椭圆柱绕流的影响
图8~图10为攻角θ=0时,平均阻力系数
¯cD 、最大升力系数|cL|max 和涡脱频率fs 随椭圆柱形状轴长比的变化规律,以及流场流线情况。由图8可以看出,随着AR的不断增大,平均阻力系数
¯cD 以及最大升力系数|cL|max 均逐渐增大。由于椭圆柱垂直于水流方向的投影面积不断增大,当流体流过椭圆柱时,对椭圆柱的冲击力增大,因此椭圆柱对流体的反作用力也会增大,椭圆柱所受到来自流体的阻力¯cD 必然会增大。但是,相比阻力而言,升力的增大幅度比较小。由图9可以看出,涡脱频率fs 随着轴长比AR的变化趋势可以分成3个阶段。初期:fs 随着轴长比AR的增大而降低,当AR=0.2时,椭圆柱的涡脱频率fs 最大,约为0.6 Hz;第二阶段:当AR进入某一区域(AR=0.5~0.7)后,fs 开始稳定下来,约为0.38 Hz;后期:当AR超出这一区间后,涡脱频率fs 将再次随着AR的增大而递减,“稳定现象”消失,在AR=0.9时,fs达到最小值。由图10可知,当AR=0.2时,椭圆柱更具流线型,流动阻力低,尾部并没有明显的分离涡,流场比较平稳;当AR=0.5时,尾部流线波动开始逐渐明显;随着轴长比的增大,流态波动更加剧烈;当AR=0.7时,在后方开始产生一个尺度较小的分离涡,其位置向上偏离椭圆柱的中心。
图11为Re=800时不同轴长比工况下,同一时刻下的瞬时涡量图(
−8≤wzD/U≤8 ,wz 为涡量;实线为正值,虚线为负值)。可以看出:当AR=0.2时,尾迹涡以S型模态脱落,呈现出上侧负漩涡与下侧正漩涡交替脱落向下游演化发展的卡门涡街现象。脱落的正负漩涡之间的纵向间距y0、周期T较小,形成明显的2列,尾涡强度(尺寸)也较小,漩涡的分离点基本位于x轴上;随着轴长比AR的增大,卡门涡街规律依然明显,但漩涡强度(尺寸)和周期T逐渐增大,漩涡分离点偏离x轴,且脱落位置x0增大,正、负漩涡中心间距y0减小逐渐趋于一列。3.3 不同攻角对椭圆柱绕流的影响
图12给出了AR=0.2,Re=800时,不同攻角
θ 下,升、阻力系数的变化规律。由图12可知,当
θ 由0∘ 向90∘ 发展时,平均阻力系数¯cD 整体呈现为上升的趋势,这是由于随着攻角的逐渐增加,垂直于来流方向的正投影面积也逐渐增加,必然导致阻力的随之增大。而最大升力系数|cL|max 的变化规律不同于平均阻力系数¯cD ,|cL|max 的发展趋势分为两个阶段,以θ=45∘ 为临界攻角,当θ<45∘ 时,最大升力系数|cL|max 随着θ 的增大而增大,而当椭圆柱摆动超出临界攻角,即当θ>45∘ 时,|cL|max 便随着θ 的增大而逐渐减小。最大升力系数与平均阻力系数的比值
|cL|max/¯cD 也是呈先增加后减小趋势,但升阻比的临界攻角则为θ=25∘ 。当攻角θ 小于临界攻角时,曲线呈现出随θ 的增大而逐渐增长的趋势,此时|cL|max/¯cD<1 ,|cL|max 数值与¯cD 的数值相差不大;当θ=25∘ 时,升阻比达到最大,为1.01;当θ 超过临界攻角后,升阻比则呈现随θ 的增加而逐渐减小的趋势,这时升力在θ=45∘ 之前虽然依旧呈现增长的趋势,但是增长速率却远小于阻力的增长速率,而超过最大升力系数的临界攻角后,由于升力开始逐渐减小而阻力却仍然呈现增长的趋势,所以升阻比必然会逐渐减小。为了讨论产生这种现象的原因,分析了位于临界攻角附近的4个工况(θ=0∘,θ=25∘,θ=45∘,θ=75∘ )的流场瞬时压力场,如图13所示。从图13中也可看出,当
θ=0∘ 时椭圆柱上下表面的压力对称,压差接近零,升力系数很小;随着攻角增大,椭圆柱的上表面的压力,不管从数值和承压面积上来看都明显比下表面大,即上下表面的压差增大,故升力系数随攻角的增大而增大。当攻角发展到临界攻角θ=45∘ 时,压差达到最大,之后升力系数的值随着攻角的增大而逐渐减小。由图14漩涡脱频率曲线可以看到,涡脱频率
fs 随着攻角θ 的不断增大,整体呈现先增加再逐渐递减的趋势。其变化趋势可以分成3个阶段:初期,θ =2°~5°阶段,fs 随着攻角θ 的增大而增加,θ=2∘ 时,涡脱频率达到最大,约为fs =0.61 Hz;第二阶段,当攻角进入θ =2°~10°区域后,fs 有所稳定,递减速率较小,在0.6 Hz左右浮动;后期,当攻角θ >10°后,涡脱频率fs 将随着攻角θ 增大而逐渐减小,后期fs 曲线逐渐平缓,减小的速率逐渐减慢。图15为Re=800时不同攻角工况下椭圆柱尾迹涡的对比,由图可见,固体表面分布的边界层、分离剪切层的卷起和椭圆柱后侧的分离区,以及涡的脱落等现象都清晰可见。与圆柱绕流一样,椭圆柱绕流尾迹涡也是呈现出明显的一对正负漩涡从上下两侧交替脱落,向下游发展的卡门涡街现象。但是,椭圆柱摆动的攻角幅值对流态起决定性作用。随着攻角的增加,涡的脱落周期、尺寸、形状和演化形式等均有较大差异,涡脱落的反对称性越发明显。类似于水轮机活动导叶大攻角运行或导弹无侧滑大攻角飞行时,背流(风)面分离涡的不对称运动。攻角越大越容易产生不希望的附加偏航力,或加大水轮机活动导叶转动轴折断的风险。
随着椭圆柱摆动攻角的增大,漩涡脱落频率降低,周期T加长,漩涡尺寸加大。当攻角
θ≤30∘ 时,在每个周期内,有两个符号相反的涡从椭圆柱上下侧脱落,在下游形成交替出现的涡对,是典型的反对称S型脱落模态,如图16(a)所示,类似于静止圆柱的涡脱落方式。当θ =45°~60°时,上下侧脱落的漩涡尾部均拖着一个方向相同的子涡(2个涡核);在向下游演化发展的过程中,上侧脱落的负漩涡的子涡与母涡融合,而下侧脱落的正漩涡的子涡则产生分裂,分裂后的次正漩涡位于最上侧,靠近负漩涡,最终尾迹涡表现为3排,在下游形成了一对涡和一个单涡交替脱落的“P+S”Ⅰ型模态,如图16(b)所示。当θ≥75∘ 时,则情况相反,上侧脱落的负漩涡的子涡发生分裂,分裂后的次负漩涡位于最下侧,靠近正漩涡,而椭圆柱下侧脱落的正漩涡的子涡则与母涡则融合,在下游形成了新的“P+S”Ⅱ型模态,如图16(c)所示。4. 结 论
采用一种锐利界面(sharp-interface)浸入边界法对二维椭圆柱绕流问题进行了数值计算,并讨论了不同轴长比AR、攻角
θ 对涡结构分布特征及升阻力系数、涡脱频率等水力不稳定现象的影响,通过对比分析,得到以下结论:(1)在模拟钝体绕流时,通过与文献和试验结果的对比,本文的锐利界面浸入边界法结果与试验结果吻合较好,验证了本方法的准确性和有效性;
(2)对于椭圆柱绕流问题,随着不同轴长比的增加,升阻力系数呈递增趋势,而流场流态波动也越来越剧烈,尾迹涡脱落位置、强度(尺寸)、间距、周期也越来越大。周期加长,涡脱频率则降低;
(3)当攻角发生变化时,平均阻力系数随攻角的增大而呈递增趋势,而最大升力系数,以及最大升力系数和平均阻力比(升阻比),均表现为前期随着攻角的增大而增大,超过临界攻角后,再降低的趋势;升阻比的临界攻角为25°,而最大升力系数的临界角为45°;
(4)当攻角发生变化时,涡脱频率总体呈随着攻角增大而逐渐减小的趋势。漩涡脱落频率降低,周期加长,漩涡尺寸加大;当攻角
θ≤30∘ 时,尾涡呈反对称S型脱落模态,当45∘≤θ≤60∘ 时,尾涡呈反对称“P+S”Ⅰ型脱落模态,当θ≥75∘ 时,尾涡呈反对称“P+S”Ⅱ型脱落模态。 期刊类型引用(1)
1. 郭涛,张纹惠,王文全,罗竹梅. 基于IBM法的低雷诺数下涡激振动高质量比效应的研究. 工程力学. 2022(03): 222-232 . 百度学术
其他类型引用(3)
-
计量
- 文章访问数: 1598
- HTML全文浏览量: 840
- PDF下载量: 64
- 被引次数: 4