一种基于超声导波的适用于任意形状截面的应力检测方法

1.本发明涉及无损检测技术领域，特别是指一种基于超声导波的适用于任意形状截面的应力检测方法。


背景技术：

2.热膨胀或收缩、扩散、相变、轧制、拉伸、锻造和焊接等过程会在零件中引入残余应力。残余应力的常见的有害影响包括加工零件的翘曲、拉制产品的开裂、淬火裂纹、疲劳引起的过早失效和应力腐蚀开裂。对结构和机械部件的应力状态进行监测，可提供有关载荷模式、结构完整性和材料退化的宝贵信息。对于结构安全、估计结构或部件的剩余使用寿命以及降低维护和维修成本至关重要。施加应力或残余应力的监测对于工业零件制造的质量控制和结构的健康管理非常重要。
3.目前，材料的应力测量可以简略分为有损检测(机械方法)和无损检测两大类。有损检测往往要对检测的试件进行机械切割或分离，从而根据测量出的应变计算出残余应力。常用的方法有环芯法、钻孔法、取条法、剥层法，这些方法理论完善、技术成熟、精度较高，但都要对工件造成不同程度的破坏。且机械方法复杂，不适合在线监测。无损检测即物理检测法正好解决了此方面的问题，不需要对被测试件进任何破坏即可测试出残余应力。多年以来，大量科研人员探索出许多无损检测的方法，其中主要有x射线法、中子衍射法、磁性法、超声波法、扫描电子声显微镜法、电子散班干涉法等。但这些方法往往成本较高、便携性低、对人体有害、适用材料的局限性高。与其它检测方法相比，超声导波受距离限制小、灵敏度高、方向性好，是这些无损检测中最具有发展前景的一种方法。
4.与传统超声检测相比，超声导波技术可在较大范围进行检测。一种结合了经典的声弹性导波理论与半解析有限元的计算方法，可以求解受轴向力的模型在各个方向的频散曲线。一种基于声弹性引起的相速度变化与导波的传播角之间的正弦依赖性的应力反演算法，可以检测出任意形状横截面的双轴力的大小与方向。一种基于多个频段下相速度的变化的反演算法。可以检测出双轴力的大小。一种结合频散曲线、半解析有限元和最小二乘法的超声导波检测方法，可以定量的对应力原位检测。


技术实现要素：

5.针对上述背景技术中存在的不足，本发明提出了一种基于超声导波的适用于任意形状截面的应力检测方法，解决了现有无损检测方法成本较高、测量范围小、检测精度低，且现有的超声波检测方法通用性较差的技术问题。
6.本发明的技术方案是这样实现的：
7.一种基于超声导波的适用于任意形状截面的应力检测方法，其步骤如下：
8.步骤一：使用声弹性导波理论与半解析有限元相结合的方法求解声弹性导波的控制方程，并将声弹性导波的控制方程应用于商业软件comsol multiphysics中，分别获得单轴受力模型、双轴受力模型和自由状态模型的频散曲线；
9.步骤二：将单轴受力模型的频散曲线与自由状态模型的频散曲线做差获得单轴的相速度变化量；将双轴受力模型的频散曲线与自由状态模型的频散曲线做差获得双轴的相速度变化量；
10.步骤三：根据单轴的相速度变化量计算单轴受力模型在多个传播角度的相速度变化的计算值，将相速度变化的计算值与应力真实值作为样本以通过最小二乘法拟合出系数k1、k2；
11.步骤四：将获得的系数k1、k2与双轴的相速度变化量输入最小二乘法，进行多角度反演与多频反演，输出双轴应力的大小与方向。
12.优选地，在步骤一中：
13.在声弹性导波理论中，研究对象存在三种状态，包括：未产生变形时的自然状态、在预应力作用下产生了变形的初始状态、存在弹性波传播的变形状态；从自然状态到初始状态、从自然状态到变形状态以及从初始状态到变形状态的变形(ui,uf,u)分别表示为：
14.ui＝x-ξ
15.uf＝x-ξ
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)；
16.u＝x-x＝u
f-ui17.其中，ui表示从自然状态到初始状态的变形，x表示质点初始状态下的位置向量，ξ表示质点自然状态下的位置向量，uf表示从自然状态到变形状态的变形，x表示质点变形状态下的位置向量，u表示从初始状态到变形状态的变形；
18.初始状态的lagrangian应变张量为：其中，为初始状态下的lagrangian应变张量，均表示从自然状态到初始状态的变形；ξ
i'
、ξj均表示质点自然状态下的位置向量；i'、j、k、l、m、n、o、p均表示索引；
19.当研究对象的材料为超弹性时，green-lagrange应变张量通过本构方程与第二piola-kirchoff应力张量相关，只保留二阶和三阶弹性常数，则：其中，为应力张量；均表示lagrangian应变张量，c
i'jkl
表示二阶弹性常数，c
i'jklmn
表示三阶弹性常数；如果初始应变很小，进一步简化方程：
[0020][0021]
其中，c
ljop
表示二阶弹性常数，为柯西应变张量：均表示从自然状态到初始状态的变形，ξ
p
、ξo均表示质点自然状态下的位置向量；
[0022]
分别以双轴力方向x
i”和导波传播方向为主轴方向建立的两个坐标系，在自然状态下，在坐标系中的初始应力张量可表示为：
[0023]
[0024]
其中，t为初始应力张量，σ
11
表示沿方向的主轴力，σ
22
表示沿方向的主轴力；
[0025]
通过导波传播角度的旋转变换，在x
i”坐标系中应力张量表示为：
[0026][0027]
其中，均表示x
i”轴和轴之间夹角的余弦值，t
mn
表示坐标系中特定方向的初始应力张量；
[0028]
声弹性波相对于初始状态的增量位移方程为：
[0029][0030]
其中，ρi是初始状态的密度，u
i'
、uk均表示从初始状态到变形状态的变形，x
l
、xj、x
p
均表示质点初始状态下的位置向量，c
i'jkl
、c
i'jpl
、c
pjkl
均表示二阶弹性常数，c
i'jklmn
表示表示三阶弹性常数，t表示时间；
[0031]
将式(2)代入等式(5)，可以得到：
[0032][0033]
其中c
ljno
、c
i'mkl
、c
mjkl
、c
i'jml
、c
i'jkm
均表示二阶弹性常数，c
i'jklmn
表示三阶弹性常数，均表示初始状态下的柯西应变张量，δ
i'k
表示克罗内克函数，表示从自然状态到初始状态的变形，xm表示质点在初始状态下的位置向量；
[0034]
在边界上的应力张量τ
i'
可以导出为：
[0035][0036]
其中，nj是在初始状态下垂直于边界的向外单位向量；从式(6)和(7)可以看出，预应力对位移增量的影响用声弹性控制方程中的系数γ
i'jkl
来表示；
[0037]
根据半解析有限元算法的特征，假设导波在x3方向上以简谐波的形式传播，波导中的位移可以写为：
[0038][0039]
其中，u
i'
为截面位移；下标i'＝1,2,3；i为虚数单位；ω表示角频率；k'表示x3方向上的波数；位移的导数可以写为：
[0040][0041]
将式(9)代入式(6)并消除相位项，声弹性导波的控制方程可以写为：
[0042][0043]
其中，uk表示截面位移，γ
i'3kl
、γ
i'jk3
、γ
i'3k3
均表示声弹性系数；
[0044]
引入一个新变量pk＝k'uk将式(10)转化为一般线性形式，则在波导域内(inω)：
[0045][0046]
ρiω2δ
i'k
p
k-k'ρiω2δ
i'k
uk＝0 in ω
ꢀꢀꢀꢀ
(12)；
[0047]
其中，ω表示波导域；
[0048]
在波导边界上，τ
i'
的振型s
i'
可以写为：
[0049][0050]
在式(10)-(13)中隐藏了对重复索引k＝1,2,3和j,l＝1,2的求和；
[0051]
通过找到导波在选定角频率ω下的波数k'和模式形状来求解特征值问题；通过选择实波数模式，可以识别出每个解中的传播模式，并且可以通过计算相速度；最后，通过在特定频率上求解特征值问题，并在每个频率步长上组合具有最相似模式形状的模式来获得频散曲线。
[0052]
优选地，在步骤三中：
[0053]
将声弹性引起的相速度变化与导波传播角度之间的正弦依赖性扩展为适用于任意形状的截面；双轴应力引起的相速度变化与导波传播角度之间的正弦关系表示为：
[0054][0055]
其中，
△cp
(
·
)表示相速度；由式(14)，根据单轴受力模型在多个传播角度的相速度变化，通过最小二乘法拟合出不同频率下对应的系数k1、k2。
[0056]
优选地，在步骤四中：
[0057]
将坐标系x
i”顺时针旋转角度θ获得坐标系x
i”'作为测试坐标系；因此，方向与x1'方向之间夹角满足在坐标系x
i”'中，相速度
△cp
的变化为：
[0058]
[0059]
其中，系数a0、a1和a2分别表示为：
[0060][0061]
根据双轴的相速度变化量以及常数k1和k2，通过最小二乘法确定式(15)中的系数a0、a1和a2；再通过转换等式(15)得到多角度反演公式，计算施加的应力σ
11
、σ
22
和角度θ：
[0062][0063]
通过转换等式(15)对多角度频散数据的计算需求转换为多个频段频散的计算需求，则可以得到轴向力大小的多频反演公式：
[0064][0065]
其中，b1、b2均为系数，且
[0066]
与现有技术相比，本发明产生的有益效果为：本发明采用两个声弹性系数来建立频散信息与应力检测结果之间的关系，利用相速度变化量与应力之间的正弦关系，通过两次系数拟合得到双轴力的大小与方向，即实现轴向应力的定量评估；本发明的反演算法简单，网络架构搭建容易；适用性强、精度高；适用于任意形状截面的应力的高精度快速定量检测。
附图说明
[0067]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0068]
图1为本发明的流程图。
[0069]
图2为本发明中的二维半解析有限元模型几何图。
[0070]
图3为本发明中的模型坐标系示意图。
[0071]
图4为本发明中s0模式相速度变化量结果图。
[0072]
图5为本发明中a0模式相速度变化量结果图。
[0073]
图6为本发明中相速度变化量与传播角正弦关系拟合结果图。
[0074]
图7为本发明中多角度反演结果图。
[0075]
图8为本发明中多频反演结果图。
具体实施方式
[0076]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有付出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0077]
如图1所示，本发明实施例提供了一种基于超声导波的适用于任意形状截面的应力检测方法，具体步骤如下：
[0078]
步骤一：将经典的声弹性导波理论与半解析有限元(semi-analytical finite element,safe)结合，使用safe方法求解声弹性导波的控制方程，并将声弹性导波的控制方程应用于商业软件comsol multiphysics中，分别获得单轴受力模型、双轴受力模型和自由状态模型在多个方向的频散曲线。如图2所示构建二维几何模型来计算轴向受力的模型的声场。模型的宽为0.64mm，厚度为3.2mm，在二维模型两侧耦合周期边界条件(periodic boundary condition,pbc)可以消除边界对模型的限制，模拟导波在无限宽板中的传播，这里以无限宽铝板为例：
[0079]
在经典的声弹性导波理论中，研究对象存在三种状态，包括：未产生变形时的自然状态、在预应力作用下产生了变形的初始状态、存在弹性波传播的变形状态(即最终状态)；同一个质点在三种状态下分别通过位置向量“ξ”、“x”、和“x”表示。三种状态下的所有物理变量和材料属性分别用上标“0”、“i”和“f”表示。从自然状态到初始状态、从自然状态到变形状态以及从初始状态到变形状态的变形分别表示为：
[0080]
ui＝x-ξ
[0081]
uf＝x-ξ
ꢀꢀꢀꢀꢀꢀꢀꢀ
(1)；
[0082]
u＝x-x＝u
f-ui[0083]
其中，ui表示从自然状态到初始状态的变形，x表示质点初始状态下的位置向量，ξ表示质点自然状态下的位置向量，uf表示从自然状态到变形状态的变形，x表示质点变形状态下的位置向量，u表示从初始状态到变形状态的变形。
[0084]
初始状态的lagrangian应变张量为：其中，为初始状态下的lagrangian应变张量，均表示从自然状态到初始状态的变形；ξ
i'
、ξj均表示质点自然状态下的位置向量；i'、j、k、l、m、n、o、p均表示索引。
[0085]
当研究对象的材料为超弹性时，green-lagrange应变张量通过本构方程与第二piola-kirchoff应力张量相关，只保留二阶和三阶弹性常数，则：其中，为应力张量；均表示lagrangian应变张量，c
i'jkl
表示二阶弹性常数，c
i'jklmn
表示三阶弹性常数；如果初始应变很小，进一步简化方程：
[0086]
[0087]
其中，c
ljop
表示二阶弹性常数，为柯西应变张量：均表示从自然状态到初始状态的变形，ξ
p
、ξo均表示质点自然状态下的位置向量。
[0088]
如图3所示，为了获得在双轴力作用下沿任意角度的导波传播，分别以双轴力方向x
i”和导波传播方向为主轴方向建立的两个坐标系，在自然状态下，在坐标系中的初始应力张量可表示为：
[0089][0090]
其中，t为初始应力张量，σ
11
表示沿方向的主轴力，σ
22
表示沿方向的主轴力。
[0091]
通过导波传播角度的旋转变换，在x
i”坐标系中应力张量表示为：
[0092][0093]
其中，均表示x
i”轴和轴之间夹角的余弦值，t
mn
表示坐标系中特定方向的初始应力张量。
[0094]
声弹性波相对于初始状态的增量位移方程为：
[0095][0096]
其中，ρi是初始状态的密度，u
i'
、uk均表示从初始状态到变形状态的变形，x
l
、xj、x
p
均表示质点初始状态下的位置向量，c
i'jkl
、c
i'jpl
、c
pjkl
均表示二阶弹性常数，c
i'jklmn
表示表示三阶弹性常数，t表示时间；将式(2)代入等式(5)，可以得到：
[0097][0098]
其中c
ljno
、c
i'mkl
、c
mjkl
、c
i'jml
、c
i'jkm
均表示二阶弹性常数，c
i'jklmn
表示三阶弹性常数，表示三阶弹性常数，均表示初始状态下的柯西应变张量，δ
i'k
表示克罗内克函数，表示从自然状态到初始状态的变形，xm表示质点在初始状态下的位置向量。
[0099]
在边界上的应力张量τ
i'
可以导出为：
[0100][0101]
其中，nj是在初始状态下垂直于边界的向外单位向量；从式(6)和(7)可以看出，预应力对位移增量的影响用声弹性控制方程中的系数γ
i'jkl
来表示。
[0102]
根据半解析有限元算法的特征，假设导波在x3方向上以简谐波的形式传播，波导
中的位移可以写为：
[0103][0104]
其中，u
i'
为截面位移；下标i'＝1,2,3；i为虚数单位；ω表示角频率；k'表示x3方向上的波数；位移的导数可以写为：
[0105][0106]
将式(9)代入式(6)并消除相位项，声弹性导波的控制方程可以写为：
[0107][0108]
其中，uk表示截面位移，γ
i'3kl
、γ
i'jk3
、γ
i'3k3
均表示声弹性系数。
[0109]
为了将该方法应用于商用有限元软件包(comsol multiphysics)，需要对声弹性导波的控制方程进行改写，使其表达形式与软件中特征值问题求解的输入公式的一般表达式相同。为此引入一个新变量pk＝k'uk将式(10)转化为一般线性形式，则在波导域内(inω)：
[0110][0111]
ρiω2δ
i'k
p
k-k'ρiω2δ
i'k
uk＝0 in ω
ꢀꢀ
(12)；
[0112]
其中，ω表示波导域。
[0113]
在边界上，τ
i'
的振型s
i'
可以写为：
[0114][0115]
在式(10)-(13)中隐藏了对重复索引k＝1,2,3和j,l＝1,2的求和；
[0116]
通过找到导波在选定角频率ω下的波数k'和模式形状来求解特征值问题；通过选择实波数模式，可以识别出每个解中的传播模式，并且可以通过计算相速度；最后，通过在特定频率上求解特征值问题，并在每个频率步长上组合具有最相似模式形状的模式来获得频散曲线。根据上述方法计算单轴、双轴受力模型(σ
11
和σ
22
的范围均为0到100mpa，步长为2.5mpa)与自由状态的模型在多个方向(的范围为0
°
到90
°
，步长为2.5
°
)的频散曲线。
[0117]
步骤二：将单轴受力模型的频散曲线与自由状态模型的频散曲线做差获得单轴的相速度变化量；将双轴受力模型的频散曲线与自由状态模型的频散曲线做差获得双轴的相速度变化量；图4和图5分别为3.2mm厚铝板在不同单轴力下，s0模式和a0模式在导波传播角度分别为0
°
、30
°
、60
°
和90
°
时的相速度变化量随频率的变化。在任何选定的频率下，两种
模式的相速度变化都随应力线性变化。在相同应力下，在1mhz以下频率范围内s0模式的相速度变化比a0模式的相速度更明显，两种模式的相速度变化均在高频范围内趋于平坦。显著而稳定的相速度变化往往意味着对应力更高的敏感性，以及更好的反演结果。因此相速度变化量的绝对值大小，是应力反演中频率段选择的重要参考因素。
[0118]
步骤三：根据单轴的相速度变化量计算单轴受力模型在多个传播角度的相速度变化的计算值，将相速度变化的计算值与应力真实值作为样本以通过最小二乘法拟合出系数k1、k2；
[0119]
将声弹性引起的相速度变化与导波传播角度之间的正弦依赖性扩展为适用于任意形状的截面；双轴应力引起的相速度变化与导波传播角度之间的正弦关系表示为：
[0120][0121]
其中，
△cp
(
·
)表示相速度；公式(14)用xi坐标系表示，即施加的轴向应力平行于坐标轴。根据s0、a0和a1模式各自的相速度变化曲线，不同模式采用不同频段进行反演：s0模式为100-1900khz，a0模式为100-1900khz，a1模式为600-2400khz，步长均为4khz。由公式(14)，根据单轴受力模型在多个传播角度的相速度变化，通过最小二乘法拟合出不同频率下对应的系数k1、k2。图6为单轴受力100mpa时，相速度变化量与传播角的正弦关系拟合结果图。对于s0，k1和k2分别为-4.63
×
10-7
和7.31
×
10-8
；对于a1，k1和k2分别为-7.54
×
10-7
和1.98
×
10-7
。
[0122]
步骤四：将获得的系数k1、k2与双轴的相速度变化量输入最小二乘法，进行多角度反演与多频反演，输出双轴应力的大小与方向。
[0123]
将坐标系x
i”顺时针旋转角度θ获得坐标系x
i”'作为测试坐标系；因此，方向与x1'方向之间夹角满足在实践中双轴力的方向，即θ通常是未知的，而随变化。在坐标系x
i”'中，相速度
△cp
的变化为：
[0124][0125]
其中，系数a0、a1和a2分别表示为：
[0126][0127]
根据双轴的相速度变化量以及常数k1和k2，通过最小二乘法确定式(15)中的系数a0、a1和a2；再通过转换等式(15)得到多角度反演公式，计算施加的应力σ
11
、σ
22
和角度θ：
[0128][0129]
通过转换等式(15)对多角度频散数据的计算需求转换为多个频段频散的计算需求，则可以得到轴向力大小的多频反演公式：
[0130][0131]
其中，b1、b2均为系数，且
[0132]
输出的反演结果与真实值之间的误差用均方误差(mean square error，mse)表示。
[0133]
图7为本发明中多角度反演结果图。反演时，角度的范围为10到80度(步长为10度)。由于双轴力σ
11
和σ
22
是等价的，所以图7(a)中只显示了σ
11
的估计误差，图7(b)为θ的误差估计。在1mhz以下，虽然只选择了8个角度的数据作为已知量进行反演，三种模式对σ
11
的反演误差都很小，而s0模式的低频段和a0模式的高频段之间的误差是最低的。各模态对角度θ的反演精度高且稳定。
[0134]
图8为本发明中多频反演结果图。图8(a)、(b)、(c)依次为s0、a0、a1模式的多频反演结果，三个模式各自选择的频率区间可以参考多角度反演的结果：s0模式为100-500khz，a0模式为1500-1900khz，a1模式为800-1200khz，步长均为80khz。当角度取一定值取一定值时，三种模式的σ
11
和σ
22
的误差分别达到最小值：对于s0，为22.5度时σ
11
误差为0.3485％，为67.5度时σ
22
的误差为0.3561％；对于a0，为32.5度时σ
11
的误差为0.5046％，为57.5度时σ
22
的误差为0.5219％；对于a1，为27.5度时σ
11
的误差为2.475％，为62.5度时σ
22
的误差为2.737％。从图8可以看出，无论是哪个模式，σ
11
和σ
22
的反演精度基本对称相同，且和互余。在σ
11
的反演精度达到最高时，σ
22
的反演精度并不理想，说明在实际应用中，需要根据所选模式分别采集两个特定角度的数据。实际应用中，多频反演法结合扫频法，可以实现利用两个方向的测量来确定双轴力的大小，这也是一种在实际应用中很有前途的应力检测方法。
[0135]
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。