CN103926315B - 一种基于单纯形法的各向同性薄板材料弹性性质获取方法 - Google Patents

Abstract

一种基于单纯形法的各向同性薄板材料弹性性质获取方法，属于超声导波无损检测与评估领域；基于声学显微镜技术，自行开发的材料弹性系数的超声测量系统，采用线聚焦PVDF探头，通过纵波和表面波波速的同时测量，可实现材料的弹性系数无损检测。通过一种新的反演算法，利用该线聚焦超声探头来获取薄板材料的弹性性质。该方法基于单纯形法诱导目标函数于频散特征方程的系数矩阵行列式中，弹性性质与式样密度都可以通过声学性质及所测密度来获得。本方法可解决各向同性薄板材料波速提取的问题，是一种先进的材料波速提取方法。

Description

一种基于单纯形法的各向同性薄板材料弹性性质获取方法

技术领域

本发明涉及一种基于单纯形法获取涂层材料弹性系数方法，属于超声导波无损检测与评估领域。

背景技术

随着材料科学的发展，各种新型材料不断涌现，许多材料表现出超常的力学性能，因此，材料的力学性能测试是材料科学中的一个重要组成部分。研究表明，材料的弹性系数与材料制备工艺密切相关。利用超声波对材料弹性性质进行测量是无损检测领域很有前景的测量方法之一。在各向同性均质薄板材料中，超声波的传播具有频散特性，而该特性又包含了大量材料机械性质的信息，因此通过波速与波长或频率的关系——即频散曲线，即可反演出薄板材料的机械性质，如厚度、密度、纵波波速、横波波速等。

基于声学显微镜技术，自行开发的材料弹性系数的超声测量系统，采用线聚焦PVDF探头，通过纵波和表面波波速的同时测量，可实现材料的弹性系数无损检测。通过一种新的反演算法，利用该线聚焦超声探头来获取薄板材料的弹性性质。该方法基于单纯形法诱导目标函数于频散特征方程的系数矩阵行列式中，弹性性质（杨氏模量E，剪切模量G）与式样密度都可以通过声学性质（纵波波速，横波波速）及所测密度来获得。

发明内容

本发明的目的是为了解决各向同性薄板材料波速提取的问题，提出一种先进的材料波速提取方法。

步骤1：确立仿真目标函数

通过实验获得的频率f与表面波波速c与改变C_L,C_T,ρ,h后的残差值进行叠加，经过叠加后的目标函数在某一组C_L,C_T,ρ,h值时最小，此时可反演出其横、纵波波速，镀层密度、镀层厚度。

${\mathrm{\Π}}_s=\underset{j=1}{\overset{N_s}{\mathrm{\Σ}}}\left[K(f_j^s,c_j^s;C_L^g,C_T^g,{\mathrm{\ρ}}^g,h^g)\right]$

步骤2：搭建测试系统。

为了方便散焦步进测量，搭建了一套进行散焦步进测量的测试系统，如图1所示。该测试系统主要包括：试样1、水槽与水2、换能器3、移动平台4、脉冲激励/接收仪5、示波器6、GPIB总线7、PXI总控制系统8、移动伺服马达9、旋转轴10。其中，在移动平台4下面安装换能器3，换能器3与脉冲激励/接收仪5相连，脉冲激励/接收仪5与示波器6相连，示波器6通过GPIB总线7与PXI总控制系统8相连，PXI总控制系统8与移动伺服马达9相连，同时PXI总控制系统8与旋转轴10相连。

步骤3：聚焦面数据采集。

将被测试样置于换能器的聚焦面，脉冲激励/接收仪5在发出一个带宽为10-200MHz的脉冲后转换为接收状态，当接收到反射信号后，将信号传输进示波器6，示波器的采样频率为f_S，f_S为0.5-5GHz，采样点数为N_s，N_s的取值范围为10000-100000点。经过示波器的低通滤波后，通过GPIB总线7存储进PXI总控制系统8。

步骤4：散焦测量。

将换能器垂直向下移动一个距离Δz₀，Δz₀的取值范围为1-50μm，待移动完成后进行数据采集，采样频率为f_S，采样点数为N_s。采集结束后再将换能器垂直向下移动Δz₀进行数据采集，如此循环往复，共移动距离z，z的取值范围为2-20mm，因此将得到M组电压数据，M由z与Δz₀共同决定，为40-20000组。

步骤5：时域傅里叶变换。

将所有数据沿散焦距离排列好，对测得的数据进行时域傅里叶变换：

$A_i\left[k\right]=\underset{n=0}{\overset{N_s-1}{\mathrm{\Σ}}}{x}_i\left[n\right]e^{-j2\mathrm{\πnk}/N_s}$

其中：A_i为时域傅里叶变换后的频谱值，x_i代表一组电压数据，i＝0,1,2…M-1，k＝0,1,2…N_s-1，j代表虚部。

步骤6：空间傅里叶变换。

为了得到精确的振荡周期Δz，需要对时域傅里叶变换的结果再进行沿散焦距离方向的空间傅里叶变换，将散焦距离z变换至z^-1域：

$B_i\left[k\right]=\underset{m=0}{\overset{M-1}{\mathrm{\Σ}}}{A}_m\left[k\right]e^{-j2\mathrm{\πmi}/M}$

其中：B_i为空间傅里叶变换后的频谱值，A_m代表沿散焦方向的时域傅里叶变换的频谱值，i＝0,1,2…M-1，k＝0,1,2…N_s-1，j代表虚部。沿z^-1域的曲线峰值即为振荡周期Δz的倒数。

步骤7：多模态追踪。

对多模态中每个模态的极大值分别进行追踪，即可得到连续的z-1值，其倒数即为Δz。

步骤8：波速提取。

V(z)曲线理论，可根据如下公式进行波速的计算：

$v_{\mathrm{Lamb}}=v_w\·{[1-{(1-\frac{v_w}{2\·f\·\mathrm{\Δz}})}^2]}^{-1/2}$

将水的超声波波速v_W、每一个极大值对应的频率f与Δz代入其中：Δz为V(z)曲线振荡周期，v_w为水中的超声波波速，f为换能器的激励频率，v_Lamb为材料的兰姆波波速。测量被测材料的V(z)曲线振荡周期是波速提取的关键。通过实验数据获得其频散波速曲线。

步骤9：单纯形法确定波速值

将试验测得数据带入目标函数中，通过单纯形法改变C_L,C_T,ρ,h四个变量值，在其保证目标函数值最小时，即确定各值。

附图说明

图1散焦测量系统示意图。

图2表面波传播示意图。

图3Lamb波在各向同性薄板中传播示意图。

图4聚焦面时域波形图。

图57.5MHz频率下V(z)振荡曲线图。

图67.5MHz频率下z^-1域曲线图。

图中：1、试样，2、水槽与水，3、换能器，4、移动平台，5、脉冲激励/接收仪，6、示波器，7、GPIB总线，8、PXI总控制系统，9、移动伺服马达，10、旋转轴。

具体实施方式

以下结合具体实例对本发明的内容做进一步的详细说明：

步骤1确立仿真目标函数

通过对各项同性材料试件内部部分波线性组合，将其定义在一个关于边界条件的特征方程中，从而得到频散关系。一块厚度为2h的板将其置于笛卡尔坐标系下，如图1所示，x1-x3平面是Lamb波传播的矢面，而坐标x1代表波的传播方向。全局坐标（x1,x2,x3）为下面的理论服务。

步骤1.1根据运动方程：（可忽略重力ρf_j影响）；及本构方程：σ_ij＝c_ijklε_kl；几何方程：推知控制波的传播方程 $c_{\mathrm{ijkl}}\·u_{k,\mathrm{li}}=\mathrm{\ρ}{\stackrel{\·\·}{u}}_j.$

$\left.\begin{array}{c}{\mathrm{\σ}}_{\mathrm{ij},i}={\mathrm{\ρ}\stackrel{\·\·}{u}}_j\\ {\mathrm{\σ}}_{\mathrm{ij}}=c_{\mathrm{ijkl}}{\mathrm{\ϵ}}_{\mathrm{kl}}\\ {\mathrm{\ϵ}}_{\mathrm{ij}}=\frac{1}{2}(u_{i,j}+u_{j,i})\end{array}\right\}\⇒\underset{\‾}{{\mathrm{\σ}}_{\mathrm{ij}}}=c_{\mathrm{ijkl}}{\mathrm{\ϵ}}_{\mathrm{kl}}=c_{\mathrm{ijkl}}\frac{1}{2}(u_{k,l}+u_{l,k})=\underset{\‾}{c_{\mathrm{ijkl}}\·u_{k,l}}\⇒c_{\mathrm{ijkl}}\·u_{k,\mathrm{li}}={\mathrm{\ρ}\stackrel{\·\·}{u}}_{j.}$

其中σ_ij为应力张量，ρ是密度，u_j是质点位移向量，t为时间，圆点代表时间的不同，指针从逗号后开始代表空间坐标的不同，重复指针的合计就如张量标志自动的被设定在自此之后。c_ijkl是弹性系数在之后我们可以将其收缩在C_IJ系数中，ε_kl是应变张量，其与质点位移有关。

根据弹性动力学理论，现假设一个平面谐波沿着x₁传播，其角频率为ω，该波相应的物理量在x₂方向上独立，波的形式u_k如下：

$(u_1,u_2,u_3)=(U,V,W)\·e^{\mathrm{j\ξ}(x_1+\mathrm{\α}{x}_3)}\·e^{-\mathrm{j\ωt}}$

其中ξ为波数，而α则表示一个未知系数，(ξ,0,αξ)为该波的传播向量，(U,V,W)为该平面波相应的振幅。将u_k带入控制波的传播方程 $c_{\mathrm{ijkl}}\·u_{k,\mathrm{li}}=\mathrm{\ρ}{\stackrel{\·\·}{u}}_j.$ 得到： ${\left[K_{\mathrm{ij}}\right]}_{3\×3}\·\left\{\begin{array}{c}U\\ V\\ W\end{array}\right\}=\left\{\begin{array}{c}0\\ 0\\ 0\end{array}\right\},$ 该线性代数形式[K_ij]_3×3其中：

K₁₁＝C₁₁+2·C₁₅·α+C₅₅·α²-ρ₀·c²

K₂₂＝C₆₆+2·C₄₆·α+C₄₄·α²-ρ₀·c²

K₃₃＝C₅₅+2·C₃₅·α+C₃₃·α²-ρ₀·c²

K₁₂＝K₂₁＝C₁₆+(C₁₄+C₅₆)·α+C₄₅·α²

K₁₃＝K₃₁＝C₁₅+(C₁₃+C₅₅)·α+C₃₅·α²

K₂₃＝K₃₂＝C₅₆+(C₃₆+C₄₅)·α+C₃₄·α²

为了使ω，ξ存在值，α的值要使K矩阵行列式为零，也就是说α可以作为K矩阵的特征值，向量(U,V,W)是其相对相应的特征向量。

在各向同性材料薄板的波场中可以利用以下的弹性系数关系。首先我们将C_IJ定义在一个各向同性材料中，得到的弹性系数为:

$C_{\mathrm{IJ}}=\left[\begin{array}{cccccc}\mathrm{\λ}+2\mathrm{\μ}& \mathrm{\λ}& \mathrm{\λ}& 0& 0& 0\\ & \mathrm{\λ}+2\mathrm{\μ}& \mathrm{\λ}& 0& 0& 0\\ & & \mathrm{\λ}+2\mathrm{\μ}& 0& 0& 0\\ & & & \mathrm{\μ}& 0& 0\\ & \mathrm{Sym}.& & & \mathrm{\μ}& 0\\ & & & & & \mathrm{\μ}\end{array}\right]$

λ与μ为拉梅系数，我们定义C_IJ＝c_ijkl，其中ij→I或J，定义11→1，22→2，33→3，23或32→4，31或13→5，12或21→6。

获得α的表示式，从而可知该波的传播方式。通常，在各向同性材料中，α的值有四个。

$\left\{\begin{array}{c}{\mathrm{\α}}_1=-{\mathrm{\α}}_3=\sqrt{{\left(\frac{c}{C_L}\right)}^2-1}\\ {\mathrm{\α}}_2=-{\mathrm{\α}}_4=\sqrt{{\left(\frac{c}{C_T}\right)}^2-1}\end{array}\right.$

将不同的α值下对应的波的形式累加，得到以下位移、应力公式来表示波场。A_q为该平面波的振幅系数

$(u_1,u_2,u_3)=\underset{q=1}{\overset{4}{\mathrm{\Σ}}}(U_q,0{,W}_q)\·A_q\·e^{\mathrm{i\ξ}(x_1+{\mathrm{\α}}_q{x}_3)}\·e^{-\mathrm{i\ωt}}$

$({\mathrm{\σ}}_{33},{\mathrm{\σ}}_{13},{\mathrm{\σ}}_{23})=\underset{q=1}{\overset{4}{\mathrm{\Σ}}}(D_{1q},D_{2q},D_{3q})\·A_q\·e^{\mathrm{i\ξ}(x_1+{\mathrm{\α}}_q{x}_3)}\·e^{-\mathrm{i\ωt}}$

获得边界条件的系数矩阵：

在无外力为边界条件的自由载荷状态下，板厚为2h，得到量化的漏表面波频散曲线，其系数矩阵为[M]_4×4：

$\left\{\begin{array}{c}{\mathrm{\σ}}_{33}{|}_{x_3=+h}\\ {\mathrm{\σ}}_{13}{|}_{x_3=+h}\\ {\mathrm{\σ}}_{33}{|}_{x_3=-h}\\ {\mathrm{\σ}}_{13}{|}_{x_3=-h}\end{array}\right\}={\left[M\right]}_{4\×4}\left\{\begin{array}{c}{A}_1\\ A_2\\ A_3\\ A_4\end{array}\right\}=\left\{\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right\}$

在各向同性的平面板中相应的位移幅值比为以下值：

$\left\{\begin{array}{c}{U}_q=1\\ W_1={\mathrm{\α}}_1,W_2=-\frac{1}{{\mathrm{\α}}_2},W_3={\mathrm{\α}}_3,W_4=-\frac{1}{{\mathrm{\α}}_4}\end{array}\right.$

α_q已经得出为:

$\left\{\begin{array}{c}{\mathrm{\α}}_1=-{\mathrm{\α}}_3=\sqrt{{\left(\frac{c}{C_L}\right)}^2-1}\\ {\mathrm{\α}}_2=-{\mathrm{\α}}_4=\sqrt{{\left(\frac{c}{C_T}\right)}^2-1}\end{array}\right.$

在各向同性的平面板中相应的应变幅值比为以下值:

$\left\{\begin{array}{c}{D}_{1q}=\mathrm{\λ}+(\mathrm{\λ}+2\mathrm{\μ}{W}_q){\mathrm{\α}}_q\\ D_{2q}=\mathrm{\μ}(W_q+{\mathrm{\α}}_q)\end{array}\right.$

λ,μ是拉梅系数,可表示横波波速C_T，纵波波速C_L为以下形式:

$\left\{\begin{array}{c}{C}_L^2=\frac{\mathrm{\λ}+2\mathrm{\μ}}{\mathrm{\ρ}}\\ C_T^2=\frac{\mathrm{\μ}}{\mathrm{\ρ}}\end{array}\right.$

最后，可以得出系数矩阵[M]_4×4:

${\left[M\right]}_{4\×4}=\left[\begin{array}{cccc}{D}_{11}{e}^{+\mathrm{i\ξh}{\mathrm{\α}}_1}& D_{12}{e}^{+\mathrm{i\ξh}{\mathrm{\α}}_2}& D_{13}{e}^{+\mathrm{i\ξh}{\mathrm{\α}}_3}& D_{14}{e}^{+\mathrm{i\ξh}{\mathrm{\α}}_4}\\ D_{21}{e}^{+\mathrm{i\ξh}{\mathrm{\α}}_1}& D_{22}{e}^{+\mathrm{i\ξh}{\mathrm{\α}}_2}& D_{23}{e}^{+\mathrm{i\ξh}{\mathrm{\α}}_3}& D_{24}{e}^{+\mathrm{i\ξh}{\mathrm{\α}}_4}\\ D_{11}{e}^{-\mathrm{i\ξh}{\mathrm{\α}}_1}& D_{12}{e}^{-\mathrm{i\ξh}{\mathrm{\α}}_2}& D_{13}{e}^{-\mathrm{i\ξh}{\mathrm{\α}}_3}& D_{14}{e}^{-\mathrm{i\ξh}{\mathrm{\α}}_4}\\ D_{21}{e}^{-\mathrm{j\ξh}{\mathrm{\α}}_1}& D_{22}{e}^{-\mathrm{i\ξh}{\mathrm{\α}}_2}& D_{23}{e}^{-\mathrm{i\ξh}{\mathrm{\α}}_3}& D_{24}{e}^{-\mathrm{i\ξh}{\mathrm{\α}}_4}\end{array}\right]$

由边界条件[M]_4×4行列式可得四个方向的频散曲线，再得到其频散特征函数K，为了得到非零解，函数K的值应为0，C_L,C_T,ρ,h都为已知，改变f,c的值使目标函数值最小保证其尽量趋近于0，频散曲线，频散特征函数，仿真的目标函数Π_s分别如下：

$\left\{\begin{array}{c}{\mathrm{\σ}}_{33}{|}_{x_3=+h}\\ {\mathrm{\σ}}_{13}{|}_{x_3=+h}\\ {\mathrm{\σ}}_{33}{|}_{x_3=-h}\\ {\mathrm{\σ}}_{13}{|}_{x_3=-h}\end{array}\right\}={\left[M\right]}_{4\×4}\left\{\begin{array}{c}{A}_1\\ A_2\\ A_3\\ A_4\end{array}\right\}=\left\{\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right\}$

K(f,c;C_L,C_T,ρ,h)＝log₁₀(|[M]_4×4|)

${\mathrm{\Π}}_s=\underset{j=1}{\overset{N_s}{\mathrm{\Σ}}}\left[K(f_j^s,c_j^s;C_L^g,C_T^g,{\mathrm{\ρ}}^g,h^g)\right]$

步骤1.2：确立波速提取的公式。

在单频激励/接收的情况下，图2所示的漏表面波传播示意图中，上表面的直接反射回波I传播的时间与漏表面波L的传播时间分别为：

$t_1=\frac{2(R-\mathrm{\Δz})}{v_w}$

$t_2=\frac{2(R-\frac{\mathrm{\Δz}}{\cos {\mathrm{\θ}}_{\mathrm{SAW}}})}{v_w}+\frac{2\·\mathrm{\Δz}\·\tan {\mathrm{\θ}}_{\mathrm{SAW}}}{v_{\mathrm{SAW}}}$

其中R为聚焦半径，Δz为散焦距离，v_w为水的超声波波速，θ_SAW为产生表面波的瑞利角，v_SAW为材料的表面波波速。因此两者的时间差为：

$\mathrm{\Δt}=t_2-t_1=\frac{2(1-\cos {\mathrm{\θ}}_{\mathrm{SAW}})}{v_w}\·\mathrm{\Δz}$

即：

$\cos {\mathrm{\θ}}_{\mathrm{SAW}}=1-\frac{v_w\·\mathrm{\Δt}}{2\·\mathrm{\Δz}}$

将Snell定律：

$\sin {\mathrm{\θ}}_{\mathrm{SAW}}=\frac{v_w}{v_{\mathrm{SAW}}}$ 或 ${\mathrm{\θ}}_{\mathrm{SAW}}=s{\mathrm{in}}^{-1}\left(\frac{v_w}{v_{\mathrm{SAW}}}\right)$

代入后，可得：

$\frac{v_w}{v_{\mathrm{SAW}}}=\sqrt{1-{(1-\frac{v_w}{2}\·\frac{\mathrm{\Δt}}{\mathrm{\Δz}})}^2}$

此时如果Δz恰为一个V(z)曲线的振荡周期时，1/Δt则为换能器的激励频率f。如果Δz能够确定，便可使用如下公式进行表面波波速的计算：

$v_{\mathrm{SAW}}=v_w\·{[1-{(1-\frac{v_w}{2\·\mathrm{f\Δz}})}^2]}^{-1/2}$

因此，测量被测材料的V(z)曲线振荡周期成为波速提取的重点。

步骤2：搭建测试系统。

为了方便散焦步进测量，搭建了一套进行散焦步进测量的测试系统，如图1所示。该测试系统主要包括：试样1、水槽与水2、换能器3、移动平台4、脉冲激励/接收仪5、示波器6、GPIB总线7、PXI总控制系统8、移动伺服马达9、旋转轴10。其中，在移动平台4下面安装换能器3，换能器3与脉冲激励/接收仪5相连，脉冲激励/接收仪5与示波器6相连，示波器6通过GPIB总线7与PXI总控制系统8相连，PXI总控制系统8与移动伺服马达9相连，同时PXI总控制系统8与旋转轴10相连。

步骤3：聚焦面数据采集。

以长方体碳化钨为被测试样，其尺寸为40mm×40mm×10mm，将换能器3聚焦到试样的上表面，通过脉冲激励/接收仪5在发出一个带宽为10-200MHz的脉冲后转换为接收状态，当接收到反射信号后，将信号传输进示波器6，示波器的采样频率f_S＝2.5GHz，采样点数N_s＝10000。经过示波器的低通滤波后，通过GPIB总线7存储进PXI总控制系统，聚焦面的时域波形如图4所示。

步骤4：散焦测量。

将换能器朝试样方向移动Δz₀＝10μm，待移动完成后进行电压数据采集，采集结束后再将换能器朝试样方向移动Δz₀＝10μm进行数据采集，采样频率f_S＝2.5GHz，采样点数N_s＝10000，如此循环往复，共移动4mm，因此将得到400组电压数据，将聚焦面的电压数据包含在内共得到M＝401组电压数据。

步骤5：时域傅里叶变换。

将测得的数据进行时域傅里叶变换。

$A_i\left[k\right]=\underset{n=0}{\overset{N_s-1}{\mathrm{\Σ}}}{x}_i\left[n\right]e^{-j2\mathrm{\πnk}/N_s}$

其中：A_i为时域傅里叶变换后的频谱值，x_i代表一组电压数据，i＝0,1,2…M-1，k＝0,1,2…N_s-1，j代表虚部，N_s＝10000。例如，所得A_i[k]，i＝0,1,2…M-1，k＝0,1,2…N_s-1。7.5MHz频率下的振荡曲线如图5所示。

步骤6：空间傅里叶变换。

为了得到精确的振荡周期Δz，需要对时域傅里叶变换的结果再进行沿散焦距离方向的空间傅里叶变换，将散焦距离z变换至z^-1域：

$B_i\left[k\right]=\underset{m=0}{\overset{M-1}{\mathrm{\Σ}}}{A}_m\left[k\right]e^{-j2\mathrm{\πmi}/M}$

其中：B_i为空间傅里叶变换后的频谱值，A_m代表沿散焦方向的时域傅里叶变换的频谱值，i＝0,1,2…M-1，k＝0,1,2…N_s-1，M＝401，j代表虚部。所得B_i[k]，i＝0,1,2…M-1，k＝0,1,2…N_s-1。例，7.5MHz频率下z^-1域的曲线如图6所示。

步骤7：模态追踪。

对2.5-22.5MHz范围内的峰值进行追踪，即可找出该频率段连续的Δz值。

步骤8：波速提取。

将水中的超声波波速v_W＝1500m/s，每一个峰值对应的频率与Δz带入公式即可得到该频率段内连续的表面波波速。

步骤9：单纯形法获得横、纵波波速

通过单纯形法使目标函数残差绝对值最接近零，单纯形法改变C_L,C_T,ρ,h使目标函数值的达到所需范围（单纯形法简介）。

单纯形搜索法通过构造单纯形来逼近极小点，每构造一个单纯形，确定其最高点和最低点，然后通过扩展或压缩、反射构造新的单纯形，目地是使得极小点能包含于单纯形内。

本方法中未知量有四个（C_L,C_T,ρ,h），为四维变量问题。

用单纯形搜索法求无约束问题minF(x)，x∈Rⁿ的算法步骤如下：

①选取初始单纯形{x⁰,x¹,…,xⁿ}，反映系数α＞1，紧缩系数θ∈(0,1)，扩展系数γ＞1，收缩系数β∈（0，1）及精度ε>0，置k=0；

②将单纯形的n+1个顶点按目标函数值的大小重新编号，使顶点的编号满足F(x⁰)≤F(x¹)≤…≤F(x^n-1)≤F(xⁿ)；

③令 $x^{n+1}=\frac{1}{n}\underset{j=0}{\overset{n-1}{\mathrm{\Σ}}}{x}^j,$ 若 ${\left\{\frac{1}{n+1}\underset{j=0}{\overset{n}{\mathrm{\Σ}}}{[F\left(x^j\right)-F\left(x^{n+1}\right)]}^2\right\}}^2\≤\mathrm{\ϵ}$ 停止迭代输出x⁰，否则转入④；

④计算xⁿ⁺²＝xⁿ⁺¹+α(xⁿ⁺¹-x），若F(xⁿ⁺²)＜F(x⁰)，转⑤，否则当F(xⁿ⁺²)＜F(x^n-1)时转⑥，当F(xⁿ⁺²)≥F(x^n-1)转⑦；

⑤计算xⁿ⁺³＝xⁿ⁺¹+γ(xⁿ⁺²-xⁿ⁺¹)，若F(xⁿ⁺³)＜F(x⁰)，令xⁿ＝xⁿ⁺³，转②，否则转⑥；

⑥令xⁿ＝xⁿ⁺²，转②

⑦令xⁿ＝{xⁱ|F(xⁱ)＝min(F(xⁿ),F(xⁿ⁺²))}，计算xⁿ⁺⁴＝xⁿ⁺¹+β(xⁿ-xⁿ⁺¹)，若F(xⁿ⁺⁴)＜F(xⁿ)，令xⁿ＝xⁿ⁺⁴，转②，否则转⑧；

⑧令x^j＝x⁰+θ(x^j-x⁰)，j＝0,1,…,n，转②

其中，N_s为通过的仿真频散曲线得到的数据点数，这样只有那些相应残差范围内的被累加计算起来，就不必再将整幅理论与实际曲线的波速平方的残差总和最小值都计算出来。并且以前的方法，必须把所有的理论频散曲线都确定，这种方式耗时且数据量冗长，因此，单纯形的方法可以加速并化简程序。

Claims (2)

1.一种基于单纯形法的各向同性薄板材料弹性性质获取方法，其特征在于：该步骤的实现过程如下，

步骤1：确立仿真目标函数

通过对各项同性材料试件内部部分波线性组合，将其定义在一个关于边界条件的特征方程中，从而得到频散关系；设一块厚度为2h的板，将其置于笛卡尔坐标系下，x1-x3平面是Lamb波传播的矢面，而坐标x1代表波的传播方向；全局坐标(x1,x2,x3)为理论服务；

步骤1.1：根据运动方程：及本构方程：σ_ij＝c_ijklε_kl；几何方程：推知控制波的传播方程

$\left.\begin{array}{c}{\mathrm{\σ}}_{ij,i}=\mathrm{\ρ}{\stackrel{\·\·}{u}}_j\\ {\mathrm{\σ}}_{ij}=c_{ijkl}{\mathrm{\ϵ}}_{kl}\\ {\mathrm{\ϵ}}_{ij}=\frac{1}{2}(u_{i,j}+u_{j,i})\end{array}\right\}\⇒\underset{\‾}{{\mathrm{\σ}}_{ij}}=c_{ijkl}{\mathrm{\ϵ}}_{kl}=c_{ijkl}\frac{1}{2}(u_{k,l}+u_{l,k})=\underset{\‾}{c_{ijkl}\·u_{k,l}}\⇒c_{ijkl}\·u_{k,li}=\mathrm{\ρ}{\stackrel{\·\·}{u}}_{j.}$

其中σ_ij为应力张量，ρ是密度，u_j是质点位移向量，t为时间，圆点代表时间的不同，指针从逗号后开始代表空间坐标的不同，重复指针的合计就如张量标志自动的被设定在自此之后；c_ijkl是弹性系数在之后可以将其收缩在C_IJ系数中，ε_kl是应变张量，其与质点位移有关；

根据弹性动力学理论，现假设一个平面谐波沿着x₁传播，其角频率为ω，该波相应的物理量在x₂方向上独立，波的形式u_k如下：

$(u_1,u_2,u_3)=(U,V,W)\·e^{j\mathrm{\ξ}(x_1+{\mathrm{\αx}}_3)}\·e^{-j\mathrm{\ω}t}$

其中ξ为波数，而α则表示一个未知系数，(ξ,0,αξ)为该波的传播向量，(U,V,W)为该平面波相应的振幅；将u_k带入控制波的传播方程得到：该线性代数形式[K_ij]_3×3其中：

K₁₁＝C₁₁+2·C₁₅·α+C₅₅α²-ρ₀·c²

K₂₂＝C₆₆+2·C₄₆·α+C₄₄·α²-ρ₀·c²

K₃₃＝C₅₅+2·C₃₅·α+C₃₃·α²-ρ₀·c²

K₁₂＝K₂₁＝C₁₆+(C₁₄+C₅₆)·α+C₄₅·α²

K₁₃＝K₃₁＝C₁₅+(C₁₃+C₅₅)·α+C₃₅·α²

K₂₃＝K₃₂＝C₅₆+(C₃₆+C₄₅)·α+C₃₄·α²

为了使ω，ξ存在值，α的值要使K矩阵行列式为零，也就是说α可以作为K矩阵的特征值，向量(U,V,W)是其相对相应的特征向量；

在各向同性材料薄板的波场中可以利用以下的弹性系数关系；首先我们将C_IJ定义在一个各向同性材料中，得到的弹性系数为:

$C_{IJ}=\left[\begin{array}{cccccc}\mathrm{\λ}+2\mathrm{\μ}& \mathrm{\λ}& \mathrm{\λ}& 0& 0& 0\\ & \mathrm{\λ}+2\mathrm{\μ}& \mathrm{\λ}& 0& 0& 0\\ & & \mathrm{\λ}+2\mathrm{\μ}& 0& 0& 0\\ & & & \mathrm{\μ}& 0& 0\\ & Sym.& & & \mathrm{\μ}& 0\\ & & & & & \mathrm{\μ}\end{array}\right]$

λ与μ为拉梅系数，我们定义C_IJ＝c_ijkl，其中ij→I或J，定义11→1，22→2，33→3，23或32→4，31或13→5，12或21→6；

获得α的表示式，从而可知该波的传播方式；通常，在各向同性材料中，α的值有四个；

$\left\{\begin{array}{c}{\mathrm{\α}}_1=-{\mathrm{\α}}_3=\sqrt{{\left(\frac{c}{C_L}\right)}^2-1}\\ {\mathrm{\α}}_2=-{\mathrm{\α}}_4=\sqrt{{\left(\frac{c}{C_T}\right)}^2-1}\end{array}\right.$

将不同的α值下对应的波的形式累加，得到以下位移、应力公式来表示波场；A_q为该平面波的振幅系数

$(u_1,u_2,u_3)=\underset{q=1}{\overset{4}{\Σ}}(U_q,0,W_q)\·A_q\·e^{i\mathrm{\ξ}(x_1+{\mathrm{\α}}_q{x}_3)}\·e^{-i\mathrm{\ω}t}$

$({\mathrm{\σ}}_{33},{\mathrm{\σ}}_{13},{\mathrm{\σ}}_{23})=\underset{q=1}{\overset{4}{\Σ}}(D_{1q},D_{2q},D_{3q})\·A_q\·e^{i\mathrm{\ξ}(x_1+{\mathrm{\α}}_q{x}_3)}\·e^{-i\mathrm{\ω}t}$

获得边界条件的系数矩阵：

在无外力为边界条件的自由载荷状态下，板厚为2h，得到量化的漏表面波频散曲线，其系数矩阵为[M]_4×4：

$\left\{\begin{array}{c}{\mathrm{\σ}}_{33}{|}_{x_3=+h}\\ {\mathrm{\σ}}_{13}{|}_{x_3=+h}\\ {\mathrm{\σ}}_{33}{|}_{x_3=-h}\\ {\mathrm{\σ}}_{13}{|}_{x_3=-h}\end{array}\right\}={\[M\]}_{4\×4}\left\{\begin{array}{c}{A}_1\\ A_2\\ A_3\\ A_4\end{array}\right\}=\left\{\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right\}$

在各向同性的平面板中相应的位移幅值比为以下值：

$\left\{\begin{array}{c}{U}_q=1\\ W_1={\mathrm{\α}}_1,W_2=-\frac{1}{{\mathrm{\α}}_2},W_3={\mathrm{\α}}_3,W_4=-\frac{1}{{\mathrm{\α}}_4}\end{array}\right.$

α_q已经得出为:

$\left\{\begin{array}{c}{\mathrm{\α}}_1=-{\mathrm{\α}}_3=\sqrt{{\left(\frac{c}{C_L}\right)}^2-1}\\ {\mathrm{\α}}_2=-{\mathrm{\α}}_4=\sqrt{{\left(\frac{c}{C_T}\right)}^2-1}\end{array}\right.$

在各向同性的平面板中相应的应变幅值比为以下值：

$\left\{\begin{array}{c}{D}_{1q}=\mathrm{\λ}+(\mathrm{\λ}+2\mathrm{\μ}{W}_q){\mathrm{\α}}_q\\ D_{2q}=\mathrm{\μ}(W_q+{\mathrm{\α}}_q)\end{array}\right.$

λ,μ是拉梅系数,可表示横波波速C_T，纵波波速C_L为以下形式:

$\left\{\begin{array}{c}{C}_L^2=\frac{\mathrm{\λ}+2\mathrm{\μ}}{\mathrm{\ρ}}\\ C_T^2=\frac{\mathrm{\μ}}{\mathrm{\ρ}}\end{array}\right.$

最后，可以得出系数矩阵[M]_4×4:

${\[M\]}_{4\×4}=\left[\begin{array}{cccc}{D}_{11}{e}^{+{\mathrm{i\ξh\α}}_1}& D_{12}{e}^{+{\mathrm{i\ξh\α}}_2}& D_{13}{e}^{+{\mathrm{i\ξh\α}}_3}& D_{14}{e}^{+{\mathrm{i\ξh\α}}_4}\\ D_{21}{e}^{+{\mathrm{i\ξh\α}}_1}& D_{22}{e}^{+{\mathrm{i\ξh\α}}_2}& D_{23}{e}^{+{\mathrm{i\ξh\α}}_3}& D_{24}{e}^{+{\mathrm{i\ξh\α}}_4}\\ D_{11}{e}^{-{\mathrm{i\ξh\α}}_1}& D_{12}{e}^{-{\mathrm{i\ξh\α}}_2}& D_{13}{e}^{-{\mathrm{i\ξh\α}}_3}& D_{14}{e}^{-{\mathrm{i\ξh\α}}_4}\\ D_{21}{e}^{-{\mathrm{i\ξh\α}}_1}& D_{22}{e}^{-{\mathrm{i\ξh\α}}_2}& D_{23}{e}^{-{\mathrm{i\ξh\α}}_3}& D_{24}{e}^{-{\mathrm{i\ξh\α}}_4}\end{array}\right]$

由边界条件[M]_4×4行列式可得四个方向的频散曲线，再得到其频散特征函数K，为了得到非零解，函数K的值应为0，C_L、C_T、ρ、h都为已知，改变f、c的值使目标函数值最小保证其尽量趋近于0，频散曲线，频散特征函数，仿真的目标函数Π_s分别如下：

$\left\{\begin{array}{c}{\mathrm{\σ}}_{33}{|}_{x_3=+h}\\ {\mathrm{\σ}}_{13}{|}_{x_3=+h}\\ {\mathrm{\σ}}_{33}{|}_{x_3=-h}\\ {\mathrm{\σ}}_{13}{|}_{x_3=-h}\end{array}\right\}={\[M\]}_{4\×4}\left\{\begin{array}{c}{A}_1\\ A_2\\ A_3\\ A_4\end{array}\right\}=\left\{\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right\}$

K(f,c；C_L,C_T,ρ,h)＝log₁₀(|[M]_4×4|)

${\mathrm{\Π}}_s=\underset{j=1}{\overset{N_s}{\Σ}}\[K(f_j^s,c_j^s;C_L^g,C_T^g,{\mathrm{\ρ}}^g,h^g)\]$

步骤1.2：确立波速提取的公式；

在单频激励/接收的情况下，上表面的直接反射回波I传播的时间与漏表面波L的传播时间分别为：

$t_1=\frac{2(R-\mathrm{\Δ}z)}{v_w}$

$t_2=\frac{2(R-\frac{\mathrm{\Δ}z}{{\mathrm{cos\θ}}_{SAW}})}{v_w}+\frac{2\·\mathrm{\Δ}z\·{\mathrm{tan\θ}}_{SAW}}{v_{SAW}}$

其中R为聚焦半径，Δz为散焦距离，v_w为水的超声波波速，θ_SAW为产生表面波的瑞利角，v_SAW为材料的表面波波速；因此两者的时间差为：

$\mathrm{\Δ}t=t_2-t_1=\frac{2(1-{\mathrm{cos\θ}}_{SAW})}{v_w}\·\mathrm{\Δ}z$

即：

${\mathrm{cos\θ}}_{SAW}=1-\frac{v_w\·\mathrm{\Δ}t}{2\·\mathrm{\Δ}z}$

将Snell定律：

或

代入后，可得：

$\frac{v_w}{v_{SAW}}=\sqrt{1-{(1-\frac{v_w}{2}\·\frac{\mathrm{\Δ}t}{\mathrm{\Δ}z})}^2}$

此时如果Δz恰为一个V(z)曲线的振荡周期时，1/Δt则为换能器的激励频率f；如果Δz能够确定，便可使用如下公式进行表面波波速的计算：

$v_{SAW}=v_w\·{\[1-{(1-\frac{v_w}{2\·f\·\mathrm{\Δ}z})}^2\]}^{-1/2}$

因此，测量被测材料的V(z)曲线振荡周期成为波速提取的重点；

步骤2：搭建测试系统；

为了方便散焦步进测量，搭建了一套进行散焦步进测量的测试系统；该测试系统主要包括试样(1)、水槽与水(2)、换能器(3)、移动平台(4)、脉冲激励/接收仪(5)、示波器(6)、GPIB总线(7)、PXI总控制系统(8)、移动伺服马达(9)、旋转轴(10)；其中，在移动平台(4)下面安装换能器(3)，换能器(3)与脉冲激励/接收仪(5)相连，脉冲激励/接收仪(5)与示波器(6)相连，示波器(6)通过GPIB总线(7)与PXI总控制系统(8)相连，PXI总控制系统(8)与移动伺服马达(9)相连，同时PXI总控制系统(8)与旋转轴(10)相连；

步骤3：聚焦面数据采集；

以长方体碳化钨为被测试样，其尺寸为40mm×40mm×10mm，将换能器3聚焦到试样的上表面，通过脉冲激励/接收仪(5)在发出一个带宽为10-200MHz的脉冲后转换为接收状态，当接收到反射信号后，将信号传输进示波器(6)，示波器的采样频率f_S＝2.5GHz，采样点数N_s＝10000；经过示波器的低通滤波后，通过GPIB总线(7)存储进PXI总控制系统(8)；

步骤4：散焦测量；

将换能器朝试样方向移动Δz₀＝10μm，待移动完成后进行电压数据采集，采集结束后再将换能器朝试样方向移动Δz₀＝10μm进行数据采集，采样频率f_S＝2.5GHz，采样点数N_s＝10000，如此循环往复，共移动4mm，因此将得到400组电压数据，将聚焦面的电压数据包含在内共得到M＝401组电压数据；

步骤5：时域傅里叶变换；

将测得的数据进行时域傅里叶变换；

$A_i\[k\]=\underset{n=0}{\overset{N_s-1}{\Σ}}{x}_i\[n\]e^{-j2\mathrm{\π}nk/N_s}$

其中：A_i为时域傅里叶变换后的频谱值，x_i代表一组电压数据，i＝0,1,2…M-1，k＝0,1,2…N_s-1，j代表虚部，N_s＝10000；即所得A_i[k]，i＝0,1,2…M-1，k＝0,1,2…N_s-1；

步骤6：空间傅里叶变换；

为了得到精确的振荡周期Δz，需要对时域傅里叶变换的结果再进行沿散焦距离方向的空间傅里叶变换，将散焦距离z变换至z^-1域：

$B_i\[k\]=\underset{m=0}{\overset{M-1}{\Σ}}{A}_m\[k\]e^{-j2\mathrm{\π}mi/M}$

其中：B_i为空间傅里叶变换后的频谱值，A_m代表沿散焦方向的时域傅里叶变换的频谱值，i＝0,1,2…M-1，k＝0,1,2…N_s-1，M＝401，j代表虚部；所得B_i[k]，i＝0,1,2…M-1，k＝0,1,2…N_s-1；

步骤7：模态追踪；

对2.5-22.5MHz范围内的峰值进行追踪，即可找出该频率段连续的Δz值；

步骤8：波速提取；

将水中的超声波波速v_W＝1500m/s，每一个峰值对应的频率与Δz带入公式即可得到该频率段内连续的表面波波速；

步骤9：单纯形法获得横、纵波波速；

通过单纯形法使目标函数残差绝对值最接近零，单纯形法改变C_L、C_T、ρ、h使目标函数值的达到所需范围；

单纯形搜索法通过构造单纯形来逼近极小点，每构造一个单纯形，确定其最高点和最低点，然后通过扩展或压缩、反射构造新的单纯形，目的是使得极小点能包含于单纯形内。

2.根据权利要求1所述的一种基于单纯形法的各向同性薄板材料弹性性质获取方法，其特征在于：本方法中未知量有四个(C_L,C_T,ρ,h)，为四维变量问题；

用单纯形搜索法求无约束问题minF(x)，x∈Rⁿ的算法步骤如下：

①选取初始单纯形{x⁰,x¹,···,xⁿ}，反映系数α＞1，紧缩系数θ∈(0,1)，扩展系数γ＞1，收缩系数β∈(0，1)及精度ε>0，置k＝0；

②将单纯形的n+1个顶点按目标函数值的大小重新编号，使顶点的编号满足F(x⁰)≤F(x¹)≤···≤F(x^n-1)≤F(xⁿ)；

③令若停止迭代输出x⁰，否则转入④；

④计算xⁿ⁺²＝xⁿ⁺¹+α(xⁿ⁺¹-xⁿ)，若F(xⁿ⁺²)＜F(x⁰)，转⑤，否则当F(xⁿ⁺²)＜F(x^n-1)时转⑥，当F(xⁿ⁺²)≥F(x^n-1)转⑦；

⑤计算xⁿ⁺³＝xⁿ⁺¹+γ(xⁿ⁺²-xⁿ⁺¹)，若F(xⁿ⁺³)＜F(x⁰)，令xⁿ＝xⁿ⁺³，转②，否则转⑥；

⑥令xⁿ＝xⁿ⁺²，转②

⑦令xⁿ＝{xⁱ|F(xⁱ)＝min(F(xⁿ),F(xⁿ⁺²))}，计算xⁿ⁺⁴＝xⁿ⁺¹+β(xⁿ-xⁿ⁺¹)，若F(xⁿ⁺⁴)＜F(xⁿ)，令xⁿ＝xⁿ⁺⁴，转②，否则转⑧；

⑧令x^j＝x⁰+θ(x^j-x⁰)，j＝0,1,···,n，转②。