基于有限元形函數的鉆孔水頭插值方法

錢武文，肖 芳

(南昌工程學院 水利與生態工程學院，江西 南昌 330029)

1 研究背景

針對水利工程的滲流場分析，廣泛采用數值分析的方法，利用已積累的水文地質資料將研究區域內的工程巖體劃分為若干個相同參數的材料分區[1]。在實際工程中，由于在關鍵位置通常會布置一定量的觀測孔用于觀測地下水位，假定一組滲透系數計算觀測孔位置處的水頭值，通過評價模型計算值與實測水頭值的匹配程度的方式反演出各分區的滲透系數[2-5]。此外，參數反演法具有比現場試驗法快速經濟、代表性好、應用性廣以及可靠性高等優點[6]。

使用有限元法對地下水系統進行數值模擬時，為了評價模型的模擬精度，往往需要計算鉆孔處的水頭以便與其實測值比較。當鉆孔點不在網格節點上時，為了避免重建模型，可由已知的網格節點的值插值計算出模型中任意點的值。Kriging插值法[7]、徑向基函數插值法[8]、有限元插值函數法[9-10]等插值方法都可實現此功能。其中，Kriging模型和徑向基函數法都具有良好的非線性函數擬合特性，但僅適用于變化平緩的情況，當相鄰節點的值發生較大變化時，二者并不適用[11-12]。此外，當有限元節點數量較多時，二者都存在空間待插值節點數量過大的問題。由于有限元插值函數法使用有限元單元形函數作為插值函數，因此具有計算簡單、易于實現以及精度高等優點。

工程實踐中廣泛使用等參單元進行有限元計算，以便于在模擬任意形狀和方位的模型時均能采用標準格式的數值積分法[13]。然而，已知的實際觀測點位置多為有限元模型的整體坐標，需要執行等參逆變換以計算觀測點對應的局部坐標，等參逆變換本質上為求解非線性方程組的解[14-16]。采用解析方法實現等參逆變換既不直觀也不容易，錢向東等[10]采用Taylor展開技術并作了較大程度的簡化處理構造了一種等參元逆變換算法，朱以文等[14]指出該算法在處理靠近單元邊界處的插值計算時難以收斂，并提出了一種等參元逆變換的高效算法，但僅適用于平面四節點單元和空間八節點單元。為了實現一般等參元的逆變換，本文使用改進差分進化算法估計插值點的局部坐標。由于反演維數很小(維數≤3)，從而保證了該反演過程具有較高的反演效率和快速計算能力。在計算局部坐標前，需要確定插值點位于哪個單元之內，以便獲取該單元的節點水頭和單元形函數。查閱文獻發現，尚無系統給出使用有限元插值函數法插值觀測點水頭的報道。

本研究使用有限元插值函數法插值觀測點水頭，給出了判斷插值點與單元的拓撲關系的方法，以及由整體坐標計算單元局部坐標的反演方法，并使用二維和三維模型測試了所提出方法的插值性能。

2 研究方法

2.1 插值原理

空間單元下，有限元單元插值過程如下式所示：

(1)

式中：n為單元的節點總數；U(x,y,z)為坐標(x,y,z)的待插值點對應的水頭；Uje為有限元中第j個單元節點的水頭值;e為單元標識；ψj(x,y,z)為第j個單元節點對應的形函數。

工程實踐中廣泛使用等參單元進行有限元計算，以便于在模擬任意形狀和方位的模型時仍能采用標準格式的數值積分法。以空間有限元問題為例，通過等參變換，單元的局部坐標(ξ,η,ζ)可轉換為全局坐標(x,y,z)：

(2)

相應地，以局部坐標(ξ,η,ζ)表示的插值公式為：

U(ξ,η,ζ)≈Ue(ξ,η,ζ)

(3)

顯然，等參元是以局部坐標表示的單元形函數，已知局部坐標即可根據公式(3)求出該位置的水頭值。一般鉆孔位置由全局坐標表示，而局部坐標難以由全局坐標顯式表示，需要由公式(2)進行等參逆變換。因此，可采用優化方法由全局坐標間接求得局部坐標?？捎煤薪^對值的表達式表示反演結果與實際值的匹配程度，以插值點的全局坐標的估計值與已知值的絕對誤差之和為目標函數(E)，即：

(4)

若已知插值點所處單元的節點坐標和形函數時，可通過改變(ξ,η,ζ)使目標函數達到最小的方式估計插值點對應的局部坐標(ξ,η,ζ)。則公式(4)在約束條件下可表示為如下優化問題：

(5)

2.2 確定插值點的歸屬單元

以上優化問題都是建立在已知插值點的單元歸屬的前提下，因為需要由插值點所處的單元確定所采用的單元形函數。確定插值點的歸屬單元，即判斷點與平面單元或空間單元的拓撲關系，屬于計算機圖形學的范疇。

在計算機圖形學中，常用射線法[17-18]、叉積判斷法[19]、角度法和面積(體積)法[20]確定點與多邊形或多面體的拓撲關系。由于有限元單元屬于凸多邊形(多面體)，這些方法均可正確判定。本文采用相對簡單的射線法和叉積判斷法分別判定點與平面單元和空間單元的拓撲關系。

由于判定點與一維線單元的拓撲關系非常簡單，本文不再贅述，僅介紹射線法和叉積判斷法。由于有限元模型中單元眾多，并非所有單元都有必要進行判定，因此，在判定插值點與單元拓撲關系前，首先應篩選處于插值點(記為P點)坐標范圍內的單元，然后再進行識別操作。

2.2.1 射線法 方法描述：從插值點引一條射線穿過平面單元，若該射線與單元所有邊的交點數目為奇數，則說明插值點在單元內部，否則在單元外部。圖1為射線法識別插值點與平面單元(三角形和四邊形單元)的拓撲關系示意圖。

圖1 射線法判定點與三角形和四邊形單元的拓撲關系示意圖

2.2.2 叉積判斷法 在數學上，兩向量點乘的結果為標量，可用于判定兩向量的空間方向，而兩叉乘的結果為垂直于這兩個向量且方向服從右手螺旋法則的向量。向量點乘和叉乘的示意圖如圖2所示。因此，可將向量點乘和叉乘組合起來判定點與多面體的拓撲關系。

圖2 兩向量叉乘、點乘的空間關系以及單元節點的編號順序

以有限元四面體和六面體單元為例，令其節點編號順序如圖2所示，則點與六面體單元和四面體單元的拓撲判定分述如下。

對于六面體單元，判定其與點P的拓撲關系的流程如下：

(3)識別步驟(2)中計算的向量內積是否存在大于0的情況，若存在，則判定點P在單元外部，否則繼續執行下一步；

(6)判斷步驟(5)中計算的向量內積是否存在大于0的情況，若存在，則判定點P在單元外部；否則在單元內部。

對于四面體單元，判定其與點P的拓撲關系的流程如下：

(3)判定步驟(2)中計算的向量內積是否存在大于0的情況，若存在，則判定點P在單元外部，否則繼續執行下一步；

(6)判定步驟(5)中計算的向量內積是否存在大于0的情況，若存在，則可判定點P在單元外部，否則在單元內部。

圖3 點與六面體單元及四面體單元的拓撲關系的兩種情況

2.3 反演方法及流程圖

該優化問題的求解方法有很多，傳統的方法如牛頓下山法，雖然該方法具有收斂速度快的優點，但由于等參單元的形函數與局部坐標多為非線性關系，初始點的選取對收斂性能影響很大。在無先驗知識的前提下，牛頓下山法極易陷入局部最優，在對高階等參元進行等參逆變換時，甚至無法完成正確計算。

智能搜索方法全局收斂性好，但其計算量遠多于傳統優化算法，考慮到各鉆孔點的局部坐標只需反演一次，本文使用文獻[21]中提出的改進差分進化算法作為估計鉆孔點局部坐標的優化算法。該算法使用包含反射變異策略[22]的3種變異策略共同用于生成變異個體，結合控制參數自適應調整機制[23]，對于多模態函數優化問題具有非常好的全局收斂性能，詳情可參照文獻[21]。差分進化算法計算插值點水頭值的流程如圖4所示。

圖4 插值計算觀測點位置處水頭值的流程

3 算 例

為了測試本文方法對于平面單元和空間單元的插值效果，使用兩個滲流場模型作為測試算例。由于4節點四邊形單元和3節點三角形單元是使用最廣的平面單元，4節點四面體單元和8節點六面體單元是使用最廣的立體單元，因此，本文僅測試該方法在這些單元中的插值效果。對于高階單元，本文方法雖然同樣適用，但需根據具體的單元編號順序編寫相應的程序。

由公式(3)可知，有限元插值方法的插值精度取決于局部坐標的反演精度。將不同反演精度(1×10-3、1×10-5和1×10-8m)下的插值結果與使用Tecplot軟件計算的插值結果進行比較，以研究本文使用的有限元插值方法的插值性能；將Tecplot軟件計算的鉆孔水頭作為參照值，與本文方法的插值結果進行比較以檢驗該方法的計算有效性。

本文方法使用Fortran語言編寫，軟件環境為Visual Studio 2017 with intel Fortran 2019，在處理器為AMD Ryzen 7 3700U with Radeon Vega Mobile Gfx、內存為8 GB的電腦上運行。

3.1 二維模型

使用修改自文獻[24]中的二維模型測試本文方法的插值性能。研究區域為一個面積為36 km2、貯藏系數為0.001的承壓含水層，其中存在3個抽水井(P1、P2和P3)，抽水流量分別為6 000、3 000和800 m3/d，研究區域劃分成3個導水系數分別為200、150和50 m2/d的材料分區(區域1、2和3)，靠近模型上邊界2 000 m的區域內有0.15×10-3m/d的入流補給。模型下邊界為Dirichlet邊界(水頭為100 m)，左邊界為Neumann邊界(單寬流量為0.25 m3/(d·m))，模型上、下邊界為不透水邊界。

假設研究區域中有18個水位觀測井。為了研究本文插值算法的性能，采用330 m×330 m的單元尺寸劃分網格，使各觀測井不在單元節點上。圖5為研究區域的模型構建示意圖、有限元網格劃分以及抽水井和觀測井位置分布。平面單元下不同精度的有限元插值與Tecplot插值的插值效果對比如圖 6所示，表1反映了平面單元不同精度下的有限元插值方法的耗時情況。

圖5 研究區域模型構建示意、有限元網格劃分以及抽水井和觀測井位置分布(標注單位：m)

由圖6和表1可知，有限元插值與Tecplot插值方法計算的觀測井水頭相差較小且耗時很少，說明了本文方法的有效性；隨著觀測孔局部坐標反演精度的提高，絕對誤差變小，可以發現，采用1.00×10-5m作為終止條件時，其取得的誤差精度與1.00×10-8m非常接近；觀測井OB7～OB13位于三角形單元內，有限元插值的水頭與Tecplot插值的誤差較小，這是由于三角形單元的形函數在單元面上為平面(如圖7(a))，而矩形單元的形函數在單元面上為直邊曲面(如圖7(b))，因此，位于3節點三角形單元內的有限元插值實則為線性插值。

圖6 平面單元下不同精度的有限元插值與Tecplot插值的插值效果對比

圖7 三角形單元和矩形單元的形函數示意圖

表1 平面單元下不同精度的有限元插值方法的耗時對比

3.2 三維模型

將某水電站樞紐工程的有限元初始滲流場分析模型作為本文的三維算例。研究區域的地質構造較為復雜，其出露的基巖屬強度高、抗風化、吸水率低的堅硬巖類，為侏羅系變質砂巖和燕山期侵入正長巖，侵入巖與變質砂巖膠結良好。根據風化程度，可分為強風化層、弱風化層和新鮮基巖，其中強風化層厚度約為2～5 m，弱風化層厚度約為12～18 m。此外，區域內存在兩條較大規模的位于河床和右岸壩肩處的斷層，分別記為F59和F2。

數值模型的建立使用六面體和四面體單元進行剖分，包括15 335個節點、35 874個單元和8個單元組。剖分后的有限元模型及單元組組成如圖 8所示，模型單元組次序及使用的參數列于表2。研究域中共有20個鉆孔，其位置如圖9所示。

圖9 三維算例模型使用的20個鉆孔點位置示意圖

表2 三維算例模型的輸入參數

圖8 三維算例模型的組成及單元剖分(標注單位：m)

圖10反映了20個鉆孔點對應的線性插值結果與不同反演精度(1.00×10-3、1.00×10-5和1.00×10-8m)所對應的有限元插值結果的比較。由圖10可知，各精度下的插值結果與Tecplot軟件插值的結果相差較小，說明了本文方法的正確性。需要注意的是，ZK103、ZK111、ZK149和ZK172都處于水頭變化較大的六面體單元內，其單元形函數屬于非線性，因此與Tecplot使用的插值具有一定的偏差。

圖10 三維算例不同收斂精度的有限元插值與Tecplot插值的結果對比

表3反映了本文方法中各主要步驟的耗時情況。由表3可看出，本文插值方法總體耗時極少，各主要步驟中耗時最多的為插值點歸屬單元的確定過程，僅需87 ms。綜上所述，由于本文方法受時間的影響較小，且各插值點的等參逆變換反演過程僅需進行1次。為了使插值結果更加準確，建議采用較高的收斂精度標準(即1.00×10-8m)。

表3 三維算例不同精度的有限元插值的耗時對比

4 結 論

本文針對觀測點不在有限元網格節點上的變量(以水頭為例)插值問題進行了研究，將點與平面、空間單元的拓撲關系判別法和智能優化方法結合，提出了一種基于等參單元形函數的有限元插值函數法。通過兩個算例插值計算了插值點的水頭值，并與Tecplot軟件的計算值進行比較，結論如下：

(1)本文方法結果可靠，計算耗時短；

(2)本文方法基于有限元插值機理，非常適合于高階非線性單元的插值；

(3)本文方法可推廣至超大規模的有限元模型插值問題，其耗時量主要取決于歸屬單元的確定，本文程序使用順序查找算法，可采用高性能查找算法進一步減少耗時。