磁電彈結(jié)構(gòu)多場耦合分析的穩(wěn)定Node-based 光滑徑向基點插值法

任書慧，孟廣偉，王吉賢，周立明

（吉林大學機械與航空航天工程學院，吉林長春 130025）

磁電彈材料是一種由壓電相（BaTiO3）與壓磁相（CoFe2O4）復合而成的智能材料，能夠?qū)C械能、電能與磁能相互轉(zhuǎn)化[1-2].磁電彈材料因其具有力電、力磁、磁電效應而被廣泛應用于智能結(jié)構(gòu)中，引起了國內(nèi)外學者的廣泛關(guān)注[3-4].Jiang 等[5]推導了在均布載荷作用下磁電彈懸臂梁響應的解析解，為未來磁電彈結(jié)構(gòu)的設(shè)計與分析奠定了基礎(chǔ).Wang 等[6]求解了三維磁電彈圓柱板的自由振動問題，得出了頻率方程.Ebrahimi 等[7-8]利用Hamilton 原理推導了磁電彈納米板的非局部控制方程，研究了納米板的屈曲行為.

由于傳統(tǒng)有限元法（FEM）在求解中存在“過剛”、體積自鎖等問題，導致結(jié)果不準確.近年來，Liu等[9]提出了廣義梯度光滑技術(shù)，并基于該方法構(gòu)造出了一系列光滑有限元法（S-FEM）[10]和光滑徑向基點插值法（S-RPIM）[11].何智成等[12-13]和陳澤聰?shù)萚14]將光滑有限元應用到聲學模擬中.周立明等[15-16]將光滑有限元擴展到了求解裂紋問題和多場耦合問題中，驗證了光滑方法的準確性.在S-RPIM 中，Node-based光滑徑向基點插值法（NS-RPIM）能夠消除FEM 中“過剛”的問題，為所求解問題提供能量范數(shù)的上界解[9].該方法使用徑向基函數(shù)對場函數(shù)進行近似，其形函數(shù)具有Kronecker Delta 函數(shù)屬性，邊界條件可以如FEM 一樣直接施加.基于伽遼金弱形式與節(jié)點積分技術(shù)，推導出系統(tǒng)方程.基于這些優(yōu)點，NS-RPIM在求解靜力學以及多場耦合問題中得到了廣泛的應用.Li 等[17]采用NS-RPIM 分析了二維、三維固體力學問題，驗證了此算法的準確性與優(yōu)越性.Zhou 等[18]將NS-RPIM 引入多場耦合問題的研究中，結(jié)果表明，NS-RPIM 對于求解磁電彈結(jié)構(gòu)的響應問題是有效且可靠的.

盡管NS-RPIM 在求解許多問題時表現(xiàn)良好，但研究表明[11，19]，由于NS-RPIM 的模型過于柔軟，會令其在求解動態(tài)問題時產(chǎn)生偽非零能量模式，導致算法存在時間不穩(wěn)定性.為增強系統(tǒng)剛度，解決時間不穩(wěn)定性，Wang 等[20]提出了一種穩(wěn)定算法，解決了Node-based 光滑有限元方法中的缺陷并且減少了求解聲學問題中的色散誤差.Feng 等[21]提出了一種穩(wěn)定的節(jié)點積分方法，分析了電磁問題.Yang 等[22]解決了Node-based 光滑有限元方法中的時間不穩(wěn)定性，更準確的求解了金屬成型問題.

本文提出了穩(wěn)定NS-RPIM（SNS-RPIM），基于傳統(tǒng)NS-RPIM 方法，引入與場變量梯度方差相關(guān)的穩(wěn)定項，消除了不確定參數(shù)，推導了求解多場耦合問題的SNS-RPIM 方程，分析了磁電彈結(jié)構(gòu)在靜力作用下的響應，并將所得結(jié)果與有限元法計算結(jié)果進行了對比.

1 基本方程

磁電彈材料平衡方程如下：

式中：σij、Dl、Bl分別為應力分量、電位移分量、磁感應強度分量.

磁電彈材料的幾何方程如下：

式中：Sij為應變分量；Ek為電場強度分量；Hk為磁場強度分量；Φ 與Ψ 為電勢與磁勢.

磁電彈材料的本構(gòu)方程如下：

式中：Cij、eki、qki分別為彈性常數(shù)、壓電系數(shù)與壓磁系數(shù)；εlk、mlk、μlk分別為介電常數(shù)、磁電耦合系數(shù)與磁導率.i，j=1，2，…，6；l=1，2，3；k=1，2，3.

邊界條件為：

式中：Γd與Γs分別為位移邊界與力邊界；Γe與Γt分別為電勢邊界與電位移邊界；Γm與Γi分別為磁勢邊界與磁通量邊界為Γd上給定的位移為Γe上給定的面力為Γe上給定的電勢為Γt上給定的電位移；為Γm上給定的磁勢為Γi上給定的磁通量.

2 穩(wěn)定Node-based 光滑徑向基點插值法

2.1 Cell-based T2L 方案

Cell-based T2L 方案[9]為計算xQ的形函數(shù)值選擇合適的局部支持節(jié)點.該方案選擇xQ周圍的兩層節(jié)點作為局部支持節(jié)點.第一層節(jié)點為xQ所在三角形單元的頂點；第二層節(jié)點為與第一層節(jié)點直接連接的那些節(jié)點，如圖1 所示.

圖1 Cell-based T2L 方案Fig.1 Cell-based T2L scheme

2.2 Node-based 光滑徑向基點插值法

二維問題域Ω 被離散為ne個三角形單元，包含nn個節(jié)點.通過將節(jié)點xi=[xi，zi]T周圍三角形的邊中點與質(zhì)心依次相連，構(gòu)造以xi為中心的光滑域Ωi，如圖2 所示.

問題域內(nèi)任一點xi處的近似位移、近似電勢與近似磁勢可表示為：

式中：Nu（xi）、NΦ（xi）與NΨ（xi）分別為NS-RPIM 的位移形函數(shù)、電勢形函數(shù)與磁勢形函數(shù)；ns為局部支持節(jié)點的數(shù)量；u、Φ 和Ψ 分別表示位移向量、電勢向量和磁勢向量.

圖2 基于節(jié)點xi 的光滑域Fig.2 Smoothing domain of node xi

通過引入梯度光滑技術(shù)，根據(jù)式（16）～（18），節(jié)點xi處的光滑應變、光滑電場強度與光滑磁場強度分別為：

通過引入標準高斯積分，式（22）（23）中的各項可表示為：

式中：nG為高斯點的數(shù)量；nseg為光滑域邊界的數(shù)量；為光滑域第p 段邊界中單位法向量矩陣的分量；為第p 段邊界上第q 個高斯點處的形函數(shù)值；為第q 個高斯點處的權(quán)值；Ai為光滑域的面積.

二維磁電彈的NS-RPIM 靜力學方程可表示為：

式中：F 為外力.各剛度矩陣為：

采用自由度凝聚技術(shù)[18]，式（25）可表示為：

式中：等效力向量Feq、等效剛度矩陣Keq以及電勢Φ和磁勢Ψ 的求解公式見參考文獻[18].

2.3 穩(wěn)定Node-based 光滑徑向基點插值法

為了在提高NS-RPIM 計算精度的同時消除時間不穩(wěn)定性，在該算法中引入與場變量梯度方差相關(guān)的穩(wěn)定項來提高模型的剛度，令其更接近真實情況.

以二維問題為例，如圖3 所示，光滑域Ωi被近似為具有相同面積的圓，將近似域進一步劃分為4個子光滑域.局部坐標系與光滑域的交點1，2，3，4；i=1，2，3，…，nn）作為補充積分點，nn為結(jié)構(gòu)包含節(jié)點數(shù).積分點與節(jié)點xi之間的距離lc相等，大小為近似域的半徑.lc的計算公式為：

圖3 SNS-RPIM 積分區(qū)域Fig.3 Integral domain of SNS-RPIM

假設(shè)場變量的梯度在光滑域Ωi中連續(xù)且一階可導，其在4 個積分點處的泰勒展開式分別為：

式中：m=1，2；n=3，4；l=x，z；上標sc 表示在近似域中的4 個補充積分點處進行計算.將式（31）～式（33）代入光滑伽遼金弱形式，式（25）中各剛度矩陣可以改寫為：

3 數(shù)值算例

3.1 算例1

磁電彈材料（BaTiO3-CoFe2O4）板在邊AD 受100 N/m2的均布載荷作用，為平面應變問題，如圖4 所示，邊長a=2.0 m.表1 給出了磁電彈板的材料參數(shù)，質(zhì)量密度為5 730 kg/m3，邊界條件為：ux=0（邊CD），uz=Ф=Ψ=0（邊BC），每個邊界的表面電荷與表面磁感應強度均為零.

采用不同節(jié)點數(shù)量（121、441 和1 681 個）求解磁電彈板的廣義位移（位移ux、uz，電勢Φ，磁勢Ψ），驗證SNS-RPIM 的正確性以及收斂性.表2 給出了文獻[23]中A 點處位移、電勢與磁勢的解析解，以及SNS-RPIM 在不同節(jié)點數(shù)量下的計算結(jié)果.可見，SNS-RPIM 的結(jié)果與解析解誤差很小，隨著節(jié)點數(shù)量的增加，誤差減小，驗證了SNS-RPIM 求解磁電彈結(jié)構(gòu)多場耦合問題的正確性、有效性.

圖4 磁電彈板Fig.4 Magneto-electro-elastic plate

3.2 算例2

磁電彈材料懸臂梁如圖5 所示，長度L=0.030 m，寬度h=0.002 m，在B 點承受200 N 的靜力，為平面應力問題.懸臂梁在固定端處滿足ux=uz=Φ=Ψ=0.懸臂梁材料屬性見表1，質(zhì)量密度為5 730 kg/m3.

圖5 磁電彈懸臂梁Fig.5 Magneto-electro-elastic cantilever beam

在證明了SNS-RPIM 正確性的基礎(chǔ)上，對磁電彈懸臂梁在靜力作用下的響應進行研究.圖6 給出了邊AB 處的廣義位移，SNS-RPIM 與FEM 采用三角形單元，節(jié)點數(shù)量為305 個.其中，參考解為FEM采用180×12 個四邊形單元的結(jié)果.圖7 給出了靜力作用下懸臂梁的云圖.由結(jié)果可知，在所用節(jié)點數(shù)量相同的情況下，SNS-RPIM 的計算結(jié)果比FEM 的結(jié)果更加接近參考解.算例結(jié)果驗證了SNS-RPIM 的高精度、正確性和有效性.

在節(jié)點數(shù)為305、637 和1 089 個時，對比了SNS-RPIM 與FEM 的能量誤差[24]，如圖8 所示.可知在節(jié)點數(shù)相同的情況下，SNS-RPIM 的能量誤差遠遠低于FEM，并且隨著所用節(jié)點數(shù)的增加，能量誤差逐漸降低.從而進一步驗證了SNS-RPIM 的精確性，高收斂性與有效性.

圖6 懸臂梁AB 邊界處ux、uz，Φ 和Ψ 變化趨勢Fig.6 Variation of ux、uz，Φ and Ψ at edge AB of cantilever beam

圖7 懸臂梁μx、μz、Φ 和Ψ 云圖Fig.7 Contour of μx、μz、Φ and Ψ of cantilever beam

圖8 穩(wěn)定Node-based 光滑徑向基點插值法與有限元法在不同節(jié)點數(shù)下的能量誤差Fig.8 Energy error of SNS-RPIM and FEM with different number of nodes

3.3 算例3

磁電彈材料傳感器為對稱結(jié)構(gòu)，在A 點受100 N 的靜力F，傳感器幾何形狀及尺寸如圖9 所示，為平面應力問題.傳感器邊界條件為：ux=uz=0 （左右固定端處），Ф=Ψ=0（下邊界處）.結(jié)構(gòu)所用材料見表1，質(zhì)量密度為5 730 kg/m3.

圖9 磁電彈傳感器（單位：mm）Fig.9 Magneto-electro-elastic sensor（unit：mm）

圖10 給出了邊BC 處的廣義位移，SNS-RPIM與FEM 采用三角形單元，節(jié)點數(shù)量為1 157 個.其中，參考解為FEM 采用4 000 個四邊形單元的結(jié)果.圖11 所示為磁電彈傳感器位移、電勢和磁勢的云圖.由圖可得：在同樣節(jié)點條件下，SNS-RPIM 的計算結(jié)果比FEM 的計算結(jié)果更加接近參考解，從而驗證了SNS-RPIM 可解決有限元系統(tǒng)剛度過硬的問題，提高結(jié)果的精度，能夠有效求解復雜磁電彈結(jié)構(gòu)在靜力作用下的響應.

圖10 傳感器BC 邊界處ux、uz、Φ 和Ψ 變化趨勢Fig.10 Variation of ux、uz、Φ and Ψ at edge BC of magneto-electro-elastic sensor

圖11 傳感器μx、μz、Φ 和Ψ 云圖Fig.11 Contour of μx、μz、Φ and Ψ of magneto-electro-elastic sensor

4 結(jié)論

本文基于場變量梯度方差構(gòu)造了穩(wěn)定項，并將其引入了傳統(tǒng)Node-based 光滑徑向基點插值法，提出了穩(wěn)定Node-based 光滑徑向基點插值法.隨后求解了磁電彈結(jié)構(gòu)在靜力作用下的響應，得出以下結(jié)論：

1）將SNS-RPIM 的結(jié)果與解析解進行對比，二者吻合良好，說明了本方法的正確性及有效性.

2）計算了SNS-RPIM 與FEM 在不同節(jié)點數(shù)量下的能量誤差，結(jié)果顯示SNS-RPIM 具有良好的收斂性與準確性.

3）SNS-RPIM 利用較少的節(jié)點能夠達到更高的精度，消除了FEM 模型剛度過硬的問題.

4）通過考慮SNS-RPIM 與FEM 對不同模型的求解結(jié)果，表明SNS-RPIM 在求解磁電彈結(jié)構(gòu)多場耦合問題時的可靠性和適用性.