[发明专利]一种基于单纯形法的各向同性薄板材料弹性性质获取方法

摘要

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

主权项

一种基于单纯形法的各向同性薄板材料弹性性质获取方法，其特征在于：该步骤的实现过程如下，步骤1确立仿真目标函数通过对各项同性材料试件内部部分波线性组合，将其定义在一个关于边界条件的特征方程中，从而得到频散关系；设一块厚度为2h的板将其置于笛卡尔坐标系下，x1‑x3平面是Lamb波传播的矢面，而坐标x1代表波的传播方向；全局坐标（x1,x2,x3）为理论服务；步骤1.1根据运动方程：及本构方程：σ_ij＝c_ijklε_kl；几何方程：推知控制波的传播方程$\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）、水槽与水（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：散焦测量；将换能器朝试样方向移动Δ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；步骤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：模态追踪；对2.5‑22.5MHz范围内的峰值进行追踪，即可找出该频率段连续的Δz值；步骤8：波速提取；将水中的超声波波速v_W＝1500m/s，每一个峰值对应的频率与Δz带入公式即可得到该频率段内连续的表面波波速；步骤9：单纯形法获得横、纵波波速；通过单纯形法使目标函数残差绝对值最接近零，单纯形法改变C_L,C_T,ρ,h使目标函数值的达到所需范围；单纯形搜索法通过构造单纯形来逼近极小点，每构造一个单纯形，确定其最高点和最低点，然后通过扩展或压缩、反射构造新的单纯形，目地是使得极小点能包含于单纯形内。