礫石單元法與SPH 耦合模型在數(shù)值波浪水槽中的應(yīng)用研究

李彥卿, 別社安, 李大鳴, 王 鑫

(1. 天津大學(xué) 水利工程仿真與安全國家重點(diǎn)實(shí)驗(yàn)室, 天津 300072; 2. 國家海洋技術(shù)中心, 天津 300112)

光滑粒子動力學(xué) (Smoothed Particle Hydrodynamics, 簡記為SPH)方法是一種處理大變形邊界運(yùn)動問題的有效方法。20 世紀(jì)70 年代為研究天體運(yùn)動的物理問題, Gingold 與Monaghan 提出了用追蹤粒子的方法描述宇宙天體中星云的運(yùn)動過程[1]。1992 年Monaghan 將SPH 方法擴(kuò)展到求解不可壓縮流體的大變形運(yùn)動問題[2]。其后許多學(xué)者在SPH 的應(yīng)用中做了多方面研究[3-6]。離散單元法(Discrete Element Method, 簡記為DEM)最早是由Peter Cundall 博士于1971 年提出的[7-8]。DEM 的基本思想是利用一系列離散的粒子來模擬固體介質(zhì), 其中每個粒子為一個單元,具有固定的大小和形狀, 粒子與粒子之間可以在小范圍內(nèi)相互擠壓重疊從而產(chǎn)生相互作用力, 在通過牛頓運(yùn)動定律計算出粒子的運(yùn)動參數(shù), 模擬固體介質(zhì)的運(yùn)動情況。近年來SPH-DEM 耦合模型也有新的研究成果[9-11]。SPH 與DEM 同時離散為粒子時, 兩者顆粒尺度相當(dāng), 只是在顆粒性質(zhì)上有區(qū)別[12]; 當(dāng)SPH 離散為粒子, DEM 用組合形式表示固體物面時, 主要描述大尺度、數(shù)量不多、結(jié)構(gòu)形式確定(可以有一定變形)、具有獨(dú)立運(yùn)動(一般為懸浮)的固體物面。對堆積礫石來說,與SPH 流體粒子尺度不同, 堆積礫石數(shù)量較多, 礫石輪廓形式復(fù)雜(可以統(tǒng)一?；?, 堆積礫石相互擠壓、遮擋運(yùn)動相互制約。為解決堆積礫石的群體運(yùn)動問題, 本文提出了礫石單元法(Gravel Element Method, 簡記為GEM), 建立了單礫石近似為力學(xué)球的概念, 用球度系數(shù)和形狀系數(shù)對礫石幾何形態(tài)進(jìn)行近似; 對于不規(guī)則的礫石輪廓, 用傅里葉級數(shù)展開礫石邊界輪廓, 將礫石表面的力學(xué)特征賦予在周期函數(shù)變化的展開函數(shù)中,在力學(xué)球形態(tài)的礫石表面增加對應(yīng)非球形礫石位置處的滾動摩擦阻力系數(shù)、滑動摩擦阻力系數(shù)和轉(zhuǎn)動慣量分布。研究了有級配非圓形狀顆粒對給定壩型的隨機(jī)填充方法, 對堆積礫石力學(xué)球的受力狀態(tài)進(jìn)行了分析,提出了擬序排列、分級求解和靜力傳遞的方法。建立了SPH-GEM 耦合波浪水槽和堆積礫石塌落、滾落的數(shù)學(xué)模型。

1 堆積礫石的受力分析

1.1 單礫石幾何形態(tài)

假定礫石顆粒內(nèi)部密度均勻, 在邊緣輪廓線段的基礎(chǔ)上, 可以求出礫石顆粒投影平面θ上的圖形形心(xc, yc), 以形心為極點(diǎn), 可以確定顆粒邊緣輪廓曲線的極角φ和極坐標(biāo)半徑R, 顆粒邊緣輪廓線的半徑以T=2π 為周期循環(huán)變化, 可以表示為R(θ, 2π+φ)=R(θ,φ)。

理論上講礫石顆粒具有三維特征, 可以在球坐標(biāo)中將礫石的任意切面以平面坐標(biāo)展開, 足夠的切面組合可以描述礫石的三維特征, 研究任意切面上礫石的輪廓則具有二維特征。

對不規(guī)則的礫石輪廓, 在θ平面上用傅里葉級數(shù)來逼近周期函數(shù)R(θ,φ), 得

式中:n=1, 2, 3, ...為級數(shù)的各項(xiàng)級數(shù);a0為常數(shù)項(xiàng), 表示顆粒投影平面上的平均粒徑;an和bn為傅里葉級數(shù)各項(xiàng)系數(shù), 表示在平均粒徑基礎(chǔ)上的半徑修正量。

傅里葉系數(shù)由以下公式確定

當(dāng)an和bn均為零時, 礫石顆粒半徑為,a0為圓形礫石。低頻諧波(n較小)反映礫石顆粒表面較大尺度的形態(tài)變化, 高頻諧波(n較大)反映礫石顆粒表面較小尺度的形態(tài)變化。諧波越多, 傅里葉級數(shù)生成的輪廓線越接近于實(shí)際礫石顆粒的輪廓線。圖1a 為礫石顆粒的圖像采集, 圖1b 為礫石顆粒的輪廓模擬。

圖1 礫石顆粒圖像和輪廓模擬Fig. 1 Image of gravel particles and contour simulation

礫石與球體的偏離程度可以用球度系數(shù)和形狀系數(shù)來描述[13], 球度系數(shù)公式為[14]

式中:a為礫石長軸;b為礫石中長軸;c為礫石短軸。

形狀系數(shù)公式為[14]

對退化的二維礫石顆粒, 可以假定c=b。對二維礫石顆粒, 引用公式(1)中傅里葉級數(shù)各項(xiàng)系數(shù)an、bn作為傅里葉級數(shù)各級礫石的長軸和中長軸, 可以將非圓顆粒形狀用數(shù)個圓顆粒部分表面的組合形狀來近似(圖2)。礫石的二維圓周矢量半徑展開曲線為圖3。

圖2 非圓顆粒的組合形狀Fig. 2 Combined shape of non-circular particles

圖3 礫石的二維圓周矢量半徑分布Fig. 3 Distribution of the 2D circumferential vector radius of the gravel

1.2 單礫石力學(xué)球概念

通過考察不同投影平面的礫石傅立葉級數(shù)各級系數(shù)半徑, 可以得到礫石的球度系數(shù)或形狀系數(shù),不僅從礫石表面整體上獲取礫石對球的近似程度,而且可以對不同位置的球面近似形態(tài)進(jìn)行描述?？梢詫⒛嗌愁w粒的球度系數(shù)延伸到礫石的球態(tài)研究中,表1 為礫石幾何形狀與球度系數(shù)的關(guān)系。

表1 礫石幾何形狀與球度系數(shù)的關(guān)系Tab. 1 Relationship between gravel geometry and sphericity coefficient

用球體近似礫石不僅是通過球度系數(shù)對礫石形態(tài)的形似, 而且要賦予球表面上礫石的力學(xué)性質(zhì)。球形礫石容易滾動, 與其他球形顆粒主要為點(diǎn)接觸,表面相對阻力較小, 具有較均勻的表面摩擦力和不變的轉(zhuǎn)動慣量。而非球形礫石的棱角有時會阻礙滾動, 有時會使顆粒迅速翻轉(zhuǎn), 在各部位轉(zhuǎn)動慣量不一致, 局部較平緩粗糙的表面會增加表面摩擦力。為達(dá)到用球形礫石模擬非球形礫石運(yùn)動的目的, 引入力學(xué)球的概念, 即當(dāng)球形礫石與非球形礫石的力學(xué)特征相同時, 這種球形礫石稱為力學(xué)球。

力學(xué)球以球的形狀存在, 用球度系數(shù)作為礫石與球體的形狀修正系數(shù), 在球形礫石表面增加對應(yīng)非球形礫石位置處的滾動摩擦阻力系數(shù)、滑動摩擦阻力系數(shù)和轉(zhuǎn)動慣量分布, 這樣就可以達(dá)到用球形礫石模擬非球形礫石運(yùn)動的效果。這在理論上完全可行, 但在實(shí)際操作中面對數(shù)量較多、表面復(fù)雜的堆石工程來說, 不可能對每一礫石顆粒進(jìn)行表面測量和輪廓模擬。但可以選擇有代表性的礫石, 采集礫石表面輪廓特征, 確定力學(xué)球表面的樣本函數(shù), 引入隨機(jī)分布的力學(xué)函數(shù)來近似解決這一問題。

力學(xué)球的滑動或滾動摩擦的阻力系數(shù)可以表示為

式中:fr0為球面滑動或滾動摩擦系數(shù);fr為非球面平均滑動或滾動摩擦系數(shù);F(θ)為隨θ變化的隨機(jī)函數(shù)。

力學(xué)球的轉(zhuǎn)動慣量可以表示為

式中:Js為非球面平均轉(zhuǎn)動慣量;Fs(θ)為隨θ變化的隨機(jī)函數(shù)。

1.3 堆積礫石的受力狀態(tài)

堆積礫石表層石塊所受到的作用力會通過接觸壓力、反力、碰撞傳遞給其它礫石, 使礫石進(jìn)入運(yùn)動狀態(tài)或達(dá)到新的平衡狀態(tài)。堆積礫石受力見圖4, 力學(xué)球礫石堆見圖5。

圖4 礫石受力示意圖Fig. 4 Gravel strength diagra m

可以看出內(nèi)部礫石除自身重力外, 還有礫石上面的壓力作用, 和下面礫石的反力。顯然要研究每一礫石的受力情況必須從外層礫石顆粒入手。為求解每一力學(xué)球上的受力, 可以由外及內(nèi)逐次遞進(jìn)對力學(xué)球上的作用力進(jìn)行求解, 求解順序顯得尤為重要。需要根據(jù)求解層次進(jìn)行重新編號, 將這種重新編號的過程稱為擬序排列過程, 圖6 為按求解層次編排的級別號, 圖7 為按級別求出的壓力分布。

當(dāng)取直徑d=1.0 m, 密度與重力加速度乘積ρg= 2 650 kg/m2s2。按球體計算, 累加球體重力為321 745.344 88 N, 通過傳遞力模式計算的地面總反力為321 745.344 08 N, 工程上并不需要如此精度的數(shù)值,這里主要是為驗(yàn)證傳遞力模式計算結(jié)果的可靠性; 同時地面給予力學(xué)球的摩擦力為193 047.2 N, 以保證堆積球體不滑落。

圖5 力學(xué)球礫石堆Fig. 5 Mechanical ball gravel pile

圖6 按求解層次編排的級別號Fig. 6 A level number organized by the solution level

圖7 按級別求出的壓力分布Fig. 7 Pressure distribution by grade

2 礫石單元法(GEM)

2.1 GEM 模式

GEM 是以獨(dú)立的礫石整體作為一個剛性單元,礫石表面的輪廓線是單元邊界, 礫石以整體運(yùn)動,具有平動和轉(zhuǎn)動的形式。礫石間可以通過接觸傳遞力和力矩, 具有運(yùn)動速度的礫石可以傳遞動量。在GEM 模式中引入力學(xué)球的概念, 不考慮礫石的形狀變化; 力學(xué)性質(zhì)可以分布在力學(xué)球表面, 使力學(xué)球具有靜、動摩擦系數(shù), 不同表面位置具有不同的轉(zhuǎn)動慣量。力學(xué)球之間不能壓入和穿越, 運(yùn)動受擬序排列控制, 將礫石剛體碰撞效果代替微彈性假定, 力用接觸壓力和反力模擬, 力矩用接觸摩擦和重力對支撐點(diǎn)之矩模擬[15-16]。運(yùn)動受力球動力傳遞模式如圖8所示。

圖8 礫石力學(xué)球動力傳遞示意圖Fig. 8 Schematic of the dynamic ball transmission of gravel mechanics

礫石的運(yùn)動模式主要有塌落、平移、滾落、碰撞和跳躍。作用力較小或礫石較大時平移和跳躍不容易發(fā)生, 塌落和滾落是兩種主要的運(yùn)動形式, 碰撞處理為力的傳遞。

2.2 GEM 主要力學(xué)公式

在碰撞過程中力的傳遞是以與位移和速度的形式表達(dá)力的傳遞, 可以表示為

式中:Sx,Sy分別為礫石傳遞的沖量;vx,vy為位移速度;mi,mj為礫石質(zhì)量。

根據(jù)力的合成原理, 礫石受到的合外力Fi以及合外力矩Ti為:

式中:n為礫石顆粒i相接觸的顆粒j的數(shù)量;fs為顆粒間的滑動摩擦力;ni和iτ 為礫石輪廓線處的法線和切線的單位向量。

力和力矩傳遞表示為

式中:f1,f2表示接觸力;M1,M2表示接觸力矩。

礫石間的不嵌入條件為

式中:Rij表示力學(xué)球心間的距離;di,dj表示兩力學(xué)球的直徑。

2.3 GEM 與DEM 比較

GEM 與DEM 主要有以下幾點(diǎn)不同:

(1) 礫石單元法以礫石為整體剛性單元, 主要以硬球的接觸模式為主, 忽略礫石表面的變形細(xì)節(jié); 離散單元法以粒子形式分布單元, 主要考慮軟球模式的接觸方式, 把顆粒間的法向力簡化為彈簧和阻尼器。

(2) 礫石單元法以堆積礫石為研究對象; 離散單元法可以模擬固體顆粒群的運(yùn)動, 一般不考慮礫石承壓和遮蔽引起的運(yùn)動次序的差別。

(3) 礫石單元法是在堆積礫石穩(wěn)定性研究基礎(chǔ)上提出來的, 在力學(xué)模式上具有外部失穩(wěn)但內(nèi)部穩(wěn)定的物理意義; 離散單元法以模擬固體顆粒群體運(yùn)動的效果為基礎(chǔ), 具有整體失穩(wěn)的物理意義。

3 SPH 與GEM 耦合波浪礫石模型

3.1 SPH 波浪水槽模型及其控制方程

采用SPH 方法, 建立了波浪水槽的數(shù)值模型, 水槽由底部和右側(cè)的固壁邊界及左側(cè)的造波推板邊界組成, 其內(nèi)部填充流體粒子。水槽初始水位為0.4 m,粒子初始間距為0.01 m, 光滑長度取3.2 倍的光滑長度取3.2 倍的粒子初始間距。

SPH 控制方程嚴(yán)格遵循質(zhì)量守恒、動量守恒及能量守恒這三條物理守恒定律, 將描述流體運(yùn)動的Navier-Stokes 方程, 通過SPH 法基本理論的思路對N-S 方程進(jìn)行空間離散化, 得到了在笛卡爾坐標(biāo)系下并且適用于廣義流體動力學(xué)的SPH 控制方程。SPH中最常用的連續(xù)方程粒子近似式為:

式中:m為粒子質(zhì)量;ρ為粒子密度;v為粒子速度;W為粒子核函數(shù); 下標(biāo)i和j為粒子編號;Nj為光滑長度內(nèi)粒子數(shù)量。

根據(jù)以上連續(xù)方程可知, 粒子密度變化率與所求粒子與其相鄰粒子的相對速度有關(guān)。

在SPH 方法中, 帶有層流黏性項(xiàng)的動量守恒方程為:

式中:p為壓強(qiáng);v0為黏性系數(shù);x為坐標(biāo)分量。

在SPH 的計算中, 流體被認(rèn)為是弱可壓縮的。這便于使用狀態(tài)方程來確定流體壓力, 這比求解泊松方程等方程要快得多。SPH 的狀態(tài)方程為:

式中:γ為常數(shù), 在大部分情況下取7; ρ0為參照密度,通常為1 000 kg/m3,B為參數(shù)用于限制密度的最大改變量, 在本文中為在參考密度下的聲速。

圖9 數(shù)值波浪水槽模擬的造波過程和壓力分布Fig. 9 Wave-making process and pressure distribution simulated using the wave numeric channel

圖9 為波浪水槽造波過程和壓力分布, 波浪周期取1.2 s, 輸出數(shù)據(jù)的時間間隔為0.01 s, 總共輸出4.5 s。模型計算由0.75 s 到1 s 的時間間隔內(nèi), 產(chǎn)生第一個大波, 此時刻波峰位置為距造波板0.5 m 左右,并持續(xù)向另一側(cè)邊壁移動。 在模型計算由1 s 至1.25 s的過程中, 造波板經(jīng)由右側(cè)最大擺幅處開始向回運(yùn)動, 此時波峰運(yùn)動到距造波板0.8 m 處。波峰左側(cè)的流體粒子都擁有向右側(cè)運(yùn)動的加速度, 而由于造波板回移導(dǎo)致其附近處水面塌落, 波谷開始形成。在模型計算的1.25 s 至1.5 s 的過程中, 造波板完成一個造波周期, 波峰位置移動到距造波板1.2 m 處, 波谷位置仍處于造波板附近。其后, 在模型計算的1.5 s至3.0 s, 造波板重復(fù)此前的運(yùn)動過程。

建立物理模型試驗(yàn), 對數(shù)值模擬結(jié)果進(jìn)行驗(yàn)證分析。數(shù)值波浪模擬結(jié)果如圖10 所示, 為周期1.2 s的波浪數(shù)值模擬變化與實(shí)測波浪數(shù)值測試。由圖10中可以看出, 當(dāng)數(shù)值波浪模擬由初始時刻開始到達(dá)穩(wěn)定狀態(tài)后, 數(shù)值模擬與實(shí)測波浪模擬結(jié)果較吻合。

3.2 SPH 與GEM 耦合模型數(shù)學(xué)基礎(chǔ)

SPH 與GEM 耦合模型由SPH 粒子和堆積礫石構(gòu)成, SPH 粒子運(yùn)動由方程(12)控制, 堆積礫石處為SPH 粒子的固體邊界, SPH 粒子運(yùn)動速度為

圖10 波浪水槽數(shù)值模擬與實(shí)測數(shù)據(jù)對比Fig. 10 Numerical simulation of the wave channel compared with the measured data

式中,Г表示礫石表面;?表示礫石內(nèi)部;n表示礫石表面法線方向。

SPH 粒子對礫石的作用力為沿礫石表面的法線方向壓力的積分, 利用勢流的伯努利方程為

式中,p∞表示礫石表面;vθ表示礫石表面切線方向速度;n表示礫石表面法線方向。

礫石運(yùn)動軌跡方程為

式中,vx0,vy0,ω0分別為礫石在計算時段初始t0的運(yùn)動速度和角速度; ∑Fx, ∑Fy,∑M表示礫石所受的合力和合力矩;J表示礫石的轉(zhuǎn)動慣量。

礫石間力的傳遞由公式(7)至公式(10)確定。

3.3 SPH 與GEM 耦合模型模擬過程

SPH 與GEM 的耦合過程主要有以下幾點(diǎn):

(1) 應(yīng)用壩體輪廓線控制將具有級配的堆積礫石隨機(jī)放置在SPH 數(shù)值水槽中, 確定每一塊礫石輪廓線內(nèi)受控制的SPH 粒子, 這些粒子稱為受控粒子,不包含SPH 粒子的礫石是水面以上礫石或尺度較小的礫石, 不包含SPH 粒子的礫石不會影響礫石運(yùn)動計算模擬。受控粒子與非受控粒子分布見圖11。

(2) 在SPH 模型邊界上使用推板造波, 礫石中的受控SPH 粒子參與核函數(shù)計算, 但運(yùn)動速度為零,相當(dāng)于將礫石對周邊SPH 粒子的影響分解為受控粒子對周邊SPH 粒子的影響, 體現(xiàn)了礫石邊界輪廓對SPH 模式的耦合效果。

圖11 受控粒子與非受控粒子分布Fig. 11 Distribution of controlled and uncontrolled particles

(3) 同步統(tǒng)計礫石的受力和周邊的流速場影響(見圖12), 礫石受到的總壓力可以用礫石中受控粒子的總壓力代替, 受控粒子(圖12 黑輪廓線內(nèi)粒子)的壓力分布決定了礫石運(yùn)動形態(tài)。礫石周邊流場取礫石輪廓以外、礫石光滑長度以內(nèi)的粒子做統(tǒng)計分析, 圖12 中黑色礫石輪廓外與紅色光滑長度圓內(nèi)的粒子, 黑色矢量表示礫石光滑長度, 本文取2 倍的力學(xué)球半徑, 這一區(qū)域的流速、流向可以處理為波浪對礫石的沖擊效果, 礫石是否進(jìn)入運(yùn)動狀態(tài), 這也是重要判別條件之一。

(4) 在堆積礫石的受力分析中, 用擬序排列、虛位移原理確定礫石的承壓力、反力、摩擦力和擠壓力, 確定礫石進(jìn)入運(yùn)動狀態(tài)的級別, 通常級別的確定會在每次計算前重新劃分, 因?yàn)槊恳淮蔚[石運(yùn)動都會改變礫石堆的排列形式, 一般一次計算中只需考慮1 級礫石的運(yùn)動, 其他級別得不到運(yùn)動的機(jī)會。

(5) 波浪力不足以撼動礫石時, 礫石保持原有狀態(tài); 波浪力撼動礫石時, 要考慮礫石的反力分布確定礫石運(yùn)動形式, 一般情況為沒有反力作用時為塌落, 有支撐點(diǎn)時要找出支撐點(diǎn)礫石。

(6) 礫石運(yùn)動和就位應(yīng)在瞬時完成, 所以在一個時間步長內(nèi)就應(yīng)當(dāng)確定出礫石運(yùn)動后的位置, 這一位置應(yīng)避免與其它礫石壓入, 并位于礫石運(yùn)動的路徑上;礫石進(jìn)入新的位置后將會放棄以前的受制粒子, 同時產(chǎn)生出新的受制粒子, 完成了一個計算的循環(huán)過程。

3.4 SPH 與GEM 耦合模型模擬結(jié)果

應(yīng)用SPH-GEM 方法耦合的數(shù)學(xué)模型來模擬波浪作用下防波堤模型由靜態(tài)平衡到動態(tài)平衡的變化過程。防波堤模型在數(shù)值波浪的作用發(fā)生改變, 邊坡部分計算顆?；? 但流體粒子仍正常穿越礫石邊界之間的孔隙參加運(yùn)動。數(shù)值模擬波浪水槽水深為0.4 m,波高為0.10 m。防波堤模型選取坡度為1︰1.5, 肩臺高0.4 m, 堤頂高0.6 m。圖13 展示了計算模型隨時間的變化。當(dāng)模型計算由初始0 時刻到3 s 的過程中, 造波程序生成波浪并傳遞至防波堤模型, 并開始于防波堤模型進(jìn)行相互作用。當(dāng)模型由3 s 至4.5 s 的計算過程中, 由于波浪作用導(dǎo)致防波堤模型的邊坡粒子部分滑落并于坡面下部堆積。當(dāng)模型計算由4.5 s 至5 s 的計算過程中, 防波堤模型達(dá)到穩(wěn)定狀態(tài)。由圖可見在計算時刻為5 s 的防波堤模型, 其斷面形狀已經(jīng)達(dá)到反S 的形狀。

圖13 防波堤模型坡面顆粒變化過程Fig. 13 Process of changing the grain slope in the breakwater model

圖14 為防波堤模型初始狀態(tài), 圖15 為防波堤模型穩(wěn)定狀態(tài), 圖16 為最終斷面平衡對比圖。防波堤模型礫石的最終動力平衡斷面為反S 形狀。

3.5 耦合模型的可靠性驗(yàn)證

為了進(jìn)一步驗(yàn)證數(shù)值模型結(jié)果的可靠性, 采用本文數(shù)值模擬方法對文獻(xiàn)[17]中的物模結(jié)果進(jìn)行驗(yàn)證分析。模型選取的水槽長度為5 m, 水槽高度為0.8 m, 對應(yīng)的水深為15 cm, 波高為5 cm, 周期選擇1.5 s。以SPH 方法生成的初始流體粒子間距為0.01 m。防波堤物理模型的對應(yīng)位置及測點(diǎn)位置如圖17 所示。防波堤模型選取粒徑為6 cm 左右的單一介質(zhì), 坡面角度為1︰1.5, 堤頂高30 cm, 堤頂寬10 cm, 模型的孔隙率約為0.35。圖18 為數(shù)值模型的計算示意圖。根據(jù)文獻(xiàn)中的物理模型實(shí)驗(yàn)結(jié)果, 分別對4 個測點(diǎn)位置的波高數(shù)值進(jìn)行驗(yàn)證分析。

對比模型的計算數(shù)據(jù)與文獻(xiàn)中的試驗(yàn)數(shù)據(jù)發(fā)現(xiàn)波高數(shù)據(jù)基本吻合, 見圖19。但是波浪模型的計算仍然存在誤差, 考慮由于水深和波高都較小導(dǎo)致流體粒子的運(yùn)算受到影響, 并且由于流體粒子的黏性系數(shù)考慮波高不能與實(shí)驗(yàn)數(shù)據(jù)完全吻合。

圖14 防波堤模型初始狀態(tài)Fig. 14 Initial state of the breakwater model

圖15 防波堤模型穩(wěn)定狀態(tài)Fig. 15 Steady state of the breakwater model

圖16 防波堤模型斷面變化對比Fig. 16 Comparison of the cross section of the breakwater model

4 波浪礫石物理模型試驗(yàn)

4.1 物理模型試驗(yàn)

物理模型試驗(yàn)在天津大學(xué)北洋園校區(qū)港口航道與海岸工程實(shí)驗(yàn)室的小水槽進(jìn)行。試驗(yàn)水槽長35 m,寬1 m, 高1 m, 模型距離造波板28 m, 造波機(jī)為天津理工大學(xué)制造的AFM-124 型造波機(jī), 該造波機(jī)可根據(jù)試驗(yàn)條件需求產(chǎn)生規(guī)則波與JONSWAP 譜的不規(guī)則波, 由計算機(jī)系統(tǒng)控制運(yùn)行。試驗(yàn)波浪參數(shù)的測量與分析, 由中國水電科學(xué)研究院制造的波高儀、壓力傳感器和SG-800 型水工模型試驗(yàn)數(shù)據(jù)采集與處理系統(tǒng)完成, 流速測量采用聲學(xué)多普勒流速儀Vectrino-II,試驗(yàn)過程由計算機(jī)進(jìn)行數(shù)據(jù)輸入、輸出、和處理。波浪試驗(yàn)水槽見圖20。

圖17 文獻(xiàn)實(shí)驗(yàn)布局示意圖(單位: cm)Fig. 17 Literature experiment layout

圖18 數(shù)值模型計算示意圖Fig. 18 Schematic of the numerical model calculation

試驗(yàn)斷面和設(shè)備分布如圖21, 堤頂高程為50 cm,堤頂寬度為20 cm, 肩臺高程為42 cm, 肩臺寬度為25 cm, 迎波堤面坡度為1︰1.5, 肩臺處及迎波面礫石護(hù)面層厚度為10 cm, 堤頂處礫石厚度為8 cm, 背波堤面坡度為1︰1.5, 背波堤面礫石護(hù)面層厚度為5 cm。

圖20 波浪試驗(yàn)水槽Fig. 20 Experimental wave tank

圖21 試驗(yàn)設(shè)計斷面和試驗(yàn)設(shè)備在水槽中的分布Fig. 21 Test design and distribution section of the test equipment in the water tank

4.2 物理模型試驗(yàn)結(jié)果

圖22 為在1 號點(diǎn)處物理模型的波浪觀測值與數(shù)學(xué)模型模擬值的比較?？梢钥闯鑫锢砟Ｐ驮觳ㄆ椒€(wěn)、波高變化較小, 數(shù)學(xué)模型造波周期穩(wěn)定, 但波高起伏較大。因數(shù)學(xué)模型選取的粒子有限, 在波高處粒子相對少一些, 會影響波峰模擬的精度。但總變化趨勢相同, 平均波高接近, 說明數(shù)值波浪水槽模擬成果基本正確。

圖23 為在極端水位+ 40 cm 高程, 周期1.2 s, 波高H=10 cm 的規(guī)則波作用下, 波浪對坡面的沖刷區(qū)間開始上移, 整個+42 cm 高程肩臺都在沖刷范圍內(nèi)。+42 cm 高程肩臺在波浪的作用下迅速被沖刷并坍落。在持續(xù)波浪作用下, 肩臺被沖刷的范圍相比設(shè)計低水位及設(shè)計高水位有明顯變化。當(dāng)波浪持續(xù)作用200 個波后, 斷面模型達(dá)到最終動力平衡狀態(tài)。整個+42 cm 高程肩臺被沖刷掉約7 cm 左右。坡面最大沖刷厚度約為4 cm, 沖刷滑落的塊石在坡面下部堆積約4 cm, 部分塊石滑落至坡腳并堆出近3 cm 左右的寬度。整個沖刷斷面呈明顯反S 形, 與數(shù)學(xué)模型結(jié)果變化趨勢一致。

圖22 波浪觀察與波浪模擬的過程比較Fig. 22 Comparison between wave observation and wave simulation

圖23 防波堤斷面1 模型試驗(yàn)工況1 動力平衡斷面Fig. 23 Section 1 of the breakwater model test Condition 1 dynamic balance section

5 結(jié)論

(1) 采用傅氏級數(shù)的方法對塊石輪廓進(jìn)行了曲線展開, 說明塊石輪廓沿曲線的矢徑分布情況, 通過分布礫石輪廓上力學(xué)性質(zhì), 可以對礫石輪廓賦予力學(xué)上的摩擦系數(shù)、轉(zhuǎn)動慣量等物理量, 提出了力學(xué)球體的概念。

(2) 分析了堆積礫石的受力特點(diǎn), 提出堆積礫石分級方法和擬序排列求解方法, 對靜態(tài)礫石堆應(yīng)用擬序排列求得的解, 與理論值誤差很小, 說明擬序排列求解方法可行。

(3) 提出了GEM 模型方法, 建立了SPH 波浪水槽與GEM 堆積礫石耦合數(shù)學(xué)模型, 模擬了礫石在波浪作用下的堤壩變形過程, 模擬結(jié)果與物理模型試驗(yàn)結(jié)果基本一致。