• 全国中文核心期刊
  • 中国科技核心期刊
  • 美国工程索引(EI)收录期刊

留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

水平圆柱强迫振荡水动力特性的数值模拟

毛鸿飞 袁剑平 庞建华 贾宝柱

毛鸿飞, 袁剑平, 庞建华, 贾宝柱. 水平圆柱强迫振荡水动力特性的数值模拟[J]. 水科学进展, 2020, 31(2): 278-286. doi: 10.14042/j.cnki.32.1309.2020.02.014
引用本文: 毛鸿飞, 袁剑平, 庞建华, 贾宝柱. 水平圆柱强迫振荡水动力特性的数值模拟[J]. 水科学进展, 2020, 31(2): 278-286. doi: 10.14042/j.cnki.32.1309.2020.02.014
MAO Hongfei, YUAN Jianping, PANG Jianhua, JIA Baozhu. Numerical investigation of hydrodynamic characteristics on horizontal circular cylinder induced by a forced oscillatory[J]. Advances in Water Science, 2020, 31(2): 278-286. doi: 10.14042/j.cnki.32.1309.2020.02.014
Citation: MAO Hongfei, YUAN Jianping, PANG Jianhua, JIA Baozhu. Numerical investigation of hydrodynamic characteristics on horizontal circular cylinder induced by a forced oscillatory[J]. Advances in Water Science, 2020, 31(2): 278-286. doi: 10.14042/j.cnki.32.1309.2020.02.014

水平圆柱强迫振荡水动力特性的数值模拟

doi: 10.14042/j.cnki.32.1309.2020.02.014
基金项目: 

国家自然科学基金资助项目 51479017

广东海洋大学科研启动费资助项目 51479017

详细信息
    作者简介:

    毛鸿飞(1985—) , 男 , 辽宁铁岭人 , 讲师 , 博士 , 主要从事波浪与结构物相互作用方面研究。E-mail : maohongfei-gdou@qq.com

    通讯作者:

    袁剑平 , E-mail:yjp_103@163.com

  • 中图分类号: TV139.2+6

Numerical investigation of hydrodynamic characteristics on horizontal circular cylinder induced by a forced oscillatory

Funds: 

The study is financially supported by the National Natural Science Foundation of China 51479017

The study is financially supported by the National Natural Science Foundation of China 51479017

图(10) / 表 (1)
计量
  • 文章访问数:  25
  • HTML全文浏览量:  8
  • PDF下载量:  5
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-08-19
  • 网络出版日期:  2020-02-14
  • 刊出日期:  2020-03-01

水平圆柱强迫振荡水动力特性的数值模拟

doi: 10.14042/j.cnki.32.1309.2020.02.014
    基金项目:

    国家自然科学基金资助项目 51479017

    广东海洋大学科研启动费资助项目 51479017

    作者简介:

    毛鸿飞(1985—) , 男 , 辽宁铁岭人 , 讲师 , 博士 , 主要从事波浪与结构物相互作用方面研究。E-mail : maohongfei-gdou@qq.com

    通讯作者: 袁剑平 , E-mail:yjp_103@163.com
  • 中图分类号: TV139.2+6

摘要: 为了研究处于自由面以下完全淹没状态的水平圆柱在强迫振荡运动时的水动力特性, 采用基于黏性流理论建立的二维两相流数值波浪水槽模型, 对不同液相黏性条件下强迫振荡水平圆柱的受力进行计算, 并对压力、黏性切力和圆柱运动之间的相位关系特征进行对比和分析, 进而结合流场分析解释黏性影响机理。结果表明:黏性切力和涡旋压力对流体作用力的贡献差别是导致不同流体黏性下流体作用力结果差异的主要原因; 涡旋运动相对圆柱振荡运动的滞后性受流体黏性影响显著, 导致不同流体黏性下压力之间有相位差; 流体水质点相对于圆柱的滞后运动在大黏性流体中更为显著, 这导致了其黏性切力的相位超前现象。

English Abstract

毛鸿飞, 袁剑平, 庞建华, 贾宝柱. 水平圆柱强迫振荡水动力特性的数值模拟[J]. 水科学进展, 2020, 31(2): 278-286. doi: 10.14042/j.cnki.32.1309.2020.02.014
引用本文: 毛鸿飞, 袁剑平, 庞建华, 贾宝柱. 水平圆柱强迫振荡水动力特性的数值模拟[J]. 水科学进展, 2020, 31(2): 278-286. doi: 10.14042/j.cnki.32.1309.2020.02.014
MAO Hongfei, YUAN Jianping, PANG Jianhua, JIA Baozhu. Numerical investigation of hydrodynamic characteristics on horizontal circular cylinder induced by a forced oscillatory[J]. Advances in Water Science, 2020, 31(2): 278-286. doi: 10.14042/j.cnki.32.1309.2020.02.014
Citation: MAO Hongfei, YUAN Jianping, PANG Jianhua, JIA Baozhu. Numerical investigation of hydrodynamic characteristics on horizontal circular cylinder induced by a forced oscillatory[J]. Advances in Water Science, 2020, 31(2): 278-286. doi: 10.14042/j.cnki.32.1309.2020.02.014
  • 在水动力学研究中, 准确地预测结构物上的流体作用力是指导工程设计和安全生产的重要问题。作为一个经典的水动力学问题, 淹没水平圆柱强迫振荡问题受到广泛关注, 早期开展的相关研究多基于势流理论。Wu[1]采用多级子展开方法对淹没水平圆柱大幅强迫振荡受力开展了理论解析研究, 给出了圆柱上流体作用力各倍频分量。Tyvand和Miloh[2]采用时间级数展开方法对水平圆柱初始冲击运动问题开展了解析研究, 求得圆柱上流体作用力和辐射波面。Wu和Eatock Taylor[3]采用边界元和有限元结合方法对淹没水平圆柱强迫振荡的非线性辐射波进行了计算。Liu等[4]采用边界积分方法对淹没水平圆柱强迫振荡的非线性辐射波浪开展了解析研究。Kent和Choi[5]采用边界条件分解显式求解方法, 对垂向强迫振荡的淹没水平圆柱受力进行了数值计算, 其一倍频力与Wu[1]的解析解吻合较好, 而高倍频力结果有一定差别。Guerber等[6-7]应用基于边界元方法建立的二维数值波浪水槽模型对淹没水平圆柱强迫振荡受力进行了数值计算研究, 所得各倍频力与Wu[1]的解析解对比表明, 二倍频力受到自由面条件非线性的影响。Teng等[8]基于黏性流理论建立了二维数值波浪水槽, 并给出了非线性波浪作用下的圆柱所受各倍频力的数值结果, 其与势流结果相比差别明显, 并通过流场分析解释了原因。Teng等[9]应用黏性流水槽模型, 考虑一种液相黏性, 对强迫振荡水平圆柱的受力进行了计算, 并与势流结果进行对比, 发现各倍频力的幅值有差别, 并通过流场分析解释了这一现象的成因, 即大振幅情况下产生的涡旋对圆柱所受压力幅值的影响。

    由于建立在流体无黏和流动无旋的基本假定下, 势流理论无法对黏性影响比较显著的物理问题给出准确预测。上述基于势流理论对淹没强迫振荡水平圆柱受力问题的研究没有考虑流体黏性对流体作用力的影响, 研究结果的准确性需要讨论。同时, 为了进一步说明流体黏性对这一黏性效应较显著问题的影响, 并深层次的揭示黏性影响机理, 需要基于黏性流理论、考虑多种流体黏性、针对淹没水平圆柱强迫振荡下的水动力特性进行预测和研究。借鉴Teng等[8-10]、Jiang等[11]、毛鸿飞和陈洪洲[12]分别对波浪作用于淹没圆柱结构物、液舱晃荡问题和水面附近圆柱受波浪冲击等问题的研究手段, 本文基于黏性流理论建立数值波浪水槽模型, 以Navier-Stokes(N-S)方程为控制方程描述流体运动;采用有限体积方法进行数值离散;采用Volume of Fluid(VOF)方法进行自由面捕捉;为了进行长时间序列的计算, 结合了基于松弛方法的造波和消波功能。为验证数值模型对结构物上流体作用力计算的准确性, 计算淹没水平圆柱在波浪作用下的受力, 并与前人的实验数据进行对比验证。应用黏性流数值模型对强迫振荡水平圆柱上的流体作用力进行数值计算, 考察不同流体黏性下圆柱所受流体作用力及其与圆柱运动之间的关系特征, 并与势流结果对比, 探明不同流体黏性下结果差别。通过分析压力、黏性切力和圆柱运动的相位关系以及圆柱周围局部黏性流流场特征, 解释不同流体黏性下数值结果之间差别形成的原因, 揭示流体黏性影响机理。

    • 控制方程为N-S方程, 考虑不可压缩黏性流体流动问题, 连续性方程和动量守恒方程可以张量形式表示为:

      $$ \frac{{\partial \rho {u_i}}}{{\partial {x_i}}} = 0 $$ (1)
      $$ \frac{{\partial \rho {u_i}}}{{\partial t}} + \frac{{\partial \rho {u_i}({u_i} - u_j^m)}}{{\partial {x_j}}} = - \frac{{\partial p}}{{\partial {x_i}}} + \mu \frac{\partial }{{\partial {x_j}}}(\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}) + \rho g $$ (2)

      式中: uiuj分别为流体质点速度在ij方向上的分量;ujm为网格运动速度;t为时间;p为流体压强;ρ为流体的密度;μ为流体动力学黏性系数;g为重力加速度。

    • 数值模型采用VOF方法对自由表面进行捕捉, α为计算单元体积分数, 表达为

      $$ \alpha = \left\{ {\begin{array}{*{20}{c}} 0&气相\\ {0 \sim 1}&自由面\\ 1&液相 \end{array}} \right. $$ (3)

      进而, 单元内流体密度和黏性系数为:

      $$ \rho {\rm{ = }}\alpha {\rho _g} + (1 - \alpha ){\rho _1} $$ (4)
      $$ \mu = \alpha \mu + (1 - \alpha ){\mu _1} $$ (5)

      式中:下标g和l分别表示气液两相。

      界面方程可表达为

      $$ \frac{{\partial \alpha }}{{\partial t}} + \frac{{\partial ({u_i} - u_j^m)\alpha }}{{\partial {x_i}}} + \frac{{\partial (u_i^r\alpha (1 - \alpha ))}}{{\partial {x_i}}} = 0 $$ (6)

      式中: uir为相对压缩速度, (uirα(1-α))/xi为人工压缩项, 该项存在自由面区域内, 其引入是为了提高界面分辨率和数值计算稳定性[13]

    • 数值波浪水槽模型及其边界如图 1所示。波浪生成边界AC定义为速度入口条件, 即采用速度边界法造波, 速度按所需波浪理论给出, 压强为p*∂ n =0(p*=p-ρgh为动态压强, h为液相中某点的深度)。消波边界BD和底边界CD均定义为不可滑移条件, 初始速度和压强分别为ui=0, p*∂ n =0。气相边界AB定义为可自由进出条件, 初始速度和压强条件分别为∂ui/∂ n =0, p*=0。水槽中结构物边界定义为不可滑移条件, 初始速度和压强条件分别为ui=0, p*∂ n =0。

      图  1  数值水槽边界和松弛区示意

      Figure 1.  Sketch of boundaries and relaxation zone for the numerical wave tank

      采用松弛方法[14]结合速度边界法实现造波和消波功能。如图 1所示, 为协助生成波浪和吸收结构物反射波浪, 松弛区Ⅰ设置在速度入口端;为消除消波边界的波浪反射影响, 松弛区Ⅱ设置在速度出口端。在松弛区内, 根据解析形式对数值求解过程进行修正, 区内流体速度和体积分数分别为

      $$ \left\{ {\begin{array}{*{20}{c}} {{u_r} = {\delta _R}{u_c} + (1 - {\delta _R}){u_a}}\\ {{a_r} = {\delta _R}{a_c} + ((1 - {\delta _R}){\alpha _a}} \end{array}} \right. $$ (7)

      式中: u为速度矢量;下标r、c和a分别表示对应物理量的目标值、计算值和解析值(松弛区Ⅰ);在消波区(松弛区Ⅱ)内, u a=0, pa=ρgh, αa为静水条件下定义值;松弛函数定义为δRχR, 表达式为

      $$ {\delta _R} = 1 - \frac{{\exp (\chi _R^{3.5}) - 1}}{{\exp (1) - 1}}\;\;\;\;\;\;{\chi _R} \in (0.1) $$ (8)
    • 结构物所受流体作用力为压力F p和黏性切力F v的合力, 表达式为

      $$ F = {F_P} + {F_V} = - \smallint pnds + \smallint \tau qds $$ (9)

      式中: n为单位法向量;q为单位切向量;τ =μdvs/dn为剪切应力, vs为切向速度, dvvs/dn为法向速度梯度。

    • 根据Chaplin[15]对淹没水平圆柱在波浪作用下受力的物理实验研究, 验证本文数值模型对结构物受流体作用力计算的准确性, 验证计算设置如图 2所示。r=5.1 cm为圆柱半径;S=5r=0.255 m为圆柱淹没深度, 即圆柱轴心到静水面距离;圆柱轴心距离造波区和消波区分别1.65L和4L(L为波长);造波区和消波区水平宽度分别为2.5L和3.0L。入射波条件为周期T=1.20 s, 波幅a=1.68~10.15 cm, 采用5阶Stokes波浪理论[16]造波。对应深水条件下的Kc数则可表示为

      图  2  波浪作用于淹没水平圆柱计算域示意

      Figure 2.  Sketch of computational domain for wave action on a submerged horizontal circular cylinder

      $$ {K_c} = \frac{{UT}}{{2r}} = \frac{{\pi A}}{r}{e^{ - kS}} $$ (10)

      式中:U为流体速度标量;k为波数。

      图 3为圆柱所受水平波浪力随Kc数变化数值结果与前人实验数据的对比。F/(ρr3ω2Kc)和a/r分别表示量纲一化的力和波幅。可见, 数值结果与实验数据吻合较好, 这验证了本文所用数值模型对结构物上流体作用力的计算具有较好的计算精度。需要说明的是, 量纲一力随着波幅增大呈现出先减小后增大的特征, 其原因是小波幅下, 圆柱所受升力与惯性力方向相反, 并随着波幅增大升力比例增加;大波幅下, 涡旋脱落对圆柱所受压力造成了影响, 升力的作用逐渐减小。

      图  3  波浪力随波幅的变化情况

      Figure 3.  Wave forces as functions of wave amplitude

    • 考虑振荡方向为水平方向, 模拟淹没水平圆柱强迫振荡。参照Wu[1]的势流解析条件, 数值计算参数设置如下:水深为d=4 m, 圆柱半径为r=0.1 m, 初始时刻圆柱位于水槽水平向中央, 轴心距离静水面S=3r=0.3 m, 水槽长度为6L, 水槽两端设有消波区, 长度均为2.5L, 如图 4所示。

      图  4  淹没水平圆柱强迫振荡计算示意

      Figure 4.  Sketch of computation for a submerged horizontal circular cylinder undergoing forced oscillation

      水平圆柱强迫振荡运动方程表示为

      $$ x(t) = Bsin(\omega t) $$ (11)

      式中:B为振幅;ω为振荡角频率。Re数则为

      $$ {\mathop{\rm Re}\nolimits} = \frac{{2Ur}}{v} = \frac{{2B\omega r}}{v} $$ (12)

      式中:υ=μ/ρ为流体运动学黏性系数。

      振荡频率以量纲一形式表示为kr=0.1和1.0;振幅以量纲一形式表示为B/r=0.20~1.75。为了探明流体黏性的影响, 运动学黏性系数分别为υ=1×10-6(m2/s)和υ=1×10-4(m2/s), 对应的计算参数如表 1所示。为便于表述, 将两组工况称为“大Re数”和“小Re数”工况, 对应真实流体可参考水和油品。

      表 1  数值计算参数

      Table 1.  Parameters for the numerical simulations

      序号 B/r υ=1×10-6(m2/s) υ=1×10-4(m2/s)
      Re (kr=0.1) Re (kr=1.0) Re (kr=0.1) Re (kr=1.0)
      1 0.20 6.26×103 1.98×104 6.26×10 1.98×102
      2 0.40 1.25×104 3.96×104 1.25×102 3.96×102
      3 0.60 1.88×104 5.94×104 1.88×102 5.94×102
      4 0.80 2.51×104 7.92×104 2.51×102 7.92×102
      5 1.00 3.13×104 9.90×104 3.13×102 9.90×102
      6 1.25 3.92×104 1.24×105 3.92×102 1.24×103
      7 1.50 4.70×104 1.49×105 4.70×102 1.49×103
      8 1.75 5.48×104 1.73×105 5.48×102 1.73×103

      数值计算中, 网格设置原则为一个波长内(根据辐射波要素)划分单元个数为100, 一个波高范围内划分单元数量为20。圆柱周围网格设置加密区, 加密半径设置为0.75r。加密区内沿圆柱法向设置等比例变宽单元(由内及外变大), 数量为24;沿圆柱轮廓设置等长度单元, 数量为120;初始时间步长为T/1 000(T为振荡周期), 并进行自适应调整, 计算总时长为40T

    • 对于水平圆柱所受流体作用力的水平分量即水平力, 本文黏性流数值结果与势流结果的对比如图 5所示。用于对比的势流结果分别为Wu[1]考虑线性自由面条件的势流结果和本文采用高阶边界元方法考虑非线性自由面条件(理论方法参见文献[17])所求得的势流结果。F/(ρBπr2ω2)为水平力的量纲一化形式。可见, 大Re数, 小振幅下, 黏性流与势流结果吻合较好;大振幅下, 随着振幅增大, 黏性流结果小于势流结果, 且二者差距逐渐增大。小Re数, 小振幅下, 黏性流结果大于势流结果;随着振幅增大, 黏性流结果先减小后增大。黏性流和势流结果之间明显差别说明流体黏性对流体作用力影响显著;两组黏性流结果不同特征表明, 流体黏性对流体作用力的影响机理不同;两组势流结果间的微小差别来源于自由面条件的不同。

      图  5  水平流体作用力随振幅变化情况

      Figure 5.  Horizontal hydrodynamic forces as functions of amplitude

      kr=0.1, B/r=0.40、1.25和1.75为例, 考察流体作用力中压力和黏性切力对总流体作用力的贡献, 对比黏性流和势流理论结果, 如图 6所示。可见, 小振幅下, 黏性流与势流下压力大小和相位均较接近;随着振幅增大, 黏性流下的压力小于势流结果;当振幅增大到一定程度, 小Re数下量纲一压力增大, 甚至又大于势流结果。需要强调的是, 大振幅下, 黏性流下的压力相比势流结果有相位滞后现象, 且小Re数下更为滞后。相比大Re数结果, 小Re数下的黏性切力相位明显超前。此外, 小Re数工况下, 黏性切力对流体作用力贡献较大, 导致了小振幅下的流体作用力大于势流结果。

      图  6  压力和黏性切力时间历程(kr=0.1)

      Figure 6.  Time-series of pressure forces and viscous shear forces (kr=0.1)

      考察当圆柱上水平流体作用力为极值时圆柱轴心的水平位置情况, 对比黏性流和势流结果, 以kr= 0.1为例, 如图 7所示。量纲一形式x/B表示水平位置, -Fmax和+Fmax分别代表力为负极值和正极值的瞬时时刻。可见, 黏性流下, 随着振幅增大, 在水平力为极值时刻, 圆柱轴心位置逐渐向振荡中心(x/B=0)偏移, 且该现象小Re数下更为显著。另外, 较小振幅下(B/r=0.20), 大Re数下圆柱轴心位置与其他结果有差别, 这是由于较大的黏性切力影响所致。

      图  7  水平力达到极值时圆柱轴心位置(kr=0.1)

      Figure 7.  Axis positions of the circular cylinder when the horizontal forces reach the extreme values (kr=0.1)

    • 对两种流体黏性下, 圆柱周围的涡量场(V=D(∂v/∂x-∂u/∂z)/U, 其中uv为速度在水平和垂向分量, D为圆形域直径, U为水质点运动速度)和动态压强场(p*)进行分析, 以kr=0.1为例, B/r=1.00、1.25和1.75工况, 水平压力为正极值时的流场为例, 如图 8所示。Fi为流体相对加速度方向, 即惯性力方向。可见, 较大振幅下, 圆柱附近产生涡旋, 并在压力为极值时涡旋出现在高压区一侧。由于涡旋压强小于外场压强, 使得此处圆柱表面压强减小, 即正压压强减小, 进而导致压力幅值小于势流结果。特殊地, 小ReB/r=1.75工况下, 涡旋压力在流体作用力占主导, 涡旋起增大压力的作用, 这导致黏性流结果大于势流结果。

      图  8  水平压力为正极值时圆柱周围涡量场和压强场(kr=0.1)

      Figure 8.  Vorticity and pressure fields around the cylinder when the horizontal pressure forces reach the positive extreme values (kr=0.1)

      kr=0.1、B/r=1.25工况为例, 一个周期内几个瞬时时刻的黏性涡量场如图 9所示。ucac分别表示圆柱运动速度和加速度方向, +B代表水平正向振荡端点位置。可见, 涡旋的运动始终滞后于圆柱运动(加速度相位), 致使黏性流下压力相位滞后于势流结果, 并且在小Re数下, 涡旋运动之后现象更明显, 对应压力相位的滞后更加显著。

      图  9  一个周期内几个瞬时时刻圆柱周围的涡量场(kr=0.1, B/r=1.25)

      Figure 9.  Vorticity and pressure fields around the cylinder when the horizontal pressure forces reach the positive extreme values (kr=0.1, B/r=1.25)

      kr=0.1、B/r=1.25工况为例, 考察黏性切力为水平正极值时圆柱周围流线分布, 对比两种流体黏性结果, 如图 10所示。C表示振荡中心位置。可见, 黏性切力为正极值时, 大Re数下圆柱轴心与振荡中心更加接近, 此时圆柱速度接近极值;而小Re数下, 由于流体运动滞后性较强, 圆柱向振荡中心运动过程中, 流体运动方向与圆柱运动方向相反, 流体与圆柱的相对速度达到极值。因此, 相比大Re数结果, 小Re数下的黏性切力有相位超前现象。

      图  10  黏性切力为正极值时圆柱周围流线分布(kr=0.1, B/r=1.25)

      Figure 10.  Streamline distribution around the circular cylinder when the horizontal viscous shear forces reach the positive extreme values (kr=0.1, B/r=1.25)

    • (1) 黏性切力对流体作用力的贡献以及涡旋压力对流体作用力的影响程度不同, 导致了不同流体黏性下圆柱受力的差别。

      (2) 涡旋运动的滞后性致使黏性流下的压力相位滞后于势流结果, 且该现象在流体黏性较大时更加显著。

      (3) 流体黏性较大时, 相比圆柱运动, 流动运动滞后性更强, 二者相对运动速度极值发生在周期内更早的瞬时时刻, 导致了对应黏性切力的相位超前现象。

参考文献 (17)

目录

    /

    返回文章
    返回