實用COMSOL后處理二次開發技術＊

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

?

實用COMSOL后處理二次開發技術＊

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

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

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

腐蝕；陰極保護；有限元；插值計算

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

1 控制方程

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

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

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

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

當外電源向船體注入電流時，由于海水是強電解質，具有導電性，因此會在船體周圍的一塊水體里產生電場，存在電位分布，如果圍繞船體周圍劃分出一塊適當大小的水體Ω，如圖2所示，則當水體中電位處于穩態時，其分布滿足Laplace方程：

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

1.1 邊界條件

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

1.1.1 電化學反應邊界條件

圖2中：邊界①上發生電化學反應（1）和（2），其電流和電勢關系遵循Tafel公式：

表1 電化學參數

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

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

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

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

1.1.2 電流邊界條件

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

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

外加電流密度是通過 DNV RP B401規范確定的：在船體裸露的情況下，要使船體表面達到保護電位，所需的電流密度為120 mA/m2。考慮到破損率，所需電流密度變為=c×120 mA/m2=16.8 mA/m2。船體總面積為11.564 m2，因此所需要的總電流為：=11.564×16.8=194.3 mA=0.194 3 A。因為有4個電極供電，所以每個電極的供電電流為0.194 3/4=0.048 568 8≈0.05 A。

1.1.3 無限遠邊界條件

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

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

圖3 計算結果

1.2 幾何模型離散并求解

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

2 后處理的二次開發

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

2.1 數據的導出

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

2.2 三角形單元被保護面積的計算

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

表2 數據結構

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 特殊情形的處理

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

圖4 三角形單元

2.2.2 一般情況的處理

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

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

因為在12邊上，電位分布為線性，因此梯度為：

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

式（12）中：1a為點1、a間距離，待求。

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

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

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

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

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

2.3.1 特殊情況

在各種相交關系中，情況1因為單元與截面無交點，因此不予考慮；情況2、3中，單元頂點和截面重合，因此不用計算，直接將單元頂點信息存儲到數據結構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為數據結構，存儲了三角形與截面相交的點的信息；count_point為結構當前存儲的交點數目。情況4、5中的三角形的邊都和截面相交，須通過一定的算法求得交點信息。截面與三角形位置關系分類如表4所示。

表3 三角形頂點電位大小比較關系

情況分類保護面積 情況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上分別截取交點、，如圖5所示。

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

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

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

其中：

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

表4 截面與三角形位置關系分類

三角形單元與截面的位置關系判據 情況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

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

2.3.3 曲線II的求解

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

2.4 算例

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

圖5 交點信息的求解

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

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

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

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

3 結論

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

圖6 精確解與程序解對比

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

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

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

［4］王海冰，王攀達，李文華.基于COMSOL Multiphysics液壓節流閥內部流場數值模擬研究［J］.液壓與氣動，2015（8）：26-29.

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

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

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

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

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

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

中央高校基礎科研基金（編號：HEUCFJ171005）

2095－6835（2018）24－0031－06

TQ150

A

10.15913/j.cnki.kjycx.2018.24.031

〔編輯：嚴麗琴〕