莊曉瑩，李 彬

（同濟(jì)大學(xué)土木工程學(xué)院，上海 200092）

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

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

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

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

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

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

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

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

1 撓曲電效應(yīng)

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

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

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

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

平衡方程可推得：

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

2 C1型三角形單元

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

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

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

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

撓曲電效應(yīng)弱形式可表達(dá)為

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

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

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

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

應(yīng)變?chǔ)?、高階應(yīng)變?chǔ)羌半妶鯡的取值為ε=(ε11ε222ε12)T，η=(η111η1222η112η211η2222η212)T，E=(E1E2)T，可計(jì)算為

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

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

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

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

3 算例

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

3.1 厚壁筒平面應(yīng)變問題

圖3為一厚壁筒，它的內(nèi)徑r0為10 μm，外徑r1為20 μm，其受到內(nèi)外表面電勢(shì)差1.0 V作用。此外，內(nèi)外徑施加指定的位移載荷，內(nèi)徑位移ur0為0.045 μm，外徑位移ur1為0.05 μm，撓曲電材料參數(shù)取值參照文獻(xiàn)［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)的對(duì)稱性，位移和電勢(shì)只與半徑有關(guān)，方程（38）的解可表示為［10-11］

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

作者貢獻(xiàn)聲明：

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