CN105930309A - 一种纤维增强树脂基复合材料层合板声线示踪算法 - Google Patents

Abstract

一种纤维增强树脂基复合材料层合板声线示踪方法，属于复合材料超声检测技术领域。该方法包括以下步骤：以单铺层为单元对计算区域分区，并利用弹性刚度矩阵及其旋转变换，定量描述FRP复合材料弹性特性空间分布；结合Christoffel方程求解，分别获得不同纤维取向铺层对应的准纵波群速度值关于传播方向角的函数关系式；计算区域网格化，利用Dijkstra最短路径搜索算法，搜寻超声波由源点传播至目标点所经过的节点并计算对应声时。该方法能够实现具有多层结构、弹性各向异性以及不同纤维铺放顺序的FRP复合材料中超声波传播路径和声时的快速、精准计算，能够为研究超声波传播行为、优化检测参数、提高超声成像质量和精度提供支持。

Description

一种纤维增强树脂基复合材料层合板声线示踪算法

技术领域

本发明涉及一种纤维增强树脂基(Fiber Reinforced Plastic，FRP)复合材料层合板声线示踪算法，属于复合材料超声检测技术领域。

背景技术

超声检测广泛应用于FRP复合材料内部缺陷的表征与定量检测。然而，FRP复合材料因其层结构以及纤维取向性排布，弹性各向异性及非均质性显著，致使超声波在其中的传播行为复杂，具体表现为声束偏折以及择优取向传播，进而导致超声成像质量差，缺陷定位及定量失准。因此，准确描述超声波的传播路径能够为超声成像相位校正提供支持，并能有效提高检测效果及缺陷检测能力。

声线示踪法通过步进式分析介质局部区域弹性特性及界面对超声波传播的影响，重建出路径信息。可用于预测超声波束的有效覆盖范围以及估算超声波传播时间。现有声线示踪算法多基于Snell定律，计算界面折射所致超声波传播路径的改变。德国联邦材料研究与测试研究所的Sanjeevareddy Kolkoori等人通过迭代使用Snell定律，提出了用于预测各向异性奥氏体钢多层结构透射声场分布的三维声线示踪算法；法国原子能委员会的SJourniac等则是借鉴地震学射线追踪理论，推导出超声波在碳纤维增强树脂基(CarbonFiber Reinforced Plastic，CFRP)复合材料中传播时声线位置矢量及慢度矢量关于传播时间的微分式。然而，上述声线示踪方案均采用迭代求解方式，计算量大，在高精度重建多条声线时显得尤为低效。而源自计算机技术的Dijkstra最短路径搜索算法则提供了一种新的声线示踪方案，能够为超声波传播行为研究检测参数优选以及超声成像相位校正提供快速、精准的分析工具。因此，开展Dijkstra最短路径搜索算法用于FRP复合材料声线示踪的研究，具有十分重要的意义。

发明内容

本发明提出一种纤维增强树脂基复合材料层合板声线示踪算法，可实现超声波传播路径的描绘以及超声波传播时间的精确计算，能够为超声波传播行为研究、检测参数优选以及超声成像相位校正提供快速、精准的分析工具。

本发明采用的技术方案是：一种纤维增强树脂基复合材料声线示踪算法，以单铺层为单元对计算区域分区，并利用弹性刚度矩阵及其旋转变换，实现FRP复合材料弹性特性空间分布的定量描述；将弹性刚度矩阵代入Christoffel方程求解，获得不同纤维取向铺层对应的准纵波群速度值关于传播方向角的函数关系式；计算区域网格化，利用Dijkstra最短路径搜索算法，搜寻超声波由源点传播至目标点所经过的路径并计算对应声时，具体执行如下步骤：

一、将FRP复合材料计算区域以铺层为单元划分为一系列分区，各分区等效为均质弹性各向异性介质，并以弹性刚度矩阵描述分区的弹性特性；

二、以0°纤维取向铺层的弹性刚度矩阵C₀为参考值，其它纤维取向铺层的弹性刚度矩阵值C(θ)由C₀进行Bond变换求得：

$C_0=\left(\begin{array}{cccccc}{c}_{11}& c_{12}& c_{12}& & & \\ c_{12}& c_{22}& c_{22}-2c_{44}& & & \\ c_{12}& c_{22}-2c_{44}& c_{22}& & & \\ & & & c_{44}& & \\ & & & & c_{55}& \\ & & & & & c_{55}\end{array}\right)\left(GPa\right)---\left(1\right)$

$M\left(\mathrm{\θ}\right)=\left(\begin{array}{cccccc}{\sin }^2\mathrm{\θ}& 0& {\cos }^2\mathrm{\θ}& 0& -2\cos \mathrm{\θ}\sin \mathrm{\θ}& 0\\ 0& 1& 0& 0& 0& 0\\ {\cos }^2\mathrm{\θ}& 0& {\sin }^2\mathrm{\θ}& 0& 2\sin \mathrm{\θ}\cos \mathrm{\θ}& 0\\ 0& 0& 0& -\sin \mathrm{\θ}& 0& -\cos \mathrm{\θ}\\ \cos \mathrm{\θ}\sin \mathrm{\θ}& 0& -\sin \mathrm{\θ}\cos \mathrm{\θ}& 0& {\sin }^2\mathrm{\θ}-{\cos }^2\mathrm{\θ}& 0\\ 0& 0& 0& \cos \mathrm{\θ}& 0& -\sin \mathrm{\θ}\end{array}\right)---\left(2\right)$

C(θ)＝M(θ)×C₀×M(θ)^T (3)

其中，θ是纤维铺放角度，M为Bond变换矩阵；

三、将弹性刚度矩阵值C(θ)以及试样密度值ρ代入Christoffel方程，求解得到不同纤维取向铺层准纵波群速度的慢度曲线；

四、借助MATLAB科学计算软件，构建各纤维取向铺层准纵波群速度值关于传播方向角的函数式

其中N表示拟合的阶数，a(n)是拟合多项式的系数序列；

五、建立FRP复合材料层合板试样的笛卡尔坐标系，以超声波源点为坐标原点O，水平向右为x轴正向，深度方向为y轴正向；定义声传播方向角为x轴逆时针方向旋转的角度；则超声波由源点O传播至目标点E的路径及所需时间，可通过执行如下步骤获得：

1.将源点O与目标点E之间的矩形计算区域离散化，其中y轴方向网格化间距等于单铺层厚度δ，x轴方向网格化间距需根据计算精度要求，选择小于单铺层厚度范围内的值；

2.判断源点O(x_O,y_O)与目标点E(x_E,y_E)之间的相对位置关系：

表1

其中，

Δx_OE＝x_E-x_O,Δy_OE＝y_E-y_O (5)

$d_{OE}=\sqrt{{(x_E-x_O)}^2+{(y_E-y_O)}^2}---\left(6\right)$

而当Δy_OE≠0且Δx_OE≠0时，需采用Dijkstra最短路径搜索算法进行声线示踪；此外，当Δy_OE<0时，源点与目标点坐标值互换；

3.在Δy_OE＞0且Δx_OE＞0时，采用Dijkstra最短路径搜索算法进行声线示踪的具体执行步骤如下：

(1)指定节点之间的连接关系：与当前节点i(x_i,y_i)能够直接连接的节点j(x_j,y_j)间的相对位置关系需满足如下条件：

Δx_ij＝x_j-x_i>0且Δy_ij＝y_j-y_i＝δ (8)

(2)计算超声波在邻接节点组合间传播所需时间，在邻接节点i(x_i,y_i)和j(x_j,y_j)，首先根据二节点的坐标值，计算得到其连线对应的传播方向角根据连线所经过铺层的慢度曲线，计算得到准纵波群速度值则超声波由节点i传播至j所需时间

(3)借助MATLAB科学计算软件，采用Dijkstra最短路径搜索算法计算超声波由源点O传播至目标点E的路径及声时，遵循以下搜索原则：

①将源点O指定为起始节点，因此其对应的声时值等于0μs，同时指定超声波由源点O传播至其它所有节点的声时为∞；

②将所有节点状态标记为“未访问”，并分配至“未访问节点集合”S；

③以起始节点O为当前节点，考察超声波经当前节点传播至其邻接节点的路径，相应的声时值为当前节点的声时加上超声波在邻接节点组合间传播的声时；若新路径的声时值小于原有值，则以此新路径代替对应邻接节点的原有路径；反之，则保持路径和声时信息不变；完成遍历后，将当前节点的状态改为已访问，并从集合S中移除；

④比较步骤③所得最新路径下的声时值，将集合S中具有最小声时值的节点作为新的起始节点；

⑤重复步骤①～④，直至集合S为空集，搜索结束；输出超声波由源点O传播至目标点E所经过路径的声线图像及声时值。

本发明的有益效果是：利用该算法可对超声波在FRP复合材料中的传播路径和声时进行快速、精确地计算和描述，能够为超声传播行为研究、检测参数优选、提高超声成像质量和精度提供支持。

附图说明

下面结合附图，以碳纤维增强树脂基(Carbon Fiber Reinforced Plastic，CFRP)复合材料层合板为实施例对本发明作进一步说明。

图1是CFRP复合材料单铺层的准纵波群速度慢度曲线，分别对应(a)0°、(b)45°、(c)90°和(d)-45°纤维取向。

图2是计算区域网格化示意图。

图3是节点间连接关系示意图。

图4是Dijkstra最短路径搜索算法用于声线示踪的原理示意图。

具体实施方式

下面结合附图，以CFRP复合材料层合板为实施例对本发明作进一步说明。其铺放顺序为[0/45/90/-45]，共4层，单铺层厚度δ＝0.15mm，密度为ρ＝1.56g/cm³。

一、将CFRP复合材料层合板试样以铺层为单元，划分为4个分区，各分区等效为均质弹性各向异性介质，并以弹性刚度矩阵描述其弹性特性。其中，0°纤维铺层对应的弹性刚度矩阵值C₀如下

$C_0=\left(\begin{array}{cccccc}126.83& 6.39& 6.39& 0& 0& 0\\ 6.39& 12.88& 6.94& 0& 0& 0\\ 6.39& 6.94& 12.88& 0& 0& 0\\ 0& 0& 0& 2.97& 0& 0\\ 0& 0& 0& 0& 4.85& 0\\ 0& 0& 0& 0& 0& 4.85\end{array}\right)\left(GPa\right)---\left(1\right)$

二、以0°纤维铺层的弹性刚度矩阵C₀为参考值，其它纤维取向铺层的弹性刚度矩阵值C(θ)可由C₀进行Bond变换求得

$M\left(\mathrm{\&theta;}\right)=\left(\begin{array}{cccccc}{\sin }^2\mathrm{\&theta;}& 0& {\cos }^2\mathrm{\&theta;}& 0& -2\cos \mathrm{\&theta;}\sin \mathrm{\&theta;}& 0\\ 0& 1& 0& 0& 0& 0\\ {\cos }^2\mathrm{\&theta;}& 0& {\sin }^2\mathrm{\&theta;}& 0& 2\sin \mathrm{\&theta;}\cos \mathrm{\&theta;}& 0\\ 0& 0& 0& -\sin \mathrm{\&theta;}& 0& -\cos \mathrm{\&theta;}\\ \cos \mathrm{\&theta;}\sin \mathrm{\&theta;}& 0& -\sin \mathrm{\&theta;}\cos \mathrm{\&theta;}& 0& {\sin }^2\mathrm{\&theta;}-{\cos }^2\mathrm{\&theta;}& 0\\ 0& 0& 0& \cos \mathrm{\&theta;}& 0& -\sin \mathrm{\&theta;}\end{array}\right)---\left(2\right)$

C(θ)＝M(θ)×C₀×M(θ)^T (3)

其中，θ是纤维铺放角度，M为Bond变换矩阵。

以45°纤维铺层为例，其弹性刚度矩阵值可由下式计算得到：

$\begin{array}{c}C\left(\frac{\mathrm{\&pi;}}{4}\right)=\left(\begin{array}{cccccc}\frac{1}{2}& 0& \frac{1}{2}& 0& -1& 0\\ 0& 1& 0& 0& 0& 0\\ \frac{1}{2}& 0& \frac{1}{2}& 0& 1& 0\\ 0& 0& \frac{1}{2}& -\frac{\sqrt{2}}{2}& 0& -\frac{\sqrt{2}}{2}\\ \frac{1}{2}& 0& -\frac{1}{2}& 0& 0& 0\\ 0& 0& 0& \frac{\sqrt{2}}{2}& 0& -\frac{\sqrt{2}}{2}\end{array}\right)\&times;\left(\begin{array}{cccccc}126.83& 6.39& 6.39& 0& 0& 0\\ 6.39& 12.88& 6.94& 0& 0& 0\\ 6.39& 6.94& 12.88& 0& 0& 0\\ 0& 0& 0& 2.97& 0& 0\\ 0& 0& 0& 0& 4.85& 0\\ 0& 0& 0& 0& 0& 4.85\end{array}\right)\\ \&times;\left(\begin{array}{cccccc}\frac{1}{2}& 0& \frac{1}{2}& 0& \frac{1}{2}& 0\\ 0& 1& 0& 0& 0& 0\\ \frac{1}{2}& 0& \frac{1}{2}& \frac{1}{2}& -\frac{1}{2}& 0\\ 0& 0& 0& -\frac{\sqrt{2}}{2}& 0& \frac{\sqrt{2}}{2}\\ -1& 0& 1& 0& 0& 0\\ 0& 0& 0& -\frac{\sqrt{2}}{2}& 0& -\frac{\sqrt{2}}{2}\end{array}\right)=\left(\begin{array}{cccccc}42.94& 33.27& 6.67& 0& 0& -28.49\\ 33.27& 42.97& 6.67& 0& 0& -28.49\\ 6.67& 6.67& 12.88& 0& 0& 0.28\\ 0& 0& 0& 3.91& -0.94& 0\\ 0& 0& 0& -0.94& 3.91& 0\\ -28.49& -28.49& 0.28& 0& 0& 31.73\end{array}\right)\left(GPa\right)\end{array}---\left(4\right)$

同理，分别求得-45°及90°纤维铺层对应的弹性刚度矩阵值如下：

$C(-\frac{\mathrm{\&pi;}}{4})=\left(\begin{array}{cccccc}42.97& 33.27& 6.67& 0& 0& 28.49\\ 33.27& 42.97& 6.67& 0& 0& 28.49\\ 6.67& 6.67& 12.88& 0& 0& -0.28\\ 0& 0& 0& 3.91& 0.94& 0\\ 0& 0& 0& 0.94& 3.91& 0\\ 28.49& 28.49& -0.28& 0& 0& 31.73\end{array}\right)\left(GPa\right)---\left(5\right)$

$C\left(\frac{\mathrm{\&pi;}}{2}\right)=\left(\begin{array}{cccccc}12.88& 6.39& 6.94& 0& 0& 0\\ 6.39& 126.83& 6.39& 0& 0& 0\\ 6.94& 6.39& 12.88& 0& 0& 0\\ 0& 0& 0& 4.85& 0& 0\\ 0& 0& 0& 0& 2.97& 0\\ 0& 0& 0& 0& 0& 4.85\end{array}\right)\left(GPa\right)---\left(6\right)$

三、分别将步骤二计算所得0°、±45°及90°纤维铺层的弹性刚度矩阵值及密度值代入Christoffel方程，求解得到各自对应的准纵波群速度慢度曲线(如图1示)。然后，借助MATLAB科学计算软件的多项式拟合功能，构造出准纵波群速度值v(km/s)关于声传播方向角的函数式：

其中N表示拟合的阶数，a(n)是拟合多项式的系数序列。

由图1可知，90°纤维铺层为横观各向同性平面，准纵波群速度值恒等于2.89km/s；同时，45°与-45°纤维铺层的准纵波群速度慢度曲线相同。因此，仅需对0°及45°铺层进行多项式拟合。这里对范围内的数据点进行拟合，阶数N＝6，得到函数式如下：

四、建立CFRP复合材料层合板试样的笛卡尔坐标系，以超声波源点为坐标原点O，水平向右为x轴正向，深度方向为y轴正向。定义声传播方向角为x轴逆时针方向旋转的角度。以A(0.6,0)、B(0,0.6)及C(0.6,0.6)三个目标点为例，超声波由源点O分别传播至目标点A、B及C的路径及所需时间，可通过执行如下步骤获得：

1.将源点O与目标点C之间的矩形计算区域离散化，其中y轴方向网格化间距等于单铺层厚度0.15mm，x轴方向网格化间距等于0.1mm，如图2所示；

2.判断源点与目标点之间的相对位置关系。显然，A点位于O点右侧同一水平位置，B点位于O点正下方，C点位于O点右下侧。因此，超声波由O点分别传播至A和B点的路径即为两点间的连线；

3.超声波由源点O分别传播至目标点A和B的路径和声时如表1所示：

表1

而对于目标点C，需采用Dijkstra最短路径搜索算法进行声线示踪。具体执行如下步骤：

(1)指定节点之间的连接关系，遵循以下原则：当前节点只能与位于其右侧的相邻下层节点直接连接；相对位置关系满足图3所示“×”标识的节点组合间不能直接连接。以1#节点为例，其邻接节点包括7/12/17/22/27/32#节点；

(2)超声波在邻接节点组合间传播的路径为两点间的连线，声时值可由两点间直线距离与对应传播方向上准纵波群速度的比值求得。因此，超声波由1#节点传播至其邻接7/12/17/22/27/32#节点的路径和声时如表2所示：

表2

(3)借助MATLAB科学计算软件，采用Dijkstra最短路径搜索算法计算超声波由源点O传播至目标点C的路径及声时，遵循以下搜索原则：

①将源点O(即1#节点)指定为起始节点，因此其对应的声时值等于0μs，同时指定超声波由源点O传播至其它所有节点的声时为∞；

②将所有节点状态标记为“未访问”，并分配至“未访问节点集合”S；

③以1#起始节点为当前节点，考察超声波经当前节点传播至其邻接节点的路径，相应的声时值为当前节点的声时加上超声波在邻接节点组合间传播的声时。若新路径的声时值小于原有值，则以此新路径代替对应邻接节点的原有路径；反之，则保持路径和声时信息不变。完成遍历后，将当前节点的状态改为“已访问”，并从集合S中移除。据此，1#当前节点的邻接7/12/17/22/27/32#节点路径、声时及访问状态变化情况如表3所示：

表3

④比较步骤③所得最新路径下的声时值，将集合S中具有最小声时值的节点作为新的起始节点。由表3可知，“未访问”节点中，7#节点的声时值为0.060μs的最小值，因此以7#节点为新的起始节点；

⑤重复步骤①～④，直至集合S为空集，搜索结束。超声波由源点O传播至目标点C的路径及声时如表4所示，声线示踪结果如图4所示。

表4

源点	目标点	路径	声时t/μs
				O	C	{1,17,23,29,35}	0.248

Claims (1)

1.一种纤维增强树脂基复合材料声线示踪算法，其特征在于：以单铺层为单元对计算区域分区，并利用弹性刚度矩阵及其旋转变换，实现FRP复合材料弹性特性空间分布的定量描述；将弹性刚度矩阵代入Christoffel方程求解，获得不同纤维取向铺层对应的准纵波群速度值关于传播方向角的函数关系式；计算区域网格化，利用Dijkstra最短路径搜索算法，搜寻超声波由源点传播至目标点所经过的路径并计算对应声时，具体执行如下步骤：

一、将FRP复合材料计算区域以铺层为单元划分为一系列分区，各分区等效为均质弹性各向异性介质，并以弹性刚度矩阵描述分区的弹性特性；

二、以0°纤维取向铺层的弹性刚度矩阵C₀为参考值，其它纤维取向铺层的弹性刚度矩阵值C(θ)由C₀进行Bond变换求得：

$C_0=\left(\begin{array}{cccccc}{c}_{11}& c_{12}& c_{12}& & & \\ c_{12}& c_{22}& c_{22}-2c_{44}& & & \\ c_{12}& c_{22}-2c_{44}& c_{22}& & & \\ & & & c_{44}& & \\ & & & & c_{55}& \\ & & & & & c_{55}\end{array}\right)\left(GPa\right)---\left(1\right)$

$M\left(\mathrm{\&theta;}\right)=\left(\begin{array}{cccccc}{\sin }^2\mathrm{\&theta;}& 0& {\cos }^2\mathrm{\&theta;}& 0& -2\cos \mathrm{\&theta;}\sin \mathrm{\&theta;}& 0\\ 0& 1& 0& 0& 0& 0\\ {\cos }^2\mathrm{\&theta;}& 0& {\sin }^2\mathrm{\&theta;}& 0& 2\sin \mathrm{\&theta;}\cos \mathrm{\&theta;}& 0\\ 0& 0& 0& -\sin \mathrm{\&theta;}& 0& -\cos \mathrm{\&theta;}\\ \cos \mathrm{\&theta;}\sin \mathrm{\&theta;}& 0& -\sin \mathrm{\&theta;}\cos \mathrm{\&theta;}& 0& {\sin }^2\mathrm{\&theta;}-{\cos }^2\mathrm{\&theta;}& 0\\ 0& 0& 0& \cos \mathrm{\&theta;}& 0& -\sin \mathrm{\&theta;}\end{array}\right)---\left(2\right)$

C(θ)＝M(θ)×C₀×M(θ)^T (3)

其中，θ是纤维铺放角度，M为Bond变换矩阵；

三、将弹性刚度矩阵值C(θ)以及试样密度值ρ代入Christoffel方程，求解得到不同纤维取向铺层准纵波群速度的慢度曲线；

四、借助MATLAB科学计算软件，构建各纤维取向铺层准纵波群速度值关于传播方向角的函数式

其中N表示拟合的阶数，a(n)是拟合多项式的系数序列；

五、建立FRP复合材料层合板试样的笛卡尔坐标系，以超声波源点为坐标原点O，水平向右为x轴正向，深度方向为y轴正向；定义声传播方向角为x轴逆时针方向旋转的角度；则超声波由源点O传播至目标点E的路径及所需时间，可通过执行如下步骤获得：

1.将源点O与目标点E之间的矩形计算区域离散化，其中y轴方向网格化间距等于单铺层厚度δ，x轴方向网格化间距需根据计算精度要求，选择小于单铺层厚度范围内的值；

2.判断源点O(x_O,y_O)与目标点E(x_E,y_E)之间的相对位置关系：

表1

其中，

Δx_OE＝x_E-x_O,Δy_OE＝y_E-y_O (5)

$d_{OE}=\sqrt{{(x_E-x_O)}^2+{(y_E-y_O)}^2}---\left(6\right)$

而当Δy_OE≠0且Δx_OE≠0时，需采用Dijkstra最短路径搜索算法进行声线示踪；此外，当Δy_OE<0时，源点与目标点坐标值互换；

3.在Δy_OE＞0且Δx_OE＞0时，采用Dijkstra最短路径搜索算法进行声线示踪的具体执行步骤如下：

(1)指定节点之间的连接关系：与当前节点i(x_i,y_i)能够直接连接的节点j(x_j,y_j)间的相对位置关系需满足如下条件：

Δx_ij＝x_j-x_i>0且Δy_ij＝y_j-y_i＝δ (8)

(2)计算超声波在邻接节点组合间传播所需时间，在邻接节点i(x_i,y_i)和j(x_j,y_j)，首先根据二节点的坐标值，计算得到其连线对应的传播方向角根据连线所经过铺层的慢度曲线，计算得到准纵波群速度值则超声波由节点i(x_i,y_i)传播至j(x_j,y_j)所需时间

(3)借助MATLAB科学计算软件，采用Dijkstra最短路径搜索算法计算超声波由源点O传播至目标点E的路径及声时，遵循以下搜索原则：

①将源点O指定为起始节点，因此其对应的声时值等于0μs，同时指定超声波由源点O传播至其它所有节点的声时为∞；

②将所有节点状态标记为“未访问”，并分配至“未访问节点集合”S；

③以起始节点O为当前节点，考察超声波经当前节点传播至其邻接节点的路径，相应的声时值为当前节点的声时加上超声波在邻接节点组合间传播的声时；若新路径的声时值小于原有值，则以此新路径代替对应邻接节点的原有路径；反之，则保持路径和声时信息不变；完成遍历后，将当前节点的状态改为已访问，并从集合S中移除；

④比较步骤③所得最新路径下的声时值，将集合S中具有最小声时值的节点作为新的起始节点；

⑤重复步骤①～④，直至集合S为空集，搜索结束；输出超声波由源点O传播至目标点E所经过路径的声线图像及声时值。