【正文】 = q （ 11） ?隔水邊界時 q=0，計算時將（ 11）式所得邊界方程與內結點方程聯(lián)立求解，即可求得各結點水頭值。 Hn??1,knjH? 11,knjH??LOGO 第三章 地下水資源數(shù)量數(shù)值計算 ?七、例題 一維承壓含水層非穩(wěn)定滲流數(shù)學模型有限差分法。 ? 承壓含水層均質等厚，各向同性，頂板和底板水平且無限延伸，含水層水頭與兩側河流水位相同。由于某種原因，兩條對稱邊界河流水位突然下降至 0－ 0，求解含水層的水頭隨時間下降的過程。 ? ? ＋ ε ＝ u* (0＜ x＜ L ， t＞ 0) ? ? H(x ,t)＝ Ha (0≤x≤ L ， t=0) ? H(x ,t)＝ Hb＝ 0 ( x＝ 0， x＝ L， t＞ 0) ? 解： ? θ [Ti+1/2（ － ）－ Ti1/2( － )]＋ ? (1 － θ )[Ti+1/2( － )－ Ti1/2( － )＋ε iΔ x178。 ? = uiΔ x178。( － )/Δ tK )Hxx????（ T Ht??11kiH?? 1kiH? 1kiH? 11kiH??1kiH? kiH kiH 1kiH?1kiH? kiHLOGO 第三章 地下水資源數(shù)量數(shù)值計算 ?當 θ =0、 1/ 1分別為顯式、中心差分和隱式差分格式。 ? ai ＋ bi ＋ ci ＝ di ? i= ...n ?區(qū)域離散化： Δ x＝ L/(n+1)，將區(qū)域（ 0～ L） n+1等份， 0結點和 n結點位于左端和右端邊界上，邊界結點 H0＝ Hn+1＝ Hb。 ?將差分方程依次用于全部內結點則有方程組： ? a1 ＋ b1 ＋ c1 ＝ d1 ? a2 ＋ b2 ＋ c2 ＝ d2 ? ? an ＋ bn ＋ ＝ dn ?將邊界水頭 和 代入第一個方程和第 n個方程得 11kiH?? 1kiH? 11kiH??10kH? 11kH? 12kH?11kH? 12kH? 13kH?11knH?? 1knH? 11knH??10kH? 11knH??LOGO 第三章 地下水資源數(shù)量數(shù)值計算 ? ? b1 ＋ c1 ＝ d1a1 ＝ d1185。 ? an ＋ bn ＝ dn ＝ dn185。 ? 代入上述方程組得： ? b1 ＋ c1 ＝ d1185。 ? a2 ＋ b2 ＋ c2 ＝ d2 ? an1 ＋ bn1 ＋ 1 ＝ dn1 ? ? an ＋ bn ＝ dn185。 ? 12kH?11kH?11knH?? 1knH? 11knH??10kH?11kH? 12kH?11kH? 12kH? 13kH?11knH?? 1knH?12knH?? 11knH?? 1knH?LOGO 第三章 地下水資源數(shù)量數(shù)值計算 ? 矩陣形式： ? ? ? ? ＝ ? 系數(shù)矩陣 A中元素都位于三條主對角線上，其他元素均為零，為三對角矩陣，可用追趕法求解方程組。 112222220... ...0nnnnnbca b ca b cab?????????????????????????1112111...KKKnknHHHH???????????????????????????????1112...nndddd?????????????????????????LOGO 第三章 地下水資源數(shù)量數(shù)值計算 ?八、例題 二維承壓含水層非穩(wěn)定滲流數(shù)學模型有限差分法 ?ai,j ＋ bi,j ＋ ci,j ＋ ei,j ＋ fi,j ＝ di,j ?將邊界條件和初始條件代入方程構成一封閉的線性代數(shù)方程組。 ?可采用常規(guī)的線性方程組求解方法求解。但上式是五對角線性方程組即系數(shù)矩陣為大型稀疏陣，為提高計算效率可采用超松弛法或交替方向隱式差分法求解。 ?九、潛水地下水非穩(wěn)定滲流數(shù)學模型有限差分法 ? 由于 T=KHp， Hp為潛水含水層平均水位，由于未知所以潛水地下水運動差分方程為非線性方程組，在求解時首先進行線性化，一般采用先給定各結點一個水位值，由其確定各相應的 T值，則差分方程化為線性方程組。采用超松弛法或交替方向隱式差分法求解，用所得結果與假定值比較，若精度滿足要求，計算結束，否則用所得結果重新計算 T值，形成一個迭代過程。 11,kijH??1,kijH? 11,kijH?? 1,1kijH ?? 1,1kijH ??LOGO 第三章 地下水資源數(shù)量數(shù)值計算 ?第三節(jié) 地下水資源數(shù)量計算有限元法 ?一、基本原理與求解思路 ? 利用 剖分插值 把描述地下水運動的偏微分方程及其求定解問題，離散為求解線性代數(shù)方程組的問題。 ? 一般將地下水區(qū)域離散化三角形網(wǎng)格狀的小單元和有限個結點，并對時間坐標離散化；采用簡單的插值函數(shù)近似表示單元上的水頭分布，以此求偏微分方程的加權積分，然后離散集合起來形成線性代數(shù)方程組；對邊界條件進行處理，并對井群處理，解代數(shù)方程組；可求解各個結點水頭。 ?（一）區(qū)域剖分并編號 ? 根據(jù)計算區(qū)域的形狀、水文地質條件、邊界條件、補排狀況以及計算精度要求，按一定的規(guī)則將計算區(qū)域化成有限個小區(qū)域。從邊界開始剖分為 P個三角形單元，三角形的公共交點為結點，共計 m個；未知水頭結點 n個（其中內結點 n1個， n2個二類邊界點），已知水頭結點 n3個。按順序依次編內結點號、二類邊界點號（邊界流量未知的結點）、一類邊界點號（水位變化已知的結點）。 LOGO 第三章 地下水資源數(shù)量數(shù)值計算 ?（二）剩余量加權積分法 ?函數(shù) L(U)＝ 0定義在二維區(qū)域 Ω上 ,設平面插值函數(shù) ? ? u＝ α iΦi 是方程的近似解，且滿足邊界條件， ? 則 L(u)＝ R≠0 ， R稱為剩余量。 ?α i為待定系數(shù)； Φi為按某種要求選定的基函數(shù)， Φ Φ2， ? ...,Φn線性獨立。 ?選用與系數(shù) α i相對應的權函數(shù) Wi，使剩余量 R在區(qū)域 Ω上 ? 的加權積分等于零。 ? ＝ 0 (i=1,2,...， n) ?伽遼金取權函數(shù)選用基函數(shù) Φi，則 ＝ 0 ? (i＝ 1,2,...， n) ?解方程可求待定系數(shù) α i,進而可求近似解。 ni?R W id xdy???R id xdy????LOGO 第三章 地下水資源數(shù)量數(shù)值計算 ?(三） 基函數(shù) Φi的構成 ?設三角形單元 e的三個結點按順時針編號 i, j, k，其坐標（ xi ,yi ) ,(xj ,yj ) ,(xk ,yk )。單元形狀函數(shù)： ? φ i ( X ,Y )=α i +β ix +γ i y ? φ j ( X ,Y )=α j +β jx +γ j y ? φ k ( X ,Y )=α k +β kx +γ k y ?單元形狀函數(shù)的特點：在自身結點上取值為 1，在對邊結點上取值為 0 ，在單元內部取值按線性變化。 ? φ i ( Xi ,Yi )＝ α i +β ixi +γ i yi＝ 1 ? φ i ( Xj ,Yj )＝ α j +β ixj +γ i yj＝ 0 ? φ i ( Xk ,Yk )＝ α j +β ixk +γ i yk＝ 0 LOGO 第三章 地下水資源數(shù)量數(shù)值計算 ? 按克萊姆法則求 α i , β i , γ i ? φ i(x ,y)＝ (ai＋ bix＋ ciy)/2Δ e ? 同理 φ j(x ,y)＝ (aj＋ bjx＋ cjy)/2Δ e ? φ k(x ,y)＝ (ak＋ bkx＋ cky)/2Δ e ? φ i(x ,y)＋ φ j (x ,y)＋ φ k (x ,y)＝ 1 ? ai、 bi、 ci為坐標輪減數(shù)， Δ e 為單元 e的面積 。 ?區(qū)域 Ω內任意結點 i可能是幾個 (p1）三角形單元的公共頂點， ?此結點 i的基函數(shù) Φi(x,y)是各三角形單元形狀函數(shù)所組成。 ? Φi(x,y)＝ (x,y) x,y∈e1 ? (x,y) x,y∈e2 ? ...... ? (x,y) x,y∈ep1 1ei?2ei?1epi?LOGO 第三章 地下水資源數(shù)量數(shù)值計算 ? 基函數(shù) Φi（ x,y）的特點：在結點 i上取值為 1 ，在相鄰結點上取值為 0 ，在 p1個三角形單元圍成的區(qū)域內按線性變化。 ?（四） 迦遼金有限元基本方程 ? 區(qū)域 Ω上的任一點水頭采用插值函數(shù)近似表達： ? (x,y,t)＝ Hj(t)Φj(x,y) ?在區(qū)域內任一單元 e內，其任一點的水頭為： ? (x,y,t)＝ Hi(t) (x,y)+Hj(t) (x,y)+Hk(t) (x,y) ?則迦遼金有限元基本方程（單元 e內）： ? ＝ ＝ ＝ 0 ?如區(qū)域上有 n個未知水頭，則可表達為 ? ＝ （ i=1,2,...,n) mj?H()eL H W id x d y??? ()eL H idxdy???? () eeL H idx dy????eH ei?ej?ek?()eL H W id x d y??? 1()mjjJeL H id x d y???????LOGO 第三章 地下水資源數(shù)量數(shù)值計算 ? 若區(qū)域上有 p個單元，則有： ? ＝ 0 （ i=1,2,...,n) ? 該式為伽遼金有限元基本方程式。 ?（五）承壓二維非穩(wěn)定流迦遼金有限元解法 ? ＋ ＋ ε ＝ u* （ 1） ? H(x,y,0)＝ Ha（ x,y) x,y∈ Ω ,t=0 ? H(x,y,t)＝ H185。(x,y,t) x,y∈ Γ 1 ,t0 ? ? ∣ Γ 2 ＝ q(x,y,t) n為外法線方向， x,y∈ Γ 2 ,t0 ? ? 1()mjjJeL H id x d y???????1pe??()HT xx????()HT yy????Ht??HTn??LOGO 第三章 地下水資源數(shù)量數(shù)值計算 ?解 ?區(qū)域剖分為 p個單元， m個結點， n1個內結點， n2個邊界點， ? 未知水