肖 瀟 ， 蔣昌波 ， 程永舟

（1長(zhǎng)沙理工大學(xué) 水利工程學(xué)院，長(zhǎng)沙410004；2湖南省水沙科學(xué)與水災(zāi)害防治重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙 410004）

水流對(duì)浮體作用的SPH方法模擬

肖 瀟1，2， 蔣昌波1，2， 程永舟1，2

（1長(zhǎng)沙理工大學(xué) 水利工程學(xué)院，長(zhǎng)沙410004；2湖南省水沙科學(xué)與水災(zāi)害防治重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙 410004）

水流與浮體相互作用時(shí)，運(yùn)動(dòng)情況十分復(fù)雜，目前很難準(zhǔn)確有效地模擬。而SPH作為一種純拉格朗日方法，無需構(gòu)建網(wǎng)格，用核函數(shù)近似粒子進(jìn)行離散，能較好地解決一些自由面大變形問題。文章利用SPH法對(duì)潰壩時(shí)引起的高速水流沖擊浮體以及水體晃動(dòng)時(shí)破損浮體的運(yùn)動(dòng)過程進(jìn)行模擬。模擬結(jié)果表明，SPH法能有效地進(jìn)行水流對(duì)浮體作用的研究。

浮體；光滑粒子流體動(dòng)力學(xué)；移動(dòng)最小二乘法；核函數(shù)；數(shù)值模擬

1 引 言

光滑粒子流體動(dòng)力學(xué)（Smoothed Particle Hydrodynamics）法[1]是模擬流體流動(dòng)的一種無網(wǎng)格拉格朗日粒子法,最初被用于解三維開放空間天體物理學(xué)問題?；舅枷胧菍⑦B續(xù)系統(tǒng)用一系列粒子來描述，這些粒子攜帶獨(dú)立的物理信息并遵循物理控制方程進(jìn)行運(yùn)動(dòng)。這種方法無需網(wǎng)格及質(zhì)點(diǎn)間相互作用的計(jì)算，使得大變形問題的處理變得容易。

近年來SPH法不斷發(fā)展、改進(jìn)和完善,逐漸趨向于成熟，并應(yīng)用于很多領(lǐng)域。Sigalotti等[2]運(yùn)用SPH法模擬了強(qiáng)沖擊過程，Bulian等[3]運(yùn)用SPH法模擬了不同角度情況下箱中水體的晃動(dòng)。Marongiu等[4]結(jié)合SPH-ALE法模擬了水輪機(jī)的自由表面流動(dòng)，Capone等[5]運(yùn)用SPH法進(jìn)行了滑坡涌浪的模擬，以及Narayanaswamy[6]將SPH與Funwave結(jié)合進(jìn)行了近岸波浪傳播過程的模擬等。在國內(nèi)也有一些學(xué)者進(jìn)行了SPH法的研究應(yīng)用，宗智等[7]進(jìn)行了二維水下爆炸問題的SPH數(shù)值模擬，崔巖等[8]進(jìn)行了二維矩形水槽晃蕩過程的SPH數(shù)值模擬，李大鳴等[9]進(jìn)行了大壩表孔泄流的SPH數(shù)值模擬，但作為一個(gè)日趨成熟和完善的方法，SPH法的優(yōu)勢(shì)還未被完全發(fā)揮。

在海洋、河流中高速水流對(duì)船舶等水上漂浮物的運(yùn)動(dòng)影響很大，使得它們的運(yùn)動(dòng)軌跡發(fā)生很大變化，因此，需要研究其運(yùn)動(dòng)規(guī)律，降低水流沖擊的破壞程度。對(duì)這類問題有不少學(xué)者進(jìn)行了相關(guān)研究。Contento等[10]簡(jiǎn)單地計(jì)算了波浪對(duì)不規(guī)則浮體的作用，Had?i等[11]用網(wǎng)格法進(jìn)行了波浪作用下浮體運(yùn)動(dòng)情況的模擬，Yan等[12]用QALE-FEM方法進(jìn)行了非線性波與浮體相互作用的模擬，王大國等[13]進(jìn)行了波浪對(duì)箱形船作用的計(jì)算，王化明等[14]進(jìn)行了船舶在淺水域橫向停靠運(yùn)動(dòng)粘性水動(dòng)力計(jì)算。這些研究都為解決實(shí)際問題提供了良好的依據(jù)?；跉W拉網(wǎng)格的方法，由于網(wǎng)格的局限性，對(duì)高速水流作用下浮體的大幅度運(yùn)動(dòng)情況下不能很好地進(jìn)行模擬，而SPH法卻能很好地解決這個(gè)問題。

本文利用SPH方法，對(duì)潰壩引起的高速水流對(duì)浮體的沖擊作用，以及水體晃動(dòng)對(duì)開口浮體影響進(jìn)行了二維數(shù)值模擬，以此分析高速水流沖擊作用下浮體的運(yùn)動(dòng)過程。

2 計(jì)算模型

2.1 SPH基本原理

在SPH方法中，任一宏觀變量都能借助于一組無序點(diǎn)上的值表示成積分插值計(jì)算得到。對(duì)于任意函數(shù)A()r[1]有：

將（1）式轉(zhuǎn)化為支持域內(nèi)所用粒子疊加求和的離散化形式，則在粒子a處的函數(shù)可表示為：

式中：mb為粒子質(zhì)量；ρb為粒子密度，Wab=ra-rb′,( )h為核函數(shù)。

SPH方法中所用的核函數(shù)應(yīng)具有歸一性、緊支性、非負(fù)性、衰減性、δ函數(shù)性質(zhì)、對(duì)稱性質(zhì)和光滑性等特性[9，15]。目前，使用較多的核函數(shù)有Gaussian函數(shù)、B樣條函數(shù)、高次核樣條函數(shù)等,本次數(shù)值模擬采用目前應(yīng)用比較廣泛的三次樣條核函數(shù)：

式中：αD應(yīng)用于二維問題時(shí)?。籷為相對(duì)長(zhǎng)度，q=。

2.2 基本方程

SPH方法引入基本方程進(jìn)行離散，采用人工粘性來產(chǎn)生剪切粘性[16]和體積粘性，得到離散形式的基本方程組為：

浮體的運(yùn)動(dòng)情況控制采用以下方程：

式中：fk是邊界粒子k單位質(zhì)量所受的作用力；fka是水粒子a對(duì)邊界粒子k的單位質(zhì)量作用力；WPs是對(duì)邊界粒子k有影響的水粒子；mk是粒子k的質(zhì)量；ma是粒子a的質(zhì)量；M是浮體的質(zhì)量；V是浮體的速度；BPs是浮體的邊界粒子；I是慣性力；Ω是浮體的旋轉(zhuǎn)速度；R0是重心的位置；uk是邊界粒子k的速度。

粒子密度采用移動(dòng)最小二乘法修正[17]后的密度求和公式進(jìn)行求解，該方法能有效地克服邊緣效應(yīng)。邊界條件采用Dalrymple（2000）[18]提出的動(dòng)力邊界條件，該方法使邊界粒子滿足流體粒子一樣的方程，但是它們需滿足固壁不可滑移的條件。

3 模型驗(yàn)證

圖1為潰壩實(shí)物模型示意圖，水箱長(zhǎng)、寬、高為4m×0.1m×3m，水體長(zhǎng)、寬、高為1m×0.1m×2m,水體前放置一擋板[19]。數(shù)模計(jì)算參照實(shí)體模型的剖面進(jìn)行相應(yīng)的二維模擬設(shè)置，取粒子的初始速度為0，粒間距為0.012m，水體粒子的總數(shù)為27 472個(gè)，進(jìn)行潰壩模擬。

圖2為典型時(shí)刻潰壩過程中水粒子速度分布圖，形象地給出了流體和氣體界面及流速隨時(shí)間的演進(jìn)過程。在t＞0時(shí)刻，擋板瞬間移開，水體的勢(shì)能釋放為動(dòng)能，將快速向右下方流動(dòng)，水體右下方底邊界附近速度較大；右上角由于水體的下落，速度也逐漸增大。水體運(yùn)動(dòng)過程中，波前峰快速覆蓋底部干邊界，接近0.8s時(shí)，波前峰抵達(dá)右側(cè)邊壁，此時(shí)動(dòng)能沒有完全耗散，水體在慣性作用下沿右側(cè)邊壁繼續(xù)向上運(yùn)動(dòng)，動(dòng)能轉(zhuǎn)化為勢(shì)能。當(dāng)水流運(yùn)動(dòng)接觸到上壁邊界時(shí)，由于壁面摩檫力和重力作用，達(dá)到極限運(yùn)動(dòng)高度，開始回落。圖3為水流前進(jìn)過程中波峰到達(dá)位置與時(shí)間的關(guān)系圖，與實(shí)測(cè)數(shù)據(jù)比較[19，20]，最大相對(duì)誤差不超過10%，結(jié)果合理，能很好的反映真實(shí)流動(dòng)。

通過對(duì)潰壩過程進(jìn)行數(shù)學(xué)模型驗(yàn)證表明：利用SPH方法進(jìn)行數(shù)值模擬結(jié)果合理可靠，在流體大變形的模擬上、體積守恒上具有相對(duì)的優(yōu)勢(shì)，自由表面流的模擬結(jié)果合理，因而，該方法能夠很好地應(yīng)用于相關(guān)問題的研究。

圖1 潰壩試驗(yàn)?zāi)Ｐ褪疽鈭DFig.1 Collapses of a water column test model diagram

圖2 潰壩過程圖Fig.2 Collapses of a water column in a tank

圖3 水流到達(dá)位置與時(shí)間的關(guān)系Fig.3 The relation of location flow arriving time

4 水流對(duì)浮體作用的數(shù)值模擬

4.1 潰壩時(shí)高速水流對(duì)浮體的沖擊作用

圖4為假定的潰壩水流對(duì)浮體沖擊的計(jì)算模型平面示意圖。浮體邊長(zhǎng)為0.25m的正方矩形,重1kg，放在水體右下角0.5m高的位置,在槽中2m的位置布置一個(gè)高1m的障礙物。數(shù)學(xué)模型參照實(shí)體模型進(jìn)行相應(yīng)的二維模擬設(shè)置，粒子和其它條件設(shè)置都與此前潰壩模型算例相同。

圖5給出了潰壩引起的高速水流對(duì)浮體的沖擊以及障礙物的阻礙作用下整個(gè)流動(dòng)過程隨時(shí)間變化的分布情況。從圖中可以看出擋板被瞬間移開后，水體向右下方塌落并帶動(dòng)浮體向右運(yùn)動(dòng)；遇障礙物后，水體在慣性作用下帶動(dòng)浮體向上運(yùn)動(dòng)，繞過障礙物，波峰處水體向前方飛濺破碎；在重力作用下，水體和浮體向下回落。

圖6給出了浮體在整個(gè)沖擊過程中所處位置的情況，結(jié)合圖5可以看出，最初浮體受重力作用下落，后經(jīng)水流帶動(dòng)移至障礙物前；在慣性作用下隨水流上移，達(dá)最高位置為1.75m；在重力作用下回落，由水流推動(dòng)達(dá)最右端，即3.25m處，后經(jīng)墻壁反射回來的水流作用移至2.5m處。

經(jīng)過對(duì)SPH方法模擬潰壩引起的高速水流對(duì)浮體的沖擊作用的流動(dòng)過程，該模擬結(jié)果接近實(shí)際情況，但是具體精確程度需要相關(guān)模型試驗(yàn)進(jìn)一步驗(yàn)證。

圖4 潰壩水流對(duì)浮體沖擊的模型示意圖Fig.4 Flow impact to floating body test model diagram

圖5 潰壩水流對(duì)浮體沖擊的過程Fig.5 Process of flow impact to floating body

圖6 浮體整個(gè)運(yùn)動(dòng)過程的位置情況圖Fig.6 Floating body position throughout the movement

圖7 晃動(dòng)水流對(duì)浮體作用的模型示意圖Fig.7 Sloshing water on floating body model diagram

4.2 水流晃蕩對(duì)浮體的作用

圖7為假定的水流晃蕩對(duì)浮體作用的計(jì)算模型示意圖。浮體邊長(zhǎng)為0.25m的上方開口的立方體,重0.5kg，水流初始布置為一傾斜的水體。數(shù)學(xué)模型參照實(shí)體模型進(jìn)行相應(yīng)的設(shè)置，粒子的初始速度都為0，粒間距為0.02m，水體粒子的總數(shù)為19 751個(gè)。

圖8給出了開口浮體在晃動(dòng)水流作用下整個(gè)流動(dòng)形態(tài)過程隨時(shí)間變化的分布情況。傾斜的水體在重力作用下，向右回落，浮體隨液面移動(dòng)；右側(cè)水體在慣性作用下，沿壁向上運(yùn)動(dòng)，水體灌入破損浮體；重力作用下，右側(cè)水體回落，浮體在水流作用下翻轉(zhuǎn)；左側(cè)水體在慣性作用下沿壁向上運(yùn)動(dòng)，浮體下沉。由模型計(jì)算結(jié)果可知，SPH可以模擬在水體晃蕩作用下，水流灌入開口浮體的情況。

圖8 水體晃動(dòng)對(duì)浮體作用的過程Fig.8 Process of the sloshing water on floating body

5 結(jié) 論

本文基于SPH方法，對(duì)潰壩引起的高速水流對(duì)浮體的沖擊作用，以及水體晃動(dòng)對(duì)開口浮體的影響進(jìn)行了二維數(shù)值模擬，分析結(jié)果得出以下結(jié)論：

（1）利用SPH法對(duì)潰壩水流進(jìn)行數(shù)模驗(yàn)證，表明SPH法可用于比較精確地模擬水流的非線性大變形流動(dòng)問題。

（2）相對(duì)于歐拉等網(wǎng)格法，SPH法作為一種拉格朗日法能有效地追蹤單個(gè)粒子的運(yùn)動(dòng)過程，保證整個(gè)流體的總體積不變，更能真實(shí)地反映實(shí)際情況。

（3）SPH方法能很好地追蹤浮體的運(yùn)動(dòng)軌跡，準(zhǔn)確模擬高速水流沖擊作用下浮體大幅度運(yùn)動(dòng)的過程。

（4）SPH方法能進(jìn)行水流灌入開口浮體的模擬，很好地反映浮體在波浪作用下的翻轉(zhuǎn)過程。

（5）SPH方法是值得研究發(fā)展的方法，今后可以做進(jìn)一步改進(jìn)，并進(jìn)行相關(guān)的實(shí)物模擬，以獲得更廣闊的研究空間。

[1]Monaghan J J.Simoothed particle hydrodynamics[J].Annual Rev.Astron.Appl.,1992,30:543-574.

[2]Sigalotti L Di G,Lopez H,Trujillo L.An adaptive SPH method or strong shocks[J].Journal of Computational Physics,2009,228:5888-5907.

[3]Bulian G,Souto-Iglesias A,Delorme L.Smoothed particle hydrodynamics simulation of a tuned liquid damper with angular motion[J].Journal of Hydraulic Research,Extra Issue,2010,48:28-39.

[4]Marongiu J C,Leboeuf F.Free surface flows simulations in Pelton turbines using an hybrid SPH-ALE method[J].Journal of Hydraulic Research,Extra Issue,2010,48:40-49.

[5]Capone T,Panizzo A.SPH modelling of water waves generated by submarine landslides[J].Journal of Hydraulic Research,Extra Issue,2010,48:80-84.

[6]Narayanaswamy M.SPHysics-FUNwave hybrid model for coastal wave propagation[J].Journal of Hydraulic Research,Extra Issue,2010,48:85-93.

[7]宗 智,鄒 麗,劉謀斌等.模擬二維水下爆炸問題的光滑粒子（SPH）方法[J].水動(dòng)力學(xué)研究與進(jìn)展,2007,22(1):61-67.

[8]崔 巖,吳 衛(wèi),龔 凱等.二維矩形水槽晃蕩過程的SPH方法模擬[J].水動(dòng)力學(xué)研究與進(jìn)展A輯,2008,23(06):618-624.

[9]李大鳴,劉江川,徐亞男等.SPH法在大壩表孔泄流數(shù)值模擬中的應(yīng)用[J].水科學(xué)進(jìn)展,2008,16(6):841-845.

[10]Contento G.Numerical wave tank computations of nonlinearmotions of wo-dimensional arbitrarily shapedfree floating bodies[J].Ocean Engineering,2000,27:531-556.

[11]Had?i I,Hennig J,Peri a M.Computation of flow-induced motion of floating bodies[J].Applied Mathematical Modelling,2005,29:1196-1210.

[12]Yan S,Ma Q W.Numerical simulation of fully nonlinear interaction between steep waves and 2D floating bodies using the QALE-FEM method[J].Journal of Computational Physics,2007,211:666-692.

[13]王大國,鄒志利,唐春安等.波浪對(duì)箱形船作用的三維耦合計(jì)算模型[J].船舶力學(xué),2007,11(4):533-544.

[14]王化明,鄒早建,田喜民等.船舶在淺水域橫向?？窟\(yùn)動(dòng)粘性水動(dòng)力計(jì)算[J].船舶力學(xué),2010,14(7):723-731.

[15]Jin Hongbin,Ding Xin.On criterions for smoothed particle hydrodynamics kernels in stable field[J].Journal of Computational Physics,2005,202:699-709.

[16]Monaghan J J.Simulating free surface flow with SPH[J].Journal of Computational Physics,1994,110:399-406.

[17]Dilts G A.Moving-Least-Squares-Particle Hydrodynamics-I.Consistency and stability[J].Int.J Numer.Meth.Engng,1999,44(8):1115-1155.

[18]Dalrymple R A,Knio O.SPH modelling of water waves[C]//Proc.Coastal Dynamics.Lund,2000:779-787.

[19]Koshizuka S,Oka Y.Moving-particle semi-implicit method for fragmentation of incompressible fluid[J].Nucl Sci Eng,1996,123:421-434.

[20]Gomez-Gesteira M,Rogers B D.State-of-the-art of classical SPH for free-surface flows[J].Journal of Hydraulic Research,Extra Issue,2010,48:6-27.

Simulation of flow-induced floating-body motion with SPH method

XIAO Xiao,JIANG Chang-bo,CHENG Yong-zhou
(1 School of Hydraulic Engineering,Changsha University of Science&Techology,Changsha 410004,China;2 Hunan Province Key Laboratory of Water,Sendiment Sciences&Flood Hazard Prevention,Changsha 410004,China)

Flow-induced floating body motions,that movement is very complicatied,can not be simulated exactly at present.SPH as a Lagrangian method without component grid,uses kernel fuction approximate to particle discretely.And it can solve some problems with strong deformation of free surface.In this paper,SPH method is applied to simulate the process of the floating body motion which results from the collapse of a water column and the movement of damaged floating body in sloshing water.Simulation results show that the SPH method can effectively study the flow indued motion of floating body.

floating body;smoothed particle hydrodynamics （SPH）;moving least squares;kernel function;numerical simulation

TV13 U661.1

A

1007-7294（2011）08-0861-06

2010-09-25 修改日期：2011-02-18

國家自然科學(xué)基金項(xiàng)目（50909009，50979008）；水文水資源與水利工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室開放研究基金項(xiàng)目（2008490911）；教育部高校博士點(diǎn)基金項(xiàng)目（20094316110002）

肖 瀟（1986-），女，碩士研究生，研究方向?yàn)樗W(xué)及河流動(dòng)力學(xué)，Email:xxxtina@163.com；蔣昌波（1970-），男，長(zhǎng)沙理工大學(xué)教授，博士生導(dǎo)師。