宿 曉 輝，張 明 亮，賀 英 芝

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

0 引 言

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

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

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

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

1 滲流計(jì)算方法和原理

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

(1)

1.1 控制方程離散

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

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

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

?Vcv?·FdV=0

(2)

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

(3)

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

ScvdS=0

(4)

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

(5)

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

(6)

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

(7)

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

(8)

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

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

(9)

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

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

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

(2)通過初始值H0，求得初始?xì)埐钕蛄縭0，并利用預(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，…進(jìn)行計(jì)算，直到收斂完成。

α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)用在土石壩滲流計(jì)算時(shí)，滲流場在邊界上還應(yīng)滿足一定的條件，主要包括以下幾種[22-23]。

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

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

(19)

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

(20)

浸潤面Γ3：

(21)

溢出邊界面Γ4：

(22)

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

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

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

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

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

(23)

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

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

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

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

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

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

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

1.4 GPU并行加速

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

2 算例驗(yàn)證

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

假設(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，計(jì)算壩體的滲流場分布，計(jì)算結(jié)果如圖4～5所示，并與GEO-Studio(2018學(xué)習(xí)版)計(jì)算結(jié)果進(jìn)行了對比。

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

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

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

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

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

圖6 不同計(jì)算方法的時(shí)間對比匯總

3 工程應(yīng)用

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

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