一種用于陰極保護數值模擬的高性能計算方法

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

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

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

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

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

C#是一種面向對象的高級程序語言，從C和C++語言演化而來[16]，吸收了C、C++與Java語言的優點，具有Visual Basic編譯語言的高效率和C++編譯語言的強大功能。C#語言開發的系統具有界面友好、易維護和升級等優點，并能有效地保護算法和數據。因此，C#語言是開發Windows數值模擬應用程序的較優選擇。但是C#語言是基于Microsoft .NET Framework(簡稱.Net)運行環境的，在.Net環境下，一維數組和多維數組作為單個托管對象(由.Net自動分配與釋放內存),其存儲能力都不超過2 GB，這限制了數值模擬的計算規模(使用雙精度浮點數，節點數或矩陣的階數不得超過16 384)，因此需要研究.Net運行環境下大規模矩陣在內存中的存儲方法，以提高數值模擬的計算規模。

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

Ax=b

(1)

本工作針對邊界元陰極保護數值模擬，使用C#語言，基于.Net運行環境和Lapack函數庫，以實現數值模擬中系數矩陣并行生成、大規模矩陣存儲和大規模線性方程組快速求解，對實際工程應用和開發我國自主的陰極保護數值模擬商業軟件具有一定的指導作用。

1 陰極保護數值模擬的基本原理

1.1 控制方程

在電解質中，陰極保護相關電場中的電流密度和電位滿足歐姆定律[11]，如式(2)所示。

(2)

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

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

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

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

(3)

當陰極保護系統相關電場達到平衡狀態且無內部場源時，由電荷守恒可知，微元體中的電流通量為0。因此通過式(3)可以得到電場域內控制方程即Laplace方程

(4)

1.2 邊界條件

陰極保護邊界元數值模擬的兩個物理量分別是結構表面的保護電位和保護電流密度。陰極保護問題的邊界條件一般有兩類：第二類邊界條件和第三類邊界條件。

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

(5)

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

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

(6)

1.3 基于邊界元的陰極保護數值模擬

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

(7)

式中：ci為節點處與邊界光滑度相關的系數;N為形函數；m為單元數；φi為節點i處的電位；Γi為第j個單元；φi*為基本解；qi*為基本解的法向導數。

(8)

(9)

對每個節點，按式(7)求積分，則可得線性方程組

(10)

其中

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

(11)

(12)

(13)

令

(14)

則式(11)可以化成

Hφ=Gq

(15)

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

Aφ=b

(16)

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

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

(17)

2 并行計算方法

本工作中單機并行計算技術的實現基于微軟.NET 4.0平臺，利用任務并行庫TPL[18](Task Parallel Library)改寫計算程序。TPL采用Task類(任務類)，它的聲明對象落實到了具體的計算任務，可以自動實現線程的創建和執行。

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

對系數矩陣并行算法進行了測試，測試環境如表1所示，測試結果如圖2所示。從圖2可以看出，當線程數小于計算機處理器數目時，增加線程數可以提升計算效率。在線程數等于處理器數目時，計算效率最高，進一步增加線程數，會降低計算效率。

表1 測試環境Tab. 1 Testing enviroment

圖2 計算時間與線程數的關系Fig. 2 Relationship between computing time and number of threads

在研究中發現,計算機同時進行計算的線程數不會超過處理器核數，當使用的線程數超過處理器核數時，超出的線程在前面的線程執行完畢后才會執行。由于每個線程的計算量是一樣的，故在不考慮其他因素的理想情況下，本例中并行計算時間為

(18)

式中：t為計算時間;C為由總任務量決定的常數。

公式(18)說明，當NCPU等于32時，計算效率最高，繼續增加NCPU會降低計算效率。

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

另一方面，使用過多的線程數會占用過多的計算機資源，申請和釋放線程也會消耗一定的時間，從而也會降低計算效率。

從測試結果和以上分析可知，在并行計算的線程數等于計算機處理器核數時，并行計算的效率最高。

3 大規模矩陣存儲方法

在.Net運行環境下，一維數組和多維數組作為單個托管對象(由.Net自動分配與釋放內存)的存儲能力都不超過2 GB，矩陣階數不得超過16 384，這限制了數值模擬的計算規模。本工作通過使用非托管內存實現了大規模矩陣的存儲。

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

AllocHGlobal申請的內存空間是連續的，將矩陣存儲在其中后，可以直接使用Lapack函數庫的矩陣運算函數進行計算，因此在.Net運行環境下，使用非托管內存進行大規模矩陣的存儲是最理想的選擇。

4 大規模線性方程組快速求解方法

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

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

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

目前MKL和ACML提供的Lapack函數庫都是采用Fortran開發的，因此其中的函數由機器語言組成。使用C、C++和Fortran語言都可以通過Lapack函數庫對應的頭文件直接調用其中的函數。使用C#語言可以直接調用基于.Net的動態鏈接庫(DLL文件)中的函數，但不能夠直接調用普通DLL的函數。在C#語言中，可以通過DllImport[19]調用基于機器語言的動態鏈接庫函數。

采用表1中的測試環境，對MKL和ACML中dgetrf和dgetrs函數的求解性能進行測試，測試結果如圖3所示。從圖3中可以看出：在本測試環境下，ACML的求解效率是MKL的兩倍左右，因此在實際使用中應根據處理器類型采用合適的函數庫。

圖3 矩陣求解時間與節點數的關系Fig. 3 Matrix solving time vs number of nodes

5 求解案例

采用兩個實際工程案例對求解性能進行了分析，測試環境如表1所示，并行計算時使用了32個線程。

5.1 張力腿平臺陰極保護數值模擬

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

圖4 張力腿平臺邊界元網格模型Fig. 4 Boundary element model of tension-leg platform

圖5 張力腿平臺犧牲陽極陰極保護的模擬結果Fig. 5 Simulated result of sacrificial anode cathodic protection for tension-leg platform

計算方案生成系數矩陣時長方程組求解時長總時長串行生成&MKL求解566.2113.6679.8并行生成&AMD求解18.057.475.4

5.2 半潛平臺陰極保護數值模擬

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

圖6 半潛平臺邊界元網格模型Fig. 6 Boundary element model of semi-submersible platform

圖7 半潛平臺外加電流陰極保護的模擬結果Fig. 7 Simulated result of impressed current cathodic protection for semi-submersible platform

計算方案生成系數矩陣時長方程組求解時長總時長串行生成&MKL求解456.961.7518.6并行生成&AMD求解14.531.345.8

6 結論

本工作使用C#編程語言，在.Net運行環境下，針對陰極保護邊界元數值模擬，從系數矩陣生成、存儲和大規模線性方程組求解三個方面對高性能計算方法進行了研究，并得到了以下結論。

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

(2) 在.Net環境下，使用非托管內存可以使用計算機空閑的最大連續內存空間，是數值模擬中存儲大規模矩陣的理想選擇；

(3) 在.Net環境下，可以使用Lapack函數庫中的高性能運算函數對線性方程組進行求解，不同處理器應采用相應的函數庫。