1、二、泥沙模型（一）、懸移質(zhì)泥沙數(shù)學模型目前，國內(nèi)外泥沙數(shù)模的研制已有了相當大的進展，二維數(shù)模的理論研究已 經(jīng)成熟，以前鑒于計算速度較慢和理論上的復雜性，有限元法在水流和泥沙耦合 計算中的應(yīng)用在國內(nèi)一直處于次要的地位，因此，發(fā)展適用于強對流、高含沙水 流的高精度有限元泥沙模型，具有一定的現(xiàn)實意義。我們知道，對于大多數(shù)水流 模型，在緩流情況下均能得到非常精確的解，當水流出現(xiàn)局部急流的情況下，只 需在該處稍加一些粘滯性，即可將其抑制，并且它對整體的求解精度并無太大的 影響；而在以擴散為主的傳質(zhì)模型中，我們也可以通過增大擴散系數(shù)和滯后因 子來避免出現(xiàn)濃度分布的不合理現(xiàn)象。但是，泥沙在河流中的運動機理非

2、常復雜， 一般情況下，無論是急流還是緩流，懸移質(zhì)的對流作用始終處處占優(yōu)，采用傳統(tǒng) 的有限元方法或者再通過增大泥沙擴散系數(shù)來保證計算結(jié)果合理是行不通的。同 時，由于它還同床面發(fā)生交換從而改變水流條件，因而泥沙輸運問題在理論上講 始終是一個非恒定問題。泥沙模型的計算除了借鑒水質(zhì)模型中的數(shù)學方法外，還 要對床面進行不斷的修正來重新計算流場，通常的外部或人工的耦合求解方法只 能適用于流速和含沙量的變化率均較小的情況。為此，本文主要研究迎風及特征 有限元法在泥沙輸移計算中的應(yīng)用，編制相應(yīng)的程序并進行算例的對比，從而得 出合理的結(jié)論。2基本方程水流控制方程如下其中u，v為流速；h，z為水深和床面高程；E為

3、紊動粘滯系數(shù)；C為謝才 系數(shù)，表達式為：C=1/nhi/6b泥沙輸移控制方程如下河床變形方程其中，為泥沙擴散系數(shù)(m2/s), Y 為淤積物十容重(kg/m3), 3為顆粒 沉速(m/s)，s和s*分別為垂線平均含沙量和水流挾沙力(kg/m3)，a為恢復飽和 系數(shù)。為簡化計算，本文采用如下假設(shè)：(1)懸移質(zhì)為均質(zhì)沙并且床面有足夠數(shù)量的床沙質(zhì)可與其交換。(2)泥沙擴散系數(shù)取一常數(shù)，挾沙力采用一維水流挾沙力公式：s*=k(U3/ghw)m。、平面二維非均勻沙數(shù)值模擬方法水流泥沙數(shù)學模型作為河流模擬的一個分枝，是定量預(yù)測水沙運動及河床演 變的重要手段.研究出一種性能良好的平面二維水沙數(shù)學模型，其著眼

4、點主要有 二1： 一是計算參數(shù)的選取及泥沙基本理論的完善；二是準確、可靠而又迅捷的 計算方法.本文力求從以上兩個方面入手，在平面二維水沙計算方法上做一些創(chuàng) 新，在模擬的內(nèi)容方面作一些擴展.1基本方程平面二維水沙數(shù)學模型由以下方程控制水流連續(xù)方程(1)水流運動方程如 du 如 dzzjJ震以+撰 襯& 目口S +一4g + g-r/+) = 0(2)dy 3ady(2)泥沙連續(xù)方a泥沙連續(xù)方a標 咨 溢程懸移質(zhì)河床變形方程”與床沙組成方程式中Z為水位；Z為河底高程；h為水深，h=Z-Z ； u、v為垂線平均流速沿x、y 方向的分量；C為謝才系數(shù)，C=(1/n)hi/6 n為糙率系數(shù)；Y為紊動粘性

5、系數(shù)； %、s*為第k組垂線平均含沙量及挾沙力；3 k為泥沙沉速；！為泥沙擴散系數(shù); a為*恢復飽和系數(shù)；Y 為泥沙十容重；Pk、Pok分別為混合層及天然河床床沙組 成，E為混合層厚度，和為標記，純淤積計算時 =0,否則， =1, 當混合層下邊界波及到原1 始河床時 2=1，否則， 2 = 0；1k為非均勻沙分組系、黃河河口海岸二維非恒定水流泥沙數(shù)學模型自本世紀七十年代以來，由于計算機技術(shù)的迅猛發(fā)展，國內(nèi)外相繼出現(xiàn)了眾 多的河口海岸泥沙數(shù)學模型IE，有力地促進了河口海岸的泥沙研究的發(fā)展。然 而，已有的河口海岸數(shù)學模型大多是模擬含沙量較低的河口海岸的泥沙運動，而 能夠模擬多沙和岸線延伸劇烈的河口

6、海岸泥沙數(shù)學模型還極為少見。近年來，已 有個別學者嘗試用泥沙數(shù)學模型模擬黃河河口海岸的泥沙運動，如張世奇開發(fā)了 一套黃河口平面二維泥沙沖淤數(shù)學模型，得到了較好的效果1820。為了全面系統(tǒng) 地反映黃河三角洲海陸動態(tài)交互影響機理和泥沙運動與濕地演替關(guān)系，本文開發(fā) 和建立了徑流、潮流和波浪作用下的黃河河口海岸平面二維動邊界非恒定流非均 勻沙不平衡輸沙數(shù)學模型。1模型結(jié)構(gòu)1.1水流運動基本方程1.1水流運動基本方程(2)式中：U、V分別為潮流速在x及y方向的垂線平均值分量；Z為潮位；Cf為謝才 系數(shù)；F為柯氏系數(shù)，F(xiàn)=23 sin中，式中3為自轉(zhuǎn)角速度，中為地理緯度；h 為水深；z為海底起始高程。1.

7、2潮流和波浪共同作用下泥沙運動方程 在河口海岸地區(qū)，潮流和波浪是泥沙 運動的最主要動力。“波浪掀沙和潮流輸沙”使河口海岸地區(qū)的泥沙運動極為活 躍，但也更為復雜。因此，在計算和預(yù)報河口海岸地區(qū)的泥沙運動和沖淤變化時， 僅僅考慮潮流的作用是不全面的，還應(yīng)考慮波浪的作用。二維非恒定流不平衡情況下懸移質(zhì)含沙量隨空間和時間的變化規(guī)律的偏微 分方程式(考慮泥沙粒徑不均勻性時，以下各方程式中含沙量和沉速是相應(yīng)于某 一粒徑組的) TOC o 1-5 h z 式中：l表示非均勻沙中的第l組泥沙，h為水深；U、V分別x、y方向垂線平 均流速；*為第l組泥沙垂線平均含沙量；S*為第l組泥沙水流挾沙力；3 為 第l組

8、泥沙沉速； s為泥沙擴散系數(shù)；a為恢復飽和系數(shù)。1為方便起見，在以下分析中省略下角標l.設(shè)U、U和V、V分別表示潮流速 和波浪質(zhì)點速度在x和y方向上的分量，則在流場任意點x和yW方向上瞬時流速 為：U=U +U(5)c wV=Vc+Vw(6)將上述二式代入式(4)即得到潮流和波浪共同作用下的懸移質(zhì)泥沙不平衡輸沙方程式如下：波浪對泥沙的作用在一個周期內(nèi)為往返輸移，其單向輸沙作用一般只表現(xiàn)在 半個周期內(nèi)。由于波浪的周期比較小，一般只有幾秒至十幾秒；而潮流的周期是 半天或一天。因此，當觀察的時間或計算的步長大于波浪周期時，就可以對上式 進行波浪周期的平均。由于潮流速與波浪質(zhì)點速度基本無關(guān)，而波浪質(zhì)點

9、速度在 一個波周期內(nèi)的平均值基本為零，對式(5)在波浪周期內(nèi)取時均后簡化為：率十竺件十竺羿2十四di dx 莎 I & j式(8)在形式上與式(4)相同，但差別在于式(8)中的挾沙能力S為潮流和波浪共 TOC o 1-5 h z 同作用下的挾沙能力。*由于懸移質(zhì)的不平衡輸沙運動引起海底地形的沖淤變化，其河床變形方程 為： HYPERLINK l bookmark22 o Current Document 空+迎*竺竺dl dx dy dxdi式中：y 為海底淤積泥沙的十容重，n為海底沖淤厚度。通過上述分析可見，在計算河口海岸大范圍泥沙沖淤變化時，只需計算潮流 流速場，無需計算波浪流速場，只是在

10、計算挾沙能力時考慮波高因子的影響，這 與竇國仁的分析結(jié)果是一致的。(四)、紊流封閉問題現(xiàn)有的紊流計算方法主要是基于渦粘性假設(shè)，渦粘性概念將紊流封閉問題轉(zhuǎn) 化為如何確定渦粘性在流場中的分布問題。有多種平均和模型方法可用于計算渦 粘性，譬如，若對基本方程進行雷諾時均，可采用任何一種紊流模型；若對基本 方程進行空間平均或濾波處理，可采用亞格子尺度模型(Sub grid scale model) o然而，對于某一具體問題，選用何種程度的紊流封閉模型尚無統(tǒng)一的標 準。以紊動輸運項在基本方程中的相對重要性作為選擇的物理依據(jù)，且綜合考慮 精度、效率及復雜性等實際因素的做法是合理的。對于很多大范圍的淺水流動問

11、題，動量方程的慣性項主要由壓力梯度項或浮 力項平衡，是否應(yīng)用精細的紊流模型就顯得并不重要。因此，出于簡化計算的目 的，很多流動和輸運三維模型6,7,11,12,15,2225或者使用一個經(jīng)驗系數(shù)或者使用一個 簡單紊流模型。此外，由于動量方程中的水平渦粘性項與垂向渦粘性項相比通常 是小量，且其影響常為數(shù)值耗散所掩蓋，很多淺水模型6,22或者忽略水平擴散項 或者簡單地使用一個恒定的水平渦粘性系數(shù)。然而，由于該二階微分項在數(shù)值計 算中具有光滑效應(yīng)，能夠增加數(shù)值計算的穩(wěn)定性，在深水區(qū)，特別當離散網(wǎng)格尺 寸較大時，或者當計算域的地形變化劇烈時，保留水平渦粘性項是必要的。當流動比較復雜或者需要考慮物質(zhì)輸運

12、時，使用精細的紊流模型就顯得非常 必要。K % 雙方程模型是最簡便的精細紊流模型，也是目前應(yīng)用最廣泛的紊流 模型，并在相當廣泛的范圍內(nèi)得到了充分的檢驗，證明是實用和有效的，很多三 維流動和輸運模型8,13,14,26,27也都使用了該模型。然而，標準的K -模型假定渦 粘性系數(shù)是各向同性的標量，即水平渦粘性與垂向渦粘性相等，而淺水流動通常 表現(xiàn)出垂向紊動尺度小于水深和水平紊動尺度遠大于水深的紊動結(jié)構(gòu)，這種雙重 紊動結(jié)構(gòu)以及強烈的各向異性導致水平和垂向擴散系數(shù)存在較大差異，因此標準 的單尺度(single length scale)K -模型將低估水平渦粘性系數(shù)。應(yīng)力-通量 方程模型、應(yīng)力-代數(shù)

13、模型以及兩尺度(two-scale) K -模型克服了單尺度 K -模型的這一缺陷。但由于計算繁瑣，較少應(yīng)用到大尺度水體的計算中，而 且基于雷諾平均的紊流模型無法模擬水平大尺度漩渦的動態(tài)演化。研究表明，紊流是由不同尺度的旋渦組成，大渦對紊流能量和雷諾應(yīng)力的產(chǎn) 生以及對各種量的紊動擴散起主要作用，其行為強烈地依賴于邊界條件，隨流動 類型而異；小渦主要對能量耗散起作用，在高雷諾數(shù)下近似于各向同性，受邊界 影響較小。因此，用雷諾平均的方法并不能真實地反映紊流的上述特點，因為此 時的雷諾應(yīng)力項包含了全部大渦引起的脈動。更合理的方法是把紊流的瞬時運動 通過某種濾波函數(shù)分解為大尺度運動和小尺度運動，大尺度

14、運動由基本方程直接 數(shù)值求解，小尺度運動對大尺度運動的影響在動量方程中表示為亞格子雷諾應(yīng) 力，使用亞格子模型模擬小渦。亞格子模型的主要優(yōu)點是大渦部分可以精確求解， 計算受人工模型的影響較小。然而，亞格子模型要求網(wǎng)格足夠細，導致計算量大 增，解決大尺度的實際工程問題尚有困難。實際的潮流的計算中，常建立垂向渦粘性與流速或剪切應(yīng)力之間的相關(guān)關(guān)系 ，水平渦粘性系數(shù)則由Samagorinsky模型28計算。試驗顯示，密度分層對水 平紊動擴散系數(shù)影響很小，可忽略不計，但垂向擴散系數(shù)需要由Munk-Anderson 公式修正。5數(shù)值計算方法淺水流動及輸運三維數(shù)值計算方法在過去二十多年的時間里取得了巨大的 進

15、展：在空間離散上，有差分法、有限元法和有限體積法；在適應(yīng)物理域的復雜 幾何形狀上，有貼體坐標變換及。坐標變換；在時間積分上，有顯式、隱式、半 隱格式以及時間分步；在求解技術(shù)上，有ADI法、迭代法、多重網(wǎng)格法以及并行 計算技術(shù)；在十濕、露灘動邊界的處理上，有固定網(wǎng)格和動態(tài)網(wǎng)格技術(shù)；在懸沙 輸運重力沉降項的處理上，有源項化和對流化的做法；在提高精度上，有歐拉- 拉格朗日法和差分有限元雜交法等。限于篇幅，本文只給出了幾種常用的比較有 效的計算方法，并分析它們的優(yōu)缺點。5.1時間分步法 時間分步法的基本思想是引入一個或多個中間變量，在每一時 間步內(nèi)把對時間的積分分解成兩個或多個子時間過程，象著名的AD

16、I差分格式、 預(yù)測-校正格式以及高階Runge-Kutta法。常用的時間分步法可分為基于空間概 念上的分步和基于物理概念上的分步，按照坐標軸方向?qū)刂品匠踢M行的分步是 空間概念上的分步法，象廣泛使用的局部一維化處理方法；根據(jù)方程中各項的物 理意義進行的分步則屬于物理概念上的分步法。空間分步的優(yōu)點是能夠帶來計算 的簡化，省時省內(nèi)存，文獻11對懸沙輸移方程分解成平面二維和垂向一維兩部 分進行求解。經(jīng)物理概念上的分步后，可以對方程的各項選用最適合的離散方法， 文獻15,29,30將控制方程分解為對流、擴散和波動三個不同的物理過程。文獻 31提出的三步有限元則屬于分步法在有限元中的應(yīng)用。5.2模式分裂

17、技術(shù) 鑒于具有自由表面的淺水流動在物理上表現(xiàn)為快速傳播的表 面重力波和緩行的內(nèi)重力波的兩重屬性，可以將模型方程相應(yīng)的分裂為內(nèi)模式和 外模式兩部分，此乃模式分裂的物理依據(jù)。內(nèi)模式方程為原始的三維方程，主要 代表動量的垂向交換過程，外模式方程通過沿水深積分內(nèi)模式方程獲得，內(nèi)外模 式耦合計算。盡管求解外模式要求較小的時間步長，但是由于不受表面重力波波 速的限制，內(nèi)模式的三維計算可以使用較大的時間步長，從而獲得較高的計算效 率，這就是進行模式分裂的目的和結(jié)果。著名的POM模型、CH3D系列模型就采 用了該技術(shù)，此外很多淺水流動和輸運模型也使用了該技術(shù)6,7,14,29。需要注意的 是，在進行模式分裂時

18、，必須確保內(nèi)外模式各物理量及底部剪切應(yīng)力的相容性， 否則可能導致數(shù)值計算的不收斂。5.3半隱格式 由于允許使用較大的時間步長且無需求解大型稀疏矩陣方程，同 全顯、全隱的差分格式相比，ADI差分格式帶來了計算效率的提高。然而，ADI 格式使用交替方向隱式求解，不能恰當考慮不同方向流動之間的相互作用，數(shù)值 解存在流速向量向某坐標軸方向偏斜的一維化趨勢，即所謂的ADI效應(yīng)32。當流 向與坐標軸夾角為45時誤差最大，且誤差隨時間步長的增大而增大，盡管限 制時間步長可以減小ADI效應(yīng)，但無疑降低了 ADI格式的效率。為克服上述ADI格式及模式分裂技術(shù)等的不足，文獻23提出了一種求解三 維淺水方程的半隱差

19、分格式。通過在每一時間步內(nèi)對動量方程中的壓力梯度項和 垂向渦粘性項以及自由水面方程的流速作隱格式差分離散，對動量方程中的對流 項、水平渦粘性項以及柯氏力項使用歐拉-拉格朗日法離散，既克服了顯式離散 的穩(wěn)定性條件受表面重力波波速和垂向擴散項限制的缺點，又避免了隱式離散需 要求解大型稀疏矩陣方程的困難，并較之歐拉法具有較高的離散精度。得益于上 述離散格式，可以采取如下的求解方法。先對由水位方程離散得到的線性、對稱、 正定、五對角且對角占優(yōu)的矩陣方程應(yīng)用共軛梯度法求解得到水位；然后對由動 量方程離散得到的線性、三對角且對角占優(yōu)的矩陣方程應(yīng)用直接求法解得到水平 流速；最后由連續(xù)方程積分得到垂向流速。穩(wěn)

20、定性分析及計算實例表明4,23，該格 式較之傳統(tǒng)的全顯、全隱差分格式有顯著的效率優(yōu)勢，在作者的San Francisco 灣三維潮流計算中，最大Courant數(shù)超過了 50，平均Courant數(shù)為20。此后， 此半隱差分格式被進一步應(yīng)用于密度分層和非靜壓的三維自由表面流的計算中 撲。Stansby將該半隱思想應(yīng)用到o坐標和有限體積法中修如。5.4有限元離散 為克服結(jié)構(gòu)化網(wǎng)格(Structured grid)對復雜邊界處理的困難， 計算中常采用貼體曲線坐標。盡管近年來各種代數(shù)與微分方程曲線坐標網(wǎng)格生成 技術(shù)取得了很大進展，但對于非常復雜的實際工程問題，特別是三維問題，曲線 坐標網(wǎng)格的生成仍然是相

21、當困難和費時的。而且坐標變換使控制方程復雜化，致 使貼體坐標變換在實際應(yīng)用中受到限制。采用非結(jié)構(gòu)化網(wǎng)格(Unstructured grid) 是處理復雜幾何邊界問題的重要方向。與差分法不同，有限元是積分形式的離散方法。與基于變分原理的經(jīng)典的結(jié) 構(gòu)分析中的有限元不同，現(xiàn)代有限元建立在函數(shù)分析或優(yōu)化理論等數(shù)學基礎(chǔ)之 上，可以看作是廣義加權(quán)余量法的控制方程弱解的統(tǒng)一形式，若選取特定的權(quán)函 數(shù)，有限元法可以退化為有限差分、有限體積等方法。由于具有任意局部加密、 邊界適應(yīng)能力強、通用性強及精度高等優(yōu)點，有限元在固體力學中獲得了巨大的 成功。然而，傳統(tǒng)的Galerkin法等價于中心差分格式，隨著雷諾數(shù)的增

22、大，N-S 方程不僅呈現(xiàn)高度的非線性而且方程的性質(zhì)也會發(fā)生退化，呈現(xiàn)混合型方程的性 質(zhì)，簡單的采用Galerkin法將導致穩(wěn)定和收斂方面的問題，這是阻礙有限元在 計算流體力學領(lǐng)域廣泛應(yīng)用的一個主要原因。隨著近年來諸如流線迎風、 Taylor-Galerkin、分步31 及高分辨率34等有限元格式的提出，有限元已經(jīng)能 夠適用于對流占優(yōu)的高雷諾數(shù)流動。此外，現(xiàn)代迭代技術(shù)使用共軛梯度法求解隱 式有限元矩陣方程也是非常高效的，且隨著并行計算的發(fā)展，顯格式結(jié)合多重網(wǎng) 格加速的有限元也將是更加高效的。文獻9給出了有限元在國外三維潮流模擬中的應(yīng)用情況。文獻35使用 Petro-Galerkin有限元法建立了河流及河口泥沙輸移三維數(shù)學模型。文獻 25,36使用半隱有限元格式對三維淺水及分層流動進行了計算。Zienkiewicz 發(fā)展的基于特征分裂的有限元淺水方程模型具有較高的計算效率。Kashiyama38 在并行計算機上實現(xiàn)了潮流和風暴潮的數(shù)值模擬。6