朱閃閃 楊 濤
(西南科技大學(xué)信息工程學(xué)院特殊環(huán)境機(jī)器人技術(shù)四川省重點(diǎn)實(shí)驗(yàn)室 綿陽(yáng) 621010)
?
基于聲線跟蹤法的空氣中超聲波聲場(chǎng)仿真*
朱閃閃楊濤
(西南科技大學(xué)信息工程學(xué)院特殊環(huán)境機(jī)器人技術(shù)四川省重點(diǎn)實(shí)驗(yàn)室綿陽(yáng)621010)
摘要針對(duì)空氣中超聲檢測(cè)存在人力計(jì)算復(fù)雜、算法驗(yàn)證困難的問(wèn)題,提出一種基于聲線跟蹤法的空氣中快速超聲檢測(cè)仿真平臺(tái)。利用CAD軟件對(duì)目標(biāo)對(duì)象三維建模,以STL文件格式導(dǎo)入OpenGL環(huán)境中,完成了STL模型的讀取和重繪。遍歷搜索目標(biāo)對(duì)象反射面上的反射點(diǎn),采用聲線跟蹤法根據(jù)超聲波在空氣中的傳播模型和斯奈爾定律確定反射聲線,接收陣元上不同聲線的聲壓疊加得到超聲波接收陣元在發(fā)射信號(hào)驅(qū)動(dòng)下的回波響應(yīng)信號(hào)。仿真結(jié)果與理論分析相吻合,為空氣中超聲檢測(cè)相關(guān)算法的研究提供了基礎(chǔ)。
關(guān)鍵詞聲線跟蹤法; 三維建模; 聲壓疊加
Simulation of Ultrasonic Sound Field in the Air Based on Ray Tracing Method
ZHU ShanshanYANG Tao
(Key Lab of Robot Technology Used for Special Environmental of Sichuan Province,School of Information Engineering, Southwest University of Science and Technology, Mianyang621010)
AbstractSince human computation is complex and algorithm verification is difficult for ultrasonic testing in the air, a rapid ultrasonic testing simulation method in the air based on tracing method is proposed. The proposed method utilizes CAD software for three-dimensional modeling of the target object, imports the OpenGL environment under the Microsoft Visual Studio C++ by STL file format, completes the reading and redrawing of STL models. Then, it traverses the reflection points of the target object reflective surfaces, adopts ray tracing law to identify the reflective sound ray according to ultrasonic wave propagation model in the air and Snell’s law, superposes the pressure of sound ray, obtains the echo response signal of ultrasonic receiving array driving by the transmitting signal from different sound ray pressure of superposition on receiving array. Simulation results coincide well with theoretical analysis, and provide the foundations to research ultrasonic testing in the air and related algorithms.
Key Wordsray tracing method, three-dimensional modeling, sound pressure of superposition
Class NumberTG115.28
1引言
超聲檢測(cè)因成本低、方向性好、穿透力強(qiáng)、不受電磁輻射影響等優(yōu)點(diǎn)廣泛用于工業(yè)生產(chǎn)和日常生活中。特別是在大霧、雨雪、強(qiáng)輻射等惡劣環(huán)境中,超聲波更表現(xiàn)出其獨(dú)有的優(yōu)勢(shì)。
隨著計(jì)算機(jī)技術(shù)的發(fā)展,超聲檢測(cè)仿真和可視化的研究也逐步深入。目前,超聲波仿真和可視化方面的研究主要是利用二維分布圖形對(duì)換能器發(fā)射聲場(chǎng)的計(jì)算仿真,對(duì)聲場(chǎng)三維仿真及可視化的研究相對(duì)比較少。通過(guò)解析法[1]和數(shù)值法[2]等數(shù)學(xué)建模方法建立超聲檢測(cè)系統(tǒng)的模型。數(shù)值法主要包括有限元(FEM)[3]、邊界元(BEM)[4]、時(shí)域有限差分(FDM)[5]等為代表的波動(dòng)聲學(xué)方法,求解精
度很高、計(jì)算規(guī)模大、效率低,更適合于小范圍區(qū)域的求解。解析法以聲線跟蹤和虛聲源幾何聲學(xué)方法為代表[6],包括聲線模型、高斯聲束模型和衍射理論模型等,采用簡(jiǎn)化物理模型,計(jì)算效率高,但是求解精度會(huì)降低,不適用于較復(fù)雜的場(chǎng)景。文獻(xiàn)[7]基于瑞利積分和Pencil建立半解析法計(jì)算模型,實(shí)現(xiàn)超聲檢測(cè)發(fā)射聲場(chǎng)和缺陷回波聲場(chǎng)的快速計(jì)算。聲線跟蹤法通用性好,基于聲線跟蹤方法的超聲波建模技術(shù),主要用于傳播路徑的計(jì)算,在回波信號(hào)的定量計(jì)算方面缺少深入的研究。文獻(xiàn)[8]基于聲線模型、超聲傳播原理和幾何光學(xué)簡(jiǎn)化實(shí)現(xiàn)了超聲傳播路徑的模擬。文獻(xiàn)[9]以聲線跟蹤法為基礎(chǔ),運(yùn)用Visual C++和OpenGL對(duì)超聲波傳播路徑進(jìn)行了模擬。
本文用3D參數(shù)化CAD軟件PTC Creo Parametric 3.0對(duì)目標(biāo)對(duì)象三維建模,以STL文件格式導(dǎo)入Microsoft Visual Studio C++下的OpenGL環(huán)境中,完成了STL模型的讀取和重繪,遍歷搜索目標(biāo)對(duì)象反射面上的反射點(diǎn),求取發(fā)射線的方向向量,采用聲線跟蹤法根據(jù)超聲波在空氣中的傳播模型和斯奈爾定律確定反射聲線,考慮接收陣元接收區(qū)域的限制及反射點(diǎn)和發(fā)射和接收陣元的遮擋問(wèn)題,確定接收點(diǎn),進(jìn)而完成超聲波發(fā)射和接收路徑的繪制。接收陣元上不同聲線的聲壓疊加得到超聲波接收陣元在發(fā)射信號(hào)驅(qū)動(dòng)下的回波響應(yīng)信號(hào)。
2超聲傳播路徑仿真的實(shí)現(xiàn)
2.1對(duì)象目標(biāo)建模及導(dǎo)入
STL文件是若干空間小三角形面片的集合,通過(guò)三維模型表面的三角網(wǎng)格化獲得,每個(gè)三角形面片用三角形的三個(gè)頂點(diǎn)和指向模型外部的三角面片的法向量組成[10]。STL文件有二進(jìn)制和ASCII兩種格式,前者更緊湊,后者更容易理解。
本文將采用PTC Creo Parametric 3.0構(gòu)建直徑10cm、高15cm的圓柱體模型,并保存為STL文件格式,該文件可對(duì)目標(biāo)對(duì)象表面自動(dòng)進(jìn)行三角剖分。根據(jù)步長(zhǎng)的不同保存為不同的STL文件格式,步長(zhǎng)越小,剖分越精細(xì)。利用C++標(biāo)準(zhǔn)庫(kù)中提供的文件流讀取及字符操作等功能結(jié)合OpenGL中三角面片繪制的編程技術(shù)實(shí)現(xiàn)對(duì)STL文件格式的讀取和顯示,將目標(biāo)對(duì)象的由局部坐標(biāo)系變換到三維場(chǎng)景的世界坐標(biāo)系中。STL文件讀取和顯示的流程如圖1所示。
2.2超聲波路徑傳播仿真
2.2.1聲線跟蹤法及建模數(shù)學(xué)基礎(chǔ)
聲線法認(rèn)為聲在空間中以聲線的方式向四周傳播,聲線碰到目標(biāo)對(duì)象后,在碰撞點(diǎn)處作鏡面反射,沿著新的方向繼續(xù)前進(jìn),判定是否被接收陣元接收。超聲波傳播路徑示意圖如圖2所示。
圖1 STL文件讀取和顯示流程圖
圖2 超聲波傳播路徑示意圖
利用OpenGL函數(shù),繪制兩個(gè)圓柱體,頂部的圓作為探頭的發(fā)射和接收陣元,左邊白色圓代表發(fā)射陣元,右邊黑色視為接收陣元。其中,發(fā)射點(diǎn)A(x1,y1,z1)為發(fā)射陣元中心,陣元半徑為R1;C(x0,y0,z0)為三角面片的重心;n(v1,v2,v3)為三角面片的法向量;B(x2,y2,z2)接收陣元的中心,陣元半徑為R2,在構(gòu)建的三維場(chǎng)景中接收平面的法向量已知設(shè)p(p1,p2,p3),發(fā)射線的方向向量r1由式(1)可得,反射線的方向向量r2(u1,u2,u3)由式(2)可得:
r1={(x0-x1),(y0-y1),(z0-z1)}
(1)
(2)
由反射線方向向量r2和反射線上點(diǎn)C的坐標(biāo),可得到反射線的方程,其中N為參數(shù)。
(3)
將平面方程寫成點(diǎn)法式方程形式,即為
p1(x-x2)+p2(y-y2)+p3(z-z2)=0
(4)
將式(3)和式(4)聯(lián)立可得
(5)
令D=p1u1+p2u2+p3u3,如果D為0,則沒(méi)有交點(diǎn);如果不為0,求出N,帶入式(3)可求得與陣元平面的交點(diǎn)O(x3,y3,z3)。交點(diǎn)O到接收陣元中心的距離為
(6)
判斷d是否小于R2,如果小于R2,并且此時(shí)三角面片的重心C(x0,y0,z0)與發(fā)射和接收陣元在同一側(cè),則此交點(diǎn)存在,反射線和反射點(diǎn)符合要求,超聲波的發(fā)射和反射的傳播路徑確定。
2.2.2超聲波傳播路徑確定流程
超聲波的傳播路徑在聲線跟蹤法基礎(chǔ)上,根據(jù)斯奈爾定律確定反射聲線。通過(guò)構(gòu)建三維場(chǎng)景和坐標(biāo)系變換,將目標(biāo)對(duì)象坐標(biāo)和發(fā)射接收陣元坐標(biāo)變換到世界坐標(biāo)系中。遍歷三角面片,得到各面片重心坐標(biāo),并采用式(2)計(jì)算出反射線方向向量,由反射線和接收陣元所在平面構(gòu)建方程組,求取接收點(diǎn)。考慮到實(shí)際情況,接收點(diǎn)應(yīng)在陣元之內(nèi)且反射點(diǎn)必須和發(fā)射和接收陣元在同側(cè),進(jìn)而確定超聲波發(fā)射和接收的傳播路徑。流程圖如圖3所示。
圖3 傳播路徑確定流程圖
3聲場(chǎng)計(jì)算分析
3.1超聲波聲場(chǎng)計(jì)算
圓形平面活塞鑲嵌在一個(gè)無(wú)限大的障板上,當(dāng)活塞以速度u=uaejωt振動(dòng)時(shí),其輻射聲場(chǎng)中的聲壓分布為[11~12]
(7)
其中a為超聲探頭的半徑;ua為振速幅值;ρ0為媒質(zhì)的靜態(tài)密度;ω為諧波頻率;k為波數(shù);J1為一階第一類貝塞爾函數(shù)。
從式(7)可以推導(dǎo)出圓形活塞式超聲波換能器的指向性。根據(jù)貝塞爾函數(shù)性質(zhì),當(dāng)x=0時(shí),J1(x)/x=1/2,超聲波換能器的聲場(chǎng)指向性可表述為
(8)
假設(shè)發(fā)射聲場(chǎng)初始聲壓為p0(t),標(biāo)定聲壓為pbd(t)。
pbd(t)=D(θ)·e-αrbdp0(t)
(9)
(10)
(11)
式中,t0為標(biāo)定信號(hào)的延時(shí),r0為標(biāo)定信號(hào)發(fā)射陣元和接收陣元的間距,c0為空氣中超聲波的速度。
在t時(shí)刻,第i條聲線的延遲時(shí)間見(jiàn)式(12):
(12)
在t時(shí)刻,第i條聲線的聲壓見(jiàn)式(13):
(13)
其中ri(t)為發(fā)射點(diǎn)到反射點(diǎn)間的距離和反射點(diǎn)到接收點(diǎn)的距離之和,θi(t)為發(fā)射線和中心軸線的夾角,tsi為延遲時(shí)間,c0為空氣中超聲波的速度,αq為目標(biāo)對(duì)象吸收系數(shù),e-αri為介質(zhì)的吸收衰減。在t時(shí)刻,則接收陣元的上的聲壓見(jiàn)式(14):
(14)
3.2聲場(chǎng)的標(biāo)定
圖4 示波器采集回波響應(yīng)波形圖
在空氣中,對(duì)于單個(gè)發(fā)射和接收探頭采集數(shù)據(jù)信息。其中,f取40kHz,ρ0取1.205kg/m3,a取0.004m,角度為0,r取1.0m,Tektronix TDS2024B數(shù)字存儲(chǔ)示波器采集獲得回波響應(yīng)波形圖如圖4所示,將示波器采集的數(shù)據(jù)導(dǎo)入Matlab采樣讀取1024個(gè)數(shù)據(jù),得到波形如圖5所示。
圖5 Matlab采樣讀取波形圖
采用標(biāo)定的方法,需要求出發(fā)射陣元的驅(qū)動(dòng)信號(hào)的數(shù)學(xué)模型,根據(jù)1024個(gè)采樣數(shù)據(jù),利用Matlab工具箱擬合得到擬合函數(shù)為
pbd(t)=1.3880sin(2.5250×105t+1.0600)
+0.65420sin(2.5750×105t-0.5357)
+0.59350sin(2.4570×105t+3.6550)
+3.9570sin(2.4070×105t-2.3190)
+0.6322sin(2.662×105t-3.4430)
+3.727sin(2.405×105t-5.3260)
+0.6086sin(2.664×105t-0.4687)
+0.04323sin(3.426×105t-1.291)
(15)
4仿真結(jié)果分析
在Microsoft Visual Studio C++編程下的OpenGL環(huán)境中,對(duì)單個(gè)發(fā)射和接收陣元的超聲波探頭的傳播路徑和回波信息進(jìn)行仿真。其中,左邊圓為發(fā)射陣元,右邊圓為接收陣元。考慮到實(shí)際環(huán)境中,超聲波傳播距離的限制及推廣到超聲相控陣中陣列的排布,選取發(fā)射和接收間距為27.5mm,陣元和目標(biāo)對(duì)象之間的距離l選擇0.8m、0.6m、0.2m。
當(dāng)距離為0.8m時(shí),間距為27.5mm,發(fā)射和接收陣元坐標(biāo)分別為(60,60,1200)和(87.5,60,1200),結(jié)果如圖6(a)和圖6(b)所示。
圖6(a) l取0.8m時(shí)傳播路徑聲線圖
圖6(b) l取0.8m時(shí)接收陣元聲壓圖
當(dāng)距離為0.6m時(shí),間距為27.5mm,發(fā)射和接收陣元坐標(biāo)為:(60,60,800)和(87.5,60,800),結(jié)果如圖7(a)和圖7(b)所示。
圖7(a) l取0.6m時(shí)傳播路徑聲線圖
圖7(b) l取0.6m時(shí)接收陣元聲壓圖
當(dāng)距離為0.2m時(shí),間距為27.5mm,發(fā)射和接收陣元坐標(biāo)為(60,60,400)和(87.5,60,400),結(jié)果如圖8(a)和圖8(b)所示。
圖8(a) l取0.2m時(shí)傳播路徑聲線圖
圖8(b) l取0.2m時(shí)接收陣元聲壓圖
由圖6(a)、圖7(a)、圖8(a)和圖6(b)、圖7(b)、圖8(b)可以看出,隨著距離的減小及聲線數(shù)的增加,接收點(diǎn)的聲壓信號(hào)越準(zhǔn)確。
5結(jié)語(yǔ)
本文通過(guò)對(duì)目標(biāo)對(duì)象建立三維幾何模型,根據(jù)聲線跟蹤法和斯奈爾定律對(duì)超聲波聲場(chǎng)的發(fā)射和接收傳播路徑進(jìn)行仿真,對(duì)聲線聲壓疊加得到接收陣元的聲壓信號(hào)。本仿真平臺(tái)可擴(kuò)展到多個(gè)發(fā)射陣元、單個(gè)接收陣元,單個(gè)發(fā)射陣元、多個(gè)接收陣元以及多發(fā)射、多接收的仿真場(chǎng)景中,實(shí)現(xiàn)超聲檢測(cè)相關(guān)算法性能的快速評(píng)估。
參 考 文 獻(xiàn)
[1] Kundu T, Placko D, Rahani E K, et al. Ultrasonic Field Modeling: A Comparison Of Analytical, Semi-Analytical, And Numerical Techniques[J]. Ultrasonics Ferroelectrics & Frequency Control IEEE Transactions on,2010,57(12):2795-2807.
[2] Liu Q, Wirdelius H. A 2D model of ultrasonic wave propagation in an anisotropic weld[J]. Ndt & E International,2007,40(3):229-238.
[3] Hill R, Forsyth S A, Macey P. Finite element modelling of ultrasound, with reference to transducers and AE waves[J]. Ultrasonics,2004,42(1-9):253-8.
[4] Maruyama T, Saitoh T, Hirose S. Simulation for Air-Coupled Ultrasound Testing Using Time-Domain BEM[J]. Aip Conference Proceedings,2014,1581(1):550-555.
[5] Yamawaki H, Saito T. Numerical calculation of ultrasonic propagation with anisotropy[J]. Ndt & E International,2000,33(7):489-497.
[6] 張偉志,剛鐵,王軍.超聲波檢測(cè)計(jì)算機(jī)模擬和仿真的研究及應(yīng)用現(xiàn)狀[J].應(yīng)用聲學(xué),2003,22(3):39-44.
ZHANG Weizhi, GANG Tie, WANG Jun. A survey of the simulation and its application in ultrasonic testing[J]. Applied Acoustics,2003,22(3):39-44.
[7] 丁輝,張俊,張益成,等.超聲檢測(cè)聲場(chǎng)計(jì)算模型的建立與仿真軟件的開(kāi)發(fā)[J].無(wú)損檢測(cè),2009,(8):614-618.
DING Hui, ZHANG Jun, ZHANG Yicheng, et al. The Establishment of Ultrasonic Field Calculation Model and the Development of Simulation Software[J]. Nondestructive Testing,2009,(8):614-618.
[8] 張偉志.基于聲線模型的超聲波傳播路徑模擬[D].哈爾濱:哈爾濱工業(yè)大學(xué),2002.
ZHANG Weizhi. Ultrasonic Simulation of propagation path Based on Ray Modeling[D]. Harbin: Harbin Institute of Technology,2002.
[9] 曲文元,熊建設(shè),秦勃,等.基于聲線跟蹤法開(kāi)發(fā)超聲無(wú)損檢測(cè)動(dòng)態(tài)仿真軟件[J].計(jì)算機(jī)仿真,2006,23(4):269-273.
QU Wenyuan, XIONG Jianshe, QIN Bo, et al. Development of Simulation System for Nonderstructive Ultrasonic Inspection Based on Ray-Tracing Method[J]. Computer Simulation,2006,23(4):269-273.
[10] 嚴(yán)梽銘,鐘艷如.基于VC++和OpenGL的STL文件讀取顯示[J].計(jì)算機(jī)系統(tǒng)應(yīng)用,2009,18(3):172-175.
YAN Zhiming, ZHONG Yanru. STL Reading and Visualization Based on VC++ and OpenGL[J]. Computer System Applications,2009,18(3):172-175.
[11] 杜功煥.聲學(xué)基礎(chǔ)[M].南京:南京大學(xué)出版社,2012.
DU Gonghuan. Foundamentals of Acoustics[M]. Nanjing: Nanjing University Press,2012.
[12] 馮若.超聲手冊(cè)[M].南京:南京大學(xué)出版社,1999.
FENG Ruo. Ultrasonics Handbook[M]. Nanjing: Nanjing University Press,1999.
中圖分類號(hào)TG115.28
DOI:10.3969/j.issn.1672-9722.2016.03.012
作者簡(jiǎn)介:朱閃閃,女,碩士研究生,研究方向:超聲檢測(cè)建模與仿真。楊濤,男,博士,教授,研究方向:強(qiáng)輻射機(jī)器人技術(shù),聲學(xué)陣列信號(hào)處理。
基金項(xiàng)目:特殊環(huán)境機(jī)器人技術(shù)四川省重點(diǎn)實(shí)驗(yàn)室開(kāi)放基金(編號(hào):13zxtk06)資助。
收稿日期:2015年9月10日,修回日期:2015年10月31日