基于离散涡方法的粒子强度交换法模拟流场

基于离散涡方法的粒子强度交换法模拟流场

0基于随机涡流的弹性涡流模拟方法,主要特点如下该市的方法是模拟非定常流的拉格朗日方法，可以忽略离散误差。该方法具有良好的强度特性。这两个特点使得涡方法成为模拟流场的有效的工具。二维涡方法发展最成熟的是由Chorin提出的随机走步涡团法,它是用随机走步来模拟涡量的粘性扩散,具有统计平均意义,该方法已经成功地运用到了桥梁断面的流场模拟和抗风分析中,如AllanLarsen和Zhou。但是它的缺点不能忽视,对长时间历程问题,误差急剧增长,随机走步模拟扩散精度较低,且用新生涡团处理无滑移边界条件有较高的数值噪声,确定性涡方法克服了随机涡方法的缺点,而不改变涡方法作为Lagrange方法的优点。过去的二十年里,在涡方法领域内,已经取得了较大的发展:快速多极子展开的应用(Greedgard和Rohklin,Barnes和Hut)用粒子强度交换法精确地模拟粘性扩散(Degond和Mas-Gallic),以及精确施加粘性边界条件(Koumoutsakos,etal,Mas-Gallic,Leonardetal,Ploumhans,etal,Benhaddouch)。所有这些使得涡方法能够更高精度地模拟具有粘性条件的流场。例如,Koumoutsakos和Leonard,Koumoutsakos和Shiels,用确定性涡方法分别模拟了突然起动的圆柱和平板的流场。当然,他们这种精确的方法也只应用到了特殊的截面(圆形和平板)。目前,粒子强度交换法尚未应用到工程中,且国内对这种方法的研究比较少,对其在工程上的应用更是空白。本文对粒子强度交换法进行了研究并程序化,将其初步应用到桥梁截面的静风分析中。1确定波方法1.1边界条件条件涡量场的二维Navier-Stokes方程可以汇总为如下数学模型∂ω∂t+u⋅∇ω=ν∇2ωinD(1a)∇2Ψ=-ωinD(1b)u=∇×ΨinD(1c)ω(x,0)=ω0(x)inD(1d)u=Ub(t)+Ω(t)ez×(xs-xb)on∂D(1e)u→U∞at∞(1f)∂ω∂t+u⋅∇ω=ν∇2ωinD(1a)∇2Ψ=−ωinD(1b)u=∇×ΨinD(1c)ω(x,0)=ω0(x)inD(1d)u=Ub(t)+Ω(t)ez×(xs−xb)on∂D(1e)u→U∞at∞(1f)其中ω为涡量,ν、ρ分别为流体的运动粘性系数和流体的密度,Ψ为流函数,固体的平动速度为Ub(t),绕其形心xb的角速度为Ω(t),U∞是无穷远处的速度,式(1b)为泊松方程,(1c)为不可压缩条件,(1d)为初始条件,(1e)为固体边界条件,(1f)为远场边界条件。以上方程和二维Navier-Stokes方程是等效的,涡量-流函数方程使得方程中压力项消失,但是对于有边界的流体域,它引入了附加的运动约束,需要将速度边界条件转变为涡量边界条件。1.2didt、dh、ni本文中的粘性扩散采用粒子强度交换法,该方法是用积分算子来估计拉普拉斯算子,∇2ω(x)≈2σ2∫ησ(x-y)(ω(y)-ω(x))dy(2)∇2ω(x)≈2σ2∫ησ(x−y)(ω(y)−ω(x))dy(2)其中ησ(x)=(1/σ2)η(|x|/σ)‚η(s)=-1sddsζ(s)‚ζ(s)=12πexp(-s22),因此dΓidt=2υσ2Ν∑j=1ησ(xi-xj)(ΓjSi-ΓiSj)(3)Γni=Γn-1i+2υσ2Ν∑j=1ησ(xi-xj)(Γn-1jSi-Γn-1Sj)(4)在实际运用中,由于η的迅速降低,只有靠近粒子i的粒子才对dΓi/dt有明显的影响,为了减少计算量,只考虑靠近粒子i的粒子。称对dΓi/dt有明显影响的粒子集为Di,如果采用高斯光滑核,Di包括离粒子xi距离小于5σ所有粒子。Γni≅Γn-1i+2υσ2∑j∈Diησ(xi-xj)(Γn-1jSi-Γn-1iSj)(5)当所有粒子的不同时,方程(5)变为,Γni≅Γn-1i+2υσ2ij∑j∈Diηij(xi-xj)(Γn-1jSi-Γn-1iSj)(6)其中ηij(xi-xj)=(1/σ2ij)η(|xi-xj|/σij),σ2ij=(σi+σj)/2。1.3几何标准指数网格的构建粒子分布的开始阶段,像线动量和角动量等全局量一般不会有误差产生,随着时间的增长流场中的粒子会失去重叠或重叠太多,由于局部粒子的扭曲不能正确显示涡结构,而粒子的聚集导致一些非物理现象产生,从而不能精确地求解光滑的涡域。因此在粒子强度交换法中,粒子强度重分布是必不可少的一步。采用四阶内插核,满足涡量守恒,涡量的一阶矩,二阶矩和三阶矩守恒,Λ3为粒子周围没有边界的情况,Λ′2为边界位于x=-12‚Λ´3边界位于x=-1,网格间距为1,单位强度的粒子位于-12≤x≤12,新粒子的强度分别为Γ1、Γ2、Γ3、Γ4,Λ3中新粒子的位置为-32、-12、12、32‚Λ´2中新粒子的位置分别为0、1、2,Λ′3中新粒子的位置分别为-12、12、32、52,其方程分别为式(7-9)。二维示意图如图1所示。在外部流场中,涡量主要位于边界层和尾流中,涡量的大小随着涡流向下游逐渐变小。因此需要在固体边界采用更高的精度而在尾流中可以降低精度。为了实现这个目标,将具有规则序数不均匀的网格(i,j)覆盖精度逐渐变化的物理域上。当然可以考虑不同的网格形式,本文采用文献中提出的指数网格形式。粒子在流体中的面积为:S=h2≃((2πmR0)⋅exp(2πm(i+12)))2(10)h≃(2πmR0)⋅exp(2πm(i+12))(11)其中m为轴向划分的个数,R0为最小圆的半径,h为网格宽度,S为粒子在流体中的面积。2表面环量与子强度交换的模拟结果对雷诺数为550的突然起动圆柱的流场进行模拟,采用指数网格形式,m取592,其粒子的分布如图2所示(图中间距是计算所取间距的四倍即m=148),时间步取0.01,表面涡片长度为0.0053,对流和粒子强度交换模拟扩散过程均采用欧拉格式,表面涡片向流场内扩散过程采用高斯四点积分,将计算的圆柱上表面产生的环量dΓup/dt=υR∫π0∂ω∂tdθ和其他文献的结果做了比较,如图3所示,其中KL为Kou-moutsakos和Leonard方法的计算结果,G、H、V分别为Ploumhans提出的一般性方法,圆柱表面加“晕环”和粒子强度重分配时使用“随机振动”,详细参考文献,从图中可以看出本文也取得了较好的结果。同时对阻力系数和Ix=∫ΩyωdΩ分别与文献和解析解做了比较,阻力系数如图4所示,Ix的比较见图5。图6给出了突然起动圆柱的涡量等势线,并于Ploumhans的结果进行了比较,从图中可以看出计算结果比较理想。3阻力系数的时程曲线为了进一步验证程序的正确性,计算了风攻角为15°的方形流场,计算方法与圆形的方法相同,只是指数网格的圆心不在圆点,与文献中一致,取在(-2.5,-0.67),m取为850,计算的阻力系数如图7所示。图7也给出了文献中阻力系数的时程曲线,尽管阻力系数的振动幅值有所不同,但是其平均值和周期还是比较接近的。图8给出了本文与文献计算的在T=30的涡量等势线的比较。从阻力系数和涡量等势线的比较均可看出本文均取得了较为理想的计算结果。4计算误差分析大海带东桥引桥主梁截面为闭口箱形截面,其尺寸如图9所示,将粒子强度交换法首次应用到桥梁抗风计算中,本文选择大海带桥进行数值模拟,为了使得网格更为贴体,对本桥梁断面的计算选用了与笛卡尔坐标平行的方形网格,网格的间距取为0.01,无量纲时间步为0.01,雷诺数为Re=UB/υ=104,其中U为风速,B为桥宽,υ为运动粘性系数。计算的三分力系数时程如图10所示,并与随机涡方法计算结果的比较,从图中可以看出,两个计算的结果吻合得很好,但是随机涡方法计算的结果振荡比较明显,确定性涡方法计算的曲线比较光滑。取阻力系数时程稳定阶段进行平均,得到平均阻力系数,且对升力系数时程进行傅立叶变换得到St=fH/U,其中f为旋涡频率,H为截面高度。表1给出了确定性涡方法计算的阻力系数和St数与其他数值方法和实验计算结果的比较,计算结果表明确定性涡方法的结果和实验结果比较接近。图11给出了确定性涡方法计算的涡量等势线。图12为确定性涡方法计算的局部速度和随机涡方法的局部