周子杰，劉英偉，張洋，張建能

?

實(shí)用COMSOL后處理二次開發(fā)技術(shù)＊

周子杰，劉英偉，張洋，張建能

（哈爾濱工程大學(xué) 材料科學(xué)與化學(xué)工程學(xué)院，黑龍江 哈爾濱 150001）

采用COMSOL軟件對船體的腐蝕與防護(hù)進(jìn)行了模擬，模擬結(jié)果需要通過后處理模塊顯示出來。一般情況下，COMSOL的后處理模塊所包含的功能能滿足大部分要求，但有些時候，對于一些特殊的要求卻難以滿足。比如計算某一等值線所包圍的面積以及沿某后天截取的曲線上顯示某物理量的分布等。針對這些不足，基于空間幾何和線性插值理論，提出一套算法，并根據(jù)該算法編制了程序進(jìn)行計算。結(jié)果表明，面積計算非常精確，曲線物理量的計算結(jié)果和精確解十分吻合，說明提出的算法正確、程序可靠。該算法也給軟件使用者提供了更多的表示形式，同時也大大地拓展了該軟件的功能。該算法具有通用性，可以移植到其他軟件上。

腐蝕；陰極保護(hù)；有限元；插值計算

COMSOL軟件全稱為COMSOL Multiphysics，是瑞典COMSOL公司于2003年正式推出的全球第一款真正的多物理場耦合分析軟件，它以其方便、易用、高效等優(yōu)點(diǎn)而被廣泛應(yīng)用于各個領(lǐng)域的科學(xué)研究和工程計算當(dāng)中[1-4]。其中的腐蝕模塊專門用來模擬和仿真水或電解質(zhì)溶液中，金屬結(jié)構(gòu)的腐蝕及其防護(hù)，很多學(xué)者使用該模塊進(jìn)行了大量卓有成效的研究[5-7]。該軟件功能強(qiáng)大，但也存在不足，比如對于一個防腐設(shè)計者而言，除了關(guān)心被腐蝕構(gòu)件的表面電位、電流分布外，還關(guān)心構(gòu)件表面腐蝕或未腐蝕面積的大小，即某一等值線所包圍的面積，而COMSOL的后處理模塊卻沒有提供這樣的功能；此外，COMSOL的后處理無法將某一曲線上的物理量以一維圖形顯示出來。對于那些幾何建模時，自然形成的曲線（先天曲線，如圖1（b）中的I）后處理可以實(shí)現(xiàn)上述功能，但對于后天曲線（如圖1（b）中的II、III，它是計算完畢后，用某個截面與船體相截形成的曲線），則現(xiàn)有的后處理模塊不能實(shí)現(xiàn)上述功能。除非在建模開始，就對感興趣的部分進(jìn)行切割以形成先天曲線，但這無疑大大增加了建模的復(fù)雜性，而且對哪些部位感興趣一般也無法確定。針對這些不足，本文開發(fā)了一套程序，該程序能將計算結(jié)果做進(jìn)一步處理，通過一定的算法，能實(shí)現(xiàn)上述功能。這為設(shè)計者提供了更多的、更直接信息，這大大有利于設(shè)計的優(yōu)化，同時也大大擴(kuò)展了該軟件的功能。

1 控制方程

帶有四個輔助陽極的船體幾何模型如圖1所示，陽極材料為惰性有色金屬鉑。當(dāng)船體浸沒在海水中時，除了鉑表面外，其他表面將遭受海水的腐蝕，發(fā)生如下電化學(xué)反應(yīng)：

2Fe→2Fe2++4e. （1）

O2+2H2O+4e→4OH-. （2）

在反應(yīng)（2）中，水中的氧不斷地從船體表面奪取電子而發(fā)生還原反應(yīng)；而反應(yīng)（1）需要不斷地為反應(yīng)（2）提供電子；失去了電子的鐵原子變成離子狀態(tài)，在水分子這種極性分子的作用下，很容易溶解到水中，從而導(dǎo)致了船體的腐蝕。為了阻止腐蝕的發(fā)生，必須抑制反應(yīng)（1）的進(jìn)行。外加電流陰極保護(hù)，就是采用外電源向船體表面注入電流，以代替反應(yīng)（1）。注入的電流使船體表面產(chǎn)生陰極極化，當(dāng)極化電位達(dá)到800 mV時（參比電極為Ag/AgCl/Seawater），表面就會得到保護(hù)。

當(dāng)外電源向船體注入電流時，由于海水是強(qiáng)電解質(zhì)，具有導(dǎo)電性，因此會在船體周圍的一塊水體里產(chǎn)生電場，存在電位分布，如果圍繞船體周圍劃分出一塊適當(dāng)大小的水體Ω，如圖2所示，則當(dāng)水體中電位處于穩(wěn)態(tài)時，其分布滿足Laplace方程：

式（3）中：為電解液電勢；，，∈Ω。

1.1 邊界條件

圖2中的水體Ω具有封閉的表面，其上存在以下幾種邊界條件。

1.1.1 電化學(xué)反應(yīng)邊界條件

圖2中：邊界①上發(fā)生電化學(xué)反應(yīng)（1）和（2），其電流和電勢關(guān)系遵循Tafel公式：

表1 電化學(xué)參數(shù)

參數(shù)鐵氧 平衡電位/V0.760.189 塔菲爾斜率/（V/m）0.410.18 交換電流度/（A/m2）7.7e-77.1e-5 修正后交換電流密度/（A/m2）1.078 7e-79.94e-6 海水電導(dǎo)率/S4

式（6）中：，分別是和涂層厚度以及船體所處環(huán)境有關(guān)的參數(shù)，分別為0.02和0.012；是使用年限，這里按10年計算，這樣計算得到破損率。

確定了破損率以后，按下式修正交換電流密度：

式（7）中：c是修正后的交換電流密度；是破損率為0時的交換電流密度，具體數(shù)據(jù)如表1所示。

1.1.2 電流邊界條件

電極和海水的接觸面②共有4個（圖2（a）只標(biāo)識了2個）。當(dāng)外電源向船體表面注入電流時，電流也會通過接觸面注入到海水中，因此存在如下關(guān)系：

式（8）中：為海水電導(dǎo)率；為單位表面法向矢量；_為注入的電流密度。

外加電流密度是通過 DNV RP B401規(guī)范確定的：在船體裸露的情況下，要使船體表面達(dá)到保護(hù)電位，所需的電流密度為120 mA/m2?？紤]到破損率，所需電流密度變?yōu)?c×120 mA/m2=16.8 mA/m2。船體總面積為11.564 m2，因此所需要的總電流為：=11.564×16.8=194.3 mA=0.194 3 A。因?yàn)橛?個電極供電，所以每個電極的供電電流為0.194 3/4=0.048 568 8≈0.05 A。

1.1.3 無限遠(yuǎn)邊界條件

所謂無限遠(yuǎn)不是空間位置的無限遠(yuǎn)，而是指電場強(qiáng)度很弱，近乎為0的地方，如圖2（a）中的表面③，此處表面電流可看作為0，公式為：

圖2 邊界條件與有限元離散

圖3 計算結(jié)果

1.2 幾何模型離散并求解

設(shè)置完邊界條件后，對幾何模型進(jìn)行有限元離散，結(jié)果如圖2（b）所示，之后就可以進(jìn)行計算。計算結(jié)果如圖3所示。由圖3可知，靠近電極的地方電位較高，電極中間區(qū)域電位低一些?？傮w來看，船體表面電位均小于保護(hù)電位0.8 V，可見整個船體都沒有得到保護(hù)，先前估算的電流偏低，應(yīng)適度增加電流供應(yīng)。

2 后處理的二次開發(fā)

對設(shè)計者而言，有時候還需要了解表面腐蝕面積的大小，即某一等值線所包圍的面積（以下簡稱“面積”），然而COMSOL后處理模塊卻沒有提供這一功能，而這一數(shù)據(jù)對設(shè)計者而言無疑是關(guān)鍵的。另外，計算完畢后，設(shè)計者常常想把某物理量沿感興趣的曲線顯示出來，對于那些建模時自然形成的曲線，COMSOL的后處理模塊能實(shí)現(xiàn)這一功能，但是對于后天通過平面和幾何體相交得到的曲線，后處理模塊卻無法實(shí)現(xiàn)這一功能，這大大限制了該軟件的功能。為此本文開發(fā)了一套程序，通過對計算結(jié)果進(jìn)行二次處理，實(shí)現(xiàn)了上述功能。

2.1 數(shù)據(jù)的導(dǎo)出

由于涉及的計算都和船體表面有關(guān)，因此需要把船體表面信息輸出到文本文件中，再導(dǎo)入程序中處理。船體表面原來是曲面，用四面體單元離散幾何體后，表面就由大量平面三角形組成，因此需將這些三角形單元的信息輸出。這些信息包括單元編號及其所包含的節(jié)點(diǎn)號、節(jié)點(diǎn)的坐標(biāo)、每個節(jié)點(diǎn)的電位等。由于計算域有很多表面，而所處理的表面只是其中一部分，如圖3（a）所示，因此需要將這部分表面單獨(dú)定義。這可以通過COMSOL 軟件中的Definitions功能實(shí)現(xiàn)。接下來要把這部分單獨(dú)定義的表面信息提取出來（COMSOL的計算結(jié)果遍布整個求解域，并不區(qū)分屬于哪個表面），方法是在后處理模塊的Results節(jié)點(diǎn)中，選擇Data Sets命令，把先前定義的表面的有關(guān)數(shù)據(jù)從總的結(jié)果中提取出來；最后，在后處理模塊Results中，選擇Export→Data命令，此時右側(cè)會出現(xiàn)界面，在界面里，可以選擇要處理的數(shù)據(jù)集(先前定義的Data Sets會出現(xiàn))；在界面的Output欄目下，選擇數(shù)據(jù)格式（這里選為Sectionwise）、存儲數(shù)據(jù)的路徑和文件名等。需要注意的是，在該界面的Advanced欄目下，不要選Include header一項(xiàng)，因?yàn)檫@會在數(shù)據(jù)文件中混入英文字符，導(dǎo)致數(shù)據(jù)處理麻煩。當(dāng)選擇、設(shè)置完畢后，按界面上方的Export按鈕，就可以把數(shù)據(jù)輸出到文本文件中了，將文本文件打開，可以看到表2所示的數(shù)據(jù)格式（片段）：前1 828行記錄的是1 828個節(jié)點(diǎn)的坐標(biāo)，其編號就是所在的行號，第1 829～5 323行中，每一行代表一個單元所屬的三個頂點(diǎn)的編號（例如第1 829行為第1行），據(jù)此可知單元數(shù)目為5 323－1 829+1=3 495.從5 324行開始到最后，每一行都代表一個節(jié)點(diǎn)的電位，其中5 324行和第1行相對應(yīng)，其余類推到結(jié)束。

2.2 三角形單元被保護(hù)面積的計算

將上述數(shù)據(jù)導(dǎo)入到程序中，并設(shè)計合適的變量將其儲存起來，然后就可以利用這些數(shù)據(jù)計算保護(hù)面積的大小。由于幾何體采用四面體單元離散，單元形狀函數(shù)為線性[8]，因此四面體的每一個三角形的單元的形狀函數(shù)也為線性，因此單元邊上的電位分布為線性（如圖4（b）所示）。另外需要注意的是，三角形頂點(diǎn)坐標(biāo)是三維的。

表2 數(shù)據(jù)結(jié)構(gòu)

1﹣5.320.867 461 8310.097 393 2﹣5.320.963 158 033﹣0.069 61 3﹣5.320.818 393 3240.089 777 ………… 1 828﹣0.412 73﹣0.014 350.3 1 829123 1 830142 1 831564 ………… 5 3232342531 827 5 3240.801 651 5 3250.802 263 ………… 7 1520.812 09

2.2.1 特殊情形的處理

如果設(shè)保護(hù)電位PROTECT=0.8 V的話，那么三角形3個頂點(diǎn)電位和保護(hù)電位相比存在幾種大小關(guān)系，如表3所示。情況1下，由于3個頂點(diǎn)電位都大于0.8 V，因此，根據(jù)插值函數(shù)的性質(zhì)，整個三角形內(nèi)任意點(diǎn)的電位均大于0.8 V，因此整個單元均得到保護(hù)，保護(hù)面積就等于三角形面積，其大小可根據(jù)相關(guān)公式[9]求得；情況2的情形和情況1相反，3個頂點(diǎn)的電位均小于0.8 V，同樣根據(jù)插值函數(shù)的性質(zhì)，三角形內(nèi)任一點(diǎn)的電位均低于0.8 V，因此整個三角形都沒有被保護(hù)，保護(hù)面積為零。情況3比情況4復(fù)雜些，需進(jìn)一步討論。

圖4 三角形單元

2.2.2 一般情況的處理

以情況3為例，圖4中三角形123的邊12，可用三維空間直線方程描述：

式（10）中：為參數(shù)且∈[0，1]。

因?yàn)樵?2邊上，電位分布為線性，因此梯度為：

這樣，邊12上任意一點(diǎn)的電位可通過線性插值確定：

式（12）中：1a為點(diǎn)1、a間距離，待求。

由于1＞0.8且2＜0.8，因此在1、2點(diǎn)之間，一定存在一點(diǎn)a，使得a=PROTECT=0.8 V；同理邊23上同樣存在，使b=PROTECT=0.8 V。這樣一來，根據(jù)插值函數(shù)的特點(diǎn)，三角形Δ2內(nèi)任意點(diǎn)的電位均低于0.8 V，因此三角形Δ2未得到保護(hù)。

將a=0.8代入式（12），就可以求得1a，再根據(jù)式（10），得到1a的另一表達(dá)式：

二者相等，就確定了a，將a代入式（10），就求出了點(diǎn)的坐標(biāo)。同理，可求出23邊上b點(diǎn)的坐標(biāo)。這樣三角形2三個頂點(diǎn)的坐標(biāo)已知，三角形的面積Δa2b也就確定了，這樣一來，三角形123中只有部分面積得到了保護(hù)，其大小為PROTECTED=Δ123－Δa2b.情況4的處理與此類似，只不過受保護(hù)面積為Δa2b而不是Δ123－Δa2b.這樣，一個三角形單元內(nèi)，受保護(hù)部分的面積就確定下來了，將所有三角形做如此處理，最后將總的保護(hù)面積累加起來就得到整個表面被保護(hù)面積的大小。

2.3 三角形截線上的物理量的計算

這里僅考慮垂直于軸、軸的特殊截面和船體表面交截后形成的曲線，如圖1（b）中的曲線II、III。曲線是由截面和船體相交形成的，由于船體離散后，船體表面實(shí)際上是由平面三角形構(gòu)成，因此截面和船體表面相交，實(shí)際上就是與大量的三角形相交，因而會截得大量直線段，這條曲線，實(shí)際上是這些直線段構(gòu)成的折線。求出線段端點(diǎn)（交點(diǎn)）的空間坐標(biāo)和物理量（電位）后，那么沿著這條折線物理量的分布就得到了。下面以曲線III（是由截面=0=0截得）為例，介紹如何求出上述信息。根據(jù)三角形單元和截面的位置關(guān)系，將二者相交情況分為5類，如表4所示。

2.3.1 特殊情況

在各種相交關(guān)系中，情況1因?yàn)閱卧c截面無交點(diǎn)，因此不予考慮；情況2、3中，單元頂點(diǎn)和截面重合，因此不用計算，直接將單元頂點(diǎn)信息存儲到數(shù)據(jù)結(jié)構(gòu)Point_Cut中，以情況2為例：

count_point=count_point+1

Point_Cut（count_point，1）=x2

Point_Cut（count_point，2）=y2=y0

Point_Cut（count_point，3）=z2

Point_Cut（count_point，4）=V2

其中，Point_Cut為數(shù)據(jù)結(jié)構(gòu)，存儲了三角形與截面相交的點(diǎn)的信息；count_point為結(jié)構(gòu)當(dāng)前存儲的交點(diǎn)數(shù)目。情況4、5中的三角形的邊都和截面相交，須通過一定的算法求得交點(diǎn)信息。截面與三角形位置關(guān)系分類如表4所示。

表3 三角形頂點(diǎn)電位大小比較關(guān)系

情況分類保護(hù)面積 情況1V1＞0.8，V2＞0.8，V3＞0.8SPROTECTED=SΔ123 情況2V1＜0.8，V2＜0.8，V3＜0.8SPROTECTED=0 情況3V1＞0.8，V2＜0.8，V3＞0.8SPROTECTED=SΔ123-SΔa2b 情況4V1＜0.8，V2＞0.8，V3＜0.8SPROTECTED=SΔa2b

2.3.2 一般情況

以情況5中的三角形①為例，截面（=0）與三角形123相交，在邊12和13上分別截取交點(diǎn)、，如圖5所示。

三角形123的邊12，用形如公式（10）的三維直線方程描述：

由于n=0，代入式（13）求得：

將1代入式（13），就可確定n、n。同理可確定點(diǎn)坐標(biāo)（m，m，m）。接下來求、點(diǎn)的物理量——電位，這可以根據(jù)式（12）確定：

其中：

1、2分別為線段12和13上電位梯度：

表4 截面與三角形位置關(guān)系分類

三角形單元與截面的位置關(guān)系判據(jù) 情況1Cut plane（y=y0）①（y1－y0）＜0∩（y2－y0）＜0∩（y3－y0）＜0；②（y1－y0）＞0∩（y2－y0）＞0∩（y3－y0）＞0 情況2Cut plane（y=y0）①（y1－y0）＜0∩（y2－y0）=0∩（y3－y0）＜0 ；②（y1－y0）＞0∩（y2－y0）=0∩（y3－y0）＞0 情況3Cut plane（y=y0）①（y1－y0）=0∩（y2－y0）=0∩（y3－y0）＜0；②（y1－y0）=0∩（y2－y0）=0∩（y3－y0）＞0 情況4Cut plane（y=y0）（y1－y0）=0∩（y2－y0）＞0∩（y3－y0）＜0 情況5Cut plane（y=y0）①（y1－y0）>0∩（y2－y0）<0∩（y3－y0）<0；②（y1－y0）>0∩（y2－y0）>0∩（y3－y0）<0

至此，三角形與某一截面相交，交點(diǎn)的坐標(biāo)以及電位均可求得，將這些結(jié)果加入到結(jié)構(gòu)Point_Cut中。每個三角形都這樣處理，最終就得到了一個由截面和船體表面相截而得到的曲線上若干交點(diǎn)的信息，經(jīng)過進(jìn)一步處理（去掉冗余、排序等），就可以把這些物理量沿曲線的分布繪制出來，如圖6（d）、6（f）所示。

2.3.3 曲線II的求解

曲線II和曲線III的不同在于二者所垂直的坐標(biāo)軸不同，但在算法上，十分類似，只需對2.3.2中的算法略作修改即可。此時曲線II是由垂直于軸的平面所截，設(shè)此平面方程為=0，將它代入式（13）中第一式，就可確定參數(shù)，進(jìn)而可確定節(jié)點(diǎn)坐標(biāo)，接下來，參照2.3.2中的余下步驟就可確定節(jié)點(diǎn)電位，這里不再贅述。

2.4 算例

由圖3（a）可知，船體表面未得到有效保護(hù)，必須增大電流密度，現(xiàn)將電流密度增加到0.1 A/m2，結(jié)果如圖6（a）、6（b）所示，此時，表面電位有所增大。用開發(fā)的程序計分別算了等值線為0.795 2 V、0.8 V（保護(hù)電位）和0 V所包圍的面積，結(jié)果如表5所示。可見，盡管0.795 2 V和0.8 V相差并不大，可所包圍的面積卻相差很大，不過由于0.795 2 V和0.8 V相差不大，因此可以認(rèn)為船體表面大部分得到了保護(hù)。另外，等值線0所包圍的面積就是船體底面總面積，計算結(jié)果為12.193 9 m2。

圖5 交點(diǎn)信息的求解

表5 不同電位等值線包圍面積

電流密度/（A/m2）0.1 等值線值/V0.795 20.80 包圍面積/m27.727 12.359 312.193 9

COMSOL后處理具有計算某一表面的總面積的功能，是在Results→Derived Values菜單里實(shí)現(xiàn)的。為了驗(yàn)證算法的正確性，本文利用COMSOL提供的這一功能計算了船體表面積，如圖3（a）所示，結(jié)果為12.194 m2，而本程序計的結(jié)果為12.193 9 m2，二者基本一致（考慮舍入的話，二者幾乎相等），這證明了算法正確性以及程序的可靠性。

圖6（c）、6（e）是在建模時，在相關(guān)位置進(jìn)行切割得到曲線(先天曲線)，經(jīng)程序計算后，得到的物理量沿曲線精確分布；6（d）、6（f）分別為沒有經(jīng)過切割，直接采用本文開發(fā)的程序，在相同部位曲線上計算得到的物理量的分布。由圖可以看出，二者符合得非常好，說明算法的正確性和程序的可靠性。

3 結(jié)論

本文開發(fā)的程序能夠計算某一等值線所包圍的面積，以及某一三維截線上的物理量分布，這大大彌補(bǔ)了該軟件后處理的不足。計算結(jié)果證明：本文提出的算法正確、程序可靠；計算結(jié)果能夠給設(shè)計者提供準(zhǔn)確的信息，大大有利于改進(jìn)設(shè)計；該算法具有可移植性。

圖6 精確解與程序解對比

［1］劉芊，曹江勇，羅勇，等.基于COMSOL Multiphysics的磁場仿真分析［J］.大學(xué)物理實(shí)驗(yàn)，2015，28（5）：106-108.

［2］周靜雷，王源.COMSOL Multiphysics在微型揚(yáng)聲器熱效應(yīng)中的應(yīng)用［J］.電子測量技術(shù)，2016，39（2）：29-32.

［3］高霞，王志斌.基于COMSOL Multiphysics壓電鈮酸鋰晶片仿真［J］.壓電與聲光，2015，37（2）：291-293.

［4］王海冰，王攀達(dá)，李文華.基于COMSOL Multiphysics液壓節(jié)流閥內(nèi)部流場數(shù)值模擬研究［J］.液壓與氣動，2015（8）：26-29.

［5］辛艷萍，梁月.基于COMSOL Multiphysics的金屬儲罐陰極保護(hù)方案的優(yōu)化［J］.石油化工高等學(xué)校學(xué)報，2016，29（4）：92-96

［6］王雷，董麗娜.基于COMSOL Multiphysics的雜散電流腐蝕仿真分析［J］.新技術(shù)新工藝，2014（1）：22-24.

［7］甘甜，仵杰，邢德鍵.基于COMSOL帶腐蝕凹陷套管井的聲場數(shù)值［J］.電聲技術(shù)，2017，21（4）：111-114．

［8］王勖成.有限單元法［M］.北京：清華大學(xué)出版社，2003.

［9］王煥定，王偉.有限單元法教程［M］.哈爾濱：哈爾濱工業(yè)大學(xué)出版社，2003.

周子杰（1995—），男，山東濟(jì)寧人，碩士研究生。

中央高?；A(chǔ)科研基金（編號：HEUCFJ171005）

2095－6835（2018）24－0031－06

TQ150

A

10.15913/j.cnki.kjycx.2018.24.031

〔編輯：嚴(yán)麗琴〕