CN103926329A - 一种基于单纯形法的半无限域基体涂层结构材料弹性性质获取方法 - Google Patents

Abstract

一种基于单纯形法的半无限域基体涂层结构材料弹性性质获取方法，材料的弹性系数与材料制备工艺密切相关；利用超声波对材料弹性性质进行测量是无损检测领域很有前景的测量方法之一。基于声学显微镜技术，自行开发的涂层材料弹性系数的超声测量系统，采用线聚焦PVDF探头，通过纵波和表面波波速的同时测量，可实现材料的弹性系数无损检测。本方法具有如下有益效果，不必再将整幅理论与实际曲线的波速平方的残差总和最小值都计算出来；不必把所有的理论频散曲线都确定，这种方式耗时且数据量冗长，单纯形的方法可以加速并化简程序。

Description

一种基于单纯形法的半无限域基体涂层结构材料弹性性质获取方法

技术领域

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

背景技术

表面工程是再造技术的重要组成部分，同时又为先进制造技术的发展提供了很好的技术支撑，表面工程是指利用物理、化学或其他方法在材料表面制备一层不同于基体材料且有一定厚度的强化、防护覆盖层。涂层材料广泛应用于航空航天、机械、石油、化工、核电及微器件与微制造等领域，其主要功能是防腐、耐摩擦、抗氧化等表面防护。鉴于涂层在工业领域应用越来越广泛，现代工业对装备可靠性、安全性的要求不断提高，涂层强度和失效分析越来越重要，同时涂层的力学性能的表征可以及时对工业镀覆技术进行改进提供参考。

声速和材料显微结构参数中一部分参数相关，超声波声速在材料中传播速度的变化通常与显微结构特性和力学性能的显著变化相关，同时由于超声波具有穿透能力强，频带宽和可实现对材料的无损检测等优势，超声波速法在材料质量检测和评价中得到广泛应用。

当超声波波长小于基底厚度时，视为半无限域基体涂层结构材料，通过小尺寸材料弹性常数测量系统，根据V(f，z)频域分析法，采用线聚焦PVDF探头，进行频散分析。将纵波和表面波波速的同时测量，可实现材料的弹性系数无损检测。经由一种新的反演算法，利用该线聚焦超声探头来获取半无限域基体涂层结构材料的表面波波速。该方法基于单纯形法诱导目标函数于频散特征方程的系数矩阵行列式中，半无限域基体涂层结构材料的弹性性质与试样密度都可以通过声学性质及所测密度来获得。

发明内容

本发明的目的是为了解决半无限域基体涂层结构材料波速提取的问题，提出一种先进的材料波速提取方法。

为实现上述目的，本发明采用的技术方案为一种基于单纯形法的半无限域基体涂层结构材料弹性性质获取方法，该方法的实现步骤如下，

步骤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四个变量值，在其保证目标函数值最小即接近0时，即确定各值。

本发明具有以下优点：（1）不必再将整幅理论与实际曲线的波速平方的残差总和最小值都计算出来。（2）不必把所有的理论频散曲线都确定，这种方式耗时且数据量冗长，单纯形的方法可以加速并化简程序。

附图说明

图1为为本方法的实施方法流程图。

图2为散焦测量系统示意图。

图3为表面波传播示意图。

图4为Lamb波在半无限域基体涂层结构材料中传播示意图。

图5为聚焦面时域波形图。

图6为7.5MHz频率下V(z)振荡曲线图。

图7为7.5MHz频率下z^-1域曲线图。

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

具体实施方式

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

步骤1确立仿真目标函数

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

步骤1.1：根据运动方程：σ_ij,i+ρf_j＝ρü_j（可忽略重力ρf_j影响）；及本构方程：σ_ij＝c_ijklε_kl；几何方程：推知控制波的传播方程c_i·,＝ρ^..。

$\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_ijkl·u_k,li＝ρü_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(1\right)}=\left[\begin{array}{c}\begin{array}{cccccc}{\mathrm{\λ}}^{\left(1\right)}+2{\mathrm{\μ}}^{\left(1\right)}& {\mathrm{\λ}}^{\left(1\right)}& {\mathrm{\λ}}^{\left(1\right)}& 0& 0& 0\end{array}\\ \begin{array}{ccccc}{\mathrm{\λ}}^{\left(1\right)}+2{\mathrm{\μ}}^{\left(1\right)}& {\mathrm{\λ}}^{\left(1\right)}& 0& 0& 0\end{array}\\ \begin{array}{cccc}{\mathrm{\λ}}^{\left(1\right)}+2{\mathrm{\μ}}^{\left(1\right)}& 0& 0& 0\end{array}\\ \begin{array}{ccc}{\mathrm{\μ}}^{\left(1\right)}& 0& 0\end{array}\\ \begin{array}{ccc}\mathrm{Sym}.& {\mathrm{\μ}}^{\left(1\right)}& 0\end{array}\\ {\mathrm{\μ}}^{\left(1\right)}\end{array}\right]$

$C_{\mathrm{IJ}}^{\left(2\right)}=\left[\begin{array}{c}\begin{array}{cccccc}{\mathrm{\λ}}^{\left(2\right)}+2{\mathrm{\μ}}^{\left(2\right)}& {\mathrm{\λ}}^{\left(2\right)}& {\mathrm{\λ}}^{\left(2\right)}& 0& 0& 0\end{array}\\ \begin{array}{ccccc}{\mathrm{\λ}}^{\left(2\right)}+2{\mathrm{\μ}}^{\left(2\right)}& {\mathrm{\λ}}^{\left(2\right)}& 0& 0& 0\end{array}\\ \begin{array}{cccc}{\mathrm{\λ}}^{\left(2\right)}+2{\mathrm{\μ}}^{\left(2\right)}& 0& 0& 0\end{array}\\ \begin{array}{ccc}{\mathrm{\μ}}^{\left(2\right)}& 0& 0\end{array}\\ \begin{array}{ccc}\mathrm{Sym}.& {\mathrm{\μ}}^{\left(2\right)}& 0\end{array}\\ {\mathrm{\μ}}^{\left(2\right)}\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{\mathrm{\ω}/c_p}{\mathrm{\ξ}}\right)}^2}-1\\ {\mathrm{\α}}_2,{\mathrm{\α}}_4=\±\sqrt{{\left(\frac{\mathrm{\ω}/c_t}{\mathrm{\ξ}}\right)}^2}-1\end{array}\right.$

c_p为纵波波速，c_t为横波波速。将不同的α值下对应的波的形式累加，得到以下位移、应力公式来表示波场。其中A_q为未知系数。

可解得特征值，这些特征值是给定ω，ξ情况下波运动的基本解，波在半无限域基体涂层结构材料中的位移场可以通过特征值的线性组合表示，不同是，基底中没有通过半无限域边界反射回来的波。可得到以下位移、应力公式来表示波场：

$u_k^{\left(1\right)}=U_{\mathrm{kq}}^{\left(1\right)}\·A_q^{\left(1\right)}\·e^{\mathrm{j\ξ}(x_1+{\mathrm{\α}}_q^{\left(1\right)}{x}_3)}\·e^{-\mathrm{j\ωt}}$ q＝1,2,3,4

$u_k^{\left(2\right)}=U_{\mathrm{kq}}^{\left(2\right)}\·A_q^{\left(2\right)}\·e^{\mathrm{j\ξ}(x_1+{\mathrm{\α}}_q^{\left(2\right)}{x}_3)}\·e^{-\mathrm{j\ωt}}$ q＝1,2

${\mathrm{\σ}}_{\mathrm{ij}}^{\left(1\right)}=c_{\mathrm{ijkl}}^{\left(1\right)}\·u_{k,l}^{\left(1\right)}$

q＝1,2,3,4

${\mathrm{\σ}}_{\mathrm{ij}}^{\left(2\right)}=c_{\mathrm{ijkl}}^{\left(2\right)}\·u_{k,l}^{\left(2\right)}$

$\left[\begin{array}{c}{\mathrm{\σ}}_{31}^{\left(2\right)}\\ {\mathrm{\σ}}_{32}^{\left(2\right)}\\ {\mathrm{\σ}}_{33}^{\left(2\right)}\end{array}\right]=\left[\begin{array}{c}{D}_{1q}^{\left(2\right)}\\ D_{2q}^{\left(2\right)}\\ D_{3q}^{\left(2\right)}\end{array}\right]\·A_q^{\left(2\right)}\·e^{\mathrm{j\ξ}(x_1+{\mathrm{\α}}_q^{\left(2\right)}{x}_3)}\·e^{-\mathrm{j\ωt}}$ q＝1,2

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

在x₃＝d为自由载荷状态（T₃₃＝T₁₃=0），根据位移与应力在x₃＝0处的边界条件，得到以下关系式：

${\mathrm{\σ}}_{31}^{\left(1\right)}{|}_{x_3=d}=0$

${\mathrm{\σ}}_{33}^{\left(1\right)}{|}_{x_3=d}=0$

$u_1^{\left(1\right)}{|}_{x_3=0}=u_1^{\left(2\right)}{|}_{x_3=0}$

$u_3^{\left(1\right)}{|}_{x_3=0}=u_3^{\left(2\right)}{|}_{x_3=0}$

${\mathrm{\σ}}_{31}^{\left(1\right)}{|}_{x_2=0}={\mathrm{\σ}}_{31}^{\left(2\right)}{|}_{x_2=0}$

${\mathrm{\σ}}_{33}^{\left(1\right)}{|}_{x_2=0}={\mathrm{\σ}}_{33}^{\left(1\right)}{|}_{x_2=0}$

并最终可得： $\left[\begin{array}{cccccc}{D}_{11}^{\left(1\right)}& D_{12}^{\left(1\right)}& D_{13}^{\left(1\right)}& D_{14}^{\left(1\right)}& 0& 0\\ D_{31}^{\left(1\right)}& D_{32}^{\left(1\right)}& D_{33}^{\left(1\right)}& D_{34}^{\left(1\right)}& 0& 0\\ U_{11}^{\left(1\right)}& U_{12}^{\left(1\right)}& U_{13}^{\left(1\right)}& U_{14}^{\left(1\right)}& {-U}_{11}^{\left(2\right)}& {-U}_{12}^{\left(2\right)}\\ U_{31}^{\left(1\right)}& U_{32}^{\left(1\right)}& U_{33}^{\left(1\right)}& U_{34}^{\left(1\right)}& {-U}_{31}^{\left(2\right)}& {-U}_{32}^{\left(2\right)}\\ D_{11}^{\left(1\right)}& D_{12}^{\left(1\right)}& D_{13}^{\left(1\right)}& D_{14}^{\left(1\right)}& {-D}_{11}^{\left(2\right)}& {-D}_{12}^{\left(2\right)}\\ D_{31}^{\left(1\right)}& D_{32}^{\left(1\right)}& D_{33}^{\left(1\right)}& D_{34}^{\left(1\right)}& {-D}_{31}^{\left(2\right)}& {-D}_{32}^{\left(2\right)}\end{array}\right]\left[\begin{array}{c}{A}_1^{\left(1\right)}\\ A_2^{\left(1\right)}\\ A_3^{\left(1\right)}\\ A_4^{\left(1\right)}\\ A_1^{\left(2\right)}\\ A_2^{\left(2\right)}\end{array}\right]=0$

$U_{\mathrm{kq}}^{\left(1\right)}=\left[\begin{array}{cccc}1& 1& 1& 1\\ 0& 0& 0& 0\\ {\mathrm{\α}}_1^{\left(1\right)}& \frac{1}{-{\mathrm{\α}}_2^{\left(1\right)}}& {\mathrm{\α}}_3^{\left(1\right)}& \frac{1}{-{\mathrm{\α}}_4^{\left(1\right)}}\end{array}\right]\begin{array}{cc}\←& U_q^{\left(1\right)}\\ \←& V_q^{\left(1\right)}\\ \←& W_q^{\left(1\right)}\end{array}$

此时： $U_{\mathrm{kq}}^{\left(2\right)}=\left[\begin{array}{cc}1& 1\\ 0& 0\\ {\mathrm{\α}}_1^{\left(2\right)}& \frac{1}{-{\mathrm{\α}}_2^{\left(2\right)}}\end{array}\right]\begin{array}{cc}\←& U_q^{\left(1\right)}\\ \←& V_q^{\left(1\right)}\\ \←& W_q^{\left(1\right)}\end{array}$

根据： $W_q=\frac{{\mathrm{\ρc}}^2-C_{11}-C_{55}{\mathrm{\α}}_q^2}{(C_{13}+C_{55}){\mathrm{\α}}_q}$

可得：

$W_1\begin{array}{c}=\frac{\mathrm{\ρ}{c}^2-C_{11}-C_{55}{\mathrm{\α}}_1^2}{(C_{13}+C_{55}){\mathrm{\α}}_1}=\frac{({\mathrm{\α}}_1^2+1)\mathrm{\ρ}{c}_p^2-(\mathrm{\λ}+2\mathrm{\μ})-\mathrm{\μ}{\mathrm{\α}}_1^2}{(\mathrm{\λ}+\mathrm{\μ}){\mathrm{\α}}_1}=\frac{({\mathrm{\α}}_1^2+1)(\mathrm{\λ}+2\mathrm{\μ})-(\mathrm{\λ}+2\mathrm{\μ})-\mathrm{\μ}{\mathrm{\α}}_1^2}{(\mathrm{\λ}+\mathrm{\μ}){\mathrm{\α}}_1}\\ ={\mathrm{\α}}_1\end{array}$

$W_2\begin{array}{c}=\frac{\mathrm{\ρ}{c}^2-C_{11}-C_{55}{\mathrm{\α}}_2^2}{(C_{13}+C_{55}){\mathrm{\α}}_2}=\frac{({\mathrm{\α}}_2^2+1)\mathrm{\ρ}{c}_s^2-(\mathrm{\λ}+2\mathrm{\μ})-\mathrm{\μ}{\mathrm{\α}}_2^2}{(\mathrm{\λ}+\mathrm{\μ}){\mathrm{\α}}_2}=\frac{({\mathrm{\α}}_2^2+1)\mathrm{\μ}-(\mathrm{\λ}+2\mathrm{\μ})-\mathrm{\μ}{\mathrm{\α}}_2^2}{(\mathrm{\λ}+\mathrm{\μ}){\mathrm{\α}}_2}\\ =-\frac{1}{{\mathrm{\α}}_2}\end{array}$

W₃＝α₃＝-α₁

$W_4=-\frac{1}{{\mathrm{\α}}_4}=\frac{1}{{\mathrm{\α}}_2}$

经由D_pq表达式，整理、去除为零部分后可得：

D_1q＝C₅₅(W_q+α_q)

D_3q＝C₃₁+(C₃₃W_q)α_q

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

$\left\{\begin{array}{c}\left(\begin{array}{c}{\mathrm{\σ}}_{31}{|}_{x_3=d}\\ {\mathrm{\σ}}_{33}^{\left(1\right)}{|}_{x_3=d}\end{array}\right)\\ \begin{array}{cc}{u}_1^{\left(1\right)}{|}_{x_3=0}& -u_1^{\left(2\right)}\end{array}{|}_{x_3=0}\\ \begin{array}{cc}{u}_3^{\left(1\right)}{|}_{x_3=0}& =u_3^{\left(2\right)}{|}_{x_3=0}\end{array}\\ \begin{array}{cc}{\mathrm{\σ}}_{31}^{\left(1\right)}{|}_{x_3=0}& ={\mathrm{\σ}}_{31}^{\left(2\right)}{|}_{x_3=0}\end{array}\\ \begin{array}{cc}{\mathrm{\σ}}_{33}^{\left(1\right)}{|}_{x_3=0}& ={\mathrm{\σ}}_{33}^{\left(1\right)}{|}_{x_3=0}\end{array}\end{array}\right\}={\left[M\right]}_{6\×6}\left\{\begin{array}{c}{A}_1\\ A_2\\ A_3\\ A_4\\ A_5\\ A_6\end{array}\right\}=\left\{0\right\}$

K(f,c;C_L,C_T,ρ,h)＝log₁₀([M]_6×6)

${\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：确立波速提取的公式。

在单频激励/接收的情况下，图3所示的漏表面波传播示意图中，上表面的直接反射回波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{\Δt}}$

将Snell定律：

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

代入(4)后，可得：

$\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\·f\·\mathrm{\Δz}})}^2]}^{-1/2}$

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

步骤2：搭建测试系统。

为了方便散焦步进测量，搭建了一套进行散焦步进测量的测试系统，如图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总控制系统，聚焦面的时域波形如图5所示。

步骤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频率下的振荡曲线如图6所示。

步骤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域的曲线如图7所示。

步骤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ⁿ)；

③令若 ${\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/1}\≤\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为通过的仿真频散曲线得到的数据点数，这样只有那些相应残差范围内的K被累加计算起来，就不必再将整幅理论与实际曲线的波速平方的残差总和最小值都计算出来。并且以前的方法，必须把所有的理论频散曲线都确定，这种方式耗时且数据量冗长，因此，单纯形的方法可以加速并化简程序。

Claims (3)

1.一种基于单纯形法的半无限域基体涂层结构材料弹性性质获取方法，其特征在于：该方法的实现步骤包括如下，

步骤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）、水槽与水（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四个变量值，在其保证目标函数值最小时，即确定各值。

2.根据权利要求1所述的一种基于单纯形法的半无限域基体涂层结构材料弹性性质获取方法，其特征在于：通过对各项同性材料试件内部部分波线性组合，将其定义在一个关于边界条件的特征方程中，从而得到频散关系；薄膜厚度为d的半无限域基体试样置于笛卡尔坐标系下，x1-x3平面是Lamb波传播的矢面，而坐标x1代表波的传播方向；全局坐标（x1,x2,x3）为下面的理论服务；

步骤1.1：根据运动方程：σ_ij,i+ρf_j＝ρü_j；及本构方程：σ_ij＝c_ijklε_kl；几何方程：推知控制波的传播方程c_ijkl·u_k,li＝ρü_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带入控制波的传播方程ci_jkl·u_k,li＝ρü_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_15·α+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(1\right)}=\left[\begin{array}{c}\begin{array}{cccccc}{\mathrm{\λ}}^{\left(1\right)}+2{\mathrm{\μ}}^{\left(1\right)}& {\mathrm{\λ}}^{\left(1\right)}& {\mathrm{\λ}}^{\left(1\right)}& 0& 0& 0\end{array}\\ \begin{array}{ccccc}{\mathrm{\λ}}^{\left(1\right)}+2{\mathrm{\μ}}^{\left(1\right)}& {\mathrm{\λ}}^{\left(1\right)}& 0& 0& 0\end{array}\\ \begin{array}{cccc}{\mathrm{\λ}}^{\left(1\right)}+2{\mathrm{\μ}}^{\left(1\right)}& 0& 0& 0\end{array}\\ \begin{array}{ccc}{\mathrm{\μ}}^{\left(1\right)}& 0& 0\end{array}\\ \begin{array}{ccc}\mathrm{Sym}.& {\mathrm{\μ}}^{\left(1\right)}& 0\end{array}\\ {\mathrm{\μ}}^{\left(1\right)}\end{array}\right]$

$C_{\mathrm{IJ}}^{\left(2\right)}=\left[\begin{array}{c}\begin{array}{cccccc}{\mathrm{\λ}}^{\left(2\right)}+2{\mathrm{\μ}}^{\left(2\right)}& {\mathrm{\λ}}^{\left(2\right)}& {\mathrm{\λ}}^{\left(2\right)}& 0& 0& 0\end{array}\\ \begin{array}{ccccc}{\mathrm{\λ}}^{\left(2\right)}+2{\mathrm{\μ}}^{\left(2\right)}& {\mathrm{\λ}}^{\left(2\right)}& 0& 0& 0\end{array}\\ \begin{array}{cccc}{\mathrm{\λ}}^{\left(2\right)}+2{\mathrm{\μ}}^{\left(2\right)}& 0& 0& 0\end{array}\\ \begin{array}{ccc}{\mathrm{\μ}}^{\left(2\right)}& 0& 0\end{array}\\ \begin{array}{ccc}\mathrm{Sym}.& {\mathrm{\μ}}^{\left(2\right)}& 0\end{array}\\ {\mathrm{\μ}}^{\left(2\right)}\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{\mathrm{\ω}/c_p}{\mathrm{\ξ}}\right)}^2}-1\\ {\mathrm{\α}}_2,{\mathrm{\α}}_4=\±\sqrt{{\left(\frac{\mathrm{\ω}/c_t}{\mathrm{\ξ}}\right)}^2}-1\end{array}\right.$

c_p为纵波波速，c_t为横波波速。将不同的α值下对应的波的形式累加，得到以下位移、应力公式来表示波场；其中A_q为未知系数；

可解得特征值，这些特征值是给定ω，ξ情况下波运动的基本解，波在半无限域基体涂层结构材料中的位移场可以通过特征值的线性组合表示，不同是，基底中没有通过半无限域边界反射回来的波；可得到以下位移、应力公式来表示波场：

$u_k^{\left(1\right)}=U_{\mathrm{kq}}^{\left(1\right)}\·A_q^{\left(1\right)}\·e^{\mathrm{j\ξ}(x_1+{\mathrm{\α}}_q^{\left(1\right)}{x}_3)}\·e^{-\mathrm{j\ωt}}$ q＝1,2,3,4

$u_k^{\left(2\right)}=U_{\mathrm{kq}}^{\left(2\right)}\·A_q^{\left(2\right)}\·e^{\mathrm{j\ξ}(x_1+{\mathrm{\α}}_q^{\left(2\right)}{x}_3)}\·e^{-\mathrm{j\ωt}}$ q＝1,2

${\mathrm{\σ}}_{\mathrm{ij}}^{\left(1\right)}=c_{\mathrm{ijkl}}^{\left(1\right)}\·u_{k,l}^{\left(1\right)}$

$\left[\begin{array}{c}{\mathrm{\σ}}_{31}^{\left(1\right)}\\ {\mathrm{\σ}}_{32}^{\left(1\right)}\\ {\mathrm{\σ}}_{33}^{\left(1\right)}\end{array}\right]=\left[\begin{array}{c}{D}_{1q}^{\left(1\right)}\\ D_{2q}^{\left(1\right)}\\ D_{3q}^{\left(1\right)}\end{array}\right]\·A_q^{\left(1\right)}\·e^{\mathrm{j\ξ}(x_1+{\mathrm{\α}}_q^{\left(1\right)}{x}_3)}\·e^{-\mathrm{j\ωt}}$ q＝1,2,3,4

${\mathrm{\σ}}_{\mathrm{ij}}^{\left(2\right)}=c_{\mathrm{ijkl}}^{\left(2\right)}\·u_{k,l}^{\left(2\right)}$

$\left[\begin{array}{c}{\mathrm{\σ}}_{31}^{\left(2\right)}\\ {\mathrm{\σ}}_{32}^{\left(2\right)}\\ {\mathrm{\σ}}_{33}^{\left(2\right)}\end{array}\right]=\left[\begin{array}{c}{D}_{1q}^{\left(2\right)}\\ D_{2q}^{\left(2\right)}\\ D_{3q}^{\left(2\right)}\end{array}\right]\·A_q^{\left(2\right)}\·e^{\mathrm{j\ξ}(x_1+{\mathrm{\α}}_q^{\left(2\right)}{x}_3)}\·e^{-\mathrm{j\ωt}}$ q＝1,2

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

在x₃＝d为自由载荷状态（T₃₃＝T₁₃=0），根据位移与应力在x₃＝0处的边界条件，得到以下关系式：

${\mathrm{\σ}}_{31}^{\left(1\right)}{|}_{x_3=d}=0$

${\mathrm{\σ}}_{33}^{\left(1\right)}{|}_{x_3=d}=0$

$u_1^{\left(1\right)}{|}_{x_3=0}=u_1^{\left(2\right)}{|}_{x_3=0}$

$u_3^{\left(1\right)}{|}_{x_3=0}=u_3^{\left(2\right)}{|}_{x_3=0}$

${\mathrm{\σ}}_{31}^{\left(1\right)}{|}_{x_2=0}={\mathrm{\σ}}_{31}^{\left(2\right)}{|}_{x_2=0}$

${\mathrm{\σ}}_{33}^{\left(1\right)}{|}_{x_2=0}={\mathrm{\σ}}_{33}^{\left(1\right)}{|}_{x_2=0}$

并最终可得： $\left[\begin{array}{cccccc}{D}_{11}^{\left(1\right)}& D_{12}^{\left(1\right)}& D_{13}^{\left(1\right)}& D_{14}^{\left(1\right)}& 0& 0\\ D_{31}^{\left(1\right)}& D_{32}^{\left(1\right)}& D_{33}^{\left(1\right)}& D_{34}^{\left(1\right)}& 0& 0\\ U_{11}^{\left(1\right)}& U_{12}^{\left(1\right)}& U_{13}^{\left(1\right)}& U_{14}^{\left(1\right)}& {-U}_{11}^{\left(2\right)}& {-U}_{12}^{\left(2\right)}\\ U_{31}^{\left(1\right)}& U_{32}^{\left(1\right)}& U_{33}^{\left(1\right)}& U_{34}^{\left(1\right)}& {-U}_{31}^{\left(2\right)}& {-U}_{32}^{\left(2\right)}\\ D_{11}^{\left(1\right)}& D_{12}^{\left(1\right)}& D_{13}^{\left(1\right)}& D_{14}^{\left(1\right)}& {-D}_{11}^{\left(2\right)}& {-D}_{12}^{\left(2\right)}\\ D_{31}^{\left(1\right)}& D_{32}^{\left(1\right)}& D_{33}^{\left(1\right)}& D_{34}^{\left(1\right)}& {-D}_{31}^{\left(2\right)}& {-D}_{32}^{\left(2\right)}\end{array}\right]\left[\begin{array}{c}{A}_1^{\left(1\right)}\\ A_2^{\left(1\right)}\\ A_3^{\left(1\right)}\\ A_4^{\left(1\right)}\\ A_1^{\left(2\right)}\\ A_2^{\left(2\right)}\end{array}\right]=0$

$U_{\mathrm{kq}}^{\left(1\right)}=\left[\begin{array}{cccc}1& 1& 1& 1\\ 0& 0& 0& 0\\ {\mathrm{\α}}_1^{\left(1\right)}& \frac{1}{-{\mathrm{\α}}_2^{\left(1\right)}}& {\mathrm{\α}}_3^{\left(1\right)}& \frac{1}{-{\mathrm{\α}}_4^{\left(1\right)}}\end{array}\right]\begin{array}{cc}\←& U_q^{\left(1\right)}\\ \←& V_q^{\left(1\right)}\\ \←& W_q^{\left(1\right)}\end{array}$

此时： $U_{\mathrm{kq}}^{\left(2\right)}=\left[\begin{array}{cc}1& 1\\ 0& 0\\ {\mathrm{\α}}_1^{\left(2\right)}& \frac{1}{-{\mathrm{\α}}_2^{\left(2\right)}}\end{array}\right]\begin{array}{cc}\←& U_q^{\left(1\right)}\\ \←& V_q^{\left(1\right)}\\ \←& W_q^{\left(1\right)}\end{array}$

根据： $W_q=\frac{{\mathrm{\ρc}}^2-C_{11}-C_{55}{\mathrm{\α}}_q^2}{(C_{13}+C_{55}){\mathrm{\α}}_q}$

可得：

$W_1\begin{array}{c}=\frac{\mathrm{\ρ}{c}^2-C_{11}-C_{55}{\mathrm{\α}}_1^2}{(C_{13}+C_{55}){\mathrm{\α}}_1}=\frac{({\mathrm{\α}}_1^2+1)\mathrm{\ρ}{c}_p^2-(\mathrm{\λ}+2\mathrm{\μ})-\mathrm{\μ}{\mathrm{\α}}_1^2}{(\mathrm{\λ}+\mathrm{\μ}){\mathrm{\α}}_1}=\frac{({\mathrm{\α}}_1^2+1)(\mathrm{\λ}+2\mathrm{\μ})-(\mathrm{\λ}+2\mathrm{\μ})-\mathrm{\μ}{\mathrm{\α}}_1^2}{(\mathrm{\λ}+\mathrm{\μ}){\mathrm{\α}}_1}\\ ={\mathrm{\α}}_1\end{array}$

$W_2\begin{array}{c}=\frac{\mathrm{\ρ}{c}^2-C_{11}-C_{55}{\mathrm{\α}}_2^2}{(C_{13}+C_{55}){\mathrm{\α}}_2}=\frac{({\mathrm{\α}}_2^2+1)\mathrm{\ρ}{c}_s^2-(\mathrm{\λ}+2\mathrm{\μ})-\mathrm{\μ}{\mathrm{\α}}_2^2}{(\mathrm{\λ}+\mathrm{\μ}){\mathrm{\α}}_2}=\frac{({\mathrm{\α}}_2^2+1)\mathrm{\μ}-(\mathrm{\λ}+2\mathrm{\μ})-\mathrm{\μ}{\mathrm{\α}}_2^2}{(\mathrm{\λ}+\mathrm{\μ}){\mathrm{\α}}_2}\\ =-\frac{1}{{\mathrm{\α}}_2}\end{array}$

W₃＝α₃＝-α₁

$W_4=-\frac{1}{{\mathrm{\α}}_4}=\frac{1}{{\mathrm{\α}}_2}$

经由D_pq表达式，整理、去除为零部分后可得：

D_1q＝C₅₅(W_q+α_q)

D_3q＝C₃₁+(C₃₃W_q)α_q

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

$\left\{\begin{array}{c}\left(\begin{array}{c}{\mathrm{\σ}}_{31}{|}_{x_3=d}\\ {\mathrm{\σ}}_{33}^{\left(1\right)}{|}_{x_3=d}\end{array}\right)\\ \begin{array}{cc}{u}_1^{\left(1\right)}{|}_{x_3=0}& -u_1^{\left(2\right)}\end{array}{|}_{x_3=0}\\ \begin{array}{cc}{u}_3^{\left(1\right)}{|}_{x_3=0}& =u_3^{\left(2\right)}{|}_{x_3=0}\end{array}\\ \begin{array}{cc}{\mathrm{\σ}}_{31}^{\left(1\right)}{|}_{x_3=0}& ={\mathrm{\σ}}_{31}^{\left(2\right)}{|}_{x_3=0}\end{array}\\ \begin{array}{cc}{\mathrm{\σ}}_{33}^{\left(1\right)}{|}_{x_3=0}& ={\mathrm{\σ}}_{33}^{\left(1\right)}{|}_{x_3=0}\end{array}\end{array}\right\}={\left[M\right]}_{6\×6}\left\{\begin{array}{c}{A}_1\\ A_2\\ A_3\\ A_4\\ A_5\\ A_6\end{array}\right\}=\left\{0\right\}$

K(f,c;C_L,C_T,ρ,h)＝log₁₀([M]_6×6)

${\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：确立波速提取的公式；

在单频激励/接收的情况下，上表面的直接反射回波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{\Δt}}$

将Snell定律：

$\sin {\mathrm{\θ}}_{\mathrm{SAW}}=\frac{v_w}{v_{\mathrm{SAW}}}$ 或 ${\mathrm{\θ}}_{\mathrm{SAW}}={\sin }^{-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\·f\·\mathrm{\Δz}})}^2]}^{-1/2}$

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

3.根据权利要求1所述的一种基于单纯形法的半无限域基体涂层结构材料弹性性质获取方法，其特征在于：通过单纯形法使目标函数残差绝对值最接近零，单纯形法改变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ⁿ)；

③令，若 ${\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/1}\≤\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，转②。