本發(fā)明涉及一種波動類型動態(tài)數(shù)據(jù)重構(gòu)方法，具體涉及一種奇異邊界法的波動類型動態(tài)數(shù)據(jù)重構(gòu)方法。
背景技術(shù)：
：波動現(xiàn)象作為自然界最常見的自然現(xiàn)象之一，在工程應(yīng)用領(lǐng)域具有廣泛應(yīng)用與影響，如聲波減噪、海浪消能、地震波預(yù)測等科研工程領(lǐng)域，均需處理大量波動類型動態(tài)數(shù)據(jù)。由于波動現(xiàn)象基本控制方程的獨特性，如波傳播的滯后性、三維無后效性等，使波動方程基本解具有與擴散方程、拉普拉斯方程基本解有明顯的不同形式，這使得波動類型動態(tài)數(shù)據(jù)的重構(gòu)極為困難。傳統(tǒng)上處理波動類型動態(tài)數(shù)據(jù)重構(gòu)，一般有有限元方法、邊界元方法、有限差分法等。但作為傳統(tǒng)區(qū)域型網(wǎng)格方法，有限元在處理如聲波散射，聲波輻射等一系列無限域問題的數(shù)據(jù)重構(gòu)時，往往會出現(xiàn)網(wǎng)格劃分困難、計算速度慢、計算荷載大等難以克服的問題。而邊界元方法作為邊界型網(wǎng)格方法，雖然在一定程度上克服了有限元方法需要劃分區(qū)域網(wǎng)格的缺點，提高了計算效率，但由于在計算過程中需要處理大量奇異與近奇異積分，嚴(yán)重影響了數(shù)據(jù)的重構(gòu)速度?；窘夥椒ㄗ鳛橐环N無網(wǎng)格方法一定程度上克服了傳統(tǒng)網(wǎng)格類算法需要劃分網(wǎng)格的缺陷，但由于在基本解方法中為克服奇異性引入了虛假邊界，使得該方法在處理一些復(fù)雜工況條件下的數(shù)據(jù)重構(gòu)問題時，常出現(xiàn)數(shù)值不穩(wěn)定的現(xiàn)象(見文獻1.YoungDL,GuMH,FanCM.Thetime-marchingmethodoffundamentalsolutionsforwaveequations[J].EngineeringAnalysiswithBoundaryElements,2009,33(12):1411-1425.)。技術(shù)實現(xiàn)要素：發(fā)明目的：本發(fā)明的目的在于針對現(xiàn)有技術(shù)的波動類型動態(tài)數(shù)據(jù)重構(gòu)技術(shù)過程復(fù)雜、計算荷載大計算時間長的缺陷，提供一種簡單、高效、無積分、無網(wǎng)格的基于奇異邊界法的波動類型動態(tài)數(shù)據(jù)重構(gòu)方法，從而達到節(jié)約重構(gòu)時間，提高重構(gòu)效率的目的。技術(shù)方案：本發(fā)明提供了一種基于奇異邊界法的波動類型動態(tài)數(shù)據(jù)重構(gòu)方法，包括以下步驟：(1)在所考察物質(zhì)的內(nèi)部和邊界配置若干測試點，獲得這些測試點的初始數(shù)據(jù)和邊界數(shù)據(jù)；(2)直接采用三維波方程的時間基本解作為核函數(shù)，建立波傳播問題相應(yīng)的插值矩陣；(3)利用經(jīng)驗公式計算處理波動數(shù)據(jù)重構(gòu)的源點強度因子；(4)將源點強度因子代入插值矩陣，計算插值矩陣的未知系數(shù)；(5)計算任意時刻任意內(nèi)點的波動數(shù)據(jù)值。進一步，步驟(1)中，初始時刻在所考察物質(zhì)的內(nèi)部配置N1個測試點，獲得初始時刻這些測試點的波動數(shù)據(jù)值ui,i＝1,...,N1；在所考察物質(zhì)的表面配置N2個測試點，隨著波動的進行，獲得不同時刻這些點的波動數(shù)據(jù)值ui,i＝N1+1,...,N1+N2。進一步，步驟(2)中采用的波動問題控制方程為：u=Δu-1c2∂2u∂t2=0,(x,y,z)∈Ω,t>0u|Γ=u‾u|t=0=u0,∂u∂t|t=0=υ1---(1)]]>其中，Ω表示所考察物質(zhì)的區(qū)域，(x,y,z)為空間坐標(biāo)，t為(x,y,z)所對應(yīng)的時刻，c代表波速，u代表勢函數(shù)，代表邊界條件，u0和υ1表示初值條件，Δ表示Laplace算子，在聲波傳播中表示聲壓值；三維波傳播問題對應(yīng)基本解為：G(r,t)=14πrδ[(t-τ)-r/c]---(2)]]>其中δ表示狄拉克函數(shù)，c表示波速，t和τ表示配點和源點的時刻，r表示距離；由于波的傳播需要時間，基本解G僅僅用延遲時刻的相應(yīng)聲壓值的相關(guān)組合便可表示當(dāng)前時刻的聲壓值，延遲時刻取決于源點與待測點之間的距離和波的傳播速度；因此，對于特定測試點，只要求得關(guān)于源點{sj}的延遲時刻的未知系數(shù){αj}，則測試點所求時刻的聲壓值即可求出；原方程u可以被看作一個初邊值問題，應(yīng)用疊加原理，將其拆分為u1和u2，如式(3)和式(4)所示，即其中u1可以被看作一個邊值問題，u2可以被看作一個初值問題：u1*=Δu1*-1c2∂2u1*∂t2=0,(x,y,z)∈Ω,t>0u1*|Γ=u‾-u‾2*=u‾1*u1*|t=0=0,∂u1*∂t|t=0=0---(3)]]>u2*=Δu2*-1c2∂2u2*∂t2=0,(x,y,z)∈Ω,t>0u2*|Γ=u‾2*u2*|t=0=0,∂u2*∂t|t=0=0,(x,y,z)∉Ωu2*|t=0=u0,∂u2*∂t|t=0=υ1,(x,y,z)∈Ω---(4)]]>式(4)用三維泊松方程直接求出：u2*(xi)=14π∂∂r∫∫SctMu0rds+14cπ∫∫SctMυ1rds,xi∈Ω---(5)]]>其中表示半徑為ct以M為圓心的球面；在式(4)中，僅需要求出測試點M影響區(qū)域中的聲壓值，所以在影響域中布置初始數(shù)據(jù)采集點；其后，考慮方程(3)由惠更斯原理“波傳播中所到達的每一點都可以作為一個新的次波源，所有這些次波所形成的包絡(luò)面構(gòu)成下一時刻的新波面”，假定在邊界存在一系列隨時間變化的點波源這些點波源所發(fā)出的次波的總和形成計算域Ω中其后任意時刻的聲壓值其中表示聲壓強度，xm表示邊界點，則計算域Ω中的聲壓可以被表示為u1*(xi)=∫∫∫Ωα(r)4πrδ(t-rc)dv,xi∈Ω---(6)]]>在邊界Γ上布置數(shù)據(jù)采集點，以Δt為時間間隔，v表示空間積分變量，tn∈((n-1)Δt,nΔt)，t表示計算時刻，t的下標(biāo)n表示時間層數(shù)，聲壓可以被基本解G的一系列線性組合來近似：u1*(xi)=u(xi)-u2*(xi)=Σj=1,j≠iNαj(tR)Gij+αi(tR)Qi,tR=tn-rc>0,i=1,2,3......---(7)]]>其中，αj(tR)表示對于源點sj的在延遲時刻tR的未知系數(shù)，Gij表示在延遲時刻tR的三維波方程基本解，Qi表示在延遲時刻tR的源點強度因子，xi表示第i個配點，xj表示第j個源點，N為數(shù)據(jù)采集點總數(shù)；在簡諧波動中，假定αj(tR)=αj*(tR)e-iωtR≈αj*(tmΔt)e-iωtR---(8)]]>其中，0≤tmΔt-tR＜Δt，αj(tR)表示延遲時刻未知系數(shù)，表示分離時間變量t后的延遲時刻未知系數(shù)，m和n表示時間層數(shù)，如果m＜n，則未知系數(shù)αj(tR)已經(jīng)被前步求出，如果tR＜0，則表示波未傳至配點處，ω表示波的頻率，k＝ω/c表示波數(shù)；式(5)、式(7)和式(8)共同構(gòu)成了奇異邊界法處理波傳播問題的差值矩陣。進一步，步驟(3)中，計算源點強度因子的經(jīng)驗公式為：Qi=14π[π425Aj+(lnπ)2S]+ik4π---(9)]]>其中Aj表示源點sj的影響區(qū)域，S表示整個計算區(qū)域的表面積。進一步，根據(jù)步驟(3)所推得源點強度因子，未知系數(shù)由式(7)和(8)求解得出：通過式(7)和(8)，求出在第n時間層，nΔt秒的未知系數(shù)進而由式(8)求得延遲時刻的未知系數(shù)αj(tR)。進一步，步驟(5)在時刻t任意內(nèi)點x的聲壓u通過無積分計算公式求得：u2*(xi)=14π∂∂r∫∫SctMu0rds+14cπ∫∫SctMυ1rds,xi∈Ω---(10)]]>u1*(xi)=u(xi)-u2*(xi)=Σj=1,j≠iNαj(tR)Gij+αi(tR)Qi,tR=tn-rc>0,i=1,2,3,......---(11)]]>u=u1*+u2*---(12)]]>有益效果：本發(fā)明僅需邊界配置邊界數(shù)據(jù)采集點，部分區(qū)域配置初始數(shù)據(jù)采集點，無數(shù)值積分、奇異積分，無網(wǎng)格劃分，直接應(yīng)用波動方程時間基本解作為核函數(shù)，無復(fù)雜的數(shù)學(xué)變換，其簡單高效的特點切合波動類型動態(tài)數(shù)據(jù)重構(gòu)需要數(shù)據(jù)重構(gòu)技術(shù)快速、穩(wěn)定、精確的特征要求；實驗對比表明，應(yīng)用本發(fā)明所提出的技術(shù)，處理波動類型動態(tài)數(shù)據(jù)重構(gòu)問題，在取得相似精度條件下，一般耗時只需傳統(tǒng)線性邊界元算法的10％左右，具有精度高，計算快，數(shù)學(xué)簡單，程序簡便的特點，可應(yīng)用于聲波減噪，海浪消能，地震預(yù)測等波動型數(shù)據(jù)重構(gòu)，圖像處理等工程
技術(shù)領(lǐng)域：
。附圖說明圖1為本發(fā)明波動類型動態(tài)數(shù)據(jù)重構(gòu)方法流程圖；圖2(a)為三維波傳播問題邊界數(shù)據(jù)采集點布置示意圖，圖2(b)為區(qū)域數(shù)據(jù)采集點布置示意圖；圖3為邊界數(shù)據(jù)采集點sj影響區(qū)域示意圖；圖4為輪胎形區(qū)域數(shù)據(jù)采集點布置示意圖；圖5為奇異邊界法數(shù)據(jù)重構(gòu)結(jié)果隨邊界數(shù)據(jù)采集點點數(shù)增加變化收斂圖；圖6為奇異邊界法數(shù)據(jù)重構(gòu)結(jié)果隨時間延續(xù)變化趨勢圖；圖7為聲輻射無量綱實部聲壓重構(gòu)數(shù)值結(jié)果隨波數(shù)變化圖；圖8為聲輻射無量綱虛部聲壓重構(gòu)數(shù)值結(jié)果隨波數(shù)變化圖。具體實施方式下面對本發(fā)明技術(shù)方案進行詳細說明，但是本發(fā)明的保護范圍不局限于所述實施例。實施例：如圖1所示，一種基于無網(wǎng)格奇異邊界法重構(gòu)波動類型動態(tài)數(shù)據(jù)，具體步驟如下：(1)初始時刻在所考察物質(zhì)的內(nèi)部配置N1個測試點，獲得初始時刻這些測試點的波動數(shù)據(jù)值ui,i＝1,...,N1；在所考察物質(zhì)的表面配置N2個測試點，隨著波動的進行，獲得不同時刻這些點的波動數(shù)據(jù)值ui,i＝N1+1,...,N1+N2。(2)根據(jù)奇異邊界法的基本思想，直接采用三維波方程的時間基本解作為核函數(shù)，建立波傳播問題相應(yīng)的奇異邊界法插值矩陣：應(yīng)用疊加原理，原數(shù)據(jù)u可分為初始數(shù)據(jù)u2和邊界數(shù)據(jù)u1，即初始數(shù)據(jù)可以用三維泊松方程直接求出u2*(xi)=14π∂∂r∫∫SctMu0rds+14cπ∫∫SctMυ1rds,xi∈Ω---(S.1)]]>其中表示半徑為ct以M為圓心的球面。在式(S.1)中，奇異邊界法僅需要求出測試點M影響區(qū)域中的聲壓值，所以奇異邊界法僅需要在影響域中布置初始數(shù)據(jù)采集點，如圖2所示。其后，考慮邊界數(shù)據(jù)奇異邊界法在邊界Γ上布置邊界數(shù)據(jù)采集點，如圖2所示，以Δt為時間間隔，t的下標(biāo)表示時間層數(shù)。聲壓tn∈((n-1)Δt,nΔt)，可以被基本解G的一系列線性組合來近似，如式(S.2)所示：u1*(xi)=u(xi)-u2*(xi)=Σj=1,j≠iNαj(tR)Gij+αi(tR)Qi,tR=tn-rc>0,i=1,2,3......---(S.2)]]>其中αj(tR)表示對于源點sj的在延遲時刻tR的未知系數(shù)，Gij表示在延遲時刻tR的三維波方程基本解，Qi表示在延遲時刻tR的源點強度因子。在簡諧波動中，我們假定αj(tR)=αj*(tR)e-iωtR≈αj*(tmΔt)e-iωtR---(S.3)]]>其中0≤tmΔt-tR＜Δt，如果m＜n，則未知系數(shù)αj(tR)已經(jīng)被前步求出，如果tR＜0，則表示波未傳至配點處，ω表示波的頻率，k＝ω/c表示波數(shù)。式(S.1)，式(S.2)和式(S.3)共同構(gòu)成了奇異邊界法處理波傳播問題的插值矩陣。(3)利用經(jīng)驗公式計算奇異邊界法的源點強度因子：三維波動方程的源點強度因子計算公式如式(S.4)所示。Qi=14π[π425Aj+(lnπ)2S]+ik4π---(S.4)]]>其中Aj表示源點sj的影響區(qū)域，如圖3所示，S表示整個計算區(qū)域的表面積。(4)將源點強度因子代入奇異邊界插值矩陣，計算奇異邊界法求解方程的未知系數(shù)。應(yīng)用步驟(3)所推得源點強度因子，步驟(4)中奇異邊界法的未知系數(shù)可由式(S.2)和(S.3)求解得出。通過式(S.2)和(S.3)，求出在第n時間層，nΔt秒的未知系數(shù)進而由式(S.3)求得延遲時刻tR的未知系數(shù)αj(tR)。(5)根據(jù)奇異邊界法公式，計算任意時刻任意內(nèi)點的波動數(shù)據(jù)值。步驟(5)中，通過步驟(4)求得的未知系數(shù)，在時刻t任意內(nèi)點x的聲壓u可通過下面奇異邊界法公式求得：u2*(xi)=14π∂∂r∫∫SctMu0rds+14cπ∫∫SctMυ1rds,xi∈Ω---(S.5)]]>u1*(xi)=u(xi)-u2*(xi)=Σj=1,j≠iNαj(tR)Gij+αi(tR)Qi,tR=tn-rc>0,i=1,2,3......(S.6)]]>u=u1*+u2*---(S.7)]]>實施例1：考慮如圖4所示的輪胎形區(qū)域波動數(shù)據(jù)重構(gòu)，區(qū)域方程為其中R＝0.8，r＝0.2，邊界數(shù)據(jù)采集點如圖3所示。波動數(shù)據(jù)控制方程及精確解為：utt=c2(uxx+uyy+uzz),(x,y,z)∈Ωu|t=0=0ut|t=0=ck(cos(kx)+cos(ky)+cos(kz))ut|t=0=(cos(kr)+cos(ky)+cos(kz))sin(ckt),(x,y,z)∈Γ]]>本例中，我們應(yīng)用時間步長Δt＝2×10-1，波速c＝10，對每個測試點布置區(qū)域數(shù)據(jù)采集點Nf＝1255，測試點被布置在半徑為0.8，t＝1s的時間層上，圖5給出了奇異邊界法隨邊界采集數(shù)據(jù)點點數(shù)增加的數(shù)值精度收斂圖，其中波數(shù)分別為(k1＝0.5，k2＝5，k3＝10)?？梢园l(fā)現(xiàn)奇異邊界法以2階收斂速度快速收斂，當(dāng)k＝5時，收斂速度甚至達到了C＝5.0。其后，應(yīng)用時間步長Δt＝2×10-1，波速c＝10，布置邊界數(shù)據(jù)采集點Ns＝972，對每個測試點布置區(qū)域數(shù)據(jù)采集點Nf＝1255，圖6給出了在不同波數(shù)情況下(k1＝0.5，k2＝5，k3＝10)奇異邊界法數(shù)據(jù)重構(gòu)最大絕對誤差隨時間變化情況，可以看到隨時間延續(xù)，奇異邊界法最大相對誤差仍與精確解保持一致，并未隨時間變化出現(xiàn)離散情況。實施例2：考慮在單位球聲壓輻射數(shù)據(jù)重構(gòu)，其中半徑a＝1，波速c＝v0，精確解可以被表示為：u(r,θ,t)=ar(ikaz0ika-1)v0eik(r-a)+iwt,]]>其中z0＝ρ0c，表示介質(zhì)阻抗，ρ0為介質(zhì)密度，c為波速，ω＝kc為波頻。圖7和圖8給出了無量綱聲壓實部解和虛部解在邊界數(shù)據(jù)采集點Ns＝400，時間步長Δt＝0.5s條件下的數(shù)值結(jié)重構(gòu)結(jié)果?？梢园l(fā)現(xiàn)奇異邊界法重構(gòu)結(jié)果和精確解精確擬合，相比之下當(dāng)虛邊界d＝0.5時，基本解方法重構(gòu)結(jié)果和精確解擬合良好，但是當(dāng)虛邊界取為d＝0.9時，數(shù)據(jù)重構(gòu)結(jié)果發(fā)生了明顯離散。同時，需要指明的是，對于聲輻射，聲散射等外域問題的波動數(shù)據(jù)重構(gòu)，本發(fā)明方法無需采集區(qū)域波動數(shù)據(jù)，大大提高了波動類型數(shù)據(jù)的重構(gòu)效率。綜上，本發(fā)明波動類型動態(tài)數(shù)據(jù)重構(gòu)方法以奇異邊界法為基礎(chǔ)，直接應(yīng)用波動方程時間基本解作為核函數(shù)，用源點強度因子代替源點奇異項，避免了數(shù)值積分和網(wǎng)格，僅需邊界布點，無網(wǎng)格、無數(shù)值積分和數(shù)學(xué)變換，提高了數(shù)據(jù)重構(gòu)效率。相較于現(xiàn)有技術(shù)，無需做復(fù)雜的快速傅里葉變換或繁瑣的多級波形迭代求解，直接在時域上運用時間依賴基本解作為插值基函數(shù)，通過應(yīng)用源點強度因子重構(gòu)波動類型數(shù)據(jù)，其簡單高效的特點切合波動類型動態(tài)數(shù)據(jù)重構(gòu)需要計算工具快速、穩(wěn)定、精確的特征要求，適用于聲波，水波，地震波等標(biāo)量波傳播動態(tài)數(shù)據(jù)的重構(gòu)問題。實驗對比表明，應(yīng)用本發(fā)明方法處理波動類型動態(tài)數(shù)據(jù)重構(gòu)問題，在取得相似精度條件下，一般耗時只需傳統(tǒng)線性邊界元算法的10％左右，具有精度高，計算快，數(shù)學(xué)簡單，程序簡便的特點。該技術(shù)為波動型動態(tài)數(shù)據(jù)的重構(gòu)提供了新的，簡單高效的技術(shù)路線。當(dāng)前第1頁1 2 3