水流作用下地形不穩(wěn)定性計算模型

王明明，鄒志利

（大連理工大學海岸和近海工程國家重點實驗室，大連 116024）

水流作用下地形不穩(wěn)定性計算模型

王明明，鄒志利

（大連理工大學海岸和近海工程國家重點實驗室，大連 116024）

河流中水流和水底泥沙之間的相互作用會產生各種復雜的地形形態(tài)，交替型沙壩是其中的一種地貌特征，研究表明這一地貌的形成機理是水底地形的不穩(wěn)定性。給出了這一不穩(wěn)定性所導致的交替型沙壩的數值模擬方法和相應結果。水流方程采用了交替方向隱式算法求解，水底變形方程采用了比傳統(tǒng)方法更加穩(wěn)定的歐拉加權基本無振蕩迎風格式求解，從而建立了水流與水底地形變化耦合計算模型。運用該模型研究了初始地形擾動產生的交替沙壩在長直河道中的演變。

地形演變；非線性；WENO；不穩(wěn)定性

Biography：WANG Ming-ming（1986-），male，master student.

近岸區(qū)域和河流中地形演變依賴于水流、波浪、泥沙之間復雜的相互作用，水流作用下地形的不穩(wěn)定性是這一復雜相互作用的突出表現形式。這一不穩(wěn)定性導致水底地形即使是在初始時刻是均勻平直的水流作用下也會產生沿空間和時間變化的地貌特征，交替型沙壩就是這種復雜地貌特征的常見典型例子。在求解這種地形和水動力相互作用所引起的水底不穩(wěn)定性問題時，如何求解地形演變方程是地形變化數值模擬中的一個難題，因為該方程是非線性的，會導致間斷解的出現。以往多采用Lax-Wendroff格式［1］或者其改進后的格式，這些格式在地形峰值后面會產生數值振蕩，從而不能夠模擬地形的長時間演變過程。Hundson［2］對比了Lax-Friedrichs格式、Lax-Wendroff格式和MacCormack格式，提出更加穩(wěn)定的Roe格式，但是很難應用于二維模型。加權基本無振蕩（Weighted Essentially Non-Oscillatory）格式是1994年由Liu和Osher等［3］在ENO（Essentially Non-Oscillatory）構造思想的基礎上提出的。ENO［4］格式具有在間斷區(qū)分辨率高等優(yōu)點。為了進一步提高該格式的計算精度，WENO格式引入了變化的加權因子，使計算結果更加穩(wěn)定［5］。文章研究了地形和水動力相互作用所引起的水底地形不穩(wěn)定性中WENO格式的應用問題。對水流方程仍然采用了常用的ADI方法［6］，但對地形演變方程采用了Euler-WENO方法，在此基礎上建立了水流與水底地形變化耦合計算模型，并與傳統(tǒng)的FTCS離散格式建立的模型的計算結果進行對比。應用該模型對交替沙壩地形的二維長直河道的形成進行了數值模擬，給出了這種地形形成的不穩(wěn)定性機理的數值模擬結果。

圖1 坐標系統(tǒng)Fig.1 Coordinate system

1 水流和水底地形的耦合數學模型

考慮一有限寬度的長直河道，取x軸和y軸分別沿河道的長度方向和寬度方向。水底沿x方向有一個微小的傾角i0，初始水流沿整個河道是均勻的。水流和地形變化滿足以下基本方程。

式中：u=（u，v）為水流流速矢量；η為水面升高，即自由水面與靜水面之間的距離；H為總水深，H=η+h*-zb，h*為靜水水深，zb為床面高程；Lx，Ly分別為河道的長度和寬度；f為科氏力系數（f=2ωsinφ，ω為地球自轉速度，φ為河流所處的當地緯度）；Cz=H1/6/nr為謝才系數，nr為水底糙率；g為重力加速度；vo是與泥沙粒徑相關的參數，一般取經驗值。若將空隙率np考慮在內，可以用v來代替vo；q=（qx，qy）為推移質體積輸沙率；b取值范圍通常在2～7；uc為推移質的起動流速；γ′▽zb為坡度項，表示泥沙順坡向下滑動的趨勢，坡度項能使計算得到的沙丘形狀更符合實際，γ′通常取1。這里泥沙輸移只考慮推移質輸沙。

文中對水動力方程（1）～（3）采用了ADI（交替方向隱式）法，網格劃分如圖2所示。計算時將每個時間步（t→t+Δt）的計算分成兩步：在t→t+Δt/2，連續(xù)方程和x方向的動量方程采用半隱式格式離散，而y方向的動量方程采用顯式格式離散；在t+Δt/2→t+Δt，連續(xù)方程和y方向的動量方程采用半隱式格式離散，而x方向的動量方程采用顯式格式離散。在每個時間分步，可形成以η為變量的三對角方程組，從而可采用追趕法快速求解。下面僅給出前半步t→t+Δt/2的離散過程，后半步的離散過程與此類似。

vn+1/2的顯示求解：

圖2 交錯網格布置圖Fig.2 Disposition of variables in staggered grid

un+1/2的半隱式求解：

x方向的動量方程采用半隱式離散，其離散格式為

對水底變形方程（4），最簡單的離散格式為FTCS（forward in time，central in space），即時間上向前差分，空間上采用中心差分。具體離散格式為

FTCS格式雖然形式簡單易懂，但對于求解水底變形方程非常不穩(wěn)定（第2節(jié)和第4節(jié)中詳細給出了FTCS格式在一維和二維地形演變計算中出現的不穩(wěn)定情況）。為了克服該格式的不足，本文除了采用該格式外，還采用如第1節(jié)所述的加權基本無振蕩WENO格式。該格式具有在間斷區(qū)分辨率高且具有五階精度［8］等特點，可以為建立一個穩(wěn)定精確的水流與水底地形變化耦合計算模型提供基礎。下面以一維情況為例，介紹WENO方法應用的具體步驟［9］。

將輸沙率q在x正方向和負方向分解成與地形傳播速度相關的兩部分

二維水底變形方程的Euler-WENO離散格式與以上處理完全類似，即將x方向和y方向的輸沙率（qx，qy）都按上面計算q的方法計算（式（15）中各量也同樣適用）。所對應的差分方程為

在第2節(jié)和第4節(jié)中將具體討論Euler-WENO格式在計算地形演變時的精度和穩(wěn)定性。

2 模型驗證

2.1 水動力模型驗證

取一平底河道，長1 320 m，寬220 m，水深5 m，坡度i0=6.11×10-5，將坐標系x取在水底平面上，z軸垂直水底向上。流速為1.0 m/s穩(wěn)定水流。在渠道兩端處施加出流邊界條件。由于沿y方向地形沒有變化，不考慮科氏力影響，因此由x方向動量方程可得

將 u=1.0 m/s，H=5.0 m，Cz=H1/6/nr代入式（18），解得糙度 nr=0.022 856。

圖3 數值模擬的流速Fig.3 Simulation of velocity

將nr=0.022 856代入數值模型，模擬結果如圖3所示。流場在2 h左右達到穩(wěn)定，得到的流速為0.999996m/s，水動力數值模擬結果與解析解吻合良好。

2.2 地形演變模型驗證

Hudson給出了在長1 000 m、水深10 m矩形水槽中高1.0 m的沙丘在單向水流作用下的演變解析解，本節(jié)將分別采用FTCS方法和Euler-WENO方法模擬該問題并與Hudson得到的解析解進行比較。計算時所采用的地形、輸沙率公式和其他參數同Hudson的一致。

初始地形為

地形方程為方程（4）的一維形式，輸沙率公式為 q=aub，a=0.001 s2/m，b=3.0，np=0.4。

圖4和圖5分別給出了4個時刻的FTCS方法和Euler-WENO方法計算結果與Hudson的解析解的對比。圖4和圖5中從左到右依次為t=0 s（初始地形）、t=93 035 s、t=186 070 s、t=238 079 s時的地形。

圖4 FTCS格式計算結果與解析解的對比Fig.4 Comparison between FTCS results and analytical solutions

圖5 Euler格式計算結果與解析解的對比Fig.5 Comparison between Euler-WENO results and analytical solutions

由圖4可知，FTCS格式在t=93 035 s和t=186 070 s2個時刻的計算結果與解析解吻合；但是隨著計算時間增長，在t=238 079 s時，沙丘頂部出現數值振蕩，與解析解不完全吻合；由圖5可知，Euler-WENO格式的計算結果在3個時刻都能與解析解吻合，并且未出現類似FTCS格式的不穩(wěn)定振蕩。這些結果表明，Euler-WENO格式在求解水底變形方程時計算精度較高，且長時間計算過程中穩(wěn)定性優(yōu)于FTCS格式，因此適合用來模擬計算水底演變。

3 河流水底地形不穩(wěn)定性

目前研究表明，雖然河流中水流是平直的，但水底地形不能保持平面形態(tài)，受到微小地形變化的擾動后，產生沙紋和較大尺寸的交替型沙壩等復雜地形形態(tài)，這一水底變形特征即為河流地形的不穩(wěn)定性。Colombini［10］對長直渠道中的交替沙壩的演變進行了線性穩(wěn)定性分析，表明非線性作用能夠抑制地形波動峰值的無限增長，從而使波動達到一個穩(wěn)定的峰值。Schielen［11］在Colombini的基礎上得到了一個用于描述不穩(wěn)定交替沙壩群包絡振幅非線性演變的Ginzburg-Landau方程，結果表明，周期性的沙壩樣式會變得不穩(wěn)定，并表現出準周期性的變化特征。本節(jié)給出有關不穩(wěn)定性線性化理論分析結果，為下一節(jié)將通過數值模擬給出非線性的地形不穩(wěn)定性數值模擬結果提供基礎。

將（21）、（22）、（23）三式代入模型方程（1）～（5），假設擾動為小量，將控制方程線性化，得到關于擾動量的控制方程。假定擾動量的解的形式為

式中：ω=ωr+iωi，ωr為地形的波動頻率，ωi為增長率，ωi＞0 表示地形擾動隨時間增長，地形是不穩(wěn)定的；k 為擾動波數。將式（24）代入控制擾動量的控制方程，可得到ω和k的關系，其數值結果如圖6所示（無因次結果）。

圖6 ω和k的關系曲線Fig.6 Relation curves ofωand k

圖 6 中無因次波數 k′=2π/λ′，λ′=λx/Ly為無因次的地形擾動波長。由圖 6可見，k′在一定范圍內（0.19～3.32），有ωi＞0，即地形在這些波數情況下會變得不穩(wěn)定；在k=1.05時，ωi的值達到最大。下一節(jié)將針對這一波數情況對正體擾動的演化進行數值模擬，以了解地形擾動的演化過程和將正體上節(jié)所介紹的數值方法應用于數值模擬這一演化過程。

4 非線性水底地形不穩(wěn)定數值模擬研究

由第2節(jié)知，水底變形方程（5）的非線性特征（方程中輸沙率對地形變化zb存在著非線性依賴關系）要求其求解的計算格式具有較高的分辨率，傳統(tǒng)的計算格式（如FTCS方法），雖方法簡單，但卻存在數值振蕩（圖5），所以不能用于地形演變的計算研究，而Euler-WENO方法可以克服這些不足。下面將針對二維水底地形演變來繼續(xù)探討這一問題。計算針對上節(jié)給出的最大不穩(wěn)定波數k=1.05的地形演變情況，將分別采用FTCS格式和Euler-WENO格式來模擬對應的二維水底地形不穩(wěn)定性演變。

圖7 初始擾動地形Fig.7 Initially perturbed bed

如同上面第1節(jié)一樣，考慮一個有限長度的具有微小傾角的長直渠道，其初始時刻水底為平直的。將坐標系x和y軸取在水底平面上，具體坐標系統(tǒng)如圖1所示。在初始地形上加入以下地形變化的擾動

該擾動在x方向和y方向都是按余弦函數變化的，zb0為擾動的幅值，kx和ky分別為x方向和y方向上的地形擾動的波數：kx=2π/λx，ky=2π/λy，λx和 λy分別為擾動波長。記 Lx和 Ly分別為計算域所包含的渠道的長和寬。計算中取Lx=λx，Ly=0.5λy。由上一節(jié)理論分析知，最大不穩(wěn)定地形波數為k′=1.05，即

該式給出的地形擾動x方向波長為λx=2πLy/k′，即最大不穩(wěn)定擾動的波長是依賴于渠道的寬度Ly。文中計算取Ly=220 m，所以可得擾動波長

因此計算時近似取λx=1 320 m，對應的計算域的長度Lx=λx=1 320 m。其他參數取值為：靜水水深h*=5.0 m，擾動幅值 zb0=1.0 m，輸沙率公式中系數 v=0.1，b=6，γ′=1，uc=0，糙率 nr=0.022 856，坡度 i0=6.11×10-5。具有這些參數的初始地形如圖7所示。

地形邊界采取周期性邊界條件 qy|y=0=0，qy|y=Ly=0，qx|x=0=qx|x=Lx，zb|x=0=zb|x=Lx。

圖8 FTCS格式地形演變三維圖Fig.8 Bed evolution of FTCS scheme in 3D

圖9 Euler-WENO格式地形演變三維圖Fig.9 Bed evolution of Euler-WENO scheme in 3D

圖8給出了FTCS格式在不同時刻t=2 h和5 h的計算結果。由圖8可知，FTCS格式計算進行2 h，左右邊界處的地形開始出現輕微振蕩；隨著計算時間的繼續(xù)增加，邊界處地形的振蕩幅值增加，并且振蕩范圍向域內擴展；整個計算過程只能維持5 h左右，這時地形振蕩明顯，幅值較大，造成計算發(fā)散。圖9中給出了對應的Euler-WENO方法的計算結果，與圖8中FTCS方法的結果不同，地形整體光滑，未出現類似圖8中的數值振蕩。為了進一步顯示該方法的計算穩(wěn)定性，圖10給出了第14 h和40 h地形結果，仍未出現振蕩，表明該方法能夠適應長時間地形演變的模擬計算。由圖10可見，地形在x方向上順水流向前移動，波動幅值隨計算時間逐漸增長，地形波動不再呈初始時刻的余弦曲線狀，而是向鋸齒形剖面形狀演化。從第14 h的三維地形圖可以看出，峰值兩側地形逐漸不對稱，地形峰值上游呈緩坡狀，下游成陡坡狀。

圖11給出了8個時刻地形的等高線圖。圖11中結果顯示，3 h過后，地形開始發(fā)生比較明顯的變化。地形向前移動的同時，沿y方向上也開始有輕微變化，等高線輪廓不再保持左右對稱，隨著計算時間增加，不對稱性更加明顯；由圖11中水流速度矢量可見，流場也隨著地形移動出現輕微的擺動。圖12給出了靠近邊界點（x，y）=（660 m，210 m）處的地形和流速的隨時間的變化歷程。由圖12可知，地形波動的最大幅值由最初的1.0 m在30 h增大到2.0 m，之后保持為常數，地形演變在30 h左右達到一種穩(wěn)定狀態(tài)。圖12中流速在30 h以后呈周期性變化，波動周期為14.5 h。

圖10 Euler-WENO格式地形演變三維圖Fig.10 Bed evolution of Euler-WENO scheme in 3D

圖11 地形演變等高線圖Fig.11 Contour ofbed evolution

綜合以上結果可知，采用Euler-WENO格式的計算地形演變，在精度和穩(wěn)定性方面都優(yōu)于傳統(tǒng)的FTCS方法，且能夠用于長時間的地形演變研究。

圖12 點（660 m，210 m）處u，v和zb的時間歷程Fig.12 Time series of u ，v and zbat（660 m，210 m）

5 結論

文章采用交替方向隱式格式（ADI）求解水動力方程，結合采用具有五階精度的Euler-WENO格式求解的水底變形方程，建立了水流與水底地形變化耦合計算模型。通過對比FTCS格式和Euler-WENO格式求解一維和二維地形演變的計算結果，得知Euler-WENO格式在精度和穩(wěn)定性上都明顯優(yōu)于傳統(tǒng)FTCS格式，能夠用于模擬長時間地形演變過程。

同時應用該模型研究了河流水底地形的不穩(wěn)定性，模擬結果表明交替沙壩地形在單向水流作用下，地形變得不穩(wěn)定，地形變化是三維的，沿縱向和橫向都呈周期性變化，并且方程的非線性特征導致波動峰值是前后不對稱的，呈鋸齒形變形。波動在一定時刻（本文計算為30 h）達到一個穩(wěn)定峰值；對應的水流變化也呈現出周期性波動。水流和地形的變化是相互耦合的。

［1］Hakeem K J，Julio A Z.Controlling spatial oscillations in bed level update schemes［J］.Coastal Engineering，2002，46：109-126.

［2］Justin H，Jesper D，Nick D，et al.Numerical approaches for 1D morphodynamic modelling［J］.Coastal Engineering，2005，52：691-707.

［3］LIU X D，Osher S，CHAN T.Weighted essentially non-oscillatory schemes［J］.Comput.Phys.，1994，115：200-212.

［4］Harten A.High resolution schemes for hyperbolic conservation laws［J］.Comput.Phys.，1983，49：357-393.

［5］趙海洋，劉偉，萬國新.加權基本無振蕩格式研究進展［J］.力學季刊，2005，26（1）：87-95.

ZHAO H Y，LIU W，WAN G X.Research on weighted essentially non-oscillatory schemes［J］.Chinese Quarterly of Mechanics，2005，26（1）：87-95.

［6］趙士清，張鏡潮.連云港潮流的數值模擬［J］.海洋學報，1981，3（3）：500-515.

ZHAO S Q，ZHANG J C.Numerical simulation of current in Lianyungang［J］.Acta Oceanologica Sinica，1981，3（3）：500-515.

［7］Van R L C.Handbook of Sediment Transport by Currents and Waves［R］.Netherlands：Delft Hydraulics，1989.

［8］JIANG G S，WU C C.A high-order WENO finite difference scheme for the equations of ideal magnetohydrodynamics［J］.Comput.Phys.，1999，150：561-594.

［9］WEN L，James T K，ZHI Y S.A numerical scheme for morphological bed level calculations［J］.Coastal Engineering，2008，55：167-180.

［10］Colombini M，Seminara G，Tubino M.Finite-amplitude alternate bars［J］.Fluid Mech，1987，181：213-232.

［11］Schielen R，Doelman A，Swart H E D.On the nonlinear dynamics of free bars in straight channels［J］.Fluid Mech，1993，252：325-356.

Numerical model of morphology instability under flow

WANG Ming-ming，ZOU Zhi-li
（State Key Laboratory of Coastal and Offshore Engineering，Dalian University of Technology，Dalian116024，China）

The interaction between flow and topography may result in different forms of complex morphologies in the river，and alternate bars are one kind of these morphologies.Researches show that the corresponding formation mechanism is the instability of the bottom topography.In this paper,numerical simulation method and results of alternate bars caused by the bed instability were presented.The ADI scheme was implemented to solve the flow equation，and the equation for bed level change was solved by Euler-WENO scheme，which was more stable than the traditional schemes.This coupled computation model of flow and bottom morphology was used to study the evolution of the straight channel with alternate bars.

bed evolution；nonlinear；WENO；instability

P 731.23；O 242.1

A

1005-8443（2011）06-0389-10

2011-05-09；

2011-06-08

國家自然科學基金項目（51079024）；遼寧省教育廳高等學校自然科學基金

王明明（1986－），男，河南省禹州人，碩士研究生，主事海岸及近海工程的研究。