[發(fā)明專利]一種基于單純形法的各向同性薄板材料彈性性質(zhì)獲取方法

有效

• 申請(qǐng)?zhí)枺?/td>: 201410135545.3
• 申請(qǐng)日: 2014-04-04
• 公開（公告）號(hào): CN103926315A
• 公開（公告）日: 2014-07-16
• 發(fā)明（設(shè)計(jì)）人: 宋國榮;劉宏實(shí);呂炎;李子木;何存富;吳斌
• 申請(qǐng)（專利權(quán)）人: 北京工業(yè)大學(xué)
• 主分類號(hào): G01N29/04
• 分類號(hào): G01N29/04;G01N29/46
• 代理公司: 北京思海天達(dá)知識(shí)產(chǎn)權(quán)代理有限公司 11203
• 代理人: 沈波

1.一種基于單純形法的各向同性薄板材料彈性性質(zhì)獲取方法，其特征在于：該步驟的實(shí)現(xiàn)過程如下，

步驟1確立仿真目標(biāo)函數(shù)

通過對(duì)各項(xiàng)同性材料試件內(nèi)部部分波線性組合，將其定義在一個(gè)關(guān)于邊界條件的特征方程中，從而得到頻散關(guān)系；設(shè)一塊厚度為2h的板將其置于笛卡爾坐標(biāo)系下，x1-x3平面是Lamb波傳播的矢面，而坐標(biāo)x1代表波的傳播方向；全局坐標(biāo)（x1,x2,x3）為理論服務(wù)；

步驟1.1根據(jù)運(yùn)動(dòng)方程：及本構(gòu)方程：σ_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.]]>

其中σ_ij為應(yīng)力張量，ρ是密度，u_j是質(zhì)點(diǎn)位移向量，t為時(shí)間，圓點(diǎn)代表時(shí)間的不同，指針從逗號(hào)后開始代表空間坐標(biāo)的不同，重復(fù)指針的合計(jì)就如張量標(biāo)志自動(dòng)的被設(shè)定在自此之后；c_ijkl是彈性系數(shù)在之后可以將其收縮在C_IJ系數(shù)中，ε_kl是應(yīng)變張量，其與質(zhì)點(diǎn)位移有關(guān)；

根據(jù)彈性動(dòng)力學(xué)理論，現(xiàn)假設(shè)一個(gè)平面諧波沿著x₁傳播，其角頻率為ω，該波相應(yīng)的物理量在x₂方向上獨(dú)立，波的形式u_k如下：

(u1,u2,u3)=(U,V,W)·ejξ(x1+αx3)·e-jωt]]>

其中ξ為波數(shù)，而α則表示一個(gè)未知系數(shù)，(ξ,0,αξ)為該波的傳播向量，(U,V,W)為該平面波相應(yīng)的振幅；將u_k帶入控制波的傳播方程cijkl·uk,li=ρu··j]]>得到：[Kij]3×3·UVW=000,]]>該線性代數(shù)形式[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)是其相對(duì)相應(yīng)的特征向量；

在各向同性材料薄板的波場中可以利用以下的彈性系數(shù)關(guān)系；首先我們將C_IJ定義在一個(gè)各向同性材料中，得到的彈性系數(shù)為:

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

λ與μ為拉梅系數(shù)，我們定義C_IJ＝c_ijkl，其中ij→I或J，定義11→1，22→2，33→3，23或32→4，31或13→5，12或21→6；

獲得α的表示式，從而可知該波的傳播方式；通常，在各向同性材料中，α的值有四個(gè)；

α1=-α3=(cCL)2-1α2=-α4=(cCT)2-1]]>

將不同的α值下對(duì)應(yīng)的波的形式累加，得到以下位移、應(yīng)力公式來表示波場；A_q為該平面波的振幅系數(shù)

(u1,u2,u3)=Σq=14(Uq,0,Wq)·Aq·eiξ(x1+αqx3)·e-iωt]]>

(σ33,σ13,σ23)=Σq=14(D1q,D2q,D3q)·Aq·eiξ(x1+αqx3)·e-iωt]]>

獲得邊界條件的系數(shù)矩陣：

在無外力為邊界條件的自由載荷狀態(tài)下，板厚為2h，得到量化的漏表面波頻散曲線，其系數(shù)矩陣為[M]_4×4：

σ33|x3=+hσ13|x3=+hσ33|x3=-hσ13|x3=-h=[M]4×4A1A2A3A4=0000]]>

在各向同性的平面板中相應(yīng)的位移幅值比為以下值：

Uq=1W1=α1,W2=-1α2,W3=α3,W4=-1α4]]>

α_q已經(jīng)得出為:

α1=-α3=(cCL)2-1α2=-α4=(cCT)2-1]]>

在各向同性的平面板中相應(yīng)的應(yīng)變幅值比為以下值：

D1q=λ+(λ+2μWq)αqD2q=μ(Wq+αq)]]>

λ,μ是拉梅系數(shù),可表示橫波波速C_T，縱波波速C_L為以下形式:

CL2=λ+2μρCT2=μρ]]>

最后，可以得出系數(shù)矩陣[M]_4×4:

[M]4×4=D11e+iξhα1D12e+iξhα2D13e+iξhα3D14e+iξhα4D21e+iξhα1D22e+iξhα2D23e+iξhα3D24e+iξhα4D11e-iξhα1D12e-iξhα2D13e-iξhα3D14e-iξhα4D21e-jξhα1D22e-iξhα2D23e-iξhα3D24e-iξhα4]]>

由邊界條件[M]_4×4行列式可得四個(gè)方向的頻散曲線，再得到其頻散特征函數(shù)K，為了得到非零解，函數(shù)K的值應(yīng)為0，C_L,C_T,ρ,h都為已知，改變f,c的值使目標(biāo)函數(shù)值最小保證其盡量趨近于0，頻散曲線，頻散特征函數(shù)，仿真的目標(biāo)函數(shù)Π_s分別如下：

σ33|x3=+hσ13|x3=+hσ33|x3=-hσ13|x3=-h=[M]4×4A1A2A3A4=0000]]>

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

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

步驟1.2：確立波速提取的公式；

在單頻激勵(lì)/接收的情況下，圖2所示的漏表面波傳播示意圖中，上表面的直接反射回波I傳播的時(shí)間與漏表面波L的傳播時(shí)間分別為：

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

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

其中R為聚焦半徑，Δz為散焦距離，v_w為水的超聲波波速，θ_SAW為產(chǎn)生表面波的瑞利角，v_SAW為材料的表面波波速；因此兩者的時(shí)間差為：

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

即：

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

將Snell定律：

sinθSAW=vwvSAW]]>或θSAW=sin-1(vwvSAW)]]>

代入后，可得：

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

此時(shí)如果Δz恰為一個(gè)V(z)曲線的振蕩周期時(shí)，1/Δt則為換能器的激勵(lì)頻率f；如果Δz能夠確定，便可使用如下公式進(jìn)行表面波波速的計(jì)算：

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

因此，測(cè)量被測(cè)材料的V(z)曲線振蕩周期成為波速提取的重點(diǎn)；

步驟2）：搭建測(cè)試系統(tǒng)；

為了方便散焦步進(jìn)測(cè)量，搭建了一套進(jìn)行散焦步進(jìn)測(cè)量的測(cè)試系統(tǒng)；該測(cè)試系統(tǒng)主要包括試樣（1）、水槽與水（2）、換能器（3）、移動(dòng)平臺(tái)（4）、脈沖激勵(lì)/接收儀（5）、示波器（6）、GPIB總線（7）、PXI總控制系統(tǒng)（8）、移動(dòng)伺服馬達(dá)（9）、旋轉(zhuǎn)軸（10）；其中，在移動(dòng)平臺(tái)（4）下面安裝換能器（3），換能器（3）與脈沖激勵(lì)/接收儀（5）相連，脈沖激勵(lì)/接收儀（5）與示波器（6）相連，示波器（6）通過GPIB總線（7）與PXI總控制系統(tǒng)（8）相連，PXI總控制系統(tǒng)（8）與移動(dòng)伺服馬達(dá)（9）相連，同時(shí)PXI總控制系統(tǒng)（8）與旋轉(zhuǎn)軸（10）相連；

步驟3：聚焦面數(shù)據(jù)采集；

以長方體碳化鎢為被測(cè)試樣，其尺寸為40mm×40mm×10mm，將換能器3聚焦到試樣的上表面，通過脈沖激勵(lì)/接收儀（5）在發(fā)出一個(gè)帶寬為10-200MHz的脈沖后轉(zhuǎn)換為接收狀態(tài)，當(dāng)接收到反射信號(hào)后，將信號(hào)傳輸進(jìn)示波器（6），示波器的采樣頻率f_S＝2.5GHz，采樣點(diǎn)數(shù)N_s＝10000；經(jīng)過示波器的低通濾波后，通過GPIB總線（7）存儲(chǔ)進(jìn)PXI總控制系統(tǒng)，聚焦面的時(shí)域波形所示；

步驟4：散焦測(cè)量；

將換能器朝試樣方向移動(dòng)Δz₀＝10μm，待移動(dòng)完成后進(jìn)行電壓數(shù)據(jù)采集，采集結(jié)束后再將換能器朝試樣方向移動(dòng)Δz₀＝10μm進(jìn)行數(shù)據(jù)采集，采樣頻率f_S＝2.5GHz，采樣點(diǎn)數(shù)N_s＝10000，如此循環(huán)往復(fù)，共移動(dòng)4mm，因此將得到400組電壓數(shù)據(jù)，將聚焦面的電壓數(shù)據(jù)包含在內(nèi)共得到M＝401組電壓數(shù)據(jù)；

步驟5：時(shí)域傅里葉變換；

將測(cè)得的數(shù)據(jù)進(jìn)行時(shí)域傅里葉變換；

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

其中：A_i為時(shí)域傅里葉變換后的頻譜值，x_i代表一組電壓數(shù)據(jù)，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，需要對(duì)時(shí)域傅里葉變換的結(jié)果再進(jìn)行沿散焦距離方向的空間傅里葉變換，將散焦距離z變換至z^-1域：

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

其中：B_i為空間傅里葉變換后的頻譜值，A_m代表沿散焦方向的時(shí)域傅里葉變換的頻譜值，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：模態(tài)追蹤；

對(duì)2.5-22.5MHz范圍內(nèi)的峰值進(jìn)行追蹤，即可找出該頻率段連續(xù)的Δz值；

步驟8：波速提??；

將水中的超聲波波速v_W＝1500m/s，每一個(gè)峰值對(duì)應(yīng)的頻率與Δz帶入公式即可得到該頻率段內(nèi)連續(xù)的表面波波速；

步驟9：單純形法獲得橫、縱波波速；

通過單純形法使目標(biāo)函數(shù)殘差絕對(duì)值最接近零，單純形法改變C_L,C_T,ρ,h使目標(biāo)函數(shù)值的達(dá)到所需范圍；

單純形搜索法通過構(gòu)造單純形來逼近極小點(diǎn)，每構(gòu)造一個(gè)單純形，確定其最高點(diǎn)和最低點(diǎn)，然后通過擴(kuò)展或壓縮、反射構(gòu)造新的單純形，目地是使得極小點(diǎn)能包含于單純形內(nèi)。

2.根據(jù)權(quán)利要求1所述的一種基于單純形法的各向同性薄板材料彈性性質(zhì)獲取方法，其特征在于：本方法中未知量有四個(gè)（C_L,C_T,ρ,h），為四維變量問題；

用單純形搜索法求無約束問題minF(x)，x∈Rⁿ的算法步驟如下：

①選取初始單純形{x⁰,x¹,…,xⁿ}，反映系數(shù)α＞1，緊縮系數(shù)θ∈(0,1)，擴(kuò)展系數(shù)γ＞1，收縮系數(shù)β∈（0，1）及精度ε>0，置k=0；

②將單純形的n+1個(gè)頂點(diǎn)按目標(biāo)函數(shù)值的大小重新編號(hào)，使頂點(diǎn)的編號(hào)滿足F(x⁰)≤F(x¹)≤…≤F(x^n-1)≤F(xⁿ)；

③令xn+1=1nΣj=0n-1xj,]]>若{1n+1Σj=0n[F(xj)-F(xn+1)]2}2≤ϵ]]>停止迭代輸出x⁰，否則轉(zhuǎn)入④；

④計(jì)算xⁿ⁺²＝xⁿ⁺¹+α(xⁿ⁺¹-x），若F(xⁿ⁺²)＜F(x⁰)，轉(zhuǎn)⑤，否則當(dāng)F(xⁿ⁺²)＜F(x^n-1)時(shí)轉(zhuǎn)⑥，當(dāng)F(xⁿ⁺²)≥F(x^n-1)轉(zhuǎn)⑦；

⑤計(jì)算xⁿ⁺³＝xⁿ⁺¹+γ(xⁿ⁺²-xⁿ⁺¹)，若F(xⁿ⁺³)＜F(x⁰)，令xⁿ＝xⁿ⁺³，轉(zhuǎn)②，否則轉(zhuǎn)⑥；

⑥令xⁿ＝xⁿ⁺²，轉(zhuǎn)②

⑦令xⁿ＝{xⁱ|F(xⁱ)＝min(F(xⁿ),F(xⁿ⁺²))}，計(jì)算xⁿ⁺⁴＝xⁿ⁺¹+β(xⁿ-xⁿ⁺¹)，若F(xⁿ⁺⁴)＜F(xⁿ)，令xⁿ＝xⁿ⁺⁴，轉(zhuǎn)②，否則轉(zhuǎn)⑧；

⑧令x^j＝x⁰+θ(x^j-x⁰)，j＝0,1,…,n，轉(zhuǎn)②。