海洋结构物在海洋环境中通常承受波-流的联合作用，而流的存在会显著改变结构物的受力，因此研究波-流-物体的相互作用对于海洋结构物的设计具有重要意义．基于时域线性理论，文献[1]模拟了波-流与三维单圆柱的相互作用问题，计算表明流的存在会对自由表面波和水动力产生显著影响；文献[2]计算了线性波-流与二维漂浮体作用的问题；文献[3]采用高阶边界元法计算了波-流作用在四个垂直圆柱上的波浪爬高和水动力，计算结果也表明流的存在对结构物上的载荷影响较大；文献[4-5]对三维单体结构物的绕射及运动问题也做了研究．关于线性波流体的研究文献还有很多，在此不多赘述．二阶波-流和结构物作用问题也有不少研究工作[6-8]，但大多是单体结构物．基于势流全非线性理论的波-流和单体结构物相互作用也有一些工作[9-11]，与线性或二阶理论相比，基于势流非线性理论的数值计算需要更多的计算资源和时间．关于波-流和多体结构物相互作用的研究，文献[12]基于时域二阶理论，采用高阶有限元法对稳态流作用下两个受迫振动的矩形截面柱体进行了计算，文献[13]基于时域二阶理论对波-流与四个座底圆柱的相互作用进行了模拟，结果表明流速对波浪爬高和水动力有着重要影响．通过文献分析发现，目前基于波-流与多体结构物相互作用的研究很有限，且对三维问题的研究均针对座底圆柱．为此本研究采用有限元法，基于时域二阶理论对处于均匀流中的双截断圆柱绕射问题进行研究，计算自由表面波浪爬高和水动力，详细分析流速、无因次波数及流向角对波浪和水动力的影响，同时分析了双圆柱对彼此的相互干扰现象．1 计算模型和控制方程如图1所示坐标系，两圆柱横剖面半径均为a，中心之间的距离为Lcy，圆柱1，2剖面中心分别位于(-Lcy/2，0)和(Lcy/2，0)．无扰动即水面静止时圆柱吃水d=a，水深h=3a．oxy平面和静水面S0重合，z轴垂直向上，坐标轴的原点位于双圆柱对称点上．Sf为瞬时自由表面，Sb为圆柱表面，单位方向矢量指向圆柱内部，Sbed为底部边界，位置在z=-h的平面上，Sc为人工截断边界．波浪沿x正轴方向传播；同时，一均匀流U沿与x轴成β方向流入．10.13245/j.hust.230174.F001图1坐标系根据势流理论，流场的速度势满足拉普拉斯方程∇2ϕ=0    (in  ∀), (1)式中∀为流体区域，in ∀为在整个流体计算域内．若流场内存在沿与x轴成夹角为β(流向角)且大小为U=(U0cos β,U0sin β)的稳态流，则流场总的速度势为Φ=U∙x+ϕ，(2)式中x=(x,y)．将速度势和波高在静水面S0上进行摄动展开到二阶，即ϕ=ϕ(0)+εϕ(1)+ε2ϕ(2)；(3)η=εη(1)+ε2η(2)，(4)式中：ε为摄动常数(一般可取为波陡)，其值为小量；ϕ(0)为由水流引起的扰动速度势，而由其引起的波高忽略不计．ϕ(r)，η(r)(r=1，2)分别为第k阶速度势和波高．由流引起的扰动势定解问题(r=0)为∇2ϕ(r)=0    (in ∀ (0))；(5)∂φ(r)/∂z=0    (on z=0)；(6)∂φ(r)/∂n=-U∙n    (on Sb)；(7)∂φ(r)/∂z=0    (on z=-h)，(8)式中：∀ (0)为由平均自由表面即静止水面S0及Sc，Sb和Sbed构成的流域；n为边界上的法向方向，n为其矢量形式，指向流域外部．对于绕射问题，一、二阶速度势ϕ(r)可进一步分解为ϕ(r)=ϕI(r)+ϕD(r) (r=1,2)，其中，下标I为入射；D为绕射．而一、二阶波也可做类似的分解，即η(r)=ηI(r)+ηD(r) (r=1,2)．一、二阶绕射势定解问题(r=1，2)为：∇2ϕD(r)=0    (in ∀  (0))；(9)∂ϕD(r)∂z-ηD(r)∂t=f1(r)    (on S0)；(10)∂ϕD(r)∂t+gηD(r)=f2(r)    (on S0)；(11)∂ϕD(r)∂n=-∂ϕI(r)∂n=-∇ϕI(r)∙n    (on Sb)；(12)∂ϕD(r)∂z=0  (r=1,2)    (on Sbed)，(13)式中g为重力加速度．对于上述各阶定解问题的详细描述及f1(r)，f2(r)(r=1，2)表达式参考文献[6，13]，而入射波为斯托克斯二阶波．此外，由于在时域内求解ϕD(r)(r=1,2)，还须提供初始条件．作用在物体上的水动力为F=∫Sbpnds, (14)式中p为压力，可由伯努利方程计算得出，p=-ρ(∂Φ/∂t+∇Φ∇Φ/2+gz)，(15)其中ρ为流体密度．将式(2)代入式(15)，再代入式(14)中并进行摄动展开得F=F(0)+F(1)+F(2)，(16)式中：F(0)为由流引起的力；F(1)，F(2)分别为一阶力和二阶力．2 有限元法和数值处理基于有限元法来进行数值模拟，采用6节点棱柱形单元(图2(a))．网格生成步骤如下：首先在静水面生成二维平面三角形网格(图2(b))，然后将该网格沿垂向向下拉伸，到达圆柱底部及以下空间须加上圆柱底部内生成的网格，这样就可生成图2(c)中的三维网格．便于显示，图2(c)中只给出了边界上的网格单元，圆柱表面的网格如图2(d)所示．10.13245/j.hust.230174.F002图2单元和网格示意图在生成网格后，便可采用有限元法进行计算．采用伽辽金加权余量法导出有限元方程，即KΦ(r)=F(r)    (r=0,1,2)，(17)式中：Φ(r)=[ϕ1(r),ϕ2(r),⋯,ϕN(r)]T为流场内节点上的k阶速度势向量，N为节点总数；K为有限元系数矩阵；F (r)为右端向量．式(17)可通过SSOR预条件共轭梯度法求解，从而获得流场内所有节点上的各阶速度势．流场内的各阶速度势ϕ(r) (r=0,1,2)获得后，通过差分法计算自由表面网格节点上的速度[13]．采用四阶Adams-Bashforth公式来更新每一时刻自由面上的一、二阶速度势和波位移．此外，为保证外传播条件，在流体区域外围设置人工阻尼区来吸收反射波[13]．3 计算结果与分析在波-流与圆柱相互作用问题研究中，关于座底单圆柱的情况较多，故本研究先模拟座底圆柱与波-流相互作用问题来验证本文数值方法的可靠性．取a=h，圆柱剖面中心位于坐标原点．弗劳德数Fr=U0/ga．图3为当无因次波数ka=1.0时计算所得自由表面与圆柱交线上点的一阶波幅值与文献[2]和[7]的计算结果对比(k为无流时一阶波波数)，图中θ为通过坐标原点与交线上的点的直线与x轴的夹角．从图中可以看出：本研究在Fr=-0.1，0.0，0.1下计算结果与文献[1]和[6]的结果均符合很好，特别是与文献[1]的结果更接近．10.13245/j.hust.230174.F003图3当ka=1.0时绕圆柱的一阶波幅值随θ/π的变化关系图4为当ka=0.5和kh=0.05时圆柱左端(x=-a，y=0)在三个不同弗劳德数下一阶和二阶合成波时间历程，图中：A为线性波波幅；T为无流时一阶波周期，图中还给出了文献[7]的计算结果．这里H=2A为线性入射波波高．从图中可以看出：本研究计算结果和文献[6]的结果符合良好，波谷和波峰值比文献[6]的结果略大，两者的相位无差异．10.13245/j.hust.230174.F004图4当ka=0.5和kh=0.05时圆柱左端一阶和二阶合成波时间历程通过图3和4中关于座底圆柱与波-流相互作用的计算结果对比表明了本文方法的可靠性．下面对两个截断圆柱与波-流相互作用问题进行数值模拟．双圆柱平面示意图如图5所示，流场区域限制在长L=8λ，宽B=6λ(λ为一阶波波长)的矩形区域内．计算中线性入射波波陡ε=0.05．为确保数值模型的有效性，对双圆柱模型作收敛性和网格敏感性分析，取无因次波数kca=1.0，Fr=0.1的情况进行验证，kc为有流时一阶波波数，考虑流存在时的波浪频率ωc2=kcgtanhkch，ωc=ω+kU0，其中ω为无流存在时波浪频率，U0为入射流的流速．选择两种不同的工况进行对比．工况一中自由表面三角形单元为2.531 3×104个，每个圆柱表面单元数为728个，流域内棱柱形单元总数为8.201 28×105，时间步长为∆t=T/200；工况二中自由表面三角形单元为2.837 1×104个，每个圆柱表面单元数为1 152个，流域内棱柱形单元总数为1.265 352×106，时间步长为∆t=T/400．由于二阶波对计算结果更加敏感，因此只对二阶波进行验证，结果如图6所示(其中A1，B1，A2和B2点位置见图5)，图中Tc为波流作用下一阶波的周期，Tc=2π/ωc．由图6可知：两种工况下的计算结果曲线几乎重合，故在工况一的网格大小和时间步长下的计算结果是收敛的．10.13245/j.hust.230174.F005图5双圆柱平面示意图10.13245/j.hust.230174.F006图6时间收敛性分析图7和8给出了当kca=1.0时Fr=-0.10，0.00，0.05和0.10下双圆柱左右两端点A1，B1，A2和B2处的一阶和二阶波时间历程．可以看出：四个点处的一阶波幅值基本上均随弗劳德数的增大而增大，而二阶波幅值只有B1，A2处的幅值随弗劳德数的增大而增大．因此，每个圆柱左、右端点上二阶波幅值并非都随弗劳德数增大而呈现规律性的变化，这与单个圆柱情况不同．单个圆柱情况下，左右两端一阶波幅值均随弗劳德数增大而增大，左端二阶波幅值随弗劳德数增大而减小，而右端二阶波幅值随弗劳德数增大而增大(单圆柱结果本研究未给出)．主要原因是二阶波受到另一个圆柱的干扰较严重，从而导致其随弗劳德数变化的规律被破坏．10.13245/j.hust.230174.F007图7当kca=1.0时一阶波时间历程10.13245/j.hust.230174.F008图8当kca=1.0时二阶波时间历程由于内侧两点B1和A2处的波浪受到圆柱相互干扰影响较大，图9给出了Fr=0.0，0.1下B1和A2两点处的一阶及一阶和二阶合成波对比．从图中可以看出：两点处的波峰和波谷表现出较明显的差异，特别是B1点，这说明两点处波的非线性比较明显，B1点处的波具有很强的非线性特性．从图9(a)和图9(c)可以看出：当Fr=0.1时，波谷的双峰更明显，显示更强的非线性特性．10.13245/j.hust.230174.F009图9当Fr=0.0，0.1时B1和A2两点处一阶及一阶和二阶合成波时间历程图10~12给出了圆柱上A1，A2，B1和B2四点处一阶波、二阶波幅值及一阶和二阶合成波峰值随kca的变化关系．从图10可知：在计算的无因次波数范围内，四个点处的一阶波幅值基本上均随弗劳德数增大而增大；B2点处的一阶波幅值随kca的增大而减小，其余三点处波幅值随kca增加呈现出上下波动即在局部出现极大值和极小值现象．而二阶波幅值(图11)，由于受到另一圆柱干扰较严重，其随Fr增大而并无变化规律．相比于一阶波幅值(图10)，一阶和二阶合成波峰值(图12)有了较明显的变化．随着弗劳德数的增加，A1点合成波峰值在整个无因次波数范围内仍然随弗劳德数增大而增大，而其他三点在kca较大时基本保持该特性，但在kca较小时随弗劳德数的增大其变化却并无规律，这主要是由二阶波幅值在小kca时其值较大且变化无规律造成(图11)．10.13245/j.hust.230174.F010图10一阶波幅值随kca变化10.13245/j.hust.230174.F011图11二阶波幅值随kca变化10.13245/j.hust.230174.F012图12一阶和二阶合成波峰值随kca变化图13给出了一阶波幅值和单个圆柱情况下的对比．由于B1受到圆柱干扰影响最大，这里只给出B1点处一阶波的对比．从图13可知：双圆柱下波的幅值和单圆柱相差很大，远大于A1点(A1结果由于篇幅限制未给出)．在单圆柱下一阶波幅值随kca基本呈线性减小，而双圆柱情况下则出现明显的波动．10.13245/j.hust.230174.F013图13B1点一阶波幅值对比图14~16给出了作用在圆柱1，2上的一阶力幅值|Fx(1)|和|Fz(1)|、二阶力幅值|Fx(2)|和|Fz(2)|及一阶和二阶合成力峰值Fx-max和Fz-max随kca变化关系．从图14可以看出：在整个无因次波数范围内，圆柱1和2的一阶水平和垂直力分量基本上均随弗劳德数增大而增大，这与一阶波类似．而二阶力(图15)，除圆柱2上的垂向力外(图15(d))，均只是在局部kca范围内才显示出随弗劳德数增大而增大的规律．而合成力峰值(图16)，除了个别kca外，也显示出与一阶力幅值相类似的规律，即随弗劳德数增大而增大．主要原因是由于在合成力成分中，一阶力起主导作用而二阶力相对较小．10.13245/j.hust.230174.F014图14一阶力幅值随kca变化10.13245/j.hust.230174.F015图15二阶力幅值随kca变化10.13245/j.hust.230174.F016图16一阶和二阶合成力峰值随kca变化上文主要计算了流向角β=0°和180°时双圆柱的波浪和水动力．下面对波-流存在一定夹角的情况进行计算讨论，夹角β设置为0°~180°范围每隔30°取一值．这里只给出水动力随流向角的变化关系．图17~19给出了各波-流夹角下圆柱1和2的水动力幅(峰)值随kca的变化关系．由图17可知：圆柱1上的水平和垂向力幅值在kca1.7范围内均随流向角增大而减小，而圆柱2上的水平和垂向力幅值在整个kca范围内均随流向角增大而减小．对二阶力(图18)而言，随着流向角的增大，可以有效减小圆柱1和2的二阶垂向力幅值，而对二阶水平力的影响相对较小．此外，圆柱1和2的一阶和二阶合成力峰值曲线(图19)的变化规律类似于一阶力(图17)．对于波浪而言也有类似的结论．10.13245/j.hust.230174.F017图17不同流向角下一阶力幅值随kca变化10.13245/j.hust.230174.F018图18不同流向角下二阶力幅值随kca变化10.13245/j.hust.230174.F019图19线性和二阶合成力峰值随kca变化图20给出了t=11.76Tc时刻三种弗劳德数下一阶波及一阶和二阶合成波的自由表面轮廓图．从图中可以看出弗劳德数对线性波及线性和二阶合成波的波形在空间上的影响．10.13245/j.hust.230174.F020图20当kca=1.0，t/Tc=11.76时自由表面轮廓图4 结论本研究基于时域有限元法对波-流与截断双圆柱的二阶相互作用问题，计算了双圆柱在波-流联合作用下的波浪和水动力．此外，通过考虑不同流向角，分析了流向角对波浪和水动力的影响．通过数值计算，可得出结论如下．a．双圆柱上内外侧四点处的一阶波基本随弗劳德数的增大而增大，而二阶波变化无明显规律．b．双圆柱内侧点受到干扰比外侧点严重．内侧点B1点处波浪的非线性特性要明显强于A1，A2和B2三个位置，且随着弗劳德数的增大其非线性特性增强．c．圆柱1，2的一阶水平和垂向力幅值基本上随着弗劳德数的增大而增大．d．波-流存在一定夹角能够较好地削弱波浪和水动力，且在多数情况下，结构物上的一阶波和一阶水动力随着流向角的增大而逐渐减小．