基于多種湍流模型的電解加工溫度多場耦合仿真

陳遠龍 林華,2 陳培譞 李回歸

(1.合肥工業(yè)大學(xué) 機械工程學(xué)院,安徽 合肥 230009；2.皖西學(xué)院 機械與車輛工程學(xué)院,安徽 六安 237000)

電解加工由于其獨特的加工優(yōu)勢，被廣泛應(yīng)用于航空、航天、兵器等國防軍工及高端裝備領(lǐng)域的難加工材料、復(fù)雜型面零件的精密加工等，有效支撐和推動了我國航空航天事業(yè)的蓬勃發(fā)展[1- 4]。電解加工過程是多物理場綜合作用的過程，加工間隙內(nèi)電解液的溫度分布不均，對電導(dǎo)率分布有重要影響，使得電場分布不均，最終影響著工件的成型精度?；趩挝锢韴龊投辔锢韴鲴詈系姆抡婕夹g(shù)在電解加工技術(shù)研究中有著越來越廣泛的應(yīng)用，有效提升了加工精度，縮短了研發(fā)周期[5- 6]。

多物理場比單物理場的仿真分析更能表達電解加工過程中多場間錯綜復(fù)雜的關(guān)系，目前電解加工過程仿真模型已由純電場或流場的單場仿真向多場耦合仿真的方向發(fā)展[7- 9]。Deconinck等[10]基于層流和多離子傳輸模型研究了電解加工的溫度分布和材料去除規(guī)律。周小超等[5]基于歐拉雙流體氣液兩相流模型進行了葉片電解加工過程多場耦合研究。陳遠龍等[11]基于湍流k-ε模型進行了高頻脈沖電化學(xué)加工過程多場耦合仿真研究。劉國強等[12]對小孔內(nèi)擴孔的電解加工過程進行了多場耦合仿真研究。劉強等[13]基于k-ε模型對航空發(fā)動機中葉片式擴壓器電解加工的流場進行了優(yōu)化設(shè)計。朱彥偉等[14]對采用不同湍流模型對渦輪葉片表面換熱計算的影響進行了比較，并指出在不同計算域需采用不同湍流模型和壁面處理函數(shù)進行計算。可見，在仿真中使用不同的湍流模型，流場求解結(jié)果也不盡相同，尤其是不同流場模型在耦合電場、溫度場等物理場時，計算的結(jié)果差別更大，很大程度上影響著整個電解加工過程和加工精度分析。

本文針對型面電解加工過程中的間隙溫度分布不均勻、難以預(yù)測和測量的問題，建立基于兩相流場模型并耦合電場、流場和溫度場的溫度多場耦合仿真模型，采用多種流場模型和壁處理方式進行間隙流速和溫度分布的求解，分析其求解精度并與試驗結(jié)果進行對比，以獲得具有較高求解精度的電解加工溫度多場耦合模型及其求解方法。

1 幾何模型

假設(shè)加工進入平衡狀態(tài)，建立型面電解加工仿真的二維幾何模型如圖1所示。采用側(cè)流式供液方式，左側(cè)為電解液入口，右側(cè)為電解液出口。工件陽極尺寸為100 mm×90 mm，工具陰極尺寸為100 mm×40 mm，厚度尺寸為12.7 mm，曲率半徑為70 mm，加工間隙為0.5 mm，Ω為電解液區(qū)域。在加工間隙中設(shè)置6個溫度觀測點1～6，各點間距18 mm，點1和點6距左右端面各5 mm，均高出陰極表面0.05 mm。O1O2為間隙內(nèi)流速觀察線。

圖1 電解加工仿真幾何模型

2 多場耦合分析

2.1 電場

電場模型如圖1所示，Γ1為工件陽極邊界，Γ2為工具陰極邊界，Γ3=U，Γ4=0，U為加工電壓。根據(jù)電荷守恒定律和歐姆定律，由電場控制方程，電解液中電位滿足[2]:

?·(κ?φ)=0

(1)

式中：κ為電解液電導(dǎo)率，φ為電解液電位。在電解液入口和出口滿足

?φ/?n=0

(2)

式中：n為電位梯度方向上的單位矢量。電解加工中電解液電導(dǎo)率會受到電解液溫升和氣泡率的影響，其模型可描述為[2]：

κ=κ0[1+ξ(T-T0)](1-β)m

(3)

式中：κ0、T0為電解液初始電導(dǎo)率和溫度；ξ為溫度系數(shù)，取0.02；β為氣泡率；m為氣泡率影響指數(shù)，取2。

2.2 流場

電解加工中，流動的電解液要足以排出間隙中的電解產(chǎn)物與所產(chǎn)生的熱量，要能減小電極附近的濃差極化，因此其必須具有一定的速度，呈湍流狀態(tài)，且流速均勻性要好。為簡化計算，假設(shè)電解液為理想狀態(tài)的液體，不含氣泡、固體顆粒等，為連續(xù)不可壓縮黏性流體。由質(zhì)量守恒定律和動量守恒定律，電解液流體的流動應(yīng)滿足[5]：

(4)

式中：ρ為電解液密度；p為電解液壓力；v為電解液流速；μ為電解液動力黏度；μT為湍流黏性系數(shù)。

2.2.1 單相流場模型

湍流的計算方法主要分為3類：雷諾時均模擬、尺度解析模擬和直接數(shù)值模擬，前兩類可以看成是非直接模擬。常用的湍流模型有SA模型、k-ε模型、k-ω模型、SST模型、雷諾應(yīng)力模型和大渦模擬等[14]。

(1)Spalart-Allmaras模型

湍流Spalart-Allmaras(簡稱SA模型)模型是單方程模型，只需求解湍流黏性的輸運方程。該模型通常用于模擬中等復(fù)雜的內(nèi)流和外流以及壓力梯度下的邊界層流流動。

(2)標(biāo)準(zhǔn)k-ε模型

湍流標(biāo)準(zhǔn)k-ε模型(簡稱k-ε模型)是兩個方程模型，要解兩個變量，即速度和長度尺度，其在現(xiàn)代工程流場計算中的應(yīng)用最為廣泛，但由于其假設(shè)流動為完全湍流，故只適用于完全湍流的流動過程模型。模型中湍動能輸運方程是通過精確的方程推導(dǎo)得到的，而耗散率方程則是通過物理推理及數(shù)學(xué)上模擬相似原型方程得到。k-ε模型液相湍動能及耗散率方程為[11]：

(5)

式中：k為湍動能；ε為湍流耗散率；σk、σε、C1ε和C2ε為模型常數(shù)；Pk為平均速度梯度產(chǎn)生的湍動能。在標(biāo)準(zhǔn)k-ε模型中，模型的常數(shù)取值為

C1ε=1.44，C2ε=1.92，σk=1.0，σε=1.3。

(3)標(biāo)準(zhǔn)k-ω模型

湍流標(biāo)準(zhǔn)k-ω模型(簡稱k-ω模型)基于Wilcoxk-ω模型，并考慮低雷諾數(shù)、可壓縮性和剪切流傳播對其修改得到[13]。模型預(yù)測了自由剪切流傳播速率，可以應(yīng)用于墻壁束縛流動和自由剪切流動。k-ω方程克服了k-ε方程在近壁流動模擬的缺陷，對于近壁流動或存在逆壓梯度流動的湍流計算時具有較大優(yōu)勢，同時還考慮到了湍流剪應(yīng)力的影響修改了湍流黏性公式。在近壁流動模型中k-ω模型精度較高，其適用于近壁低雷諾數(shù)區(qū)域的模擬，但比k-ε模型收斂困難，計算結(jié)果對初始條件較敏感。

(4)SSTk-ω模型

湍流SSTk-ω模型[15](簡稱SST模型)綜合了k-ω模型在近壁處邊界層內(nèi)部的模擬和k-ε模型在外部區(qū)域高雷諾數(shù)流動模擬的優(yōu)點，比標(biāo)準(zhǔn)k-ω模型應(yīng)用更為廣泛。

2.2.2 氣液兩相流場模型

分散氣泡流模型的理論依據(jù)是歐拉-歐拉模型，主要求解3個方程[16- 17]，其中流體動量方程為

Φlρlg+F

(6)

混合物連續(xù)性方程為

(7)

氣相輸運方程為

(8)

式中：下標(biāo)為l和g表示液相和氣相；Φ為體積分?jǐn)?shù)；mgl表示液相到氣相的傳質(zhì)速率。

氣泡流速度為

(9)

陰極產(chǎn)生氫氣的質(zhì)量通量為

(10)

式中,NH2為氫氣質(zhì)量通量，M為氫氣的摩爾質(zhì)量，i為局部電流密度，η為電流效率，F(xiàn)為法拉第常數(shù)。

2.3 溫度場

工具陰極和工件陽極的電阻很小，可忽略其加工中所產(chǎn)生的焦耳熱，電解加工過程中產(chǎn)生的熱量主要由電解液焦耳熱和電解液/電極邊界處的電化學(xué)反應(yīng)熱組成。電解液溫度分布可由流體對流擴散方程描述[10- 11]：

(11)

Qbulk=Qj+Qr

(12)

式中,cp為電解液常壓比熱容,kt為電解液熱導(dǎo)率,T為電解液的溫度,Qbulk為電解加工產(chǎn)生的熱量,Qj為焦耳熱,Qr為電化學(xué)反應(yīng)熱。

電解液中產(chǎn)生的焦耳熱為

Qj=κ(?φ)2

(13)

工具陰極和工具陽極與電解液的邊界電化學(xué)反應(yīng)熱為

Qr=Ukik

(14)

式中：Uk為過電位，ik為雙電層的電流密度。

假設(shè)工具陰極和陽極的邊界Γ4和Γ3與空氣自然對流，其余邊界均為熱絕緣。

2.4 多物理場耦合分析

在電解加工中，陽極的溶解速率直接取決于加工間隙內(nèi)各點法向電流密度。當(dāng)施加的外電壓一定時，間隙電場分布受到電解液電導(dǎo)率和幾何結(jié)構(gòu)的影響，電解液的電導(dǎo)率受溫度和氣泡的共同影響，而溫度分布又受到兩相流場和電場等的作用。因此，電解加工過程間隙的溫度分布是多物理場耦合的結(jié)果[12]。

在多場耦合仿真中，通過流體傳熱的對流擴散方程中的對流項來耦合流場和溫度場的影響；通過熱源項(電解液焦耳熱和電化學(xué)反應(yīng)熱)來耦合電場的影響。在模型中可通過設(shè)置陰極/電解液邊界的氣體通量公式實現(xiàn)電場和氣相流場的耦合；通過設(shè)置電解液電導(dǎo)率公式來實現(xiàn)電場、溫度場和兩相流場的耦合。溫度多場耦合仿真流程如圖2所示。

圖2 溫度多場耦合仿真流程圖

2.5 仿真參數(shù)

基于多物理場有限元仿真軟件COMSOL Multiphysics進行求解。將幾何模型導(dǎo)入仿真軟件中，添加多相流場、電場和流體傳熱等模塊，設(shè)置邊界條件并進行網(wǎng)格劃分，單元尺寸選擇極細化。工具陰極、工件陽極和電解液的物理特性設(shè)置如表1所示。

表1 工具陰極、工件陽極和電解液的物理特性

3 仿真分析

近壁區(qū)流速的求解精度影響著電流密度分布的求解精度，流場分布最終影響著溫度場分布。近壁區(qū)流場的求解通常有兩種方法：一種是使用壁函數(shù)將壁面上的物理量與湍流核心區(qū)內(nèi)的相應(yīng)物理量聯(lián)系起來，能大大提高計算速度和收斂性，但由于簡化了邊界層內(nèi)流動的計算，未模擬緩沖區(qū)中的流動，因此該方法對邊界層流場計算精度不高，但在工程中是很實用，應(yīng)用廣泛；另一種是使用低雷諾數(shù)壁處理模型對近壁區(qū)黏性底層和過渡層進行求解，精度較高，但計算成本也相對較高。

3.1 不同湍流模型的間隙流動分析

設(shè)置入口流速為1 m/s，出口壓力為0 MPa(相對于1個大氣壓)，分別在湍流SA、k-ε、k-ω、SST和低雷諾數(shù)k-ε等5個模型中計算間隙內(nèi)流速分布。在仿真中，求解近壁區(qū)流動時，k-ε模型使用壁函數(shù)，另4個模型使用自動壁處理和低雷諾數(shù)壁處理，求得O1O2的流速分布如圖3所示。

從圖3(a)和圖3(b)中可以看出，5種流場模型采用不同壁處理方式時仿真得到的流速分布有區(qū)別，尤其是在近壁區(qū)。在近壁區(qū)的流速值，k-ε模型的仿真值最高，而SA、SST和低雷諾數(shù)k-ε模型的仿真值則很接近。在湍流中心區(qū)的流速值，k-ε模型和k-ω模型比另3個模型要低。從圖3(a)可知，采用自動壁處理的近壁區(qū)流速仿真值差別較

(a)自動壁處理

(b)低雷諾數(shù)壁處理

大，其中k-ω模型為4 m/s，SA、SST和低雷諾數(shù)k-ε模型為1.2 m/s，k-ε模型為11 m/s。從圖3(b)中可知，除k-ε模型以外，其它4個采用低雷諾數(shù)壁求解的速度在邊界層內(nèi)均是從0 m/s開始，一直到湍流中心區(qū)均有求解，其中SA、SST和低雷諾數(shù)k-ε模型計算結(jié)果相近。

造成圖3(a)和圖3(b)中的模型求解結(jié)果的差異是因為選擇自動壁處理時，仿真軟件會根據(jù)網(wǎng)格精細程度自動選擇壁函數(shù)或低雷諾數(shù)壁處理，當(dāng)網(wǎng)格足夠精細時自動調(diào)用求解精度高的低雷諾數(shù)壁處理。k-ε模型是高雷諾數(shù)模型，只對充分發(fā)展的湍流有效，對流動情況變化大的近壁區(qū)黏性底層和緩沖過渡層不進行精確求解，只使用壁函數(shù)進行簡化計算。

SA、SST、k-ω和低雷諾數(shù)k-ε等4個模型是低雷諾數(shù)模型，均可對邊界黏性底層和緩沖過渡層進行求解，在近壁區(qū)的計算精度相對較高[18]。這是因為：SA模型增加了一個額外的無衰減運動學(xué)渦流粘度變量，它可求解實體壁之內(nèi)的整個流場；ω模型是近壁湍流模型，求解的是動能耗散的具體速率，對近壁區(qū)邊界黏性底層低雷諾數(shù)流動模擬較好；SST模型在近壁處采用k-ω模型對近壁區(qū)邊界層流速較低的黏性底層和緩沖過渡層進行計算；低雷諾數(shù)k-ε模型是對標(biāo)準(zhǔn)k-ε模型適應(yīng)于求解低雷諾數(shù)流動而進行修正的模型。在本模型中劃分網(wǎng)絡(luò)時尺寸選擇的是極細化，因而使用自動壁和低雷諾數(shù)壁處理的仿真結(jié)果相差很小。

對于電解加工過程的流體傳熱研究而言，提高近壁區(qū)的流場求解精度無疑將提高溫度場分布求解結(jié)果的可信度。SA模型常用于空氣動力學(xué)流動分析，對墻壁束縛流動模擬效果較好，而對電解加工中流動尺度變換比較大的情況并不太適合；k-ω模型非線性程度大，對初值敏感且難以收斂。SST模型結(jié)合了k-ε模型和k-ω模型的優(yōu)點，求解精度高，求解成本不高，其與低雷諾數(shù)k-ε模型一樣是電解加工仿真中值得研究的模型。

3.2 間隙溫度演變分析

選擇低雷諾數(shù)k-ε湍流模型，采用低雷諾數(shù)壁處理，進行溫度多場耦合仿真，設(shè)置入口流速為1 m/s，出口壓力為0 MPa，加工電壓為20 V。在0.1 s時，模型的溫度場分布如圖4所示。間隙內(nèi)各測量點在0～0.1 s時間段內(nèi)的溫升變化如圖5所示，為更清晰展示溫度隨時間的變化規(guī)律，采用對數(shù)橫坐標(biāo)。

圖4 溫度分布

從圖4可以看出，熱量在電解液和邊界處產(chǎn)生，大部分熱量由電解液帶離加工間隙，在靠近出口處溫升最大，部分熱量在電極內(nèi)進行熱傳導(dǎo)。

由圖5可知，在同一時間點上，加工間隙內(nèi)沿程1～6點的溫度均呈上升態(tài)勢。隨著加工時間的增加，各點溫度急劇上升后又急速下降，經(jīng)過一段時間后逐漸趨于穩(wěn)定。在加工剛開始時，電流密度很大，瞬間溫升很快，在流場的作用下，溫度沿程累積(圖5中點6較為明顯)，各點在0.02 s后溫升值逐漸穩(wěn)定。若記進入穩(wěn)態(tài)溫升的±(2%～5%)范圍內(nèi)為進入準(zhǔn)穩(wěn)態(tài)，在當(dāng)前加工參數(shù)下，點6在0.05 s左右進入準(zhǔn)穩(wěn)態(tài)，其余各點的時間更短。因此，后文可取0.05 s時的物理場分布進行分析。

圖5 點1～6的溫升變化

3.3 不同湍流模型下間隙流動傳熱分析

基于不同湍流模型(k-ε、SST和低雷諾數(shù)k-ε模型)，以及耦合氣泡的低雷諾數(shù)k-ε模型對間隙溫度進行多場耦合仿真。設(shè)置入口流速為1 m/s，出口壓力為0 MPa，加工電壓為20 V。在0.05 s時，點1～6的溫升值如圖6所示。

圖6 不同流場模型的溫度分布

由圖6可以看出，不同流場模型對間隙溫度分布求解結(jié)果不盡相同，沿程各點溫度值呈上升趨勢。SST和低雷諾數(shù)k-ε模型求解值近似相等，且高于k-ε模型求解值，這是因為前兩個模型的計算流速值近似相等且低于后者(如圖3(b)所示)，電解液流速大能在相同時間內(nèi)帶走更多的熱量，因而溫升相對較小。根據(jù)式(3)可知，電解液中氣泡的存在會造成電解液電導(dǎo)率、電流密度的下降，因此耦合氣泡的低雷諾數(shù)k-ε模型求解的沿程各點溫升會有所減小。

3.4 不同加工電壓和入口流速時的溫度分布

選擇耦合氣泡的低雷諾數(shù)k-ε模型進行電解加工溫度多場耦合仿真。圖7為加工電壓分別為16、20和24 V時的沿程溫度分布。圖8為入口流速分別為0.6、1.0和1.4 m/s時的沿程溫度分布。

圖7 不同加工電壓時沿程溫度分布

圖8 不同入口流速時沿程溫度分布

從圖7和圖8中可以看出，在電解加工中，加工電壓越大，單位時間內(nèi)產(chǎn)生的熱量越多，在電解液的帶動下，溫度沿程累積，在出口處最大溫升分別為16、26和40 K左右。增加入口電解液流速或壓力，可使得加工間隙內(nèi)電解液的流速增加，帶走熱量的能力增強，因而沿程溫升變小。

4 試驗驗證

為更進一步驗證前文溫度模型求解的精度以及不同模型求解的差異性，與相關(guān)試驗結(jié)果進行對比分析。試驗數(shù)據(jù)來源見文獻[19]。仿真中，調(diào)整幾何模型與文獻所述的試驗裝置一致，流場模型選擇低雷諾數(shù)k-ε模型(低雷諾數(shù)壁處理)、k-ε模型(壁函數(shù)壁處理)，以及耦合氣泡的低雷諾數(shù)k-ε模型等3種模型進行溫度多場耦合仿真。仿真參數(shù)的設(shè)置與試驗值一致：間隙入口質(zhì)量流速為11.3×10-5m3/s，出口背壓為0.17 MPa，加工電流密度為16 A/cm2，電解液質(zhì)量分?jǐn)?shù)為10%的NaCl溶液，初始溫度為293 K，電導(dǎo)率為13 S/m，陰極進給速度為0.35 mm/min，加工時間為3 min。1～6點仿真溫升數(shù)據(jù)與試驗數(shù)據(jù)對比如圖9所示。

從圖9可以看出，仿真值與試驗值的變化趨勢一致，其中低雷諾數(shù)k-ε模型仿真值比k-ε模型低，這是因為前者在測量點處的流速相對較高。還可以看出，耦合氣泡的低雷諾數(shù)k-ε模型的仿真值整體上與試驗值更為接近，在出口處的溫度近似相等，但仿真和試驗值的差值在流程的前半程比后半程大，其主要原因是試驗裝置中測溫傳感器的安裝是伸出陰極表面的，在間隙入口處氫氣泡在測溫頭附近更容易累積，使得測溫頭附近氣液混合流體的電阻率更大，以致測得的溫度會比實際值高，試驗誤差較大，而在后半程由于氣泡的充分累積，試驗誤差變小。

5 結(jié)論

(1)在流場模型中，采用低雷諾數(shù)模型以及低雷諾數(shù)壁處理能對邊界層流動更精確地求解，且求解精度高于k-ε模型。SST模型和低雷諾數(shù)k-ε模型求解的近壁區(qū)流速值相近，基于該模型的溫度分布也近似相同。

(2)電解加工開始后，間隙內(nèi)的溫度場分布會在一個較短時間內(nèi)進入一個準(zhǔn)穩(wěn)態(tài)。

(3)基于耦合氣泡的低雷諾數(shù)k-ε模型的溫度多場耦合模型對電解加工間隙溫度分布的求解精度高，與溫度分布的實際情況更接近。