楊 杰，張義佐，黃 錚，馬夢杰，劉延波

(1.中交上海航道局有限公司，上海 200002；2.天津天科工程管理有限公司，天津 300450；3.哈爾濱工程大學 船舶工程學院，黑龍江省 哈爾濱 150001)

拋石作業(yè)被廣泛應用于我國沿海沿江地區(qū)的河道治理、護岸、防汛搶險、橋墩防沖等工程中。在拋石施工過程中，要重點把握好“拋準、拋足、拋勻”這3個指標。其中“拋準”是拋石作業(yè)中最為關(guān)鍵的評價指標，即應該準確地將拋石投在規(guī)定的區(qū)域范圍內(nèi)，完全覆蓋塌陷區(qū)域[1]。在實際的拋石施工中，多采用網(wǎng)兜、液壓反鏟式挖掘機和抓斗等方式進行機械拋填，在拋石量較大時，還會采用底開駁船、側(cè)拋船等特殊作業(yè)船舶進行直接拋投。在實際工程中，拋石多為混合石料群拋，而非單顆粒拋石，因此有必要對多粒徑塊石群拋后的沉降軌跡及落點位置進行研究。

1 研究背景

拋石漂移距是指拋投的塊石從進入水中落到設計位置，中間因為水流作用而產(chǎn)生橫向或者縱向水平移動的距離。若拋投位置不準確，可能會造成拋石超過預期的拋投范圍，不能達到設計要求的相關(guān)高程，最終會大幅減弱治理的效果；也可能會導致塊石拋投不均勻，主要維護的地段會出現(xiàn)大片空白，甚至形成工程薄弱區(qū)域。這對于工程整體的作用是一種削弱，甚至可能會影響到工程的安全性。

為探究水下拋石施工過程中的塊石漂移過程，為工程提出指導意義，梁潤[2]、詹義正[3]、姚仕明等[4]推導出拋石漂移距的相關(guān)公式，但尚未找到適用于多粒徑水上塊石群拋的漂距估算公式，且經(jīng)驗公式得出的塊石漂距往往與實測值偏差較大。韓海騫等[5]、李小超等[6]設計相關(guān)的水槽試驗，并在施工標段做了現(xiàn)場試驗，但試驗需要不斷調(diào)整，成本較高。相比之下，數(shù)值模擬成本低廉，只需要在計算機上進行模擬和數(shù)據(jù)處理，在不能進行試驗的惡劣天氣下也能采用數(shù)值模擬估算出漂距。

拋石運動過程是一個典型的流-固耦合現(xiàn)象，塊石顆粒為固體相，水流為流體相，為精準模擬這一流固耦合過程，張強強[7]、陳凱華等[8]、劉卡等[9]、李衛(wèi)國[10]均建立了數(shù)值模型對水下拋石的運動過程進行模擬，他們的研究證明數(shù)值模擬的方法可以對拋石工程達到良好的模擬效果，得到的預測結(jié)果較為可靠，可以指導相關(guān)的施工作業(yè)。但現(xiàn)存數(shù)值研究多集中于對單顆粒塊石漂距的研究，針對多粒徑混合塊石群拋擴散過程的研究較少。為了更好地探究塊石粒徑對水下拋石的漂移距離與沉降時間的關(guān)系，本文采用計算流體力學-離散元(CFD-DEM)方法建立針對水下拋石作業(yè)的耦合分析模型，并以塊石粒徑與來流速度為變量對水下拋石的漂移距離與沉降時間進行研究。

在針對拋石施工過程的數(shù)值仿真中，為模擬塊石碰撞中的接觸、壓縮與回彈過程，采用離散元方法(DEM)對塊石碰撞的實際情況進行模擬。在DEM計算過程中，采用牛頓第二定律和接觸力-位移方程完成塊石漂移運動和塊石碰撞的模擬。水流的運動變化情況則可采用經(jīng)典的計算流體力學方法(CFD)進行求解。最早在1992年Tsuji等[11]提出了CFD-DEM耦合的方法，基于歐拉-拉格朗日顆粒軌道模型，將流體相和顆粒相分別用計算流體力學和離散元法處理計算，用于研究一維水平管的無黏性節(jié)涌流化過程。近年來，隨著CFD和DEM不斷發(fā)展完善，CFD-DEM耦合方法在研究固液兩相流方面也具有更強的優(yōu)勢，能夠滿足大部分條件下的固-液兩相流模擬。因此本文將計算流體力學程序OpenFOAM和離散元程序LIGGGHTS進行耦合，建立CFD-DEM耦合模型研究拋石群的沉降過程及河床分布情況。

2 數(shù)值模型建立

2.1 離散單元法理論

在CFD-DEM耦合方法中，流體被看作是含顆粒的連續(xù)介質(zhì)，采用含體積分數(shù)的Navier-Stokes方程描述；顆粒被看作是離散單元，通過檢測單元間接觸量，代入剛度模型由牛頓第二定律描述，流相-離散相間耦合通過求解每一迭代步中流場對每個顆粒的動量與能量傳遞實現(xiàn)。

在CFD-DEM耦合方法中的離散單元部分，顆粒之間的運動遵循牛頓第二定律，并在動力守恒方程中考慮了流體對顆粒的作用力，在給定的時間步長中，通過顯式迭代顆粒相關(guān)位置進行求解。在DEM部分顆粒平動與轉(zhuǎn)動控制方程分別如下：

(1)

(2)

式中：m為固體質(zhì)量；t為時間；δ為平動位移；Fe為外力；Fc為接觸力；Ff為流體施加在固體上的力；c為阻尼系數(shù)；IM為慣性矩張量；θ為角位移；M為總扭矩，包括外力、流固相互作用力的扭矩。

在DEM系統(tǒng)中，顆粒間的運動模擬通過考慮顆粒間相互接觸量與接觸處的本構(gòu)關(guān)系實現(xiàn)。顆粒接觸本構(gòu)關(guān)系又被稱為接觸剛度，考慮到顆粒間的平動與轉(zhuǎn)動，顆粒的接觸剛度通常被分解為接觸法向割線剛度與接觸切向切線剛度，并由不同的接觸剛度模型給出。顆粒的法向接觸剛度與切向剛度定義式為：

(3)

為模擬塊石碰撞時的連續(xù)接觸，本文使用Hertz-Mindlin無滑移“軟球”模型，作為DEM部分的接觸剛度模型以計算顆粒的法向與切向的接觸剛度。Hertz正應力和切應力的控制方程為：

(4)

式中：δ為兩顆粒接觸時的重疊范圍，其值為D-d，其中D為直徑，d為兩顆粒之間距離；ri、rj分別為顆粒i、j的半徑；meff為顆粒i、j的等效質(zhì)量，其值為mimj/(mi+mj)，其中mi、mj分別為顆粒i、j的質(zhì)量；kn、kt分別為法向、切向彈性系數(shù)；γn、γt分別為法向、切向黏彈性阻力系數(shù)；Δst為兩顆粒之間的切向位移矢量；nij為沿著連接兩個顆粒中心線的單位矢量；vs，n、vs，t為兩顆粒相對速度的法向、切向分量。

2.2 計算流體力學理論

在CFD部分，由于顆粒占據(jù)了相應的流體空間，因此需要引入含體積分數(shù)的Navi-Stokes方程質(zhì)量守恒式描述含顆粒流體的運動。由于需要計算流體與顆粒間的相互交換動量，須在Navi-Stokes方程動量守恒式末尾加入流體交換的顆粒動量項，使用CFDEM程序自身的耦合求解器進行計算求解。本文使用cfdemSolverPiso求解器，該求解器采用含體積分數(shù)的Navi-Stokes守恒式為：

(5)

式中：αl為在單個計算網(wǎng)格內(nèi)的流體含量；ρl為流體密度；ul、us分別為流體、固體速度；p為壓強；Ksl為固液相間的隱含動量交換系數(shù)；τ為流體剪切應力；g為重力加速度；f為固相交換至液相的動量。

2.3 CFD-DEM耦合理論

(6)

式中：fd，i為流體與顆粒相對運動產(chǎn)生的拖曳力；f?p，i為由流體壓力梯度產(chǎn)生的壓力梯度力；f?·τ，i為黏性流體在顆粒表面產(chǎn)生的相應剪切應力。

近年來，針對流體與顆粒相對運動產(chǎn)生的拖曳力fd，i，許多學者采用模型試驗與數(shù)值方法等手段提出許多計算模型。其中Hill、Koch等利用格子玻爾茲曼方法研究流體與固體的相互作用，進一步提出包含顆粒雷諾數(shù)Rep和顆粒體積分數(shù)φ的拖曳力分段計算公式。該拖曳力模型后由Benyahia拓展到整個體積分數(shù)、顆粒雷諾數(shù)空間上[12-13]。其中fd，i可由式(7)表達，f?p，j與f?·τ，i可以由式(8)表示：

(7)

(8)

式中：εs為顆粒體積分數(shù)；Rep為顆粒雷諾數(shù)；F0(εs)、F1(εs)與F3(εs)是由εs與Rep確定的相應拖曳力分量；Vp為顆粒體積；p為顆粒所受壓力；τ′為流體黏度。

流體與顆粒的動量交換模型包括在顆粒邊界上對流體作用力積分的解析模型和采用經(jīng)驗拖曳力模型的非解析模型。由于引入了大量顆粒模擬堆石漂移群的運動，因此本文采用PISO算法以求解非解析模型下流體施加于顆粒上的經(jīng)驗拖曳力以提高求解效率。PISO(pressure implicit split operator，求解壓力的隱式算子分裂)算法是典型的2步校正算法，主要實施步驟包括預估步、第一校正步、第二校正步。PISO算法具有相鄰校正，即在每個迭代步中，相鄰網(wǎng)格的速度值都會采用最新的速度預測值，其具體計算步驟為：1)假定壓強，隱式求解運動方程，得到速度預測值u*，即速度預測；2)利用u*求解連續(xù)方程，得到壓強值p*，即壓強求解；3)利用p*再求解運動方程，但改為顯示求解，得到u**，即速度校正；4)回到第2步，迭代求解，直至允許的迭代公差。

2.4 模型參數(shù)設置

在CFD工況設計部分中，流場尺寸設定為15 m×10 m×15 m(長×寬×高)，入口設置為沿x正方向速度分別為1.00、1.25、1.50、2.00 m/s的入口，出口設置為與入口流速相同。迭代時間步為0.5 ms，并采用PISO算法對流場的速度與壓力進行二次修正。網(wǎng)格劃分采用結(jié)構(gòu)化六面體網(wǎng)格進行，網(wǎng)格尺寸根據(jù)實際塊石粒徑確定。根據(jù)不同網(wǎng)格尺寸的試算結(jié)果，網(wǎng)格尺寸與塊石顆粒直徑的比值確定為2～3。

在DEM工況設計部分中，在不同粒徑顆粒群拋投仿真中，設置顆粒材料的泊松比為0.45，彈性模量為50 GPa，法向恢復系數(shù)為0.3，顆粒轉(zhuǎn)動阻尼為0.5。仿真中設置的拋石群顆粒粒徑級配參數(shù)見表1。

表1 拋石群顆粒級配參數(shù)

3 拋石下落規(guī)律分析

3.1 不同粒徑顆粒下落位置隨時間的變化

以拋石粒徑為單一變量，通過對5種粒徑分別為0.075、0.110、0.160、0.260、0.300 m的顆粒在水中運動過程的模擬，得到了顆粒的運動軌跡，其沉降時間-水平移距關(guān)系曲線如圖1所示?？梢钥闯觯谄渌麠l件不變的情況下，改變顆粒的粒徑對水平移距的影響相對較大。粒徑越小者，水平漂移距離越遠，運動停止也越慢，而當塊石粒徑大于等于0.26 m時，塊石在下落過程中漂移現(xiàn)象并不顯著。

圖1 不同粒徑顆粒水平移距隨時間變化

3.2 不同流速下不同粒徑顆粒沉降時間的變化

分別以顆粒粒徑為0.075、0.110、0.160、0.260、0.300 m，水流流速為1.00、1.25、1.50、2.00 m/s進行模擬。5種不同粒徑顆粒在4種水流流速下的沉降時間變化如圖2所示。可以看出，若控制其他變量，只改變水流的初始速度，則相同粒徑的顆粒隨流速的增加，其沉降時間越久；而在相同的流速條件下，顆粒的粒徑越大，則沉降時間越短。

圖2 不同流速下5種粒徑顆粒的沉降時間變化

3.3 不同粒徑拋石群擴散過程及流場分布

在CFD-DEM耦合算法中，一定迭代步數(shù)后會更新顆粒位置以更新流場網(wǎng)格中的體積分數(shù)，流場的流線也隨著顆粒的運動發(fā)生變化。由于顆粒間相互碰撞與摩擦，初始階段分布在生成區(qū)域前端的顆粒由于后部顆粒的碰撞會在拋投過程中不斷向前擴散，而后部顆粒由于在碰撞中將動量傳遞至前部顆粒，漂距相較于前部顆粒較小。在2.0 m/s流速下，不同粒徑顆粒組成的拋石群拋投后在水下的擴散過程及流場形態(tài)如圖3所示?？梢悦黠@看出粒徑大的顆粒首先開始沉降且水平漂距較近，粒徑較小顆粒沉降時間較久并且水平漂距相對較遠，同時會產(chǎn)生一定的堆積現(xiàn)象。

圖3 拋石群拋投后顆粒擴散過程及流場側(cè)視圖

3.4 不同流速下拋石群的落點分布

模擬不同粒徑塊石在水流條件分別為1.00、1.25、1.50及2.00 m/s下進行群拋，得到的塊石顆粒落至底床后的分布形態(tài)如圖4所示?？梢钥闯觯煌降膲K石顆粒群拋投后落至床底后的分布形狀較為不規(guī)則。在流速相同的條件下，粒徑越小者漂移的距離越遠，而粒徑大者落地距離較近。通過圖4a)～d)對比可看出，流速大的情況下各個粒徑顆粒的漂距都較遠，流速小時則漂距相對較近。

圖4 不同流速下不同粒徑顆粒落點分布

4 結(jié)論

1)在同一流場中，多粒徑塊石顆粒群拋過程中，塊石落點分布不規(guī)則。仿真結(jié)果表明，由于前部顆粒與后部顆粒之間的相互碰撞以及塊石粒徑差異造成的重力分量不同，塊石的漂距表現(xiàn)為由最小漂距與最大漂距間的隨機分布。

2)在相同來流條件下，多粒徑塊石群中的粗粒徑塊石漂距顯著小于細粒徑塊石，粒徑大于等于0.26 m的塊石在下落過程中漂移現(xiàn)象不顯著。

3)相同的流速條件下，由于塊石運動過程中塊石所受外力中重力分量隨塊石粒徑增加而增加，塊石沉降時間隨塊石粒徑增加而顯著增加。因此，在實際工程中應結(jié)合相關(guān)規(guī)范，盡量選取粗粒徑塊石，以提高塊石群拋的施工精度。