模擬多孔介質(zhì)中二維水流運(yùn)動(dòng)的高效多尺度有限元方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明屬于水力學(xué)技術(shù)領(lǐng)域，具體涉及一種模擬多孔介質(zhì)中二維水流運(yùn)動(dòng)的高效 多尺度有限元方法。
【背景技術(shù)】
[0002] 地下水資源是水資源的重要組成部分，是工業(yè)、農(nóng)業(yè)和城市用水的重要來源之一。 在水文地質(zhì)學(xué)中地下水位能夠反映地下水所具有的機(jī)械能的大小。地下水的分布情況與工 程的實(shí)施方案、施工方式、施工時(shí)間、工程資金等因素關(guān)系密切;所以，研究地下水位相關(guān)的 數(shù)值計(jì)算方法，對于分析地下水分布、運(yùn)動(dòng)情況是非常必要，具有重要的研究價(jià)值。
[0003] 傳統(tǒng)有限單元法是目前常用的地下水?dāng)?shù)值計(jì)算方法之一，在水文地質(zhì)領(lǐng)域應(yīng)用十 分廣泛。然而，該方法受到單元內(nèi)部巖性必須相同的限制，在處理非均質(zhì)多孔介質(zhì)中的地下 水問題時(shí)，必須通過精細(xì)剖分的方式保證精度，在研究區(qū)域較大時(shí)，需要大量的計(jì)算時(shí)間和 空間，效率較低。因此，國外數(shù)學(xué)家提出了多尺度有限單元法(Hou and Wu 1997)用于解決 這一問題，是一種既可以減少剖分單元數(shù)又能保證計(jì)算結(jié)果精度的高效方法。該方法突破 了有限單元法的單元內(nèi)部巖性必須相同的限制，在求解非均質(zhì)多孔介質(zhì)中的地下水問題 時(shí)，該方法可以使用較大的單元。同時(shí)，該方法還能夠通過基函數(shù)滿足粗網(wǎng)格單元的局部橢 圓問題來抓住細(xì)尺度的巖性性質(zhì)，從而節(jié)約計(jì)算消耗并保證精度(Hou and Wu 1997，Ye et al.2004，Xie et al2014)。鑒于上述傳統(tǒng)有限單元法的局限性和多尺度有限單元法的優(yōu)越 性，開展多尺度有限單元法的改進(jìn)算法具有重要的理論和實(shí)際意義。近年來，經(jīng)濟(jì)、科技發(fā) 展迅速，人們對于地下水問題越來越關(guān)注，也謀求求解更加復(fù)雜的地下水問題，相應(yīng)研究面 積、水文周期也越來越大，如地面沉降問題、海水入侵問題等。在求解此類問題時(shí)，多尺度有 限單元法的基函數(shù)構(gòu)造消耗過高，需要進(jìn)一步提高。

【發(fā)明內(nèi)容】

[0004] 針對于上述現(xiàn)有技術(shù)的不足，本發(fā)明的目的在于提供一種模擬多孔介質(zhì)中二維水 流運(yùn)動(dòng)的高效多尺度有限元方法，該方法運(yùn)用區(qū)域分解技術(shù)改進(jìn)了基函數(shù)的構(gòu)造算法和剖 分方式，將基函數(shù)的構(gòu)造問題分解為若干個(gè)子問題，分批求解未知項(xiàng)，能夠大幅降低構(gòu)造基 函數(shù)所需的計(jì)算消耗，以解決現(xiàn)有技術(shù)中求解復(fù)雜的地下水問題時(shí)基函數(shù)構(gòu)造消耗過高等 問題。
[0005] 為達(dá)到上述目的，本發(fā)明的一種模擬多孔介質(zhì)中二維水流運(yùn)動(dòng)的高效多尺度有限 元方法，包括步驟如下：
[0006] (1)根據(jù)所要模擬的研究區(qū)域確定邊界條件，設(shè)定粗網(wǎng)格單元尺度，剖分該研究區(qū) 域，得到粗網(wǎng)格單元；
[0007] (2)設(shè)定中網(wǎng)格單元尺度，剖分上述粗網(wǎng)格單元，得到中網(wǎng)格單元；
[0008] (3)設(shè)定細(xì)網(wǎng)格單元尺度，剖分上述中網(wǎng)格單元，得到細(xì)網(wǎng)格單元；
[0009] (4)根據(jù)滲透系數(shù)K以及基函數(shù)的邊界條件，以中網(wǎng)格單元為最小子單元，在粗網(wǎng) 格單元上求解退化的橢圓型問題，確定所有中網(wǎng)格單元頂點(diǎn)處的基函數(shù)值；
[0010] (5)運(yùn)用區(qū)域分解技術(shù)將上述粗網(wǎng)格上的局部橢圓問題分解為每個(gè)中網(wǎng)格單元上 的子問題；
[0011] (6)根據(jù)滲透系數(shù)K、中網(wǎng)格單元頂點(diǎn)處的基函數(shù)值以及改進(jìn)的基函數(shù)邊界條件得 到所有子問題的邊界條件，以細(xì)網(wǎng)格單元為最小子單元，在每個(gè)中網(wǎng)格單元上求解子問題 得到基函數(shù)在每個(gè)中網(wǎng)格單元中所有節(jié)點(diǎn)上的值；
[0012] (7)計(jì)算各粗網(wǎng)格單元的剛度矩陣，相加得總剛度矩陣；根據(jù)研究區(qū)域的邊界條 件、源匯項(xiàng)，計(jì)算右端項(xiàng)，形成有限元方程；
[0013 ] (8)采用cho 1 e shy分解法，求得研究區(qū)域上每個(gè)節(jié)點(diǎn)的水頭。
[0014] 優(yōu)選地，上述的步驟(1)中，采用三角形單元剖分研究區(qū)域，以形成粗網(wǎng)格單元。
[0015] 優(yōu)選地，上述的步驟(2)中，采用三角形單元剖分粗網(wǎng)格單元，以形成中網(wǎng)格單元。
[0016] 優(yōu)選地，上述的步驟(3)中，采用放射狀的三角形單元剖分方式剖分中網(wǎng)格單元， 以形成細(xì)網(wǎng)格單元。
[0017]優(yōu)選地，上述的步驟(4)中，所述的研究區(qū)域最小子單元上的滲透系數(shù)K取這個(gè)單 元的所有頂點(diǎn)上的滲透系數(shù)平均值。
[0018] 優(yōu)選地，上述的步驟(6)中，所述的研究區(qū)域最小子單元上的滲透系數(shù)K取這個(gè)單 元的所有頂點(diǎn)上的滲透系數(shù)平均值。
[0019] 優(yōu)選地，上述的步驟(7)中，細(xì)網(wǎng)格單元上的源匯項(xiàng)值取這個(gè)單元的所有頂點(diǎn)上的 源匯項(xiàng)的平均值。
[0020] 優(yōu)選地，上述的步驟(7)中，粗網(wǎng)格單元上的源匯項(xiàng)值取這個(gè)單元中所有細(xì)網(wǎng)格單 元的源匯項(xiàng)的平均值。
[0021] 本發(fā)明的有益效果：
[0022] 1.令多尺度有限單元法可以將其每個(gè)粗網(wǎng)格單元上的基函數(shù)構(gòu)造問題轉(zhuǎn)換為若 干個(gè)子問題，從而分批求解未知項(xiàng)，大幅降低基函數(shù)的構(gòu)造消耗；
[0023] 2.可以將整個(gè)研究區(qū)剖分為粗、中、細(xì)三重網(wǎng)格單元，單元靈活度更高，單元的抗 畸形能力較強(qiáng)；
[0024] 3.改進(jìn)了基函數(shù)的邊界條件，不要求邊界端點(diǎn)必須分別為0和1;
[0025] 4.在研究區(qū)粗網(wǎng)格數(shù)目相同且粗網(wǎng)格的細(xì)網(wǎng)格數(shù)目相同時(shí)，本發(fā)明獲得的精度和 傳統(tǒng)多尺度有限單元法相近，計(jì)算時(shí)間大幅降低；
[0026] 5.能夠有效求解非均質(zhì)多孔介質(zhì)地下水穩(wěn)定流和非穩(wěn)定流問題，如二維穩(wěn)定流的 連續(xù)介質(zhì)問題、二維非穩(wěn)定流的漸變介質(zhì)問題，二維穩(wěn)定流的突變介質(zhì)問題，二維穩(wěn)定流的 多尺度介質(zhì)模型(兩種不同尺度），二維潛水流模型的數(shù)值模擬，且原理簡單，高效精確；
[0027] 6.在求解大尺度，長時(shí)間等高計(jì)算量問題時(shí)，本發(fā)明的效率更高。
【附圖說明】
[0028] 圖1為高效多尺度有限單元法的粗網(wǎng)格單元剖分方式示意圖。
[0029] 圖2為高效多尺度有限單元法的中網(wǎng)格單元剖分方式示意圖。
[0030]圖3為二維穩(wěn)定流的連續(xù)介質(zhì)模型，子例1.1中各數(shù)值方法在y=100m截面處的水 頭絕對誤差值的示意圖。
[0031] 圖4為二維非穩(wěn)定流的漸變介質(zhì)模型，子例2.1中各數(shù)值方法在y = 5200m截面處的 水頭曲線的不意圖。
[0032] 圖5為二維非穩(wěn)定流的漸變介質(zhì)模型，子例2.2中各數(shù)值方法所計(jì)算的在y = 5200m 截面處的水頭曲線的示意圖。
[0033] 圖6為二維潛水流模型，各數(shù)值方法在y = 0.5截面處水頭的相對誤差值的示意圖。
【具體實(shí)施方式】
[0034] 為了便于本領(lǐng)域技術(shù)人員的理解，下面結(jié)合實(shí)施例與附圖對本發(fā)明作進(jìn)一步的說 明，實(shí)施方式提及的內(nèi)容并非對本發(fā)明的限定。
[0035] 參照圖1所示，本發(fā)明的一種模擬多孔介質(zhì)中二維水流運(yùn)動(dòng)的高效多尺度有限元 方法，其基函數(shù)構(gòu)造過程如下：
[0036]在粗網(wǎng)格單元Δ,訴上考慮構(gòu)造方程，即退化的橢圓方程：
[0037] V · Α'νψ, = 0,
[0038] 其中Κ為滲透系數(shù)，Ψi為粗網(wǎng)格單元Δ ijk(圖1)在i點(diǎn)的基函數(shù)。
[0039] 步驟一:在粗網(wǎng)格單元△#(圖1)上，以中網(wǎng)格單元為最小單元，求解構(gòu)造方程，可 以得到基函數(shù)在0點(diǎn)（即Mo點(diǎn)）的值;通過基函數(shù)的邊界條件可以得到基函數(shù)在其余中網(wǎng)格 頂點(diǎn)（即i，j，k，e，f，p，q)上的值。
[0040]步驟二:運(yùn)用區(qū)域分解技術(shù)將Δ uk上的構(gòu)造問題轉(zhuǎn)換為中網(wǎng)格單元Δ Qgl、Δ _、 Δ。#、AQfg、Δ<^、AqkgApfe、Apej上的子問題；中網(wǎng)格單元的某一邊界ξτι上的改進(jìn)的基函數(shù) 線性、振蕩邊界條件如下：
[0043] 步驟三:在每個(gè)中網(wǎng)格單元上（圖2)，以細(xì)網(wǎng)格單元為最小單元，求解子問題即可 得到基函數(shù)在內(nèi)點(diǎn)Μι_Μ8的值。
[0044] 其中：
[0045] 在步驟一中，在通過基函數(shù)的邊界條件獲得基函數(shù)在中網(wǎng)格頂點(diǎn)上的值時(shí)，可以 對邊界加密提升所獲精度。
[0046] 在步驟二中，中網(wǎng)格單元數(shù)目并不僅限于圖1所示的8個(gè)，可以根據(jù)需要增多。
[0047] 完成上述基函數(shù)的構(gòu)造后，通過變分形式即可求得單元?jiǎng)偠染仃?，相加得總剛?矩陣，運(yùn)用cholesky分解法求解總剛度矩陣和右端項(xiàng)的方程組即可獲得水頭值。
[0048] 下面結(jié)合具體實(shí)施例對本發(fā)明做進(jìn)一步的解釋，其中涉及一些簡寫符號，以下為 注解:LFEM:有限單元法；
[0049] LFEM-F:有限單元法(精細(xì)剖分）；
[0050] MSFEM-L:傳統(tǒng)多尺度有限單元法，使用線性邊界條件；
[0051 ] MSFEM-0:傳統(tǒng)多尺度有限單元法，使用振蕩邊界條件；
[0052] EMSFEM-L:高效多尺度有限單元法，使用線性邊界條件；
[0053] EMSFEM-0:高效多尺度有限單元法，使用振蕩邊界條件。
[0054]本發(fā)明所有實(shí)施例在基函數(shù)構(gòu)造的步驟一中，利用粗網(wǎng)格的基函數(shù)邊界條件得到 中網(wǎng)格頂點(diǎn)基函數(shù)值時(shí)，將粗網(wǎng)格每邊分為8份，并使用振蕩邊界條件。
[0055] 實(shí)施例1:二維穩(wěn)定流的連續(xù)介質(zhì)模型
[0056] 研究區(qū)為正方形區(qū)域：Ω =[50m，150m] X [50m，150m]，滲透系數(shù)K(x，y)=x2m/d， 研究方程為穩(wěn)定流方程：
[0058] 邊界條件為定水頭邊界條件//|iu = V -3/ <源匯項(xiàng)為〇，此模型有解析解:!1 = ￥- 3y2〇
[0059] 子例 1 · 1:采用 LFEM，LFEM-F，MSFEM-L，MSFEM-0，EMSFEM-L 和 EMSFEM-0 求解。LFEM-F 將研究區(qū)剖分為88200個(gè)單元，其他