礫石單元法與SPH 耦合模型在數值波浪水槽中的應用研究

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

(1. 天津大學 水利工程仿真與安全國家重點實驗室, 天津 300072; 2. 國家海洋技術中心, 天津 300112)

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

1 堆積礫石的受力分析

1.1 單礫石幾何形態

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

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

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

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

傅里葉系數由以下公式確定

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

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

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

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

形狀系數公式為[14]

對退化的二維礫石顆粒, 可以假定c=b。對二維礫石顆粒, 引用公式(1)中傅里葉級數各項系數an、bn作為傅里葉級數各級礫石的長軸和中長軸, 可以將非圓顆粒形狀用數個圓顆粒部分表面的組合形狀來近似(圖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 單礫石力學球概念

通過考察不同投影平面的礫石傅立葉級數各級系數半徑, 可以得到礫石的球度系數或形狀系數,不僅從礫石表面整體上獲取礫石對球的近似程度,而且可以對不同位置的球面近似形態進行描述。可以將泥沙顆粒的球度系數延伸到礫石的球態研究中,表1 為礫石幾何形狀與球度系數的關系。

表1 礫石幾何形狀與球度系數的關系Tab. 1 Relationship between gravel geometry and sphericity coefficient

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

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

力學球的滑動或滾動摩擦的阻力系數可以表示為

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

力學球的轉動慣量可以表示為

式中:Js為非球面平均轉動慣量;Fs(θ)為隨θ變化的隨機函數。

1.3 堆積礫石的受力狀態

堆積礫石表層石塊所受到的作用力會通過接觸壓力、反力、碰撞傳遞給其它礫石, 使礫石進入運動狀態或達到新的平衡狀態。堆積礫石受力見圖4, 力學球礫石堆見圖5。

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

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

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

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

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

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

2.2 GEM 主要力學公式

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

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

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

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

力和力矩傳遞表示為

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

礫石間的不嵌入條件為

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

2.3 GEM 與DEM 比較

GEM 與DEM 主要有以下幾點不同:

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

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

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

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

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

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

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

式中:m為粒子質量;ρ為粒子密度;v為粒子速度;W為粒子核函數; 下標i和j為粒子編號;Nj為光滑長度內粒子數量。

根據以上連續方程可知, 粒子密度變化率與所求粒子與其相鄰粒子的相對速度有關。

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

式中:p為壓強;v0為黏性系數;x為坐標分量。

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

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

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

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

建立物理模型試驗, 對數值模擬結果進行驗證分析。數值波浪模擬結果如圖10 所示, 為周期1.2 s的波浪數值模擬變化與實測波浪數值測試。由圖10中可以看出, 當數值波浪模擬由初始時刻開始到達穩定狀態后, 數值模擬與實測波浪模擬結果較吻合。

3.2 SPH 與GEM 耦合模型數學基礎

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

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

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

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

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

礫石運動軌跡方程為

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

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

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

SPH 與GEM 的耦合過程主要有以下幾點:

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

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

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

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

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

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

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

3.4 SPH 與GEM 耦合模型模擬結果

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

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

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

3.5 耦合模型的可靠性驗證

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

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

圖14 防波堤模型初始狀態Fig. 14 Initial state of the breakwater model

圖15 防波堤模型穩定狀態Fig. 15 Steady state of the breakwater model

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

4 波浪礫石物理模型試驗

4.1 物理模型試驗

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

圖17 文獻實驗布局示意圖(單位: cm)Fig. 17 Literature experiment layout

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

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

圖20 波浪試驗水槽Fig. 20 Experimental wave tank

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

4.2 物理模型試驗結果

圖22 為在1 號點處物理模型的波浪觀測值與數學模型模擬值的比較。可以看出物理模型造波平穩、波高變化較小, 數學模型造波周期穩定, 但波高起伏較大。因數學模型選取的粒子有限, 在波高處粒子相對少一些, 會影響波峰模擬的精度。但總變化趨勢相同, 平均波高接近, 說明數值波浪水槽模擬成果基本正確。

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

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

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

5 結論

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

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

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