层合压电材料冲击问题的时域间断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 方法，所发展的方法在冲击作用下压电动力学问题求解中，能较好地捕捉位移波、电波传播过程中间断；同时能较好消除波传播过程中波前波后数值振荡；能较好消除功能梯度层合材料界面处虚假的波前波后数值振荡。问题的求解为冲击作用下压电材料电子元器件的研究提供了高精度的数值依据，为进一步广义热压电问题的研究奠定了基础。