層合壓電材料沖擊問題的時(shí)域間斷Galerkin有限元方法求解

郭 攀， 衛(wèi)洪濤， 武文華， 徐廣濤， 趙 軍

(1.鄭州大學(xué) 力學(xué)與工程科學(xué)學(xué)院，鄭州 450001；2.大連理工大學(xué) 運(yùn)載與力學(xué)學(xué)部工程力學(xué)系工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室，大連 116024； 3.鄭州大學(xué) 機(jī)械工程學(xué)院，鄭州 450001)

由于壓電材料優(yōu)良的機(jī)電耦合效應(yīng)，近幾十年來，受到機(jī)械、土木、海洋、安全防護(hù)等工程及科學(xué)領(lǐng)域?qū)W者的普遍關(guān)注[1-3]，壓電材料自身性能得到了學(xué)者們廣泛的研究。這類材料目前廣泛應(yīng)用于航空航天、軍事[4]、建筑[5]、超聲[6]等領(lǐng)域內(nèi)各類換能器、傳感器、敏感器的制備。電子監(jiān)測(cè)設(shè)備及物聯(lián)網(wǎng)領(lǐng)域的快速發(fā)展，材料結(jié)構(gòu)及信號(hào)傳輸穩(wěn)定性變得日趨重要[7-10]，使得高靈敏度、集成化、高精度成為電子元器件高科技相關(guān)產(chǎn)業(yè)快速發(fā)展的驅(qū)動(dòng)力[11-12]。另一方面，當(dāng)荷載作用時(shí)間達(dá)到微納尺度時(shí)，材料結(jié)構(gòu)往往表現(xiàn)出復(fù)雜的物理現(xiàn)象。因此，研究沖擊作用下極限時(shí)間尺度壓電材料動(dòng)力學(xué)及波傳播行為顯得尤為必要。

隨著計(jì)算方法以及計(jì)算機(jī)技術(shù)的發(fā)展，數(shù)值計(jì)算已經(jīng)成為科學(xué)研究的重要工具[13-19]。采用數(shù)值方法對(duì)壓電耦合效應(yīng)進(jìn)行研究也是目前學(xué)術(shù)界的熱點(diǎn)。目前采用較多的數(shù)值方法主要包括：有限元方法、有限體積方法、等幾何方法等。在沖擊問題的分析中，眾多研究表明[20-24]，間斷 Galerkin 有限元法(Discontinuous Galerkin Finite Element Method, DGFEM) 具備自動(dòng)引入數(shù)值耗散和濾去虛假的高階模式和數(shù)值振蕩效應(yīng)的能力，進(jìn)而得到了固體力學(xué)界廣泛的注意和研究。區(qū)別于時(shí)域連續(xù)Galerkin有限元方法，DGFEM有效的消除了數(shù)值模擬中波傳播過程中波的虛假的非物理數(shù)值振蕩。在動(dòng)力學(xué)以及波傳播問題時(shí)域求解過程中，DGFEM的重要特征是在空間域和時(shí)間域同時(shí)采用有限元離散，對(duì)問題的半離散(空間離散)控制方程中節(jié)點(diǎn)基本未知數(shù)向量及其時(shí)間導(dǎo)數(shù)向量在時(shí)間域中獨(dú)立分片插值，并允許它們?cè)陔x散的時(shí)間段之間存在間斷，其間斷值通過變分原理確定。

Huges等最早將DGFEM引入固體領(lǐng)域?qū)Y(jié)構(gòu)動(dòng)力學(xué)以及波傳播問題進(jìn)行求解。Wiberg等通過對(duì)基本未知變量和未知變量的導(dǎo)數(shù)在時(shí)間步內(nèi)做分段線性插值提出了一種結(jié)構(gòu)雙曲型動(dòng)力學(xué)方程的DGFEM，其在離散的時(shí)刻處允許變量值的間斷，并通過單自由度和多自由度問題進(jìn)行了驗(yàn)證計(jì)算，結(jié)果顯示該方法有著較高的計(jì)算精度，但這種方法需要求解較大的方程組，計(jì)算的內(nèi)存需求大、效率較低。Mancuso等針對(duì)DGFEM需要求解較大的方程組而帶來的計(jì)算效率較低的問題，對(duì)方程進(jìn)行了降階處理，通過隱式迭代的方法在保證穩(wěn)定性和可控的數(shù)值耗散的前提下提高了計(jì)算效率。Li等在傳統(tǒng)DGFEM的基礎(chǔ)上做了進(jìn)一步的發(fā)展，僅使位移基本向量的時(shí)間導(dǎo)數(shù)向量存在間斷，節(jié)點(diǎn)位移則在離散的時(shí)刻自動(dòng)保證連續(xù)，方程組得到了解耦和縮減，在保持求解精度的同時(shí)，極大地提高了計(jì)算的效率。研究發(fā)現(xiàn)，這類方法在沖擊作用下間斷波模擬過程中，波前會(huì)出現(xiàn)虛假的數(shù)值振蕩，進(jìn)而降低了問題求解的精度，針對(duì)這類現(xiàn)象進(jìn)行研究，進(jìn)一步發(fā)展了這類方法，較好地消除了這類問題模擬出現(xiàn)的波前、波后數(shù)值振蕩強(qiáng)間斷沖擊作用下的壓電動(dòng)力學(xué)控制方程，本質(zhì)上同結(jié)構(gòu)動(dòng)力學(xué)問題和單向延遲熱傳播問題是一致的，屬于雙曲型波動(dòng)方程。因此，在沖擊作用下的壓電動(dòng)力學(xué)及波傳播問題求解中，采用現(xiàn)有的DGFEM仍舊存在計(jì)算效率低的問題，抑或是存在波前及界面間數(shù)值振蕩，精度低、不能較好捕捉波間斷特征等問題。針對(duì)這類問題，本文擬構(gòu)建出壓電材料動(dòng)力學(xué)及波傳播問題的改進(jìn)的DGFEM，并針對(duì)沖擊作用下的層合壓電材料動(dòng)力學(xué)及波傳播過程進(jìn)行求解。進(jìn)而，為壓電元器件的設(shè)計(jì)及優(yōu)化提供技術(shù)支持和精確的數(shù)值依據(jù)。

1 壓電動(dòng)力學(xué)基本控制方程及時(shí)域間斷Galerkin有限元方法

不計(jì)體力的運(yùn)動(dòng)方程和靜電學(xué)平衡方程

(1)

壓電場(chǎng)的本構(gòu)方程

σij=Cijklεkl-ekijEk

(2)

Di=ekijεkl-pikEk

(3)

幾何方程為

(4)

式中：σ為應(yīng)力；ρ為密度；Di為電位移分量；Cijkl為彈性模量；εkl為應(yīng)變分量；ekij為壓電常數(shù)；Ek為電場(chǎng)強(qiáng)度；ρik為介電常數(shù)。

對(duì)于壓電動(dòng)力學(xué)方程，空間域有限元離散后可得

(5)

壓電動(dòng)力學(xué)方程的時(shí)域可以離散為

0<…

(6)

對(duì)任一時(shí)刻tn，其時(shí)刻處的階躍函數(shù)可以表示為

(7)

(8)

(9)

其中，

(10)

基本未知函數(shù)的導(dǎo)數(shù)vn在時(shí)間步內(nèi)采用線性插值，具體形式為

(11)

省略上標(biāo)，式(11)可簡(jiǎn)化為

v=vnλ1+vn+1λ2

(12)

(13)

式(13)分別對(duì)基本未知量及其時(shí)域一階導(dǎo)數(shù)進(jìn)行變分，可以得到

(14)

式(14)解耦后，可得

(15)

其中，

為了消除DGFEM方法波前虛假數(shù)值振蕩，基于高頻振蕩敏感于比例剛度阻尼、低頻振蕩敏感于比例質(zhì)量阻尼的原理，引入比例剛度阻尼Ck=βK，對(duì)原有DGFEM方法進(jìn)行進(jìn)一步發(fā)展。式(14)可以改寫為

(16)

2 算例及分析

算例一 為了驗(yàn)證方法的可靠性，首先考慮具有理論解的一維問題。結(jié)構(gòu)長度為0.1 m，截面積A=1.0×10-4m2，壓電材料為PZT-4，其中密度ρ=7 500 kg/m3，介電常數(shù)ζ33=59.2×10-10F/m2，彈性常數(shù)C33=62.1 GPa，壓電常數(shù)e33=18.58 C/m2。具體模型如圖1所示，結(jié)構(gòu)右端固定接地，左端為開路，并受到荷載作用，具體形式為

F(t)=-360sin(125 660t)Nt∈[0.25 μs]

(17)

所發(fā)展的DGFEM計(jì)算時(shí)，結(jié)構(gòu)劃分成400個(gè)單元，時(shí)間步長取為0.5 μs。理論計(jì)算和數(shù)值計(jì)算得到的時(shí)程結(jié)果如圖2和圖3所示，從圖中可以看出，所發(fā)展DGFEM的位移和電勢(shì)計(jì)算結(jié)果與理論解基本一致，表明方法在力電耦合動(dòng)力學(xué)問題計(jì)算中具有較好的可靠性。

圖1 算例一模型示意圖Fig.1 Schematic of configuration of example 1

圖2 左端點(diǎn)位移時(shí)程的解析結(jié)果與DGFEM結(jié)果比較圖Fig.2 Comparison of displacement distributions between analytical solution and DGFEM over time

圖3 左端點(diǎn)電勢(shì)時(shí)程的解析結(jié)果與DGFEM結(jié)果比較圖Fig.3 Comparison of electric potential distributions between analytical solution and DGFEM over time

算例二 為了驗(yàn)證方法的優(yōu)越性，考慮一維層合壓電材料問題。結(jié)構(gòu)分為三層，材料依次為水泥、PZT-4壓電材料、水泥，厚度依次為0.03 m，0.07 m，0.1 m，截面積為1.0×10-4m2。其中壓電材料屬性同算例一，水泥材料密度ρ=2 500 m3，彈性常數(shù)C33=34.5 GPa。具體模型如圖4所示，壓電層左端開路、右端接地，結(jié)構(gòu)右端固定，左端并受到荷載作用，具體形式為

(18)

計(jì)算時(shí)，結(jié)構(gòu)劃分成200個(gè)單元，時(shí)間步長取為0.5 μs。圖5和圖6分別為5 μs，15 μs，30 μs時(shí)刻結(jié)構(gòu)的位移、電勢(shì)分布圖。從圖5可以看出，5 μs時(shí)刻N(yùn)ewmark方法在波后會(huì)出現(xiàn)虛假的數(shù)值振蕩，15 μs，30 μs時(shí)，波動(dòng)前沿跨過了0.03 m，0.1 m界面，Newmark方法在界面前后均會(huì)出現(xiàn)虛假的數(shù)值振蕩，這類數(shù)值振蕩進(jìn)一步加劇了原有的波前數(shù)值振蕩。從圖6可以看出，相比發(fā)展的DGFEM，Newmark方法位移計(jì)算結(jié)果的數(shù)值振蕩導(dǎo)致了電勢(shì)計(jì)算結(jié)果的較大誤差。從計(jì)算結(jié)果不難發(fā)現(xiàn)，所發(fā)展的DGFEM則較好地消除了上述非物理數(shù)值振蕩，進(jìn)一步較好地捕捉了間斷波及層合界面處波陣面的間斷現(xiàn)象，得到了具有較高精度和穩(wěn)定性的結(jié)果，驗(yàn)證了該方法的優(yōu)越性。

圖4 算例二模型示意圖Fig.4 Schematic of configuration of example 2

圖5 Newmark方法和DGFEM不同時(shí)刻結(jié)構(gòu)位移分布比較圖Fig.5 Comparison of displacement distributions between Newmark and DGFEM at different time

圖6 Newmark方法和DGFEM不同時(shí)刻結(jié)構(gòu)電勢(shì)分布比較圖Fig.6 Comparison of electric potential distributions between Newmark and DGFEM at different time

算例三 為了表明所發(fā)展方法在實(shí)際問題較好的應(yīng)用前景，本算例針對(duì)二維層合壓電材料結(jié)構(gòu)進(jìn)行分析。結(jié)構(gòu)長度和寬度假定為1 m，材料依次為水泥、PZT-6B壓電材料[25]、水泥，不同材料厚度依次為0.25 m，0.5 m，0.25 m。

壓電材料密度ρ=7 550 kg/m3，彈性常數(shù)C11=16.8 GPa，C13=60.0 GPa，C31=60.0 GPa，C33=16.3 GPa,C44=27.1 GPa，壓電常數(shù)e31=-0.9 C/m2，e33=7.1 C/m2，e15=4.6 C/m2，介電常數(shù)ζ11=36.2×10-10F/m2，ζ33=34.0×10-10F/m2。 水泥材料密度ρ=2 500 m3，彈性常數(shù)C33=34.5 GPa，泊松比為v=0.19。如圖7所示，左端開路、右端接地，結(jié)構(gòu)右端固定，左端中點(diǎn)以下部分并受到荷載作用，具體形式為

(19)

計(jì)算時(shí)，結(jié)構(gòu)劃分成20×20個(gè)單元，時(shí)間步長取為1.0 μs。圖8和圖9分別是40 μs，80 μs時(shí)刻N(yùn)ewmark方法和DGFEM計(jì)算得到的位移和電勢(shì)分布圖。從圖8位移分布圖中可以看出，在40 μs時(shí)刻，位移波動(dòng)的前沿剛好到達(dá)0.75 m界面處，Newmark方法僅在波陣面內(nèi)部會(huì)出現(xiàn)虛假的數(shù)值振蕩。80 μs時(shí)刻，位移波動(dòng)跨過了0.75 m處界面，Newmark方法在界面附近前后均會(huì)出現(xiàn)虛假的數(shù)值振蕩，這類數(shù)值振蕩進(jìn)一步加劇了原有的波后數(shù)值振蕩。從圖9不同時(shí)刻電勢(shì)計(jì)算結(jié)果可以看出，Newmark方法位移計(jì)算的較大誤差，導(dǎo)致了電勢(shì)計(jì)算結(jié)果的較大誤差。類似算例二，相比Newmark方法，所發(fā)展的DGFEM則能較好地消除二維壓電動(dòng)力學(xué)問題中表現(xiàn)出來非物理數(shù)值振蕩現(xiàn)象，保證了位移和電勢(shì)計(jì)算結(jié)果具有較高的精度，進(jìn)一步能較好地捕捉間斷波及層合界面處波陣面的間斷現(xiàn)象。算例進(jìn)一步表明算法在高梯度、強(qiáng)間斷沖擊作用下力電耦合動(dòng)力學(xué)問題分析中，具備較好的計(jì)算精度和穩(wěn)定性等優(yōu)越性能，可以為這類問題的精確求解提供有效的技術(shù)支持。

圖8 Newmark(左側(cè))和DGFEM(右側(cè))不同時(shí)刻位移分布圖Fig.8 Comparison of displacement distributions between Newmark and DGFEM at different time

圖9 Newmark(左側(cè))和DGFEM(右側(cè))不同時(shí)刻電勢(shì)分布圖Fig.9 Comparison of electric potential distributions between Newmark and DGFEM at different time

3 結(jié) 論

本文構(gòu)建了壓電動(dòng)力學(xué)問題的改進(jìn)時(shí)域間斷Galerkin有限元方法，針對(duì)力電耦合問題進(jìn)行了求解。算例計(jì)算表明，方法具有較高的精度和可靠性。相比傳統(tǒng)的時(shí)域連續(xù)Galerkin 有限元方法如Newmark 方法，所發(fā)展的方法在沖擊作用下壓電動(dòng)力學(xué)問題求解中，能較好地捕捉位移波、電波傳播過程中間斷；同時(shí)能較好消除波傳播過程中波前波后數(shù)值振蕩；能較好消除功能梯度層合材料界面處虛假的波前波后數(shù)值振蕩。問題的求解為沖擊作用下壓電材料電子元器件的研究提供了高精度的數(shù)值依據(jù)，為進(jìn)一步廣義熱壓電問題的研究奠定了基礎(chǔ)。