莊曉瑩，李 彬

（同濟大學土木工程學院，上海 200092）

力電耦合效應廣泛存在與各種材料之中，包括人工和天然材料。自然界的很多生物功能都是力電耦合的，如聽覺感知、與動作有關(guān)的神經(jīng)元突起等。這些現(xiàn)象激發(fā)了對此相關(guān)機理的研究，并開發(fā)可以模仿它們的新材料。在所有力電耦合效應中，研究最多的就是壓電效應，它是應變和極化之間的力電耦合效應。

式中：pi是電極化；eijk為三階壓電系數(shù)；εjk為應變。壓電效應是由晶體內(nèi)原子的相對位移引起的，只存在于非中心對稱晶體中。

不同于壓電效應的應變與極化的耦合，撓曲電效應是應變梯度與極化或者應變與極化梯度的耦合。

式中：fijkl為四階撓曲電系數(shù)；ηjkl為高階應變。撓曲電效應產(chǎn)生的原因是晶體材料局部對稱性的變化。傳統(tǒng)宏觀尺度上應變梯度較小，撓曲電效應可以忽略不計，但在微納尺度上結(jié)構(gòu)可以產(chǎn)生很大的應變梯度，使得撓曲電效應在納米尺度上和壓電效應相當。撓曲電效應存在于所有晶體類型中，并在超過居里溫度時仍然出現(xiàn)。相比壓電效應，撓曲電效應更為普遍，而且它還能產(chǎn)生壓電不具備的其他力電功能，關(guān)于撓曲電效應理論及相關(guān)應用可見以下文獻綜述［1-7］。

固體材料結(jié)構(gòu)隨著特征尺度的減小表現(xiàn)出明顯的尺度效應，但經(jīng)典連續(xù)介質(zhì)力學本構(gòu)關(guān)系中不包含任何尺度參量，無法描述微構(gòu)件力學性能的尺寸效應，因此需要引入長度標度，如梯度彈性理論。Mindlin在1965年提出了一個梯度彈性理論［8］，該理論中，應變能密度被認為是應變和一階、二階應變梯度的函數(shù)，從而應力σ、應變ε的本構(gòu)關(guān)系可以寫成［9］：

式中：c1、c2、c3是三個具有長度平方量綱的獨立梯度參數(shù)；λ和G為Lamé常數(shù)；I為單位張量；?為微分算子；?2為拉普拉斯算子；tr為跡運算；?為張量積運算。當c2=0，c1=c3=l2時，式（3）就是在撓曲電效應分析中常見的拉普拉斯型梯度彈性本構(gòu)方程：

式中：l即為材料單參數(shù)的內(nèi)稟尺度。由于撓曲電效應中涉及應變梯度，要求插值函數(shù)至少C1連續(xù)，傳統(tǒng)的有限元方法采用拉格朗日插值往往僅具備C0連續(xù)。常見的應用于撓曲電效應分析的方法有混合有限元方法［10-11］、無網(wǎng)格方法［12-13］、等幾何方法［14-15］?；旌嫌邢拊椒ㄖ?，位移和位移導數(shù)是相互獨立的變量，它們的運動學關(guān)系需要額外的自由度施加約束，從而使得格式較為復雜。對于無網(wǎng)格方法和等幾何方法，它們不是插值類型的單元，因此不容易處理復雜的邊界條件，在商用軟件中也不容易擴展。C1型有限單元開始用于平板的受彎分析，相關(guān)學者開發(fā)了一系列的C1型有限單元［16］，其中Argyris三角形單元也用于撓曲電效應的分析［17］，但對于C1型有限單元在撓曲電效應分析的研究還很少。

本文主要研究Bell三角形單元，它是由Argyris單元簡化而來，其節(jié)點自由度包含位移及全部的一階和二階導數(shù)，Bell單元已經(jīng)被用于許多梯度彈性問題的求解［18-20］?；谔荻葟椥岳碚?，給出了Bell單元分析撓曲電效應的一般格式，并通過數(shù)值算例驗證了Bell單元的準確性和收斂性，分析了內(nèi)稟尺度對結(jié)構(gòu)變形的影響。此外，利用撓曲電結(jié)構(gòu)的非對稱性設(shè)計壓電結(jié)構(gòu)，闡述了相關(guān)原理，并說明了撓曲電的尺寸效應。

1 撓曲電效應

不考慮壓電效應，對于各向同性撓曲電材料，電焓密度H［10］可以表示為

式中：εij為應變，定義為其中u為位移；Ei為電場強度，定義為Ei=-?，i，其中?為電勢；f1和f2為撓曲電材料兩個相互獨立的系數(shù)；κ為介電系數(shù)；l為材料內(nèi)稟長度，它在梯度彈性理論（式（4））中引入。

本構(gòu)方程可表示為

式中：σjk為應力；τijk為高階應力；Di為電位移；ηijk為高階應變，定義為ηijk=εjk，i。

平衡方程可推得：

式中：bk為體積力；ρ為體自由電荷。

2 C1型三角形單元

Argyris三角形單元（圖1）和Bell三角形單元（圖2）是兩種重要的C1型三角形單元。對于Argyris單元，單元自由度包括節(jié)點位移w及其一階、二階導數(shù)，以及三邊中點法向?qū)?shù)。位移插值為五次函數(shù)，邊界上的法向?qū)?shù)插值為四次函數(shù)。然而，邊中節(jié)點法向?qū)?shù)自由度會顯著增大單元剛度矩陣的帶寬，這三個自由度可通過角節(jié)點的位移導數(shù)消去，即為Bell單元。對于Bell單元，位移的插值依然是五次函數(shù)，而邊界上的法向?qū)?shù)插值變?yōu)槿魏瘮?shù)。撓曲電效應是四階偏微分方程問題，故Argyris和Bell單元均滿足完備性和協(xié)調(diào)性要求，當劃分單元足夠小，有限元解趨近于精確解。

圖1 Argyris三角形單元Fig.1 Argyris triangle element

圖2 Bell三角形單元Fig.2 Bell triangle element

Bell三角形單元是基于插值函數(shù)的協(xié)調(diào)單元，便于建模和施加復雜的邊界條件，也很容易在商用軟件中擴展。對于撓曲電效應分析，Bell單元每個節(jié)點有18個自由度，分別是位移、電勢以及它們的一階、二階導數(shù)。

撓曲電效應弱形式可表達為

式中：b、t、ρ和ω分別為體積力、面積力、體自由電荷和表面電荷密度。

離散后，單元內(nèi)部的位移場和電勢場可表示為

式（12）中：ve和φe分別為單元節(jié)點位移和電勢；Nu和N?分別為單元位移和電勢形函數(shù)，定義為

式中：xj和yj為坐標值；i、j、k的取值采用指標輪換，即1→2，2→3，3→1。N7～N18的取值也采用同樣的指標輪換。

應變ε、高階應變η及電場E的取值為ε=(ε11ε222ε12)T，η=(η111η1222η112η211η2222η212)T，E=(E1E2)T，可計算為

式中：Bu、Hu和B?分別為位移形函數(shù)的梯度矩陣、海森矩陣和電勢形函數(shù)的梯度矩陣。

將式（12）和（24）代入弱形式方程（11），可得如下單元節(jié)點平衡方程：

式中：v和φ分別為結(jié)構(gòu)整體位移和電勢；結(jié)構(gòu)整體耦合剛度矩陣和節(jié)點載荷矩陣表達為

式中：材料參數(shù)C、Q、f、κ所對應的矩陣為

3 算例

Bell三角形單元通過Abaqus用戶單元子程序UEL實現(xiàn)，采用TECPLOT進行結(jié)果的后處理。本文首先通過厚壁筒內(nèi)外表面受壓下變形的數(shù)值解與理論值的對比，驗證Bell單元的準確性，討論內(nèi)稟尺度對結(jié)果的影響。然后，利用撓曲電結(jié)構(gòu)的非對稱性設(shè)計產(chǎn)生壓電效應，并研究撓曲電的尺寸效應。

3.1 厚壁筒平面應變問題

圖3為一厚壁筒，它的內(nèi)徑r0為10 μm，外徑r1為20 μm，其受到內(nèi)外表面電勢差1.0 V作用。此外，內(nèi)外徑施加指定的位移載荷，內(nèi)徑位移ur0為0.045 μm，外徑位移ur1為0.05 μm，撓曲電材料參數(shù)取值參照文獻［10］，見表1。

圖3 厚壁筒模型Fig.3 Model of cylinder

表1 撓曲電材料參數(shù)Tab.1 Material properties of flexoelectricity

不考慮體自由電荷，將本構(gòu)方程（6）—（8）代入平衡方程（9）、（10），可得：

由于厚壁筒結(jié)構(gòu)的對稱性，位移和電勢只與半徑有關(guān)，方程（38）的解可表示為［10-11］

其中：I1和K1為一階修正的第一類和第二類貝塞爾（Bessel）函數(shù)，方程（41）和（42）有C1～C6共6個未知數(shù)，可以通過如下6個邊界條件求得：

厚壁筒的切向位移和電勢均為零，沿筒壁厚度方向徑向位移的理論解以及不同網(wǎng)格密度的數(shù)值解見圖4。對于比較稀疏的網(wǎng)格（環(huán)向48個單元，徑向5個單元），數(shù)值解也與理論解吻合得很好。呈環(huán)狀的位移云圖也反映了位移僅與半徑有關(guān)。圖5為不同網(wǎng)格密度下，徑向位移數(shù)值解與理論解誤差絕對值的對數(shù)坐標圖（兩端點由于施加位移約束，誤差過小，因此略去）?？傮w來講，隨著網(wǎng)格密度的增加，數(shù)值解與理論解的誤差減少。當網(wǎng)格比較稀疏，誤差在對數(shù)坐標下的曲線呈現(xiàn)比較明顯的振蕩。徑向網(wǎng)格數(shù)量為20和40的誤差曲線大致重合。圖6為筒壁厚度方向電勢的理論解和數(shù)值解對比，稀疏網(wǎng)格會造成數(shù)值解與理論值略有偏差，當徑向網(wǎng)格數(shù)量為20時，數(shù)值解與理論值位移曲線幾乎重合。圖7為不同網(wǎng)格密度下，徑向電勢數(shù)值解與理論解誤差絕對值的對數(shù)坐標圖（端點值略去）。隨著網(wǎng)格的加密，誤差同樣減少，當徑向網(wǎng)格數(shù)量加密到20時，誤差曲線已較為光滑，數(shù)值解收斂。從圖5和圖7中可以看出，此時徑向位移和電勢數(shù)值解與理論解的誤差絕對值分別小于10-5μm和10-2V。位移和電勢的數(shù)值解與理論值對比說明了Bell三角形單元的準確性與解的收斂性。

圖4 位移ur的理論解和不同網(wǎng)格密度的數(shù)值解Fig.4 Theory solution and numerical solutions of different mesh grids of ur

圖5 位移ur的理論解和不同網(wǎng)格密度的數(shù)值解誤差絕對值Fig.5 Absolute value of ur error between theory solution and numerical solutions of different mesh grids

圖6 電勢?r的理論解和不同網(wǎng)格密度的數(shù)值解Fig.6 Theory solution and numerical solutions of different mesh grids of ?r

圖7 電勢?r的理論解和不同網(wǎng)格密度的數(shù)值解誤差絕對值Fig.7 Absolute value of ?r error between theory solution and numerical solutions of different mesh grids

材料的內(nèi)稟尺度l對徑向位移的影響見圖8，網(wǎng)格密度設(shè)為188×20，即環(huán)向188個單元，徑向20個單元。從圖中可以看出，數(shù)值解同樣與理論解完全重合。此外，內(nèi)稟尺度l越小，徑向位移曲線彎曲程度越大，即徑向應變（?ur/?r）的梯度越大，這也可以從圖9中更直觀地看出。材料系數(shù)Q反映了高階應力與高階應變的關(guān)系，其表達式中含有內(nèi)稟尺度的平方項，因此較大的內(nèi)稟尺度會減小結(jié)構(gòu)的應變梯度。在撓曲電材料計算中，內(nèi)稟尺度大都直接設(shè)定，其取值方法需要進一步研究。

圖8 不同內(nèi)稟尺度下的位移ur的解Fig.8 Solutions of ur for different intrinsic lengths

圖9 不同內(nèi)稟尺度下的應變梯度ur''Fig.9 Strain gradient ur''for different intrinsic lengths

3.2 利用撓曲電材料結(jié)構(gòu)的非對稱性產(chǎn)生壓電效應

相比壓電材料，撓曲電效應存在于所有電介質(zhì)中，而且不受居里溫度的限制。撓曲電效應的一個重要應用就是利用撓曲電材料（非壓電材料）設(shè)計結(jié)構(gòu)產(chǎn)生壓電效應［21-22］。在正方形中心開孔，當受到均勻拉力時，結(jié)構(gòu)內(nèi)部由于應力分布不均勻就會產(chǎn)生應變梯度。但如果開孔是上下左右均對稱，如圓形等，結(jié)構(gòu)的平均電極化仍為零。因此，產(chǎn)生非零電極化需要利用開孔的不對稱性，如圖10所示的三角形結(jié)構(gòu)。結(jié)構(gòu)邊長S為1 μm，開孔三角形的內(nèi)接圓半徑r為0.4 μm，左邊位移和電勢均為零，右邊施加1.02×10-4N·μm-1的線荷載F，此問題設(shè)為平面應變問題。

圖10 三角形開孔的正方形結(jié)構(gòu)Fig.10 Square structure with a triangular hole

不同內(nèi)稟尺度下，結(jié)構(gòu)右邊電位移D1（水平方向）如圖11所示。由于結(jié)構(gòu)的對稱性，電位移曲線沿y=0.5 μm軸對稱。當內(nèi)稟尺度較大時，電位移很小，而當內(nèi)稟尺度較小時，電位移則比較顯著，這是由于較大的內(nèi)稟尺度減弱了結(jié)構(gòu)的應變梯度。右側(cè)電位移D2（豎直方向）如圖12所示，此時電位移曲線沿中軸反對稱，電位移D2沿邊線的積分為零。由于結(jié)構(gòu)是上下對稱、左右非對稱的，電位移才會出現(xiàn)上述的對稱和反對稱性，這也說明了開孔需是非對稱的才可以產(chǎn)生壓電效應。與圖11類似，當內(nèi)稟尺度較小時，結(jié)構(gòu)中會產(chǎn)生更大的應變梯度，電位移D2會更為顯著。撓曲電材料一個重要特點是尺寸效應，即撓曲電效應在小的尺度下尤為顯著。對于許多材料，撓曲電系數(shù)很小，但在納米尺度下，撓曲電效應非常強，而不能忽略。當逐漸縮小結(jié)構(gòu)邊長S，由微米尺度到納米尺度（內(nèi)稟尺度和荷載以相同比例縮?。?，結(jié)構(gòu)的電位移D1如圖13所示。隨著結(jié)構(gòu)尺寸的減小，結(jié)構(gòu)電位移隨之增大。結(jié)構(gòu)尺寸每減小10倍，電位移相應增大10倍左右。這一特點有利于微納尺度下?lián)锨姴牧辖Y(jié)構(gòu)設(shè)計。

圖11 不同內(nèi)稟尺度下的電位移D1Fig.11 Electric displacement D1 for different intrinsic lengths

圖12 不同內(nèi)稟尺度下的電位移D2Fig.12 Electric displacement D2 for different intrinsic lengths

圖13 不同結(jié)構(gòu)大小的電位移D1Fig.13 Electric displacement D1 for different structure sizes

4 結(jié)語

本文基于Bell三角形單元進行了撓曲電效應的分析。Bell三角形單元是C1型協(xié)調(diào)單元，便于處理復雜邊界條件，也很容易在商用軟件中擴展。給出了在撓曲電分析中Bell單元的一般格式，并在數(shù)值算例中驗證了它的準確性和收斂性。梯度彈性理論中引入的內(nèi)稟尺度，可以用來預測微小尺度材料和結(jié)構(gòu)中的尺度效應。較大的內(nèi)稟尺度會減小結(jié)構(gòu)的應變梯度，減弱正向撓曲電效應產(chǎn)生的電位移。此外，本文還利用結(jié)構(gòu)的非對稱性產(chǎn)生壓電效應，通過電位移的分布闡述了相應的原理。撓曲電效應會隨著結(jié)構(gòu)尺度的減小顯著增大，這有利于其在微納尺度中的應用。

Bell三角形單元的一個劣勢是其節(jié)點自由度較多，另外它擴展到三維單元還存在一定難度。本文采用了單參數(shù)的內(nèi)稟尺度，其取值大都直接給出，此外還有學者采用多個內(nèi)稟尺度參數(shù)進行撓曲電效應的分析［23］。Bell單元的計算效率、擴展問題，以及內(nèi)稟尺度的選取，都需要進一步研究。

作者貢獻聲明：

莊曉瑩：學術(shù)指導，研究方法，撰寫論文，論文修改。李彬：理論推導，數(shù)值計算，撰寫論文。