在众多导波模态中,S0模态导波作为典型对称型模态[1-3],具有衰减低、传播稳定等优点,对板内部损伤亦有很高的检测敏感性,因此被视为复合材料板早期损伤识别的理想载体[4-6]。然而,由于复合材料各向异性特性显著,其导波传播规律更加复杂,这为精确建模与信号分析带来挑战。
传统有限元法(finite element method,FEM)[7-8]尽管适应性强,但在处理高频导波或大规模模型时,需要满足每个导波波长范围里至少要有10~12个节点的网格划分要求[9],导致系统自由度数以百万计,刚度矩阵组装与存储开销巨大,从而造成计算效率低下与内存资源消耗严重。相比之下,谱元法(spectral element method,SEM)[10-15]凭借高阶插值函数与稀疏单元的特性,可将单波长所需节点数减少一个数量级,显著提升计算效率。然而,现有研究多集中于串行计算,针对复合材料板导波传播的GPU并行化研究仍较为薄弱[16]。如何在保持计算精度的同时,突破计算效率瓶颈,依然是亟须解决的关键问题。
基于上述背景,本文提出了一种并行时域谱元法(parallel time-domain SEM)。该方法基于高阶谱元离散构建复合材料板导波传播模型,并结合GPU并行计算架构,实现大规模模型的高效求解。同时,通过信号处理方法分离时域响应,提取S0模态特征并计算波速分布曲线。在数值模拟基础上,进一步搭建压电激励和激光接收实验系统,对T300复合材料板进行实验验证。本文方法为复合材料板结构无损检测及健康监测提供了有效的技术路径。
1 基本理论
1.1 谱单元建模
在主结构和压电传感器结构中,弹性导波传播的平衡方程可以表示为
$\nabla \boldsymbol{\sigma}=\rho \frac{\partial^{2} \boldsymbol{u}}{\partial t^{2}}-f$。
式中:u表示位移向量;t表示该位移的时间变量;ρ为质量密度;f为体载荷;σ为应力矢量。在压电结构耦合场中还需要满足电学平衡方程
$\nabla \boldsymbol{U}^{\mathrm{D}}+\rho_{\mathrm{e}}=0 \text { 。 }$
式中:UD是电位移矩阵;ρe为自由电荷密度。压电本构方程为
$\begin{array}{l}\boldsymbol{\sigma}=\boldsymbol{c}^{\mathrm{E}} \boldsymbol{\varepsilon}-\boldsymbol{d} E, \\\boldsymbol{U}^{\mathrm{D}}=\boldsymbol{d} \boldsymbol{\varepsilon}+\boldsymbol{g}^{\boldsymbol{\varepsilon}} E 。\end{array}$
式中:ε为应变矢量;cE为恒定电场强度下的弹性常数矩阵(上标E表示电场强度固定时测得);d为压电常数矩阵;gε为应力恒定时的介电常数张量(上标ε表示应变固定时测得);E为电场强度。
对大面积复合材料板进行网格划分时,采用36节点谱单元(如图1所示)。谱单元节点定义为 Gauss-Lobatto-Legendre(GLL)方程的根:
$\left(1-\xi^{2}\right) \frac{\mathrm{d}}{\mathrm{~d} \xi} P_{n}(\xi)=0$。
式中:Pn(ξ)为n阶Legendre多项式;ξi(i=0,1,2,…,n)为谱单元节点坐标,相应的积分权重为
$\left(\begin{array}{l} u^{(c)} \\ \varphi_{x}^{(c)} \\ v^{(c)} \\ \varphi_{y}^{(c)} \\ w^{(c)} \end{array}\right)=\left(\begin{array}{cccccccccccc} S_{1}(\xi) S_{1}(\eta) & 0 & 0 & 0 & 0 & \cdots & S_{n}(\xi) S_{n}(\eta) & 0 & 0 & 0 & 0 \\ 0 & S_{1}(\xi) S_{1}(\eta) & 0 & 0 & 0 & \cdots & 0 & S_{n}(\xi) S_{n}(\eta) & 0 & 0 & 0 \\ 0 & 0 & S_{1}(\xi) S_{1}(\eta) & 0 & 0 & \cdots & 0 & 0 & S_{n}(\xi) S_{n}(\eta) & 0 & 0 \\ 0 & 0 & 0 & S_{1}(\xi) S_{1}(\eta) & 0 & \cdots & 0 & 0 & 0 & S_{n}(\xi) S_{n}(\eta) & 0 \\ 0 & 0 & 0 & 0 & S_{1}(\xi) S_{1}(\eta) & \cdots & 0 & 0 & 0 & 0 & S_{n}(\xi) S_{n}(\eta) \end{array}\right)\left(\begin{array}{c} \hat{u}_{1}^{(c)} \\ \hat{\varphi}_{x, 1}^{(c)} \\ \hat{v}_{1}^{(c)} \\ \hat{\varphi}_{y, 1}^{(c)} \\ \hat{w}_{1}^{(c)} \\ \vdots \\ \hat{u}_{n}^{(c)} \\ \hat{\varphi}_{x, n}^{(c)} \\ \hat{v}_{n}^{(c)} \\ \hat{\varphi}_{y, n}^{(c)} \\ \hat{w}_{n}^{(c)} \end{array}\right)=\boldsymbol{S}^{(c)} \boldsymbol{U}^{(c)} 。$
式中:S(c)为单元形状函数;u(c)、v(c)和w(c)表示各方向位移; ${\varphi }_{x}^{\left(c\right)}$ 和 ${\varphi }_{y}^{\left(c\right)}$ 为复合板中间平面法线分别对x轴和y轴的旋转;U(c)为单元内节点位移向量;Sn(ξ)Sn(η)是基于GLL节点构造的谱单元形状函数;Sn(ξ)是在参考坐标ξ∈[-1,1]上定义的一维插值函数,其中n表示第n个节点在ξ方向的插值函数。
$\begin{aligned}\boldsymbol{m}^{(c)} & =\int_{-1}^{1} \int_{-1}^{1}\left[\boldsymbol{S}^{(c)}(\xi, \eta)\right]^{\mathrm{T}} \boldsymbol{H}^{(c)}(\xi, \eta) \boldsymbol{S}^{(c)}(\xi, \eta) \operatorname{det} \boldsymbol{J}^{(c)} \mathrm{d} \xi \mathrm{~d} \eta=\sum_{i=1}^{n} w_{i} \sum_{j=1}^{n} w_{j}\left[\boldsymbol{S}^{(c)}\left(\xi_{i}, \eta_{j}\right)\right]^{\mathrm{T}} \boldsymbol{H}^{(c)}\left(\xi_{i}, \eta_{j}\right) \boldsymbol{S}^{(c)}\left(\xi_{i}, \eta_{j}\right) \operatorname{det} \boldsymbol{J}^{(c)} ; \\\boldsymbol{k}^{(c)} & =\int_{-1}^{1} \int_{-1}^{1}\left[\boldsymbol{B}^{(c)}(\xi, \eta)\right]^{\mathrm{T}} \boldsymbol{D}^{(c)}(\xi, \eta) \boldsymbol{B}^{(c)}(\xi, \eta) \operatorname{det} \boldsymbol{J}^{(c)} \mathrm{d} \xi \mathrm{~d} \eta=\sum_{i=1}^{n} w_{i} \sum_{j=1}^{n} w_{j}\left[\boldsymbol{B}^{(c)}\left(\xi_{i}, \eta_{j}\right)\right]^{\mathrm{T}} \boldsymbol{D}^{(c)}\left(\xi_{i}, \eta_{j}\right) \boldsymbol{B}^{(c)}\left(\xi_{i}, \eta_{j}\right) \operatorname{det} \boldsymbol{J}^{(c)}。\end{aligned}$
式中:J表示全局与局部坐标映射的雅可比矩阵;B(c)是应变-位移矩阵;D(c)为弹性矩阵;H(c)(ξ,η)为相关的系数矩阵。相应的面密度H0、一阶质量矩H1、二阶质量矩H2分别定义为
式中:m为复合板铺层层数;ρr为第r层材料密度;lr为中性面至第r层上表面的距离;lr-1为中性面至第r层下表面的距离。
运动方程为
$\boldsymbol{M} \ddot{U}+\boldsymbol{C} \dot{U}+\boldsymbol{K} \boldsymbol{U}=\boldsymbol{F}_{\circ}$。
式中:M为全局质量矩阵;K为全局刚度矩阵;C是全局阻尼矩阵;U、 $\dot{U}$ 和Ü分别为节点位移、速度和加速度向量;F为随时间变化的激励向量。零初始条件下,采用中心差分法求解上述运动方程,可表示为
$\begin{array}{l}\left(\frac{1}{\Delta t^{2}} \boldsymbol{M}+\frac{1}{2 \Delta t} \boldsymbol{C}\right) \boldsymbol{U}_{t+\Delta t}= \\\quad \boldsymbol{F}-\boldsymbol{K} \boldsymbol{U}+\frac{2}{\Delta t^{2}} \boldsymbol{M} \boldsymbol{U}_{t}-\left(\frac{1}{\Delta t^{2}} \boldsymbol{M}-\frac{1}{2 \Delta t} \boldsymbol{C}\right) \boldsymbol{U}_{t-\Delta t} 。\end{array}$
式中:t为时间;Δt表示时间步长。
1.2 内力矢量并行计算
在显式时间积分计算中,结构动力学方程的求解效率高度依赖于内力项的计算速度。传统方法通常需组装全局刚度矩阵并与位移向量进行乘法运算,以得到结构内力;但这一过程在处理大规模模型时面临严重的内存与计算瓶颈。尤其在导波传播模拟中,高频激励要求更细致的空间离散,导致系统自由度激增、刚度矩阵极其庞大,组装与存储成本高昂,严重限制计算规模与效率。
针对上述问题,本文采用了无矩阵化的并行策略,将内力计算过程由传统的“矩阵-向量乘法”转化为“单元级内力逐点叠加”,避免了全局刚度矩阵的显式构建与存储。该方法基于谱单元的局部性计算特征,充分利用GPU的高并行性,在单元层面独立完成内力计算,并通过索引映射机制高效集成各单元力分量至全局力向量。
实现该方法的关键在于节点编号重排与单元解耦。由于谱单元在共享节点处存在自由度重叠,直接并行计算易导致数据覆盖。因此,本文引入重复节点分离机制,在预处理阶段对全局节点编号进行重新排序,使得每个单元在逻辑上皆有独立的局部节点集合,确保单元间计算互不干扰。具体操作为,首先建立全局节点与局部节点的索引映射关系:
{Li}=Q({Gi})。
式中,Q表示第i个单元的全局节点集Gi到局部节点集Li的映射。在谱单元离散中,共享节点的编号映射遵循“全局-局部”解耦原则:任一被多单元共用的节点在全局坐标系下具有唯一编号,但在各单元局部坐标系中则分别赋予不同编号,体现了“单全局编号”到“多局部编号”的映射方式。根据索引映射,单元级变量(如形函数、雅可比矩阵、节点位移等)被重新以局部索引排列,具体为
$\boldsymbol{S}_{\xi}^{(p)}=\left(\begin{array}{cccc}S_{\xi}^{(e=1)} & 0 & \cdots & 0 \\0 & S_{\xi}^{(e=2)} & \cdots & 0 \\\vdots & \vdots & \ddots & \vdots \\0 & 0 & \cdots & S_{\xi}^{(e=n)}\end{array}\right)$;
$\left(\boldsymbol{J}^{(p)}\right)^{-1}=\left(\left(\boldsymbol{J}^{(e=1)}\right)^{-1}\left(\boldsymbol{J}^{(e=2)}\right)^{-1} \quad \cdots \quad\left(\boldsymbol{J}^{(e=n)}\right)^{-1}\right) ;$
$\boldsymbol{u}_{i}^{(p)}=\left(\begin{array}{llll}\boldsymbol{u}_{i}^{(e=1)} & \boldsymbol{u}_{i}^{(e=2)} & \cdots & \boldsymbol{u}_{i}^{(e=n)}\end{array}\right) 。$
式中:上标(p)表示单元拆分后用于并行计算的矩阵;(e=n)表示第n号单元。在此基础上,内力计算可在GPU上实现高度并行化:每个线程负责一个谱单元的内力积分任务,针对高斯积分点循环计算应力-应变响应,并借助形函数插值获得单元节点力。所有单元的内力计算完成后,通过预设映射索引将单元节点力合并为全局力向量
FG= ${F}_{L1}^{\left(p\right)}$ + ${F}_{L2}^{\left(p\right)}$ +…+ ${F}_{Ln}^{\left(p\right)}$ 。
式中: ${F}_{Li}^{\left(p\right)}$ 为并行计算得到的单元节点力向量;FG为组装后的全局力向量。将内力FG代入公式(10),即可获得时间步t+Δt时的位移Ut+Δt,其表达式为
Ut+Δt= $\frac{1}{\frac{1}{\Delta {t}^{2}}M+\frac{1}{2\Delta t}C}\left[{F}_{t}-{F}_{G}+\left(\frac{2}{\Delta {t}^{2}}M\right){U}_{t}+\left(-\frac{1}{\Delta {t}^{2}}M+\frac{1}{2\Delta t}C\right){U}_{t-\Delta t}\right]$ 。
式中:M为全局质量矩阵;K为全局刚度矩阵;C是全局阻尼矩阵;U为节点位移;Ft为随时间变化的激励向量;t为时间;Δt为时间步长。
2 数值模拟与实验分析
2.1 数值模拟
为了验证并行时域谱元法在CFRP层压板中导波传播模拟的精度,本文选取尺寸为1 000 mm×1 000 mm×2.4 mm的T300复合层压板作为研究对象。该板共16层,每层厚度0.15 mm,纤维铺层顺序为[0\90\45\-45\45\-45\90\0]s。单层材料力学参数见表1。
表1 T300/M914碳纤维板单层材料参数Tab.1 Material properties of a single ply of T300/M914 carbon fiber composite |
| 参数 | 数值 |
|---|---|
| ρ/(kg·m-3) | 1 560 |
| E1/GPa | 139.9 |
| E2/GPa | 10.5 |
| E3/GPa | 10.5 |
| G12/GPa | 5.7 |
| G13/GPa | 5.7 |
| G23/GPa | 3.38 |
| v12 | 0.313 |
| v13 | 0.313 |
| v23 | 0.478 |
当激励中心频率为400 kHz时,板内S0模态导波的波长λ约为3.9 mm。为保证每个波长范围内包含至少12个节点,时域谱元法模拟时将单元尺寸设置为1.5 mm,对应模型共454 657个谱单元、11 373 096个节点和56 865 480个自由度。谱元法模型中考虑了压电耦合效应,通过在压电片区域施加实际压电激励信号,模拟实验中PZT的激励过程。
作为对比,采用ABAQUS有限元软件对相同工况进行分析。为了使总自由度数量基本一致,将单元类型设置为四边形(4节点),单元尺寸为0.3 mm,网格划分后得到11 108 889个线性单元、11 115 556个节点及55 577 780个自由度。有限元模型未考虑压电耦合效应,因为若在ABAQUS模型中引入压电耦合与PZT激励,则需采用隐式求解器,并组装与求解规模庞大的电-机耦合全局方程组,模型规模会进一步增加,自由度剧增,计算成本远超可接受范围。因此,本对比模型未引入压电耦合,而以等效单点力模拟压电片激励,以保证导波场在可计算前提下产生。
尽管两种方法的自由度规模处于同一数量级,且ABAQUS模型采用16核并行计算提升效率,其求解时间仍长达1 023 min。而在考虑压电耦合、物理过程更完整的情况下,时域谱元法的求解时间仅为834 min,明显优于ABAQUS。可见,即便模拟更为复杂的压电耦合激励,时域谱元法仍能凭借高度并行的单元级运算策略及无全局刚度矩阵组装的优势,展现出显著的计算效率提升。
为测量CFRP层压板的波速分布曲线,设计了CFRP层压板超声导波波速测量系统。如图4所示,以中心 PZT 传感器为圆心,在板表面划定半径100 mm的采样圆,沿圆周每隔10°设置一个信号采样点。图5为并行时域谱元法与ABAQUS有限元法在不同角度(0°、60°、90°等)采样点处归一化时域信号的对比。结果显示,两种方法所得S0模态波包在多个角度的到达时间高度一致,波包峰值对应的时间几乎重合。由于ABAQUS采用等效力加载方式,而谱元法为压电耦合激励,S0波包形态略有差别。总体来看,并行时域谱元法能在导波传播模拟的时域传播特性上与有限元一致。图6展现了不同时刻导波在复合板中的传播云图,可清晰观察A0和S0两种模态在板内的传播过程,波阵面近似为圆形。
图5 并行时域谱元法与ABAQUS有限元法在不同传播方向上的归一化时域信号对比注:网络版为彩图。 Fig.5 Comparison of normalized time-domain signals in different propagation directions between the parallel time-domain spectral element method and the ABAQUS finite element method |
由于各测点到激励点的距离已知,在获得导波到达各测点的时间后,可进一步计算导波在板内任意方向的传播速度v(θ)。将0°~360°范围内离散测点获得的v(θ)按角度排序并插值处理,便可得到连续的波速分布曲线。图7展示了S0模态群速度随传播方向(0°~360°)的方向性变化规律,该规律反映了T300/M914复合材料板因材料各向异性导致的速度差异。从图中可以看出,时域谱元法模拟结果与ABAQUS有限元法在群速度变化趋势上高度一致,总体吻合较好。此外,为进一步验证群速度计算的准确性,将实验采集到的信号进行傅里叶变换,结果表明实际激励的中心频率约为360 kHz。以0°方向为例,选取360 kHz的理论频散群速度值(参见图3)进行对比,发现时域谱元法计算的S0模态群速度与理论值间的差异约为3%,进一步验证了时域谱元法在导波模拟中的数值精度和可靠性。
2.2 实验验证
为进一步验证数值模拟的准确性,搭建了压电激励-激光接收的超声导波实验平台。该系统主要由Tektronix-AFG31102任意函数发生器、Tegam-2350双通道功率放大器、Tempo-2D激光接收仪、运动控制平台以及激光超声成像软件(Laser Ultrasound System V5.2.3)组成。Tempo-2D激光接收仪可同时测量任意表面的离面位移和面内位移。实验装置与试件如图8所示。
实验材料为T300/M914碳纤维复合板。压电传感器采用直径5 mm、厚度0.4 mm的PZT压电陶瓷片,固定于试件中心作为激励源。图9为实验系统的局部放大图,包括中心贴附的PZT激励片、Tempo-2D激光探头布局以及激光光斑在板表面的位置。为确保信号接收质量,激光探头垂直对准采样点,并通过精密运动平台实现逐点扫描定位。
实验中采用5周期正弦调制信号作为激励信号,激励频率为400 kHz。同步与触发方面,函数发生器的触发输出同时送至Tempo-2D的触发输入,作为时间基准,确保激励时刻与激光接收时序的一致。Tempo-2D的采样率设为250 MHz,采样点数为50 000,采集窗口覆盖了从触发起始后的完整传输过程。为降低随机噪声并提升信噪比,每个采样点信号均通过64次重复累计后取平均。激励-接收路径的有效接收时间窗和带通滤波器带宽在数据处理中统一设置(中心400 kHz,带宽±100 kHz),以便有效提取S0波包并抑制其他模态及高频干扰。采样点沿以激励点为圆心、半径100 mm的圆周均匀分布,方位间隔为10°,该布点方式与数值模型中的采样点相对应(见图4)。在每一采样点处,通过运动控制平台精确定位激光探头至目标位置,以接收对应的导波信号。
基于各采样点的导波信号,计算S0模态波包的到达时间。以中心PZT压电片激励时刻作为时间零点,对不同采样点的S0波包峰值时刻进行对比,结合传播距离和时间差,计算S0模态波速,并进一步拟合得到全方位的S0模态波速分布曲线,最后将其与数值模拟结果进行对比,如图10所示。
以中心压电片激励触发时刻为统一时间零点,对各测点采集的导波信号进行带通滤波和Hilbert包络分析,获得S0模态波包峰值对应的到达时间tS0。结合激光测距标定得到的传播距离r,可以按下式计算各角度下的波速:
vS0(θ)= $\frac{r}{{t}_{{S}_{0}}}$ 。
对0°~360°范围内以10°为间隔的36组数据进行插值,可获得连续、光滑的S0模态波速分布曲线。图10对比了实验结果与时域谱元法数值预测结果,可以看到2条波速曲线几乎重合,实验测试值与时域谱元法模拟值的误差在1%~5%,平均相对误差约为3.5%。其中纤维主方向(0°和90°)误差最小,低于2%。这表明时域谱元法能够准确模拟高频导波的传播,对S0模态传播速度的计算有较高精度。
2.3 并行计算加速效率验证
为评估所提出的并行时域谱元法在导波传播模拟中的计算性能,本节对GPU并行计算与传统CPU串行计算在相同模型规模下的运行效率进行了对比分析。实验测试平台的配置为:Intel Core i9-13900K处理器、NVIDIA RTX4090显卡和128 GB DDR5内存。通过记录不同自由度规模下完整时间步计算所需运行时间,对两种计算方式的加速比进行评估,并分析其随模型规模变化的趋势。
图11给出了GPU与CPU在相同模型规模下的计算耗时对比。结果显示,随着模型单元数量和节点数的增加,2种计算方式的耗时均近似线性递增,但CPU计算时间显著高于GPU,且两者差距随模型规模扩大。当模型节点数超过100万时,GPU的计算时间仅为CPU的约1/25,实现了25倍以上的加速。
3 结论
1)本文建立了基于并行时域谱元法的各向异性复合材料板导波传播数值模拟模型。通过引入逐元并行计算策略,实现单元级任务的均衡分配与高效协同,提升了计算性能。模型采用“无矩阵化”单元级内力并行组装策略,避免了全局刚度矩阵存储,有效降低了内存占用,为大规模导波传播数值模拟提供了稳定的基础保障。
2)数值结果表明,该方法对S0模态导波的模拟精度与传统有限元法高度吻合,在400 kHz激励条件下,仅需1.5 mm单元尺寸即可满足单波长≥12节点的精度要求。同时,GPU并行计算有效克服了传统谱元法内存占用大、串行计算效率低的限制,验证了并行时域谱元法在大规模高频导波模拟中的高效性与可行性。
3)与实验结果对比,时域谱元法并行计算的结果与实测数据高度吻合,S0模态波速计算精度优异,能够准确捕获S0模态波包的时域特征,从而证明了该方法在表征各向异性复合材料中S0模态传播规律方面的准确性,为导波检测提供了可靠的数值支撑。