基于有限體積法的三維穩(wěn)態(tài)滲流計算方法及其應(yīng)用

宿 曉 輝，張 明 亮，賀 英 芝

(大連理工大學(xué) 建設(shè)工程學(xué)部水利工程學(xué)院，遼寧 大連 116033)

0 引 言

滲流分析和滲流控制是大壩工程中極為重要的一項研究內(nèi)容，直接關(guān)系到工程的費用和安全。國內(nèi)外許多破壞或潰決的工程多數(shù)都是由于滲流控制不當所引起的。據(jù)不完全統(tǒng)計，中國擁有水庫大壩近10萬座，其中95%以上為土石壩。這些水庫大多興建于20世紀50～70年代，并且限于當時的技術(shù)和經(jīng)濟條件，許多工程均是邊勘測、邊設(shè)計、邊施工，這給大壩的安全運行埋下了巨大的隱患。在土石壩失事事故的原因統(tǒng)計中，約有1/2是由于滲流問題所引起的，所以對病險壩進行滲流計算就顯得尤為重要[1-4]。

隨著計算機技術(shù)的高速發(fā)展，數(shù)值解法在求解滲流問題領(lǐng)域逐漸占據(jù)主導(dǎo)地位。常用的數(shù)值求解方法主要是有限差分法[5-6]和有限單元法[7-10]，其他方法，如邊界元法[11]、無單元法[12]等也有應(yīng)用。目前，在這些數(shù)值計算方法當中，較為流行的是有限單元法的應(yīng)用[8]，但是其計算工作量大，應(yīng)用到求解大型工程滲流計算時，耗時較長。有限體積法與有限元相比，計算工作量小，并且其控制體積內(nèi)局部絕對守恒的特性對于處理流動問題具有更好的性能[13-14]。

近年來GPU并行加速技術(shù)發(fā)展迅速且已相對成熟，浮點計算能力是同代CPU的10倍以上。NVIDIA等公司發(fā)布的較為完善的CUDA和OpenAcc并行編程標準打破了GPU技術(shù)應(yīng)用的瓶頸[15-16]，使得GPU技術(shù)應(yīng)用到多體系統(tǒng)模擬、計算流體力學(xué)、核科學(xué)等各個計算領(lǐng)域。在水利計算中，GPU并行優(yōu)化技術(shù)已經(jīng)應(yīng)用到水動力學(xué)模型[17-19]、分布式水文模型[20]、梯級水庫調(diào)度[16]的求解中，取得了較好效果，而在土石壩滲流計算領(lǐng)域中，目前尚未有研究明確實現(xiàn)GPU并行加速。

基于此，本文提出了基于非結(jié)構(gòu)化網(wǎng)格的有限體積法滲流模型。模型采用預(yù)處理共軛梯度算法快速求解代數(shù)方程組，使用固定網(wǎng)格變滲透系數(shù)法求解計算域內(nèi)含有浸潤面的問題，創(chuàng)建無系數(shù)矩陣數(shù)值迭代求解方法對滲流控制方程進行離散，應(yīng)用OpenACC-GPU并行技術(shù)加速求解。通過與商用軟件GEO-Studio進行結(jié)果對比，證明本文提出的模型滿足工程滲流實際需要。

1 滲流計算方法和原理

非均質(zhì)各向異性穩(wěn)態(tài)飽和滲流問題的控制方程為

(1)

1.1 控制方程離散

采用非結(jié)構(gòu)化網(wǎng)格的有限體積法對滲流控制方程進行離散，采用預(yù)處理共軛梯度法(PCG)加速對離散方程求解。采用Cell-Vertex方法形成控制體，并在控制體上對公式(1)進行積分。該方法中四面體單元頂點的控制體如圖1所示。對于每一個頂點的控制體，用四面體的中值對偶來實現(xiàn)；四面體單元的頂點分別為A，P，B，C，點O是四面體單元APBC的中心，點a，b，c分別是AP，BP，CP的中點；1，2，3分別為面APC，CBP和ABP的形心；在Cell-Vertex模式中，所有的計算變量值均存儲在頂點A，P，B，C上；三角形O1a，O3a，O3b，O2b，O1c和O2c組成了四面體單元中節(jié)點P的控制體表面；同樣的，構(gòu)建了剩余節(jié)點A，B，C的控制體表面[21]?；谶@種思想構(gòu)造的控制體會在一定程度上降低網(wǎng)格奇異性的影響，提高計算的準確度。

圖1 圍繞節(jié)點P的控制體

對滲流控制方程在上述控制體上進行積分(以P為例)：

?Vcv?·FdV=0

(2)

采用高斯散度定理對公式(2)積分：

(3)

式中：ncell代表與控制節(jié)點P相關(guān)聯(lián)的單元數(shù)目，ΔScj表示控制體在四面體單元j中的對應(yīng)邊界面。對于給定的四面體單元，可以采用閉合曲面積分公式計算ΔScj：

ScvdS=0

(4)

在控制體所有邊界面法向量積分為0，可以得到：

(5)

將其代入式(3)可以得到：

(6)

四面體單元的中心梯度是通過單元頂點變量值采用格林公式計算得來：

(7)

四面體單元的中心滲透系數(shù)是通過單元頂點滲透系數(shù)值采用體積加權(quán)法計算得來：

(8)

式中：Hi為四面體單元j中頂點i處的變量值；Si為四面體單元j中頂點i所對應(yīng)的面的法向量，其大小為頂點i所對應(yīng)面的面積；V為四面體單元體積；(Ki)m為四面體單元j中頂點m處的滲透系數(shù)值。

由式(6)推導(dǎo)可得：

(9)

式中：R(H)為殘余誤差，當R(H)<ε(ε為收斂控制精度)時，認為計算結(jié)果已經(jīng)收斂，達到所需精度要求。公式(9)同時表明，本文提出的穩(wěn)態(tài)滲流計算方法無需建立二維系數(shù)矩陣A，而將AH整合為一列向量，節(jié)省了存儲空間，減少了內(nèi)存占用。

本文采用預(yù)處理共軛梯度法快速求解代數(shù)方程組，其求解步驟為：

(1)給定求解域初始值H0。

(2)通過初始值H0，求得初始殘差向量r0，并利用預(yù)處理對角陣M-1得到P0。

r0=b-AH0=0-AH0=-R(H0)

(10)

z0=M-1r0

(11)

p0=z0

(12)

式中：M=diag(a11，a22，...，ann)。

(3)對k=0，1，…進行計算，直到收斂完成。

αk=(rk，zk)/(Apk，pk)

(13)

Hk+1=Hk+αkpk

(14)

rk+1=rk-αkApk

(15)

zk+1=M-1rk+1

(16)

βk=(zk+1，zk+1)/(rk，zk)

(17)

pk+1=zk+1+βkpk

(18)

1.2 邊界條件

滲流模型應(yīng)用在土石壩滲流計算時，滲流場在邊界上還應(yīng)滿足一定的條件，主要包括以下幾種[22-23]。

第一類邊界Γ1，即已知常水頭邊界：

H(x，y，z)=Hfixed(x，y，z)

(19)

第二類邊界Γ2，即流量邊界條件q0=qfixed，在滲流計算分析中設(shè)定q0=0得到不透水邊界：

(20)

浸潤面Γ3：

(21)

溢出邊界面Γ4：

(22)

溢出邊界面的位置采用迭代方法獲得[22，24]，其計算基本過程如下：

(1)假設(shè)溢出面范圍?？梢栽O(shè)定求解區(qū)域下游面上的下游水位以上區(qū)域全部為滲流溢出面。

(2)將假設(shè)的溢出面給定第一類邊界條件。將溢出邊界面上節(jié)點的總水頭給為位置高程，即H(x，y，z)=y。

(3)通過求解滲流方程得到計算域內(nèi)的總水頭分布。

(4)溢出面修正。在計算收斂之后，根據(jù)溢出邊界面的另一個控制條件即式(23)進行更正：

(23)

(5)將不滿足式(22)的溢出面節(jié)點從溢出邊界面剔除，并將下游面內(nèi)溢出面以上H(x，y，z)>y的節(jié)點歸入溢出面范圍。

(6)根據(jù)修正后的溢出面返回步驟(2)重新進行迭代計算，直到兩次迭代計算溢出面節(jié)點數(shù)之差小于某一給定值，則迭代結(jié)束，完成計算。

1.3 固定網(wǎng)格變節(jié)點滲透系數(shù)

本文使用固定網(wǎng)格變節(jié)點滲透系數(shù)方法計算域內(nèi)浸潤線存在問題。將浸潤線以上節(jié)點滲透系數(shù)設(shè)置為近似不透水情況時，以滿足式(21)所設(shè)定的控制條件，其主要步驟如下[25]：

(1)假定全區(qū)域均為飽和區(qū)即H(x，y，z)≥y，各節(jié)點滲透系數(shù)均為飽和滲透系數(shù)，進行初次滲流迭代計算。

(2)初次滲流計算完成后，通過比較節(jié)點水頭H(x，y，z)與節(jié)點坐標y值的大小，將區(qū)域分為飽和區(qū)與無壓區(qū)，無壓區(qū)即浸潤面以上的區(qū)域，將無壓區(qū)所有節(jié)點的滲透系數(shù)降低1 000倍，進入下次迭代計算。

(3)通過不斷更新無壓區(qū)并改變其滲透系數(shù)大小，尋找浸潤線位置，當無壓區(qū)節(jié)點數(shù)量小于控制收斂精度時，得到滿足收斂精度下的穩(wěn)態(tài)飽和滲流計算結(jié)果。

1.4 GPU并行加速

OpenACC是目前主流的GPU并行加速技術(shù)之一，是一種以編譯器導(dǎo)語形式實現(xiàn)的GPU并行標準，其優(yōu)勢為減少底層代碼的修改實現(xiàn)加速手段。本文通過支持OpenACC加速標準的PGI編譯器，采用kernel構(gòu)件實現(xiàn)本文模型的并行加速計算功能。

2 算例驗證

本文采用下面二維均質(zhì)土壩滲流問題進行模型評估與驗證。

假設(shè)一二維均質(zhì)黏土壩[23]，壩底相對高程0 m，壩高50.0 m，頂寬10 m，上游壩坡坡比1∶2，下游壩坡坡比1∶1.5，飽和滲透系數(shù)為3.4×10-8m/s，上游水位45 m，下游無水，底面為不透水邊界，大壩剖面圖見圖2，其網(wǎng)格剖分見圖3，計算壩體的滲流場分布，計算結(jié)果如圖4～5所示，并與GEO-Studio(2018學(xué)習(xí)版)計算結(jié)果進行了對比。

圖2 均質(zhì)黏土壩剖面圖(尺寸單位：m)

圖3 驗證模型網(wǎng)格(節(jié)點9 148個；單元17 795個)

圖4 驗證算例壓力水頭等值線(單位：m)

圖5 驗證算例總水頭等值線(單位：m)

對比本文模型與GEO-Studio(2018學(xué)習(xí)版)結(jié)果，發(fā)現(xiàn)兩者計算結(jié)果誤差在±2.5%內(nèi)，壓力水頭與總水頭等值線的分布吻合，浸潤線的位置接近，能滿足工程計算需要。圖6匯總了本文采用的各種計算方法及GEO-Studio的計算時間，通過對比時間結(jié)果發(fā)現(xiàn)，當采用預(yù)處理共軛梯度法(PCG)時，比只應(yīng)用共軛梯度法(CG)速度提升了3.1倍，因為評估算例為均質(zhì)壩，在一定程度上削減了預(yù)處理共軛梯度法的優(yōu)勢，當將其應(yīng)用在非均質(zhì)壩時，相信預(yù)處理共軛梯度法對于共軛梯度法的效率提升將會更加明顯。當OpenACC與預(yù)處理共軛梯度法聯(lián)合運用時，加速效果極為顯著，加速比達到70，計算時間也遠小于網(wǎng)格數(shù)量較多的GEO-Studio模型計算所需時間。

圖6 不同計算方法的時間對比匯總

3 工程應(yīng)用

甘肅省代古寺水庫混凝土面板堆石壩位于大“S”形河道內(nèi)，河谷狹窄，呈深切“V”字形，為橫向或斜向河谷，兩岸山體雄厚，岸坡高陡直立。河床底高程為1 667.00 m，正常蓄水位為1 804.00 m，下游尾水位為1 675.00 m，壩頂高程為1 809.00 m。表1為三維滲流計算的各材料分區(qū)滲透系數(shù)。圖7為本模型采用的有限體積網(wǎng)格(節(jié)點725 036個；單元4 040 769個)，圖8～9分別為大壩最大典型斷面(壩0+115)的壓力水頭等值線與總水頭等值線分布結(jié)果。由圖可知，本文模擬結(jié)果與商業(yè)軟件GEO-Studio結(jié)果具有較高一致性，驗證了本文計算方法的正確性和可靠性。

表1 面板堆石壩各材料分區(qū)滲透系數(shù)

圖7 應(yīng)用算例有限體積網(wǎng)格

圖8 工程應(yīng)用算例最大典型斷面(壩0+115)壓力水頭等值線(單位：m)

圖9 工程應(yīng)用算例最大典型斷面(壩0+115)總水頭等值線(單位：m)

4 結(jié) 語

本文提出了一種基于非結(jié)構(gòu)化網(wǎng)格控制體積法滲流計算模型。模型采用預(yù)處理共軛梯度算法求解代數(shù)方程組，使用固定網(wǎng)格變滲透系數(shù)法求解域內(nèi)含浸潤面問題，發(fā)展了無系數(shù)矩陣數(shù)值迭代求解方法。驗證算例與工程實際應(yīng)用表明，本文所提出的土石壩三維穩(wěn)態(tài)滲流計算方法有效可行，滿足工程的實際需要。本文采用的OpenACC-GPU并行策略具有較高的加速比，為提高大型工程滲流計算效率提供了新的手段。