一種用于陰極保護(hù)數(shù)值模擬的高性能計(jì)算方法

(1. 中海石油(中國(guó))有限公司 深圳分公司 深水工程建設(shè)中心，深圳 518054; 2. 大連理工大學(xué) 船舶工程學(xué)院，大連 116024)

陰極保護(hù)(Catholic Protection，簡(jiǎn)稱CP)是通過電化學(xué)方法降低金屬腐蝕速率的一種技術(shù)[1]，有犧牲陽(yáng)極陰極保護(hù)(Galvanic Anode Catholic Protection，簡(jiǎn)稱GACP)和外加電流陰極保護(hù)(Impressed Current Catholic Protection，簡(jiǎn)稱ICCP)。犧牲陽(yáng)極陰極保護(hù)設(shè)計(jì)[2]，大多根據(jù)經(jīng)驗(yàn)公式確定需要的陽(yáng)極數(shù)量，然后再根據(jù)設(shè)計(jì)人員的經(jīng)驗(yàn)確定陽(yáng)極布置方案，因此該方法的誤差較大。采用該方法時(shí)一般選取較大的安全系數(shù)，但這會(huì)導(dǎo)致陽(yáng)極用量增加和經(jīng)濟(jì)成本升高，而且增加了安裝施工難度和降低結(jié)構(gòu)強(qiáng)度。外加電流陰極保護(hù)系統(tǒng)設(shè)計(jì)[3]，采用的輔助陽(yáng)極的單體輸出電流較大，輔助陽(yáng)極的數(shù)量較少，結(jié)構(gòu)表面的保護(hù)電位分布范圍較大，且基于規(guī)范經(jīng)驗(yàn)公式的設(shè)計(jì)方法，往往會(huì)導(dǎo)致結(jié)構(gòu)表面發(fā)生欠保護(hù)或過保護(hù)現(xiàn)象，從而降低結(jié)構(gòu)強(qiáng)度。

為了彌補(bǔ)經(jīng)驗(yàn)公式的不足，科研人員將有限差分法(Finite Difference Method，簡(jiǎn)稱FDM)[4-5]、邊界元法(Boundary Element Method，簡(jiǎn)稱BEM)[6-8]和有限元法(Finite Element Method，簡(jiǎn)稱FEM)[9-11]等數(shù)值模擬方法應(yīng)用到陰極保護(hù)設(shè)計(jì)中。在犧牲陽(yáng)極陰極保護(hù)設(shè)計(jì)中，通過數(shù)值模擬可以確定最優(yōu)的陽(yáng)極布置方案，即采用最少的陽(yáng)極達(dá)到最佳的保護(hù)效果，從而降低經(jīng)濟(jì)成本和施工難度。在外加電流陰極保護(hù)設(shè)計(jì)中，也可通過數(shù)值模擬確定最優(yōu)的輔助陽(yáng)極數(shù)量、布置和電流輸出方案，避免過保護(hù)問題。

英國(guó)Beasy公司針對(duì)陰極保護(hù)數(shù)值模擬，已開發(fā)出比較成熟的基于邊界元法的商業(yè)軟件。目前在國(guó)內(nèi)，基于邊界元法的陰極保護(hù)數(shù)值模擬也有了較多的研究[12-15]，但還沒有成熟的商業(yè)軟件。采用邊界元法進(jìn)行陰極保護(hù)數(shù)值模擬，無需考慮域內(nèi)電位分布，只需要對(duì)邊界進(jìn)行離散，且可以直接得到結(jié)構(gòu)外表面的保護(hù)電位分布，因此，邊界元法是目前最理想的陰極保護(hù)數(shù)值模擬方法。

然而，隨著海洋結(jié)構(gòu)物的結(jié)構(gòu)形式日益復(fù)雜，尺寸越來越大，數(shù)值模擬的節(jié)點(diǎn)規(guī)模也越來越大，大規(guī)模模型的數(shù)值模擬計(jì)算效率很低，嚴(yán)重影響了模擬的進(jìn)度和效率。當(dāng)前典型海洋工程結(jié)構(gòu)物邊界元模型的節(jié)點(diǎn)數(shù)通常在104～105數(shù)量級(jí)，對(duì)應(yīng)的雙精度型浮點(diǎn)系數(shù)矩陣將占用計(jì)算機(jī)0.7～74.5 GB內(nèi)存空間，對(duì)計(jì)算機(jī)硬件資源耗用巨大。通過常規(guī)串行計(jì)算方式生成如此規(guī)模的系數(shù)矩陣，計(jì)算效率很低。邊界元數(shù)值模擬中，系數(shù)矩陣中每行數(shù)據(jù)的計(jì)算相互獨(dú)立且沒有交叉調(diào)用，易于實(shí)現(xiàn)并行計(jì)算，理論上通過并行算法提高的計(jì)算效率倍數(shù)為計(jì)算機(jī)的CPU數(shù)，因此在CPU數(shù)目較多的服務(wù)器上通過并行算法可以大幅提升系數(shù)矩陣的生成效率。

C#是一種面向?qū)ο蟮母呒?jí)程序語(yǔ)言，從C和C++語(yǔ)言演化而來[16]，吸收了C、C++與Java語(yǔ)言的優(yōu)點(diǎn)，具有Visual Basic編譯語(yǔ)言的高效率和C++編譯語(yǔ)言的強(qiáng)大功能。C#語(yǔ)言開發(fā)的系統(tǒng)具有界面友好、易維護(hù)和升級(jí)等優(yōu)點(diǎn)，并能有效地保護(hù)算法和數(shù)據(jù)。因此，C#語(yǔ)言是開發(fā)Windows數(shù)值模擬應(yīng)用程序的較優(yōu)選擇。但是C#語(yǔ)言是基于Microsoft .NET Framework(簡(jiǎn)稱.Net)運(yùn)行環(huán)境的，在.Net環(huán)境下，一維數(shù)組和多維數(shù)組作為單個(gè)托管對(duì)象(由.Net自動(dòng)分配與釋放內(nèi)存),其存儲(chǔ)能力都不超過2 GB，這限制了數(shù)值模擬的計(jì)算規(guī)模(使用雙精度浮點(diǎn)數(shù)，節(jié)點(diǎn)數(shù)或矩陣的階數(shù)不得超過16 384)，因此需要研究.Net運(yùn)行環(huán)境下大規(guī)模矩陣在內(nèi)存中的存儲(chǔ)方法，以提高數(shù)值模擬的計(jì)算規(guī)模。

對(duì)于線性方程組，見式(1)，不同求解方法的計(jì)算效率會(huì)相差幾倍、幾十倍，甚至更高，在計(jì)算大規(guī)模問題時(shí)，計(jì)算效率會(huì)嚴(yán)重影響模擬的進(jìn)度和效率。因此需要對(duì)大規(guī)模線性方程組的快速求解方法進(jìn)行研究，以提高陰極保護(hù)數(shù)值模擬的計(jì)算效率。

Ax=b

(1)

本工作針對(duì)邊界元陰極保護(hù)數(shù)值模擬，使用C#語(yǔ)言，基于.Net運(yùn)行環(huán)境和Lapack函數(shù)庫(kù)，以實(shí)現(xiàn)數(shù)值模擬中系數(shù)矩陣并行生成、大規(guī)模矩陣存儲(chǔ)和大規(guī)模線性方程組快速求解，對(duì)實(shí)際工程應(yīng)用和開發(fā)我國(guó)自主的陰極保護(hù)數(shù)值模擬商業(yè)軟件具有一定的指導(dǎo)作用。

1 陰極保護(hù)數(shù)值模擬的基本原理

1.1 控制方程

在電解質(zhì)中，陰極保護(hù)相關(guān)電場(chǎng)中的電流密度和電位滿足歐姆定律[11]，如式(2)所示。

(2)

式中：q為電流密度矢量;ρ為海水電阻率;φ為電位。

取電解質(zhì)中一立方微元體作為研究對(duì)象，則垂直于立方微元體六個(gè)面的電流密度如圖1所示。

圖1 三維微元體Fig. 1 Cubic volume element

任意時(shí)刻，圖1中的立方微元體中的電流變化量為

(3)

當(dāng)陰極保護(hù)系統(tǒng)相關(guān)電場(chǎng)達(dá)到平衡狀態(tài)且無內(nèi)部場(chǎng)源時(shí)，由電荷守恒可知，微元體中的電流通量為0。因此通過式(3)可以得到電場(chǎng)域內(nèi)控制方程即Laplace方程

(4)

1.2 邊界條件

陰極保護(hù)邊界元數(shù)值模擬的兩個(gè)物理量分別是結(jié)構(gòu)表面的保護(hù)電位和保護(hù)電流密度。陰極保護(hù)問題的邊界條件一般有兩類：第二類邊界條件和第三類邊界條件。

(1) 第二類邊界條件，表面電流密度為0，例如：海水與空氣的交界面和結(jié)構(gòu)防腐蝕涂層完好部分。

(5)

式中：n為邊界處的法線方向。

(2) 第三類邊界條件，防腐蝕涂層破損部位，裸鋼和陽(yáng)極表面，電位和電流密度的關(guān)系滿足極化曲線。

(6)

1.3 基于邊界元的陰極保護(hù)數(shù)值模擬

對(duì)于Laplace方程，其離散形式的邊界積分方程為[17]

(7)

式中：ci為節(jié)點(diǎn)處與邊界光滑度相關(guān)的系數(shù);N為形函數(shù)；m為單元數(shù)；φi為節(jié)點(diǎn)i處的電位；Γi為第j個(gè)單元；φi*為基本解；qi*為基本解的法向?qū)?shù)。

(8)

(9)

對(duì)每個(gè)節(jié)點(diǎn)，按式(7)求積分，則可得線性方程組

(10)

其中

C=diag(c1,c2,…,cn)

(11)

(12)

(13)

令

(14)

則式(11)可以化成

Hφ=Gq

(15)

帶入1.2中的邊界條件，將式(15)變換成線性方程組

Aφ=b

(16)

由于邊界條件一般是非線性的，需要反復(fù)迭代求解方程組(16)，得到各節(jié)點(diǎn)處的保護(hù)電位。

對(duì)于節(jié)點(diǎn)數(shù)為N的模型，在求解過程中，需要首先生成H和G矩陣，生成H和G的時(shí)間復(fù)雜度為O(N2)。通常情況下，對(duì)于節(jié)點(diǎn)數(shù)N=10 000的模型，生成系數(shù)矩陣需要10 min左右；而對(duì)于節(jié)點(diǎn)數(shù)N=100 000的模型，生成系數(shù)矩陣則需要大概17 h。A是N×N的滿陣，在求解線性方程組(14)時(shí)，時(shí)間復(fù)雜度為O(N3)，空間復(fù)雜度為O(N2)，使用雙精度浮點(diǎn)數(shù)(一個(gè)雙精度浮點(diǎn)數(shù)占8個(gè)字節(jié)的內(nèi)存)存儲(chǔ)矩陣時(shí)，占用的內(nèi)存M可通過下式計(jì)算。

(17)

2 并行計(jì)算方法

本工作中單機(jī)并行計(jì)算技術(shù)的實(shí)現(xiàn)基于微軟.NET 4.0平臺(tái)，利用任務(wù)并行庫(kù)TPL[18](Task Parallel Library)改寫計(jì)算程序。TPL采用Task類(任務(wù)類)，它的聲明對(duì)象落實(shí)到了具體的計(jì)算任務(wù)，可以自動(dòng)實(shí)現(xiàn)線程的創(chuàng)建和執(zhí)行。

針對(duì)陰極保護(hù)數(shù)值模擬，從式(7)可知每個(gè)節(jié)點(diǎn)處的積分計(jì)算量基本相同，因此可以依據(jù)節(jié)點(diǎn)數(shù)和線程數(shù)劃分任務(wù)。本研究通過TPL中的Parallel.For函數(shù)開發(fā)了邊界元系數(shù)矩陣并行生成算法。

對(duì)系數(shù)矩陣并行算法進(jìn)行了測(cè)試，測(cè)試環(huán)境如表1所示，測(cè)試結(jié)果如圖2所示。從圖2可以看出，當(dāng)線程數(shù)小于計(jì)算機(jī)處理器數(shù)目時(shí)，增加線程數(shù)可以提升計(jì)算效率。在線程數(shù)等于處理器數(shù)目時(shí)，計(jì)算效率最高，進(jìn)一步增加線程數(shù)，會(huì)降低計(jì)算效率。

表1 測(cè)試環(huán)境Tab. 1 Testing enviroment

圖2 計(jì)算時(shí)間與線程數(shù)的關(guān)系Fig. 2 Relationship between computing time and number of threads

在研究中發(fā)現(xiàn),計(jì)算機(jī)同時(shí)進(jìn)行計(jì)算的線程數(shù)不會(huì)超過處理器核數(shù)，當(dāng)使用的線程數(shù)超過處理器核數(shù)時(shí)，超出的線程在前面的線程執(zhí)行完畢后才會(huì)執(zhí)行。由于每個(gè)線程的計(jì)算量是一樣的，故在不考慮其他因素的理想情況下，本例中并行計(jì)算時(shí)間為

(18)

式中：t為計(jì)算時(shí)間;C為由總?cè)蝿?wù)量決定的常數(shù)。

公式(18)說明，當(dāng)NCPU等于32時(shí)，計(jì)算效率最高，繼續(xù)增加NCPU會(huì)降低計(jì)算效率。

需要注意的是，公式(18)是在理想狀態(tài)下給出的，不能用于求解并行計(jì)算時(shí)間，僅用來定性分析計(jì)算效率。

另一方面，使用過多的線程數(shù)會(huì)占用過多的計(jì)算機(jī)資源，申請(qǐng)和釋放線程也會(huì)消耗一定的時(shí)間，從而也會(huì)降低計(jì)算效率。

從測(cè)試結(jié)果和以上分析可知，在并行計(jì)算的線程數(shù)等于計(jì)算機(jī)處理器核數(shù)時(shí)，并行計(jì)算的效率最高。

3 大規(guī)模矩陣存儲(chǔ)方法

在.Net運(yùn)行環(huán)境下，一維數(shù)組和多維數(shù)組作為單個(gè)托管對(duì)象(由.Net自動(dòng)分配與釋放內(nèi)存)的存儲(chǔ)能力都不超過2 GB，矩陣階數(shù)不得超過16 384，這限制了數(shù)值模擬的計(jì)算規(guī)模。本工作通過使用非托管內(nèi)存實(shí)現(xiàn)了大規(guī)模矩陣的存儲(chǔ)。

在C#語(yǔ)言中，Marshal.AllocHGlobal(IntPtr cb)函數(shù)可以申請(qǐng)連續(xù)的非托管內(nèi)存(通過代碼申請(qǐng)與釋放的內(nèi)存)，IntPtr為指針類型的數(shù)據(jù)，可以表示內(nèi)存中的所有地址，即AllocHGlobal函數(shù)可以申請(qǐng)計(jì)算機(jī)空閑的最大連續(xù)內(nèi)存空間。內(nèi)存申請(qǐng)完畢后，可以通過非安全指針進(jìn)行訪問。

AllocHGlobal申請(qǐng)的內(nèi)存空間是連續(xù)的，將矩陣存儲(chǔ)在其中后，可以直接使用Lapack函數(shù)庫(kù)的矩陣運(yùn)算函數(shù)進(jìn)行計(jì)算，因此在.Net運(yùn)行環(huán)境下，使用非托管內(nèi)存進(jìn)行大規(guī)模矩陣的存儲(chǔ)是最理想的選擇。

4 大規(guī)模線性方程組快速求解方法

目前，求解線性方程組主要借助Lapack函數(shù)庫(kù)，Intel和AMD公司分別開發(fā)了針對(duì)自己處理器的Lapack函數(shù)庫(kù)：MKL(Math Kernel Library)和ACML(AMD Core Math Library)。使用其中的矩陣運(yùn)算函數(shù)，比自行編寫的矩陣運(yùn)算函數(shù)效率高很多。

求解線性方程組(16)，一般使用Lapack函數(shù)庫(kù)中dgetrf和dgetrs兩個(gè)函數(shù)，其中dgetrf函數(shù)對(duì)矩陣A進(jìn)行三角分解，dgetrs利用分解后的矩陣進(jìn)行求解。

目前Lapack函數(shù)庫(kù)都封裝在DLL(Dynamic Link Library)中，DLL分為兩種：一種是直接由機(jī)器語(yǔ)言組成(由C、C++和Fortran等語(yǔ)言編譯的DLL)；另一種是由.Net中間語(yǔ)言組成，在被調(diào)用時(shí)，再由.Net翻譯為機(jī)器語(yǔ)言。

目前MKL和ACML提供的Lapack函數(shù)庫(kù)都是采用Fortran開發(fā)的，因此其中的函數(shù)由機(jī)器語(yǔ)言組成。使用C、C++和Fortran語(yǔ)言都可以通過Lapack函數(shù)庫(kù)對(duì)應(yīng)的頭文件直接調(diào)用其中的函數(shù)。使用C#語(yǔ)言可以直接調(diào)用基于.Net的動(dòng)態(tài)鏈接庫(kù)(DLL文件)中的函數(shù)，但不能夠直接調(diào)用普通DLL的函數(shù)。在C#語(yǔ)言中，可以通過DllImport[19]調(diào)用基于機(jī)器語(yǔ)言的動(dòng)態(tài)鏈接庫(kù)函數(shù)。

采用表1中的測(cè)試環(huán)境，對(duì)MKL和ACML中dgetrf和dgetrs函數(shù)的求解性能進(jìn)行測(cè)試，測(cè)試結(jié)果如圖3所示。從圖3中可以看出：在本測(cè)試環(huán)境下，ACML的求解效率是MKL的兩倍左右，因此在實(shí)際使用中應(yīng)根據(jù)處理器類型采用合適的函數(shù)庫(kù)。

圖3 矩陣求解時(shí)間與節(jié)點(diǎn)數(shù)的關(guān)系Fig. 3 Matrix solving time vs number of nodes

5 求解案例

采用兩個(gè)實(shí)際工程案例對(duì)求解性能進(jìn)行了分析，測(cè)試環(huán)境如表1所示，并行計(jì)算時(shí)使用了32個(gè)線程。

5.1 張力腿平臺(tái)陰極保護(hù)數(shù)值模擬

模型包括張力腿平臺(tái)的水下主體部分、鉆井立管、張力筋腱和犧牲陽(yáng)極，共有63 024個(gè)節(jié)點(diǎn)，115 336個(gè)三角形單元，如圖4所示。求解方程組時(shí)，共迭代4次，模擬結(jié)果如圖5所示，求解耗時(shí)如表2所示。從表2可以看出：采用并行生成系數(shù)矩陣和AMD Lapack進(jìn)行求解，可以大幅提升計(jì)算效率。

圖4 張力腿平臺(tái)邊界元網(wǎng)格模型Fig. 4 Boundary element model of tension-leg platform

圖5 張力腿平臺(tái)犧牲陽(yáng)極陰極保護(hù)的模擬結(jié)果Fig. 5 Simulated result of sacrificial anode cathodic protection for tension-leg platform

計(jì)算方案生成系數(shù)矩陣時(shí)長(zhǎng)方程組求解時(shí)長(zhǎng)總時(shí)長(zhǎng)串行生成&MKL求解566.2113.6679.8并行生成&AMD求解18.057.475.4

5.2 半潛平臺(tái)陰極保護(hù)數(shù)值模擬

模型包括半潛平臺(tái)的水下主體部分、鉆井立管和張力筋腱，共有56 594個(gè)節(jié)點(diǎn)，111 204個(gè)三角形單元，如圖6所示。求解方程組時(shí)，共迭代3次，模擬結(jié)果如圖7所示，求解耗時(shí)如表3所示。從表3可以看出：采用并行生成系數(shù)矩陣和AMD Lapack進(jìn)行求解，可以大幅提升計(jì)算效率。

圖6 半潛平臺(tái)邊界元網(wǎng)格模型Fig. 6 Boundary element model of semi-submersible platform

圖7 半潛平臺(tái)外加電流陰極保護(hù)的模擬結(jié)果Fig. 7 Simulated result of impressed current cathodic protection for semi-submersible platform

計(jì)算方案生成系數(shù)矩陣時(shí)長(zhǎng)方程組求解時(shí)長(zhǎng)總時(shí)長(zhǎng)串行生成&MKL求解456.961.7518.6并行生成&AMD求解14.531.345.8

6 結(jié)論

本工作使用C#編程語(yǔ)言，在.Net運(yùn)行環(huán)境下，針對(duì)陰極保護(hù)邊界元數(shù)值模擬，從系數(shù)矩陣生成、存儲(chǔ)和大規(guī)模線性方程組求解三個(gè)方面對(duì)高性能計(jì)算方法進(jìn)行了研究，并得到了以下結(jié)論。

(1) 使用單機(jī)并行計(jì)算方法可以顯著提升邊界元數(shù)值模擬中系數(shù)矩陣的生成效率，線程數(shù)等于處理器數(shù)目時(shí)，計(jì)算效率最高；

(2) 在.Net環(huán)境下，使用非托管內(nèi)存可以使用計(jì)算機(jī)空閑的最大連續(xù)內(nèi)存空間，是數(shù)值模擬中存儲(chǔ)大規(guī)模矩陣的理想選擇；

(3) 在.Net環(huán)境下，可以使用Lapack函數(shù)庫(kù)中的高性能運(yùn)算函數(shù)對(duì)線性方程組進(jìn)行求解，不同處理器應(yīng)采用相應(yīng)的函數(shù)庫(kù)。