不同粒徑組合塊石群水下漂移數(shù)值模擬

李 小 超， 張 蓋， 李 衛(wèi) 國， 張 偉， 李 明 益， 常 留 紅

( 1.長沙理工大學(xué) 水利工程學(xué)院， 湖南 長沙 410114；2.水沙科學(xué)與水災(zāi)害防治湖南省重點(diǎn)實(shí)驗室， 湖南 長沙 410114；3.洞庭湖水環(huán)境治理與生態(tài)修復(fù)湖南省重點(diǎn)實(shí)驗室， 湖南 長沙 410114；4.湖南省環(huán)境保護(hù)河湖污染控制工程技術(shù)中心， 湖南 長沙 410114；5.湖南省水運(yùn)建設(shè)投資集團(tuán)有限公司， 湖南 長沙 410011；6.中交天津航道局有限公司， 天津 300461 )

0 引 言

水上拋石施工作業(yè)時塊石在水流作用下會發(fā)生漂移，塊石落點(diǎn)位置不易確定，拋石不到位將造成石料浪費(fèi)及作業(yè)效率降低．塊石在動水中的漂距與拋投方式、塊石自身特性、水流條件等因素有關(guān)，既有一定的規(guī)律性，又呈現(xiàn)出一定的隨機(jī)性，目前很多學(xué)者針對單顆粒塊石漂距提出了不同的經(jīng)驗公式[1-5]，但這些公式大多適用于較小水深和較小流速條件，隨著水深和流速的增加，經(jīng)驗公式得到的塊石漂距與實(shí)測值偏差逐漸增大[5]．此外，實(shí)際拋石施工中多采用抓斗、網(wǎng)兜等方式進(jìn)行機(jī)械拋投，更大型的拋投方式還有底開駁船、側(cè)拋船等船舶直接拋投，即實(shí)際工程中拋石為石料群拋，而非單顆粒拋石．對于群拋石料，下落過程中塊石之間的碰撞增加了確定落點(diǎn)位置的難度．目前對于混合石料群拋漂移特性的研究還非常少，且大部分采用實(shí)驗室水槽進(jìn)行物理模型試驗研究[6-7]．由于影響因素眾多，難以用物理模型試驗對各種影響因素條件下的塊石漂移特性進(jìn)行研究，且物理模型試驗的時間周期和成本均較高，需要建立較為實(shí)用的數(shù)值模型用于群拋石料的漂移計算．

拋石水下漂移運(yùn)動過程可以看作一種顆粒流現(xiàn)象，其中塊石顆粒為固體相，水流為液體相．在顆粒流數(shù)值模擬方面，Sun等[8]將程序OpenFOAM和LAMMPS進(jìn)行耦合，開發(fā)了用于求解固液兩相流的SediFoam開源求解器，流場和顆粒運(yùn)動的計算分別在OpenFOAM和LAMMPS中進(jìn)行，其研究重點(diǎn)關(guān)注水流作用下的泥沙輸移特性．張強(qiáng)強(qiáng)[9]將Fluent和EDEM進(jìn)行耦合用于塊石水中沉降分析，計算信息通過UDF和API接口在兩個軟件之間傳遞．劉卡等[10]采用MFIX-DEM軟件對水下拋石運(yùn)動進(jìn)行了數(shù)值模擬．以上關(guān)于塊石漂移的數(shù)值模擬研究主要針對單顆粒塊石進(jìn)行，且缺乏足夠的試驗數(shù)據(jù)驗證．

航道整治工程中很多拋石工程面臨工期緊、任務(wù)重、級配種類多、驗收標(biāo)準(zhǔn)高等問題，迫切需要對混合石料群拋漂移特性開展系統(tǒng)研究，以指導(dǎo)精準(zhǔn)拋投施工，提升水下拋石質(zhì)量控制和成本控制．本文采用計算流體力學(xué)(CFD)與離散元方法(DEM)耦合的方法對塊石群水下漂移進(jìn)行計算，考慮到OpenFOAM和LAMMPS的開源性，顆粒粒徑組合、亞格子尺度應(yīng)力模型、塊石碰撞模型等更易于實(shí)現(xiàn)，計算在SediFoam求解器中進(jìn)行，重點(diǎn)對大水深、大流速條件下的石料群拋漂移特性及級配的影響進(jìn)行探討．

1 數(shù)值方法

1.1 水流運(yùn)動方程

不可壓縮流體的連續(xù)性方程和動量方程為

?·(εsus+εfuf)=0

(1)

ρfεfg+Ffp)

(2)

式中：εs為固體體積分?jǐn)?shù)；εf=1-εs，為流體體積分?jǐn)?shù)；us為固體顆粒速度矢量；uf為流體速度矢量；p為壓力；τ為流體黏性應(yīng)力張量；ρf為流體密度；g為重力加速度矢量；Ffp為平均化的水流-顆粒相互作用力．

采用大渦模擬對流體方程進(jìn)行求解，考慮塊石漂移具有較高的雷諾數(shù)，亞格子尺度應(yīng)力模型采用單方程亞格子渦黏模型[11]，其表達(dá)式為

(3)

1.2 顆粒運(yùn)動方程

顆粒的平移和轉(zhuǎn)動方程為

(4)

(5)

式中：m為顆粒質(zhì)量；u為顆粒運(yùn)動速度；fcol為相鄰顆粒之間的接觸力，塊石之間的接觸具有較強(qiáng)的非線性，采用考慮黏彈性影響的非線性彈簧-阻尼模型[12]計算；flub是由于流體從顆粒之間的間隙中擠出而產(chǎn)生的潤滑力，采用Ball-Melrose模型[13]計算，該力與顆粒間的相對速度成正比，與顆粒間的距離成反比；ffp是水流-顆粒相互作用力，包括浮力fp、拖曳力fd、升力fl、附加質(zhì)量力fa、Basset歷史應(yīng)力fB，即ffp=fp+fd+fl+fa+fB，拖曳力、升力、附加質(zhì)量力和Basset歷史應(yīng)力的計算方法可分別參見文獻(xiàn)[14]、[15]、[16]、[17]；I為顆粒轉(zhuǎn)動慣性矩；ω為顆粒角速度；Mcol為接觸力扭矩；Mlub為潤滑力扭矩；Mfp為水流-顆粒相互作用力扭矩．

1.3 數(shù)值求解

計算域如圖1所示，為一長方體區(qū)域，尺寸為18 m×10 m×15 m(長×寬×高)，原點(diǎn)設(shè)置于0號點(diǎn)處，水流沿X方向流動．采用有限體積法對計算域進(jìn)行離散，網(wǎng)格劃分采用結(jié)構(gòu)化六面體網(wǎng)格進(jìn)行，根據(jù)不同網(wǎng)格尺寸下的顆粒沉速和漂距試算結(jié)果，網(wǎng)格尺寸與顆粒直徑的比值確定為2～3．河流的尺度非常大，若對整個區(qū)域進(jìn)行計算，計算量將非常大．塊石漂距最遠(yuǎn)只有十幾米，與整個河流相比，只是非常微小的一個河段．對于此問題，本文數(shù)值模擬在X方向上設(shè)置為周期性邊界條件，利用長度小很多(10～20 m)的河段代替實(shí)際河流，以減小模擬尺度，可在不失精確度的前提下，大大降低計算量．流速的垂向分布是影響塊石漂移的主要因素，計算采用工程上常用的指數(shù)分布：

圖1 數(shù)值水槽Fig.1 Numerical tank

u=um(Y/H)n

(6)

式中：um為水流表面流速，Y為距離床面的垂直高度，H為水深，n為指數(shù)．底部為無滑移邊界，頂部為滑移邊界．在Z方向上亦設(shè)置為周期性邊界．

(7)

式中：r為顆粒半徑，ρs為顆粒密度，G為剪切模量，vRay為Rayleigh波波速．上式確定的時間步長小于接觸模型中彈簧-質(zhì)量振子的半個周期．河道治理工程所用石料一般為石灰?guī)r和花崗巖，塊石質(zhì)地較為均勻，根據(jù)文獻(xiàn)[19]取塊石密度為2 650 kg/m3，泊松比為0.3，彈性模量為3.755×104MPa，剪切模量為1.444 2×104MPa，Rayleigh 波波速為2 000 m/s．塊石質(zhì)量范圍為1～300 kg，根據(jù)最大塊石質(zhì)量對應(yīng)的等容粒徑確定最大計算時間步長為2.467×10-6s．

2 結(jié)果與分析

2.1 數(shù)值模型驗證

塊石形狀千變?nèi)f化，很難對復(fù)雜的塊石形狀進(jìn)行模擬，計算時將塊石模擬為球體，采用原型塊石進(jìn)行，粒徑為塊石質(zhì)量對應(yīng)的等容粒徑．水的密度為1 000 kg/m3，運(yùn)動黏度系數(shù)為1×10-6m2/s．本文重點(diǎn)關(guān)注大水深、大流速下的塊石群漂移特性，計算水深為15 m、表面流速為1.5 m/s，指數(shù)n為1/12．圖2將計算流速與水槽模型試驗測試數(shù)據(jù)進(jìn)行了比較．圖3為質(zhì)量25 kg的單顆粒塊石漂移位移隨時間的變化，可以看出，塊石在初始階段有一個短暫的加速過程，之后由于流體與塊石之間的相互作用，塊石漂移速度基本保持穩(wěn)定．漂距為塊石初始下落點(diǎn)與塊石順?biāo)飨鲁恋降状仓潼c(diǎn)之間的X方向水平距離．圖4給出了塊石漂距數(shù)值模擬與模型試驗結(jié)果[5]，兩者吻合較好，數(shù)值模擬得到的漂距變化規(guī)律與模型試驗結(jié)果整體上一致，均表現(xiàn)為隨著塊石質(zhì)量(或塊石粒徑)的增大而逐漸減?。?/p>

圖2 垂向流速分布Fig.2 Vertical flow velocity distribution

圖3 質(zhì)量25 kg的單顆粒塊石漂移位移隨時間的變化Fig.3 Drift displacement against time for single particle rock with mass of 25 kg

圖4 單顆粒塊石漂距Fig.4 Drift distance for single particle rock

文獻(xiàn)[7]中作者在實(shí)驗室水槽內(nèi)開展了混合石料群拋試驗，采用一次拋投的塊石總體積來表征鏟斗、網(wǎng)兜和抓斗等不同機(jī)械拋投方式．本文對抓斗群體拋石漂移進(jìn)行計算，一次拋投塊石總體積為5 m3．初始時刻塊石群為一球體顆粒矩形空間陣列，最上層的顆粒緊貼水面，在X、Y、Z3個方向上顆粒等中心間距排列，間距大小根據(jù)試算確定，且不小于組合中顆粒最大直徑．實(shí)際工程中塊石質(zhì)量范圍比較大，但根據(jù)比尺(質(zhì)量比尺為27 000)換算為模型后，其范圍變得非常?。送?，在從碎石堆中選取塊石時，實(shí)驗人員會習(xí)慣性地參考之前已稱取好的塊石，使得選取到的塊石質(zhì)量大都比較接近，因此計算采用等粒徑顆粒進(jìn)行模擬，塊石個數(shù)與試驗時塊石個數(shù)保持一致，塊石群參數(shù)見表1，表中塊石規(guī)格為實(shí)際工程中使用的拋投石料規(guī)格，排列方式為X、Y、Z3個方向上的塊石排列個數(shù)．圖5為塊石群在底床上的落點(diǎn)分布形態(tài)，圖中m*為塊石平均質(zhì)量．塊石落至底床后沒有出現(xiàn)堆積現(xiàn)象，且分布形狀較為不規(guī)則．下落過程中水流紊動以及塊石之間的碰撞使塊石運(yùn)動方向或速度大小發(fā)生改變，使得集中在一起的塊石逐漸分散開．可以看出，在一次拋投塊石總體積不變的情況下，平均質(zhì)量越大，塊石數(shù)量越少，塊石落至底床后越分散．圖6給出了相應(yīng)的模型試驗塊石群落點(diǎn)分布照片，對比圖5和圖6可見，數(shù)值模擬和模型試驗得到的塊石群落點(diǎn)分布形態(tài)較為接近．

表1 等粒徑塊石群參數(shù)Tab.1 Parameters for rock group with equal particle diameter

(a) m*=33.129 kg

(a) m*=33.129 kg

漂距和塊石群在底床上的落點(diǎn)分布面積是施工過程中比較關(guān)心的兩個參數(shù)．圖7為塊石群平均漂距計算值與試驗值的比較，平均漂距為塊石群中所有塊石漂距的平均值．?dāng)?shù)值模擬與模型試驗結(jié)果較為一致，均隨平均質(zhì)量的增大而減?。畨K石群落點(diǎn)分布的形狀較為不規(guī)則，密集區(qū)域與空白區(qū)域交錯分布，分布面積存在一定的不確定性，本文采用如圖8所示不規(guī)則圖形的面積對塊石群落點(diǎn)分布面積進(jìn)行估算．選取塊石群落點(diǎn)最外層的A、B、C、D4個點(diǎn)，分布面積為(ZC-ZA)×(XB-XD)/2，ZA、ZC、XB、XD為各點(diǎn)的坐標(biāo)值，相當(dāng)于圖中虛線圍成的三角形ABC和三角形ACD的面積之和．計算時主要考慮較為集中的范圍，對于個別離分布集中區(qū)域較遠(yuǎn)的塊石，計算不予考慮，例如圖5(c)中右上角最遠(yuǎn)處的塊石，計算時沒有考慮．塊石群落點(diǎn)分布面積計算結(jié)果見表2，表中模型試驗值為多次拋投的平均值，偏差為相對偏差，即模擬計算值和試驗值之差的絕對值與試驗值的比值．在拋投總體積以及塊石數(shù)量一定的情況下，落點(diǎn)分布面積隨著平均質(zhì)量的增大而減小．表中，除1～200 kg的塊石計算值與試驗值存在一定偏差外，其他計算值均與試驗值吻合較好，偏差的產(chǎn)生與模型試驗塊石形狀不規(guī)則等因素有關(guān)．

圖7 塊石群平均漂距Fig.7 Average drift distance for rock group

圖8 塊石群落點(diǎn)分布面積計算示意圖Fig.8 Calculation diagram for rock group falling point distribution area

表2 塊石群落點(diǎn)分布面積計算值與試驗值的比較Tab.2 Comparison of rock group falling point distribution area between simulations and experiments

2.2 級配的影響

為了考察級配的影響，對4種不同規(guī)格的塊石分別采用4種不同粒徑組合進(jìn)行模擬計算．對于每種塊石規(guī)格，進(jìn)行不同粒徑組合時選取了5種塊石，分別表示為塊石1、塊石2、塊石3、塊石4、塊石5，其質(zhì)量和粒徑見表3，塊石1～5質(zhì)量是以平均質(zhì)量(塊石3的質(zhì)量)為中心，相應(yīng)地增加或減少相同的質(zhì)量得到，以保證平均質(zhì)量不變．對于不同塊石規(guī)格，塊石1、塊石2、塊石3、塊石4、塊石5所代表塊石并不相同．計算時，對于同一規(guī)格塊石，塊石群的總體積、總質(zhì)量、塊石數(shù)量、塊石平均質(zhì)量及排列方式均保持不變．計算采用的粒徑組合方式及各組合中采用的塊石數(shù)量見表4，組合1、組合2、組合3、組合4分別為2、3、5、1種塊石的組合，例如塊石規(guī)格1～50 kg的組合1為100個塊石1和100個塊石5組成的塊石群，組合4為上文討論的等粒徑組合．從組合1到組合4，其不均勻系數(shù)逐漸變小，見表5，即塊石群級配組成逐漸變得更均勻；從連續(xù)性上看，從組合1到組合3，級配組成逐漸變得更連續(xù)，組合4由于缺少其他粒徑的塊石，其連續(xù)性最差．

表3 塊石質(zhì)量和粒徑Tab.3 Mass and particle diameters of rocks

表4 不同粒徑組合中塊石數(shù)量Tab.4 Number of rocks for different particle diameter combinations

表5 不同塊石組合的不均勻系數(shù)Tab.5 Non-uniformity coefficients for different rock combinations

計算均在水深15 m、表面流速1.5 m/s、指數(shù)為1/12的條件下進(jìn)行．圖9～12為4種不同規(guī)格的塊石采用4種不同粒徑組合計算得到的塊石群在底床上的落點(diǎn)分布形態(tài)．從圖中可以看出，不管是何種組合方式，塊石群落點(diǎn)分布形態(tài)都呈現(xiàn)為不規(guī)則圖形，隨著組合方式的變化，分布形態(tài)的變化沒有呈現(xiàn)出任何規(guī)律性．

表6為不同粒徑組合方式下不同塊石規(guī)格的塊石群落至底床上的擴(kuò)散帶寬，水流向擴(kuò)散帶寬和橫向擴(kuò)散帶寬分別為塊石群落點(diǎn)在水流方向和垂直于水流方向上最大坐標(biāo)和最小坐標(biāo)之間的差值．對于相同規(guī)格的塊石群，隨著級配的變化，水流向擴(kuò)散帶寬和橫向擴(kuò)散帶寬均沒有出現(xiàn)單調(diào)的變化規(guī)律．

(a) 組合1

(a) 組合1

(a) 組合1

(a) 組合1

表6 擴(kuò)散帶寬Tab.6 Diffusion breadth

塊石群平均漂距及在底床上的落點(diǎn)分布面積計算結(jié)果分別見圖13、14．與水流向擴(kuò)散帶寬的變化規(guī)律類似，對于相同規(guī)格的塊石群，從組合1變化到組合3，隨著塊石群級配組成的均勻性和連續(xù)性變好，平均漂距和分布面積逐漸增大．但從組合3變化到組合4，由于連續(xù)性變差，平均漂距和分布面積有所減?。畬τ诓煌牧浇M合方式，平均漂距隨平均質(zhì)量的變化趨勢是相似的，均表現(xiàn)為隨著平均質(zhì)量(塊石規(guī)格)的增大而減?。?/p>

圖13 不同粒徑組合的塊石群平均漂距Fig.13 Average drift distance of rock groups with different particle diameter combinations

圖14 不同粒徑組合方式下塊石群落點(diǎn)分布面積Fig.14 Falling point distribution area for rock group with different particle diameter combinations

3 結(jié) 語

本文基于CFD-DEM耦合方法，建立了用于計算拋石水下漂移的數(shù)值模型．模型采用大渦模擬方法對流場進(jìn)行計算，對塊石之間的接觸力、塊石與水流之間相互作用時的各種力進(jìn)行了詳細(xì)模擬，采用物理模型試驗結(jié)果對數(shù)值模型的正確性進(jìn)行了驗證．對不同規(guī)格、不同粒徑組合的塊石群水下漂移進(jìn)行了計算，分析了塊石群在底床上的落點(diǎn)分布形態(tài)和擴(kuò)散范圍、漂距、分布面積等拋石施工過程中比較關(guān)心的塊石漂移特性．結(jié)果表明，級配對塊石群的漂移與擴(kuò)散有較大的影響，塊石群級配組成的均勻性和連續(xù)性越好，平均漂距和分布面積越大．拋石施工過程中抓取得到的塊石粒徑和質(zhì)量分布具有很強(qiáng)的隨機(jī)性，針對具體工程計算時，可參考所有塊石(從石料場運(yùn)來的塊石)的級配構(gòu)造塊石組合．