河口海岸水沙動力環(huán)境多尺度時空耦合數(shù)值模擬研究進展

王芳宇，陳汨梨,，張啟博，謝守鵬，張繼生*

(1.河海大學 港口海岸與近海工程學院，南京 210024；2.中國港灣工程有限責任公司，北京 100027)

河口海岸是陸海相互作用的集中地帶，各種物理過程耦合多變，演變機制復(fù)雜，生態(tài)環(huán)境敏感脆弱。同時，該地帶經(jīng)濟較為發(fā)達、人口集居，日益增多的人類活動會加劇河口海岸地區(qū)環(huán)境壓力。因此，河口海岸的研究、開發(fā)和保護是當前世界沿海國家和科學家關(guān)注的熱點問題。

學者們通常采用物理模型試驗、數(shù)值模擬和現(xiàn)場觀測等手段開展河口海岸水沙動力環(huán)境的多尺度時空耦合作用研究。隨著計算機科學技術(shù)的快速發(fā)展，數(shù)值模擬憑借成本低、速度快、預(yù)測性強、結(jié)果可視化等優(yōu)勢成為當前河口海岸水沙動力環(huán)境研究的重要手段，并在工程實踐中得到廣泛應(yīng)用。因此，本文總結(jié)了波浪、潮流、泥沙同一物理場的多尺度時空耦合及不同物理場間的單一時空耦合模擬研究進展，為今后河口海岸水沙動力環(huán)境數(shù)值模擬提供研究思路及方向。

1 同一物理場多尺度空間耦合模擬研究

隨著河口海岸資源的持續(xù)開發(fā)與環(huán)境保護地不斷推進，實現(xiàn)小規(guī)模沿海水沙動力的準確預(yù)測愈發(fā)重要。當使用大范圍水沙動力方法模擬大陸坡及小規(guī)模近海動力環(huán)境時，由于大小范圍網(wǎng)格間的空間尺度差異，在描述垂向網(wǎng)格變化時較為困難[1]。近年來，填海工程造成的區(qū)域泥沙運動特征的改變以及對港口、航道和海岸生態(tài)系統(tǒng)的影響也愈來愈受到廣泛的關(guān)注，并且泥沙淤積帶來的問題尤其復(fù)雜。雖然目前同一物理場大范圍模擬的方法相對成熟，能準確刻畫大范圍水沙動力特性，但在局部水沙特性的刻畫上仍不夠精細。因此，如何考慮同一物理場不同空間尺度的相互作用，已成為準確模擬河口海岸水沙動力特性的關(guān)注重點。

1.1 單向嵌套

近年來，學者們發(fā)現(xiàn)雖計算流體力學方法能求解完整的Navier-Stokes方程，捕捉不同尺度的水流流動現(xiàn)象。但是小尺度模擬對計算資源需求過大，仍不適用于工程實踐。另外，由于空間小尺度現(xiàn)象和大尺度現(xiàn)象間的相互影響不可忽視，迫切需要同步實現(xiàn)大范圍水動力場的預(yù)測及局部動力場的精確模擬。

基于此，學者開始著眼于同一物理場的多尺度空間嵌套模擬，以大范圍模型的模擬結(jié)果作為邊界條件驅(qū)動小范圍模型。吳文挺等[2]通過嵌套模擬的方法搭建了高分辨率的模型網(wǎng)格，將臨近海域的海岸動力過程降尺度傳遞到小尺度河口地區(qū)。隨后M?hring等[3]發(fā)展了流體動力平面定位系統(tǒng)的多尺度建模方式，基于有限元軟件Ansys，使用Flotran求解器實現(xiàn)模型嵌套，證明了嵌套模型相較于簡化模型的優(yōu)越性。Nash和Hartnett[4]利用MSN_Flood平臺開發(fā)了多尺度嵌套模型，該模型在精度和計算量間取得良好平衡且在復(fù)雜嵌套邊界上有獨特的處理方法(如圖1)。圖2展示了Nash等采用四級聯(lián)嵌套網(wǎng)格(90 m、30 m、6 m和2 m四個不同空間尺度)對較粗網(wǎng)格PG數(shù)據(jù)進行空間內(nèi)插，以填充下一級較細子網(wǎng)格CG邊界數(shù)組(P為父網(wǎng)格，C為子網(wǎng)格)。嵌入PG90中的第一層子網(wǎng)格CG30以3:1(△x_P/△x_C)的嵌套比在重點區(qū)域縮小，在5:1的嵌套比下CG30為CG06提供東部邊界條件，最高分辨率的CG02則以3:1嵌套比嵌入CG06中[5]。由于小尺度水流常存在復(fù)雜邊界及源、匯項等問題，且模型精度受差分網(wǎng)格幾何形態(tài)限制[6]，因此在實現(xiàn)多尺度空間水流過程的高效模擬方面仍有欠缺[7]。

圖1 嵌套過程示意圖

泥沙數(shù)學模型作為研究和解決近海水域泥沙沖淤問題的有效工具已得到較為普遍的運用[8]。一維模型計算省時但無法給出各物理量的平面范圍分布，雖二維模型能較準確模擬但計算量較大經(jīng)濟性欠缺[9]。兼顧一、二維模型的優(yōu)缺點，Zhang[10]應(yīng)用一、二維嵌套技術(shù)，將河海連接處視為節(jié)點，由此成功地模擬了黃河口的演變。李玉婷等[11]考慮河口懸沙運動沖淤的變化，采用大小尺度雙重嵌套的方法建立了河口海岸泥沙二維模型，高峰等[12]在分析鴨綠江河口實測水文泥沙資料的基礎(chǔ)上，采用MIKE21 FM模塊建立了基于多重嵌套網(wǎng)格模式的泥沙場數(shù)學模型，驗證表明該嵌套模型能較好地反映泥沙的運動特征，可以為今后開展進一步的相關(guān)工程規(guī)劃和設(shè)計提供有效工具。張修忠等[13]也建立了一種河口海岸一、二維嵌套的泥沙數(shù)學模型，通過交界面實現(xiàn)數(shù)據(jù)在每一迭代步進行耦合計算，能用較少的機時復(fù)演、預(yù)測重點部分的細部變形情況，驗證結(jié)果表明在河口海岸處建立的一、二維嵌套泥沙數(shù)學模型是一種解決某些實際工程問題的可靠和高效的工具。

1.2 雙向耦合

跨空間尺度完成小范圍與大范圍網(wǎng)格之間的數(shù)據(jù)交互，準確模擬多尺度空間耦合現(xiàn)象，對研究各物理場的本質(zhì)至關(guān)重要。因此，為更精細刻畫小尺度復(fù)雜水流過程，實現(xiàn)不同空間尺度下各流場信息的雙向互饋，曾季才等[14]建立多尺度空間數(shù)值離散的三維水流運動耦合模型，用以獲取全區(qū)域大尺度解及多子域小尺度高精度解。該模型采用了三維Darcy水頭插值算法[15]及尺度分離策略，在獲取區(qū)域大尺度粗網(wǎng)格解的同時精細刻畫局部小尺度水流運動過程。將不同尺度模型的空間信息關(guān)聯(lián)傳遞，以實現(xiàn)多尺度空間耦合模擬。相比傳統(tǒng)局部嵌套模型，該模型能更精確刻畫局部尺度的復(fù)雜水流運動過程。

在大規(guī)模沿海水流運動模擬中，Agnieszka等[16]首次提出POM-MSN_Flood建模耦合系統(tǒng)，并將其用于研究復(fù)雜潮汐引起的河口海岸水動力。驗證結(jié)果表明該模型能有效解決與水流特征尺度相對應(yīng)的流體運動，包括東北大西洋的大規(guī)模流動和沿海水域精細環(huán)流。Zhang等[17]在建立近岸波流耦合的非結(jié)構(gòu)化四叉樹網(wǎng)格數(shù)值模型時，以樹狀數(shù)據(jù)結(jié)構(gòu)實現(xiàn)網(wǎng)格連通性和數(shù)據(jù)傳遞性。Wu等[18]為準確模擬小范圍沿海水流的水動力特性，將局部計算流體力學模型[19]與大范圍非結(jié)構(gòu)化網(wǎng)格海洋模型耦合，前者用于捕獲小尺度流，后者用于預(yù)測大尺度流。在局部模型中將小尺度流場劃分為相互重疊的雙向耦合子域，基于Schwartz交替迭代[20]和研發(fā)的質(zhì)量守恒算法[21]，實現(xiàn)子域間解的銜接過渡。如圖3所示，大范圍非結(jié)構(gòu)化網(wǎng)格海洋模型采用內(nèi)外模式求解離散化流動，將局部模型與大范圍模型的三維內(nèi)模式相耦合，在網(wǎng)格交界面上交換速度分量，完成數(shù)據(jù)交互。

3-a 模型示意 3-b 網(wǎng)格示意

隨后Tang等[22]基于該耦合模型成功模擬了污染流從直徑為cm級的出口排放到寬度為km級的河口中并流向海域，并將模擬結(jié)果與僅采用局部模型方式計算的結(jié)果比較，證明該耦合方法的優(yōu)越性。與此同時，Tang等[23]還提出一種預(yù)測沿海海流的多尺度方法，將基于HM(Hybrid method)和DDM(Domain decomposition method)的計算流體力學和地球物理流體力學模型與Chimera交疊網(wǎng)格[24]相結(jié)合，實現(xiàn)了三維不可壓縮計算流體力學模型與大范圍非結(jié)構(gòu)化網(wǎng)格海洋模型的雙向耦合，并提供了數(shù)值示例證明其可行性，指出多尺度空間耦合模擬有很大的發(fā)展空間和應(yīng)用前景。

河口海岸泥沙輸運是陸海物質(zhì)循環(huán)的重要過程，對物質(zhì)入海通量和全球物質(zhì)循環(huán)具有重要研究意義。劉紀根等[25]認為當前泥沙輸移過程隨流域尺度變化的研究缺乏應(yīng)有的廣度和深度，總結(jié)了泥沙沖刷輸運的空間尺度相關(guān)問題，探討了不同尺度下的泥沙輸運現(xiàn)象與規(guī)律。朱志夏等[26]針對三維泥沙數(shù)值模擬難題，采用FVCOM模型及非結(jié)構(gòu)化網(wǎng)格二維、三維嵌套方法，建立了三維泥沙數(shù)學模型。Bai等[27]對河口海岸地區(qū)泥沙輸移現(xiàn)象展開研究，采用有限元法FEM和有限差分法FDM相結(jié)合的方法提出海河河口模型，并在三維指標系統(tǒng)中引入多層模型，最后將模型應(yīng)用于渾濁的華北海河河口。通過對比數(shù)值結(jié)果與實測資料，證明該多空間尺度泥沙模型與美國CH3D-SED模型具有可比性。

本節(jié)回顧了學者針對不同物理場展開的多尺度空間耦合模擬，由于小尺度和大尺度水、沙相互作用，當進行多尺度空間耦合模擬時，因大小網(wǎng)格的空間尺度差異性，通常會出現(xiàn)大范圍水動力的一個網(wǎng)格涵蓋了十至百個精細網(wǎng)格的現(xiàn)象。隨著數(shù)值模擬手段的不斷深入，各物理量的空間尺度跨度是逐漸增大的，例如，泥沙的大規(guī)模輸運、海工建筑物的局部沖淤、泥沙顆粒運動規(guī)律等物理現(xiàn)象的模擬，均從宏觀到微觀不斷深入，并且在該過程中空間尺度量級差是不可避免的。因此，如何在數(shù)值模擬中有效解決空間尺度量級差帶來的問題，已成為真實模擬各物理場的研究重點及難點。

2 同一物理場多尺度時間耦合模擬研究

將不同時間尺度的物理場進行耦合模擬，預(yù)測未來時間內(nèi)該物理場的真實變化，可在某種程度上滿足當前對長時間變化物理現(xiàn)象的現(xiàn)實需求。雖前人曾采用統(tǒng)計學模型的方式對河口海岸演變進行預(yù)測并得到?jīng)_淤平衡點[27-30]，但研究性質(zhì)僅限于規(guī)律，無法從機理上對長期演變進行很好地解釋，并且也無法對海岸變化的具體位置進行預(yù)測。因此多時間尺度模擬預(yù)測未來長時間變化是有必要的，如泥沙輸移的多尺度時間數(shù)值模擬被普遍認為是研究河口海岸地貌演變的重要工具?；诖?，學者們逐漸關(guān)注各物理場的多尺度時間耦合模擬，期盼能以當前時間尺度預(yù)測未來時間內(nèi)物理場變化，對厘清河口海岸水沙動力環(huán)境及相互作用機制具有重要意義[31]。

2.1 加速因子

當模擬中長期河口海岸動力地貌演變時，通常對地貌演變的時間進行加速以快速模擬動力地貌的中長期變化?；诖擞^念，Lesser[32]提出地貌加速因子法，這是一種用來解決地貌發(fā)展過程中出現(xiàn)時間尺度差異的有效方法。這一技術(shù)與Latteux[33]提出的“延長潮汐”概念相類似，旨在用短時間泥沙運動模擬實現(xiàn)長時間地貌演變模擬。

地貌加速因子法通過引入一個地貌加速因子(morphological acceleration factor)，有效擴展地貌時間步長以完成地貌變化的快速計算。通過乘上一個常量因子增加地形變化率，經(jīng)過1個潮周期的計算，模型實際上已對n(n=fMOR)個周期的地貌變化進行了模擬，如式(1)

△tmorphology=fMOR△thyfordynamic

(1)

式中：fMOR為地貌加速因子；△t為時間步長。

由于地貌演變作用的復(fù)雜性與時空尺度多樣性，當前對河口海岸動力地貌演變的理解與預(yù)測仍非常有限。Ganju等[34]首先利用ROMS模型展開河口地區(qū)多尺度時間嵌套的地貌模擬，概化了入海輸沙量并引入地貌加速因子。Ganju等[35-36]也對美國Suisun灣1867年—1887年地貌沖淤開展了模擬，與實際地形變化對比預(yù)測了Suisun灣2010年—2030年際地貌演變[37]，并指出模擬年代際尺度上的地形沖淤變化難度較大。

荷蘭學者們[38]利用Delft3D模型分別模擬研究了Scheldt河口西部的長時間地貌演變，實現(xiàn)多時間尺度耦合模擬研究。董程等[39]從黃河潮灘地貌演化角度出發(fā)，認為采用地貌加速因子后一個潮周期的地貌變化等價于10個潮周期的地貌，表明以多尺度時間方式模擬地貌變化是非常必要的。蔣超[40]系統(tǒng)分析了黃河入海水沙通量變化特征，揭示了黃河口不同時間尺度地貌演變的時空變化特征，為不同類型動力地貌的多尺度時間耦合模擬建立重要紐帶。Luan等[41-42]闡述了主導(dǎo)長江口年代際地貌沖淤演變的影響因素，基于Delft3D模型建立長江口年代際地貌演變數(shù)值模型，預(yù)測了未來年代際演變的動態(tài)過程。在珠江口方面，Xing和Wu[43]基于動力地貌數(shù)值模型模擬了珠江口長時間沖淤變化過程，指明河口地形的長期演變與復(fù)雜的河網(wǎng)格局密切相關(guān)。吳超羽等[44]綜合考慮影響珠江口動力地貌演變的驅(qū)動因子，建立長時間動力地貌演變數(shù)值耦合模型，模擬了珠江口幾千年地貌、河網(wǎng)演變過程?？紤]工程實際情況，竇希萍等[45]認為在工程初期地形沖淤變化往往較劇烈，此時不能選擇大的加速因子。隨著水沙動力與河床不斷調(diào)整，地形沖淤變化趨緩，可選擇較大加速因子以減小地貌長周期模擬計算時間。當模擬工程影響下的地貌演變時，采用隨時間變化的加速因子更能提高模擬精度和計算效率。但在中長期地貌演變模擬過程中使用加速因子，模型需要設(shè)置對應(yīng)尺度的上游和外海邊界條件。如果仍采用實測邊界條件，那么水動力過程在時間尺度上會被人為“拉長”，一個潮周期對應(yīng)的地貌時間可能為數(shù)天甚至數(shù)月，顯然不符合實際物理意義，因此在選擇相應(yīng)邊界條件需要考量得當。

2.2 長期運算

由于計算資源有限，以往的數(shù)值模擬多為短時間模擬，并且大多數(shù)基于物理機制的海岸演化模型為短時間尺度下“真實性演變過程”的一個描述，難以刻畫中長時間尺度的海岸演化過程和機制[46-48]。隨著計算技術(shù)進步，長期運算模擬手段逐漸完善。如前所述，為預(yù)測未來長周期河口海岸水沙動力變化，常用手段是使用加速因子加快時間進程，從而縮短模擬時長。但是，Roelvink等[49]指出使用加速因子的前提是水動力計算時間步長內(nèi)泥沙凈輸運量和地形沖淤變化很小，乘以加速因子后的地形沖淤變化不會對未來年際內(nèi)地貌格局產(chǎn)生較大影響。

基于上述觀點，學者們嘗試采用長期運算手段模擬河口海岸水沙動力的長期演變，以期預(yù)測未來更長時間尺度內(nèi)各物理量的變化過程。Roelvink等[50]基于二維(水平)模型模擬復(fù)雜海岸演變，認為該模型可數(shù)值再現(xiàn)多年海岸形態(tài)，但計算成本高且需要豐富專業(yè)知識。為克服現(xiàn)有海岸線模型局限性，他們開發(fā)了一種全新的海岸線模擬模型(Shorelines)，旨在預(yù)測數(shù)年到數(shù)世紀的海岸線演變。研究結(jié)果證明該模型不僅能模擬漸進的海岸變化，并且能預(yù)測未來時間尺度內(nèi)海岸的根本性轉(zhuǎn)變。該模型適用于局部到區(qū)域甚至全球領(lǐng)域，以合理的計算成本在完整生命周期中穩(wěn)健跟蹤海岸特征，捕捉海岸演化的復(fù)雜性，具有顯著的優(yōu)越性。

發(fā)展中長時間尺度的海岸演化預(yù)測方法與構(gòu)建機理模型是目前國際海岸演化研究領(lǐng)域關(guān)注的焦點和極具挑戰(zhàn)性的主題。Hanson等[47]對目前已有的基于短時間尺度發(fā)展的20個不同類型年代際尺度海岸演化模型進行對比分析和討論，結(jié)果表明這類模型大多不能進行長時間尺度海岸演化過程的準確模擬?；诖耍顕鴦俚萚51]提出將一線模型(one-line model)與三維海洋動力學數(shù)值模型耦合，經(jīng)過參數(shù)化后形成真正意義上的機理模型，能在長時間尺度上對河口海岸演化的驅(qū)動機制和影響機理進行深入地研究。錢繼春等[52]基于潮流泥沙數(shù)學模型和灘涂動態(tài)平衡的經(jīng)驗關(guān)系，建立了長歷時灘涂動態(tài)演變的經(jīng)驗動力數(shù)值模型，結(jié)果表明數(shù)值結(jié)果與實測基本吻合，能真實反映灘涂的動態(tài)演變過程。王寧舸[53]隨后對江蘇中部流域的河口潮灘開展1 a現(xiàn)場觀測，建立平面二維泥沙輸運及潮灘演變概化數(shù)學模型，模擬2.5 a內(nèi)江蘇中部潮灘變化，預(yù)測了未來年際內(nèi)潮灘演變規(guī)律。程梁秋等[54]基于多年海圖數(shù)據(jù)建立了相對應(yīng)的數(shù)字高程模型，對伶仃洋西槽地貌展開長達15 a的數(shù)值模擬，定量分析其長期演變特征。龔政等[55]以淤長型潮灘為研究對象，建立了潮灘地貌的演變概化數(shù)學模型，預(yù)測了未來16 a內(nèi)潮灘剖面的發(fā)育和演變過程，并與經(jīng)典潮灘一維數(shù)值結(jié)果對比，證明建立的長期運算模型具有較高可靠性。鄭釗[56]考慮黃河沖淤演變建立了數(shù)學模型以此預(yù)測未來50 a黃河的沖淤演變，并且驗證了預(yù)測成果的可靠性，該預(yù)測成果對未來黃河河口治理具有重要的實際價值。

石小翠[57]從預(yù)防海洋災(zāi)害角度出發(fā)，運用SWAN模型嵌套模擬計算了中國東部海區(qū)22 a波要素，以此構(gòu)成中國東部海區(qū)的波浪數(shù)據(jù)庫。Liang等[58]也采用SWAN模型以當前時間尺度模擬了渤海、黃海和東海21 a的波浪要素。Langodan等[59]采用ERA-Interim全球再分析初始化WRF模式生成高分辨率區(qū)域驅(qū)動WAVEWATCH III模式，利用30 a波浪要素追算了紅海波浪變化趨勢。徐佳麗等[60]采用TOMAWAC模型模擬中國近海40 a的波浪要素，建立相對完整的數(shù)據(jù)庫以全面描述中國海區(qū)波浪情況，以此為依據(jù)探究河口海岸地區(qū)動力演變過程。

3 多物理場單一尺度時空耦合模擬研究

河口海岸水沙動力環(huán)境的演變是個復(fù)雜的多時空尺度過程，不同時空尺度的演化過程及動力機制存在顯著差異，相關(guān)研究的理論基礎(chǔ)及技術(shù)方法也各不相同[61]。自20世紀以來，學者致力于數(shù)值模擬各物理現(xiàn)象，但已有研究主要聚焦于各種相對孤立的物理效應(yīng)。在通常情況下，聚焦主要物理場時會采用概化方法模擬其他物理場[62]。對比學者們的模擬手段發(fā)現(xiàn)，聚焦主要物理場并概化次要物理場是常見且實用的研究手段。然而，現(xiàn)實中的物理現(xiàn)象并不是孤立發(fā)生的，只重點關(guān)注單一物理場會導(dǎo)致模擬結(jié)果與真實情況有出入。

隨著計算能力不斷提升，為更全面真實模擬河口海岸水沙動力環(huán)境，學者們開始在數(shù)學模型中考慮多物理場的相互耦合作用。

3.1 單向嵌套

在河口海岸地區(qū)波浪和潮流是主要驅(qū)動力，其過程復(fù)雜且呈非線性[63]。近岸區(qū)波浪在傳播時會受到水流的影響，波浪自身也會產(chǎn)生近岸水流驅(qū)動力，兩者存在耦合作用，并且波浪與潮流在時間、空間尺度上差異顯著，因此研究河口海岸地區(qū)的波流相互作用具有一定意義[64]。當前波流相互作用的研究主要采用兩種概化方法：一是把周期變化的潮流概化為具有某一特征流疊加到波浪運動方程中，用以模擬短時間內(nèi)波流運動的變化；另一種方法是把波浪作用過程概化為潮周期中具有平均意義的波浪要素疊加到潮流運動方程中，用以描述波對流的影響。

波流模型經(jīng)歷了從一維到二維再到三維、由單向嵌套再到雙向耦合等由簡單到復(fù)雜的發(fā)展歷程[65]，Longue-Higgins和Stewart[66]首先基于輻射應(yīng)力概念闡明了波流相互作用機理以及波流的能量轉(zhuǎn)換關(guān)系；Janssen[67]隨后提出波的阻力關(guān)系，進一步豐富波流相互作用理論。在過去幾十年，特別是近20 a，學者們從不同角度對波流相互作用展開細致模擬研究。Teles等[68]基于RANS方程將水流模塊嵌套至波浪模塊，采用任意拉格朗日-歐拉(ALE)方法模擬時變自由表面動力學，在波浪時空尺度上同時求解波流運動數(shù)值解，分析自由表面流動中波流的聯(lián)合效應(yīng)。許忠厚等[69]基于Delft3D建立了三維波流耦合模型，模擬了波浪作用下出水離岸堤附近的波浪與波生流間的關(guān)系，Zhang等[17]通過求解二維深度平均非線性淺水方程，在沿海水域的水動力模塊加入波浪誘導(dǎo)的沖浪區(qū)近水平流，波流相互作用模擬結(jié)果與實測數(shù)據(jù)吻合良好。Sanchez等[70]采用二維深度平均方程求解波流作用下潮流從十幾公里寬的Colombia河口流向幾十萬平方公里海域的地形變化，實現(xiàn)了基于波流相互作用的地貌演變模擬。

在水沙動力環(huán)境下，泥沙常隨水流一同運動。由于水流是周期性的往復(fù)運動流，其作用下的泥沙運動也是往復(fù)運動，運動機理要比河流中單向恒定作用下的泥沙運動機理復(fù)雜得多[71]。因此，研究近岸區(qū)泥沙輸移及岸灘演變具有十分重要的意義[72]。受限于計算流體力學理論不成熟，在過去的幾十年中水沙模型多為非耦合模型。近20 a來得益于計算效率的提高，吳偉明等[73]通過處理連接交界面讓一維和二維水沙數(shù)學模型滿足水流、泥沙運動連接的邊界條件，從而構(gòu)建了水沙運動一維、二維嵌套模型。張華慶等[74]采用正交網(wǎng)格建立了珠江河口平面二維潮流泥沙數(shù)學模型，將珠江三角洲一維河網(wǎng)模型嵌套至珠江河口二維模型，在不同空間尺度上構(gòu)建了珠江河口水沙數(shù)值模擬系統(tǒng)的核心計算模塊。張世奇[75]在河口采用一維模式并在海岸運用二維模式，將入海口作為河海連接點，既作為一維模擬的下邊界，又作為二維模擬的上邊界。根據(jù)水流連續(xù)原理交替使用一維、二維尺度手段求解水位流量，采用擬合的方法使其趨于一致，由此建立一維、二維連接的沖淤嵌套模型，結(jié)果表明該模型能較準確模擬黃河口沖淤地形和輸沙過程。

3.2 雙向耦合

隨著科學進步以及嵌套技術(shù)的成熟，信息的雙向傳遞已經(jīng)成為河口海岸水沙動力環(huán)境數(shù)值模擬研究的關(guān)鍵，多物理場的數(shù)值模擬方式也從單向嵌套逐步轉(zhuǎn)到雙向耦合。Dietrich等[76]采用區(qū)域分解并行仿真技術(shù)，在同一非結(jié)構(gòu)化網(wǎng)格上交換網(wǎng)格頂點的解以實現(xiàn)SWAN和ADCIRC耦合，有效解決了波流相互作用的模擬難題。該模型具有高度可擴展性和通信性，在較大空間梯度上提高了局部分辨率，且能準確地模擬波流相互作用。王平等[77]基于非結(jié)構(gòu)化網(wǎng)格同時模擬潮流和波浪，并通過參數(shù)的同步傳遞實現(xiàn)波流的耦合計算，結(jié)果表明模型對模擬近岸波流的耦合作用有著很好的精度和適用性。王文鼎等[78]采用大尺度和小尺度兩種尺度的網(wǎng)格建立波流耦合數(shù)值模型，并將研究結(jié)論應(yīng)用于部分工程實踐，驗證了模型的符合性、有效性和適用性。Roland等[79]基于三維水動力模型SELFE和波譜模型WWM-II，提出一種基于非結(jié)構(gòu)網(wǎng)格的波流全耦合模型。SELFE和WWM-II模型間的耦合是在源代碼級別進行的，在共享相同子域的基礎(chǔ)上并行MPI和域分解，從根本上保證并行效率和避免插值。Roland等[79]隨后將模擬結(jié)果和現(xiàn)場資料比較，結(jié)果表明即使在復(fù)雜的實際環(huán)境下，耦合模型也能捕捉到不同尺度下波流相互作用過程。

在水流問題的研究上，王智勇等[80]提出耦合邊界水位預(yù)測校正法，由此建立一維和二維耦合水動力模型。在任一時間步計算時先在耦合連接處預(yù)測一維和二維水位邊界條件，再根據(jù)耦合邊界處的流量計算結(jié)果及當?shù)厮?shù)對水位邊界條件反復(fù)校正，直至耦合連接條件滿足計算容差。顧杰等[81]基于MIKE-FLOOD模塊建立秦皇島入海河流與近岸海域間的一維、二維耦合河流-海岸水動力模型，系統(tǒng)研究了秦皇島海域入海河流的水動力特性。來志剛[82]在珠江口實現(xiàn)了一維河網(wǎng)、三維河口的水動力耦合模擬，成功將一維、二維水動力及泥沙連接模型基本條件和計算方法擴展至一維、三維水動力模型的耦合模擬。Zhou等[83]基于質(zhì)量守恒準則，耦合了一維河網(wǎng)模型與三維河口模型，為珠江三角洲的水環(huán)境監(jiān)控提供可靠依據(jù)。同時，Lai等[84]著眼于長江中游復(fù)雜河湖大水系，提出一種新的耦合水動力模型CHAM，為河口海岸地區(qū)的復(fù)雜流態(tài)提供了有效模擬工具。

準確描述河口海岸的水流運動是合理描述泥沙運動的基礎(chǔ)，目前相關(guān)學者在泥沙數(shù)值模擬方面已開展了相關(guān)研究并取得一定成果[85]，張修忠等[86]指出對于解決河口海岸這類具有多空間尺度的水沙問題，采用單一模型不僅困難、費時，且不易掌握不同區(qū)域矛盾的主要特征?；诖藦埿拗业萚86]利用有限元的非結(jié)構(gòu)化網(wǎng)格特點和子結(jié)構(gòu)疊加的概念，建立了一種虛擬過度區(qū)間的一、二維耦合算法。蘇東升等[87]基于CFD-DEM 方法建立水沙運動耦合模型，張華慶等[88]考慮水沙相互關(guān)系建立了一維、二維耦合的水沙模型，將河口地區(qū)視為重點區(qū)域進行計算。張鵬飛和陸建剛[89]研究了通江湖泊與外部江河的水沙交換過程，基于二維非穩(wěn)態(tài)水流和泥沙耦合數(shù)學模型，模擬了鄱陽湖與外部江河的水沙交換過程。Li等[90]基于三維水動力模型與泥沙運動模型，耦合模擬了北部灣海岸地區(qū)沙波的遷移過程，模擬結(jié)果與實測基本吻合，表明多場耦合模型可以模擬現(xiàn)實中的沙波遷移規(guī)律。在驗證耦合模型能較好模擬水沙運動后，Shen等[91]考慮河口和沿海地區(qū)水域懸浮沉積物濃度的高變異性特征，建立了泥沙模型與開源TELEMAC系統(tǒng)耦合的平衡模型。葉濤焱[92]為探究杭州灣人類活動影響下的懸沙濃度，建立了三維水沙數(shù)值耦合模型ESed，分析了杭州灣灣口懸沙動力與潮灘變化的互饋機理。

由此可見，在多物理場單一尺度時空耦合數(shù)值模擬中，通常在求解控制方程時考慮多因子的相互影響。本節(jié)總結(jié)了學者在河口海岸地區(qū)波浪、潮流以及泥沙方面的研究成果，說明了考慮不同物理場的相互作用能更真實反應(yīng)河口海岸的水沙動力環(huán)境的變化。但目前數(shù)值模擬多為單一尺度時空耦合，在準確刻畫不同時空尺度下各物理場間的相互作用上仍有所欠缺。

4 模擬難題與建議

在河口海岸水沙動力環(huán)境的數(shù)值模擬中，模擬手段從單獨模擬到單向整合再發(fā)展至雙向耦合。以往的模型大多僅實現(xiàn)信息單向傳遞，即給定初始、邊界條件時A模型先計算，再將A模擬結(jié)果傳至B模型，此時B模型結(jié)果并未反饋給A模型。但是，較為理想的是結(jié)果相互反饋：A將結(jié)果傳遞至B，B同時反饋至A，循環(huán)迭代直到各物理場的計算達到穩(wěn)定。在整個模擬過程中，數(shù)據(jù)在每一時間步都進行交換，據(jù)此實現(xiàn)兩個模型間的實時雙向耦合。

實際河口海岸水沙動力環(huán)境是多物理場相互作用的結(jié)果，早在1993年Orton和Reading[93]就提出河口系統(tǒng)每一個地貌單元并不是獨立的個體，而存在頻繁的物質(zhì)交換和相互作用過程?？紤]不同物理場相互作用時，數(shù)據(jù)相互反饋使耦合模擬結(jié)果更精確。隨著社會發(fā)展及行業(yè)需求，研究重點、熱點也更加聚焦于多物理場的多尺度時空耦合模擬，河口海岸地區(qū)多物理場間的相互作用問題也越來越廣泛受到人們的重視，例如：海冰、波浪、水流的相互作用[94]，波浪、水流、泥沙的相互作用[95]，風、波浪、水流的相互作用[96]，大氣、海浪、海流的相互作用等[97]。與此同時，多物理場相互作用現(xiàn)象越發(fā)復(fù)雜，水沙動力與結(jié)構(gòu)物響應(yīng)、波浪-海床-結(jié)構(gòu)物[98]相互作用等問題也逐漸突出，多尺度計算已經(jīng)發(fā)展成為一個具有自身特征的學科，是構(gòu)成人們探索和理解物質(zhì)世界多尺度問題的基礎(chǔ)[99]。因此，在求解多物理場相互作用時，如何客觀描述控制方程是提高數(shù)值模擬精度、實現(xiàn)多尺度時空耦合的研究重點，也是亟待解決的關(guān)鍵問題。

受限于計算技術(shù)及耦合算法的難度，多物理場的全耦合模擬尚未完全實現(xiàn)。結(jié)合現(xiàn)今技術(shù)手段，采用以下三種方法實現(xiàn)模型間耦合較可行。第一種方法是將兩種模型代碼直接合并，基于內(nèi)存實現(xiàn)模型間的數(shù)據(jù)傳遞，但通常兩個模型是獨立開發(fā)，因此數(shù)據(jù)的交互具有一定難度。第二種方法是利用現(xiàn)有通信工具或協(xié)議如MPI，但該方法需模型開發(fā)人員具備專業(yè)知識。第三種方法是使用模型耦合工具提供自動傳輸數(shù)據(jù)接口，這種方法能極大減少模型開發(fā)期間修改代碼的工作量，并且非專業(yè)技術(shù)人員也能實現(xiàn)模型接口配置，相比下具有更優(yōu)性。雖第三種方法更具優(yōu)越性，但關(guān)于耦合模型接口的搭建，實現(xiàn)數(shù)據(jù)的自動傳輸難度很大。因此如何開發(fā)模型耦合工具及使用好該工具，實現(xiàn)信息交互、耦合的自動化以及模型結(jié)果的相互反饋是學者們非常期待的。

5 結(jié)語

河口海岸是地球四大圈層交匯、能量流和物質(zhì)流的重要聚散地帶，既是流域物質(zhì)的歸宿，又是海洋的開始，生態(tài)環(huán)境極其敏感、脆弱。高精度數(shù)值模擬手段是未來河口海岸水沙動力耦合作用機理及其環(huán)境影響效應(yīng)評估等研究的主要手段和發(fā)展趨勢。本文回顧了河口海岸地區(qū)潮流、波浪及泥沙等多物理場間的耦合作用及同一物理場下多尺度空間、時間耦合作用，探討了多尺度時空耦合過程中周期量級差、尺度量級差、各物理場間的相互耦合、多尺度時空耦合技術(shù)以及多物理場多尺度時空耦合等關(guān)鍵難題，指出現(xiàn)階段實現(xiàn)河口海岸多物理場多尺度時空數(shù)值耦合的模擬仍面臨許多挑戰(zhàn)，急需展開深入系統(tǒng)研究和工程實踐。