張 偉，王春磊，李衛(wèi)國(guó)，王瀚銳，李小超，常留紅

(1.中交天航港灣建設(shè)工程有限公司，天津 300450；2.湖南省水運(yùn)建設(shè)投資集團(tuán)有限公司，長(zhǎng)沙 410011；3.長(zhǎng)沙理工大學(xué) 水利工程學(xué)院，長(zhǎng)沙 410114)

拋石群水下運(yùn)動(dòng)擴(kuò)散是河道治理[1]、護(hù)岸、潛丁壩、橋墩防沖等實(shí)際工程中經(jīng)常遇到的問題之一[2]。拋石施工受環(huán)境因素影響較明顯，既有一定的規(guī)律性，又呈現(xiàn)出一定的隨機(jī)性，一直是工程界關(guān)注的熱點(diǎn)問題之一，目前尚難找到適應(yīng)于各種水下拋石的施工工藝。

以往拋石運(yùn)動(dòng)的研究主要關(guān)注相關(guān)經(jīng)驗(yàn)公式，如早期蘇聯(lián)學(xué)者伊茲巴斯通過平堵截流模型試驗(yàn)提出目前工程最常用的伊茲巴斯公式[3]、考慮水深的沙漠夫公式[4]以及姚仕明[5]、毛佩郁[6]、韓海騫[7]、姚仕明[8]等學(xué)者推導(dǎo)的拋石漂移的相關(guān)公式。此外，還有不少學(xué)者研究了拋石水下的運(yùn)動(dòng)過程。陳凱華等[9]分析了水下拋石運(yùn)動(dòng)規(guī)律以及拋石基床的堆積情況。劉曉青[10]等通過物理模型試驗(yàn)研究散拋塊石種類、尺寸及水流速度對(duì)塊石漂移量的影響。張宏千[11]等通過對(duì)水流運(yùn)動(dòng)中的透水框架受力分析，推導(dǎo)了不同搭接和拋投方式的透水框架在不同拋投高度的落距公式。劉卡等[12]基于計(jì)算流體力學(xué)與離散元法耦合的方法，從拋石初始速度和拋石粒徑等方面模擬了拋石水下運(yùn)動(dòng)過程。王茂枚等[13]通過河床結(jié)構(gòu)的異質(zhì)性,從粒子隨機(jī)運(yùn)動(dòng)角度解釋了拋石在河道中反常輸運(yùn)的物理機(jī)理。李小超等[14]通過開展模型試驗(yàn)，分析了不同因素對(duì)拋石漂移距離的影響，并導(dǎo)出了相應(yīng)的計(jì)算公式。還有學(xué)者通過開展現(xiàn)場(chǎng)試驗(yàn)，探究漲、落潮條件下拋石水下運(yùn)動(dòng)的影響因素，并推導(dǎo)拋石漂移距離公式。但鮮有關(guān)于拋石群水下擴(kuò)散規(guī)律的研究。

拋石群水下運(yùn)動(dòng)擴(kuò)散是影響拋石成型及其穩(wěn)定的重要因素。依托長(zhǎng)江干線武漢至安慶段6 m水深航道整治工程，基于OpenFOAM建立拋石水下運(yùn)動(dòng)三維流固耦合數(shù)學(xué)模型，研究不同粒徑組合方式拋石群水流方向漂移距離、拋石群落點(diǎn)分布，揭示拋石群水下運(yùn)動(dòng)擴(kuò)散規(guī)律。對(duì)水下拋石質(zhì)量控制和成本控制具有重要的理論和現(xiàn)實(shí)意義，研究成果可在類似工程推廣應(yīng)用。

1 工程概況

工程河段上起小孤山，下至華陽河口，長(zhǎng)30 km，為分汊河型，南汊為主汊，較順直，為主航道所在。北汊為支汊，又稱為馬當(dāng)圓水道，分流比很小，河道彎曲、狹窄，進(jìn)口已淤淺至航行基面以上。南汊分為3個(gè)水道：馬當(dāng)磯以上為馬當(dāng)南水道，江心有心灘——棉外洲；馬當(dāng)磯—娘娘廟為馬當(dāng)阻塞線水道(以下簡(jiǎn)稱馬阻水道)，被瓜子號(hào)洲再度分為兩汊，右汊為主汊及主航道；娘娘廟以下為東流直水道。馬當(dāng)南水道右岸為江西省彭澤縣，左岸為江西省彭澤縣棉船鎮(zhèn)，上起小孤山，下至馬當(dāng)磯，長(zhǎng)約13 km，為受節(jié)點(diǎn)控制的二級(jí)分汊河段，江中棉外洲分水道為左右兩槽，右槽為主航槽較為窄深，左槽相對(duì)寬淺。馬當(dāng)南水道上接?xùn)|北水道，下連馬阻水道。東流水道位于長(zhǎng)江下游九江—安慶河段，上起華陽河口，與馬當(dāng)水道相連，下迄吉陽磯，與官洲水道相接，全長(zhǎng)約31 km。工程的主要建設(shè)內(nèi)容包括馬當(dāng)河段棉外洲順壩加高工程、棉外洲中部護(hù)灘工程、左槽中段潛壩工程、瓜子號(hào)洲右緣護(hù)灘帶工程以及疏浚工程和東流水道天玉串溝控制工程。工程主要工作量包括：棉外洲順壩加高工程拋石50 570 m3；左槽中段潛壩工程拋石280 646 m3；瓜子號(hào)洲右側(cè)中部護(hù)灘帶工程8#護(hù)灘帶拋石89 002 m3；瓜子號(hào)洲右側(cè)中部護(hù)灘帶工程9#護(hù)灘帶拋石80 544 m3；瓜子號(hào)洲右側(cè)中部護(hù)灘帶工程10#護(hù)灘帶拋石100 808 m3；天玉串溝控制工程1#護(hù)底帶拋石69 720 m3；天玉串溝控制工程2#護(hù)底帶拋石70 948 m3。

根據(jù)實(shí)測(cè)水文資料，本工程所處水域大多為主流頂沖水域，實(shí)測(cè)流速一般在1.7 m/s左右，最大為1.98 m/s。施工過程中拋石建筑物附近的局部水流結(jié)構(gòu)不斷被改變，如果斷面拋投成型質(zhì)量不能保證，在高流速、大水深條件下很容易出現(xiàn)穩(wěn)定性不足、發(fā)生水毀等問題?；旌鲜蠏佂冻尚瓦^程影響因素復(fù)雜，拋石群體下落過程的漂移和分層等現(xiàn)象使混合石料拋投斷面形態(tài)及級(jí)配組成隨之變化，直接影響斷面穩(wěn)定性，因此有必要對(duì)混合石料群體拋石水下運(yùn)動(dòng)過程及落至底床上的分布形態(tài)進(jìn)行研究，為拋投施工穩(wěn)定性分析提供依據(jù)。本文采用CFD-DEM耦合方法對(duì)群體拋石水下運(yùn)動(dòng)過程進(jìn)行模擬，探討拋石群落至底床上的分布形態(tài)和級(jí)配組成。

2 數(shù)值模型的建立

2.1 流體運(yùn)動(dòng)控制方程

流體運(yùn)動(dòng)的控制方程為局部平均的不可壓縮Navier-Stokes方程

(1)

(2)

2.2 顆粒運(yùn)動(dòng)控制方程

固體運(yùn)動(dòng)采用牛頓第二定律進(jìn)行求解，其控制方程如下

(3)

(4)

2.3 固體顆粒間的作用控制方程

考慮塊石之間的接觸具有較強(qiáng)的非線性，采用軟球模型對(duì)塊石顆粒之間碰撞過程進(jìn)行模擬，塊石顆粒之間的接觸力采用非線性彈簧-阻尼模型計(jì)算，該模型在Hertez彈性接觸理論的基礎(chǔ)上考慮粘彈性影響，將塊石之間的接觸力沿法向和切向分解為兩個(gè)分量，其表達(dá)式分別為

(5)

(6)

(7)

(8)

(9)

(10)

式中：Ei和Ej分別為塊石i和塊石j的彈性模量，vi和vj為相應(yīng)塊石的泊松比。

顆粒間的潤(rùn)滑力采用Ball-Melrose模型進(jìn)行計(jì)算，該力與顆粒間的相對(duì)速度成正比，與顆粒間的距離成反比。

2.4 流體與顆粒間作用力控制方程

(11)

拖曳力根據(jù)下式計(jì)算

(12)

升力由下式計(jì)算

(13)

附加質(zhì)量力根據(jù)下式計(jì)算

(14)

Basset歷史應(yīng)力為兩相流中顆粒與流體存在相對(duì)加速度時(shí)所產(chǎn)生的一種非恒定流體力，其表達(dá)式

(15)

2.5 參數(shù)選擇及網(wǎng)格劃分

計(jì)算域尺寸為16 m×10 m×15 m(長(zhǎng)×寬×高)，原點(diǎn)設(shè)置于0號(hào)點(diǎn)處，水流沿X方向流動(dòng)，如圖1所示。采用有限體積法對(duì)計(jì)算域進(jìn)行離散，網(wǎng)格劃分采用結(jié)構(gòu)化六面體網(wǎng)格進(jìn)行，網(wǎng)格尺寸根據(jù)實(shí)際塊石粒徑確定。根據(jù)不同網(wǎng)格尺寸的試算結(jié)果，網(wǎng)格尺寸與塊石顆粒直徑的比值確定為2～3的范圍。流體邊界在OpenFOAM內(nèi)設(shè)置為：入流邊界(面0374)、出流邊界(面5621)、前側(cè)邊界(面0123)、后側(cè)邊界(面7654)，底部為無滑移邊界，頂部為滑移邊界。固體顆粒邊界在LAMMPS內(nèi)設(shè)置為：入流邊界、出流邊界、前側(cè)邊界、后側(cè)邊界，頂部和底部邊界為固定邊界。流體和顆粒的入流邊界與出流邊界互為連續(xù)性周期邊界，前側(cè)邊界和后側(cè)邊界互為連續(xù)性周期邊界。塊石密度2 650 kg/m3，泊松比為0.3，彈性模量為3.755×104MPa，剪切模量為1.444 2×104MPa。為了驗(yàn)證數(shù)值模型的正確性，圖2將計(jì)算得到的單一塊石漂距值與實(shí)驗(yàn)室水槽試驗(yàn)得到的塊石漂距值[14]進(jìn)行了比較，數(shù)值計(jì)算得到的塊石漂距變化規(guī)律與實(shí)驗(yàn)結(jié)果整體上一致，均表現(xiàn)為隨著塊石質(zhì)量(或塊石粒徑)的增大而逐漸減小，計(jì)算值較實(shí)驗(yàn)平均值略微偏大。

圖1 數(shù)值水槽示意圖 圖2 單一塊石漂距計(jì)算值與模型實(shí)驗(yàn)[14]的比較

3 拋石群擴(kuò)散運(yùn)動(dòng)規(guī)律分析

為了考察塊石顆粒群中顆粒級(jí)配的影響，采用兩種粒徑組合方式塊石顆粒組合進(jìn)行模擬計(jì)算。塊石顆粒排列矩陣參數(shù)見表1，表1中塊石1代表粒徑較大的塊石，塊石2代表粒徑較小的塊石。

表1 不等直徑顆粒組合塊石參數(shù)

3.1 拋石群水流方向漂距分析

塊石顆粒群在水流方向上漂距隨時(shí)間變化如圖3所示，塊石顆粒群的平均漂移速度隨塊石平均質(zhì)量的減小而增大，且質(zhì)量越小，塊石群的下落時(shí)間越長(zhǎng)。

圖3 不同粒徑組合方式塊石群水流方向漂距隨時(shí)間變化

3.2 拋石群落點(diǎn)分布分析

圖4為不同時(shí)刻塊石群在水中狀態(tài)，塊石群前端的下落速度比后端慢，個(gè)別塊石的下落速度比塊石群整體速度大。圖5為不同粒徑組合方式塊石組合顆粒群落至底床上的落點(diǎn)分布，圖中um為表面流速，流速垂向分布為指數(shù)為1/12的指數(shù)分布。塊石顆粒群落至床底后的分布形狀較不規(guī)則。雖然質(zhì)量不同，但由于顆粒間的相互影響，并非質(zhì)量小的漂距遠(yuǎn)，質(zhì)量大的漂距近。塊石顆粒群分散之后，并未出現(xiàn)單一質(zhì)量塊石集中在一起的現(xiàn)象，而是呈現(xiàn)大小塊石夾雜在一起的分布規(guī)律。受塊石運(yùn)動(dòng)的影響，塊石群所在區(qū)域的流速減小，這是由于流體質(zhì)點(diǎn)將動(dòng)能傳遞給塊石引起的。

4-a t=0 s4-b t=0.5 s4-c t=1.5

5-a m平均=33.129 kg5-b m平均=97.61 kg5-c m平均=147.502 kg5-d m平均=239.276 kg

4 結(jié)論

基于計(jì)算流體力學(xué)程序OpenFOAM和離散元程序LAMMPS建立水下拋石三維流固耦合數(shù)學(xué)模型，分析拋石群擴(kuò)散運(yùn)動(dòng)規(guī)律，得到以下主要結(jié)論：(1)不同粒徑組合方式塊石顆粒群的平均漂移速度隨塊石平均質(zhì)量的減小而增大，且質(zhì)量越小，塊石群的下落時(shí)間越長(zhǎng)。(2)不同粒徑組合方式塊石組合顆粒落至床底后，并未出現(xiàn)單一質(zhì)量塊石集中在一起的現(xiàn)象，而是呈現(xiàn)大小塊石夾雜在一起的分布規(guī)律。