翟必垚，劉 璐，張寶森，沈洪道，季順迎

（1.大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024；2.黃河水利委員會(huì) 黃河水利科學(xué)研究院，河南 鄭州 450003；3.美國(guó)Clarkson大學(xué) 土木與環(huán)境工程系，波茲坦 紐約 13699-5710）

1 研究背景

河冰在我國(guó)高寒地區(qū)的河流中非常常見，如黃河寧蒙及山東河段、黑龍江上游、松花江依蘭以下河段、嫩江上游等［1］。河冰的動(dòng)力過程非常復(fù)雜，其包含流體力學(xué)、固體力學(xué)，以及流體和固體的相互作用，其中流凌和冰壩是比較常見并備受關(guān)注的河冰現(xiàn)象［2-4］。流凌不僅會(huì)影響航運(yùn)，其巨大的動(dòng)能還會(huì)撞擊破壞橋墩和其他水工構(gòu)筑物。冰壩的形成會(huì)減小河道的正常過流斷面，加之封凍冰層對(duì)水的摩擦阻力，必然阻隔上游來水導(dǎo)致槽蓄水量增加，不斷抬高上游水位，可在極短時(shí)間內(nèi)形成凌汛險(xiǎn)情，造成嚴(yán)重的自然災(zāi)害［5-6］。

現(xiàn)場(chǎng)觀測(cè)、物模試驗(yàn)和數(shù)值計(jì)算一直是河冰研究的三個(gè)主要手段，而隨著研究認(rèn)識(shí)的深入和計(jì)算機(jī)的發(fā)展，數(shù)值模擬越來越得到人們的重視。國(guó)外對(duì)冰塞冰壩的研究最早開始于Pariset 和Hausse（1961）對(duì)冰蓋前緣水動(dòng)力條件和冰蓋穩(wěn)定時(shí)平衡狀態(tài)的受力分析，并由此建立了判斷冰蓋是否向上游發(fā)展的弗勞德數(shù)條件以及穩(wěn)定冰塞的平衡方程［7］。后來其他研究者在其基礎(chǔ)上建立一維靜態(tài)冰塞冰壩模型［8-10］，包括RIVJAM、HEC-RAS等［11］。Beltaos和Burrell對(duì)冰壩的形成、演變和潰決過程中的水位和冰厚進(jìn)行了大量的觀測(cè)和分析，并與模型結(jié)合得到了很好的模擬結(jié)果［12-13］。這些模型能夠計(jì)算冰塞冰壩的形狀和水位狀態(tài)，其不足之處是無法計(jì)算冰塞形成的動(dòng)力過程、冰水相互影響的過程，因此也無法預(yù)測(cè)冰塞發(fā)生的概率、位置和時(shí)間［14-15］。之后，一維和二維河冰動(dòng)力學(xué)模型相繼發(fā)展起來。并可以描述河冰的輸移和初始冰塞的形成［16］，如DynaRICE［14，17］、SPIKI［18］等。同時(shí)，國(guó)內(nèi)對(duì)河冰的研究也取得了很大進(jìn)展，從河冰力學(xué)性質(zhì)研究［19-20］到數(shù)學(xué)模型研究［21-25］都做了大量的工作，并結(jié)合我國(guó)河流的具體情況，在工程應(yīng)用中取得了理想的研究成果［26-28］。

然而，對(duì)于以上提到的河冰數(shù)值模型，大部分都是將河冰作為連續(xù)體而忽略了冰塊擠壓破碎的動(dòng)力過程。離散元法則可以很好地表征冰塊運(yùn)動(dòng)時(shí)的接觸過程和黏結(jié)破碎過程，近年來在海冰領(lǐng)域被廣泛采用［29-31］。在海冰離散元研究中，因?yàn)樗畡?dòng)力不如河流中那樣強(qiáng)烈，因此大多沒有考慮海冰對(duì)水力條件的影響。在河冰領(lǐng)域中，離散元法的應(yīng)用還較為少見。Hopkins 等［32-34］開展了前期初步研究，并發(fā)現(xiàn)離散元法在冰塞冰壩形成、冰塊與橋墩、攔冰柵等結(jié)構(gòu)相互作用等方面具有很強(qiáng)的優(yōu)勢(shì)。離散元法中最早采用球形顆粒和圓盤顆粒，近年來不同類型的非規(guī)則單元已相繼發(fā)展起來［35-37］。擴(kuò)展多面體單元作為一種新的離散單元構(gòu)造方法，在一定程度上避免球形離散單元中細(xì)觀計(jì)算參數(shù)的經(jīng)驗(yàn)性［38］。針對(duì)冰蓋武開河產(chǎn)生的碎冰以及由冰與河岸、冰與橋墩等水工結(jié)構(gòu)物的相互作用而破碎形成的碎冰，擴(kuò)展多面體單元可更加真實(shí)地反映碎冰的幾何形態(tài)。

為此，本文將采用擴(kuò)展多面體單元構(gòu)造河冰離散元方法，采用有限元方法求解河冰影響下的水動(dòng)力學(xué)，將河冰離散元方法與水動(dòng)力學(xué)有限元方法相耦合以建立河冰的動(dòng)力學(xué)模型。采用以上方法對(duì)明流冰封過程的水力條件變化過程進(jìn)行數(shù)值模擬，驗(yàn)證本文計(jì)算模型的有效性，并對(duì)冰壩的形成過程進(jìn)行數(shù)值模擬。

2 離散元方法與水動(dòng)力學(xué)耦合的河冰動(dòng)力學(xué)模型

本文建立的河冰動(dòng)力學(xué)模型主要包括四個(gè)部分，即：河冰的動(dòng)力學(xué)方程求解、擴(kuò)展多面體單元求解冰塊受到的接觸力、二維水動(dòng)力模型求解水力條件，以及耦合模型中冰水信息的傳遞和冰水相互影響的過程。

2.1 河冰動(dòng)力學(xué)方程河冰輸移主要受到風(fēng)和水流的拖曳力、水面梯度力、冰塊間的相互作用、河冰與河岸河床及水工建筑物的相互作用，如圖1所示。由此，河冰運(yùn)動(dòng)的動(dòng)力方程式為［14］：

式中：Mi為河冰單元的質(zhì)量；vi為河冰單元的速度；t為時(shí)間；G為重力沿水面坡度的分量；Fa為風(fēng)對(duì)冰塊的拖曳力；Fw為水對(duì)冰塊的作用力，包括水的拖曳力、浮力；Ri為河冰接觸力，包括河冰冰塊間、河冰與河岸河床及水工建筑物的接觸力，其主要通過離散元方法求得。

圖1 河冰動(dòng)力過程的受力分析

2.2 河冰擴(kuò)展多面體離散元方法擴(kuò)展多面體離散元方法是基于多面體離散元和球體離散元發(fā)展起來的，通過Minkowski Sum理論進(jìn)行構(gòu)造［39］，可寫作：

圖2 基于Minkowski Sum理論擴(kuò)展多面體的構(gòu)造

式中：A、B為兩個(gè)單元的空間體；x、y分別為基礎(chǔ)多面體A和擴(kuò)展球體B內(nèi)的坐標(biāo)矢量［39-40］，即按照Minkowski Sum理論求和后，會(huì)得到具有光滑表面的擴(kuò)展多面體，如圖2所示。接觸搜索判斷也由角點(diǎn)與棱邊的搜索判斷轉(zhuǎn)化為球面與柱面的搜索判斷，由此有效提高計(jì)算效率。

在離散元求解過程中需要確定每個(gè)時(shí)間步接觸單元間的接觸點(diǎn)以及單元重疊量。單元的接觸判斷是基于包絡(luò)函數(shù)的優(yōu)化求解方法，引入二階多面體擴(kuò)展函數(shù)g，其定義的幾何形態(tài)與擴(kuò)展多面體具有很高的相似性，可寫作：

式中：g（x、y、z）為多面體擴(kuò)展函數(shù)；N為多面體的表面數(shù)量；ai、bi和ci分別表示多面體第i個(gè)面單位外法向的三個(gè)分量；di是第i個(gè)面到坐標(biāo)原點(diǎn)的距離；r是函數(shù)的擴(kuò)展半徑；是Macaulay括號(hào)，且有即求點(diǎn)(x，y，z)到某個(gè)面的距離。

通過二階多面體擴(kuò)展函數(shù)與球面函數(shù)加權(quán)求和的函數(shù)形式可構(gòu)造具有多面體特征的光滑顆粒，即：

式中：f（x、y、z）為擴(kuò)展多面體外表面函數(shù)；k為顆粒光滑度系數(shù)，0 ＜k≤1。當(dāng)k越小時(shí)，單元越接近多面體形態(tài)；R為球面函數(shù)的半徑。由該函數(shù)定義擴(kuò)展多面體的包絡(luò)函數(shù)并建立相應(yīng)的優(yōu)化方程：min(fA+fB)，s.t.fA-fB=0，fA和fB分別是兩個(gè)顆粒的形態(tài)方程。通過光滑度系數(shù)k由大到小漸變的方法計(jì)算迭代初始點(diǎn)，并在基本多面體上進(jìn)行迭代確定最終接觸點(diǎn)和接觸的重疊量δ=min‖XA-XB‖-(rA+rB)，XA和XB分別為擴(kuò)展多面體內(nèi)部基本多面體表面上的點(diǎn)，rA和rB是擴(kuò)展半徑［41］。

擴(kuò)展多面體單元間的作用力采用非線性接觸模型，可寫作：

式中：Fn為法向力，由法向彈性力和法向阻尼力組成；Ft為切向力，由切向彈性力F te和切向阻尼力組成；kn為法向接觸剛度，E *和R*分別是等效彈性模量和等效顆粒半徑，νA、νB為顆粒的泊松比；EA、EB為顆粒的彈性模量；RA、RB為顆粒的半徑；δn、δt分別是法向和切向重疊量；分別是法向和切向相對(duì)速度；cn、ct分別是法向和切向阻尼系數(shù)，其中mAB為等效質(zhì)量mA、mB為顆粒質(zhì)量；ζn為無量綱阻尼系數(shù)，可表示為e是顆?；貜椣禂?shù)；δtmax為最大切向重疊量，表征最大靜摩擦來限制切向接觸力，即δtmax=μδn( 2-ν)( 2-2ν)，這里μ為顆粒間摩擦函數(shù)，ν為泊松比。

2.3 二維河流水動(dòng)力模型河冰動(dòng)力模型中的水動(dòng)力學(xué)采用DynaRICE 模型中的水力方程與求解方法［14］。河流水動(dòng)力學(xué)方程是基于動(dòng)量守恒的二維淺水方程推導(dǎo)而來并考慮了表面冰的影響［14］。整個(gè)控制方程包括連續(xù)性方程和動(dòng)量方程，其分別寫作：

式中：x、y和t分別是空間和時(shí)間變量；H、Ht分別是總的水深和冰蓋下的水深；η是水面高程；qtx、qty是總的流量沿x方向和y方向的分量；Nice是冰的密集度；t′是冰在水中的厚度；ρw為水的密度；τs、τb分別是冰對(duì)水流及河床對(duì)水流的摩擦力；τsx、τsy、τbx、τby分別是其沿x和y方向的分量；這里εxy是廣義渦流黏性系數(shù)。

河水流量分為冰層流量和冰下流量?jī)刹糠?，即qt=qu+ql。冰層中的水流一方面會(huì)隨著冰一起運(yùn)動(dòng)，一方面也會(huì)因?yàn)楹恿鞅砻娴乃狡露榷l(fā)生滲流。

單位河寬上冰層部分的流量可以寫成：

式中：qi為隨冰塊一起運(yùn)動(dòng)的河水流量；vi為冰量；qs=VsΔH為在冰層中運(yùn)動(dòng)的滲流流量；Vs=usi+vs j為滲濾流速，us和vs是其沿x和y方向的分量；冰層中的滲流是由河流的水力坡降產(chǎn)生；滲流流速Vs與水面坡降J的關(guān)系可以表示為［42］：

式中：λ為滲流系數(shù)，與冰的孔隙率、形狀、冰塊的大小等因素有關(guān)，一般取λ=1.0～1.6m/s。

2.4 擴(kuò)展多面體離散元方法與水動(dòng)力學(xué)的耦合模型河冰動(dòng)力過程中的強(qiáng)水動(dòng)力作用使得冰與水流間的相互作用不可忽略。本文建立的河冰動(dòng)力學(xué)模型對(duì)水和冰的耦合處理通過以下幾個(gè)方面進(jìn)行：在水力模型方程中考慮表面冰對(duì)水流的影響，基于擴(kuò)展多面體離散單元與三角形有限單元的相對(duì)位置通過插值平均進(jìn)行信息交換，水對(duì)冰的作用通過對(duì)離散單元上的受力計(jì)算完成。

（1）水與冰間的計(jì)算參數(shù)交換。河冰的離散元計(jì)算是在拉格朗日坐標(biāo)下進(jìn)行的，水動(dòng)力的有限元法計(jì)算則是基于歐拉坐標(biāo)系。河冰的厚度、速度等信息與水動(dòng)力學(xué)信息之間的傳遞是通過對(duì)離散單元與有限元網(wǎng)格節(jié)點(diǎn)相對(duì)位置進(jìn)行插值計(jì)算得到。假設(shè)P為某一離散單元的質(zhì)心，其所在的三角形網(wǎng)格的三節(jié)點(diǎn)分別為A、B、C，河冰計(jì)算所需要的流速、水深等水力信息由所在網(wǎng)格節(jié)點(diǎn)的信息加權(quán)求得，如圖3所示，fP=∑Φi fi，其中Φi=Si Stotal（i=a，b，c），Si是第i個(gè)三角形的面積，Stotal是整個(gè)三角形單元的總面積，fP可以為流速、水面高度、水深等水力變量。

圖3 冰與水間計(jì)算參數(shù)傳遞計(jì)算

然而，河流水力計(jì)算所需要的河冰信息則由該有限單元網(wǎng)格內(nèi)的所有離散元顆粒信息的平均值求得，均值結(jié)果作為該有限單元網(wǎng)格內(nèi)的值，即其中，下標(biāo)E是該三角形單元編號(hào)，i表示三角單元中的離散單元編號(hào)，NP是三角單元中的離散顆粒數(shù)總和，f表示河冰離散元顆粒的速度、剖面厚度等參數(shù)。網(wǎng)格節(jié)點(diǎn)的值通過其所在所有單元網(wǎng)格的值進(jìn)行平均求得，如圖3中節(jié)點(diǎn)A的參數(shù)

（2）浮力、拖曳力及水面梯度力。水對(duì)冰的作用力主要有浮力、拖曳力以及水面梯度力，如圖4所示。浮力大小與冰塊浸沒在水下的體積有關(guān)，通過多面體質(zhì)心與多面體水下部分每個(gè)三角面組成的四面體體積Vi求和計(jì)算，即Vsub=∑Vi。冰塊的總浮力Fb=Vsub ρw gk，k為豎直單位矢量；冰塊受到的浮力矩為：Mb=rb×Fb，其中rb為由質(zhì)心指向水下浮心的向量。浮心的確定可由公式計(jì)算，這里Oi是水下四面體的質(zhì)心。

圖4 河冰離散單元受浮力、拖曳力與水面梯度力計(jì)算

水流對(duì)冰塊的拖曳力Fd和拖曳力矩Md可分別寫作［31］：

式中：、C dM是拖曳力和拖曳力矩系數(shù)；vw、vi分別為流速和冰速；Asub為冰塊與水的接觸面積在水流方向的投影；reff是冰塊的有效半徑，即塊體所有頂點(diǎn)到質(zhì)心距離的平均值；ωi是塊體轉(zhuǎn)速。

風(fēng)對(duì)冰塊的拖曳力與流對(duì)冰的拖曳力原理與計(jì)算方式相同，這里不再說明。通過這種方法對(duì)水流拖曳力和拖曳力矩進(jìn)行計(jì)算，考慮了冰塊與水的接觸面在水流方向的投影面積對(duì)水流拖曳力的影響，拖曳力造成冰塊的平移，拖曳力矩造成冰塊的翻轉(zhuǎn)。然而，由于水動(dòng)力計(jì)算采用了二維的水力模型，忽略了流速沿水深方向的變化，因此公式中的流速實(shí)際上是此處沿水深方向的平均流速。并且計(jì)算中未考慮在豎直方向的流速對(duì)冰塊的擾動(dòng)，冰塊豎直方向的運(yùn)動(dòng)只受到浮力和浮力矩的影響。在河冰輸移過程中，豎直方向的流速相對(duì)于水流流動(dòng)方向的速度較小，對(duì)河冰的影響也比較小，一定程度上可以忽略。但是在形成冰塞或冰蓋后，冰蓋前緣處的水流豎向擾動(dòng)會(huì)促使前緣冰塊的翻轉(zhuǎn)下潛或者豎直下潛，對(duì)冰塞的形成有一定的影響，這會(huì)在之后的模型發(fā)展中考慮。

河冰的水面梯度力是由水面坡降造成的冰塊重力沿水面下降坡度的分量，即G=Gx i+Gy j，其中，ηice是該河冰顆粒對(duì)應(yīng)位置的水面高度，和是沿x方向和y方向的水面梯度。

（3）冰、水計(jì)算時(shí)間步耦合。數(shù)值計(jì)算中的時(shí)間步是影響計(jì)算穩(wěn)定性、計(jì)算精度和計(jì)算效率的重要參數(shù)。離散元方法假定一個(gè)計(jì)算時(shí)間步內(nèi)，顆粒受力和加速度保持恒定不變，這里采用了瑞利波法時(shí)間步進(jìn)行計(jì)算［43］

式中：rmin為計(jì)算顆粒中最小擴(kuò)展半徑；ρi為河冰離散元顆粒的密度；剪切模量E為彈性模量。在實(shí)際計(jì)算中，要依據(jù)顆粒運(yùn)動(dòng)的劇烈程度選取合適的時(shí)間步長(zhǎng)以保證計(jì)算的穩(wěn)定性，這里取進(jìn)行計(jì)算。

在水動(dòng)力學(xué)有限元計(jì)算中，臨界時(shí)間步長(zhǎng)為：

式中：ΔL為計(jì)算域中三角形網(wǎng)格的最小邊長(zhǎng)；Hmax為最大水深。在實(shí)際計(jì)算中，根據(jù)流場(chǎng)水力變化的劇烈情況會(huì)稍微調(diào)整步長(zhǎng)，一般取

水與冰間的計(jì)算參數(shù)交換頻率由耦合時(shí)間步?jīng)Q定，耦合時(shí)間步根據(jù)實(shí)際冰水作用的劇烈程度進(jìn)行人為設(shè)置，一般取10～60 s。對(duì)于冰水動(dòng)力過程較為劇烈的情況采用小的耦合時(shí)間步，對(duì)于冰水動(dòng)力過程相對(duì)緩和的情況采用大耦合時(shí)間步。每當(dāng)模擬時(shí)間為耦合時(shí)間步的整數(shù)倍時(shí)，即可交換水冰間的計(jì)算參數(shù)。具體過程如圖5所示。

圖5 計(jì)算耦合時(shí)間步

3 矩形河道內(nèi)封河水位變化的數(shù)值模擬

將以上河冰動(dòng)力學(xué)模型應(yīng)用到規(guī)則河道進(jìn)行簡(jiǎn)單工況的數(shù)值模擬，并通過與DynaRICE的計(jì)算結(jié)果和理論分析結(jié)果對(duì)比來檢驗(yàn)?zāi)Ｐ偷挠行耘c準(zhǔn)確性。DynaRICE 模型由美國(guó)Clarkson大學(xué)沈洪道教授建立，已在多條河流動(dòng)力學(xué)模擬中得到驗(yàn)證［14，17］并在國(guó)際上被廣泛認(rèn)可。

開闊水域在冬季經(jīng)常遇到封河的現(xiàn)象，河流上方的冰蓋會(huì)改變河道的綜合糙率及水力半徑，從而使河流的水力條件發(fā)生很大變化。先分析開闊水域產(chǎn)生冰蓋后水力條件的變化，然后再對(duì)該過程進(jìn)行數(shù)值模擬。

圖6 明渠與暗渠水力參數(shù)示意

曼寧公式可反映水流與河床的部分關(guān)系以及河床內(nèi)部各因素之間的相互關(guān)系。這里假設(shè)封河情況下，冰蓋不向下游運(yùn)動(dòng)，則通過曼寧公式可求得規(guī)則河道內(nèi)的均勻流水深，其基本表達(dá)式為：

式中：Q為水文斷面的流量，m3/s；Ac為過水?dāng)嗝娴拿娣e，m2；n為糙率；Rh為水力半徑，m；Sf為水力坡降。冰蓋的出現(xiàn)改變了明渠流的水力條件，如圖6所示，其在曼寧公式的使用上有一定的區(qū)別。

在規(guī)則矩形河道明渠均勻流中，糙率等于河底糙率，即n=nb；水力半徑Rh=，這里Ac=Bcdw；濕周Pc=Pb=2dw+Bc，Bc和dw分別是河寬和水深。式（16）整理后為h為水深，這里等于dw。

在河道有冰蓋的條件下，濕周增加了上層冰蓋部分，即Pc=Pb+Pi=( 2dw+Bc)+Bc=2(dw+Bc)，糙率等于河底糙率nb與冰蓋底部糙率ni的綜合糙率。綜合糙率的計(jì)算公式形式多樣，被廣泛應(yīng)用的有Belokon-Sabaneev 公式、Einstein 公式和Larsen 公式［44］。由于在模擬中采用的是寬深比較小的矩形明渠，因此這里采用誤差較小的Einstein公式計(jì)算綜合糙率，即：

式中α為床面區(qū)和冰蓋區(qū)的濕周比，此處由此整理后得到最終水深需考慮冰蓋中的部分，即這里hi為冰厚。通過曼寧公式計(jì)算的正常水深，將作為模擬結(jié)果準(zhǔn)確性檢驗(yàn)的參照。當(dāng)表面冰存在向下游流動(dòng)的速度時(shí)，以上分析則需要重新考慮冰速的影響［17］。

這里規(guī)則矩形河道的長(zhǎng)度L=2000 m，寬度Bc=10 m，如圖7（a）所示。河底坡降Sf=0.4×10-3，河床糙率nb=0.035，冰蓋糙率ni=0.02，上游邊界設(shè)定恒定流量Q=3.17m3/s，下游邊界設(shè)定正常水深h=0.59 m，主要計(jì)算參數(shù)列于表1。圖8為均勻恒定明渠流初始時(shí)刻的縱向水力參數(shù)。在明渠流的穩(wěn)定狀態(tài)下，沿河道方向各處流速基本穩(wěn)定在0.55 m/s，水深穩(wěn)定在0.59 m。

在穩(wěn)定明渠流的水面上生成0.2 m厚的冰蓋，冰塊單元分布如圖7（b）所示，上下游的邊界條件恒定不變，計(jì)算初始的水位流速保持與明渠流一致。隨著冰蓋的出現(xiàn)，附加的冰邊界增加了水流流動(dòng)的阻力并導(dǎo)致過流能力減小。采用本文模型計(jì)算了2個(gè)小時(shí)的河冰動(dòng)力學(xué)過程，結(jié)果如圖9所示。這里給出2 min、10 min、1 h和2 h時(shí)的冰水狀態(tài)。

表1 規(guī)則矩形河道明渠與封河的模擬參數(shù)

圖7 明渠流與有冰暗渠流模擬

圖8 明渠流的初始水力參數(shù)

在保持上游流量以及下游水位邊界條件恒定的情況下，當(dāng)河道存在冰蓋時(shí)，水位先由上游開始壅高并逐漸往下游抬高，最終在下游水位邊界條件的影響范圍以外基本都達(dá)到了正常水深。采用Dy?naRICE對(duì)這一過程進(jìn)行了數(shù)值模擬，計(jì)算結(jié)果與本文方法的對(duì)比情況如圖10所示。圖10（a）和圖10（b）分別是流速和水深的對(duì)比。在初始階段2 min時(shí)，流速和水深較DynaRICE 的結(jié)果都有一些偏差；在10 min時(shí)，偏差顯著縮?。辉谥蟮?～2 h時(shí)，兩個(gè)模型的結(jié)果幾乎完全吻合。這在一定程度上驗(yàn)證了本文計(jì)算模型的有效性。

以上計(jì)算中的前期偏差主要是離散元模型中水面的壅高過程稍稍滯后于DynaRICE。這是由于離散元的冰塊顆粒在水中是三維的運(yùn)動(dòng)形式，在初期水位快速抬升的過程中冰塊會(huì)有不穩(wěn)定的上下浮動(dòng)，冰下厚度隨之不斷變化，而冰對(duì)水流的糙率與冰的水下厚度有關(guān)。這就使得水流受到的冰下摩擦力發(fā)生了變化，因此計(jì)算得到的流速和水深隨之改變。DynaRICE是二維模型，所以不存在豎向的運(yùn)動(dòng)。這使得兩個(gè)模型在初始階段出現(xiàn)結(jié)果上的偏差。當(dāng)水位逐漸穩(wěn)定后，這種情況便會(huì)消失，結(jié)果自然與之相吻合。

根據(jù)前文對(duì)明渠和有冰蓋暗渠兩種情況下曼寧公式的理論分析，可以求得在明渠的狀態(tài)下正常水深hw=0.59 m；在有冰蓋的情況下求得的冰下水深dw=0.73 m，加上冰蓋部分的影響，最終的正常水深hw=0.91 m。由本文模型計(jì)算得到的最終正常水深為0.86 m，與理論分析對(duì)比得誤差為5.49%。對(duì)于流速分布，由于冰蓋底部及河床處的摩擦力的影響，冰下水流速度沿水深方向并非均勻分布，并且在冰蓋層內(nèi)也有水流通過，而本文模型中的二維水力模型部分是將豎向的流速平均化處理，因此這里只能將模擬結(jié)果中沿豎向平均流速（包括冰蓋內(nèi)部的流速和冰蓋下的流速）與理論分析結(jié)果中的平均流速做對(duì)比。模型計(jì)算結(jié)果中正常水深位置的平均流速為0.37 m/s，理論分析結(jié)果中的平均流速為0.35 m/s，對(duì)比得到誤差為5.71%。由于本文模型考慮了冰蓋中的水流流量，而在理論分析中較難考慮到總流量在冰蓋中及冰蓋下的分配比例而忽略了冰蓋中的流量，因此計(jì)算模型和理論分析存在較小誤差是在可接受的范圍之內(nèi)。

圖9 封河水位抬升過程模擬結(jié)果

圖10 本文模擬的流速和水深變化同DynaRICE結(jié)果的對(duì)比

4 河冰輸移、堆積及冰壩形成過程的數(shù)值模擬

河冰的動(dòng)力過程往往伴隨著水力條件的快速變化，特別是在冰塊堆積、堵塞河道并形成冰壩時(shí)，過水?dāng)嗝婕眲p小，上游水位迅速升高。通過對(duì)河冰輸移、堆積、阻塞形成冰壩過程的數(shù)值模擬，檢驗(yàn)本文河冰動(dòng)力學(xué)模型的有效性。

將河道調(diào)整為長(zhǎng)度L=2000 m，寬度Bc=6 m。在河道x=700 m 的位置處設(shè)置攔冰柵，初始時(shí)刻在攔冰柵前方距上游邊界600～700 m的位置生成了100 m長(zhǎng)5.8 m寬的連續(xù)流冰區(qū)域。冰塊的平均尺寸為0.1 m2，面密集度為0.8，冰厚為0.1 m，冰塊之間的間隔在0.08 m左右，冰塊的總數(shù)量為4060個(gè)。為保證冰水作用引起的水力變化不會(huì)影響到河道邊界，在冰區(qū)距離上下游邊界都留有一定長(zhǎng)度的開闊水域。計(jì)算參數(shù)采用表1中的數(shù)值，上游邊界設(shè)定恒定流量Q=9.45m3/s，下游邊界設(shè)定恒定水深h=1.60m，流速為1.12m/s左右。初始的冰塊分布及相關(guān)變量如圖11所示。

河冰在外力作用下向下游運(yùn)動(dòng)，遇到冰柵后開始聚集堆積。數(shù)值模擬的河冰動(dòng)力過程為20 分鐘，對(duì)整個(gè)過程的模擬結(jié)果如圖12所示。河冰的初始密集度為0.8，相互間有一定間隙，如圖12（a）所示。當(dāng)冰塊遇到攔冰柵后，開始在攔冰柵前方聚集。隨著上游來冰的增加，下游冰塊承受的壓力越來越大并開始相互擠壓。靠近攔冰柵的冰塊由平鋪水面被壓至傾斜或直立。冰塊在水流動(dòng)力作用以及冰塊相互擠壓作用下發(fā)生進(jìn)一步擠壓，開始形成水下塞體，河道斷面過水能力下降，上下游水位開始變化，如圖12（b）。隨著來冰量的增加，對(duì)冰壩頭部的作用力進(jìn)一步加大，更多的浮冰塊被擠入水下，冰壩的厚度快速增加，如圖12（c）。由于過水?dāng)嗝孢M(jìn)一步縮小，河道蓄水量進(jìn)一步增加，上游水位也明顯抬高，下游水位明顯降低，如圖12（d）（e）所示，此時(shí)冰壩堆積厚度已基本達(dá)到穩(wěn)定狀態(tài)，冰壩開始向上游延伸發(fā)展。

圖11 初始河冰分布以及相關(guān)變量

圖12 河冰堆積過程的離散元模擬結(jié)果

離散元法在計(jì)算中考慮了各個(gè)冰塊所受的作用力及運(yùn)動(dòng)狀態(tài)的變化。在冰壩形成過程的模擬中，模擬結(jié)果重現(xiàn)了冰壩形成穩(wěn)定后河道表面冰的三種不同的分布狀態(tài)，如圖13所示?？拷鼣r冰柵的區(qū)域，稱為立封區(qū)，冰塊呈豎直交錯(cuò)、雜亂重疊的堆積狀態(tài)。該區(qū)域承受著上游冰區(qū)的擠壓，是整個(gè)河道內(nèi)力最高的區(qū)域，冰塊間相互作用最嚴(yán)重。該位置水流速度最大，水流對(duì)冰的拖曳作用最強(qiáng)，同時(shí)此處也是上游高水位區(qū)向下游低水位區(qū)過渡的區(qū)域。水面坡度大，水面梯度力也最大，從而加劇了冰塊之間的擠壓。在冰區(qū)的中部，冰塊平鋪于河道，排列緊密，相互接觸但不重疊或者略微重疊的狀態(tài)，稱為平封區(qū)。此處屬于冰內(nèi)力較小的區(qū)域，外力使冰塊緊密排列但又不足以使冰塊重疊翻轉(zhuǎn)和堆積。在冰壩尾端，冰塊密集度比較低，呈自由漂移狀態(tài)，排列較為稀疏隨機(jī)，稱為冰塊的稀疏區(qū)。此處的水面由上游的正常水位過渡到冰區(qū)的水位稍稍有升高的趨勢(shì)，即水面梯度指向河道上游，重力分量指向上游在一定程度上平衡了水流的拖曳作用，使得冰自由漂移在冰壩尾部。圖14是穩(wěn)定冰壩的冰壩剖面與水力要素的沿程分布情況。由圖可見在冰壩體最低處過水?dāng)嗝婷娣e急劇減小，流速急劇增大，達(dá)到3 m/s，弗勞德數(shù)也達(dá)到了最大值1.5。在冰蓋前緣處的弗勞德數(shù)和流速分布為0.15和0.7 m/s。在冰壩上游端至中部x=0.68 km處，弗勞德數(shù)處于0.17左右，流速穩(wěn)定在0.75 m/s，河面冰區(qū)分布狀態(tài)為冰壩頭部的自由漂浮和冰壩中部的平封狀態(tài)。在x=0.68 km 至下游冰壩頭部冰厚逐漸增大，流速和弗勞德數(shù)也快速增加至最大值，這時(shí)的表面冰狀態(tài)也由平封狀態(tài)過渡至雜亂交錯(cuò)的立封狀態(tài)。由于本文尚未考慮冰蓋前緣的復(fù)雜流場(chǎng)以及使冰塊下潛翻轉(zhuǎn)的作用力，因此這里平封立封過渡的臨界弗勞德數(shù)和臨界流速與前人研究結(jié)果相比數(shù)值偏大［45-46］，只能反映冰壩形成后沿程水力要素的變化趨勢(shì)。

圖13 河面上河冰的三種分布狀態(tài)

圖14 冰壩剖面與水力要素的沿程分布情況

模型計(jì)算了冰壩形成過程中攔冰柵受到的冰壓力隨時(shí)間的變化過程，如圖15所示。在冰壩形成初始階段，冰柵由于阻攔了大量來冰，其所受冰力快速增大。隨著冰量的增加，冰柵前的冰塊變得密集，并擠壓河道兩側(cè)，河道兩側(cè)邊界開始承擔(dān)部分冰力。隨著冰壩長(zhǎng)度的增加，河道兩側(cè)承載冰力隨之增大。上游來冰造成的冰力由兩側(cè)邊界承擔(dān)的部分增加，傳遞到攔冰柵的部分減少，因而攔冰柵所受冰力的增加速度漸漸變緩，最終冰柵所受冰力隨冰壩形成穩(wěn)定后趨于穩(wěn)定。

圖15 冰壩形成過程中攔冰柵承受冰壓力隨時(shí)間的變化

圖16 冰壩形成過程中水位隨時(shí)間的變化

對(duì)整個(gè)冰壩形成過程中水位的變化情況也進(jìn)行了數(shù)值模擬，計(jì)算結(jié)果如圖16所示。圖中列出攔冰柵上游0.2 km、0.05 km處以及下游距離攔冰柵0.05 km、0.2 km處四個(gè)水位測(cè)點(diǎn)（具體位置如圖11所示）的水位變化。從中可以看出，在冰壩的影響下，上游的水位抬升近0.24 m，而下游的水位降低近0.11 m，并且攔冰柵上游兩個(gè)測(cè)點(diǎn)的初始水位差近0.15 m，水位抬升后水位差變?yōu)榻?.08 m。由此可以看出，在冰壩的影響下，上游的水面坡降減小，水面變得緩和，而下游的水面坡降未發(fā)生很大變化。

5 結(jié)論

本文建立了河冰離散元方法與二維水動(dòng)力學(xué)耦合的河冰動(dòng)力學(xué)計(jì)算模型，河冰的運(yùn)動(dòng)過程采用離散元方法求解，河流的水動(dòng)力學(xué)采用有限元方法對(duì)河冰影響下的二維淺水方程進(jìn)行求解。同時(shí)采用了擴(kuò)展多面體離散單元對(duì)開河時(shí)期冰塊幾何形狀隨機(jī)分布的特點(diǎn)進(jìn)行表征。通過對(duì)明渠出現(xiàn)冰蓋之后的水位變化過程進(jìn)行數(shù)值模擬，計(jì)算結(jié)果很好地吻合了DynaRICE 的模擬結(jié)果和理論分析結(jié)果，從而驗(yàn)證了河冰動(dòng)力學(xué)模型的有效性。在此基礎(chǔ)上將河冰動(dòng)力學(xué)模型應(yīng)用于河冰輸移、堆積和冰壩形成過程的模擬中。計(jì)算結(jié)果準(zhǔn)確地描述了冰壩形成過程中河冰的運(yùn)動(dòng)狀態(tài)以及上游水位雍抬、下游水位降低、冰下流速分布變化過程以及冰柵所受冰壓力變化過程。同時(shí)再現(xiàn)了冰壩形成穩(wěn)定后河面上冰塊的分布狀態(tài)，由冰壩頭部的豎直交錯(cuò)、雜亂重疊，冰壩中部的平坦密集，到冰壩尾部的稀疏漂浮，與現(xiàn)場(chǎng)觀察現(xiàn)象相吻合。結(jié)果表明，該模型可以準(zhǔn)確地從細(xì)觀角度表征河冰在輸移聚集堆積過程中復(fù)雜的運(yùn)動(dòng)情況，并可準(zhǔn)確地獲得河渠內(nèi)水力條件的變化過程。本文所建立的河冰動(dòng)力學(xué)模型為分析河冰的動(dòng)力演變過程提供了一種有效的數(shù)值方法。然而，現(xiàn)階段計(jì)算模型僅僅通過與理論分析和DynaRICE模擬結(jié)果對(duì)比的驗(yàn)證，還缺少與原型試驗(yàn)或現(xiàn)場(chǎng)實(shí)測(cè)數(shù)據(jù)的對(duì)比，這將在以后的工作中著重考慮。