（浙江大學 港口海岸與近海工程研究所，浙江舟山316000）

泥沙異重流常發(fā)生在水庫和高含沙河流等區(qū)域，如長江深水航道區(qū)浮泥沿河床運動形成的異重流現(xiàn)象[1].異重流在人類的生產生活中有著廣泛的應用，合理利用異重流排沙可以減小水庫淤積和航道淤塞.國內外學者通過野外觀測、實驗和數值模擬對異重流取得了豐碩的研究成果[2－5].特別是異重流數值模擬，隨著計算機科學的進步，取得了飛速的發(fā)展.層平均模型在異重流的數值模擬研究中得到了廣泛的應用[6－11].層平均模型不可避免地要引用泥沙侵蝕、水卷吸和泥沙沉速等經驗關系式.其中，水卷吸經驗關系式的精確度對數學模型的計算結果影響很大.這是因為，異重流的驅動力來自其與環(huán)境流體之間的密度差，而水卷吸經驗式的計算結果代表環(huán)境水體進入異重流的水量，若卷吸水較多，則使得含沙量較低，反之則含沙量高，并因此影響密度差和驅動力.

早在1959年，Ellison等[12]建立了異重流卷吸水量與理查德森數的經驗式.在此基礎上，范家驊[2]認為急流時泥沙異重流的水卷吸同鹽水異重流類似.Parker等[13]結合少量的持續(xù)入流式異重流水卷吸實驗數據，對Ellison經驗式進一步修正，修正公式在異重流層平均模擬中得到了廣泛應用.然而，異重流發(fā)生在水下，不易觀測到，且其與背景流體的交界面環(huán)境非常復雜（湍流、開爾文－赫姆霍茨（K－H）不穩(wěn)定性），要量化這個區(qū)域是個挑戰(zhàn)，需要獲取大量高精度野外及實驗數據來進一步率定水卷吸經驗式.當實測數據增加或存在誤差時，采用傳統(tǒng)最小二乘法擬合得到的經驗系數值極可能發(fā)生變化，這使得異重流水卷吸經驗式具有較大的不確定性.比如，Ottolenghi等[14]指出，Parker的水卷吸經驗式[15]可能會低估開閘式異重流水卷吸量一個量級.

Traer等[11]采用基于貝葉斯定理的 Monte Carlo方法（以下簡稱概率法），首次嘗試了對異重流水卷吸經驗式進行不確定性分析，并探討了其不確定性對異重流經典四方程模型的影響.Traer所采用的概率法相較于以最小二乘法為代表的傳統(tǒng)擬合方法有很大優(yōu)勢，因為通過經驗系數樣本與實測數據擬合優(yōu)劣的比較篩選過程，可以得到經驗系數的大量樣本，用以分析其均值和標準差等統(tǒng)計特性，而傳統(tǒng)擬合方法得到的是固定的經驗系數取值.然而，Traer所采用的異重流四方程模型具有恒定流和忽略地形沖淤過程的假設.一方面，恒定流的假設顯然不能全面地描述異重流隨時間的變化；另一方面，異重流運動過程中造成的地形沖淤，以及地形沖淤對異重流演化的反作用不能輕易忽略[9－10].Traer等采用四方程模型的可能原因是：異重流三方程模型（包括水沙混合體質量守恒、動量守恒和泥沙質量守恒）模擬自加速異重流時高估了泥沙侵蝕量，消耗過多能量，產能小于耗能從而不滿足能量平衡，與異重流自加速本質相反.然而，Hu等[16]通過理論推導和數值實驗發(fā)現(xiàn)，三方程模型并沒有破壞能量平衡.實際上，泥沙侵蝕量越多，異重流傳播時有越多勢能轉化為動能，增加產能，滿足能量平衡.

因此，本文考慮異重流非恒定特性、地形沖淤和地形沖淤對異重流演化的影響，研究異重流水卷吸經驗式對三方程模型框架下異重流層平均水沙耦合數學模型計算結果的影響.具體而言，本文結合現(xiàn)有異重流水卷吸實驗數據，應用基于貝葉斯定理的Gibbs隨機采樣方法，采樣獲取了異重流水卷吸經驗式中經驗系數組合的大量樣本；然后，將經驗系數樣本值代入異重流層平均水沙耦合異重流數學模型（以下簡稱耦合模型），模擬陡坡上異重流演化過程，從而討論經驗系數不確定性對模型計算結果的影響.

1 水卷吸經驗式不確定性研究方法

1.1 水卷吸經驗式的簡介

Parker等[13]基于異重流實驗數據，得到了水卷吸經驗式：

式中：ew為水卷吸系數；E1＝2.4，E2＝0.5，為Parker率定的經驗系數取值；Ri＝Rgφh／u2為理查德森數，R為有效重力，g為重力加速度，φ為層平均泥沙的體積分數，h為層平均厚度，u為層平均速度.

1.2 基于貝葉斯定理的Gibbs采樣法

用Data指代實測值（即Parker等[13]在可侵蝕底床上所做持續(xù)入流式異重流實驗中獲得的相關Ri和ew的數據，共53組），用m指代經驗系數組合，m＝ （E1，E2），用P（m｜Data）表征給定實測數據時，經驗系數某組取值的概率.根據貝葉斯定理[17]有

式中：下標j表示經驗系數樣本的序號；nT＝2×105為總樣本數；P（m）為先驗分布（即無實測數據條件下，經驗系數取值的概率）；P（Data｜m）為似然函數，即對于給定的經驗系數取值，式（1）計算值與實測值之間的擬合程度.由于每個經驗系數的取值受到實驗條件和測量數據準確性等很多因素的影響，符合中心極限定理，所以本文假設經驗系數取值符合正態(tài)分布，采用 Gibbs采樣法[18－19]，分別對E1和E2進行循環(huán)采樣.采樣過程敘述如下.

任意選取一組經驗系數值 （E1，0，E2，0）作為初始值；先固定E2，0不變，以E1，0為均值，以E1，0σE1為標準差（σE1為E1，0先驗概率的標準差），隨機取得一個新的備選系數樣本值，將其記為E＊1，則E＊1的接受概率[20]為

式中：f（E1，0｜E＊1）為以E＊1為均值，E＊1σE1為標準差的正態(tài)分布在E1，0處的值.當ra大于一個[0，1]之間的隨機數時，接受E＊1并以（E＊1，E2，0）為新的起點.固定E＊1，使用同樣的方法對E2，0進行采樣.如此循環(huán)，在采樣的過程中，不斷采取以下方法對樣本頻次進行統(tǒng)計：以每個樣本點為圓心，0.2為半徑作圓，統(tǒng)計圓內點的頻次.結果表明，當采樣所得的（E1，E2）樣本數達到2×105后，繼續(xù)增加樣本數量對樣本集合的頻次分布不再產生影響，本文認為采樣已經收斂，停止采樣.

Gibbs法的關聯(lián)性采樣體現(xiàn)在兩個方面：① 由于每個樣本的入選僅與上一個樣本相關聯(lián)，所以在采樣過程中，如果待選（E1，E2）樣本能夠比上一個樣本更好地擬合實測水卷吸數據（似然函數值較高），則該樣本能夠以更高的概率被接受；② 每組（E1，E2）樣本之間相關聯(lián)，如果一組（E1，E2）樣本使得水卷吸計算值和實驗值差距較大（似然函數?。?，則該組樣本的接受概率就小.基于上述認識可以預期的是，在通過足夠數量的采樣后，能夠使得水卷吸實驗值和計算值擬合越好的（E1，E2）樣本頻次越高.

2 經驗系數的樣本分布

得到樣本點在（E1，E2）平面的分布后，進一步得到25%nT、50%nT、75%nT和95%nT的（E1，E2）樣本所在的區(qū)間范圍，如圖1所示，其中黑色十字為概率最大樣本（E1，E2）＝（2.38，0.52）.顏色由深至淺，依次為25%nT、50%nT、75%nT、95%nT樣本的區(qū)間范圍.占總樣本95%的E1、E2取值組合中，E1和E2的區(qū)間范圍分別是[0.96Em1，2.72Em1]和[3×10－4Em2，6.59Em2]，Em1和Em2分別為概率最大的E1和E2值（即以該點為圓心，0.2為半徑的圓內，樣本頻次最高的值）.

圖1 E1和E2不同比例樣本所在的區(qū)間范圍Fig.1 Interval of different proportions of E1and E2

圖2顯示了概率最大系數組合與最小二乘法擬合結果（E1，E2）＝（2.38，0.52）的對比情況，圖2（a）為代入概率最大系數組合的式（1）與實測數據的擬合圖像，圖2（b）為原經驗式與實測數據的擬合圖像.可以看出，概率最大系數組合與原公式擬合結果相近.

圖2 概率法擬合的經驗式與原經驗式對比Fig.2 Comparison between the empirical formula fitted by probability method and the original formula

3 系數不確定性對數學模擬的影響

3.1 層平均水沙耦合異重流數學模型

在層平均假設下，非恒定的水沙耦合異重流數學模型可寫成如下守恒形式[9－10]：

式中：U為守恒變量向量；F為通量向量；Sb為幾何源項向量；Sf為底床摩擦、底床變形、水體卷吸等源項向量；z為底床高程；E＝ωEs和D＝ωφb分別為泥沙上揚通量和沉降通量，ω為泥沙沉速，Es＝8.9×10－9Z5m／（1＋8.9×10－9Z5m／0.3）為泥沙侵蝕系數[21]為泥沙雷諾數，R＝ （ρs－ρw）／ρ為有效重力，g為重力加速度，Ds為泥沙粒徑，ν為水的運動黏滯系數，φb＝r0φ為近底泥沙的體積分數，r0為近底與水深平均泥沙體積分數的比值；x為水平坐標；t為時間；g′＝Rgφ為水下有效重力加速度為摩阻流速，cD為拖曳系數，rw為異重流上下界面阻力比值；ρ＝ρw（1－φ）＋ρsφ為水沙混合體密度，ρ0＝ρwp＋ρs（1－p），ρw為清水密度，ρs為泥沙密度；p為床沙孔隙率.控制方程使用有限體積法進行離散[9－10]，并采用正方形網格，則有

式中：Δx為空間步長；Δt為時間步長；上標n和 ＊代表時間節(jié)點號；i代表空間節(jié)點號；Fi＋1／2和Fi－1／2為體積單元間的通量向量；Sb，i為底坡源項在第i個空間單元的計算值.時間步長需滿足CFL收斂條件.

本文將Parker推導的陡坡上異重流自加速臨界“點火條件”作為邊界條件，在這些初值條件下的異重流，從底床和環(huán)境水體獲得的能量高于其耗散，在沿程運動時會不斷加速，故稱為“自加速”異重流，這種異重流在實際生活中應用廣泛[8－9，11]：x＝0處，h＝2.0m，u＝0.801m／s，φ＝0.006 09.其他參數：g＝9.8m／s2，ω＝ 0.75cm／s，r0＝ 1.6[13]，cD＝0.004，Ds＝0.01cm，ν＝1×10－6m2／s，ρw＝1 000 kg／m3，ρs＝2 650kg／m3，p＝ 0.4，rw＝0.43.計算區(qū)域為6 000m，空間步長1m，時間步長1s，出口為自由出流.

3.2 （E1，E2）組合對異重流模型的影響分析

由于不同的（E1，E2）樣本值代表了不同的環(huán)境流體進入異重流的清水量，進而影響異重流與環(huán)境流體的密度差及其運動速度，為了使所有單個異重流算例均能流過最后一個統(tǒng)計斷面（1 000m），且盡可能提高計算效率，現(xiàn)隨機選擇5 000組樣本輸入模型并記錄其運動時間和運動距離.通過測試，當計算時間為2 200s，計算區(qū)域為6 000m時，所有算例均符合要求.

由于篇幅所限，現(xiàn)選擇概率最大經驗系數值所對應的異重流在4個典型時刻（100，200，300，400 s）的厚度h、速度u、泥沙體積分數φ、底床形變Δzb（初始時刻時為0），如圖3所示.其中，x表示與進口邊界的距離.

圖3 不同時刻異重流水力參數Fig.3 Turbidity current hydraulic parameters at different times

從圖3中可以看出：隨著異重流運動時間的增加，厚度和速度增加明顯，且增速在加快；泥沙體積分數先增加后減??；而底床被侵蝕得越來越多.

將所有（E1，E2）樣本值輸入上述異重流耦合模型，計算異重流過程，以x＝100m斷面為例，將計算所得水力參數范圍（0.53m≤h≤6.19m，0.98 m／s≤u≤2.30m／s，0.01≤φ≤0.11，－6.18m≤Δzb≤－2.91m）分為50等份，對每個區(qū)間作頻次（n）統(tǒng)計，結果如圖4所示.若計算水力參數的范圍越大，（E1，E2）不確定性對異重流計算結果的影響就越大.從圖4可以看出，厚度最小值頻次較高，速度和泥沙體積分數的較大值頻次高.比如，厚度區(qū)間0.53～0.64m的頻次為69 981，占比約37%，隨著厚度的增大，頻次逐漸降低，計算速度和泥沙體積分數有類似但相反的規(guī)律.底床形變的最小值和中間值（－4m左右）頻次都較高.－6.19～－6.12m區(qū)間的頻次為14 159，占比7%；－4.22～4.15m區(qū)間的頻次為33 947，占比18%；其他區(qū)間占比普遍低于3%.

為解釋這一現(xiàn)象，并了解樣本空間內不同樣本輸入模型計算所得的異重流水力參數分布情況，將E1樣本范圍[－3，7]以及E2樣本范圍[0，4]等分為100份，得到10 201組（E1，E2）.首先輸入耦合模型，得到10 201組計算結果，并取x＝100m斷面處計算結果作等值線，結果如圖5所示；然后輸入式（1）得到的10 201組ew計算值，作等值線，結果如圖6所示.

以圖5（a）為例，根據0.9m厚度等值線，將E1－E2平面分為兩個區(qū)域，區(qū)域I中ew較小，區(qū)域II中ew較大（圖6）.當（E1，E2）取區(qū)域I樣本值時（該區(qū)間樣本頻次為97 357，占比達51%），ew較小，計算的厚度相應較??；當（E1，E2）取區(qū)域II樣本值時（頻次93 214，占比達48%），ew較大，計算所得厚度較大，但是由于區(qū)域I的厚度范圍（0.53～0.9m）較區(qū)域II的厚度范圍（0.9～6.19m）小，所以頻率分布圖4（a）中出現(xiàn)了最小值頻次最高的現(xiàn)象.以異重流速度等值線圖5（b）為例，根據2.118m／s速度等值線，可將平面分為3個區(qū)域.當（E1，E2）取區(qū)域I樣本值時（頻次22 658，占比12%），速度區(qū)間范圍為2.113～2.118m／s；當（E1，E2）取區(qū)域II樣本值時（頻次80 235，占比42%），速度區(qū)間范圍為2.118～2.3m／s；當（E1，E2）取區(qū)域 III 樣本值時（頻次87 079，占比45%），速度區(qū)間范圍為0.98～2.113 m／s，雖然區(qū)間III的樣本頻次較高，但是其速度的范圍較大，所以圖4（b）中出現(xiàn)了較大速度頻次高的現(xiàn)象，泥沙體積分數和底床形變的頻率分布圖可以使用同樣的方法分析予以分析，不再贅述.

圖4 異重流水力參數的頻次統(tǒng)計（x＝100m）Fig.4 Frequency statistics of the hydraulic parameters（x＝100m）

圖5 E1－E2平面上計算異重流水力參數的等值線圖（x＝100m）Fig.5 Equivalent line of the hydraulic parameters on the E1－E2plane（x＝100m）

圖6 E1－E2 平面ew 的等值線與（E1，E2）樣本分布圖（Ri＝0.3）Fig.6 Equivalent line of ewon E1－E2plane and （E1，E2）sample distribution（Ri＝0.3）

將概率最大（E1，E2）取值輸入數學模型，計算所得各斷面水力參數分別記為hr、ur、φr和 Δzbr（以下簡稱標準值）.將（E1，E2）樣本值輸入數學模型計算所得水力參數（h、u、φ、Δzb）分別除以各斷面水力參數（hr、ur、φr、Δzbr），比較該倍數的變化情況，可以定量分析經驗系數不確定性對數學模型計算結果的影響.圖7給出了4個典型斷面（100，200，500，1 000m）上h／hr、u／ur、φ／φr、Δzb／Δzbr的范圍，并通過顏色的深淺區(qū)分了25%、50%、75%、95%樣本對應的該比值的范圍.圖中黑色虛線代表標準值隨空間的變化.

圖7 不同比例（E1，E2）樣本計算的異重流水力參數在4個斷面的區(qū)間范圍Fig.7 Interval of turbidity current hydraulic parameters calculated by different proportions of（E1，E2）

由圖7可知，隨著距進口距離的增大，計算厚度和速度的標準值迅速增大、泥沙體積分數標準值迅速減小、底床形變標準值變化不大.95%（E1，E2）樣本計算得到的比值（h／hr、u／ur、φ／φr、Δzb／Δzbr）范圍也隨著距離增大而增大.這說明，（E1，E2）不確定性對計算水力參數的影響隨著與進口邊界距離增大而增大.在x＝500m處，95%樣本計算的厚度是標準值的0.06～5.35倍、計算速度是標準值的0.49～1.24倍、計算泥沙體積分數是標準值的0.37～4.35倍，計算底床形變是標準值的0.1～1.11倍；到x＝1 000m處，計算水力參數與標準值之間倍數擴大至0.03～5.73（厚度）、0.46～1.23（速度）、0.39～5.31（泥沙體積分數）和0.01～1.08（底床形變）.由圖7可見，隨著距離的增加，95%樣本計算所得水力參數的范圍遠大于25%樣本計算所得水力參數范圍，這說明異重流數學模型對水卷吸經驗式中經驗系數（E1，E2）的取值十分敏感.比如，在x＝500m處，95%樣本計算的厚度是標準值的0.1～3.16倍，比25%樣本的范圍擴大了72%，計算速度、泥沙體積分數和床面形變也有類似現(xiàn)象.通過上述分析可以看出：經驗系數（E1，E2）對數學模型的計算結果影響很大，且隨著距離的增加，經驗系數不確定性對模型的影響在逐漸增大；概率最大（E1，E2）值代入數學模型可能低估異重流厚度和泥沙體積分數，高估異重流速度和床面形變.

4 結論

（1）本文應用基于貝葉斯定理的Gibbs隨機采樣法得到 概率最大 系數組 合 （E1，E2）＝ （2.38，0.52），與原公式的取值（E1，E2）＝（2.4，0.5）相近，該方法可以作為最小二乘法等傳統(tǒng)經驗式擬合方法的重要補充.

（2）比較（E1，E2）不同比例樣本計算結果在不同斷面的區(qū)間范圍，可以看出：相同比例樣本輸入模型計算得到異重流水力參數區(qū)間范圍隨運動距離的增加逐漸增大；經驗系數不確定性對數學模型的計算結果影響逐漸增大；采用概率最大經驗系數值可能低估異重流厚度和泥沙體積分數，高估速度和底床形變.

（3）比較不同比例（E1，E2）樣本輸入數學模型得到的同一斷面計算水力參數區(qū)間范圍可以發(fā)現(xiàn)，在x＝1 000m 處，95% （E1，E2）樣本計算得到的厚度范圍比25%樣本計算所得厚度范圍擴大了83%，說明模型對（E1，E2）系數組合的取值十分敏感.

由于概率最大經驗系數取值輸入模型預測的異重流水力參數可能出現(xiàn)偏差，所以下一步將在兩個方面開展工作：① 增加實驗數據獲取，以便對經驗系數的取值進行更加準確的率定；② 開展不同水卷吸經驗式不確定性對數學模型影響的比較，為經驗式選擇提供依據.