崔喜風(fēng),張紅亮,鄒 忠,徐宇杰,李 劼,張翮輝,賴延清
(中南大學(xué) 冶金科學(xué)與工程學(xué)院,湖南 長沙,410083)
有限元法是隨著電子計(jì)算機(jī)的發(fā)展而迅速發(fā)展起來的一種現(xiàn)代計(jì)算方法,是以計(jì)算機(jī)和矩陣運(yùn)算作為工具,對(duì)復(fù)雜工程問題或結(jié)構(gòu)進(jìn)行計(jì)算和分析的數(shù)值方法[1-2]。有限元方法的基本求解思想是把計(jì)算域劃分為有限個(gè)互不重疊的單元,在每個(gè)單元內(nèi),選擇一些合適的節(jié)點(diǎn)作為求解函數(shù)的插值點(diǎn),將微分方程中的變量改寫成由各變量或其導(dǎo)數(shù)的節(jié)點(diǎn)值與所選用的插值函數(shù)組成的線性表達(dá)式,借助于變分原理或加權(quán)余量法,將微分方程離散求解。在有限元方法中,在每個(gè)單元內(nèi)選擇基函數(shù),用單元基函數(shù)的線性組合來逼近單元中的真解。整個(gè)計(jì)算域上總體的基函數(shù)可以認(rèn)為由每個(gè)單元基函數(shù)組成的,因此,整個(gè)計(jì)算域內(nèi)的解可以看作是由所有單元上的近似解構(gòu)成[3]。由此可見,單元尺寸越小,計(jì)算結(jié)果越精確[4-5]。當(dāng)網(wǎng)格數(shù)量增加到一定程度后,再繼續(xù)細(xì)化,網(wǎng)格精度提高很小,而計(jì)算時(shí)間卻大幅度增加。所以,應(yīng)注意增加網(wǎng)格數(shù)量的經(jīng)濟(jì)性。實(shí)際應(yīng)用時(shí),可以比較2種網(wǎng)格劃分的計(jì)算結(jié)果,若2次計(jì)算結(jié)果相差較大,則可以繼續(xù)增加網(wǎng)格,若相反,則停止計(jì)算[6]。鋁電解槽熔體溫度高達(dá)960 ℃,腐蝕性強(qiáng),很難進(jìn)行溫度和電流分布情況的測(cè)定。隨著現(xiàn)代數(shù)學(xué)物理理論、數(shù)值模擬方法、計(jì)算機(jī)技術(shù)的發(fā)展,有限元計(jì)算方法以其適用于求解具有復(fù)雜幾何形狀和復(fù)雜邊界條件的問題而得到迅速發(fā)展,目前在鋁電解槽多物理場仿真計(jì)算中得到廣泛的應(yīng)用[7-8]。針對(duì)大型預(yù)焙鋁電解槽這種大型模型,單元數(shù)高達(dá)百萬數(shù)量級(jí),很容易造成單元數(shù)超過硬件甚至軟件的處理能力而無法計(jì)算。因此,在現(xiàn)有計(jì)算條件下,尋求既能夠滿足計(jì)算精度要求,又能夠保證計(jì)算經(jīng)濟(jì)性的網(wǎng)格劃分策略,顯得尤為重要。在此,本文作者在HP高性能計(jì)算集群的計(jì)算節(jié)點(diǎn)HP ProLiant BL460c(CPU為 Intel E5430×2,內(nèi)存DDR2-667 8GB)的硬件平臺(tái)上,應(yīng)用軟件ANSYS12.0,建立300 kA鋁電解槽1/4電熱模型,并分別采用不同網(wǎng)格劃分策略,分析電熱計(jì)算結(jié)果與網(wǎng)格劃分尺寸的關(guān)系,并以此來確定鋁電解槽電熱計(jì)算模型所采用的網(wǎng)格劃分策略。
圖1 模型網(wǎng)格劃分圖Fig.1 Sketch of model mesh
所研究的實(shí)驗(yàn)對(duì)象為某300 kA預(yù)焙鋁電解槽,長度為7.536 m,寬度為2.036 m,高度為1.465 m。單元尺寸為70 mm時(shí)模型的網(wǎng)格劃分圖如圖1所示。電解槽導(dǎo)電導(dǎo)熱結(jié)構(gòu)由陽極、電解質(zhì)熔體、鋁液、陰極和鋼棒組成,采用SOLID69單元進(jìn)行網(wǎng)格劃分;其他部分如氧化鋁覆蓋料、電解質(zhì)結(jié)殼、側(cè)部炭塊、槽幫伸腿、底部保溫材料以及槽殼等僅導(dǎo)熱但不導(dǎo)電部分采用SOLID70單元進(jìn)行網(wǎng)格劃分。陽極鋼爪和陽極導(dǎo)桿結(jié)構(gòu)簡單,直接采用 LINK68線單元建立。定義LINK68單元和炭陽極 SOLID69單元間的電約束方程,以連接2種不同方式建立的結(jié)構(gòu)實(shí)現(xiàn)電流的連續(xù)流動(dòng)。
本模型將鋁電解槽自上至下分為:氧化鋁覆蓋料層(含陽極),電解質(zhì)結(jié)殼層(含陽極),含陽極電解質(zhì)層,極間電解質(zhì)層,鋁液層,陰極炭塊無鋼棒層,陰極炭塊含鋼棒層以及底部保溫結(jié)構(gòu)層,如圖2所示。劃分網(wǎng)格時(shí),定義單元尺寸,按照一定的順序調(diào)用ANSYS結(jié)構(gòu)化設(shè)計(jì)語言,自上往下分層劃分網(wǎng)格。大部分體均可劃分為六面體網(wǎng)格,但極少部分形狀不規(guī)則體,無法掃過則劃分為四面體網(wǎng)格。
圖2 模型分層示意圖Fig.2 Sketch map of model delamination
電流自陽極進(jìn)入電解槽,流經(jīng)電解質(zhì)熔體和鋁液,再經(jīng)陰極炭塊和鋼棒流出。各導(dǎo)電部分特別是電解質(zhì)熔體產(chǎn)生大量熱量作為熱源,再通過電解槽外表面與外部空氣的熱對(duì)流和熱輻射作用將熱量散失[9-10]。
在數(shù)值計(jì)算中,其電場和熱場的控制方程分別如下:
式中:V為電場強(qiáng)度;σ為電導(dǎo)率;T為熱力學(xué)溫度;k為導(dǎo)熱系數(shù);qs為熱源強(qiáng)度;對(duì)于導(dǎo)電部分,qs等于控制單元的熱量,對(duì)于非導(dǎo)電部分,qs等于 0。通過電場產(chǎn)生的熱量將電場和熱場耦合起來進(jìn)行計(jì)算。所采用電熱計(jì)算方法見文獻(xiàn)[11],本文所考察的重點(diǎn)為網(wǎng)格密度與計(jì)算結(jié)果(熔體溫度、電壓等)及計(jì)算效率(迭代次數(shù)、計(jì)算時(shí)間)的關(guān)系。
整體網(wǎng)格尺寸為70,60,50,40,35和30 mm時(shí),計(jì)算得到的熔體溫度和槽電壓(槽內(nèi)歐姆壓降)以及計(jì)算時(shí)間如表1所示。
表1 采用不同全局網(wǎng)格尺寸的計(jì)算結(jié)果比較Table1 Comparison of calculated results for different overall element sizes
從表1可見:采用不同的單元尺寸進(jìn)行網(wǎng)格劃分所得的計(jì)算結(jié)果相差較大;隨著單元尺寸變小,網(wǎng)格密度變大,網(wǎng)格數(shù)目急劇增加,計(jì)算得到的熔體溫度與槽電壓更加精確,但計(jì)算時(shí)間也大幅度增加。當(dāng)單元尺寸為30 mm時(shí),由于網(wǎng)格數(shù)量太多在所采用的軟件及硬件平臺(tái)上已經(jīng)無法進(jìn)行計(jì)算(硬件系統(tǒng)的內(nèi)存分配及CPU資源已經(jīng)無法滿足要求,若進(jìn)一步計(jì)算,則需要更加先進(jìn)的軟硬件平臺(tái)),因此,暫時(shí)將35 mm確定為計(jì)算結(jié)果最接近實(shí)際值的單元尺寸。此時(shí),電解質(zhì)溫度為956 ℃,槽電壓為2.785 V,計(jì)算共耗時(shí)約4 h。
在計(jì)算數(shù)據(jù)變化梯度較大的部位要采用較大的網(wǎng)格密度[12],而鋁電解槽底部、頂部和側(cè)部溫度梯度比較大。因此,分別在整體單元尺寸為70 mm和40 mm的基礎(chǔ)上細(xì)化底部、頂部和側(cè)部單元尺寸計(jì)算得到的結(jié)果如表2所示。從表2可以看出:細(xì)化保溫結(jié)構(gòu)前后計(jì)算結(jié)果相差很小。細(xì)化保溫結(jié)構(gòu)網(wǎng)格,僅僅使該部分的溫度分布更加精確,但對(duì)該部分的散熱量影響極小,并不是造成熔體溫度和槽電壓差異較大的主要原因。
為了考察導(dǎo)電部分的網(wǎng)格細(xì)化對(duì)于仿真結(jié)果的影響,采用4個(gè)方案進(jìn)行網(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。各方案的計(jì)算結(jié)果如表3所示。
表2 保溫結(jié)構(gòu)網(wǎng)格細(xì)化前后計(jì)算結(jié)果比較Table2 Comparison of calculated results before and after element refining in insulation constructions
表3 導(dǎo)電部分網(wǎng)格細(xì)化計(jì)算結(jié)果比較Table3 Comparison of calculated results different meshing strategies in conductive parts
比較方案1和2的計(jì)算結(jié)果可知:x和y方向網(wǎng)格細(xì)化對(duì)計(jì)算結(jié)果有很大的影響;比較方案3和4的計(jì)算結(jié)果可知:導(dǎo)電部分z方向網(wǎng)格細(xì)化對(duì)計(jì)算結(jié)果也有一定影響,但影響較小,保溫結(jié)構(gòu)網(wǎng)格的細(xì)化對(duì)計(jì)算結(jié)果影響比較小,為保證計(jì)算時(shí)間的經(jīng)濟(jì)性,該部分可以忽略。因此,在確定了總體的網(wǎng)格密度之后,便可以對(duì)局部的網(wǎng)格密度(主要是整體 x,y方向以及導(dǎo)電部分z方向單元尺寸)進(jìn)行細(xì)致分析,以得到最優(yōu)化的計(jì)算效率。
保持導(dǎo)電部分z方向單元尺寸為35 mm,保溫結(jié)構(gòu)部分z方向?yàn)?0 mm,分別取水平的x和y方向單元尺寸為35,30,25和20 mm,計(jì)算結(jié)果如表4所示。比較發(fā)現(xiàn),當(dāng)取單元尺寸為20 mm時(shí)由于網(wǎng)格數(shù)量過大計(jì)算機(jī)無法處理,而單元尺寸為30 mm與25 mm,熔體溫度相差2 ℃,槽電壓相差6 mV,都比較小。綜合考慮計(jì)算精度以及計(jì)算時(shí)間經(jīng)濟(jì)性,水平x和y方向取單元尺寸為30 mm為宜。
表4 不同水平方向單元尺寸計(jì)算結(jié)果比較Table4 Results comparison of different xy coordinate element sizes
保持x和y方向單元尺寸30 mm,底部保溫結(jié)構(gòu)層z方向單元尺寸70 mm,分別取上部導(dǎo)電部分z方向單元尺寸30,35,40和50 mm,計(jì)算結(jié)果如表5所示。從表 5可見:單元尺寸為 35 mm和 40 mm時(shí),計(jì)算結(jié)果相差不大,與單元尺寸為30 mm時(shí)的計(jì)算結(jié)果相比,熔體溫度相差5 ℃,槽電壓相差0.030 V,相對(duì)較小。因此,導(dǎo)電部分的豎直方向單元尺寸為40 mm。
表5 不同導(dǎo)電部分豎直方向單元尺寸計(jì)算結(jié)果比較Table5 Comparison of calculated results different z coordinate element sizes in conductive parts
比較以上計(jì)算結(jié)果,取x和y方向整體單元尺寸30 mm,導(dǎo)電部分z方向單元尺寸40 mm,保溫結(jié)構(gòu)部分z方向單元尺寸70 mm。此時(shí),共迭代10次,計(jì)算總耗時(shí)約4 h。
對(duì)以上確定的網(wǎng)格劃分策略,雖然保證了計(jì)算精度,但計(jì)算時(shí)間比較長。由于網(wǎng)格劃分的繼承性,自上至下x和y方向單元尺寸不可變更。因此,可以采用分層稀疏化z方向網(wǎng)格的方法,以降低計(jì)算時(shí)間。
逐步將鋁液層、陰極炭塊含鋼棒層、陰極炭塊無鋼棒層、電解質(zhì)含陽極層、結(jié)殼層和上部氧化鋁覆蓋料層這些局部區(qū)域進(jìn)行網(wǎng)格稀疏化,將其z方向單元尺寸設(shè)為70 mm。計(jì)算結(jié)果和計(jì)算時(shí)間如表6所示。
從表6可見:將鋁液層、陰極層、電解質(zhì)陽極層z方向網(wǎng)格逐步稀疏化之后計(jì)算結(jié)果變化很小,但計(jì)算時(shí)間降低了近2 h。結(jié)殼層單元尺寸變大后仍然劃分為2層網(wǎng)格,所以計(jì)算結(jié)果及計(jì)算時(shí)間相同。上部氧化鋁層網(wǎng)格稀疏化之后,計(jì)算時(shí)間增至 216 min,且計(jì)算精度降低,不可取。
極間電解質(zhì)層厚度為50 mm,電勢(shì)梯度大,網(wǎng)格劃分時(shí)卻僅為單層網(wǎng)格。將電解質(zhì)層z方向網(wǎng)格細(xì)化為10 mm,計(jì)算得熔體溫度957 ℃,槽電壓2.809 V,細(xì)化電解質(zhì)層前后計(jì)算結(jié)果變化很小。由此可見,電勢(shì)梯度也不是網(wǎng)格細(xì)化與否的決定性因素。極間電解質(zhì)層仍劃分單層網(wǎng)格。
通過分析鋁電解槽電流密度分布可知:在電流由陽極鋼爪流向陽極炭塊處,即上部氧化鋁覆蓋料層,是線單元與體單元的連接區(qū)域,有集中電流,電流密度很大。以上實(shí)驗(yàn)結(jié)果證明,計(jì)算結(jié)果主要受有電流集中的區(qū)域的網(wǎng)格密度的影響。因此,可以確定鋁電解槽電熱場計(jì)算精度與網(wǎng)格密度的關(guān)系為:電流分布受網(wǎng)格劃分密度的影響很大,電流密度越大,越應(yīng)將網(wǎng)格細(xì)化,所得電場精度越高,進(jìn)而由電場影響其產(chǎn)生的熱量,最終使得電解槽的熱場精度越高,此時(shí),熱場又反過來通過材料非線性影響到電場,最終得到精度高的計(jì)算結(jié)果。
綜合考慮計(jì)算精度、計(jì)算時(shí)間及硬件所能承受極限條件,最終將實(shí)驗(yàn)電解槽的底部保溫結(jié)構(gòu)部分以及陰極層、鋁液層、電解質(zhì)層和結(jié)殼層等無電流集中的區(qū)域z方向單元尺寸定為70 mm,氧化鋁覆蓋料層有電流集中的區(qū)域z方向單元尺寸設(shè)置為40 mm。此時(shí),計(jì)算耗時(shí)127 min。
表6 豎直方向局部網(wǎng)格稀疏化后計(jì)算結(jié)果比較Table6 Comparison of calculated results before and after reducing local element density
(1) 溫度梯度比較大的各部分保溫結(jié)構(gòu)以及電勢(shì)梯度比較大的極間電解質(zhì)層,網(wǎng)格劃分尺寸對(duì)計(jì)算結(jié)果影響較小,無需劃分過密的網(wǎng)格;而電流密度比較大的部位則受網(wǎng)格密度的影響很大,如鋼爪與陽極的結(jié)合部位,此時(shí),為獲得準(zhǔn)確的計(jì)算結(jié)果,應(yīng)提高這部分的網(wǎng)格密度。因此,對(duì)于鋁電解槽電熱場計(jì)算,其計(jì)算精度主要受其電流集中區(qū)域的網(wǎng)格密度的影響,在網(wǎng)格劃分上可遵循該思路適當(dāng)調(diào)整全局和局部網(wǎng)格的密度。
(2) 對(duì)于本文的算例,上部氧化鋁層處于陽極鋼爪與陽極炭塊的體單元的過渡區(qū)域,有很大的電流集中,故細(xì)化了該部分的網(wǎng)格,其x和y方向單元尺寸取30 mm,z方向取40 mm。由于網(wǎng)格劃分的繼承性,為劃分六面體網(wǎng)格,其他部分x,y方向單元尺寸均取30 mm,z方向尺寸則可取70 mm,此時(shí),在保證了計(jì)算精度的前提下,迭代次數(shù)最少,計(jì)算時(shí)間最短。
(3) 該網(wǎng)格劃分思路及方法的確定可為提高鋁電解槽電熱仿真的精度和效率提供參考。
[1] 王瑞, 陳海霞, 王廣峰. ANSYS有限元網(wǎng)格劃分淺析[J]. 天津工業(yè)大學(xué)學(xué)報(bào), 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]. 上海電機(jī)學(xué)院學(xué)報(bào),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]. 計(jì)算機(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ā)展趨勢(shì)[J]. 計(jì)算機(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)格密度對(duì)全鋼載重子午線輪胎有限元分析的影響[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ì)算[J]. 計(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é)報(bào):自然科學(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]. 機(jī)械工程與自動(dòng)化, 2009, 152(1): 185-186.
ZHU Xiu-juan. Pivotal skills in grid division in FEA[J].Mechanical Engineering & Automation, 2009, 152(1): 185-186.