计算流体力学(computational fluid dynamics，CFD)是船舶与海洋工程领域的研究热点之一，如水下爆炸、波浪冲击、流固耦合等，研究手段分为网格法和粒子法两大类．网格法技术发展非常成熟，对于常见水动力学问题都可以提供非常好的计算结果，但无法很好地解决流体破碎、大变形等问题．粒子法采用粒子表述物理模型，不存在网格的约束，对于模拟自由液面、动边界及大变形问题具有先天优势[1]．为了解决经典连续介质力学在模拟裂纹问题上的困难，Silling等[2]提出用非局部的近场动力学理论(Peridynamic theory，PD)处理不连续问题，如今PD方法已经在固体力学、断裂力学等领域发挥重要作用．Madenci等[3]在近场动力学理论中的粒子相互作用概念及泰勒展开思想的基础上提出近场动力学微分算子(peridynamic differential operator，PDDO)，推导出任意阶微分方程与积分运算的等价关系；Bazazzadeh等[4]在势流理论基础上，采用PDDO模拟不可压缩理想流体的液舱晃荡问题；Cong等[5]采用PDDO模拟欧拉框架下不可压缩黏性流体的计算模型，研究不同雷诺数下的圆柱绕流问题，PDDO结果与ANSYS结果二者之间的误差较小；Chang等[6]采用近场动力学微分算子建立欧拉框架下的内部流动计算模型，引入压力校正及粒子细化等策略模拟方腔驱动及低雷诺数圆柱绕流，模拟效果较好；Gao等[7-8]采用PDDO建立了低雷诺数层流的非局部更新拉格朗日计算模型，分别对Poiseuille流、Couette流、泰勒-格林涡及溃坝问题进行研究，进一步将PDDO应用于流固耦合研究中，将计算结果与已有数据进行对比，验证了该计算模型的有效性．为验证近场动力学微分算子在计算流体力学领域的适用性，采用Fortran语言编制近场动力学微分算子数值计算程序，建立近场动力学微分算子的二维黏性流体计算模型，借助二维方腔驱动、溃坝及圆柱绕流三个经典流体算例对比不同时刻流体粒子的位置分布、速度云图及速度分布曲线，验证PDDO流体计算模型对简单流动、自由表面流动及绕流问题的可行性．1 基本理论1.1　PDDO理论在近场动力学微分算子中，对于一个二维标量函数f(x)=f(x,y)，其二阶泰勒展开形式为      f(x+ξ)=∑n1=0N∑n2=0N-n11n1!n2!ξ1n1ξ2n2∂n1+n2f(x)∂xn1∂yn2+R(N,x)，(1)式中：N=2为泰勒展开阶数；n1和n2分别为对x和y的求导阶数；ξ为粒子的相对位置，ξ1=x'-x，ξ2=y'-y；R(⋅)为泰勒展开余项．假设余量的贡献可忽略不计，将式(1)每一项乘以近场动力学函数gNp1p2(ξ)，在点x族内进行积分，并调用近场动力学函数的正交性质，1n1!n2!∫ξ1,ξ2ξ1p1ξ2p2gNp1p2(ξ)dξ1dξ2=δn1p1δn2p2，(2)可得到任意阶偏导数的近场动力学非局部表达式，∂p1+p2f(x)∂xp1∂yp2=∫H(x)(f(x+ξ)-f(x))gNp1p2(ξ)dV', (3)式中：δn1p1和δn2p2为Kronecker符号；V'为粒子体积；p1和p2分别为对x和y的求导阶数．将近场动力学函数重新构造为多项式形式gNp1p2(ξ)=∑q1=0N∑q2=0N-q1aq1q2p1p2wq1q2(ξ)ξ1q1ξ2q2，(4)式中：ξ为粒子间距；wq1q2(ξ)为粒子相互作用的权函数[3]，一般为wq1q2(ξ)=e-(2ξ/δ)2；q1，q2为对x，y的求解导数；δ为领域半径．将式(4)代入式(3)，可得∑q1=0N∑q2=0N-q1A(n1n2)(q1q2)aq1q2p1p2=bn1n2p1p2；(5)A(n1n2)(q1q2)=∫V'wq1q2(ξ)ξ1n1+q1ξ2n2+q2dV'；(6)bn1n2p1p2=n1!n2!δn1p1δn2p2，(7)式中：A(n1n2)(q1q2)为形状函数；aq1q2p1p2为求解的系数矩阵；bn1n2p1p2为已知的系数矩阵．1.2　控制方程使用近场动力学微分算子对流体连续性方程和动量方程进行离散求解．流体假设为弱可压缩，补充额外的状态方程求解压强与密度之间的关系，拉格朗日格式∂ρi∂t=-ρi∫Vj[(uj'-ui')g210+(vj'-vi')g201]dVj，(8)ρi∂ui∂t=ρig-∫Vj[(pj-pi)g210]dVj+μ∫Vj[(uj-ui)(g220+g202)]dVj; (9)欧拉格式    ∂ρi∂t=-ρi∫Vj[(uj'-ui')g210+(vj'-vi')g201]dVj-∫Vj(ρj-ρi)(ui'g210+vi'g201)dVj, (10)ρi∂ui∂t=ρig-∫Vj[(pj-pi)g210]dVj+μ∫Vj[(uj-ui)(g220+g202)]dVj-∫Vj(uj-ui)(ui'g210+vi'g201)dVj; (11)状态方程pi=c02(ρi-ρ0)，(12)式中：粒子i与j为相邻粒子；u=[ui，uj]为速度矢量，u'和v'为u和v在x和y方向的速度分量；p为压强；t为时间；g为重力加速度；μ为动力黏性系数；ρ0为初始时刻的密度；c0为人工声速，通常为最大粒子速度的10倍及以上；g210为近场动力学函数，下标2表示泰勒展开阶数，上标10表示对于x的一阶导数和y的零阶倒数．采用式(12)形式的状态方程可以获得更稳定的压强场[1]．1.3　边界条件模拟流体运动时，通常会存在自由面边界、固壁边界及周期性边界．在边界处粒子的作用域存在截断现象，导致计算误差，影响计算精度，因此须对该部分的粒子进行特殊处理，见图1，图中：蓝色表示流体；绿色为固壁边界；红色为流入边界；黄色为流出边界；v0为初始速度；v为流体粒子的速度；p为流体粒子的压强；ρ为流体粒子的密度；下标in表示入口边界粒子，out表示出口边界粒子，f表示流体．10.13245/j.hust.240470.F001图1流入流出边界条件示意图1.3.1　固壁边界对于固壁边界(绿色粒子)，主要有排斥力法[9]和虚拟粒子法[10]．在PDDO模型中，采用虚拟粒子法和排斥力法结合的方法具有更好的模拟效果，因此对于所模拟的流体运动都会在固定虚粒子的基础上，额外布置作用力较小的排斥力F作用，保证严格的粒子不可穿透．边界宽度根据粒子作用域进行确定，固壁粒子的位置保持不变，速度、压强和加速度根据周围流体粒子的物理量进行插值计算，并参与流体控制方程的计算，pw=∑pfwq1q2(ξ)+∑ρf(aw-g)ξwq1q2(ξ)∑wq1q2(ξ)，式中：pw和aw为固壁粒子压强和加速度；pf和ρf为流体粒子压强和密度；1.3.2　流入流出边界对于圆柱绕流问题，通常须在计算域左右两侧分别设置流入流出边界．在Tafuni等[11]提出的流入流出边界算法基础上实现流入流出边界．在流体计算域左右两侧分别设置固定虚粒子(红色粒子和黄色粒子)，物理量根据内部流体粒子进行镜像插值获得，根据计算的模型不同设置不同的进出口边界．1.4　移动最小二乘法在无网格方法中移动最小二乘法(MLS)是一种广泛使用的近似算法，通常被用于边界条件、光滑压力场及降噪等作用[12]．提出的计算模型中压力振荡现象较为明显，严重影响到数值计算的精度，因此当求解流体控制方程时，每30个时间步采用MLS算法平滑流场的密度场、压力场及速度场，有效减少数值震荡带来的误差，但过于频繁地应用MLS算法，会导致计算中无法及时捕捉到流体细节，特别是对于圆柱绕流问题尾部涡的模拟，应当适当放大时间间隔更有效．对于方腔驱动和溃坝问题，每30个时间步采用一次MLS算法；对于圆柱绕流问题，每150个时间步采用一次MLS算法，WjMLS=wq1q2(ξ)(β0+β1Δx+β2Δy)；β=[β0,β1,β2]T=A-1[1,0,0]T；A=∑j=1Niwq1q2(ξ)1ΔxΔyΔxΔx2ΔxΔyΔyΔxΔyΔy2；φi'=∑j=1NiφiWjMLS/∑j=1NiWjMLS，式中：WjMLS为光滑系数；Δx=xi-xj；Δy=yi-yj；φi和φi'为光滑前后的物理量；A和β为系数矩阵．2 基础算例与验证为研究所建立的基于近场动力学微分算子流体模型的可行性，选取方腔驱动流、溃坝流动及二维圆柱绕流3种算例验证流体计算框架的有效性，并与已有试验数据和数值结果进行对比分析．方腔驱动和溃坝采用拉格朗日形式的控制方程，圆柱绕流采用欧拉形式的控制方程．2.1　数值方法误差研究PDDO是在泰勒展开思想上建立起来的计算手段，泰勒展开阶数是影响数值精度的一大因素，因此借助函数g=x2+2y2研究其对数值精度的影响．如表1所示，分别给出PDDO泰勒展开至一阶、二阶及三阶导时解析解与数值解的平均误差值．10.13245/j.hust.240470.T001表1各阶导数误差阶数g'xg'ygxx″gyy″11.748 2×10-21.166 6×10-224.784 0×10-164.623 7×10-164.696 7×10-154.081 9×10-1532.533 3×10-151.932 2×10-153.190 9×10-142.660 5×10-14随着展开阶数越多，误差值越小，数值解与解析解基本一致，但过多展开项反而不会进一步提高精度，如展开三阶导时，gxx″和gyy″误差会略微增加．本文算例只须求得二阶导数即可，因此选取PDDO函数展开至二阶导基础上研究，领域半径范围[3]2.0 dx≤δ≤4.0 dx，dx为粒子间距．2.2　方腔驱动对于方腔驱动，流体均匀分布在1 mm×1 mm的正方形区域中，上壁面以0.001 m/s的速度沿x方向移动，其余3个壁面静止不动，流体在该作用下将形成稳态．根据物理条件，流体密度为1 000 kg/m3，运动黏性系数为1×10-6 m2/s-1，总计算时间步3 000步，时间步长dt=5×10-5 s．在流体粒子周围预先布置固定虚粒子模拟固壁边界，虚粒子的压强、速度及密度等物理量按照1.3节所述进行更新．为减小压力震荡对于计算结果的影响，每隔30个时间步长使用MLS算法计算得到新的密度场和速度场，压强采用状态方程进行计算；为减少粒子的不均匀性对计算结果的影响，采用粒子位移修正算法(PST)[13]对每一时刻流场位置以及物理量进行修正．为检验PDDO流体计算模型的精度及收敛性，分别设置粒子数分别为40×40，60×60，100×100，粒子间距分别为0.025，0.016 7，0.01 mm，计算结果与SPH[14]结果进行对比．图2给出了同一时刻的水平速度vx和垂直速度vy分布曲线．图中PDDO结果与SPH结果[14]的趋势及数值基本符合，相同位置速度的相对误差在10%左右，随着粒子数的增加，二者的误差逐渐减小，认为PDDO计算模型可行．10.13245/j.hust.240470.F002图2水平速度vx和垂直速度vy分布曲线1—SPH结果；2—dx=0.025 mm；3—dx=0.0167 mm；4—dx=0.01 mm．2.3　溃坝问题溃坝问题作为一种经典的自由表面流动，存在着复杂的大变形问题，对于无网格方法来说，不存在网格的限制，对于解决该类问题有着较大的优势．该问题一方面是要解决固壁边界问题，为保证流体粒子和固体粒子之间不发生穿透现象，采用1.3节的固定虚粒子方法实现壁面条件；另一方面对于自由面边界，一般只须将靠近自由面的流体粒子的压强设置为0 Pa即可满足运动条件[1]．采用较为简单的密度判别法[15]，即粒子密度与流体密度比值小于一特征值αi(ρi/ρ0≤αi)后，认为该粒子为自由面粒子．进行模拟计算时，在流体与右壁面接触前，αi可选取较小数值，能够更为准确地捕捉到自由表面粒子，在接触后适当提高αi以提高系统的稳定性．对比文献[16]数值模拟及试验结果进行验证．初始时刻流体静止于左壁面，底部宽度0.146 m，高度0.292 m，初始压强根据静水压强进行赋值，采用式(12)状态方程更新每一时间步的压强．流体初始密度为1 000 kg/m3，动力黏性系数为0.5 kg/(m⋅s)，g=9.8 m/s2，c0=15.0Vmax，整个流体模拟的区域布置60×30个粒子，dx=0.004 5，时间步长dt=0.05 ms，总模拟时间0.5 s，每隔30个步采用MLS算法光滑流体速度、密度及压强．图3给出4个时刻下流体粒子的运动形式及垂向速度分布云图．t=0.15 s时流体在重力的驱使下开始倒塌向右方流动，流动状态较为稳定；t=0.275 s时流体与右壁面发生接触，流体开始产生回流现象并且流体开始沿壁面上爬；t=0.4 s时流体沿着右壁面攀爬直至速度变为0 m/s开始下落，部分区域出现破碎现象．10.13245/j.hust.240470.F003图34个时刻下流体粒子的运动形式及垂向速度分布云图(色标单位：m/s)为进一步验证，选择流体底部最右端粒子为观测点，记录数值模拟过程中该粒子到右壁面的距离随时间的变化，绘制如图4所示时间-距离曲线，横、纵坐标皆进行无量纲化，d为观测点的x坐标，L为初始时刻溃坝流体的宽度，L=0.146 m．选取MPS结果[16]及模型试验数据[17]进行对比，PDDO与MPS结果和试验数据都较为符合，呈现近似线性趋势上升，验证PDDO流体计算模型对于自由表面流动模拟具有一定的可靠性．10.13245/j.hust.240470.F004图4时间-距离曲线1-MPS结果；2-PDDO结果；3-试验结果2.4　圆柱绕流圆柱绕流问题是船舶与海洋工程中较为经典的流固耦合问题，其几何结构简单，但流动现象较为复杂，特别是对于尾部涡的形成更是研究的重点，通常是研究涡激振动、绕流问题的基础．在之前研究的基础上，将PDDO流体模型应用于圆柱绕流的问题中．通过改变流场速度实现不同雷诺数，设置了Re=10和Re=100两种工况，研究稳态时流场的速度场及尾部涡的变化，验证PDDO方法的有效性．与FEM[5]计算结果进行对比，采用的计算模型如图5所示，流体初始密度为1.0 kg/m3，动力黏性系数为1.0 kg/(m∙s)，整个计算域长15 m，宽6 m，圆柱半径0.5 m，圆柱中心点(5，3) m，整个计算域总共布置300×120个粒子．为保证流场的稳定性和计算效率，总计算时间步为4×105．Re=10工况时间步长dt=0.1 ms，Re=100工况时间步长dt=0.05 ms．流入流出边界采用1.3.2节周期性边界进行设置，远方均匀来流以vx=10 m/s和vx=100 m/s速度进口，尾部采用压力出口，即进口处的压强和密度采用镜像插值求解，速度为初始值，出口处流体的速度保持不变，压力初始值设置为0 Pa．壁面和圆柱均采用非滑移壁面vwall=0，c0=15.0Vmax．为保证流场稳定性，每150时间步采用MLS算法计算流场的密度、压力和速度，并且在连续性方程中引入密度耗散项[1]．10.13245/j.hust.240470.F005图5计算模型图6为PDDO方法模拟Re=100时的流场分布云图．当Re=10时，流场的速度分布基本对称，圆柱两侧速度最大，圆柱尾部会出现规则的旋涡并且稳定地附着在圆柱表面，不会发生涡脱落现象；随着雷诺数增加，黏性力对流场的影响逐渐增加，圆柱尾部流场开始出现波动，伴随着涡的生成和脱落，最终形成图7中的反向旋转、排列有序的涡列，与FEM[5]结果对比，PDDO预测结果同样能够较好地体现流场的特征．10.13245/j.hust.240470.F006图6流场分布云图(色标单位：m/s)10.13245/j.hust.240470.F007图7流场涡量图(色标单位：s-1)图8给出了不同Re下流场达到稳定状态时的速度分布曲线．可以看出：PDDO方法与FEM方法同样具有较高的一致性，速度场的分布及变化趋势基本符合．当Re=100时，伴随着涡的不断生成和脱落，圆柱尾部流场速度曲线呈现波动式上升，进一步验证PDDO计算框架是有效的．10.13245/j.hust.240470.F008图8不同Re下流场达到稳定状态时的速度分布曲线1，2—Re=10；3，4—Re=100图9进一步给出Re=100时计算得到的圆柱升力系数(Cl)和阻力系数(CD)曲线，流动达到稳定之后，升力和阻力呈现规律周期性变化，并且平均升力系数在0.0，阻力系数峰值在1.5附近．10.13245/j.hust.240470.F009图9Re=100时计算得到的圆柱升力系数和阻力系数曲线3 结论a．采用PDDO求解流体连续性方程及动量方程实现流体运动时，各阶导数的选取基于泰勒展开思想，物理场数值振荡现象明显，采用移动最小二乘法(MLS)可以有效平缓流场中的振荡现象．b．基于流体弱可压缩假设，通过对比不同类型流体运动，与已有结果对比，建立的PDDO流体计算模型具有较好的可行性，对于计算流体力学领域的研究具有一定的意义．c．对于圆柱绕流问题，PDDO模型计算结果较好，能够准确模拟尾部涡的产生和脱落现象；而对于自由表面流动的模拟，建立的计算模型只是采用较为简单的密度判别，存在一定的误差，后续有待进一步改进．