基于半隱歐拉－拉格朗日法鹽水入侵數(shù)學(xué)模型

匡翠萍，黃 靜，陳思宇，2，劉曙光

（1.同濟大學(xué) 土木工程學(xué)院，上海 200092；2.北卡羅來納州立大學(xué) 土木、結(jié)構(gòu)和環(huán)境工程系，羅利 27606）

河口地區(qū)鹽水入侵會使河口淡水資源的水質(zhì)惡化，從而對兩岸人民的生活用水及工農(nóng)業(yè)用水造成嚴重的影響.數(shù)學(xué)模型常被用來研究長江口鹽水入侵現(xiàn)象.肖成猷等［1］建立了長江口垂向平均的二維數(shù)學(xué)模型對北支鹽水倒灌過程及其影響范圍和影響時間尺度進行了模擬研究.GU等［2］運用Delft3DFLOW軟件建立了長江口二維潮流鹽度數(shù)學(xué)模型，利用枯季大潮時北支束窄方案對長江口的鹽度進行了計算分析.羅小峰等［3］通過建立長江口平面二維水流鹽度數(shù)學(xué)模型對4種不同水文組合條件下徑流和潮流對鹽水入侵的影響進行了計算還建立了長江口北槽三維水流鹽度數(shù)值模型，研究了鹽度梯度引起的斜壓項對水流的影響.匡翠萍等［4－5］通過建立長江口鹽水入侵三維數(shù)學(xué)模型對長江口鹽度分層現(xiàn)象的特征進行了計算分析.此外，還采用了干濕判別法使ECOMSED模型實現(xiàn)了動邊界的處理，并在此基礎(chǔ)上建立了長江口三維潮流數(shù)學(xué)模型，對大潮時長江口的水動力進行了計算模擬.朱建榮等［6］、馬鋼峰等［7］基于ECOM模型亦分別建立了長江口三維鹽度數(shù)學(xué)模型對長江口鹽水入侵進行了模擬研究.

在平面坐標(biāo)系下，矩形網(wǎng)格的系統(tǒng)數(shù)據(jù)結(jié)構(gòu)組織簡單，適用于各種算法，處理效率高，但不能準(zhǔn)確地刻畫復(fù)雜的幾何形狀.正交曲線網(wǎng)格雖然能貼合計算區(qū)域邊界，但無法刻畫海岸帶離散分布的眾多島嶼、潮溝以及工程區(qū)域.而無結(jié)構(gòu)網(wǎng)格不僅能很好地擬合邊界，而且能控制網(wǎng)格的密度，易進行適應(yīng)性調(diào)整，從而提高了邊界模擬的精度.

近年來，許多專家學(xué)者基于無結(jié)構(gòu)網(wǎng)格建立了合適的數(shù)學(xué)模型對河口的水動力和鹽水入侵進行了計算研究.王志力等［8－10］采用無結(jié)構(gòu)網(wǎng)格的有限體積方法對二維河道進行了數(shù)值模擬，并在此基礎(chǔ)上分別建立了物質(zhì)輸運的高精度無數(shù)值震蕩數(shù)學(xué)模型和物質(zhì)輸運對流擴散方程離散模式，最后都應(yīng)用于甌江河口鹽度輸運的模擬研究.潘存鴻等［11－12］在無結(jié)構(gòu)網(wǎng)格下分別建立了求解二維淺水流動方程的具有二階精度的KFVS（kinetic flux vector splitting）和諧格式和Godunov格式，并分別利用典型算例和錢塘江涌潮對模型進行了計算驗證.陳昞睿等［13］建立了一個無結(jié)構(gòu)三角形網(wǎng)格二維河口海岸水動力數(shù)值模式，并將該模式應(yīng)用于長江口水動力的計算研究.

Casulli和 Walters［14］于2000年提出了一個基于正交非結(jié)構(gòu)網(wǎng)格的水動力學(xué)數(shù)值模型——UnTRIM模型，并得到了應(yīng)用.例如，SHEN等［15］運用UnTRIM模型模擬研究了美國切薩皮克灣（Chesapeake Bay）由颶風(fēng)伊莎貝爾（Isabel）引起的風(fēng)暴潮.CHEN等［16］基于UnTRIM 模型建立了二維潮流數(shù)學(xué)模型，網(wǎng)格模塊是通過多元最小二乘法重構(gòu)的無結(jié)構(gòu)三角網(wǎng)格，潮流模塊基于消除了穩(wěn)定性條件限制的半隱式的歐拉－拉格朗日法，動量方程中流通量對時間的導(dǎo)數(shù)項和對流項寫成全導(dǎo)數(shù)形式，采用拉格朗日法求解，從而避免了Courant數(shù)限制，并將模型應(yīng)用到了曹妃甸圍墾工程的潮流模擬.LIU等［17］基于UnTRIM模型建立了一個時變的三維正壓水流鹽度模型，并將其運用于臺灣北部淡水河河口地區(qū)鹽度的計算分析.

本文在文獻［16］研究的基礎(chǔ)上增加鹽度輸運模塊建立了二維潮流鹽度數(shù)學(xué)模型.鹽度輸運模塊采用有限體積法進行離散，通過調(diào)整界面通量取值以及用多元最小二乘法重構(gòu)的濃度分布值代替單元平均濃度值的方法得到了與連續(xù)性離散方程相協(xié)調(diào)的二階對流擴散離散方程，并且分別通過對流和擴散數(shù)值測試對模型進行了驗證，最后將該模型成功應(yīng)

用于模擬長江口的鹽水入侵計算.

1 數(shù)學(xué)模型

1.1 基本方程

二維潮流鹽度數(shù)學(xué)模型的控制方程為沿水深垂向平均的淺水方程.連續(xù)性方程為

式中：η為水位；t為時間；q為界面流通量，q＝DU，D為總水深，D＝H＋η，H為水深，U為流速，U＝（u，v），u，v分別為x方向和y方向的流速.動量方程為

式中：S為鹽度；Kh為鹽度擴散系數(shù).狀態(tài)方程為

式中：ρ0為淡水密度；λ為常數(shù)，λ＝0.75.

1.2 數(shù)值方法

1.2.1 方程離散

（1）潮流方程的離散.連續(xù)性方程采用半隱有限體積法，動量方程在各單元界面上進行離散，流通量對時間的導(dǎo)數(shù)項和對流項寫成全導(dǎo)數(shù)形式，采用拉格朗日法求解從而避免了Courant數(shù)限制，底部摩擦項采用隱格式離散，而擴散項、科氏力項等均采用顯格式，最終得到的離散方程如下.離散的連續(xù)性方程為

式中：δAi為單元i的面積分別為單元i在n和n＋1時刻的水位；δlfi為界面fi的長度；分別為界面fi在n和n＋1時刻的流通量矢量；n為界面fi的法向量；θ為半隱權(quán)重因子.離散的動量方程為

其中，離散方程的求解分為顯式預(yù)估和隱式校正二步求解［18］；水位余量方程采用雅可比共軛梯度法求解.詳細的離散過程和求解方法參考文獻［16］.

（2）鹽度輸運方程的離散.為了使方程既能保證質(zhì)量守恒的特性，又能在離散后滿足離散極大極小原則，本模型在離散的守恒形式的對流擴散方程中引入了離散極大極小原則.根據(jù)Casulli等［19］研究，若要使守恒形式的對流擴散方程的離散格式滿足極大極小原則，那么它必須與離散的連續(xù)性方程協(xié)調(diào).因而，改寫離散的連續(xù)性方程如下：

可見直接離散的守恒形式對流擴散方程與離散的連續(xù)性方程界面通量的取值是不協(xié)調(diào)的，那么，為了使守恒形式對流擴散方程的離散形式與離散的連續(xù)性方程式協(xié)調(diào)，調(diào)整界面通量取值可得一階顯式迎風(fēng)格式下離散的鹽度輸運方程為

因為該離散格式采用了有限體積法，離散方程以界面通量的形式來表示，該格式仍具有質(zhì)量守恒的特性，而且如果時間步長足夠小，式（9）右邊實際上是的加權(quán)平均值，而它們的加權(quán)系數(shù)非負，且和為1，因而滿足離散極大極小原則，格式是穩(wěn)定的.該離散方程式對時間步長的要求滿足［20］

由于一階迎風(fēng)格式存在較大的數(shù)值耗散，因而精度不高［21］.為了提高精度，可以采用二階格式，即采用網(wǎng)格重構(gòu)的方法在單元內(nèi)以重構(gòu)的鹽度分布值代替單元平均的鹽度值，即

為消除時間步長對潮流模塊的限制，可在鹽度輸運模塊采用子時間步.每個子時間的水深通過n和n＋1時刻的水深由線性插值求得，這樣離散的連續(xù)性方程在任何子時間步內(nèi)仍滿足.

1.2.2 邊界條件

對于入流邊界，假定鹽度的輸運以對流為主，需要給出鹽度時間變化曲線；而對于出流邊界即鹽度傳出計算區(qū)域，那么可給定鹽度梯度為零的條件.

2 數(shù)值測試

鹽度輸運包含了對流和擴散2種不同的重要特性，它們的運輸形式完全不同，因而分別對對流和擴散2種輸運形式進行測試.

2.1 對流測試：彎曲河道對流

物質(zhì)在純對流作用下輸運的測試采用了最初由Casulli和Zanolli［19］設(shè)計的U形河道試驗來測試對流格式的耗散性.該U形河道長L＝8 570m（其中直線段長3 440m，彎曲段外徑560m，內(nèi)徑440m），寬W＝120m，水深恒定D＝10m.該U型河道離散成幾何尺度為20m的三角形網(wǎng)格，如圖1所示.

圖1 U形河道無結(jié)構(gòu)三角網(wǎng)格Fig.1 Unstructured triangular mesh in the U-shaped river channel

U形河道兩端給定恒定的水頭差，入流邊界給定水位0.09m，出流邊界水位為－0.09m，從而驅(qū)動一個恒定流，并設(shè)整個渠道的曼寧系數(shù)為0.04.為了保證純對流作用，計算中設(shè)置擴散系數(shù)為零.潮流模塊時間步長取300s，而鹽度輸運模塊子時間步長為30s.模擬得水動力穩(wěn)定后，在入流邊界處連續(xù)釋放10min濃度C＝1的示蹤劑，繼續(xù)模擬到所有示蹤劑隨對流流出河道.

圖2a，2b分別顯示了一階和二階迎風(fēng)格式下模擬的示蹤劑釋放75，150和225min后的示蹤劑濃度場，可見一階迎風(fēng)格式的耗散性非常強，經(jīng)過75 min后示蹤劑就已經(jīng)發(fā)生了較大的擴散，產(chǎn)生較大的數(shù)值誤差，最大濃度已由原來的1.0減小到0.4左右，而二階迎風(fēng)格式大幅降低了數(shù)值耗散，經(jīng)過U型彎道后仍能保持較小的數(shù)值擴散，精度相比一階迎風(fēng)有了很大提高，說明二階格式可以明顯減小數(shù)值耗散，但仍有微弱的數(shù)值耗散.文獻［19］認為二階格式產(chǎn)生的數(shù)值耗散可能與網(wǎng)格的扭曲有關(guān).

2.2 擴散測試：一維鹽水入侵

采用一維鹽水入侵模型測試物質(zhì)輸運模塊中鹽度的擴散輸運.一維鹽水入侵模型的控制方程為

式中：U為恒定流流速；S為鹽度；x為橫坐標(biāo)；AH為擴散系數(shù).

其邊界條件定義為

通過分析可得其解析解為

其中，S0為開邊界恒定的鹽度.

為了生成一個一維恒定流，本文采用一個均勻水深10m的5km長、1km寬的矩形河道，該河道被離散為網(wǎng)格長度為100m的三角形網(wǎng)格.在上邊界和下邊界均給定恒定流200.0m3·s－1.鹽度擴散系數(shù)取 10.0m2·s－1，而下游開邊界鹽度S0＝35.0‰，潮流模塊時間步長取60s，而鹽度輸運模塊子時間步長為10s.

由圖3可見，模型能夠很好地處理鹽水入侵即鹽度的擴散輸運.

圖3 一維鹽水入侵模擬結(jié)果與解析解的比較Fig.3 Comparison of simulated salinity and analytical solution

3 長江口鹽水入侵的數(shù)值模擬

3.1 計算范圍及網(wǎng)格

模型計算范圍西起江陰，東至外海－40m等深線，北至連興港的北側(cè)，南至南匯嘴以南，包括了南北支、南北港和南北槽在內(nèi)的整個長江口水域.其中，東西跨度大約為278km，南北跨度大約為188 km.由于長江口具有三級分汊、四口入海的不規(guī)則幾何邊界，本文采用無結(jié)構(gòu)三角網(wǎng)格擬合岸線.計算區(qū)域共生成9 202個三角形網(wǎng)格、5 035個網(wǎng)格節(jié)點，水平空間網(wǎng)格長為200～5 200m（見圖4）.

3.2 邊界條件以及相關(guān)參數(shù)選取

模型外海開邊界采用潮位過程控制，由8個主要分潮（M2，K2，S2，N2，K1，P1，O1，Q1）的調(diào)和常數(shù)計算得到，上邊界采用實測潮位過程控制.

上游開邊界鹽度取零，外海邊界鹽度由實測值插值得到.初始潮位取上游開邊界和外海東邊界潮位的平均值，初始流速和鹽度取零.糙率按水深插值取n＝0.01＋0.01／D，鹽度紊動擴散系數(shù)取100 m2·s－1，潮流模塊時間步長取60s，而鹽度輸運模塊子時間步長為20s.

圖4 計算區(qū)域及網(wǎng)格Fig.4 Computational domain and grids

3.3 模型的驗證

選用2004年5月5日5：00至5月6日11：00的實測潮位、流速和鹽度資料［7，22］對模型進行了驗證.潮位、流速和鹽度測站的具體位置如圖5所示，圖中bg，ng和cs分別為北港、南港和北槽中的測站.圖6～8分別列出了部分測站的潮位、流速和鹽度的驗證結(jié)果.從驗證結(jié)果看，潮位的計算值比實測值稍大，兩者之間的平均誤差小于0.1m；流速的平均相對誤差小于5%；流向的平均誤差小于4°；鹽度的平均誤差小于0.12%.

圖9為長江口漲落急時刻的流場.由圖可見，漲潮時，崇明東灘、橫沙東灘、南匯東灘以及九段沙部分灘地被潮水淹沒，落潮時，它們又露出水面.

圖10為長江口漲落憩時刻的鹽度場，整體上看，模型能夠基本反映長江口鹽水入侵的實際情況：鹽度隨漲潮向上游入侵，隨落潮退向海洋；北支主要為鹽水控制，而南支為河川徑流控制.

4 結(jié)論

建立了與連續(xù)性方程協(xié)調(diào)的鹽度輸運模型，既滿足了守恒性，在一定時間步長內(nèi)又不違反離散極大極小原則，保證了鹽度輸運模型的穩(wěn)定性.同時，鹽度輸運模塊也采用了網(wǎng)格模塊的重構(gòu)機制，實現(xiàn)了高精度的數(shù)值模擬.該模型分別采用彎曲河道和矩形河道測試了鹽度輸運的純對流作用和純擴散作用，驗證結(jié)果表明該數(shù)學(xué)模型計算穩(wěn)定、計算結(jié)果精度高.最后以該模型應(yīng)用于長江口鹽水入侵的模擬計算.

［1］ 肖成猷，朱建榮，沈煥庭.長江口北支鹽水倒灌的數(shù)值模型研究［J］.海洋學(xué)報，2000，22（5）：124.

XIAO Chengyou，ZHU Jianrong，SHEN Huanting.Study on numerical modeling about salt water flow backward in the Changjiang Estuary north branch［J］.Acta Oceanologica Sinica，2000，22（5）：124.

［2］ GU Jie，HUANG Jing，HAN Bing，et al.Numerical study to choose the optimum narrowing schemes for protecting water source in the south branch of the changjiang estuary［C］／／2009 3rd International Conference on Bioinformatics and Biomedical Engineering（iCBBE2009）.Beijing：IEEE Eng in Medicine and Biology Society of USA，2009：1-4.

［3］ 羅小峰，陳志昌.長江口水流鹽度數(shù)值模擬［J］.水利水運工程學(xué)報，2004（2）：29.

LUO Xiaofeng，CHEN Zhichang.Numerical simulation of salinity in Yangtze River Estuary［J］.Hydro Science and Engineering，2004（2）：29.

［4］ 匡翠萍.長江口鹽水入侵三維數(shù)值模擬［J］.河海大學(xué)學(xué)報：自然科學(xué)版，1997，25（4）：54.

KUANG Cuiping. A 3-D numerical model for saltwater intrusion in the Changjiang Estuary［J］.Journal of Hohai University：Natural Sciences，1997，25（4）：54.

［5］ KUANG Cuiping，SUN Bo，LIU Shuguang，et al.Approach of moving boundary and its application in 3D tidal current simulation of Yangtze River Estuary based on ECOMSED Model［J］.Acta Oceanologica Sinica，2009，28（4）：78.

［6］ 朱建榮，朱首賢.ECOM模式的改進及在長江河口、杭州灣及鄰近海區(qū)的應(yīng)用［J］.海洋與湖沼，2003，34（4）：364.

ZHU Jianrong，ZHU Shouxian.Improvement of the ECOM with application to the Yangtze River Estuary，Hangzhou Bay and adjacent waters［J］.Oceanologia Et Limnologia Sinica，2003，34（4）：364.

［7］ 馬鋼峰，劉曙光，戚定滿.長江口鹽水入侵數(shù)值模型研究［J］.水動力學(xué)研究與進展：A輯，2006，21（1）：53.

MA Gangfeng， LIU Shuguang， QI Dingman. Three dimensional hydronamic model of Yangtze River Estuary［J］.Journal of Hydrodynamics：Ser A，2006，21（1）：53.

［8］ 王志力，耿艷芬，金生.具有復(fù)雜計算域和地形的二維淺水流動數(shù)值模擬［J］.水利學(xué)報，2005，36（4）：439.

WANG Zhili，GENG Yanfen，JIN Sheng.Numerical modeling of 2D shallow water flow with complicated geometry and topography［J］.Journal of Hydraulic Engineering，2005，36（4）：439.

［9］ 王志力，陸永軍，耿艷芬.基于非結(jié)構(gòu)網(wǎng)格有限體積法的二維高精度物質(zhì)輸運模擬［J］.水科學(xué)進展，2008，19（4）：531.

WANG Zhili，LU Yongjun，GENG Yanfen.High resolution twodimensional numerical model of scalar transport based on unstructured grid finite volume method［J］.Advances in Water Science，2008，19（4）：531.

［10］ 耿艷芬，王志力，陸永軍.基于無結(jié)構(gòu)網(wǎng)格單元中心有限體積法的二維對流擴散方程離散［J］.計算物理，2009，26（1）：17.

GENG Yanfen，WANG Zhili，LU Yongjun.Discretization of two-dimensional advection-diffusion equation with unstructured cell center finite volume method［J］.Chinese Journal of Computational Physics，2009，26（1）：17.

［11］ 潘存鴻，徐昆.三角形網(wǎng)格下求解二維淺水方程的KFVS格式［J］.水利學(xué)報，2006，37（7）：858.

PAN Cunhong，XU Kun.Kinetic flux vector splitting scheme for solving 2D shallow water equations with triangular mesh［J］.Journal of Hydraulic Engineering，2006，37（7）：858.

［12］ 潘存鴻.三角形網(wǎng)格下求解二維淺水方程的和諧Godunov格式［J］.水科學(xué)進展，2007，18（2）：204.

PAN Cunhong. Well-balanced Godunov-type scheme for 2D shallow water flow with triangle mesh［J］.Advanced in Water Science，2007，18（2）：204.

［13］ 陳昞睿，朱建榮，吳輝，等.無結(jié)構(gòu)網(wǎng)格二維河口海岸水動力數(shù)值模式的建立及其應(yīng)用［J］.海洋學(xué)報，2010，32（2）：31.

CHEN Bingrui， ZHU Jianrong， WU Hui， et al. An unstructured mesh two-dimensional coastal and estuarine hydrodynamic numerical model and its application［J］.Acta Oceanologica Sinica，2010，32（2）：31.

［14］ Vincenzo Casulli，Roy A Walters.An unstructured grid，threedimensional model based on the shallow water equations［J］.International Journal for Numerical Methods in Fluids，2000，32：331.

［15］ SHEN Jian，WANG Harry，Mac Sisson，et al.Storm tide simulation in the Chesapeake Bay using an unstructured grid model［J］.Estuarine，Coastal and Shelf Science，2006，68：1.

［16］ CHEN Siyu，KUANG Cuiping，LIU Shuguang，et al.A semiimplicit Eulerian-Lagrangian Finite-Volume Method for coastal circulation ［C］／／Proceedings of the Fifth International Conference on Asian and Pacific Coasts，Singapore：World Scientific，2009：244-250.

［17］ LIU Wencheng，CHEN Weibo，CHENG T Ralph，et al.Modelling the impact of wind stress and river discharge on Danshui River plume［J］.Applied Mathematical Modelling，2008，32：1255.

［18］ 柯國益.河口海岸水沙運動非結(jié)構(gòu)有限體積方法數(shù)值模擬研究［D］.上海：上海交通大學(xué)船舶海洋與建筑工程學(xué)院，2008.

KE Guoyi.Numerical simulation study on finite volume method for estuarine and coastal flow and sediment movements on the unstructured grids ［D］. Shanghai： Shanghai Jiaotong University School of Naval Architecture，Ocean and Civil Engineering，2008.

［19］ Vincenzo Casulli，Paola Zanolli.High resolution methods for multidimensional advection-diffusion problems in free-surface hydrodynamics［J］.Ocean Modelling，2005，10：137.

［20］ Greenspan D，Casulli V.Numerical analysis for applied mathematics，science and engineering［M］.Boston：Addison-Wesley Longman Publishing Co.Inc.，1988.

［21］ LeVeque R J.Numerical methods for conservation laws［M］.Berlin：Birkhauser Verlag，1992.

［22］ 孫波.三峽與南水北調(diào)工程對長江口鹽水楔影響的數(shù)值研究［D］.上海：同濟大學(xué)土木工程學(xué)院，2009.

SUN Bo.Numerical study on the impact of Three Gorges and South-North Water Transfer Project on the saline wedge of the Yangtze River Estuary［D］.Shanghai：Tongji University.College of Civil Engineering，2009.