黃懷緯 黃海博

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

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

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

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

近年來，納米尺度下?lián)锨娏阂蚱鋺?yīng)變梯度大、電響應(yīng)亦較強(qiáng)等優(yōu)點(diǎn)，在能量采集、傳感研究領(lǐng)域備受關(guān)注［18］。Hosseini 等［19］分析了深、淺彎曲功能梯度納米梁的自由振動特性。Arefi 等［20］研究了壓電雙彎曲納米殼的非局部高階電彈性彎曲。Sladek等［21］建立了鐵木辛柯梁模型用于分析開路狀態(tài)下小曲率彎曲納米梁的撓曲電效應(yīng)。

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

通過這一工作，將撓曲電材料本構(gòu)模型的二次開發(fā)代碼內(nèi)化到Abaqus 程序中，實現(xiàn)了對撓曲電結(jié)構(gòu)響應(yīng)行為高效、準(zhǔn)確的模擬。

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

1 基本理論

1.1 撓曲電本構(gòu)方程

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

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

根據(jù)電焓密度表達(dá)式，本構(gòu)方程可改寫為

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

結(jié)合平面應(yīng)變假設(shè)，將正交各向異性材料的本構(gòu)方程改寫為矩陣形式：

1.2 數(shù)值方程

Hu 等［23］根據(jù)虛功原理推導(dǎo)出梯度理論中有限元的變分公式為

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

機(jī)械位移和電勢可以用形函數(shù)Na表示：

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

在有限元V內(nèi)，全局坐標(biāo)中的梯度可由雅克比矩陣J表示為

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

其中，

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

電場強(qiáng)度被視為獨(dú)立變量。因此，

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

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

因此，結(jié)合式（11）和式（12），推導(dǎo)出：

為位移模式矩陣。

因此，

式中為有限元節(jié)點(diǎn)坐標(biāo)轉(zhuǎn)換矩陣。

將式（14）代入式（13）推導(dǎo)得

式中，為電場強(qiáng)度轉(zhuǎn)換矩陣，Si(ξ)為電場強(qiáng)度梯度轉(zhuǎn)換矩陣。

電場強(qiáng)度矢量梯度可推導(dǎo)為

其中，

應(yīng)變張量可以通過類似的推導(dǎo)得到。它的表達(dá)式可以從單元V內(nèi)部的幾何關(guān)系中獲得：

應(yīng)變張量也設(shè)為獨(dú)立變量，因此，

類似地

應(yīng)變的近似值可以表示為

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

最后，推導(dǎo)出應(yīng)變場

應(yīng)變梯度的近似值可以表示為

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

因為

式中，(ξ)、(ξ)為位移應(yīng)變梯度轉(zhuǎn)換矩陣Ψa的子矩陣。

應(yīng)變梯度場可以表示為

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

UEL 子程序中，以上方程的代碼編譯流程如圖1所示。在開始高斯循環(huán)后，根據(jù)形函數(shù)和在高斯點(diǎn)上的配點(diǎn)循環(huán)推導(dǎo)雅可比矩陣。接下來，重新定義本構(gòu)方程，并根據(jù)生成的應(yīng)力和應(yīng)變矩陣更新應(yīng)力和應(yīng)變值。最后，利用計算過程中的狀態(tài)變量導(dǎo)出剛度矩陣，結(jié)束高斯循環(huán)。

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

2 驗證及網(wǎng)格收斂性分析

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

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

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

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

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

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

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

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

3 數(shù)值案例

3.1 邊界條件對撓曲電矩形梁力電耦合響應(yīng)的影響

本節(jié)對不同約束條件下的矩形梁和弧形梁進(jìn)行模擬，以研究不同約束條件下不同梁的機(jī)電耦合響應(yīng)特性和能量采集效率。所有模型均在梁的上表面施加1 nN的均布載荷。首先研究支承方式對開路電壓的影響，在左側(cè)邊緣為鉸鏈支座的情況下，將活動鉸鏈支架從X=0.05D移動到X=D，如圖6 所示。當(dāng)活動鉸鏈支座非?？拷髠?cè)邊緣時，結(jié)構(gòu)逐漸轉(zhuǎn)變?yōu)閹缀慰勺兿到y(tǒng)，當(dāng)X/D接近0時，存在結(jié)構(gòu)失效的風(fēng)險，故研究從X=0.05D開始。

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

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

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

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

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

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

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

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

3.2 圓心角對撓曲電弧形梁力電耦合響應(yīng)的影響

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

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

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

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

4 結(jié)語

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