基于谱单元并行计算的结构弹性波传播特性模型构建方法

本发明涉及基于谱单元并行计算的结构弹性波传播特性模型构建方法，属于土木工程结构监测与仿真。

背景技术：

1、近年来，基于弹性波的无损检测技术在土木工程结构健康监测领域得到了广泛研究，适用于结构关键构件的损伤精准识别。但弹性波在实际结构中的传播路径十分复杂，存在多次反射、衍射和频散现象，甚至与既有损伤发生耦合，为了实现损伤精准识别，需要深刻理解弹性波在各种复杂结构中的传播行为及其与损伤的相互作用机理。因此，有效地模拟弹性波在结构中的传播特征、构建弹性波传播特性模型以服务于基于弹性波的结构损伤定位，对发展基于弹性波的复杂结构健康监测方法至关重要。

2、弹性波传播的数值模拟可以采用有限元法、有限差分法和边界元法等方法。有限元法对弹性波的求解过程虽然简单，但由于弹性波频率高、波长短，有限元网格尺寸大小需要满足不大于最小波长的十分之一才能够确保计算精度，这使计算效率大大降低，对尺寸稍大结构的模拟所耗费的计算资源会大到不切实际的地步。有限差分法是一种简单有效的数值解法，但当材料结构复杂时不能保证求解精度与准确性。边界元法更容易处理边界较为复杂的结构，但在计算无限介质中的弹性波发散时才更为有效。

3、在土木工程桥梁结构中，正交异性钢桥面板凭借整体性好、承载能力大、自重轻、施工快、结构美观等优点，已成为国内外大跨度桥梁普遍采用的桥面结构形式，但其极易产生疲劳裂纹，给桥梁结构的运营维护带来极大的安全隐患。正交异性钢桥面板由顶板、u型纵肋和横肋等组成，空间构造复杂，有效地模拟弹性波在正交异性钢桥面板中的传播特征并构建弹性波传播特性模型以服务于结构损伤定位，对理解弹性波的传播路径和波场分布、发展基于声发射弹性波的复杂结构健康监测方法至关重要。

技术实现思路

1、技术问题：正交异性钢桥面板几何构造极为复杂，板件厚度变化且立体交叉，其内部弹性波难以保持单一的传播路径，且在板件之间容易发生反射及衍射，此外，弹性波在板类结构中存在明显的频散现象，不同的频率成分具有不同的传播速度和能量耗散程度，这些原因使得基于声发射弹性波的无损检测技术无法为正交异性钢桥面板提供精准的损伤识别结果。

2、谱单元法采用gauss-lobatto-legendre多项式节点分布，兼具有限元法处理边界的灵活性和谱方法的快速收敛性，是一种快速、高效和高精度的高阶有限元法，它克服了有限元网格尺寸需要小于十分之一最小波长的精度要求，大大降低了对计算资源的需求，能够胜任超高频率的声发射弹性波在大型结构中传播过程的模拟，目前基于谱单元法的弹性波模拟主要集中在平板结构中，针对具有复杂几何构造的三维空间结构的数值模拟研究还有待加强。本发明旨在采用三维时域谱单元法模拟高频声发射弹性波在正交异性钢桥面板中的传播过程，并通过并行计算加快求解速度，分析弹性波传播特性，构建结构的弹性波传播特性模型以服务于基于弹性波的结构损伤定位。其对理解弹性波的传播路径和波场分布、辅助提升声发射技术的工程应用潜力具有重要意义。

3、技术方案：为了解决上述技术问题，本发明提出基于谱单元并行计算的结构弹性波传播特性模型构建方法，该方法包括以下步骤：

4、s1，使用基于legendre正交多项式的三维时域谱单元将土木工程结构或土木工程结构构件离散为有限个单元；

5、s2，利用弹性动力学控制方程建立离散后有限个三维时域谱单元的运动方程，计算单元质量矩阵、单元刚度矩阵和单元阻尼矩阵，对形成的严格对角线形式的单元质量矩阵和单元阻尼矩阵以向量形式存储并参与后续计算，组装有限个三维时域谱单元的单元质量矩阵、单元刚度矩阵和单元阻尼矩阵获得系统的控制方程；

6、s3，利用中心差分法和逐元法，将单元刚度矩阵组装成总体刚度矩阵的过程及大规模总体刚度矩阵的运算简化为有限个三维时域谱单元内的等效节点力的计算及向量组装，使系统的控制方程求解简化为向量运算以实施gpu并行加速计算求解系统的动力学响应；

7、s4，对步骤s1中土木工程结构或土木工程结构构件的不同激励位置的工况分配不同的cpu计算核心，通过cpu并行计算进行不同激励位置工况的弹性波传播模拟，获得弹性波响应及结构波场；

8、s5，计算响应点处弹性波的到达时间，构建步骤s1中土木工程结构或土木工程结构构件的弹性波传播特性时差图模型。

9、进一步的，步骤s1中，使用基于legendre正交多项式的三维时域谱单元将土木工程结构或土木工程结构构件离散为有限个单元的具体步骤如下：

10、s11，根据土木工程结构或土木工程结构构件的几何特征与尺寸，将整体空间结构划分为若干个六面体部分；

11、s12，根据若干个六面体部分的几何尺寸分别建立在同一三维坐标系下的空间模型，采用三维时域谱单元对若干个六面体部分进行离散，依次对若干个六面体部分产生的所有三维时域谱单元和节点进行编号；

12、s13，将若干个六面体部分间几何要素连接处的相同节点修改为同一编号以此组装成土木工程结构或土木工程结构构件的整体空间模型，且要保证组装后所有三维时域谱单元和节点编号连续且不重复；

13、s14，所有三维时域谱单元均采用基于legendre正交多项式的三维时域谱单元，全局坐标系下的任意形状三维时域谱单元都可以通过局部坐标系(ξ,η,γ)下的标准单元进行坐标映射得到，标准域内每个主方向内插节点的坐标由下式确定

14、

15、式中，pn′(ξ)、pn′(η)和pn′(γ)分别为局部坐标系下三个主方向的n阶legendre多项式的一阶倒数，其中，legendre多项式的递推公式p(x)为

16、

17、由上式确定三维时域谱单元内所有节点的坐标，三维时域谱单元内部的位移场则通过节点位移与形函数插值得到。

18、进一步的，步骤s2中的步骤如下：

19、s21，从土木工程结构或土木工程结构构件中取一个微元体进行受力分析，根据达朗贝尔原理可知微元体的平衡方程为：

20、

21、几何方程为：

22、

23、物理方程为：

24、σ＝dε

25、式中，f为体积力向量，ρ为质量密度，u为位移向量，t为时间，ε为应变向量，σ为应力向量，d为弹性矩阵，为的转置，表述如下：

26、

27、在不考虑体积力f的情况下可得到三维问题的波动方程为：

28、

29、s22，三维时域谱单元的单元质量矩阵me和单元刚度矩阵ke可表示为：

30、

31、

32、式中，ωi、ωj和ωk分别表示ξ、η和γ方向的gauss-lobatto-legendre积分权重因子，积分权重因子大于0，单元质量矩阵为严格的对角线形式，(ξi,ηj,γk)表示积分点坐标，n表示积分点数量，n为形函数矩阵，b为几何矩阵，je表示全局坐标(x,y,z)与局部坐标(ξ,η,γ)映射的jacobian矩阵，单元阻尼矩阵ce由质量阻尼比例系数η与单元质量矩阵定义如下：

33、ce＝ηme

34、对形成的严格对角线形式的单元质量矩阵和单元阻尼矩阵以向量形式存储并参与后续计算；

35、s23，组装各单元质量矩阵、单元刚度矩阵和单元阻尼矩阵可以得到全局的质量矩阵、刚度矩阵和阻尼矩阵，根据有限元理论，系统运动方程为：

36、

37、式中，m、c、k、f、q、和分别为结构的总体质量矩阵、总体阻尼矩阵、总体刚度矩阵、等效节点力向量、节点位移向量、节点速度向量和节点加速度向量。

38、步骤s2中，形成的三维谱单元的单元质量矩阵me和单元阻尼矩阵ce为严格的对角线形式，以向量形式存储并参与后续计算可以减少计算资源的耗费，有效提升计算效率。

39、进一步的，步骤s3中，基于组装得到的系统运动方程求解方法如下：

40、采用显示动力学计算方法，t时刻的速度和t时刻的加速度通过t+δt时刻的位移qt+δt及t-δt时刻的位移qt-δt插值得到：

41、

42、

43、式中δt表示积分的时间步长，在零初始条件下，步骤s23中系统运动方程迭代格式为：

44、

45、上式中ft为t时刻的等效节点力向量，qt为t时刻的节点位移向量，通过步骤s22中单元阻尼矩阵和单元质量矩阵组装成总体阻尼矩阵和总体质量矩阵后有：

46、c＝ηm

47、将上式代入系统运动方程迭代格式后表示为：

48、

49、其中

50、

51、

52、上式中总体质量矩阵m为严格对角线形式，在求解过程中矩阵以向量的形式存储和参与计算；

53、上式中总体刚度矩阵k为稀疏的大型矩阵，采用逐元法以实施图形处理器gpu并行加速计算，并行的做法是：在三维时域谱单元上进行单元刚度矩阵ke与单元节点位移向量的乘积运算，并将所有m个三维时域谱单元的此运算平均分配至每个gpu计算核心，使分配到每个计算核心的三维时域谱单元数量相同、单元刚度矩阵数量相同，等待所有gpu计算核心任务结束，将结果输出汇总并组装。逐元法gpu并行加速计算思想表示为下式：

54、

55、上式中为t时刻单元等效节点静力向量。

56、步骤s3中，总体刚度矩阵的总装过程及大规模的矩阵运算被简化为各单元内的等效节点力的计算及向量组装，系统的动力学方程求解被简化为向量运算，在迭代过程中有效降低了内存需求，加之采用并行处理，最大化提高了求解速度，因此本方法能够以极小的计算耗费求解弹性波在结构中的传播问题。

57、进一步的，步骤s4中，对步骤s1中土木工程结构或土木工程结构构件的多个激励位置的荷载工况，每个激励位置的荷载工况都需要进行三维时域谱单元法求解，利用中央处理器cpu并行为每一个激励位置的荷载工况分配一个cpu计算核心进行该工况的三维时域谱单元法求解，通过cpu并行计算同时获得多个工况下弹性波响应及结构波场，等待所有工况计算结束，将结果输出汇总，绘制对应工况下各测点波形图与各时刻波场图，分析弹性波传播特性。

58、步骤s4中，通过并行计算进行结构不同激励位置的弹性波传播模拟，同时求解出对应工况下系统动力学响应，最大化缩减计算总耗时。

59、进一步的，步骤s5中，采用赤池信息准则计算步骤s1中土木工程结构或土木工程结构构件在不同激励位置的荷载工况下，各响应点处弹性波的到达时间，构建结构的弹性波传播特性时差图模型，用于结构损伤的定位识别；赤池信息准则的计算公式如下：

60、aic(t)＝tlog10(var(r(1,t)))+(t-t-1)log10(var(r(t+1,t)))

61、其中，整个时间序列t被分为具有最大似然估计的前半部分(1,t)和后半部分(t+1,t)两个时间窗口，var(r)表示r的方差，r(1,t)为从1到t的时间窗口，r(t+1,t)为从t+1到t的时间窗口，该方程的最小值对应的时间t为弹性波的到达时间。

62、有益效果：与现有技术相比，本发明的技术方案具有以下有益技术效果：

63、(1)离散模型的拟合与计算精度高：本发明采用基于legendre正交多项式的三维时域谱单元对结构进行离散与计算，能够很好地拟合复杂空间形状，高阶内插值节点在单元边界处分布得更加密集，能够提高单元边界处的插值精度，有效抑制runge现象带来的误差，从而实现精确的高阶位移场插值，保证了离散后模型的精确性，兼具有限元法处理边界的灵活性和谱方法的快速收敛性。

64、(2)内存需求低：本发明通过高阶位移场插值，克服了有限元网格尺寸需要小于十分之一最小波长的精度要求，使得结构离散后产生的未知量减少，大大降低了内存需求；本发明形成的质量和阻尼矩阵为严格对角线形式，通过以向量的形式存储和参与计算也降低了内存需求；本发明采用逐元法避免大规模的矩阵运算，亦有效降低了内存需求。

65、(3)求解速度快：本发明将系统的动力学方程求解简化为向量运算，加之采用并行处理，最大化提高了求解速度；对不同激励位置工况，通过并行计算进行结构不同激励位置的弹性波传播模拟，同时求解出对应工况下系统动力学响应，最大化缩减计算总耗时。

66、(4)适用范围广：本发明适用于土木工程结构或土木工程结构构件中弹性波传播仿真，尤其能胜任超高频率的声发射弹性波在大型结构中传播过程的模拟。

67、(5)应用潜力大：本发明能够以极小的计算耗费模拟高频弹性波在土木工程结构或土木工程结构构件中的传播过程以构建结构的弹性波传播特性模型，对理解弹性波的传播路径和波场分布、服务于基于弹性波的结构损伤定位、发展基于声发射弹性波的结构健康监测方法至关重要。