崔喜風(fēng)，張紅亮，鄒 忠，徐宇杰，李 劼，張翮輝，賴延清

(中南大學(xué) 冶金科學(xué)與工程學(xué)院，湖南 長沙，410083)

有限元法是隨著電子計算機的發(fā)展而迅速發(fā)展起來的一種現(xiàn)代計算方法,是以計算機和矩陣運算作為工具，對復(fù)雜工程問題或結(jié)構(gòu)進行計算和分析的數(shù)值方法[1-2]。有限元方法的基本求解思想是把計算域劃分為有限個互不重疊的單元，在每個單元內(nèi)，選擇一些合適的節(jié)點作為求解函數(shù)的插值點，將微分方程中的變量改寫成由各變量或其導(dǎo)數(shù)的節(jié)點值與所選用的插值函數(shù)組成的線性表達式，借助于變分原理或加權(quán)余量法，將微分方程離散求解。在有限元方法中，在每個單元內(nèi)選擇基函數(shù)，用單元基函數(shù)的線性組合來逼近單元中的真解。整個計算域上總體的基函數(shù)可以認為由每個單元基函數(shù)組成的，因此，整個計算域內(nèi)的解可以看作是由所有單元上的近似解構(gòu)成[3]。由此可見，單元尺寸越小，計算結(jié)果越精確[4-5]。當網(wǎng)格數(shù)量增加到一定程度后，再繼續(xù)細化，網(wǎng)格精度提高很小，而計算時間卻大幅度增加。所以，應(yīng)注意增加網(wǎng)格數(shù)量的經(jīng)濟性。實際應(yīng)用時，可以比較2種網(wǎng)格劃分的計算結(jié)果，若2次計算結(jié)果相差較大，則可以繼續(xù)增加網(wǎng)格，若相反，則停止計算[6]。鋁電解槽熔體溫度高達960 ℃，腐蝕性強，很難進行溫度和電流分布情況的測定。隨著現(xiàn)代數(shù)學(xué)物理理論、數(shù)值模擬方法、計算機技術(shù)的發(fā)展，有限元計算方法以其適用于求解具有復(fù)雜幾何形狀和復(fù)雜邊界條件的問題而得到迅速發(fā)展，目前在鋁電解槽多物理場仿真計算中得到廣泛的應(yīng)用[7-8]。針對大型預(yù)焙鋁電解槽這種大型模型，單元數(shù)高達百萬數(shù)量級，很容易造成單元數(shù)超過硬件甚至軟件的處理能力而無法計算。因此，在現(xiàn)有計算條件下，尋求既能夠滿足計算精度要求，又能夠保證計算經(jīng)濟性的網(wǎng)格劃分策略，顯得尤為重要。在此，本文作者在HP高性能計算集群的計算節(jié)點HP ProLiant BL460c(CPU為 Intel E5430×2，內(nèi)存DDR2-667 8GB)的硬件平臺上，應(yīng)用軟件ANSYS12.0，建立300 kA鋁電解槽1/4電熱模型，并分別采用不同網(wǎng)格劃分策略，分析電熱計算結(jié)果與網(wǎng)格劃分尺寸的關(guān)系，并以此來確定鋁電解槽電熱計算模型所采用的網(wǎng)格劃分策略。

圖1 模型網(wǎng)格劃分圖Fig.1 Sketch of model mesh

1 實驗?zāi)Ｐ?/h2>

所研究的實驗對象為某300 kA預(yù)焙鋁電解槽，長度為7.536 m，寬度為2.036 m，高度為1.465 m。單元尺寸為70 mm時模型的網(wǎng)格劃分圖如圖1所示。電解槽導(dǎo)電導(dǎo)熱結(jié)構(gòu)由陽極、電解質(zhì)熔體、鋁液、陰極和鋼棒組成，采用SOLID69單元進行網(wǎng)格劃分；其他部分如氧化鋁覆蓋料、電解質(zhì)結(jié)殼、側(cè)部炭塊、槽幫伸腿、底部保溫材料以及槽殼等僅導(dǎo)熱但不導(dǎo)電部分采用SOLID70單元進行網(wǎng)格劃分。陽極鋼爪和陽極導(dǎo)桿結(jié)構(gòu)簡單，直接采用 LINK68線單元建立。定義LINK68單元和炭陽極 SOLID69單元間的電約束方程，以連接2種不同方式建立的結(jié)構(gòu)實現(xiàn)電流的連續(xù)流動。

本模型將鋁電解槽自上至下分為：氧化鋁覆蓋料層(含陽極)，電解質(zhì)結(jié)殼層(含陽極)，含陽極電解質(zhì)層，極間電解質(zhì)層，鋁液層，陰極炭塊無鋼棒層，陰極炭塊含鋼棒層以及底部保溫結(jié)構(gòu)層，如圖2所示。劃分網(wǎng)格時，定義單元尺寸，按照一定的順序調(diào)用ANSYS結(jié)構(gòu)化設(shè)計語言，自上往下分層劃分網(wǎng)格。大部分體均可劃分為六面體網(wǎng)格，但極少部分形狀不規(guī)則體，無法掃過則劃分為四面體網(wǎng)格。

圖2 模型分層示意圖Fig.2 Sketch map of model delamination

電流自陽極進入電解槽，流經(jīng)電解質(zhì)熔體和鋁液，再經(jīng)陰極炭塊和鋼棒流出。各導(dǎo)電部分特別是電解質(zhì)熔體產(chǎn)生大量熱量作為熱源，再通過電解槽外表面與外部空氣的熱對流和熱輻射作用將熱量散失[9-10]。

在數(shù)值計算中，其電場和熱場的控制方程分別如下：

式中：V為電場強度；σ為電導(dǎo)率；T為熱力學(xué)溫度；k為導(dǎo)熱系數(shù)；qs為熱源強度；對于導(dǎo)電部分，qs等于控制單元的熱量，對于非導(dǎo)電部分，qs等于 0。通過電場產(chǎn)生的熱量將電場和熱場耦合起來進行計算。所采用電熱計算方法見文獻[11]，本文所考察的重點為網(wǎng)格密度與計算結(jié)果(熔體溫度、電壓等)及計算效率(迭代次數(shù)、計算時間)的關(guān)系。

2 全局網(wǎng)格密度與仿真精度的關(guān)系

整體網(wǎng)格尺寸為70，60，50，40，35和30 mm時，計算得到的熔體溫度和槽電壓(槽內(nèi)歐姆壓降)以及計算時間如表1所示。

表1 采用不同全局網(wǎng)格尺寸的計算結(jié)果比較Table1 Comparison of calculated results for different overall element sizes

從表1可見：采用不同的單元尺寸進行網(wǎng)格劃分所得的計算結(jié)果相差較大；隨著單元尺寸變小，網(wǎng)格密度變大，網(wǎng)格數(shù)目急劇增加，計算得到的熔體溫度與槽電壓更加精確，但計算時間也大幅度增加。當單元尺寸為30 mm時，由于網(wǎng)格數(shù)量太多在所采用的軟件及硬件平臺上已經(jīng)無法進行計算(硬件系統(tǒng)的內(nèi)存分配及CPU資源已經(jīng)無法滿足要求，若進一步計算，則需要更加先進的軟硬件平臺)，因此，暫時將35 mm確定為計算結(jié)果最接近實際值的單元尺寸。此時，電解質(zhì)溫度為956 ℃，槽電壓為2.785 V，計算共耗時約4 h。

在計算數(shù)據(jù)變化梯度較大的部位要采用較大的網(wǎng)格密度[12]，而鋁電解槽底部、頂部和側(cè)部溫度梯度比較大。因此，分別在整體單元尺寸為70 mm和40 mm的基礎(chǔ)上細化底部、頂部和側(cè)部單元尺寸計算得到的結(jié)果如表2所示。從表2可以看出：細化保溫結(jié)構(gòu)前后計算結(jié)果相差很小。細化保溫結(jié)構(gòu)網(wǎng)格，僅僅使該部分的溫度分布更加精確，但對該部分的散熱量影響極小，并不是造成熔體溫度和槽電壓差異較大的主要原因。

為了考察導(dǎo)電部分的網(wǎng)格細化對于仿真結(jié)果的影響，采用4個方案進行網(wǎng)格劃分：方案1令整體單元尺寸為70 mm；方案2令x和y方向單元尺寸為35 mm，z方向70 mm；方案3令x和y方向單元尺寸為35 mm，z方向?qū)щ娊Y(jié)構(gòu)部分單元尺寸為35 mm，保溫結(jié)構(gòu)部分單元尺寸為70 mm；方案4令整體單元尺寸為35 mm。各方案的計算結(jié)果如表3所示。

表2 保溫結(jié)構(gòu)網(wǎng)格細化前后計算結(jié)果比較Table2 Comparison of calculated results before and after element refining in insulation constructions

表3 導(dǎo)電部分網(wǎng)格細化計算結(jié)果比較Table3 Comparison of calculated results different meshing strategies in conductive parts

比較方案1和2的計算結(jié)果可知：x和y方向網(wǎng)格細化對計算結(jié)果有很大的影響；比較方案3和4的計算結(jié)果可知：導(dǎo)電部分z方向網(wǎng)格細化對計算結(jié)果也有一定影響，但影響較小，保溫結(jié)構(gòu)網(wǎng)格的細化對計算結(jié)果影響比較小，為保證計算時間的經(jīng)濟性，該部分可以忽略。因此，在確定了總體的網(wǎng)格密度之后，便可以對局部的網(wǎng)格密度(主要是整體 x，y方向以及導(dǎo)電部分z方向單元尺寸)進行細致分析，以得到最優(yōu)化的計算效率。

3 局部網(wǎng)格密度與仿真精度的關(guān)系

3.1 水平方向網(wǎng)格劃分策略

保持導(dǎo)電部分z方向單元尺寸為35 mm，保溫結(jié)構(gòu)部分z方向為70 mm，分別取水平的x和y方向單元尺寸為35，30，25和20 mm，計算結(jié)果如表4所示。比較發(fā)現(xiàn)，當取單元尺寸為20 mm時由于網(wǎng)格數(shù)量過大計算機無法處理，而單元尺寸為30 mm與25 mm，熔體溫度相差2 ℃，槽電壓相差6 mV，都比較小。綜合考慮計算精度以及計算時間經(jīng)濟性，水平x和y方向取單元尺寸為30 mm為宜。

表4 不同水平方向單元尺寸計算結(jié)果比較Table4 Results comparison of different xy coordinate element sizes

3.2 豎直方向?qū)щ姴糠终w網(wǎng)格劃分策略

保持x和y方向單元尺寸30 mm，底部保溫結(jié)構(gòu)層z方向單元尺寸70 mm，分別取上部導(dǎo)電部分z方向單元尺寸30，35，40和50 mm，計算結(jié)果如表5所示。從表 5可見：單元尺寸為 35 mm和 40 mm時，計算結(jié)果相差不大，與單元尺寸為30 mm時的計算結(jié)果相比，熔體溫度相差5 ℃，槽電壓相差0.030 V，相對較小。因此，導(dǎo)電部分的豎直方向單元尺寸為40 mm。

表5 不同導(dǎo)電部分豎直方向單元尺寸計算結(jié)果比較Table5 Comparison of calculated results different z coordinate element sizes in conductive parts

比較以上計算結(jié)果，取x和y方向整體單元尺寸30 mm，導(dǎo)電部分z方向單元尺寸40 mm，保溫結(jié)構(gòu)部分z方向單元尺寸70 mm。此時，共迭代10次，計算總耗時約4 h。

3.3 豎直方向局部網(wǎng)格劃分策略

對以上確定的網(wǎng)格劃分策略，雖然保證了計算精度，但計算時間比較長。由于網(wǎng)格劃分的繼承性，自上至下x和y方向單元尺寸不可變更。因此，可以采用分層稀疏化z方向網(wǎng)格的方法，以降低計算時間。

逐步將鋁液層、陰極炭塊含鋼棒層、陰極炭塊無鋼棒層、電解質(zhì)含陽極層、結(jié)殼層和上部氧化鋁覆蓋料層這些局部區(qū)域進行網(wǎng)格稀疏化，將其z方向單元尺寸設(shè)為70 mm。計算結(jié)果和計算時間如表6所示。

從表6可見：將鋁液層、陰極層、電解質(zhì)陽極層z方向網(wǎng)格逐步稀疏化之后計算結(jié)果變化很小，但計算時間降低了近2 h。結(jié)殼層單元尺寸變大后仍然劃分為2層網(wǎng)格，所以計算結(jié)果及計算時間相同。上部氧化鋁層網(wǎng)格稀疏化之后，計算時間增至 216 min，且計算精度降低，不可取。

極間電解質(zhì)層厚度為50 mm，電勢梯度大，網(wǎng)格劃分時卻僅為單層網(wǎng)格。將電解質(zhì)層z方向網(wǎng)格細化為10 mm，計算得熔體溫度957 ℃，槽電壓2.809 V，細化電解質(zhì)層前后計算結(jié)果變化很小。由此可見，電勢梯度也不是網(wǎng)格細化與否的決定性因素。極間電解質(zhì)層仍劃分單層網(wǎng)格。

通過分析鋁電解槽電流密度分布可知：在電流由陽極鋼爪流向陽極炭塊處，即上部氧化鋁覆蓋料層，是線單元與體單元的連接區(qū)域，有集中電流，電流密度很大。以上實驗結(jié)果證明，計算結(jié)果主要受有電流集中的區(qū)域的網(wǎng)格密度的影響。因此，可以確定鋁電解槽電熱場計算精度與網(wǎng)格密度的關(guān)系為：電流分布受網(wǎng)格劃分密度的影響很大，電流密度越大，越應(yīng)將網(wǎng)格細化，所得電場精度越高，進而由電場影響其產(chǎn)生的熱量，最終使得電解槽的熱場精度越高，此時，熱場又反過來通過材料非線性影響到電場，最終得到精度高的計算結(jié)果。

綜合考慮計算精度、計算時間及硬件所能承受極限條件，最終將實驗電解槽的底部保溫結(jié)構(gòu)部分以及陰極層、鋁液層、電解質(zhì)層和結(jié)殼層等無電流集中的區(qū)域z方向單元尺寸定為70 mm，氧化鋁覆蓋料層有電流集中的區(qū)域z方向單元尺寸設(shè)置為40 mm。此時，計算耗時127 min。

表6 豎直方向局部網(wǎng)格稀疏化后計算結(jié)果比較Table6 Comparison of calculated results before and after reducing local element density

4 結(jié)論

(1) 溫度梯度比較大的各部分保溫結(jié)構(gòu)以及電勢梯度比較大的極間電解質(zhì)層，網(wǎng)格劃分尺寸對計算結(jié)果影響較小，無需劃分過密的網(wǎng)格；而電流密度比較大的部位則受網(wǎng)格密度的影響很大，如鋼爪與陽極的結(jié)合部位，此時，為獲得準確的計算結(jié)果，應(yīng)提高這部分的網(wǎng)格密度。因此，對于鋁電解槽電熱場計算，其計算精度主要受其電流集中區(qū)域的網(wǎng)格密度的影響，在網(wǎng)格劃分上可遵循該思路適當調(diào)整全局和局部網(wǎng)格的密度。

(2) 對于本文的算例，上部氧化鋁層處于陽極鋼爪與陽極炭塊的體單元的過渡區(qū)域，有很大的電流集中，故細化了該部分的網(wǎng)格，其x和y方向單元尺寸取30 mm，z方向取40 mm。由于網(wǎng)格劃分的繼承性，為劃分六面體網(wǎng)格，其他部分x，y方向單元尺寸均取30 mm，z方向尺寸則可取70 mm，此時，在保證了計算精度的前提下，迭代次數(shù)最少，計算時間最短。

(3) 該網(wǎng)格劃分思路及方法的確定可為提高鋁電解槽電熱仿真的精度和效率提供參考。

[1] 王瑞, 陳海霞, 王廣峰. ANSYS有限元網(wǎng)格劃分淺析[J]. 天津工業(yè)大學(xué)學(xué)報, 2001, 21(4): 8-11.

WANG Rui, CHEN Hai-xia, WANG Guang-feng. Analysis of ANSYS finite element mesh dividing[J]. Journal of Tianjin Institute of Textile Science and Technology, 2001, 21(4): 8-11.

[2] 李慶齡. ANSYS中網(wǎng)格劃分方法研究[J]. 上海電機學(xué)院學(xué)報,2006, 9(5): 28-30.

LI Qing-ling. Study of mesh dividing methods in ANSYS[J].Journal of Shanghai Dianji University, 2006, 9(5): 28-30.

[3] 陶建華. 水波的數(shù)值模擬[M]. 天津: 天津大學(xué)出版社, 2005:7-9.

TAO Jian-hua. Numerical simulation of water waves[M]. Tianjin:Tianjin University Press, 2005: 7-9.

[4] 金晶, 吳新躍. 有限元網(wǎng)格劃分相關(guān)問題分析研究[J]. 計算機輔助工程, 2005, 14(2): 75-78.

JIN Jing, WU Xin-yue. Analysis of the finite element mesh problems[J]. Computer Aided Engineering, 2005, 14(2): 75-78.

[5] 古成中, 吳新躍. 有限元網(wǎng)格劃分及發(fā)展趨勢[J]. 計算機科學(xué)與探索, 2008(3): 248-259.

GU Cheng-zhong, WU Xin-yue. A review of FEM and trend of development[J]. Journal of Frontiers of Computer Science &Technology, 2008(3): 248-259.

[6] 應(yīng)世洲, 張建, 王國林, 等. 網(wǎng)格密度對全鋼載重子午線輪胎有限元分析的影響[J]. 輪胎工業(yè), 2008, 28(10): 579-582.

YING Shi-zhou, ZHANG Jian, WANG Guo-lin, et al. Effects of mesh density on finite element analysis for TBR tire[J]. 2008,28(10): 579-582.

[7] 李劼, 程迎軍, 賴延清, 等. 大型預(yù)焙鋁電解槽電、熱場的有限元計算[J]. 計算物理, 2003, 20(4): 351-355.

LI Jie, CHENG Ying-jun, LAI Yan-qing, et al. Numerical simulation of current and temperature fields of aluminum reduction cells based on ANSYS[J]. Chinese Journal of Computation Physics, 2003, 20(4): 351-355.

[8] 李劼, 鄧星球, 賴延清, 等. 160 kA預(yù)焙鋁電解槽在低分子比和低溫條件下的三維電熱場[J]. 中南大學(xué)學(xué)報:自然科學(xué)版,2004, 35(6): 875-879.

LI Jie, DENG Xing-qiu, LAI Yan-qing, et al. 3D Thermo-electric of 160 kA prebaked aluminum reduction cell at low cryolitic ratio and temperature[J]. Journal of Central South University of Technology: Natural Science, 2004, 35(6): 875-879.

[9] Arkhipov G V. The mathematical modeling of aluminum reduction cells[J]. JOM, 2006, 58(2): 54-56.

[10] Dupuis M. Thermo-electric design of a 400 kA cell using mathematical models: A tutorial[C]//Peterson R D. Light Metals 2000. Nashville: TMS, 2000: 297-302.

[11] CUI Xi-feng, ZHANG Hong-liang, ZOU Zhong, et al. 3D freeze shape study of the aluminum electrolysis cells using finite element method[C]// Johnson J A. Light Metals 2010. San Diego:TMS, 2010: 447-452.

[12] 朱秀娟. 有限元分析網(wǎng)格劃分的關(guān)鍵技巧[J]. 機械工程與自動化, 2009, 152(1): 185-186.

ZHU Xiu-juan. Pivotal skills in grid division in FEA[J].Mechanical Engineering & Automation, 2009, 152(1): 185-186.