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

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

Info

Publication number
CN114739546A
CN114739546A CN202210564728.1A CN202210564728A CN114739546A CN 114739546 A CN114739546 A CN 114739546A CN 202210564728 A CN202210564728 A CN 202210564728A CN 114739546 A CN114739546 A CN 114739546A
Authority
CN
China
Prior art keywords
stress
state
phase velocity
initial state
equation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202210564728.1A
Other languages
English (en)
Other versions
CN114739546B (zh
Inventor
刘洋
赵春雨
李健
曾周末
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Tianjin University
Original Assignee
Tianjin University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Tianjin University filed Critical Tianjin University
Priority to CN202210564728.1A priority Critical patent/CN114739546B/zh
Publication of CN114739546A publication Critical patent/CN114739546A/zh
Application granted granted Critical
Publication of CN114739546B publication Critical patent/CN114739546B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01L—MEASURING FORCE, STRESS, TORQUE, WORK, MECHANICAL POWER, MECHANICAL EFFICIENCY, OR FLUID PRESSURE
    • G01L1/00—Measuring force or stress, in general
    • G01L1/25—Measuring force or stress, in general using wave or particle radiation, e.g. X-rays, microwaves, neutrons
    • G01L1/255—Measuring force or stress, in general using wave or particle radiation, e.g. X-rays, microwaves, neutrons using acoustic waves, or acoustic emission

Landscapes

  • Physics & Mathematics (AREA)
  • Acoustics & Sound (AREA)
  • Health & Medical Sciences (AREA)
  • Toxicology (AREA)
  • General Physics & Mathematics (AREA)
  • Investigating Or Analyzing Materials By The Use Of Ultrasonic Waves (AREA)

Abstract

本发明提出了一种基于超声导波的适用于任意形状截面的应力检测方法，其步骤为：首先，根据半解析有限元结合声弹性理论进行正演模拟获取频散信号；其次，分别将单轴、双轴受力模型的频散曲线与自由状态的模型的频散曲线做差获得各自的相速度变化量；然后，将单轴受力模型的相速度变化量与应力真实值作为样本进行数据预处理以拟合系数；最后，将预处理后获得的系数与双轴受力模型的相速度变化量输入最小二乘法进行训练与验证，模型预测输出即为双轴应力的大小与方向预测结果。本发明实现了应力的大小与方向的定量评估；且本发明的反演算法简单，架构搭建容易；检测精度高；适用于任意截面形状的研究对象的应力的高精度快速定量检测。

Description

一种基于超声导波的适用于任意形状截面的应力检测方法
技术领域
本发明涉及无损检测技术领域，特别是指一种基于超声导波的适用于任意形状截面的应力检测方法。
背景技术
热膨胀或收缩、扩散、相变、轧制、拉伸、锻造和焊接等过程会在零件中引入残余应力。残余应力的常见的有害影响包括加工零件的翘曲、拉制产品的开裂、淬火裂纹、疲劳引起的过早失效和应力腐蚀开裂。对结构和机械部件的应力状态进行监测，可提供有关载荷模式、结构完整性和材料退化的宝贵信息。对于结构安全、估计结构或部件的剩余使用寿命以及降低维护和维修成本至关重要。施加应力或残余应力的监测对于工业零件制造的质量控制和结构的健康管理非常重要。
目前，材料的应力测量可以简略分为有损检测(机械方法)和无损检测两大类。有损检测往往要对检测的试件进行机械切割或分离，从而根据测量出的应变计算出残余应力。常用的方法有环芯法、钻孔法、取条法、剥层法，这些方法理论完善、技术成熟、精度较高，但都要对工件造成不同程度的破坏。且机械方法复杂，不适合在线监测。无损检测即物理检测法正好解决了此方面的问题，不需要对被测试件进任何破坏即可测试出残余应力。多年以来，大量科研人员探索出许多无损检测的方法，其中主要有X射线法、中子衍射法、磁性法、超声波法、扫描电子声显微镜法、电子散班干涉法等。但这些方法往往成本较高、便携性低、对人体有害、适用材料的局限性高。与其它检测方法相比，超声导波受距离限制小、灵敏度高、方向性好，是这些无损检测中最具有发展前景的一种方法。
与传统超声检测相比，超声导波技术可在较大范围进行检测。一种结合了经典的声弹性导波理论与半解析有限元的计算方法，可以求解受轴向力的模型在各个方向的频散曲线。一种基于声弹性引起的相速度变化与导波的传播角之间的正弦依赖性的应力反演算法，可以检测出任意形状横截面的双轴力的大小与方向。一种基于多个频段下相速度的变化的反演算法。可以检测出双轴力的大小。一种结合频散曲线、半解析有限元和最小二乘法的超声导波检测方法，可以定量的对应力原位检测。
发明内容
针对上述背景技术中存在的不足，本发明提出了一种基于超声导波的适用于任意形状截面的应力检测方法，解决了现有无损检测方法成本较高、测量范围小、检测精度低，且现有的超声波检测方法通用性较差的技术问题。
本发明的技术方案是这样实现的：
一种基于超声导波的适用于任意形状截面的应力检测方法，其步骤如下：
步骤一：使用声弹性导波理论与半解析有限元相结合的方法求解声弹性导波的控制方程，并将声弹性导波的控制方程应用于商业软件COMSOL Multiphysics中，分别获得单轴受力模型、双轴受力模型和自由状态模型的频散曲线；
步骤二：将单轴受力模型的频散曲线与自由状态模型的频散曲线做差获得单轴的相速度变化量；将双轴受力模型的频散曲线与自由状态模型的频散曲线做差获得双轴的相速度变化量；
步骤三：根据单轴的相速度变化量计算单轴受力模型在多个传播角度的相速度变化的计算值，将相速度变化的计算值与应力真实值作为样本以通过最小二乘法拟合出系数K1、K2；
步骤四：将获得的系数K1、K2与双轴的相速度变化量输入最小二乘法，进行多角度反演与多频反演，输出双轴应力的大小与方向。
优选地，在步骤一中：
在声弹性导波理论中，研究对象存在三种状态，包括：未产生变形时的自然状态、在预应力作用下产生了变形的初始状态、存在弹性波传播的变形状态；从自然状态到初始状态、从自然状态到变形状态以及从初始状态到变形状态的变形(ui,uf,u)分别表示为：
ui＝X-ξ
uf＝x-ξ (1)；
u＝x-X＝uf-ui
其中，ui表示从自然状态到初始状态的变形，X表示质点初始状态下的位置向量，ξ表示质点自然状态下的位置向量，uf表示从自然状态到变形状态的变形，x表示质点变形状态下的位置向量，u表示从初始状态到变形状态的变形；
初始状态的Lagrangian应变张量为：
Figure BDA0003657440860000021
其中，
Figure BDA0003657440860000022
为初始状态下的Lagrangian应变张量，
Figure BDA0003657440860000023
均表示从自然状态到初始状态的变形；ξi'、ξj均表示质点自然状态下的位置向量；i'、j、k、l、m、n、o、p均表示索引；
当研究对象的材料为超弹性时，Green-Lagrange应变张量通过本构方程与第二Piola-Kirchoff应力张量相关，只保留二阶和三阶弹性常数，则：
Figure BDA0003657440860000024
其中，
Figure BDA0003657440860000025
为应力张量；
Figure BDA0003657440860000026
均表示Lagrangian应变张量，ci'jkl表示二阶弹性常数，ci'jklmn表示三阶弹性常数；如果初始应变很小，进一步简化方程：
Figure BDA0003657440860000031
其中，cljop表示二阶弹性常数，
Figure BDA0003657440860000032
为柯西应变张量：
Figure BDA0003657440860000033
均表示从自然状态到初始状态的变形，ξp、ξo均表示质点自然状态下的位置向量；
分别以双轴力方向xi”和导波传播方向
Figure BDA0003657440860000034
为主轴方向建立的两个坐标系，在自然状态下，在
Figure BDA0003657440860000035
坐标系中的初始应力张量可表示为：
Figure BDA0003657440860000036
其中，T为初始应力张量，σ11表示沿
Figure BDA0003657440860000037
方向的主轴力，σ22表示沿
Figure BDA0003657440860000038
方向的主轴力；
通过导波传播角度
Figure BDA0003657440860000039
的旋转变换，在xi”坐标系中应力张量表示为：
Figure BDA00036574408600000310
其中，
Figure BDA00036574408600000311
均表示xi”轴和
Figure BDA00036574408600000312
轴之间夹角的余弦值，Tmn表示
Figure BDA00036574408600000313
坐标系中特定方向的初始应力张量；
声弹性波相对于初始状态的增量位移方程为：
Figure BDA00036574408600000314
其中，
Figure BDA00036574408600000315
ρi是初始状态的密度，ui'、uk均表示从初始状态到变形状态的变形，Xl、Xj、Xp均表示质点初始状态下的位置向量，ci'jkl、ci'jpl、cpjkl均表示二阶弹性常数，ci'jklmn表示表示三阶弹性常数，t表示时间；
将式(2)代入等式(5)，可以得到：
Figure BDA00036574408600000316
其中
Figure BDA00036574408600000317
cljno、ci'mkl、cmjkl、ci'jml、ci'jkm均表示二阶弹性常数，ci'jklmn表示三阶弹性常数，
Figure BDA00036574408600000318
均表示初始状态下的柯西应变张量，δi'k表示克罗内克函数，
Figure BDA00036574408600000319
表示从自然状态到初始状态的变形，Xm表示质点在初始状态下的位置向量；
在边界上的应力张量τi'可以导出为：
Figure BDA0003657440860000041
其中，nj是在初始状态下垂直于边界的向外单位向量；从式(6)和(7)可以看出，预应力对位移增量的影响用声弹性控制方程中的系数Γi'jkl来表示；
根据半解析有限元算法的特征，假设导波在X3方向上以简谐波的形式传播，波导中的位移可以写为：
Figure BDA0003657440860000042
其中，Ui'为截面位移；下标i'＝1,2,3；I为虚数单位；ω表示角频率；k'表示X3方向上的波数；位移的导数可以写为：
Figure BDA0003657440860000043
将式(9)代入式(6)并消除相位项，声弹性导波的控制方程可以写为：
Figure BDA0003657440860000044
其中，
Figure BDA0003657440860000045
Uk表示截面位移，Γi'3kl、Γi'jk3、Γi'3k3均表示声弹性系数；
引入一个新变量Pk＝k'Uk将式(10)转化为一般线性形式，则在波导域内(inΩ)：
Figure BDA0003657440860000046
ρiω2δi'kPk-k'ρiω2δi'kUk＝0 in Ω (12)；
其中，Ω表示波导域；
在波导边界
Figure BDA0003657440860000047
上，τi'的振型Si'可以写为：
Figure BDA0003657440860000048
在式(10)-(13)中隐藏了对重复索引k＝1,2,3和j,l＝1,2的求和；
通过找到导波在选定角频率ω下的波数k'和模式形状来求解特征值问题；通过选择实波数模式，可以识别出每个解中的传播模式，并且可以通过
Figure BDA0003657440860000051
计算相速度；最后，通过在特定频率上求解特征值问题，并在每个频率步长上组合具有最相似模式形状的模式来获得频散曲线。
优选地，在步骤三中：
将声弹性引起的相速度变化与导波传播角度之间的正弦依赖性扩展为适用于任意形状的截面；双轴应力引起的相速度变化与导波传播角度
Figure BDA0003657440860000052
之间的正弦关系表示为：
Figure BDA0003657440860000053
其中，△cp(·)表示相速度；由式(14)，根据单轴受力模型在多个传播角度的相速度变化，通过最小二乘法拟合出不同频率下对应的系数K1、K2。
优选地，在步骤四中：
将坐标系xi”顺时针旋转角度θ获得坐标系xi”'作为测试坐标系；因此，
Figure BDA0003657440860000054
方向与x1'方向之间夹角
Figure BDA0003657440860000055
满足
Figure BDA0003657440860000056
在坐标系xi”'中，相速度△cp的变化为：
Figure BDA0003657440860000057
其中，系数a0、a1和a2分别表示为：
Figure BDA0003657440860000058
根据双轴的相速度变化量以及常数K1和K2，通过最小二乘法确定式(15)中的系数a0、a1和a2；再通过转换等式(15)得到多角度反演公式，计算施加的应力σ11、σ22和角度θ：
Figure BDA0003657440860000059
通过转换等式(15)对多角度频散数据的计算需求转换为多个频段频散的计算需求，则可以得到轴向力大小的多频反演公式：
Figure BDA0003657440860000061
其中，b1、b2均为系数，且
Figure BDA0003657440860000062
与现有技术相比，本发明产生的有益效果为：本发明采用两个声弹性系数来建立频散信息与应力检测结果之间的关系，利用相速度变化量与应力之间的正弦关系，通过两次系数拟合得到双轴力的大小与方向，即实现轴向应力的定量评估；本发明的反演算法简单，网络架构搭建容易；适用性强、精度高；适用于任意形状截面的应力的高精度快速定量检测。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明的流程图。
图2为本发明中的二维半解析有限元模型几何图。
图3为本发明中的模型坐标系示意图。
图4为本发明中S0模式相速度变化量结果图。
图5为本发明中A0模式相速度变化量结果图。
图6为本发明中相速度变化量与传播角正弦关系拟合结果图。
图7为本发明中多角度反演结果图。
图8为本发明中多频反演结果图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有付出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
如图1所示，本发明实施例提供了一种基于超声导波的适用于任意形状截面的应力检测方法，具体步骤如下：
步骤一：将经典的声弹性导波理论与半解析有限元(Semi-Analytical FiniteElement,SAFE)结合，使用SAFE方法求解声弹性导波的控制方程，并将声弹性导波的控制方程应用于商业软件COMSOL Multiphysics中，分别获得单轴受力模型、双轴受力模型和自由状态模型在多个方向的频散曲线。如图2所示构建二维几何模型来计算轴向受力的模型的声场。模型的宽为0.64mm，厚度为3.2mm，在二维模型两侧耦合周期边界条件(PeriodicBoundary Condition,PBC)可以消除边界对模型的限制，模拟导波在无限宽板中的传播，这里以无限宽铝板为例：
在经典的声弹性导波理论中，研究对象存在三种状态，包括：未产生变形时的自然状态、在预应力作用下产生了变形的初始状态、存在弹性波传播的变形状态(即最终状态)；同一个质点在三种状态下分别通过位置向量“ξ”、“X”、和“x”表示。三种状态下的所有物理变量和材料属性分别用上标“0”、“i”和“f”表示。从自然状态到初始状态、从自然状态到变形状态以及从初始状态到变形状态的变形分别表示为：
ui＝X-ξ
uf＝x-ξ (1)；
u＝x-X＝uf-ui
其中，ui表示从自然状态到初始状态的变形，X表示质点初始状态下的位置向量，ξ表示质点自然状态下的位置向量，uf表示从自然状态到变形状态的变形，x表示质点变形状态下的位置向量，u表示从初始状态到变形状态的变形。
初始状态的Lagrangian应变张量为：
Figure BDA0003657440860000071
其中，
Figure BDA0003657440860000072
为初始状态下的Lagrangian应变张量，
Figure BDA0003657440860000073
均表示从自然状态到初始状态的变形；ξi'、ξj均表示质点自然状态下的位置向量；i'、j、k、l、m、n、o、p均表示索引。
当研究对象的材料为超弹性时，Green-Lagrange应变张量通过本构方程与第二Piola-Kirchoff应力张量相关，只保留二阶和三阶弹性常数，则：
Figure BDA0003657440860000074
其中，
Figure BDA0003657440860000075
为应力张量；
Figure BDA0003657440860000076
均表示Lagrangian应变张量，ci'jkl表示二阶弹性常数，ci'jklmn表示三阶弹性常数；如果初始应变很小，进一步简化方程：
Figure BDA0003657440860000077
其中，cljop表示二阶弹性常数，
Figure BDA0003657440860000078
为柯西应变张量：
Figure BDA0003657440860000079
均表示从自然状态到初始状态的变形，ξp、ξo均表示质点自然状态下的位置向量。
如图3所示，为了获得在双轴力作用下沿任意角度的导波传播，分别以双轴力方向xi”和导波传播方向
Figure BDA00036574408600000710
为主轴方向建立的两个坐标系，在自然状态下，在
Figure BDA00036574408600000711
坐标系中的初始应力张量可表示为：
Figure BDA0003657440860000081
其中，T为初始应力张量，σ11表示沿
Figure BDA0003657440860000082
方向的主轴力，σ22表示沿
Figure BDA0003657440860000083
方向的主轴力。
通过导波传播角度
Figure BDA0003657440860000084
的旋转变换，在xi”坐标系中应力张量表示为：
Figure BDA0003657440860000085
其中，
Figure BDA0003657440860000086
均表示xi”轴和
Figure BDA0003657440860000087
轴之间夹角的余弦值，Tmn表示
Figure BDA0003657440860000088
坐标系中特定方向的初始应力张量。
声弹性波相对于初始状态的增量位移方程为：
Figure BDA0003657440860000089
其中，
Figure BDA00036574408600000810
ρi是初始状态的密度，ui'、uk均表示从初始状态到变形状态的变形，Xl、Xj、Xp均表示质点初始状态下的位置向量，ci'jkl、ci'jpl、cpjkl均表示二阶弹性常数，ci'jklmn表示表示三阶弹性常数，t表示时间；将式(2)代入等式(5)，可以得到：
Figure BDA00036574408600000811
其中
Figure BDA00036574408600000812
cljno、ci'mkl、cmjkl、ci'jml、ci'jkm均表示二阶弹性常数，ci'jklmn表示三阶弹性常数，
Figure BDA00036574408600000813
Figure BDA00036574408600000814
均表示初始状态下的柯西应变张量，δi'k表示克罗内克函数，
Figure BDA00036574408600000815
表示从自然状态到初始状态的变形，Xm表示质点在初始状态下的位置向量。
在边界上的应力张量τi'可以导出为：
Figure BDA00036574408600000816
其中，nj是在初始状态下垂直于边界的向外单位向量；从式(6)和(7)可以看出，预应力对位移增量的影响用声弹性控制方程中的系数Γi'jkl来表示。
根据半解析有限元算法的特征，假设导波在X3方向上以简谐波的形式传播，波导中的位移可以写为：
Figure BDA0003657440860000091
其中，Ui'为截面位移；下标i'＝1,2,3；I为虚数单位；ω表示角频率；k'表示X3方向上的波数；位移的导数可以写为：
Figure BDA0003657440860000092
将式(9)代入式(6)并消除相位项，声弹性导波的控制方程可以写为：
Figure BDA0003657440860000093
其中，
Figure BDA0003657440860000094
Uk表示截面位移，Γi'3kl、Γi'jk3、Γi'3k3均表示声弹性系数。
为了将该方法应用于商用有限元软件包(COMSOL Multiphysics)，需要对声弹性导波的控制方程进行改写，使其表达形式与软件中特征值问题求解的输入公式的一般表达式相同。为此引入一个新变量Pk＝k'Uk将式(10)转化为一般线性形式，则在波导域内(inΩ)：
Figure BDA0003657440860000095
ρiω2δi'kPk-k'ρiω2δi'kUk＝0 in Ω (12)；
其中，Ω表示波导域。
在边界
Figure BDA0003657440860000096
上，τi'的振型Si'可以写为：
Figure BDA0003657440860000097
在式(10)-(13)中隐藏了对重复索引k＝1,2,3和j,l＝1,2的求和；
通过找到导波在选定角频率ω下的波数k'和模式形状来求解特征值问题；通过选择实波数模式，可以识别出每个解中的传播模式，并且可以通过
Figure BDA0003657440860000098
计算相速度；最后，通过在特定频率上求解特征值问题，并在每个频率步长上组合具有最相似模式形状的模式来获得频散曲线。根据上述方法计算单轴、双轴受力模型(σ11和σ22的范围均为0到100MPa，步长为2.5MPa)与自由状态的模型在多个方向(
Figure BDA0003657440860000099
的范围为0°到90°，步长为2.5°)的频散曲线。
步骤二：将单轴受力模型的频散曲线与自由状态模型的频散曲线做差获得单轴的相速度变化量；将双轴受力模型的频散曲线与自由状态模型的频散曲线做差获得双轴的相速度变化量；图4和图5分别为3.2mm厚铝板在不同单轴力下，S0模式和A0模式在导波传播角度
Figure BDA0003657440860000109
分别为0°、30°、60°和90°时的相速度变化量随频率的变化。在任何选定的频率下，两种模式的相速度变化都随应力线性变化。在相同应力下，在1MHz以下频率范围内S0模式的相速度变化比A0模式的相速度更明显，两种模式的相速度变化均在高频范围内趋于平坦。显著而稳定的相速度变化往往意味着对应力更高的敏感性，以及更好的反演结果。因此相速度变化量的绝对值大小，是应力反演中频率段选择的重要参考因素。
步骤三：根据单轴的相速度变化量计算单轴受力模型在多个传播角度的相速度变化的计算值，将相速度变化的计算值与应力真实值作为样本以通过最小二乘法拟合出系数K1、K2；
将声弹性引起的相速度变化与导波传播角度之间的正弦依赖性扩展为适用于任意形状的截面；双轴应力引起的相速度变化与导波传播角度
Figure BDA0003657440860000101
之间的正弦关系表示为：
Figure BDA0003657440860000102
其中，△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。
步骤四：将获得的系数K1、K2与双轴的相速度变化量输入最小二乘法，进行多角度反演与多频反演，输出双轴应力的大小与方向。
将坐标系xi”顺时针旋转角度θ获得坐标系xi”'作为测试坐标系；因此，
Figure BDA0003657440860000103
方向与x1'方向之间夹角
Figure BDA0003657440860000104
满足
Figure BDA0003657440860000105
在实践中双轴力的方向，即θ通常是未知的，而
Figure BDA0003657440860000106
随
Figure BDA0003657440860000107
变化。在坐标系xi”'中，相速度△cp的变化为：
Figure BDA0003657440860000108
其中，系数a0、a1和a2分别表示为：
Figure BDA0003657440860000111
根据双轴的相速度变化量以及常数K1和K2，通过最小二乘法确定式(15)中的系数a0、a1和a2；再通过转换等式(15)得到多角度反演公式，计算施加的应力σ11、σ22和角度θ：
Figure BDA0003657440860000112
通过转换等式(15)对多角度频散数据的计算需求转换为多个频段频散的计算需求，则可以得到轴向力大小的多频反演公式：
Figure BDA0003657440860000113
其中，b1、b2均为系数，且
Figure BDA0003657440860000114
输出的反演结果与真实值之间的误差用均方误差(Mean Square Error，MSE)表示。
图7为本发明中多角度反演结果图。反演时，角度
Figure BDA0003657440860000115
的范围为10到80度(步长为10度)。由于双轴力σ11和σ22是等价的，所以图7(a)中只显示了σ11的估计误差，图7(b)为θ的误差估计。在1MHz以下，虽然只选择了8个角度的数据作为已知量进行反演，三种模式对σ11的反演误差都很小，而S0模式的低频段和A0模式的高频段之间的误差是最低的。各模态对角度θ的反演精度高且稳定。
图8为本发明中多频反演结果图。图8(a)、(b)、(c)依次为S0、A0、A1模式的多频反演结果，三个模式各自选择的频率区间可以参考多角度反演的结果：S0模式为100-500kHz，A0模式为1500-1900kHz，A1模式为800-1200kHz，步长均为80kHz。当角度
Figure BDA0003657440860000116
取一定值
Figure BDA0003657440860000117
Figure BDA0003657440860000118
时，三种模式的σ11和σ22的误差分别达到最小值：对于S0，
Figure BDA0003657440860000119
为22.5度时σ11误差为0.3485％，
Figure BDA00036574408600001110
为67.5度时σ22的误差为0.3561％；对于A0，
Figure BDA00036574408600001111
为32.5度时σ11的误差为0.5046％，
Figure BDA00036574408600001112
为57.5度时σ22的误差为0.5219％；对于A1，
Figure BDA00036574408600001113
为27.5度时σ11的误差为2.475％，
Figure BDA00036574408600001114
为62.5度时σ22的误差为2.737％。从图8可以看出，无论是哪个模式，σ11和σ22的反演精度基本对称相同，且
Figure BDA0003657440860000121
和
Figure BDA0003657440860000122
互余。在σ11的反演精度达到最高时，σ22的反演精度并不理想，说明在实际应用中，需要根据所选模式分别采集两个特定角度
Figure BDA0003657440860000123
的数据。实际应用中，多频反演法结合扫频法，可以实现利用两个方向的测量来确定双轴力的大小，这也是一种在实际应用中很有前途的应力检测方法。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种基于超声导波的适用于任意形状截面的应力检测方法，其特征在于，其步骤如下：
步骤一：使用声弹性导波理论与半解析有限元相结合的方法求解声弹性导波的控制方程，并将声弹性导波的控制方程应用于商业软件COMSOL Multiphysics中，分别获得单轴受力模型、双轴受力模型和自由状态模型的频散曲线；
步骤二：将单轴受力模型的频散曲线与自由状态模型的频散曲线做差获得单轴的相速度变化量；将双轴受力模型的频散曲线与自由状态模型的频散曲线做差获得双轴的相速度变化量；
步骤三：根据单轴的相速度变化量计算单轴受力模型在多个传播角度的相速度变化的计算值，将相速度变化的计算值与应力真实值作为样本以通过最小二乘法拟合出系数K1、K2；
步骤四：将获得的系数K1、K2与双轴的相速度变化量输入最小二乘法，进行多角度反演与多频反演，输出双轴应力的大小与方向。
2.根据权利要求1所述的基于超声导波的适用于任意形状截面的应力检测方法，其特征在于，在步骤一中：
在声弹性导波理论中，研究对象存在三种状态，包括：未产生变形时的自然状态、在预应力作用下产生了变形的初始状态、存在弹性波传播的变形状态；从自然状态到初始状态、从自然状态到变形状态以及从初始状态到变形状态的变形(ui,uf,u)分别表示为：
Figure FDA0003657440850000011
其中，ui表示从自然状态到初始状态的变形，X表示质点初始状态下的位置向量，ξ表示质点自然状态下的位置向量，uf表示从自然状态到变形状态的变形，x表示质点变形状态下的位置向量，u表示从初始状态到变形状态的变形；
初始状态的Lagrangian应变张量为：
Figure FDA0003657440850000012
其中，
Figure FDA0003657440850000013
为初始状态下的Lagrangian应变张量，
Figure FDA0003657440850000014
均表示从自然状态到初始状态的变形；ξi'、ξj均表示质点自然状态下的位置向量；i'、j、k、l、m、n、o、p均表示索引；
当研究对象的材料为超弹性时，Green-Lagrange应变张量通过本构方程与第二Piola-Kirchoff应力张量相关，只保留二阶和三阶弹性常数，则：
Figure FDA0003657440850000015
其中，
Figure FDA0003657440850000016
为应力张量；
Figure FDA0003657440850000017
均表示Lagrangian应变张量，ci'jkl表示二阶弹性常数，ci'jklmn表示三阶弹性常数；如果初始应变很小，进一步简化方程：
Figure FDA0003657440850000021
其中，cljop表示二阶弹性常数，
Figure FDA0003657440850000022
为初始状态下的柯西应变张量：
Figure FDA0003657440850000023
Figure FDA0003657440850000024
Figure FDA0003657440850000025
均表示从自然状态到初始状态的变形，ξp、ξo均表示质点自然状态下的位置向量；分别以双轴力方向xi”和导波传播方向
Figure FDA0003657440850000026
为主轴方向建立的两个坐标系，在自然状态下，在
Figure FDA0003657440850000027
坐标系中的初始应力张量可表示为：
Figure FDA0003657440850000028
其中，T为初始应力张量，σ11表示沿
Figure FDA0003657440850000029
方向的主轴力，σ22表示沿
Figure FDA00036574408500000210
方向的主轴力；
通过导波传播角度
Figure FDA00036574408500000211
的旋转变换，在xi”坐标系中应力张量表示为：
Figure FDA00036574408500000212
其中，
Figure FDA00036574408500000213
均表示xi”轴和
Figure FDA00036574408500000214
轴之间夹角的余弦值，Tmn表示
Figure FDA00036574408500000215
坐标系中特定方向的初始应力张量；
声弹性波相对于初始状态的增量位移方程为：
Figure FDA00036574408500000216
其中，
Figure FDA00036574408500000217
ρi是初始状态的密度，ui'、uk均表示从初始状态到变形状态的变形，Xl、Xj、Xp均表示质点初始状态下的位置向量，ci'jkl、ci'jpl、cpjkl均表示二阶弹性常数，ci'jklmn表示表示三阶弹性常数，t表示时间；
将式(2)代入等式(5)，可以得到：
Figure FDA00036574408500000218
其中
Figure FDA00036574408500000219
cljno、ci'mkl、cmjkl、ci'jml、ci'jkm均表示二阶弹性常数，ci'jklmn表示三阶弹性常数，
Figure FDA00036574408500000220
Figure FDA00036574408500000221
均表示初始状态下的柯西应变张量，δi'k表示克罗内克函数，
Figure FDA00036574408500000222
表示从自然状态到初始状态的变形，Xm表示质点在初始状态下的位置向量；
在边界上的应力张量τi'可以导出为：
Figure FDA0003657440850000031
其中，nj是在初始状态下垂直于边界的向外单位向量；从式(6)和(7)可以看出，预应力对位移增量的影响用声弹性控制方程中的系数Γi'jkl来表示；
根据半解析有限元算法的特征，假设导波在X3方向上以简谐波的形式传播，波导中的位移可以写为：
Figure FDA0003657440850000032
其中，Ui'为截面位移；下标i'＝1,2,3；I为虚数单位；ω表示角频率；k'表示X3方向上的波数；位移的导数可以写为：
Figure FDA0003657440850000033
将式(9)代入式(6)并消除相位项，声弹性导波的控制方程可以写为：
Figure FDA0003657440850000034
其中，
Figure FDA0003657440850000035
Uk表示截面位移，Γi'3kl、Γi'jk3、Γi'3k3均表示声弹性系数；
引入一个新变量Pk＝k'Uk将式(10)转化为一般线性形式，则在波导域内(inΩ)：
Figure FDA0003657440850000036
ρiω2δi'kPk-k'ρiω2δi'kUk＝0 in Ω (12)；
其中，Ω表示波导域；
在波导边界
Figure FDA0003657440850000037
上，τi'的振型Si'可以写为：
Figure FDA0003657440850000038
在式(10)-(13)中隐藏了对重复索引k＝1,2,3和j,l＝1,2的求和；
通过找到导波在选定角频率ω下的波数k'和模式形状来求解特征值问题；通过选择实波数模式，可以识别出每个解中的传播模式，并且可以通过
Figure FDA0003657440850000041
计算相速度；最后，通过在特定频率上求解特征值问题，并在每个频率步长上组合具有最相似模式形状的模式来获得频散曲线。
3.根据权利要求2所述的基于超声导波的适用于任意形状截面的应力检测方法，其特征在于，在步骤三中：
将声弹性引起的相速度变化与导波传播角度之间的正弦依赖性扩展为适用于任意形状的截面；双轴应力引起的相速度变化与导波传播角度
Figure FDA0003657440850000042
之间的正弦关系表示为：
Figure FDA0003657440850000043
其中，△cp(·)表示相速度；由式(14)，根据单轴受力模型在多个传播角度的相速度变化，通过最小二乘法拟合出不同频率下对应的系数K1、K2。
4.根据权利要求3所述的基于超声导波的适用于任意形状截面的应力检测方法，其特征在于，在步骤四中：
将坐标系xi”顺时针旋转角度θ获得坐标系xi”'作为测试坐标系；因此，
Figure FDA0003657440850000044
方向与x1'方向之间夹角
Figure FDA0003657440850000045
满足
Figure FDA0003657440850000046
在坐标系xi”'中，相速度△cp的变化为：
Figure FDA0003657440850000047
其中，系数a0、a1和a2分别表示为：
Figure FDA0003657440850000048
根据双轴的相速度变化量以及常数K1和K2，通过最小二乘法确定式(15)中的系数a0、a1和a2；再通过转换等式(15)得到多角度反演公式，计算施加的应力σ11、σ22和角度θ：
Figure FDA0003657440850000049
通过转换等式(15)对多角度频散数据的计算需求转换为多个频段频散的计算需求，则可以得到轴向力大小的多频反演公式：
Figure FDA0003657440850000051
其中，b1、b2均为系数，且
Figure FDA0003657440850000052
CN202210564728.1A 2022-05-23 2022-05-23 一种基于超声导波的适用于任意形状截面的应力检测方法 Active CN114739546B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210564728.1A CN114739546B (zh) 2022-05-23 2022-05-23 一种基于超声导波的适用于任意形状截面的应力检测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210564728.1A CN114739546B (zh) 2022-05-23 2022-05-23 一种基于超声导波的适用于任意形状截面的应力检测方法

Publications (2)

Publication Number Publication Date
CN114739546A true CN114739546A (zh) 2022-07-12
CN114739546B CN114739546B (zh) 2022-12-02

Family

ID=82287865

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210564728.1A Active CN114739546B (zh) 2022-05-23 2022-05-23 一种基于超声导波的适用于任意形状截面的应力检测方法

Country Status (1)

Country Link
CN (1) CN114739546B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116305665A (zh) * 2023-05-15 2023-06-23 广东工业大学 一种工件横截面残余应力分布的分析方法及相关装置

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101441199A (zh) * 2008-12-26 2009-05-27 北京工业大学 钢绞线预应力的高阶纵向导波测量方法
CN106768529A (zh) * 2017-01-24 2017-05-31 清华大学 具有预应力的薄壁软材料或软组织材料力学特性分析方法
CN111678630A (zh) * 2020-06-18 2020-09-18 哈尔滨工业大学（深圳） 基于超声导波应力敏感度分析的钢绞线单轴应力检测方法
WO2020233359A1 (zh) * 2019-05-20 2020-11-26 北京工业大学 一种用于金属薄板中应力分布测量的非线性Lamb波混频方法
CN112903157A (zh) * 2021-01-19 2021-06-04 吉林大学 基于纵向模态超声导波的圆管型结构的应力监测方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101441199A (zh) * 2008-12-26 2009-05-27 北京工业大学 钢绞线预应力的高阶纵向导波测量方法
CN106768529A (zh) * 2017-01-24 2017-05-31 清华大学 具有预应力的薄壁软材料或软组织材料力学特性分析方法
WO2020233359A1 (zh) * 2019-05-20 2020-11-26 北京工业大学 一种用于金属薄板中应力分布测量的非线性Lamb波混频方法
CN111678630A (zh) * 2020-06-18 2020-09-18 哈尔滨工业大学（深圳） 基于超声导波应力敏感度分析的钢绞线单轴应力检测方法
CN112903157A (zh) * 2021-01-19 2021-06-04 吉林大学 基于纵向模态超声导波的圆管型结构的应力监测方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
SARAVANAN, TJ: "Investigation on elastic guided wave propagation in a parallel axisymmetric stressed viscoelastic cylindrical waveguide", 《MECHANICS OF ADVANCED MATERIALS AND STRUCTURES》 *
刘飞: "波导结构的特征频率法及其超声导波声弹性效应研究", 《中国博士学位论文全文数据库工程科技Ⅰ辑》 *
杨正岩: "异型截面结构超声导波监测的半解析有限元分析方法研究", 《中国博士学位论文全文数据库基础科学辑》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116305665A (zh) * 2023-05-15 2023-06-23 广东工业大学 一种工件横截面残余应力分布的分析方法及相关装置
CN116305665B (zh) * 2023-05-15 2023-08-11 广东工业大学 一种工件横截面残余应力分布的分析方法及相关装置

Also Published As

Publication number Publication date
CN114739546B (zh) 2022-12-02

Similar Documents

Publication Publication Date Title
Hein et al. Computationally efficient delamination detection in composite beams using Haar wavelets
Gravenkamp et al. The computation of dispersion relations for three-dimensional elastic waveguides using the scaled boundary finite element method
Pawar et al. Damage detection in beams using spatial Fourier analysis and neural networks
Yan et al. Structural Health Monitoring Using High‐Frequency Electromechanical Impedance Signatures
Hosseini et al. Numerical simulation of Lamb wave propagation in metallic foam sandwich structures: a parametric study
Fouaidi et al. Nonlinear bending analysis of functionally graded porous beams using the multiquadric radial basis functions and a Taylor series-based continuation procedure
Barouni et al. A layerwise semi-analytical method for modeling guided wave propagation in laminated and sandwich composite strips with induced surface excitation
Du et al. Electromechanical impedance temperature compensation and bolt loosening monitoring based on modified Unet and multitask learning
CN114739546B (zh) 一种基于超声导波的适用于任意形状截面的应力检测方法
Yu et al. A new procedure for exploring the dispersion characteristics of longitudinal guided waves in a multi-layered tube with a weak interface
Li et al. Modelling of the high-frequency fundamental symmetric Lamb wave using a new boundary element formulation
Teter et al. On using load-axial shortening plots to determine the approximate buckling load of short, real angle columns under compression
Yang et al. Transient wave propagation of isotropic plates using a higher-order plate theory
Kapuria et al. AC 1‐continuous time domain spectral finite element for wave propagation analysis of Euler–Bernoulli beams
Yuan et al. Theoretical and experimental study on interface stiffness measurement of rough surface using improved acoustic model
Shi et al. An online stress monitoring strategy based on Wigner–Ville time–frequency energy extraction of single-frequency dual mode Lamb waves
Nag et al. Identification of delamination in a composite beam using a damaged spectral element
dos Santos Jr et al. Comparison of acoustoelastic methods to evaluate stresses in steel plates and bars
Hernández et al. Numerical solution of a wave propagation problem along plate structures based on the isogeometric approach
Gao et al. Stable and accurate computation of dispersion relations for layered waveguides, semi-infinite spaces and infinite spaces
Jain et al. C1-continuous time-domain spectral finite element for modeling guided wave propagation in laminated composite strips based on third-order theory
Zhao et al. Stress inversion in waveguides with arbitrary cross sections with acoustoelastic guided waves
Lopes et al. Longitudinal wave scattering in thin plates with symmetric damage considering oblique incidence
Da et al. A rapid and accurate technique with updating strategy for surface defect inspection of pipelines
Lomazov et al. Mathematical modeling of diagnostics of residual stresses in a layered medium

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant