王明明，鄒志利

（大連理工大學(xué)海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，大連 116024）

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

王明明，鄒志利

（大連理工大學(xué)海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，大連 116024）

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

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

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

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

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

1 水流和水底地形的耦合數(shù)學(xué)模型

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

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

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

vn+1/2的顯示求解：

圖2 交錯(cuò)網(wǎng)格布置圖Fig.2 Disposition of variables in staggered grid

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

x方向的動(dòng)量方程采用半隱式離散，其離散格式為

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

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

將輸沙率q在x正方向和負(fù)方向分解成與地形傳播速度相關(guān)的兩部分

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

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

2 模型驗(yàn)證

2.1 水動(dòng)力模型驗(yàn)證

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

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

圖3 數(shù)值模擬的流速Fig.3 Simulation of velocity

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

2.2 地形演變模型驗(yàn)證

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

初始地形為

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

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

圖4 FTCS格式計(jì)算結(jié)果與解析解的對(duì)比Fig.4 Comparison between FTCS results and analytical solutions

圖5 Euler格式計(jì)算結(jié)果與解析解的對(duì)比Fig.5 Comparison between Euler-WENO results and analytical solutions

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

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

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

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

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

圖6 ω和k的關(guān)系曲線Fig.6 Relation curves ofωand k

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

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

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

圖7 初始擾動(dòng)地形Fig.7 Initially perturbed bed

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

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

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

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

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

圖8 FTCS格式地形演變?nèi)S圖Fig.8 Bed evolution of FTCS scheme in 3D

圖9 Euler-WENO格式地形演變?nèi)S圖Fig.9 Bed evolution of Euler-WENO scheme in 3D

圖8給出了FTCS格式在不同時(shí)刻t=2 h和5 h的計(jì)算結(jié)果。由圖8可知，F(xiàn)TCS格式計(jì)算進(jìn)行2 h，左右邊界處的地形開(kāi)始出現(xiàn)輕微振蕩；隨著計(jì)算時(shí)間的繼續(xù)增加，邊界處地形的振蕩幅值增加，并且振蕩范圍向域內(nèi)擴(kuò)展；整個(gè)計(jì)算過(guò)程只能維持5 h左右，這時(shí)地形振蕩明顯，幅值較大，造成計(jì)算發(fā)散。圖9中給出了對(duì)應(yīng)的Euler-WENO方法的計(jì)算結(jié)果，與圖8中FTCS方法的結(jié)果不同，地形整體光滑，未出現(xiàn)類似圖8中的數(shù)值振蕩。為了進(jìn)一步顯示該方法的計(jì)算穩(wěn)定性，圖10給出了第14 h和40 h地形結(jié)果，仍未出現(xiàn)振蕩，表明該方法能夠適應(yīng)長(zhǎng)時(shí)間地形演變的模擬計(jì)算。由圖10可見(jiàn)，地形在x方向上順?biāo)飨蚯耙苿?dòng)，波動(dòng)幅值隨計(jì)算時(shí)間逐漸增長(zhǎng)，地形波動(dòng)不再呈初始時(shí)刻的余弦曲線狀，而是向鋸齒形剖面形狀演化。從第14 h的三維地形圖可以看出，峰值兩側(cè)地形逐漸不對(duì)稱，地形峰值上游呈緩坡?tīng)?，下游成陡坡?tīng)睢?/p>

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

圖10 Euler-WENO格式地形演變?nèi)S圖Fig.10 Bed evolution of Euler-WENO scheme in 3D

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

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

圖12 點(diǎn)（660 m，210 m）處u，v和zb的時(shí)間歷程Fig.12 Time series of u ，v and zbat（660 m，210 m）

5 結(jié)論

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

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

［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］趙海洋，劉偉，萬(wàn)國(guó)新.加權(quán)基本無(wú)振蕩格式研究進(jìn)展［J］.力學(xué)季刊，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］趙士清，張鏡潮.連云港潮流的數(shù)值模擬［J］.海洋學(xué)報(bào)，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

國(guó)家自然科學(xué)基金項(xiàng)目（51079024）；遼寧省教育廳高等學(xué)校自然科學(xué)基金

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