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

萬　能　杜　珂　高瑾宇　陳　濤

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

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

萬能杜珂高瑾宇陳濤

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

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

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

0　引言

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

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

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

(1)

工件陽極邊界Γa條件為

φΓa=U

(2)

陰極邊界Γc條件為

φΓc=0

(3)

在邊界Γb上邊界條件為

(4)

式中,φ為電場中各點電勢值;U為工件陽極表面電勢值;n為工件陽極表面各處的法向坐標[8](圖1)。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

采用加權(quán)余量法推導式(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)

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

(9)

Ki j=∫ΩNi·NjdΩ

(10)

Fi=0

(11)

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

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

假設葉片電解加工加工區(qū)間對應U、V、W方向上參數(shù)域的節(jié)點矢量為(0，ε1,ε2,…,εi,εi+1,…,1)、(0,η1,η2,…,ηj,ηj+1,…,1)和(0,ζ1,ζ2,…,ζk,ζk+1,…,1),基于這些節(jié)點矢量所構(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)中的積分運算就不用在整個參數(shù)域Ω內(nèi)進行。考慮參數(shù)域中所有測度不為零的間隔：

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

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

定義單元剛度矩陣Ke和單元載荷向量Fe：

Ke=∫Ωe

Fe=0

(12)

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

4.2加工間隙建模

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

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

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

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

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

u0=u1=…=up=0

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

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

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

其中,Pijk為邊界上的采樣點。

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

(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)最終得到葉片電解加工間隙體的控制頂點為

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

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

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

任意一點處的法矢可記為

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

4.3Dirichlet邊界約束處理

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

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

(13)

5　實驗驗證

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

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

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

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

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

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

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

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

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

6　結(jié)語

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

[1]李志永,朱荻,孫春華,等.發(fā)動機葉片電解加工陰極設計有限元數(shù)值解法研究[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]王蕾.葉片電解加工成形過程仿真與預測研究[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]朱棟,朱荻,徐正揚.航空發(fā)動機葉片電解加工陰極數(shù)字化修正模型及其試驗研究[J].機械工程學報,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].計算機集成制造系統(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邊界條件的配點施加方法[J].機械工程學報,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

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

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

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