模擬二維水流運動的三次樣條多尺度有限元方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明屬于水力學(xué)技術(shù)領(lǐng)域，具體涉及一種模擬多孔介質(zhì)中二維水流運動的三次 樣條多尺度有限元方法。
【背景技術(shù)】
[0002] 地下水是十分重要的水體，是水資源的重要組成部分。世界各國有很多城市的大 部分用水都取自地下水。同時，在地質(zhì)工程活動中，地下水的分布和運移是一項重要因素。 因此，研究地下水位和滲流速度的計算方法和數(shù)值模擬，對于考察地下水分布規(guī)律及預(yù)測 具有十分重要的意義。
[0003] 多尺度有限單元法（Hou and Wu 1997)是由科學(xué)工作者提出的一種求解非均質(zhì)介 質(zhì)中地下水流問題的方法。多尺度有限單元法的核心思想是在網(wǎng)格單元上令基函數(shù)滿足局 部微分算子，從而通過基函數(shù)抓住介質(zhì)的非均質(zhì)性質(zhì)。和常規(guī)的有限單元方法或者有限差 分方法不同，多尺度有限單元法不要求單元內(nèi)部的滲透系數(shù)必須近似為常數(shù)。在處理非均 質(zhì)問題時，多尺度有限單元法可以使用較少的網(wǎng)格單元數(shù)，無需對研究區(qū)精細(xì)剖分。因此， 多尺度有限單元法比有限單元方法或有限差分方法的計算效率更高，在進(jìn)行大尺度研究區(qū) 域中的水流模擬時能夠節(jié)約大量計算成本。同時，許多科學(xué)工作者已經(jīng)從數(shù)學(xué)推導(dǎo)和數(shù)值 模擬上證明了多尺度有限單元法能夠很好的求解橢圓和拋物型問題，并且高效、收斂、精確 (Hou and Wu 1997, Hou et al. 1999, Deng et al. 2008, Ye et al. 2004)。然而，和有限單 元法一樣，多尺度有限單元法無法獲得連續(xù)的節(jié)點水頭導(dǎo)數(shù)，導(dǎo)致其無法通過達(dá)西定律獲 得連續(xù)、精確的達(dá)西滲流速度場，不能對地下水運動狀態(tài)進(jìn)行精確描述。

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

[0004] 針對于上述問題，本發(fā)明的目的在于提供一種模擬二維水流運動的三次樣條多尺 度有限元方法，該方法將三次樣條法（Zhang et al. 1994)與多尺度有限單元法有機(jī)結(jié)合， 以解決現(xiàn)有技術(shù)中無法通過達(dá)西定律獲得連續(xù)、精確的達(dá)西滲流速度場，不能對地下水運 動狀態(tài)進(jìn)行精確描述的問題。
[0005] 為達(dá)到上述目的，本發(fā)明的模擬二維水流運動的三次樣條多尺度有限元方法，包 括步驟如下：
[0006] (1)根據(jù)所要模擬的研究區(qū)域確定邊界條件，設(shè)定粗網(wǎng)格單元尺度h，剖分該研究 區(qū)域，得到粗網(wǎng)格單元，對每一個粗網(wǎng)格單元進(jìn)行細(xì)剖分得到細(xì)網(wǎng)格單元；
[0007] (2)根據(jù)滲透系數(shù)K以及基函數(shù)的邊界條件，求解退化的橢圓型問題，確定基函 數(shù)，形成有限元空間；
[0008] (3)計算各粗網(wǎng)格單元的剛度矩陣，相加得總剛度矩陣；根據(jù)研究區(qū)域的邊界條 件、源匯項，計算右端項，形成有限元方程；
[0009] (4)提供有限元方程的有效解法，求得研究區(qū)域上每個節(jié)點的水頭；
[0010] (5)在研究區(qū)每條網(wǎng)格單元線上運用三次樣條技術(shù)求解水頭的一階導(dǎo)數(shù)值；
[0011] (6)通過達(dá)西定律求解研究區(qū)節(jié)點上的達(dá)西滲流速度Γ;
[0012] (7)通過研究區(qū)域每個節(jié)點的水頭及基函數(shù)得到每一個粗網(wǎng)格單元上每一個節(jié)點 水頭；
[0013] (8)在粗網(wǎng)格每條單元線上運用三次樣條技術(shù)求解水頭的一階導(dǎo)數(shù)值；
[0014] (9)通過達(dá)西定律求解每一粗網(wǎng)格單元上的達(dá)西滲流速度Vf。
[0015] 進(jìn)一步地，所述的步驟（1)中，采用矩形單元剖分研究區(qū)域，以形成粗網(wǎng)格單元。
[0016] 進(jìn)一步地，所述的步驟（1)中，采用三角形單元剖分粗網(wǎng)格單元，以形成細(xì)網(wǎng)格單 J L· 〇
[0017] 進(jìn)一步地，所述的步驟（2)、（3)中，細(xì)網(wǎng)格單元上的滲透系數(shù)K、源匯項值近似取 這個單元的所有頂點上的滲透系數(shù)、源匯項的平均值。
[0018] 進(jìn)一步地，所述的步驟（5)、（6)與步驟（7)、（8)、（9)相互獨立。
[0019] 本發(fā)明的有益效果：
[0020] 1.本發(fā)明令多尺度有限單元法可以獲得其研究區(qū)及粗網(wǎng)格單元節(jié)點上的達(dá)西滲 流速度，并保證它們的連續(xù)性；
[0021] 2.與現(xiàn)有計算達(dá)西滲流速度的常用方法相比，在研究區(qū)網(wǎng)格節(jié)點數(shù)目相同時，本 發(fā)明獲得的精度更高；
[0022] 3.本發(fā)明可以將計算整個研究區(qū)達(dá)西滲流速度的問題轉(zhuǎn)化為求解各個粗網(wǎng)格單 元節(jié)點的達(dá)西滲流速度的問題，從而節(jié)約大量計算成本；
[0023] 4.本發(fā)明能有效求解非均質(zhì)多孔介質(zhì)地下水流場和達(dá)西滲流速度場，且原理簡 單，尚效精確；
[0024] 5.在求解大范圍，長時間等高計算量問題時，本發(fā)明的效率更高。
【附圖說明】
[0025] 圖1繪示本發(fā)明中三次樣條法示意圖。
[0026] 圖2繪示二維連續(xù)介質(zhì)模型，在研究區(qū)的每邊被剖分為10、20、30等份時各數(shù)值法 獲得的G的平均相對誤差。
[0027] 圖3繪示二維高振蕩水頭模型，各數(shù)值法獲得的水頭在截面y = 0. 7處的相對誤 差。
[0028] 圖4繪示二維高振蕩水頭模型，各數(shù)值法獲得的「在截面y = 0. 7處的相對誤差。
【具體實施方式】
[0029] 為了便于本領(lǐng)域技術(shù)人員的理解，下面結(jié)合實施例與附圖對本發(fā)明作進(jìn)一步的說 明，實施方式提及的內(nèi)容并非對本發(fā)明的限定。
[0030] 本發(fā)明的模擬二維水流運動的三次樣條多尺度有限元方法求解粗網(wǎng)格單元上的 達(dá)西滲流速度Vf的過程如下：
[0031] 先通過多尺度有限單元法得到粗網(wǎng)格單元頂點水頭，再通過基函數(shù)的插值公式獲 得所有粗網(wǎng)格單元上的節(jié)點（即細(xì)尺度節(jié)點）上的水頭。假設(shè)[a，b]是粗網(wǎng)格單元口 1]H上的一條網(wǎng)格線（參照圖1所示），有n+1個節(jié)點，使用樣條函數(shù)S近似估計水頭H，樣條函 數(shù)S為三次多項式，S在[a，b]上任一節(jié)點的值等于該點的水頭值；S的二階導(dǎo)數(shù)在[a，b] 上連續(xù)，S在[a, b]上某區(qū)間[χτ χτ]的表達(dá)式如下：
[0033] 由于S的二階導(dǎo)數(shù)在[a，b]上連續(xù)，可以得到一個關(guān)于S-階導(dǎo)數(shù)的方程組，求解 可以得到S的一階導(dǎo)數(shù)，該導(dǎo)數(shù)在節(jié)點上是連續(xù)的。利用該導(dǎo)數(shù)替代水頭的一階導(dǎo)數(shù)，再通 過達(dá)西定律，可以得到粗網(wǎng)格單元節(jié)點上的連續(xù)的達(dá)西滲流速度vf。若需要求解研究區(qū)域 網(wǎng)格節(jié)點（即粗尺度節(jié)點）上的達(dá)西滲流速度V%只需利用多尺度有限單元法求得的所有 研究區(qū)網(wǎng)格節(jié)點上的水頭值，再通過三次樣條函數(shù)得到連續(xù)的水頭導(dǎo)數(shù)，通過達(dá)西定律求 解即可。具體步驟和上述求解，的步驟類似。通過對多孔介質(zhì)下的二維連續(xù)介質(zhì)模型，二 維高振蕩水頭模型，二維振蕩介質(zhì)模型，二維多尺度介質(zhì)模型（兩種不同尺度），二維漸變 介質(zhì)模型，二維穩(wěn)定流潛水介質(zhì)模型的數(shù)值模擬，發(fā)現(xiàn)本發(fā)明與解析解吻合的很好；和常用 求解達(dá)西滲流速度的方法相比，在研究區(qū)網(wǎng)格節(jié)點數(shù)目相同時，本發(fā)明獲得的精度更高；在 總節(jié)點數(shù)目（包括細(xì)尺度節(jié)點）相同，本發(fā)明精度和現(xiàn)有常用方法相近，并可以將計算整個 研究區(qū)達(dá)西滲流速度的問題轉(zhuǎn)化為求解各個粗網(wǎng)格單元節(jié)點的達(dá)西滲流速度的問題，節(jié)約 大量計算成本。
[0034] 下面結(jié)合具體實施例對本發(fā)明做進(jìn)一步的解釋，其中涉及一些簡寫符號，以下為 注解：
[0035] M:研究區(qū)邊界被等分份數(shù)
[0036] Vx:X方向上的達(dá)西滲流速度
[0037] 巧:研究區(qū)網(wǎng)格節(jié)點上的達(dá)西滲流速度Vx，即粗尺度達(dá)西滲流速度
[0038] F/ :粗單元網(wǎng)格節(jié)點上的達(dá)西滲流速度Vx，即細(xì)尺度達(dá)西滲流速度
[0039] LFEM:傳統(tǒng)有限單元法，采用三角形網(wǎng)格單元
[0040] LFEM-F:精細(xì)剖分的傳統(tǒng)有限單元法，采用三角形網(wǎng)格單元
[0041] CMSFEM :三次樣條多尺度有限單元法，采用矩形網(wǎng)格單元
[0042] CST-AS :CMSFEM應(yīng)用解析的水頭求解細(xì)尺度達(dá)西滲流速度
[0043] CMSFEM-tr :三次樣條多尺度有限單元法，采用三角形網(wǎng)格單元
[0044] Method_Yeh:Yeh的線性伽遼金有限元模型，采用LFEM求解水頭
[0045] Method-Yeh-F: Yeh的線性伽遼金有限元模型（精細(xì)剖分），采用LFEM-F求解水頭
[0046] Method-Zhang:三次樣條法，采用LFEM求解水頭
[0047] CMSFEM和CMSFEM-tr的基函數(shù)均使用振蕩邊界條件，由于三次樣條技術(shù)要求應(yīng)用 區(qū)域為矩形，CMSFEM-tr僅能夠計算粗尺度上的達(dá)西滲流速度。在求解時，CMSFEM是在 研究區(qū)網(wǎng)格上求解的，而計算F/時，CMSFEM是在每個粗網(wǎng)格單元上求解的。在計算人的相 對誤差時，若在某節(jié)點Vx，且數(shù)值解Vx< 10 則將該節(jié)點的相對誤差記為0%，否則記為 100%〇
[0048] 實施例1 :二維連續(xù)介質(zhì)模型
[0049] 研究區(qū)為一正方形單元：Ω = [0, lm] X [0, lm]，滲透系數(shù) K(x，y) = (1+x) (1+y) m/d，水流方程為：
[0051] 邊界條件為定水頭邊界條件
^此模型有解析解：H = Xy(l-X) (l_y)〇
[0052] 求解粗尺度達(dá)西滲流速度：
[0053] 采用CMSFEM、Method-Yeh和Method-Zhang求解此模型。為了保證研究區(qū)上網(wǎng) 格節(jié)點數(shù)目相同，CMSFEM、Method-Yeh和Method-Zhang均將研究區(qū)的每邊剖分為相同份 數(shù)。三種方法均將研究區(qū)每邊剖分為1〇、20、30等份（M= 10、20、30)。Method-Yeh和 Method-Zhang采用三角形剖分，則將研究區(qū)剖分為200、800、1800個三角形單元；CMSFEM采 用正方形剖分，將研究區(qū)剖分為1〇〇、400、900個矩形粗網(wǎng)格單元，并將每個粗網(wǎng)格單元剖 分為50個三角形單元。
[0054] Method-Yeh和Method-Zhang通過有限單元法（LFEM)求解水頭，當(dāng)研究區(qū)的每邊 被分為10、20、30等份時，水頭的平均相對誤差分別為2. 29 %、0. 58 %和0. 26%，CPU時間 分別為1秒、1秒、2秒。CMSFEM采用多尺度有限單元法（MSFEM)求解水頭，當(dāng)研究區(qū)的每 邊被分為10、20、30等份時，水頭的