基于时域谱单元的三维压电耦合结构波传播分析

鱼则行， 徐 超

(西北工业大学 航天学院，西安 710072)

结构健康监测(Structural Health Monitoring, SHM)是指通过获取响应信号及系统特性等信息评估工程结构受损情况、预测结构寿命的一项智能技术。其中，基于弹性导波的健康监测方法因其响应速度快、对结构影响小、对微小损伤敏感等特点，在实际工程中受到越来越多的关注[1]。弹性波在结构中传播会发生反射、干涉、衍射、相速度频散等现象，这给信号测量和信息提取带来了困难与挑战。因此，有效模拟弹性波在结构中的传播行为和数值研究损伤对波传播的影响现象等对发展基于弹性导波的结构健康监测方法至关重要。

压电片因为体积质量小、频率带宽、易集成等特点，作为导波激发和敏感元件被广泛应用于结构健康监测领域。研究弹性波在含压电元件的结构中传播时，须考虑压电元件的耦合效应对波传播行为的影响。邓庆田[2]采用有限元方法研究了在纵向载荷下层合压电杆的频域动力学响应。李鹏[3]研究了界面、梯度功能材料等因素对波在压电元件中传播的影响，其对复合材料压电设备的设计分析以及测量有重要意义。杜朝亮等[4]则采用回传矩阵-射线方法对横观各向同性的压电结构进行了频域的波传播分析。Wang[5]结合二维单层单元与三维多层压电单元，建立了双压电晶片悬臂梁的有限元方法，研究了结构在压电效应下的静力学与动力学响应。Ippolito等[6-9]基于经典有限元数值模拟方法研究了压电耦合结构中的波传播文献，但需要注意的是传统有限元方法在求解结构弹性波传播问题时存在计算耗费大、精度低等问题，很难应用于复杂结构的大规模波传播模拟。

谱单元法最早应用于求解计算流体动力学问题[10]。这种方法结合了谱方法的快速收敛性与有限元法的复杂几何适应性，特别适合模拟高频、瞬态的弹性波传播行为。Kudela等[11]将时域谱单元推广至波传播问题，并研究了一维结构中弹性波的传播行为。Komatitsch等[12-13]将这种方法推广至二维与三维结构中波传播行为的模拟。徐超等[14]推导了任意四边形的谱单元用于功能梯度材料结构中的波传播问题。冯朝辉等[15]将谱单元方法与压电效应相结合，研究了弹性波在压电耦合结构中的传播。但该文章仅研究了二维结构，忽略了导波在三维结构中的反射与散射，也并未结合实验对所做的理论工作进行验证。李富才等[16]进一步研究了压电耦合结构波传播分析的谱单元方法，并初步进行了实验验证。

近年来，宏纤维压电复合材料(Macro Fiber Composite, MFC)作动器和传感器因其柔性好、易于黏接于曲面结构上等优点在结构健康监测中的应用日益增多[17]。实施结构健康监测时，通常在结构的离散位置布置若干MFC压电作动和传感器。压电器件激发的弹性波在结构中传播本质上是三维的。同时，即使对薄板类结构，其裂纹、腐蚀等损伤也多发生在沿板厚度方向的局部位置，必须要对结构进行三维建模和分析才能有效地模拟弹性导波在结构中的传播行为及其与损伤的相互作用关系。最后，针对同时考虑压电元件和待测结构的耦合问题，对压电、黏接层和结构组合体采用统一的三维单元来描述具有形式简单、易于建模等优点。

因此，本文针对MFC压电耦合结构，推导了一种三维时域压电谱单元，并将其用于模拟弹性波在压电耦合结构中的传播行为。首先给出了三维压电耦合谱单元的推导过程；然后将谱单元计算结果与有限元解进行对比，以验证单元的可靠性与快速收敛性；最后结合实验，验证了谱单元方法对模拟弹性波在结构中传播行为的有效性，并研究了导波在多压电耦合结构中的传播特性。

1 三维压电耦合谱单元推导

含压电元件的结构，其运动平衡方程可写为

(1)

D=0 inΩp

(2)

对于压电材料，其物理方程表示为

σ=Cε-eTE

(3)

D=eε+gE

(4)

式中： 上标h和p分别表示主结构区域和压电结构区域。σ，f，ρ，q和D分别表示应力、体积力、密度、位移场和电位移矩阵。ε，e，g，E和C分别表示应变、耦合系数、介电常数、电场强度和弹性矩阵。

如图1所示，与传统有限元方法不同，谱单元采用非等距内插节点以减少龙格效应的影响，从而实现高阶插值。单元内各方向的插节点由Gauss-Lobatto-Legendre多项式的零点确定，称为GLL点。例如，在x方向，该多项式为

(5)

ξ1=-1
ξ2=-0.654 653 670 707 978
ξ3=0
ξ4=0.654 653 670 707 978
ξ5=1

各GLL插值点在局部坐标系下插值函数如图2所示。

图1 3维单元对比Fig.1 Comparison of 3D elements between

图2 GLL节点在局部坐标系下插值函数Fig.2 Interpolation function of GLL points atlocal coordinates system

若将所得节点进行拉格朗日多项式插值，则位移场可表示为

(6)

根据小变形假设可知，应变-位移关系可表示为

ε=Buq=[B1B2…Bn]q

(7)

在电场内，电场强度与电势的关系可表示为

E=-Bφφ

(8)

式中： 下标u和φ分别表示位移场和电场物理量，φ表示电势场。

结合上述公式，根据Hamilton原理，可得到三维压电耦合动力学问题的谱单元运动方程为

(9)

(10)

其中，

式中：Je为Jacobi矩阵表示单元从局部坐标系到整体坐标系的映射。ωi,ωj,ωk表示GLL积分法则在坐标方向的权重系数，由下式确定

(11)

可知，ωi>0独立于单元，并有如下正交性质

(12)

同样的性质也存在于另外两个方向。因此，谱单元的质量矩阵为对角线形式，这对降低计算耗费，提高计算速度有重要的意义。

对于典型压电材料，Kuu数量级约为108， 而Kφφ数量级约为10-11。 由式(9)～式(10)可知，若直接求解系统的控制方程，会由于两者参数差异巨大导致与Kφφ有关的项在计算过程中由于太小而被湮没。这会引起计算上不可忽略的误差。为克服这种缺点，本文引入静力凝聚法，将Kφφ有关的项凝聚掉以位移场向量表示。因此，与式(9)～式(10)可写为

(13)

(14)

式中：KI和fA分别表示由压电耦合引起的诱导刚度矩阵和等效节点力。

KI诱导刚度矩阵的计算中，Kφφ矩阵由定义可知为非正定矩阵。因此，KI矩阵的计算则与压电材料的电学边界条件相关。

(1) 压电元件作为传感器处于电学闭路：压电元件上下表面接地，电势为0，元件内不含有自由电子。可知，式(10)在电学闭路条件下可写为

(15)

式中： 各矩阵上标为：0表示压电元件下表面；i表示压电元件中间层；n表示压电元件上表面。则压电结构中未知电势可以表示为

(16)

则可知， 诱导刚度矩阵KI可表示为

(17)

(2) 压电元件作为传感器处于电学开路：压电元件下表面接地，电势为0，可知此时式(10)可写为

(18)

压电结构中未知电势表示为

(19)

此时，诱导刚度矩阵KI可表示为

(20)

(3) 压电元件作为驱动器：压电元件下表面接地，上表面施加电压V，式(10)可写为

(21)

未知电势可表示为

(22)

诱导刚度矩阵可表示为

(23)

其中，由上表面电压引起的等效节点力可以表示为

fA=Kuφφ

(24)

根据线性假设，忽略由于形变带来压电单元上各节点的电势变化。因此，φ矩阵中各节点的电势由上下表面电势场线性插值而来。Ke=Kuu+KI,Fe=f+fA分别称为三维压电单元的广义刚度矩阵与广义力矩阵。

通过组装各单元矩阵，三维压电耦合结构总体平衡方程可表示为

(25)

(26)

其中Δt表示积分时间步长，该算法稳定的条件为

(27)

式中：Tn表示结构系统的基本固有周期。

2 数值验证

2.1 单片压电耦合结构

为了验证所推导的三维压电谱单元能有效地模拟弹性波在压电耦合结构中传播行为，考虑贴压电片的薄板结构，其几何尺寸如图3所示。薄板为铝制品，压电片采用MFC材料，长25 mm，宽10 mm，厚0.3 mm。各材料属性如表1所示。

图3 单片压电耦合结构Fig.3 Thin palte with one MFC piezoeletric transducers

薄板结构Eνρ70 GPa0.332 700 kg/m3压电材料C11=3.94×1010 PaC22=2.03×1010 PaC33=3.25×1010 PaC12=1.29×1010 PaC13=0.83×1010 PaC23=0.53×1010 PaC44=0.55×1010 PaC55=0.55×1010 PaC66=1.31×1010 Pae31=13.62 cm-2e33=-4.1 cm-2e32=0.55 cm-2e24=-17.03 cm-2e15=-17.03 cm-2g11=141.2g0g11=141.2g0g11=141.2g0ρ=7 000 kg/m3g0=8.854 19×10-12CV-1m-1

针对压电元件作为驱动器与传感器两种情况分别分析：

(1) 工况一：压电元件作为驱动器，将其下表面接地，上表面接入如图4所示的电压，激励波形为频率100 kHz的5周期的汉宁窗调制正弦波。采用谱单元将结构全局划分为612个单元，单元内插节点数为4×4×4，积分时间步长选取0.01 μs。输出监测点R(0.15,0.015,0.002)处位移-时间历程曲线，并与传统有限元方法对比，结果如图5所示。

图4 激励波形Fig.4 The exciting waveform

图5(a)给出了监测点R处的x方向的位移响应，图5(b)和图5(c)分别给出了监测点y方向和z方向的位移响应。由图可知，在各个方向上，谱单元求得位移-时间历程曲线与传统有限元解能够较好吻合，验证了本文所建立单元的有效性，但是需要说明的是使用两种方法所采用的网格数量有很大差别。表2给出了两种方法在相同精度结果下所采用的计算规模。通过对比可知，谱单元方法能以极小的计算代价获得相当高精度的解。

图5 R点处位移-时间历程曲线Fig.5 Time histories of displacement at point R

SEMFEM总单元数 612241 200总自由度数68 529910 515

在模拟高频导波在结构中的传播行为时，有限元方法对网格密度有很高的要求。通过对比可知，在相同精度的结果下，谱单元方法将结构全局划分的网格数远小于传统有限元,两种方法对应的总自由度数也相差较大。因此，谱单元在研究弹性波在结构中传播规律时，能够大大降低计算耗费与必要的存储空间。

(2) 工况二：压电元件用作传感器，将其下表面接地，在薄板末端的中点(0.3,0.025,0.002)施加一垂直冲击载荷，激励波形仍如图4所示的100 kHz调制正弦波。网格划分与积分步长的选择与工况一保持一致。输出压电片上表面平均电势，并与传统有限元方法对比，结果如图6所示。

图6 压电片上表面平均电势Fig.6 The average voltage of the top surface ofpiezoelectric sensor

由图6可知，当压电元件作为传感器时，谱单元解与传统有限元解能够较好地吻合，说明了所建立谱单元的有效性。但两种方法所得结果在A0模式存在微小的相位偏差，且谱单元解波速略快于有限元解。这是由于在求解高频响应的A0波时，有限元方法对厚度方向的网格尺寸提出了更高的要求，即在每个波长内，传统有限元法需要更多的节点数才能使结果趋近于解析解，而谱单元法则可以较少的节点快速收敛于精确解。如在Kin等对这一现象也有所提及，根据Richardson外插公式与有限元理论方法，在每个A0波长内，有限元方法需至少80个节点才能近似于解析解，但谱单元法仅需10个节点即可收敛于精确解。

综上所述，本文所建立的三维时域谱单元能够有效地模拟弹性波在压电耦合结构中的传播行为，且相较于传统有限元方法，谱单元法在计算规模上具有显著的优势，并且精度更高。

3 实验验证

3.1 压电耦合实验结构

在结构健康监测中，通常采用一个压电片作驱动器，激励起弹性波在主结构中传播，另一个压电片作传感器，检测弹性波的传播并生成相关的响应信号。这种“一收一发”的压电耦合结构通常称为pitch-catch结构。为了研究弹性波在pitch-catch结构中的传播，考虑一黏有两压电片的铝制方板，其几何尺寸如图7所示。

图7 pitch-catch结构几何尺寸(m)Fig.7 Geometry of the pitch-catch structure (m)

两个压电元件均为宏压电纤维复合材料传感器(MFC, SMART MATERIAL 公司)，压电片长28 mm，宽14 mm，高0.3 mm。各结构材料属性与2.1节一致。实验装置平台及试件如图8所示。铝板处于自由-自由状态，采样频率为2 MHz，其输入信号如图4所示，为峰值50 V、频率100 kHz的汉宁窗调制正弦波。在实验室环境下对该pitch-catch结构在同样工况进行32次重复测试，测试结果重复性较好，之间差异很小。并对实验结果进行平均，将其与数值模拟结果进行对比。

图8 实验平台示意图Fig.8 Experiment platform

在内插节点数为5×5×5情况下，计算不同网格规模的波传播结果并进行收敛性分析，结果如图9所示。可知，采用谱单元方法将结构全局离散为3 265个单元时，响应趋于收敛，计算结果与实验记录结果的对比如图10所示。

图9 谱单元收敛性分析结果Fig.9 The convergence study of SEM

图10 pitch-catch结构传感器上表面平均电势Fig.10 The average voltage of the top surface of sensor

由图10可知，通过与实验结果的对比，谱单元所得结果与实际情况符合较好。谱单元在模拟高频导波传播时，对于波速的计算尤为准确。由于在实验过程中，采用柔性支撑来模拟试件自由-自由的边界条件以及试件之间的胶层等因素可能造成了两种方法的响应幅值上的差异。相较于有限元方法，本文所推导的谱单元能够更高效、准确地模拟弹性波在压电耦合结构中的传播情况，对于损伤监测、冲击识别等结构健康监测技术具有重要意义。

4 结 论

本文针对压电元件在结构健康监测领域中的广泛应用，推导建立了一种三维实体压电耦合谱单元，分别研究了压电元件用作激励器与传感器两种不同工况下，结构中的波传播行为。并结合有限元解与实验结果，研究谱单元方法在研究波传播问题的应用。主要结论如下：

(1) 通过与有限元解结果的比较，说明了谱单元方法能够有效地模拟压电元件在单独作动荷单独传感两类功能下的力学行为。通过对比两种方法的计算规模，谱单元方法能够大大降低计算规模，提高运算速度。相较于传统有限元方法，谱单元法在模型弹性导波的A0模式时，能够快速收敛到精确解。

(2) 与实验结果的比较，说明了谱单元方法在模拟多压电元件结构中波传播的有效性，验证了其用于结构健康监测的潜力与优势。