譚先濤，楊斌堂，孟 光，徐彭有，楊德華

(1. 上海交通大學(xué)機(jī)械系統(tǒng)與振動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200240；2.中國(guó)科學(xué)院國(guó)家天文臺(tái)南京天文光學(xué)技術(shù)研究所，南京 210042)

由于大口徑鏡面制造困難，目前世界上的大口徑天文光學(xué)望遠(yuǎn)鏡方案均采用了拼接子鏡主動(dòng)光學(xué)技術(shù)，采用六面體或者扇形子鏡拼接而成，并通過(guò)單個(gè)鏡面的精密控制來(lái)實(shí)現(xiàn)整體鏡面的協(xié)調(diào)一致?？刂茊蝹€(gè)鏡面的微位移驅(qū)動(dòng)器需要具有高精度、大行程以及大負(fù)載的技術(shù)要求。目前廣泛使用的微位移驅(qū)動(dòng)包括電機(jī)驅(qū)動(dòng)、液壓驅(qū)動(dòng)、壓電驅(qū)動(dòng)等驅(qū)動(dòng)方式。電機(jī)驅(qū)動(dòng)的機(jī)械式微位移驅(qū)動(dòng)器存在著間隙、傳動(dòng)誤差、摩擦損耗及爬行等現(xiàn)象，難以到達(dá)高精度的要求；液壓驅(qū)動(dòng)的微位移驅(qū)動(dòng)器分辨率較低、輸出負(fù)載能力較小、反應(yīng)速度較慢；壓電驅(qū)動(dòng)器具有高分辨率大承載等優(yōu)點(diǎn)，但由于壓電應(yīng)變量較小，難以實(shí)現(xiàn)大行程，并且需要高電壓驅(qū)動(dòng)，發(fā)熱嚴(yán)重。超磁致伸縮驅(qū)動(dòng)器(Giant Magnetostrictive Actuator，GMA)具有輸出位移和輸出力大、機(jī)械響應(yīng)速度快等優(yōu)點(diǎn)，在大型天文望遠(yuǎn)鏡拼接子鏡的精密驅(qū)動(dòng)控制方面具有潛在優(yōu)勢(shì)。

GMA的核心部件為超磁致伸縮材料(Terfenol-D)棒。Terfenol-D作為一種稀土超磁致伸縮材料，具有非常復(fù)雜的電磁-結(jié)構(gòu)-熱多物理場(chǎng)耦合特性，并且其材料屬性如相對(duì)磁導(dǎo)率、彈性模量等隨外加磁場(chǎng)、預(yù)壓應(yīng)力以及溫度的變化而變化。建立精確的數(shù)學(xué)模型，并選用有限元計(jì)算方法對(duì)GMA進(jìn)行仿真分析，從而優(yōu)化GMA的結(jié)構(gòu)和電磁性能，對(duì)提高GMA的輸出性能具有重要意義。

針對(duì)GMA的研究，Engdahl、Perez和Kannan[1-3]分別建立了描述其耦合特性的模型，但是都沒(méi)有考慮到材料屬性參數(shù)的非線性變化；Benatar[4]采用FEMLAB計(jì)算了耦合三維有限元模型，但也沒(méi)有考慮參數(shù)的非線性；趙章榮[5]引入Jiles-Atherton磁滯模型建立了三維動(dòng)態(tài)有限元模型，但由于模型的強(qiáng)非線性，三維模型計(jì)算復(fù)雜，不易得到收斂結(jié)果?；谝陨戏治?，并結(jié)合所設(shè)計(jì)驅(qū)動(dòng)器實(shí)際工作環(huán)境，本研究采用Jiles-Atherton模型計(jì)算磁化強(qiáng)度M、磁致伸縮量λ和磁場(chǎng)強(qiáng)度H、預(yù)壓應(yīng)力σ0之間的非線性關(guān)系，建立了二維軸對(duì)稱非線性有限元模型，并嘗試使用商用有限元軟件COMSOL 3.5實(shí)現(xiàn)了整個(gè)驅(qū)動(dòng)器的電磁-結(jié)構(gòu)全耦合分析。

1 超磁致伸縮驅(qū)動(dòng)器二維軸對(duì)稱非線性位移模型

直動(dòng)型超磁致伸縮驅(qū)動(dòng)器的工作原理是由勵(lì)磁線圈通入電流產(chǎn)生電磁場(chǎng)，超磁致伸縮棒(Terfenol-D)在磁場(chǎng)作用下，由于內(nèi)部磁疇的偏轉(zhuǎn)而產(chǎn)生伸縮變化，并由輸出頂桿輸出位移和力，驅(qū)動(dòng)負(fù)載。

基于驅(qū)動(dòng)器的軸對(duì)稱性質(zhì)，采用二維軸對(duì)稱方式建立驅(qū)動(dòng)器建立有限元模型，不僅可以保證求解的精度不下降，同時(shí)可以大大節(jié)省計(jì)算機(jī)資源。如圖1，為超磁致伸縮驅(qū)動(dòng)器的結(jié)構(gòu)模型。

圖1 GMA結(jié)構(gòu)示意圖Fig.1 Structure of a giant magnetostrictive actuator(GMA)

圖2 超磁致伸縮驅(qū)動(dòng)器電磁場(chǎng)和機(jī)械應(yīng)力場(chǎng)耦合關(guān)系Fig.2 Illustration of the couplings between electromagnetic field and mechanical stress field in a GMA

如圖2，驅(qū)動(dòng)器中包含電磁—機(jī)械的雙向耦合。利用有限元方法分析超磁致伸縮驅(qū)動(dòng)器的性能，實(shí)際上就是求解描述耦合物理場(chǎng)偏微分方程邊值問(wèn)題。以下過(guò)程將首先建立電磁場(chǎng)和機(jī)械應(yīng)力場(chǎng)的偏微分方程及其邊界條件，然后引入Jiles-Atherton模型描述Terfenol-D棒內(nèi)的電磁—機(jī)械耦合特性。

1.1 超磁致伸縮驅(qū)動(dòng)器的電磁場(chǎng)

驅(qū)動(dòng)器工作在電磁場(chǎng)中，采用Maxwell方程組描述電磁場(chǎng)特性[6]：

(1)

(2)

▽·D=ρ,

(3)

▽·B=0,

(4)

針對(duì)驅(qū)動(dòng)器的軸對(duì)稱特性，使用磁矢量勢(shì)描述電磁場(chǎng)，則只有圓周方向分量Aφ≠0。

根據(jù)(1)可知：

(5)

對(duì)上式兩邊同時(shí)乘以δAφ(δAφ表示Aφ的變分)，并在整個(gè)域內(nèi)積分，可得積分弱解形式方程為：

▽×δAφ)tHdΩ+

(6)

電磁場(chǎng)方程求解的初始和邊界條件為：

A|t=0=A0， ?ΩA∶A=A*

A|t=0=V0， ?ΩV∶V=V*

(7)

1.2 超磁致伸縮驅(qū)動(dòng)器的機(jī)械應(yīng)力場(chǎng)

根據(jù)牛頓第二定律，Terfenol-D棒內(nèi)部力平衡方程為：

(8)

其中T為應(yīng)力；b為體積力；u=[u,w]T；u、w分別表示徑向和軸向的位移。

應(yīng)力-應(yīng)變關(guān)系根據(jù)胡克定律有：

T=CS

(9)

其中C為剛度矩陣；S為應(yīng)變向量。應(yīng)變-位移關(guān)系為：

(10)

對(duì)式(8)兩端同時(shí)乘以δu，并在整個(gè)域內(nèi)積分，可得積分弱解形式方程：

(11)

應(yīng)力張量T由3部分組成[2]：機(jī)械應(yīng)力分量Tmech；Maxwell應(yīng)力張量Tm；預(yù)壓應(yīng)力張量T0。其中：

(12)

機(jī)械應(yīng)力場(chǎng)方程求解的初始和邊界條件為：

(13)

1.3 GMM內(nèi)的磁-結(jié)構(gòu)耦合場(chǎng)

GMM中磁場(chǎng)H由線圈的源磁場(chǎng)Hs以及磁化產(chǎn)生的磁場(chǎng)HM組成[7]：

H=Hs+HM

(14)

式中HM=-▽?duì)?φ表示簡(jiǎn)化磁標(biāo)量位，則GMM中的磁場(chǎng)強(qiáng)度可表示為：

H=-▽?duì)?Hs

(15)

A E Clark在1980年提出了描述超磁致伸縮效應(yīng)的線性壓磁方程：

S=T/YH+dH,

(16)

B=dT+μTH.

(17)

其中YH、d和μT分別表示Terfenol-D棒的彈性模量矩陣、壓磁系數(shù)矩陣和磁導(dǎo)率矩陣；上標(biāo)H,T分別表示恒定磁場(chǎng)和恒定應(yīng)力條件。

線性壓磁方程是在一定的外加偏置磁場(chǎng)和預(yù)壓力作用下，忽略磁滯非線性特性的基礎(chǔ)上得到的，能夠描述磁致伸縮和逆磁致伸縮效應(yīng)，但不能反映材料的磁滯回特性。

Jiles和Atherton建立了J-A磁滯模型[8-9]，該模型基于微磁學(xué)理論和Weiss磁籌理論，表達(dá)式簡(jiǎn)單，參數(shù)較少。本文采用J-A模型描述Terfenol-D棒的磁化過(guò)程。包含外加可變預(yù)應(yīng)力σ0的J-A模型表達(dá)式為：

(18)

表1 預(yù)壓應(yīng)力為8MPa，最大電流為8A時(shí)通過(guò)遺傳算法辨識(shí)出的J-A模型參數(shù)

壓磁方程中，總的應(yīng)變是由彈性應(yīng)變T/YH和磁致伸縮應(yīng)變dH兩部分構(gòu)成。磁致伸縮應(yīng)變?chǔ)?dH根據(jù)二次疇轉(zhuǎn)模型可寫為：

(19)

式中M可以通過(guò)J-A模型，采用四階Runge-Kutta法計(jì)算求得。計(jì)算得到的M-H,λ-H，曲線如圖3。

圖3 根據(jù)Jiles-Atherton磁滯模型計(jì)算得到的磁化強(qiáng)度、磁致伸縮量和磁場(chǎng)強(qiáng)度的函數(shù)關(guān)系Fig.3 The variations of magnetization and magnetostriction displacement with magnetic field intensity according to the Jiles-Atherton model

將上述獲得的M-H、λ-H關(guān)系寫成代數(shù)插值表達(dá)式，帶入(16)中可以得:

T=EHS-EHλ

(20)

磁致伸縮棒內(nèi)的磁感應(yīng)強(qiáng)度B可以通過(guò)B、H、M之間的關(guān)系獲得：

B=μ0(M+H)

(21)

(20)、(21)即為包含磁滯的非線性本構(gòu)關(guān)系。

根據(jù)以上分析，(6)、(11)分別為電磁場(chǎng)和機(jī)械應(yīng)力場(chǎng)的控制方程；(7)、(13)分別為兩個(gè)物理場(chǎng)的邊界條件。同時(shí)考慮(20)、(21)的本構(gòu)方程式，采用有限元方法進(jìn)行求解。

2 COMSOL Multiphysics 3.5求解及實(shí)驗(yàn)結(jié)果比較

COMSOL Multiphysics是一款大型的高級(jí)數(shù)值仿真軟件。它基于偏微分方程組定義物理模型，具有很好的交互開發(fā)環(huán)境界面。COMSOL Multiphysics成功實(shí)現(xiàn)了任意多物理場(chǎng)、直接、雙向?qū)崟r(shí)耦合。既可以直接選擇封裝的模塊來(lái)模擬某一類型的物理模型，也可以采用PDE(Partial Differential Equations)模塊來(lái)直接定義特殊的物理模型。基于超磁致伸縮材料獨(dú)特的磁-結(jié)構(gòu)非線性耦合特性，在比較了常用的數(shù)值仿真軟件如ANSYS之后，最終選用了COMSOL Multiphysics 3.5作為解決超磁致伸縮材料非線性多物理場(chǎng)耦合問(wèn)題的計(jì)算分析平臺(tái)。

在COMSOL Multiphysics3.5中，首先選擇模型的空間維度為2D軸對(duì)稱，然后設(shè)置3個(gè)物理環(huán)境：2D軸對(duì)稱應(yīng)力-應(yīng)變模型，變量{uor,w,p};2D軸對(duì)稱磁場(chǎng)模型，變量{Aphi}；2D弱解模型，變量{u1,w1,Vm}。輸入幾何模型之后，即可以設(shè)置參數(shù)、表達(dá)式、函數(shù)來(lái)定義物理環(huán)境常量及變量(材料參數(shù)、求解域參數(shù)、邊界設(shè)定等)。對(duì)于應(yīng)力-應(yīng)變模型和磁場(chǎng)模型，可以直接設(shè)置這些參數(shù)，而對(duì)于弱解模型，自定義的弱解表達(dá)式如表2。

表2 自定義的弱解形式方程式

其中，sr、sz、er、ez分別表示徑向和軸向的應(yīng)力、應(yīng)變；Br、Bz、Hr、Hz分別表示徑向和軸向的磁通量密度和磁場(chǎng)強(qiáng)度；rho_rod表示超磁致伸縮棒的密度；_w,_test分別表示弱解模式的名稱和試探函數(shù)。電磁場(chǎng)和應(yīng)力場(chǎng)相互聯(lián)系的本構(gòu)關(guān)系通過(guò)設(shè)置等式表達(dá)式的變量來(lái)實(shí)現(xiàn)。

如前所述，選擇了實(shí)驗(yàn)最佳預(yù)壓應(yīng)力σ0=8MPa，激勵(lì)源為一個(gè)760匝，線徑為1mm的銅線圈，電流施加方式為隨時(shí)間斜線增大到8A。

基于以上分析，對(duì)GMA進(jìn)行全耦合計(jì)算，計(jì)算結(jié)果如圖4～6。由圖4、5可以看出，在電流強(qiáng)度為8A情況下，超磁致伸縮棒內(nèi)部磁場(chǎng)強(qiáng)度大約為60KA/m，應(yīng)變大約為1250ppm，并且在磁致伸縮棒兩端的應(yīng)變要大于棒中間的應(yīng)變，這主要是由于磁致伸縮棒兩端的磁場(chǎng)強(qiáng)度要略大于中心位置的磁場(chǎng)強(qiáng)度。

圖4 超磁致伸縮驅(qū)動(dòng)器磁場(chǎng)強(qiáng)度分布計(jì)算結(jié)果Fig.4 Calculation result of the magnetic field intensity distribution in the GMA

圖5 磁致伸縮棒內(nèi)部應(yīng)變分布計(jì)算結(jié)果Fig.5 Calculation result of the strain distribution in the GMM rod

分別選取磁致伸縮棒幾何中心點(diǎn)和驅(qū)動(dòng)器輸出端面中點(diǎn)，通過(guò)COMSOL Multiphysics 3.5計(jì)算整個(gè)時(shí)間歷程的磁致伸縮棒內(nèi)的磁場(chǎng)強(qiáng)度和驅(qū)動(dòng)器的位移量，計(jì)算結(jié)果和實(shí)驗(yàn)結(jié)果的對(duì)比如圖6。

圖6 (a)激勵(lì)磁場(chǎng)強(qiáng)度和輸入電流強(qiáng)度之間的關(guān)系；(b)輸出位移和輸入電流強(qiáng)度之間的關(guān)系Fig.6 (a)Relation between excitation magnetic field intensity and input electric current; (b)Relation between output displacement and input electric current

由于所設(shè)計(jì)的超磁致伸縮驅(qū)動(dòng)器工作在準(zhǔn)靜態(tài)環(huán)境中，在此沒(méi)有考慮渦流損耗，所以電流強(qiáng)度和磁場(chǎng)強(qiáng)度呈線性關(guān)系；輸出位移量和試驗(yàn)值具有相同的變化趨勢(shì)，最大誤差不大于10%。

3 結(jié) 論

本文根據(jù)Jiles-Atherton磁滯模型、二次疇轉(zhuǎn)模型以及應(yīng)力-應(yīng)變模型建立了不考慮渦流損耗的超磁致伸縮驅(qū)動(dòng)器的二維軸對(duì)稱非線性有限元模型。通過(guò)分別建立機(jī)械應(yīng)力場(chǎng)，電磁場(chǎng)的積分弱解方程，并采用Jiles-Atherton模型來(lái)描述兩場(chǎng)之間的耦合關(guān)系，在COMSOL Multiphysics 3.5 軟件中進(jìn)行求解計(jì)算，得到的驅(qū)動(dòng)器輸出位移和實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了對(duì)比，發(fā)現(xiàn)誤差較小，說(shuō)明了該方法的正確性和實(shí)用性。由于驅(qū)動(dòng)電流的頻率變化對(duì)驅(qū)動(dòng)器的輸出特性有巨大的影響，下一步的工作是研究在計(jì)及渦流損耗的情況下，驅(qū)動(dòng)電流的頻率變化對(duì)輸出特性的影響。

[1]Engdahl G,Svensson L. Simulation of the magnetostrictive performance of Terfenol-D in mechanical devices[J].Appl Phys,1988,63(8)：3924-3926.

[2]J L Perez-Aparicio,H Sosa. A continuum three-dimensional,fully coupled,dynamic,non-linear finite element formulation for magnetostrictive materials[J].Smart Material and Structure,2004,13(3):493-502.

[3]Kannan Kidambi Srinivasan.Galerkin Finite Element Scheme for Magnetostrictive Structures and Composites[M].Department of Mechanical Engineering, University of Maryland,1997.

[4]Benatar J G,Flatau A B.FEM implementation of a magnetostrictive transducer[J].Smart Structures and Materials,Proceedings of SPIE,2005,5764:482-493.

[5]趙章榮，鄔義杰，顧建新，等.超磁致伸縮執(zhí)行器的三維非線性動(dòng)態(tài)有限元模型[J].浙江大學(xué)學(xué)報(bào)(工學(xué)版)，2008，42(2)：204-208.

Zhao Zhangrong, Wu Yijie,Gu Jianxin,et al.Three-dimensional nonlinear dynamic finite element model for giant magnetostrictive actuators[J].Journal of Zhejiang University(Engineering Science),2008，42(2)：204-208.

[6]趙凱華，陳熙謀.電磁學(xué)下冊(cè)[M].高等教育出版社，1985，6：793-795.

[7]MAGSOFT Corporation.ATILA Finite Element Analysis for Piezoelectric and Magnetostrictive Structures.User’s Manual,2003,34.

[8]Jiles D C,Atherton D L.Theory of ferromagnetic hysteresis[J].Journal of Magnetism and Magnetic Materials,1986,61(1-2):48-60.

[9]Jiles D C,Atherton D L.Theory of the magnetization process in ferromagnets and its application to the magnetomechanical effect[J].Phys D:Appl Phys,1984,17(6):1265-1281.