微納撓曲電結構力電耦合響應的數值模擬方法

黃懷緯 黃海博

（華南理工大學 土木與交通學院，廣東 廣州 510640）

機電耦合指的是機械運動和電磁現象之間的相互作用。其中撓曲電效應是現代高性能電子線路和微機電系統（MEMS）等領域的一個重要功能模塊。它描述的是應變梯度誘導的電極化現象（正撓曲電效應）和電場梯度誘導的機械應變現象（逆撓曲電效應）。梯度在撓曲電效應中起控制作用。材料尺寸減小且應變不變時，梯度與材料尺寸成反比，這意味著微小尺寸的材料具有更大梯度，故而撓曲電效應也更顯著。由于應變梯度與電場梯度的尺度效應，撓曲電材料在微納米尺度下的高靈敏度傳感［1］及精準致動［2］極具應用潛力。另一方面，相比壓電材料，撓曲電材料普遍不需要經歷苛刻的極化過程，也不會因溫度達到居里點而喪失機電耦合效應，這有效增長了材料的使用壽命［3］。例如，將Gohari 等［4］提出的智能圓柱體壓電輻射聲音控制器改為使用撓曲電材料，可以在極端條件下正常工作。具有撓曲電效應的材料非常廣泛，所有介電材料都可以表現出撓曲電性。本研究提供了有效分析撓曲電結構的方法。

固體材料中撓曲電理論框架完善的過程是漫長的。Mindlin［5］針對電介質提出考慮極化梯度效應的變分原理；Maugin［6］針對介電體建立了電場梯度理論；而后，Maranganti 等［7］給出了非壓電材料力電耦合效應的一般理論框架。Sahin等［8］給出了考慮應變梯度與極化梯度效應的介電材料變分原理。Kalpakides等［9］發展了包含應變梯度和電場梯度的計算理論。Shen等［10］提出針對納米電介質的電焓變分原理，為撓曲電材料的有限元計算提供依據。

撓曲電研究領域的迅速發展使數值模擬的需求日益加劇。然而由于應變梯度和極化梯度的存在，固體電介質中的撓曲電效應由四階偏微分方程控制。傳統僅基于位移的有限元法不適用于處理撓曲電固體。與應變梯度彈性一樣，撓曲電本構方程在有限元中建模困難的問題由位移場近似值中C1單元連續性（函數及其一階導數的連續性）造成。這個問題首先由Abdollahi等［11］使用局部最大熵（LME）無網格方法解決。等幾何分析（IGA）對高階連續性直接處理，是另一種可靠的撓曲電計算方法［12］。另外，如Yvonnet 等［13］的工作一樣，Argyris 三角形單元也可用作C1單元連續性要求的補救措施。類似于應變梯度彈性，Mao等［14-16］采用混合有限元法來解決撓曲電問題。Tian等［17］提出用配點法處理混合元。

近年來，納米尺度下撓曲電梁因其應變梯度大、電響應亦較強等優點，在能量采集、傳感研究領域備受關注［18］。Hosseini 等［19］分析了深、淺彎曲功能梯度納米梁的自由振動特性。Arefi 等［20］研究了壓電雙彎曲納米殼的非局部高階電彈性彎曲。Sladek等［21］建立了鐵木辛柯梁模型用于分析開路狀態下小曲率彎曲納米梁的撓曲電效應。

MFEM 與當前大多數有限元代碼的框架兼容?；贔ORTRAN語言編寫UEL子程序，在Abaqus中實現撓曲電材料的本構行為。用戶自定義變量（UVARM）子程序用于傳遞計算過程中的狀態變量，實現計算結果可視化。本文編寫CMFEM通過C0單元來捕獲C1 單元的連續性，相對于傳統拉格朗日方法使用的9 節點高階單元，簡化了計算的同時，又保持了C1 單元的連續性。它避免了引入拉格朗日乘子的額外自由度，并將單元自由度從文獻［14］中提及的87 個減少到12個，從而極大地提升了求解器的運算效率。

通過這一工作，將撓曲電材料本構模型的二次開發代碼內化到Abaqus 程序中，實現了對撓曲電結構響應行為高效、準確的模擬。

對于撓曲電納米結構已有一些討論，然而商業有限元軟件中撓曲電仿真方法仍然有待探索。主流有限元仿真軟件中撓曲電相關的模塊需要通過二次開發寫入撓曲電本構方程。為此，本研究于Abaqus中開發了一個使用CMFEM 的UEL 子程序模擬撓曲電納米結構，推導了包含撓曲電效應的力電耦合本構方程和數值方程。通過與現有成果對比，驗證了UEL子程序并分析了網格收斂性；并以撓曲電納米梁為數值案例，在Abaqus 中數值模擬了不同邊界條件下的矩形納米梁和弧形納米梁的力電耦合響應。

1 基本理論

1.1 撓曲電本構方程

根據Shen 等［10］提出的考慮應變梯度和應力梯度影響的機電耦合理論，線性電介質體中的電焓密度表示為

式中，aij為二階介電常數張量，cijkl為四階彈性常數張量，ekji和fijkl分別為壓電系數和直接撓曲電系數，張量gjklmni為應變梯度彈性的高階彈性系數，dklij和hijkl分別表示逆撓曲電系數和高階電參數，εij和εkl為應變張量，Ei、Ej和Ek為電場矢量，ηjkl、ηmni為應變梯度張量，Ei，j和Ek，l為電場梯度張量。

根據電焓密度表達式，本構方程可改寫為

其中σij、Di、τjkl和Qij分別為應力張量、電位移張量、高階應力張量和高階電位移張量。高階彈性參數gjkmni=l2cjkmnδli，與一般彈性剛度系數cjkmn、內部長度材料參數l和克羅內克函數δli相關［22］。同樣，高階電參數hijkl=q2aikδjl，由介電常數aik、另一個材料長度參數q和克羅內克函數δli組成。撓曲電系數fijkl=f1δjkδil+f2(δijδkl+δikδjl)，僅由兩個獨立分量f1和f2表征。與直接撓曲電系數fijkl類似，逆撓曲電系數dklij=d1δijδkl+(δikδjl+δilδjk) (d2δkiδlj+d3δkjδli)由3個獨立的系數d1、d2和d3組成。

結合平面應變假設，將正交各向異性材料的本構方程改寫為矩陣形式：

1.2 數值方程

Hu 等［23］根據虛功原理推導出梯度理論中有限元的變分公式為

式（5）右端為外力功，其中，tˉi、Rˉi、Sˉ、Zˉ分別表示力的密度、應力密度、電荷和高階電荷，δui、δsi、δ?、δp表示位移的變分、應變的變分、電勢的變分、電場的變分，Γt、ΓR、ΓS、ΓZ為對應邊界。

機械位移和電勢可以用形函數Na表示：

其中，x表示有限元V上的節點，ξ1、ξ3為有限元V上的局部坐標系，a為高斯積分點，{qa}和?a分別是節點位移和電勢。

在有限元V內，全局坐標中的梯度可由雅克比矩陣J表示為

式中，矩陣Y為雅克比矩陣J的逆矩陣。因此，

其中，

根據麥克斯韋方程，電場可以由式（11）給出：

電場強度被視為獨立變量。因此，

其中：(x) 為配點x上i方向的電場強度；{P(ξ)}T=(1ξ1ξ2ξ1ξ2)，為位移模式向量；為廣義坐標。

由于電勢和電場強度的電學約束通過在高斯點處配點得到滿足，故由使在局部坐標系(，)上單元V中選高斯點的配點xc處的兩個近似值相等來確定，即

因此，結合式（11）和式（12），推導出：

為位移模式矩陣。

因此，

式中為有限元節點坐標轉換矩陣。

將式（14）代入式（13）推導得

式中，為電場強度轉換矩陣，Si(ξ)為電場強度梯度轉換矩陣。

電場強度矢量梯度可推導為

其中，

應變張量可以通過類似的推導得到。它的表達式可以從單元V內部的幾何關系中獲得：

應變張量也設為獨立變量，因此，

類似地

應變的近似值可以表示為

其中，{L(ξ)}T={P(ξ)}TA-1，L(ξ)為位移應變轉換矩陣。

最后，推導出應變場

應變梯度的近似值可以表示為

其中，{Sk(ξ)}T={P，k(ξ)}TA-1。

因為

式中，(ξ)、(ξ)為位移應變梯度轉換矩陣Ψa的子矩陣。

應變梯度場可以表示為

結合式（15）、（16）、（23）、（28），變分公式式（5）可以簡化為

UEL 子程序中，以上方程的代碼編譯流程如圖1所示。在開始高斯循環后，根據形函數和在高斯點上的配點循環推導雅可比矩陣。接下來，重新定義本構方程，并根據生成的應力和應變矩陣更新應力和應變值。最后，利用計算過程中的狀態變量導出剛度矩陣，結束高斯循環。

圖1 UEL子程序編譯流程圖Fig.1 Compilation flow chart of UEL

2 驗證及網格收斂性分析

本研究選用撓曲電梁驗證UEL子程序。計算模型如圖2所示，梁的左端邊界固支，右端受一豎向集中力Q作用，且右側邊緣的電學邊界條件為接地。通過在Abaqus 中進行幾何建模、材料定義和網格劃分并編寫好INP文件完成前期工作，隨后調用UEL子程序完成數值仿真。值得注意的是，使用UEL子程序計算后，Abaqus可視化模塊中只顯示積分點的相關信息。另編寫UVARM 子程序，將積分點的結果映射到幾何節點以便分析和獲得結果。因此，本研究得到的結果是后處理的結果。

圖2 納米懸臂梁Fig.2 Nano cantilever beam

懸臂梁的幾何尺寸為長度D=500 nm，厚度W=20 nm，寬度b=10 nm，集中力Q=1 nN。材料選取PZT-5H，其材料參數為E=126 GPa，ν=0.2，l=2 nm，f1=f2=1.0×10-7C/m，a=13×10-9C2/（N·m-2）。圖3為懸臂梁在集中荷載下的機電耦合響應。

圖3 懸臂梁的力電耦合響應Fig.3 Electromechanical coupling response of cantilever beams

觀察圖3（a）中的撓度分布，懸臂梁在端部荷載作用下，撓度沿著長度方向逐漸增加，并呈現立方變化，且梁在自由端取得最大撓度40 nm。將仿真結果與Zhang 等［24］計算所得解析解進行對比，結果如圖4（a）所示，撓曲電懸臂梁的撓度結果高度一致。

圖4 本研究有限元法與解析解［24］、Tian等［17］的有限元解對比Fig.4 Comparison of the FEM in this research with the analytical solution［24］ and the Tian’s results［17］

圖3（b）所示懸臂梁的電勢沿縱向對稱分布，最大電勢為38 mV。此外，在靠近固定端的橫截面上，電位分布沿厚度方向呈線性變化。懸臂梁模型在固定端的彎矩最大，而在自由端的彎矩為零。根據截面彎矩與應變的計算關系ε=My/EIz（y為距離橫截面中心軸的距離，Iz為橫截面的慣性矩），彎矩在固定邊緣處達到最大，應變梯度達到最大，導致極化最強。圖3（c）顯示，固定邊緣處的電場強度最大，達到3.8 MV/m，這與電勢分布結果一致。將電場強度仿真結果與Zhang 等［24］計算所得解析解進行對比，如圖4（b）所示?？梢?，結果的一致性略微差于撓度。由圖4可知，本研究采用的基于Abaqus平臺開發的有限元法，與其他有限元法的撓度與電場強度指標相比，比其他方法都更接近解析解。這驗證了本研究采用的數值方法的有效性。

在Abaqus的有限元模擬中，需要足夠精細的網格以確保計算結果的精度，網格過粗會導致結果不準確。涉及高階彈性和高階電場的計算需要更密集的網格。然而，隨著單元數ρ的增加，用于求解的計算機資源也將增加。因此，在收斂的前提下，控制單元數量來控制計算成本和結果精度之間的平衡非常重要。lnρ可用于表征網格密度，中央處理器（CPU）計算總時間可用于表征計算資源消耗。從圖5可以看出，當單元數為3 000時，計算結果趨于收斂，計算資源消耗適中?？梢哉f，此時計算結果的有效性和計算速度之間達到了平衡。因此以下計算均選擇這個網格密度進行模擬分析。

圖5 不同網格密度下的電場E3和CPU總時間Fig.5 Electric field E3 and total CPU times under different mesh density

3 數值案例

3.1 邊界條件對撓曲電矩形梁力電耦合響應的影響

本節對不同約束條件下的矩形梁和弧形梁進行模擬，以研究不同約束條件下不同梁的機電耦合響應特性和能量采集效率。所有模型均在梁的上表面施加1 nN的均布載荷。首先研究支承方式對開路電壓的影響，在左側邊緣為鉸鏈支座的情況下，將活動鉸鏈支架從X=0.05D移動到X=D，如圖6 所示。當活動鉸鏈支座非?？拷髠冗吘墪r，結構逐漸轉變為幾何可變系統，當X/D接近0時，存在結構失效的風險，故研究從X=0.05D開始。

圖6 均布荷載下的簡支梁Fig.6 Simply supported beam under uniformly distributed loads

圖7示出了支座移動時梁的開路電壓、最大撓度以及x方向應變在y上的梯度絕對值|η113|的關系。在矩形梁中應變梯度|η113|對撓曲電效應起控制作用，故本研究取應變梯度|η113|進行量化分析。當X/D接近0時，梁的撓度迅速增加。此時簡支梁逐漸轉變為幾何可變系統，梁的開路電壓迅速增加。這是由于此時支座附近的應變場急劇變化，導致應變梯度增大，最后提高撓曲電梁的開路電壓。當X/D=0.690時，梁的撓度達到最小，為0.19 nm。當X/D=0.750時，梁取得最小開路電壓4.25 mV，此時梁的應變梯度|η113|也取得最小值26.9 mm-1。此時，梁的彎矩與應變梯度均為一個較小值?？梢苿鱼q鏈支座向右移動，梁的開路電壓、撓度以及應變梯度|η113|增加。當X/D=1.000時，納米梁為典型簡支梁。此時，納米梁的開路電壓、撓度和應變梯度|η113|分別為9.63 mV、1.58 nm 和642.9 mm-1。

圖7 鉸鏈支座移動時的開路電壓、撓度以及應變梯度|η113|Fig.7 Open-circuit voltage，deflection，and strain gradient|η113| when the hinge support is moving

一般情況下，簡支梁的彎矩及應變梯度比起懸臂梁都更小。類似于壓電懸臂梁，撓曲電懸臂梁的力電耦合效應也更顯著。如圖8所示的懸臂梁的左邊緣是滑動支撐的，而活動鉸鏈支架的位置從X=0變為X=D。

圖8 均布荷載下的懸臂梁Fig.8 Cantilever beam under uniform load

圖9 示出了鉸鏈支座從0 移動到D時，梁的開路電壓、最大撓度和x方向應變在y上的梯度的絕對值|η113|。當活動鉸鏈支座設置在X/D=0.550時，梁撓度最小，為0.66 nm。當活動鉸鏈支座設置在X/D=0.575時，梁的開路電壓最小，為7.86 mV，此時梁的應變梯度|η113|為60.5 mm-1。當活動鉸鏈支座設置在X/D=0 和X/D=1.000時，整個梁的最大撓度分別為15.23 nm和25.30 nm，開路電壓分別為40.37 mV和38.72 mV，應變梯度|η113|為272.2 mm-1和243.1 mm-1。X/D=0.100 以左和X/D=0.900 以右的曲線趨勢不同。這是因為，納米梁即將完全轉變為懸臂狀態。納米梁的彎曲狀態轉變為一側完全受拉，另一側受壓，此時豎向的應變的梯度最大。雖然當X/D=0 和X/D=1.000時，梁的最大彎矩相同，但梁的彎矩圖與梁形成的圖形中，X/D=1.000 時有較大的面積。根據結構力學知識可知，此時梁的撓度更大。因此，當X/D超過0.900時，撓度增加得更快。當X/D小于0.100時，梁的撓度增加較慢。當X/D=0和X/D=1.000時，納米梁的電勢云圖分布相似。X/D=0和X/D=1.000時，梁的最大彎矩相同，但梁的受彎狀態并不完全相同，故應變梯度不同，這導致了最后開路的不同。與圖7的簡支梁情況一樣，此時在應變梯度稍大的X/D=0 一側開路電壓更大。

圖9 鉸鏈支座從X=0 移動到X=D 時的開路電壓、撓度和應變梯度|η113|Fig.9 Open-circuit voltage，deflection，and strain gradient|η113| when the hinge support moves from X=0 to X=D

由此可知，矩形梁如果需要俘獲較大的開路電壓，邊界條件可以是懸臂，也可以是一端滑動另一端簡支。然而一端滑動另一條端簡支條件下支座邊緣最大位移比普通懸臂梁最大位移大66%。這是由于撓曲電效應的本質是應變梯度誘導的極化效應。當彎矩較大，沿垂直方向的應變具有較大的梯度，因此可以獲得更強的極化。懸臂邊界條件與圖6中左側邊緣簡支的邊界相比，在確保結構正常工作的前提下，可以獲得更強的力電耦合響應。圖7 和圖9 的結果與撓曲電效應本構方程式（2）相符，即極化由撓曲電系數和應變梯度大小決定。由此，在使用同一材料設計撓曲電梁時，可通過控制變形梯度來控制輸出電壓的大小。

3.2 圓心角對撓曲電弧形梁力電耦合響應的影響

異形梁如弧形梁因其結構不規則，應變分布不規則，應變梯度更大，撓曲電效應更為顯著。，對矩形梁進行彎曲，可以有效增大開路電壓。Sladek等［21］對小曲率弧形梁的撓曲電響應進行了理論分析，其推導結果表明，改變梁的曲率可有效地增大撓曲電懸臂梁的開路電壓。根據Sladek等［21］的計算，梁的E3從矩形梁到圓形角為0.95°的弧形梁，增加了四分之一。對懸臂梁和一端滑動、另一端簡支的矩形梁進行類似彎曲，并研究大曲率的撓曲電弧形梁的力電耦合響應，結構示意圖如圖10所示。

圖10 豎向均布荷載下的弧形梁Fig.10 Curved beam under vertical uniformly distributed load

當弧形梁的邊界條件為左端固定，研究弧形梁在1 nN 豎向均布荷載作用下，圓心角取值-90°至90°的撓曲電響應。模擬不同圓心角懸臂梁的撓曲電響應，匯總得到圖11（a）。其中0°到90°表示梁向下彎曲，-90°到0°表示梁向上彎曲。當圓形角為0°時，梁處于直梁狀態，梁的撓度最大，但與彎曲情況相比，開路電壓較小。這表明改變梁的彎曲程度，可有效提高梁輸出的開路電壓，減小梁的撓度。相同彎曲程度下，弧形懸臂梁向下彎曲的開路電壓大于弧形懸臂梁向上彎曲的開路電壓。當弧形懸臂向下彎曲且圓心角為52°時，獲得的開路電壓最大，其值為102.47 mV，比相同條件下矩形梁的40.37 mV大153.8%。

圖11（b）示出了弧形梁在左端滑動約束和右端鉸接的邊界條件下，向上和向下彎曲的撓曲電響應。與弧形懸臂梁一樣，圓心角的正負也代表弧形梁的彎曲方向。同樣，當梁為直梁時，梁的撓度最大，但開路電壓遠小于弧形梁。在一定范圍內增加弧形梁的彎曲程度也可以增加開路電壓，減小撓度。然而，與懸臂梁不同的是，無論是向上彎曲還是向下彎曲，弧形梁在左端滑動約束和右端鉸接的邊界條件下產生的最大開路電壓都更大。如圖11（b）所示，當弧形梁向上彎曲且圓心角為38°時，取得最大開路電壓214.07 mV。這是弧形懸臂梁邊界條件下最大開路電壓102.47 mV 的兩倍多，是矩形懸臂梁的5倍以上。此外，若將撓曲電系數設置為0 C/m，材料PZT-5H 僅考慮其壓電效應，在該模型中，沒有撓曲電效應的情況下，開路電壓將降低89.3%。

4 結語

基于UEL子程序二次開發將撓曲電本構模型內化到有限元軟件Abaqus中，實現了對撓曲電效應的數值模擬，并將單元的DOFs 減少至12。計算出保證網格收斂和計算效率的最優網格密度。數值模擬了納米矩形梁和納米弧形梁在不同邊界條件下的撓曲電響應。結果表明，在均布荷載作用下，以中心角38°向上彎曲的受一邊滑動和一邊鉸接約束的弧形梁，取得弧形梁的開路電壓最大值。本研究通過數值模擬，直觀地展現了撓曲電梁中應變梯度對梁力電耦合響應的影響。在設計撓曲電結構時，為獲取更大的輸出電壓，應從增大變形梯度的角度入手。本研究提供了一種在通用有限元軟件中分析撓曲電微納結構的有效方法，可用于Abaqus 中大型撓曲電結構建模和計算。在大變形的情況下，上述數值方法也適用，例如撓曲電蜂窩結構的模擬。通過處理狀態變量并通過UVARM 子程序輸出，還可以自由輸出其他撓曲電結構的場量，例如強度因子和J積分等參數，這些參數是材料斷裂計算中的重要參數。通過顯式動力分析用戶自定義單元子程序（VUEL）進一步開發，還可以在Abaqus 中實現撓曲電結構的動力響應分析。