光滑粒子流体动力学(smoothed particle hydrodynamics，SPH)对模拟带自由液面的流动和动边界问题具有强大的优势．为处理计算过程中产生的非物理性振荡，SPH方法模拟流体运动时往往会施加人工黏性来保证数值模拟的稳定性．但人工黏性并不是流体黏性项的近似，只是数值稳定和处理高马赫数下冲击效应的一种手段．在许多问题中人工黏性的高耗散性影响了流体的剪切力，给使用SPH法捕捉湍流结构造成较大困难．另外，在低雷诺数流动问题中[1]，用人工黏性替代流体的真实黏性无法得到正确的速度场；在自由表面流动或者界面流动问题中[2]，人工黏性不能完全防止流场中产生压力波动，且取值过高会导致流场中动能衰减过快而影响数值模拟的准确性．在SPH法动量方程中，建立既能体现流体的物理黏性又能保证数值的稳定性的黏性模型具有重要的理论意义和工程意义，也是当前SPH方法的一个研究难点．对于黏性流体，层流黏性无法准确捕捉湍流流体的运动，到目前为止还没有为SPH方法设计具体的湍流模型．若将纳维-斯托克斯方程的黏性项展开使用流体的物理黏性，SPH方法则须将所有粒子的应变率进行求解，计算量较大．因此很多学者将SPH与不同的湍流模型结合进行了相关研究．Gotoh等[3]将SPH和大涡模拟方法相结合跟踪了自由液面流动，并模拟了波浪与防波堤相互作用时的涡流，模拟波剖面和实验波轮廓符合较好．Shao[4]将k-ε模型运用到SPH方法中模拟了破断条件下斜坡上的消弧形波的破断．Dalrymple和Rogers[5]采用大涡模拟方法，通过引入亚粒子尺度技术[6]，提出了标准的SPH黏性公式．Bishoy等[7]将SPS方程引入二维不可压缩SPH方法以捕捉流场中的复杂流动．陈善群等[8]分别采用SPS-SPH和经典的SPH方法对二维溃坝问题进行了模拟，通过比较计算结果得出SPS-SPH对于自由表面流动问题具有较高的模拟精度．本研究基于Dalrymple和Rogers提出的SPH黏性公式[6]，开发了SPS-δ-SPH算法，分别采用人工黏性、层流黏性和SPS湍流黏性对流体黏性进行处理并对不同雷诺数下的二维圆柱绕流进行数值模拟，系统分析了不同黏性处理方法下的水动力系数的变化并与实验数据进行比较．在此基础上，采用本文算法对三维圆柱绕流进行了模拟，计算了尾流场的速度、水动力系数及自由液面的变化，并与相关文献的计算结果进行了比较，为海洋工程领域相关工程问题的研究提供参考．1 δ-SPH控制方程δ-SPH模型[9-10]假定流体正压、弱可压缩，以纳维-斯托克斯方程为参考方程，通过在连续性方程中施加适当的密度耗散(人工扩散)项以消除模拟数值压力场中的高频振荡．忽略流体黏性的δ-SPH控制方程为dρidt=-ρi∑j(uj-ui)∇WijVj+δhc0∑jψij∇WijVj；(1)duidt=-∑jmjpj+piρjρi∇iWij+g；(2)dri/dt=ui，(3)式中：粒子i和j为相邻粒子；Wij为粒子i对粒子j产生影响的核函数，∇i为粒子i的梯度算子；h为光滑长度；∇W为核函数W的梯度；Vj和mj分别为粒子i支持域中粒子j的体积和质量；ρi为粒子i的密度；ui为t时刻粒子i在位置ri处的速度；pi和pj分别为第i和j个粒子的压强；c0为人工声速，为确保密度波动低于1%ρ0，通常令c0大于10倍以上的最大粒子速度；ρ0为静止时流体的密度，也为参考密度；系数δ和α分别控制密度和速度耗散的强度；g为重力加速度；系数δ控制密度耗散的强度．式(1)中密度耗散项的系数ψij=2(ρj-ρi)rji /rij2-∇ρiL+∇ρjL，式中：rij为粒子i与j的间距；∇ρL为重整化的密度梯度，采用的是修正的粒子近似法计算，即∇ρiL=∑j(ρj-ρi)Mi-1∇WijΔVj，Mi=∑j∇Wij⊗(rj-ri)ΔVj，其中，上标L表示拉普拉斯算子，δ取0.1．2 不同黏性处理方法的动量方程2.1　人工黏性Monaghan [11]提出的人工黏性是SPH流体模拟中最广泛使用的一种方法，人工黏性项Πij具体表达式为Πij=-αc¯ijμij/ρ¯ij     (uijrij0),0     (uijrij≥0);μij=huijrij/(rij2+0.01h2), (4)式中：uij为粒子i与j相对速度；c¯ij为平均声速；ρ¯ij 为平均密度；α为控制速度耗散强度的系数，通常取值范围为0~0.5，α=0.01已经被证实是波浪水槽的最佳取值[12]，并且适用于波浪传播和沿海结构波浪荷载的研究[13]．在SPH表示法中，人工黏性项通常施加在压力项上，引入人工黏性的动量方程(2)可以写成duidt=-∑jmjpj+piρjρi∇iWij+1ρiΠij+g. (5)人工黏性不仅可以消耗高频能量和数值振荡引起的流场突变，还能在一定程度上阻止粒子相互接近时的非物理穿透，对于自由表面流动，它在数值上也起到了稳定数值的作用．2.2　层流黏性和SPS湍流黏性在黏性流SPH方法中，纳维-斯托克斯方程采用粒子近似法进行离散并基于弱可压缩假设进行求解，流体的物理黏性考虑层流黏性和SPS湍流模型，对边界层的精细流动进行模拟．动量方程中的层流黏滞应力可表示为[14]υ0∇2ui=∑jmj4υ0rij∇iWij(ρi+ρj)(rij2+η2)uij，(6)式中υ0为流体的运动黏度(通常水的运动黏度1×10-6 m2/s)．参量η=0.1h用于防止粒子过于靠近时可能发生的数值离散，SPH的离散形式为duidt=-∑jmjpj+piρjρi∇iWij+g+∑jmj4υ0rij∇iWij(ρi+ρj)(rij2+η2)uij. (7)SPS湍流模型动量守恒方程定义为dudt=-∇p/ρ+g+υ0∇2u+1ρ∇τ，(8)式中：层流项υ0∇2u按式(6)处理，ρ为流体密度；τ为SPS的应力．在Dalrymple和Rogers[5]的弱可压缩SPH方法中，使用涡流黏度假设对SPS应力张量进行建模，应力张量用爱因斯坦表示法表示，τ可以表示为τmn/ρ=νt2Smn-2kδmn/3-2CIΔl2δmnSmn2/3，式中：m和n分别为坐标轴x和y方向；νt=(CSΔl)22S为湍流涡黏性系数，其中，CS为Smagorinsky常数[15]，通常取0.12，Δl为粒子间距，S=2SmnSmn，Smn为SPS应变张量的一个元素；δmn为Kronecker符号函数；系数CI=0.006 6；k为SPS湍动能．运用Favre平均，式(7)可以写为duidt=-∑jmjpj+piρjρi∇Wij+g+∑jmj4υ0rij∇Wij(ρi+ρj)(rij2+η2)uij+∑jmjτmnjρj2+τmniρi2∇Wij.3 结果验证与分析开展Re=200和9 500时的圆柱绕流算例，分别采用人工黏性、层流黏性和SPS湍流模型对流体黏性进行处理，进行不同黏性处理结果的升阻力系数计算结果对比．阻力系数Cd和升力系数CL分别为Cd=2Fd/(ρU2D)；CL=2FL/(ρU2D)，式中：Fd和FL分别为圆柱受到的阻力和升力；U为来流速度；D为圆柱直径．3.1　数值计算模型直径为D的圆柱固定在无限流场中，均匀来流从左至右，左边界和上边界为流入边界，下边界和右边界为流出边界，圆柱表面为无滑移壁面边界．计算域为长方形，沿来流方向的长度20D，垂直于来流方向的宽度为12D．初始粒子间距Δl=D/100，实粒子总数为2.392 076×106，时间步长Δt=0.1Δl/U(即CFL为1)，运用粒子打包算法实现初始粒子均匀分布，如图1所示．10.13245/j.hust.24010401.F001图1近圆柱区域初始粒子分布3.2　Re=200时二维圆柱绕流升阻力系数对比以Wile等[16]开展的圆柱绕流实验和Liu等[17]、Ng等[18]、Rossi等[19]的计算结果为参考，比较了Re=200时圆柱所受阻力系数和升力系数，具体结果如表1所示．10.13245/j.hust.24010401.T001表1Re=200时圆柱绕流升阻力系数和相关文献的对比来源CdCLWile等[16]1.300 0Liu等[17]1.310 0±0.049 0±0.690Ng等[18]1.373 0±0.050 0±0.724Rossi等[19]1.354 0±0.050 0±0.680人工黏性模型1.764 0±0.476 0±1.712层流黏性模型1.327 0±0.055 8±0.805SPS湍流模型1.327 0±0.053 3±0.804层流黏性模型和SPS湍流模型计算的阻力系数均值一致而幅值略有差异，二者接近于Wile等[16]的实验值和Liu等[17]的计算值，而升力系数幅值略有增大；人工黏性模型升阻力计算值的均值和幅值均未取得良好的数值精度，系数α对该模型计算结果的影响较大，在数值模拟过程中往往须通过多次试算调整α取值，一定程度上限制了模型的适用性．图2为Re=200时采用不同黏性处理方法的阻力系数和升力系数时历曲线．可见：层流黏性模型和SPS湍流模型的升阻力曲线分别在tU/D=57，54后表现为正弦分布；相对于层流模型，SPS湍流模型计算结果更早趋于稳定，而人工黏性模型一直表现为不规则分布．10.13245/j.hust.24010401.F002图2Re = 200时不同黏性处理方法的圆柱绕流升阻力系数图3为tU/D=200时不同黏性处理方法的湍流涡黏性系数对比．层流黏性和SPS湍流模型的涡黏性系数整体分布形态相近，再现了近壁面区域尾流场涡街交替脱落引起的涡黏性系数变化，而人工黏性模型近壁面区域已变现为湍流特征，效果最差．10.13245/j.hust.24010401.F003图3Re=200时不同黏性处理方法瞬时湍流涡黏性系数3.3　Re=9 500时二维圆柱绕流升阻力系数对比为进一步验证不同黏性处理方法在较高雷诺数下的有效性，计算了Re=9 500时二维圆柱绕流升阻力系数，如图4所示．由图4可见：随着雷诺数的增大，层流和SPS湍流模型升阻力系数在湍流阶段也不再保持周期性变化．10.13245/j.hust.24010401.F004图4Re = 9 500时圆柱绕流升阻力系数进一步选取了tU/D=180~200时间段内3种模型升阻力进行对比，如图5所示．可见：在湍流阶段，3种模型阻力系数均值接近，由于边界层流动的不稳定性使得数值模拟结果存在不同程度的数值脉动；相对于层流计算结果，人工黏性和SPS湍流模型的升力系数接近．10.13245/j.hust.24010401.F005图5Re=9 500时不同黏性处理方法的圆柱绕流升阻力系数图6为tU/D=200时不同黏性处理方法下涡量对比，可见SPS湍流模型可有效捕捉到涡街脱落引起的二次涡，进一步验证了算法的有效性．10.13245/j.hust.24010401.F006图6不同黏性处理方法下瞬时涡量对比综上可知：SPS湍流模型在升阻力预报和涡量场模拟方面都有不错的表现，将其进一步应用于三维圆柱绕流有助于解决后续更多相关工程问题．4 穿透自由液面的三维圆柱绕流4.1　计算域及边界条件设置计算模型的笛卡尔坐标系统定义，原点为静止水位的圆柱体中心，分别取顺流向、横流向和垂直方向为x，y和z轴方向．来流方向从左向右，与x轴的正方向一致．计算域如图7所示，圆柱直径D，圆柱高度为6D，水深为4D，计算域来流方向的长度10D，垂直于来流方向的宽度为6D．边界条件设置如下：左边界设置为流入边界，来流速度为U，右边界为流出边界，两侧面及底面设置为对称边界，圆柱表面采用无滑移边界条件．10.13245/j.hust.24010401.F007图7三维圆柱绕流计算域示意图4.2　升阻力系数对比表2为Re=27 000，弗劳德数Fr=0.8时的阻力系数均值及升力系数均方根与文献值的对比，可见：本文阻力系数均值与无限长圆柱实验值接近，受圆柱长度影响升力系数均方根与实验值存在一定的差异．10.13245/j.hust.24010401.T002表2升阻力系数均值及升力系数均方根与文献值的对比来源阻力系数均值升力系数均方根Inoue等[20]实验值1.2000.450Kawamura等[21]计算值0.970~1.1200.240~0.320本文计算值1.2970.6394.3　自由液面预报结果对比当流体绕过圆柱时，自由液面会出现波动，在柱体前驻点，流体速度为0而出现自由液面的抬升．对于理想流体，根据伯努利方程可知其自由液面的抬升高度约为Fr2D/2．表3为自由液面抬升高度与理论值的对比，计算值与理论值非常接近．10.13245/j.hust.24010401.T003表3自由液面抬升高度与理论值的对比Fr计算值理论值0.80.350D0.32D1.00.527D0.50D1.20.700D0.72D图8为Fr=0.8时不同水深(z/D=-0.5、z/D=-1.5)处涡量分布云图．受自由液面变形影响，近自由液面区域尾流场表现为大量小尺度涡，随水深的增加尾流场逐渐出现涡街脱落现象，可见本文方法能够有效捕捉受自由液面变形引起的涡量变化，对自由液面的模拟是准确有效的．10.13245/j.hust.24010401.F008图8不同水深处涡量分布云图5 结论a．与传统的SPH方法相比，SPS-δ-SPH方法在圆柱绕流升阻力预报和涡量场模拟中表现出良好的模拟精度．算例表明该方法能够真实地体现出流体的物理黏性，适用于二维钝体绕流问题SPH方法模拟．b．采用SPS-δ-SPH模拟穿透自由液面的三维圆柱绕流模拟，自由液面的抬升预报结果与理论值接近，表明SPS-SPH方法是一种可靠且有效的预报自由表面黏性流动的研究方法．