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

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

1 研究背景

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

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

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

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

2 離散元方法與水動力學耦合的河冰動力學模型

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

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

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

圖1 河冰動力過程的受力分析

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

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

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

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

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

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

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

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

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

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

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

河水流量分為冰層流量和冰下流量兩部分，即qt=qu+ql。冰層中的水流一方面會隨著冰一起運動，一方面也會因為河流表面的水平坡度而發(fā)生滲流。

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

式中：qi為隨冰塊一起運動的河水流量；vi為冰量；qs=VsΔH為在冰層中運動的滲流流量；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 擴展多面體離散元方法與水動力學的耦合模型河冰動力過程中的強水動力作用使得冰與水流間的相互作用不可忽略。本文建立的河冰動力學模型對水和冰的耦合處理通過以下幾個方面進行：在水力模型方程中考慮表面冰對水流的影響，基于擴展多面體離散單元與三角形有限單元的相對位置通過插值平均進行信息交換，水對冰的作用通過對離散單元上的受力計算完成。

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

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

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

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

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

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

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

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

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

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

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

在水動力學有限元計算中，臨界時間步長為：

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

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

圖5 計算耦合時間步

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

將以上河冰動力學模型應用到規(guī)則河道進行簡單工況的數(shù)值模擬，并通過與DynaRICE的計算結(jié)果和理論分析結(jié)果對比來檢驗模型的有效性與準確性。DynaRICE 模型由美國Clarkson大學沈洪道教授建立，已在多條河流動力學模擬中得到驗證［14，17］并在國際上被廣泛認可。

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

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

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

式中：Q為水文斷面的流量，m3/s；Ac為過水斷面的面積，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的綜合糙率。綜合糙率的計算公式形式多樣，被廣泛應用的有Belokon-Sabaneev 公式、Einstein 公式和Larsen 公式［44］。由于在模擬中采用的是寬深比較小的矩形明渠，因此這里采用誤差較小的Einstein公式計算綜合糙率，即：

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

這里規(guī)則矩形河道的長度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，主要計算參數(shù)列于表1。圖8為均勻恒定明渠流初始時刻的縱向水力參數(shù)。在明渠流的穩(wěn)定狀態(tài)下，沿河道方向各處流速基本穩(wěn)定在0.55 m/s，水深穩(wěn)定在0.59 m。

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

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

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

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

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

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

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

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

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

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

河冰的動力過程往往伴隨著水力條件的快速變化，特別是在冰塊堆積、堵塞河道并形成冰壩時，過水斷面急劇減小，上游水位迅速升高。通過對河冰輸移、堆積、阻塞形成冰壩過程的數(shù)值模擬，檢驗本文河冰動力學模型的有效性。

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)論

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