基于SPH方法的三维液体晃动数值模拟

第４２卷第１期　２０１０年２月　南京航空航天大学学报　Ｊｏｕｒｎａｌ　ｏｆ　Ｎａｎｊｉｎｇ　Ｕｎｉｖｅｒｓｉｔｙ　ｏｆ　Ａｅｒｏｎａｕｔｉｃｓ　８Ｌ　Ａｓｔｒｏｎａｕｔｉｃｓ　Ｖｏ１．４２　Ｎｏ．１　Ｆｅｂ．２０１０　

基于ＳＰＨ方法的三维液体晃动数值模拟　

刘　富　童明波　陈建平　

（南京航空航天大学先进设计技术国防重点学科实验窒，南京，２１００１６）　

摘要：采用ＳＰＨ方法对棱形液舱在外加激励作用下，不同充液比工况所对应的舱内液体晃动进行了三维数值模　拟，并将其与实验进行对比，两者吻合较好，同时成功地模拟出液体晃动产生的波浪翻卷和破碎，进而验证了　ＳＰＨ方法在求解三维液体晃动问题的准确性。分析结果表明：液舱的充液比及晃动周期对于液体的晃动特性及　舱壁的压力有重要的影响，在大充液比情况下，迟滞特性所表现的双峰特征消失。　关键词：光滑粒子流体动力学；数值模拟；晃动；液舱　中图分类号：Ｏ３５２　文献标识码：Ａ　文章编号：１００５—２６１５（２０１０）０１—０１２２—０５　

Ｎｕｍｅｒｉｃａｌ　Ｓｉｍｕｌａｔｉｏｎ　ｏｆ　Ｔｈｒｅｅ－Ｄｉｍｅｎｓｉｏｎａｌ　Ｌｉｑｕｉｄ　Ｓｌｏｓｈｉｎｇ　

Ｂａｓｅｄ　ｏｎ　ＳＰＨ　Ｍｅｔｈｏｄ　

Ｌｉｕ　Ｆｕ，Ｔｏｎｇ　Ｍｉｎｇｂｏ，Ｃｈｅｎ　Ｊｉａｎｐｉｎｇ　（Ｋｅｙ　Ｌａｂｏｒａｔｏｒｙ　ｏｆ　Ｆｕｎｄａｍｅｎｔａｌ　Ｓｃｉｅｎｃｅ　ｆｏｒ　Ｎａｔｉｏｎａｌ　Ｄｅｆｅｎｓｅ－Ａｄｖａｎｃｅｄ　Ｄｅｓｉｇｎ　Ｔｅｃｈｎｏｌｏｇｙ　ｏｆ　Ｆｌｉｇｈｔ　Ｖｅｈｉｃｌｅ，Ｎａｎｊｉｎｇ　Ｕｎｉｖｅｒｓｉｔｙ　ｏｆ　Ａｅｒｏｎａｕｔｉｃｓ＆Ａｓｔｒｏｎａｕｔｉｃｓ，Ｎａｎｊｉｎｇ，２１００１６，Ｃｈｉｎａ）　

Ａｂｓｔｒａｃｔ：Ｗｉｔｈ　ｔｈｅ　ｅｘｃｉｔａｔｉｏｎ，ｌｉｑｕｉｄ　ｓｌｏｓｈｉｎｇ　ｉｎ　ｔｈｅ　ｔｈｒｅｅ—ｄｉｍｅｎｓｉｏｎａｌ　ｐｒｉｓｍａｔｉｃ　ｔａｎｋ　ａｔ　ｄｉｆｆｅｒｅｎｔ　ｆｉｌｌ　ｒａ～　ｔｉｏｓ　ｏｆ　ｌｉｑｕｉｄ　ｃｏｎｄｉｔｉｏｎｓ　ａｒｅ　ｔｈｒｅｅ—ｄｉｍｅｎｓｉｏｎａｌｌｙ　ｓｉｍｕｌａｔｅｄ　ｂｙ　ｔｈｅ　ｓｍｏｏｔｈｅｄ　ｐａｒｔｉｃｌｅ　ｈｙｄｒｏｄｙｎａｍｉｃｓ　

（ＳＰＨ）ｍｅｔｈｏｄ．Ｔｈｅ　ｃａｌｃｕｌａｔｉｏｎ　ｒｅｓｕｌｔｓ　ａｒｅ　ｉｎ　ａｇｒｅｅｍｅｎｔ　ｗｉｔｈ　ｔｈｅ　ｅｘｐｅｒｉｍｅｎｔａｌ　ｒｅｓｕｌｔｓ．Ｔｈｅ　ｗｅｌｔｅｒ　ａｎｄ　

ｔｈｅ　ｆｒａｇｍｅｎｔａｔｉｏｎ　ｏｆ　ｔｈｅ　ｗａｖｅ　ａｒｅ　ｓｉｍｕｌａｔｅｄ．Ａｓ　ａ　ｒｅｓｕｌｔ，ｔｈｅ　ａｃｃｕｒａｃｙ　ｏｆ　ＳＰＨ　ｍｅｔｈｏｄ　ｉｎ　ｓｏｌｖｉｎｇ　ｔｈｅ　ｐｒｏｂｌｅｍ　ｏｆ　ｔｈｒｅｅ—ｄｉｍｅｎｓｉｏｎａｌ　ｌｉｑｕｉｄ　ｓｌｏｓｈｉｎｇ　ｉｓ　ｖａｌｉｄａｔｅｄ．Ａｎａｌｙｓｉｓ　ｒｅｓｕｌｔｓ　ｓｈｏｗ　ｔｈａｔ　ｔｈｅ　ｆｉｌｌ　ｒａｔｉｏ　ｏｆ　ｔｈｅ　

ｌｉｑｕｉｄ　ａｎｄ　ｔｈｅ　ｓｌｏｓｈｉｎｇ　ｃｙｃｌｅ　ｈａｓ　ｉｍｐｏｒｔａｎｔ　ｉｎｆｌｕｅｎｃｅ　ｏｎ　ｔｈｅ　ｃｈａｒａｃｔｅｒｉｓｔｉｃｓ　ｏｆ　ｔｈｅ　ｌｉｑｕｉｄ　ｓｌｏｓｈｉｎｇ　ａｎｄ　ｔｈｅ　

ｗａｌｌ　ｐｒｅｓｓｕｒｅ．Ａｔ　ｔｈｅ　ｌａｒｇｅ　ｆｉｌｌ　ｒａｔｉｏ　ｏｆ　ｔｈｅ　ｌｉｑｕｉｄ，ｔＷＯ—ｐｅａｋ　ｆｅａｔｕｒｅ　ｒｅｖｅａｌｅｄ　ｂｙ　ｔｈｅ　ｈｙｓｔｅｒｅｓｉｓ　ｃｈａｒａｃｔｅｒ—　

ｉｓｔｉｃ　ｄｉｓａｐｐｅａｒｓ．　Ｋｅｙ　ｗｏｒｄｓ：ｓｍｏｏｔｈｅｄ　ｐａｒｔｉｃｌｅ　ｈｙｄｒｏｄｙｎａｍｉｃｓ；ｎｕｍｅｒｉｃａｌ　ｓｉｍｕｌａｔｉｏｎ；ｓｌｏｓｈｉｎｇ；ｌｉｑｕｉｄ　ｔａｎｋ　

贮箱内液体的晃动是指液体自由表面由于外　加扰动或激励而产生的波运动。液体晃动问题广泛　

存在于航天飞行器推进剂贮箱、飞机燃油箱设计、　

载液货船液舱结构设计等领域，从２Ｏ世纪中期开　始，已得到各国研究人员的广泛关注。　

对液体晃动问题的早期研究大多基于液体简　

化模型等线性理论进行，即将液体等效为单摆或者　弹簧一脉冲质量系统。液体简化模型对于小幅晃动　

是有效的，但在大幅晃动情况下，非线性效应会非　常明显，线性理论难以真实地预测自由液面的位　置。对于液体大幅晃动，一般需要采用数值解法进　行求解，如Ｈａｒ—ｌｏｗ提出的ＭＡＣ法，Ｈｉｒｔ提出的　

ＶＯＦ法，以及有限元法、边界元法等。但是ＭＡＣ法　

对硬件要求高，计算时间长；ＶＯＦ法需要高精度的　偏微分方程离散格式，有限差分方法、有限元法、边　

界元法在涉及到大变形如液面破碎或者翻卷的问　题上，会出现网格的缠结和扭曲，严重影响到计算　

精度。光滑粒子流体动力学（Ｓｍｏｏｔｈｅｄ　ｐａｒｔｉｃｌｅ　ｈｙ—　ｄｒｏｄｙｎａｍｉｃｓ，ＳＰＨ）方法作为无网格方法的一种，　

是由Ｇｉｎｇｏｌｄ和ＭｏｎａｇｈａｎＬ１　于１９７７年提出的。　

收稿日期：２００９—０２—２３；修订日期：２００９—０９—０２　作者简介：刘富，男，博士研究生，１９８３年生；童明波（联系人），男，教授，博士生导师，Ｅ—ｍａｉｌ：ｔｏｎｇｗ＠ｎｕａａ．ｅｄｕ．ｃｎ。

　第１期　刘　富，等：基于ＳＰＨ方法的三维液体晃动数值模拟　１２３　

ＳＰＨ方法最初是用于解决三维开放空间的天体物　理学问题，当Ｍｏｎａｇｈａｎ首次将ＳＰＨ方法应用于自　

由表面水波的流动模拟　后，该方法被广泛应用于　

具有材料强度的动态响应和具有大变形的流体动　

力学问题　引。由于ＳＰＨ方法的计算对象为空间运　

动的粒子，在此方法中，每一个时间步内都要应用　当前支持域内的粒子进行粒子计算，且所用的粒子　

近似式是由拉格朗日描述下的偏微分方程组推导　而得。因此自适应、元网格、粒子形式以及拉格朗日　

性质的ＳＰＨ方法避免了大变形时，网格发生畸变　

的问题，非常适用于液体大幅晃动问题的研究。　在国内，近年来逐渐有学者采用ＳＰＨ方法进　行液体晃动的研究，但基本都是基于二维晃动的数　值计算［　］，目前关于液体三维晃动的研究却很　

少。本文基于ＳＰＨ方法的基本理论，采用ＳＰＨ程　

序，进行液体三维晃动的数值计算，为今后进行飞　

机油箱燃油晃动，航天器推进剂晃动等方面的研究　奠定基础。　

１　ＳＰＨ数值方法　

１．１　ＳＰＨ基本原理　ＳＰＨ方法通过构造一个近似场函数（厂（ｚ））来　

表征某段时间内粒子的速度和能量，表达式如下　

＜厂（ｚ））＝Ｉ　ｆ（ｘ　）　（ｚ—　，ｈ）ｄｘ　（１）　

式中：　（ｚ—ｚ　，　）为光滑核函数；ｈ为光滑长度，　

它是粒子　与Ｘ　的间距，决定了核函数支持域的尺　

寸。　常用的核函数是Ｍｏｎａｇｈａｎ提出的分段三次　

样条函数　

ｌ号一Ｒ　＋寺Ｒ。０≤Ｒ＜１　Ｉ　

（Ｒ，＾）一ａａ１　（２一Ｒ）ｓ　１≤Ｒ≤２　ｌ　＼　＾　＼　ｂ　

。　尺＞２　（２）　式中尺为粒子ｚ和ｚ　之间的相对距离。在一维、二　

维、三维空间中，Ｏｆｄ分别为１／ｈ，１５／７　。，３／２￣ｒｈ。。本　

文选取该样条函数作为光滑核函数。　若设粒子Ｊ的密度为Ｐｘ（ｊ一１，２，…，Ⅳ），质量　

为ｍ』，其中Ⅳ为粒子ｉ支持域里的粒子总量，对式　

（１）离散化可得粒子ｉ处的粒子近似式　

＜厂（ｚ　）＞一∑　厂（ｚ　）・　（五一　，　）（３）　

利用核函数的可微性，函数空间导数的粒子近　似式为　（　一　ＶＷ（一　）　

（４）　从式（３，４）可以看出，粒子近似式将函数及导　

数的连续积分表达式转换成对任意排列的粒子的　离散化求和，从而实现了空间的无网格化。　

１．２　ＳＰＨ流体运动方程　

描述连续介质流体动力学的质量守恒，动量守　恒方程分别表示为　

一一ＩＤ　．比　（５）　ｄ￡　

ｄｔ一一丢　一　（６）　Ｄ　Ｖ　＼ｕ　

式中：ＩＤ为流体密度；ｕ为速度向量；　为应力张量；　

ｔ为时间。通过对上述Ｎ—Ｓ方程的转换，可以得到方　程的ＳＰＨ表达式如下ｍ］　

ｄ　Ｐ＿＿２＝一　（“　—　・　』　（７）　

警一砉　（ａ　ｔ．＋　㈣　

式中Ｗ　一　（　—ｚ　，＾）。　对于不可压缩流，通过人工压缩性口　可以将　

其视为可压缩流，粒子ｉ的压强Ｐ　由流体状态方程　求得　

一　）　㈤　

式中：　为流体初始密度；对于水，系数７＝－７；Ｃ　为　

表征粒子声速的参数，Ｃ　一√１００ｐ。ｇＨ，Ｈ为贮箱　

内液面高度。　１．３固壁边界条件　

为了阻止固壁边界被粒子穿透，在ＳＰＨ方法　中，Ｍｏｎａｇｈａｎ在固定边界上分布一组虚粒子用于　

对邻近边界的粒子作用排斥力，从而阻止邻近边界　

的这些粒子非物理穿透边界。此时在虚粒子和相邻　

实粒子之间的作用力沿着两粒子的中心线，力的大　小为　

一　一　］考　≤１　

ｌ　０　ｒｏ＞　

（１０）　式中：Ｄ，　和ｎ　为可调参数。文中Ｄ的取值为　

ｌＯｇｈ，ｈ为液面高度，ｎ　和　一般取１２和４。　为当　

前时刻虚粒子和邻近边界实粒子的间距，ｒ０为截止　半径：若太大，会造成初始阶段，实粒子就受到排斥　

力的作用，扰乱粒子的分布；若太小，粒子可能穿过　１２４　南京航空航天大学学报　第４２卷　

边界也感受不到排斥力作用，其取值一般与粒子初　

始间距的大小相近，本文取其值为粒子的初始间距　

０．００５　ｍ。　

２数值算例　

液舱内液体的晃动是大型载液货船发展的关　

键技术之一，当货船的晃动频率接近于液体的固有　

频率时，会带来液舱内液体的大幅晃动，其引起的　

晃动冲击载荷会对液舱壁及内部构件产生破坏性　的影响，带来一系列的强度和疲劳问题。　棱形液舱是大型载液货船储液系统设计常采　用的一种液舱形式。Ｍｉｋｅｌｉｓ等［１ｓ３在实验中通过对　

棱形液舱施加横摇及俯仰激励，测量液体晃动对液　

舱侧壁及顶部产生的冲击压力。液舱的晃荡实验满　

足模型相似律，即几何相似、运动相似、动力相似、　边界条件相似及初始条件相似。实验模型按１：４０　

的比例制作，在液舱的中心剖面沿宽度方向侧壁上　布置压力测试点，图１显示了具体的液舱尺寸及文　

中计算工况所用到的压力测试点布置图。Ｒ　，Ｒ。，　

尺。，Ｒ　分别距离液舱底面７５，１１３．５，２０５和　

２９６．５　ｍｍ。充液￣Ｌｈ／Ｈ分别为０．３３和０．６１，Ｈ为　

液舱的高度。外加激励为绕液舱底部中心线的横摇　运动，运动方程为　一０。ｓｉｎ（ｃｏｔ）。　

本文利用ＳＰＨ程序，对实验中的两个工况（表　

１）进行计算，并与实验进行对比。计算中不考虑液　

舱的变形，设置为刚性模型。　表１计算工况　

图ｌ棱形液舱模型　

当液舱晃动的时候，载液深度和晃动周期不同　会产生不同的晃荡波浪。常见的舱内液体晃荡现象　

有：驻波、行进波、水跃以及组合波，另外还有冲击　作用下的液体飞溅等强非线性现象。图２直观显示　

了液舱在外加横摇激励作用下自由液面的晃动模　

拟及典型晃荡现象。从图中可以清楚看出液面的翻　卷以及应用有限元方法难以模拟的碎波。　

…力　

瓤力　

（ｅ）驻波（工况２）　（ｆ）２￣（Ｔｂ２２）　

图２自由液面的晃动模拟　

图３，４分别给出了ＳＰＨ程序所计算的压力测　

试点的压力曲线图，并且和实验进行了对比。从图　

３可以看到测试点的压力呈现的双峰特征，这是液　体晃动过程中迟滞特性的一种典型表现。其中前一　

个峰值是由于液体对舱壁之间的冲击而产生的高　脉冲压力，时间短，具有明显的局部特性，主要发生　

在水跃或者行进波过程中（图２（ｃ））。后一个峰值也