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

郭 攀， 衛洪濤， 武文華， 徐廣濤， 趙 軍

(1.鄭州大學 力學與工程科學學院，鄭州 450001；2.大連理工大學 運載與力學學部工程力學系工業裝備結構分析國家重點實驗室，大連 116024； 3.鄭州大學 機械工程學院，鄭州 450001)

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

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

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

1 壓電動力學基本控制方程及時域間斷Galerkin有限元方法

不計體力的運動方程和靜電學平衡方程

(1)

壓電場的本構方程

σij=Cijklεkl-ekijEk

(2)

Di=ekijεkl-pikEk

(3)

幾何方程為

(4)

式中：σ為應力；ρ為密度；Di為電位移分量；Cijkl為彈性模量；εkl為應變分量；ekij為壓電常數；Ek為電場強度；ρik為介電常數。

對于壓電動力學方程，空間域有限元離散后可得

(5)

壓電動力學方程的時域可以離散為

0<…

(6)

對任一時刻tn，其時刻處的階躍函數可以表示為

(7)

(8)

(9)

其中，

(10)

基本未知函數的導數vn在時間步內采用線性插值，具體形式為

(11)

省略上標，式(11)可簡化為

v=vnλ1+vn+1λ2

(12)

(13)

式(13)分別對基本未知量及其時域一階導數進行變分，可以得到

(14)

式(14)解耦后，可得

(15)

其中，

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

(16)

2 算例及分析

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

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

(17)

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

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

圖2 左端點位移時程的解析結果與DGFEM結果比較圖Fig.2 Comparison of displacement distributions between analytical solution and DGFEM over time

圖3 左端點電勢時程的解析結果與DGFEM結果比較圖Fig.3 Comparison of electric potential distributions between analytical solution and DGFEM over time

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

(18)

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

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

圖5 Newmark方法和DGFEM不同時刻結構位移分布比較圖Fig.5 Comparison of displacement distributions between Newmark and DGFEM at different time

圖6 Newmark方法和DGFEM不同時刻結構電勢分布比較圖Fig.6 Comparison of electric potential distributions between Newmark and DGFEM at different time

算例三 為了表明所發展方法在實際問題較好的應用前景，本算例針對二維層合壓電材料結構進行分析。結構長度和寬度假定為1 m，材料依次為水泥、PZT-6B壓電材料[25]、水泥，不同材料厚度依次為0.25 m，0.5 m，0.25 m。

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

(19)

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

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

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

3 結 論

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