萬　能　杜　珂　高瑾宇　陳　濤

西北工業(yè)大學(xué)現(xiàn)代設(shè)計(jì)與集成制造技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，西安，710072

面向葉片電解加工分析的等幾何方法

萬能杜珂高瑾宇陳濤

西北工業(yè)大學(xué)現(xiàn)代設(shè)計(jì)與集成制造技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，西安，710072

為了解決傳統(tǒng)數(shù)值分析方法在電解加工間隙間電場(chǎng)分布邊界敏感問題上計(jì)算精度不足的問題,提出利用等幾何法原理提高葉片電解加工數(shù)值仿真精度的思路,采用NURBS基函數(shù)替代原有拉格朗日基函數(shù)建立加工間隙物理場(chǎng)的求解方程組,解決由于NURBS基函數(shù)在邊界處非插值特性引起的Dirichlet邊界條件施加誤差問題。并通過實(shí)驗(yàn)證明，等幾何方法能提高葉片電解加工數(shù)值分析的計(jì)算精度和收斂速度。

等幾何法;葉片電解加工;NURBS基函數(shù);非插值特性

0　引言

電解工藝加工是對(duì)復(fù)雜型面(例如航空發(fā)動(dòng)機(jī)葉片)加工的一種重要方法,對(duì)這類復(fù)雜工藝的預(yù)測(cè)往往采用數(shù)值仿真方法。傳統(tǒng)的數(shù)值分析方法包括有限元法、有限差分法、邊界元法等,南京航空航天大學(xué)朱荻院士團(tuán)隊(duì)長(zhǎng)期從事葉片類電解加工的研究,采用有限單元法仿真加工間隙區(qū)域的加工狀態(tài)[1-3]。Pattavanitch等[4]利用邊界元法建立了電解加工過程模型,Marius等[5]利用邊界元法分析了陽極工件的腐蝕狀態(tài),Bieniasz[6]采用有限差分法建立了電化學(xué)動(dòng)力學(xué)仿真模型。而這些方法都存在共同的不足:①分析模型準(zhǔn)備時(shí)間長(zhǎng),常出現(xiàn)網(wǎng)格劃分質(zhì)量不好的現(xiàn)象;②仿真結(jié)果表現(xiàn)為單元網(wǎng)格節(jié)點(diǎn)的位置變化,需要通過擬合網(wǎng)格節(jié)點(diǎn)重構(gòu)復(fù)雜型面;③采用多項(xiàng)式基函數(shù)的網(wǎng)格單元逼近表示邊界,從原理上不能精確表達(dá)求解區(qū)域邊界處的約束,不適于流體分析這類對(duì)邊界敏感問題的求解。因此本文借用近年來新興的等幾何思想建立統(tǒng)一的幾何建模與分析仿真的數(shù)學(xué)模型,研究基于NURBS基函數(shù)的電解加工間隙幾何參數(shù)化方法與非插值特性的邊界條件施加方法,形成支持葉片電解加工這類具有復(fù)雜敏感邊界問題的數(shù)值求解方法[7]。

1　電解加工的數(shù)學(xué)模型

一般認(rèn)為,當(dāng)電解加工過程處于平衡狀態(tài)時(shí),加工間隙內(nèi)的電場(chǎng)屬于穩(wěn)恒電場(chǎng)，電位分布符合拉普拉斯方程:

(1)

工件陽極邊界Γa條件為

φΓa=U

(2)

陰極邊界Γc條件為

φΓc=0

(3)

在邊界Γb上邊界條件為

(4)

式中,φ為電場(chǎng)中各點(diǎn)電勢(shì)值;U為工件陽極表面電勢(shì)值;n為工件陽極表面各處的法向坐標(biāo)[8](圖1)。

圖1　葉片電解加工間隙的平面模型

2　等幾何法的原理及基函數(shù)特點(diǎn)

2.1基函數(shù)統(tǒng)一的幾何與分析模型

在葉片電解加工的計(jì)算機(jī)輔助分析中，葉片工件與陰極工具之間的加工間隙一般采用CAD系統(tǒng)建立其幾何模型。目前商用CAD系統(tǒng)通常采用非均勻有理B樣條(non-uniformrationalB-splines,NURBS)表示葉片這類具有復(fù)雜曲面的幾何模型。等幾何法的基本思想就是采用同一套NURBS基函數(shù)統(tǒng)一表達(dá)幾何模型和數(shù)值分析模型。對(duì)于加工間隙采用樣條體(NURBSvolumes)表示,它是使用三個(gè)節(jié)點(diǎn)矢量定義的張量積樣條:

i≠i′j≠j′k≠k′

式中,{I,J,K}為節(jié)點(diǎn)矢量的索引集;{Ri,p,Rj,q,Rk,r}為各節(jié)點(diǎn)矢量對(duì)應(yīng)的單變量B樣條基函數(shù);ωijk、ωi′j′k′為權(quán)因子；p、q、r為基函數(shù)的次數(shù)。

2.2自然劃分加工間隙的參數(shù)域

NURBS基函數(shù)是由節(jié)點(diǎn)矢量和樣條次數(shù)定義的，而張量積樣條的節(jié)點(diǎn)矢量正好張成加工間隙參數(shù)域上的規(guī)則網(wǎng)格。等幾何法即采用節(jié)點(diǎn)矢量張成的規(guī)則網(wǎng)格做自然劃分,無需后續(xù)的網(wǎng)格剖分工作。借用經(jīng)典有限元法中單元和節(jié)點(diǎn)的概念，把等幾何法中的單元類比為測(cè)度不等于零的節(jié)點(diǎn)間隙，而節(jié)點(diǎn)類比為間隙內(nèi)非零基函數(shù)對(duì)應(yīng)的控制頂點(diǎn)。對(duì)于三維張量積樣條，單元可用直積表示為

Ωe=[εi,εi+1]?[ηj,ηj+1]?[ζk,ζk+1]

εi<εi+1ηj<ηj+1ζk<ζk+1

其中,εi、ηj、ζk表示三個(gè)參數(shù)方向上第i、j、k個(gè)參數(shù)值?？梢园l(fā)現(xiàn)單元內(nèi)的非零基函數(shù)共有(p+1)×(q+1)×(r+1)個(gè)。

2.3NURBS基函數(shù)的非插值性

相對(duì)于經(jīng)典有限元的多項(xiàng)式基函數(shù),NURBS基函數(shù)具有很多優(yōu)點(diǎn),例如它可以精確表示任意的幾何模型,在單元邊界處可以獲得更高的連續(xù)性,但它缺少一個(gè)重要的性質(zhì)即在節(jié)點(diǎn)處的插值性,即Ni(εj)≠δi j(εj為節(jié)點(diǎn)處的參數(shù)值)。單變量基函數(shù)在首末端點(diǎn)處滿足插值條件，但是二維張量積樣條基函數(shù)除四個(gè)角點(diǎn)外，在其他各節(jié)點(diǎn)處都不具有插值性。因此等幾何法不能像傳統(tǒng)有限元法那樣對(duì)葉片電解加工間隙的節(jié)點(diǎn)處場(chǎng)變量進(jìn)行插值以表示Dirichlet邊界條件[9]。

3　葉片電解加工的常微分方程組形式

采用加權(quán)余量法推導(dǎo)式(1)的等效積分弱形式。式(1)兩邊同乘以權(quán)函數(shù)w,得到：

∫ΩwΔφdΩ=0

?Ω=Γa+Γb+Γc

由于在邊界Γb上,在本質(zhì)邊界Γa和Γc上權(quán)函數(shù)w=0,由格林第一公式

∫Ωw·φdΩ-∮?Ωw(φ·n)ds=0

可得到

∫Ωw·φdΩ=0

(5)

上述等效積分形式可以寫成與特定問題無關(guān)的一般形式：

a(u,v)=l(v)

(6)

u∈S,v∈V,S={u:u∈H(Ω),u|Γa=U,u|Γc=0}

V={v:v∈H(Ω),v|?Ω=0}

其中，a(u,v)和l(v)分別為定義在區(qū)間Ω上的雙線性和線性泛函。

(7)

取權(quán)函數(shù)w為NURBS基函數(shù)族,再代入到式(6)中可以得到方程組形式如下:

∫ΩNjφj)dΩ=0i=1,2,…,n

(8)

整理得到電位場(chǎng)值問題的一階常微分方程組形式如下：

(9)

Ki j=∫ΩNi·NjdΩ

(10)

Fi=0

(11)

4　程序?qū)崿F(xiàn)

4.1剛度矩陣和載荷向量的裝配

假設(shè)葉片電解加工加工區(qū)間對(duì)應(yīng)U、V、W方向上參數(shù)域的節(jié)點(diǎn)矢量為(0，ε1,ε2,…,εi,εi+1,…,1)、(0,η1,η2,…,ηj,ηj+1,…,1)和(0,ζ1,ζ2,…,ζk,ζk+1,…,1),基于這些節(jié)點(diǎn)矢量所構(gòu)建的基函數(shù)為Ni、Nj、Nk。由于NURBS基函數(shù)的局部支撐性,即基函數(shù)Ni、Nj、Nk只在區(qū)間[εi，εi+p+1]、[ηj,ηj+p+1]、[ζk,ζk+p+1]內(nèi)有非零值,其中基函數(shù)取工程中常見的3次函數(shù),即p=3,式(10)中的積分運(yùn)算就不用在整個(gè)參數(shù)域Ω內(nèi)進(jìn)行?？紤]參數(shù)域中所有測(cè)度不為零的間隔：

Ωe=[εi,εi+1]?[ηj,ηi+1]?[ζk,ζk+1]

可以把Ωe看成為電解加工間隙劃分的等幾何分析單元,因此顯然有結(jié)論：

定義單元?jiǎng)偠染仃嘖e和單元載荷向量Fe：

Ke=∫Ωe

Fe=0

(12)

類似于有限元方法,等幾何分析法也可以看成是劃分了NURBS樣條體單元,但這種單元在加工間隙建模完成時(shí)即已完成。同樣,等幾何法也有一個(gè)單元?jiǎng)偠染仃嚭洼d荷向量的裝配過程，通過單元?jiǎng)偠染仃囇b配得到全局剛度矩陣。等幾何單元和經(jīng)典有限元單元的區(qū)別在于：①基函數(shù)不具備插值性質(zhì)，節(jié)點(diǎn)(控制定點(diǎn))有可能不在單元區(qū)域上；②等幾何單元有更高的單元邊界連續(xù)性[10]。

4.2加工間隙建模

葉片電解加工間隙的邊界Γa與邊界Γc為自由曲面,邊界Γb為平面。為了建立加工間隙的參數(shù)化幾何模型,可以對(duì)工件表面和陰極工具表面進(jìn)行采樣,獲得邊界Γa和Γc上的采樣點(diǎn),另外依據(jù)等參條件給定Γb上的采樣點(diǎn)。

建立葉片電解加工間隙參數(shù)化模型(圖2)的步驟如下：

圖2　葉片電解加工間隙的三維參數(shù)化模型

(1)設(shè)加工間隙的a、b、c邊上分別存在U、V、W方向上的采樣點(diǎn)Oi、Pj和Qk。其中，i=0,1,…,m；j=0,1,…,n;k=0,1,…,l。

(3)構(gòu)建三個(gè)方向的節(jié)點(diǎn)矢量，即令節(jié)點(diǎn)矢量中的首末參數(shù)值需要滿足:

u0=u1=…=up=0

um=um+1=…=um+p=1

(4)三個(gè)方向優(yōu)化后的節(jié)點(diǎn)矢量為ui′、vj′和wk′,其中,i′=0,1,…,m+p+1;j′=0,1,…,n+p+1;k′=0,1,…,l+p+1。以u(píng)i′、vj′和wk′作為節(jié)點(diǎn)矢量構(gòu)建邊界Γ在U、V、W方向上的NURBS基函數(shù)Ni、Nj和Nk。

(5)分別對(duì)Γa、ΓbF、ΓbB、ΓbL、ΓbR、Γb6個(gè)邊界面反求其控制頂點(diǎn)。若記加工間隙邊界上的控制頂點(diǎn)為Cctrl(Γ),則可得到：

其中,Pijk為邊界上的采樣點(diǎn)。

(6)通過對(duì)控制頂點(diǎn)Vij0、Vij1、Vi0k、Vi1k、V0jk、V1jk超限插值可以得到整個(gè)加工間隙體的控制頂點(diǎn)：

(1-εi)ζkV0jl(0,ηj,1)+εi(1-ζk)Vmj0(1,ηj,0)+

εiζk·Vmjl(1,ηj,1)]

0,ζk)+(1-εi)ηjV0nk(0,1,ζk)+

εi(1-ηj)Vm0k(1,0,ζk)+εiηjVmnk(1,1,ζk)]

(1-ηj)ζkVi0l(εi,0,1)+

ηj(1-ζk)Vin0(εi,1,0)+ηjζkVinl(εi,1,1)]

V000(0,0,0)+(1-εi)(1-ηj)ζkV00l(0,0,1)+

(1-εi)ηj(1-ζk)V0n0(0,1,0)+εi(1-ηj)(1-

ζk)Vm00(1,0,0)+(1-εi)ηjζkV0nl(0,1,1)+

εi(1-ηj)ζkVm0l(1,0,1)+εiηj(1-ζk)Vmn0(1,

1,0)+εiηjζkVmnl(1,1,1)]

(7)最終得到葉片電解加工間隙體的控制頂點(diǎn)為

Vijk=U+V+W-UW-UV-VW+UVW

其中,εi、ηj、ζk∈[0,1]。因此，葉片電解加工間隙的幾何模型可表示為

(8)陽極型面的表達(dá)式為

任意一點(diǎn)處的法矢可記為

則式(11)中的cosθ=a·b/(|a|·|b|)。

4.3Dirichlet邊界約束處理

由于采用NURBS基函數(shù)表達(dá)葉片電解加工間隙的電場(chǎng)分布,因此不能像傳統(tǒng)分段多項(xiàng)式有限元單元一樣插值Dirichlet邊界上采樣點(diǎn)的電場(chǎng)值，所以采用強(qiáng)施加方法對(duì)式(9)施加Dirichlet邊界條件[11-12]。假設(shè)弱解φ∈S可以表示為兩部分之和,即u=e+g,其中g(shù)∈S,e∈V,代入到式(6)中得到：

a(e,v)=l(v)-a(g,v)

(13)

5　實(shí)驗(yàn)驗(yàn)證

設(shè)置電解加工仿真的電勢(shì)差為15V,電解加工初始間隙設(shè)置為0.5mm,陰極進(jìn)給速度為0.5mm/min。電解液成分為NaNO3,質(zhì)量分?jǐn)?shù)為10%,初始溫度為25～30 ℃,流速為15m/s,工件材料為2Cr13鋼。

利用三維產(chǎn)品設(shè)計(jì)平臺(tái)NX建立葉片電解加工間隙的幾何模型,如圖3a所示。反求出加工間隙幾何的控制頂點(diǎn)，并利用超限插值得到加工間隙的參數(shù)化模型，如圖3b所示。

(a)加工間隙的幾何模型(b)加工間隙參數(shù)化模型的控制頂點(diǎn)圖3　葉片電解加工間隙幾何模型與參數(shù)化模型

利用邊界配點(diǎn)法施加陰陽極邊界電勢(shì)約束條件后,得到等幾何方法分析的電勢(shì)分布、電場(chǎng)矢量分布，并分別與經(jīng)典有限元法進(jìn)行比較,如圖4所示。

(a)等幾何分析電勢(shì)分布云圖(b)有限元分析電勢(shì)分布云圖

(c)等幾何分析電場(chǎng)矢量圖(d)有限元分析電場(chǎng)矢量圖圖4　葉片電解加工間隙的電勢(shì)分析比較

采用經(jīng)典有限元和等幾何分析對(duì)葉片電解加工間隙的收斂速度進(jìn)行比較，經(jīng)典有限元法采用線性、二次、三次拉格朗日單元(p=1，2，3),而等幾何分析采用了二次和三次樣條函數(shù),得到曲線如圖5所示。圖5中，fDOF為自由度，e為分析誤差。

圖5　葉片電解加工間隙的收斂速度比較

雖然樣條函數(shù)也可理解為定義在參數(shù)域內(nèi)的分段(有理)多項(xiàng)式,但它通?？梢垣@得比經(jīng)典有限元更高的單元邊界連續(xù)性。經(jīng)典有限元在單元邊界處通常是C0連續(xù),而樣條函數(shù)可以獲得Cp-r(r為節(jié)點(diǎn)重復(fù)次數(shù))連續(xù)。因此,從理論上講,采用NURBS樣條基函數(shù)能夠更精確地表達(dá)葉片這類具有復(fù)雜自由曲面的邊界幾何形狀，從而等幾何分析方法可以獲得更高的分析精度。另一方面,在相同的網(wǎng)格自由度情況下,FEM和等幾何法都能達(dá)到收斂，但由圖5中相同自由度情況下的誤差比較可發(fā)現(xiàn)，等幾何法的誤差較傳統(tǒng)有限元法精度更高。相對(duì)來說等幾何法的收斂速度明顯要快于FEM,那么在相同網(wǎng)格單元數(shù)量下,等幾何法體現(xiàn)出了更高的分析精度。

6　結(jié)語

在傳統(tǒng)基于有限元法的葉片電解加工數(shù)值分析中，幾何模型與分析模型所采用的數(shù)學(xué)描述方法不同，兩者之間需要相互轉(zhuǎn)換，帶來分析模型的準(zhǔn)備時(shí)間長(zhǎng)，轉(zhuǎn)換常出現(xiàn)模型質(zhì)量不高的缺點(diǎn)。另外由于有限元法采用多項(xiàng)式基函數(shù)網(wǎng)格單元逼近表示邊界，所以原理上不能精確表達(dá)電極邊界與電解液邊界的場(chǎng)變量分布。因此本文采用NURBS基函數(shù)取代多項(xiàng)式基函數(shù)，實(shí)現(xiàn)葉片電解加工間隙幾何建模與數(shù)值分析共用相同的基函數(shù)，即建模完成同時(shí)網(wǎng)格劃分完成。同時(shí)利用NURBS基函數(shù)構(gòu)建的樣條體單元能準(zhǔn)確表達(dá)加工間隙自由曲面邊界的特性，實(shí)現(xiàn)更加精確的加工間隙內(nèi)電場(chǎng)分析。最終通過實(shí)例驗(yàn)證了等幾何法對(duì)葉片電解加工數(shù)值分析這類具有復(fù)雜幾何邊界問題的有效性。

[1]李志永,朱荻,孫春華,等.發(fā)動(dòng)機(jī)葉片電解加工陰極設(shè)計(jì)有限元數(shù)值解法研究[J].中國機(jī)械工程,2004,15(13):1151-1154.

Li Zhiyong,Zhu Di,Sun Chunhua,et al.Study on Finite-element Arithmetic in Electrochemical Machining for Turbine Blades[J].China Mechanical Engineering,2004,15(13):1151-1154.

[2]王蕾.葉片電解加工成形過程仿真與預(yù)測(cè)研究[J].機(jī)械設(shè)計(jì)與制造,2009,29(11):101-103.Wang Lei.Study on Simulation and Prediction of Blade’s Shape in Electrochemical Machining[J].Machinery Design & Manufacture,2009,29(11):101-103.

[3]朱棟,朱荻,徐正揚(yáng).航空發(fā)動(dòng)機(jī)葉片電解加工陰極數(shù)字化修正模型及其試驗(yàn)研究[J].機(jī)械工程學(xué)報(bào),2011,47(7)：191-198.

Zhu Dong,Zhu Di,Xu Zhengyang.Experimental Study on the Catode Digital Modification of Turbine Blade in Electrochemical Machining[J].Journal of Mechanical Engineering,2011,47(7)：191-198.

[4]Pattavanitch J,Hinduja S,Atkinson J.Modelling of the Electrochemical Machining Process by the Boundary Element Method[J].CIRP Annals-manufacturing Technology,2010,59:243-246.

[5]Marius P,Leslie B.3D Electrochemical Machining Computer Simulations[J].Journal of Materials Processing Technology,2004,149:472-478.

[6]Bieniasz L K.Finite-difference Electrochemical Kinetic Simulations Using the Rosenbrock Time Integration Scheme[J].Journal of Electroanalytical Chemistry,1999,469:97-115.

[7]Xiang Y,Mo R,Wan N,et al.The High Precision Blade Electrochemical Machining Simulation and Cathode Optimization Based on Isogeometric Method[J].Applied Mechanics and Materials,2013,339:489-494.

[8]Sun C H,Zhu D,Li Z H,et al.Application of FEM to Tool Design for Electrochemical Machining Freeform Surface[J].Finite Elements in Analysis and Design,2006,43:168-172.

[9]Huges T,Cottrell J,Bazilevs Y.Isogeometric Analysis:CAD,Finite Elements,NURBS,Exact Geometry and Mesh Refinement[J].Computer Methods in Applied Mechanics and Engineering,2005,194(39/41):4135-4195.

[10]陳濤,莫蓉,張欣.固體介質(zhì)瞬態(tài)傳熱問題的等幾何分析[J].計(jì)算機(jī)集成制造系統(tǒng),2011,17(9)：1988-1996.

Chen Tao,Mo Rong,Zhang Xin.Isogeometric Analysis for Transient Heat Conduction of Solid Medium[J].Computer Integrated Manufacturing Systems,2011,17(9):1988-1996.

[11]陳濤,莫蓉,萬能.等幾何分析中Dirichlet邊界條件的配點(diǎn)施加方法[J].機(jī)械工程學(xué)報(bào),2012,48(5):1-8.

Chen Tao,Mo Rong,Wan Neng.Imposing Dirichlet Boundary Conditions with Point Collocation Method in Isogeometric Analysis[J].Journal of Mechanical Engineering,2012,48(5):1-8.

[12]Chen T,Mo R,Wan N.NURBS-based Isogeometric Finite Element Method for Analysis of Two-dimensional Piezoelectric Device[J].Procedia Engineering,2011,15:3562-3566.

(編輯王艷麗)

Electrochemical Machining Analysis of Aero-engine Blade Based on Isogeometric Method

Wan NengDu KeGao JingyuChen Tao

The Key Laboratory of Contemporary Design and Integrated Manufacturing Technology,Northwestern Polytechnical University,Xi’an,710072

For solving the low precision problem for sensitive boundary which caused by traditional numerical analysis method,this paper proposed a method of promoting the analysis precision by using the isogeometric analysis.The NURBS basis functions were used to replace the Lagrange basis function for establishing the solving equations of the electrochemical machining gap.The problem of the non-interpolation feature belonged to the NURBS basis functions was settled,which could bring the errors for imposing Dirichlet boundary condition.At last,the superiority,including precision and the rate of convergence,of the isogeometric analysis was proved by the comparison tests.

isogeometric analysis;electrochemical machining;NURBS basis function;non-interpolation feature

2014-06-11

國家自然科學(xué)基金資助項(xiàng)目(51205320);陜西省自然科學(xué)基金資助項(xiàng)目(2012JQ7002)

TP391.72DOI：10.3969/j.issn.1004-132X.2015.10.016

萬能,男,1979年生。西北工業(yè)大學(xué)機(jī)電學(xué)院副教授。主要研究方向?yàn)镃AD/CAPP/CAM。發(fā)表論文10余篇。杜珂,男,1991年生。西北工業(yè)大學(xué)機(jī)電學(xué)院碩士研究生。高瑾宇,男,1990年生。西北工業(yè)大學(xué)機(jī)電學(xué)院碩士研究生。陳濤,男,1981年生。西北工業(yè)大學(xué)機(jī)電學(xué)院博士研究生。