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

权利要求书

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

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

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

Πs=Σj=1Ns[K(fjs,cjs;CLg,CTg,ρg,hg)]]]>

表示第j个点，模拟值为s的频率；表示第j个点，模拟值为s的表面波波速；

步骤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：时域傅里叶变换；

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

Ai[k]=Σn=0Ns-1xi[n]e-j2πnk/Ns]]>

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

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

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

Bi[k]=Σm=0M-1Am[k]e-j2πmi/M]]> 1

其中，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)曲线理论，根据如下公式进行波速的计算：

vLamb=vw·[1-(1-vw2·f·Δz)2]-1/2]]>

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

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

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

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

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

σij,i=ρu··jσij=cijklϵklϵij=12(ui,j+uj,i)⇒σij‾=cijklϵkl=cijkl12(uk,l+ul,k)=cijkl·uk,l‾⇒cijkl·uk,li=ρu··j]]>

其中σ为应力张量，ρ是密度，u是质点位移向量，u上的双圆点代表时间的不同，ε是应变张量，与质点位移有关，下标代表方向；c_ijkl是弹性系数，在之后将其收缩在C_IJ系数中；

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

(u1,u2,u3)=(U,V,W)·ejξ(x1+αx3)·e-jω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定义在一个半无限域基体涂层结构材料中，得到的弹性系数为：

CIJ(1)=λ(1)+2μ(1)λ(1)λ(1)000λ(1)+2μ(1)λ(1)000λ(1)+2μ(1)000μ(1)00Sym.μ(1)0μ(1)]]>

CIJ(2)=λ(2)+2μ(2)λ(2)λ(2)000λ(2)+2μ(2)λ(2)000λ(2)+2μ(2)000μ(2)00Sym.μ(2)0μ(2)]]>

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

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

α1,α3=±(ω/CLξ)2-1α2,α4=±(ω/CTξ)2-1]]>

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

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

uk(1)=Ukq(1)·Aq(1)·ejξ(x1+αq(1)x3)·e-jωtq=1,2,3,4]]>

uk(2)=Ukq(2)·Aq(2)·ejξ(x1+αq(2)x3)·e-jωtq=1,2]]>

σij(1)=cijkl(1)·uk,l(1)]]>

σ31(1)σ32(1)σ33(1)=D1q(1)D2q(1)D3q(1)·Aq(1)·ejξ(x1+αq(1)x3)·e-jωtq=1,2,3,4]]>

σij(2)=cijkl(2)·uk,l(2)]]>

σ31(2)σ32(2)σ33(2)=D1q(2)D2q(2)D3q(2)·Aq(2)·ejξ(x1+αq(2)x3)·e-jωtq=1,2]]>

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

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

σ31(1)|x3=d=0]]>

σ33(1)|x3=d=0]]>

u1(1)|x3=0=u1(2)|x3=0]]>

u3(1)|x3=0=u3(2)|x3=0]]>

σ31(1)|x2=0=σ31(2)|x2=0]]>

σ33(1)|x2=0=σ33(1)|x2=0]]>

并最终可得：

Ukq(1)=11110000α1(1)1-α2(1)α3(1)1-α4(1)←Uq(1)←Vq(1)←Wq(1)]]>

此时：

根据：

可得：

W1=ρc2-C11-C55α12(C13+C55)α1=(α12+1)ρCL2-(λ+2μ)-μα12(λ+μ)α1=(α12+1)(λ+2μ)-(λ+2μ)-μα12(λ+μ)α1=α1]]>

W2=ρc2-C11-C55α22(C13+C55)α2=(α22+1)ρCL2-(λ+2μ)-μα22(λ+μ)α2=(α22+1)μ-(λ+2μ)-μα22(λ+μ)α2=-1α2]]>

W₃＝α₃＝-α₁

W4=-1α4=1α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分别如下：

σ31|x3=dσ31(1)|x3=du1(1)|x3=0-u1(2)|x3=0u3(1)|x3=0=u3(2)|x3=0σ31(1)|x3=0=σ31(2)|x3=0σ33(1)|x3=0=σ33(1)|x3=0=[M]6×6A1A2A3A4A5A6={0}]]>

K(f,c;CL,CT,ρ,h)=log10(|[M]6×6|)Πs=Σj=1Ns[K(fjs,cjs;CLg,CTg,ρg,hg)]]]>

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

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

t1=2(R-Δz)vw]]>

t2=2(R-ΔzcosθSAW)vw+2·Δz·tanθSAWvSAW]]>

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

Δt=t2-t1=2(1-cosθSAW)vw·Δz]]>

即：

cosθSAW=1-vw·Δt2·Δz]]>

将Snell定律：

或

代入后，可得：

vwvSAW=1-(1-vw2·ΔtΔz)2]]>

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

vSAW=vw·[1-(1-vw2·f·Δz)2]-1/2]]>

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

2.根据权利要求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ⁿ)；

③令若停止迭代输出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，转②。