奇異邊界法中的兩種反插技術(shù)

-.z奇異邊界法中的兩種反插技術(shù)1)1)10672051.陳文2)2):.〔河海大學(xué)工程力學(xué)系，****210098〕摘要奇異邊界法直接使用根本解做為插值基函數(shù)，且源點和配點為同一組物理邊界上的離散點，是一種真正的無網(wǎng)格邊界離散方法。奇異邊界法的核心是通過反插值技術(shù)，計算源點強度因子，即插值矩陣的對角線元素。而非對角線元素可以使用根本解直接求得。本文考察了奇異邊界法中消除根本解原點奇異性問題的兩種反插值技術(shù)，對它們做了比較研究；以Laplace和Helmholtz方程為例檢驗了這些方法，并做了一些數(shù)值分析。關(guān)鍵詞奇異邊界法，根本解，反插值技術(shù)引言與基于網(wǎng)格劃分的有限元法等方法相比，無網(wǎng)格法由于不需要劃分網(wǎng)格，自上世紀(jì)90年代興起了廣泛研究的熱潮ADDINNE.Ref.{08E9EFFA-1903-45BE-AFB5-47BD40BFA351}[1-3]。無網(wǎng)格法可分為區(qū)域型無網(wǎng)格法和邊界型無網(wǎng)格法，邊界型無網(wǎng)格法只需要一組邊界節(jié)點來離散求解區(qū)域邊界，直接借助于邊界離散點來構(gòu)造近似函數(shù)，因此受到眾多學(xué)者的青睞。目前備受關(guān)注的邊界型無網(wǎng)格法主要包括：根本解法〔methodoffundamentalsolutions:MFS〕ADDINNE.Ref.{3A1602D1-EE10-4CB1-90E0-EEC6FCD5176A}[4-8]、邊界節(jié)點法〔boundaryknotmethod:BKM〕ADDINNE.Ref.{B080150E-0FA6-41F6-A255-31CA74811DC5}[9,10]、奇異無網(wǎng)格法〔singularmeshlessmethod:SMM〕ADDINNE.Ref.{B0E72E48-CB7D-446F-B8CC-65DD62C0D88A}[11,12]、修正根本解法〔modifiedmethodoffundamentalsolutions:MMFS〕ADDINNE.Ref.{5300B502-A419-4A35-AA10-4B2B0F603E1E}[13]等。根本解法由于數(shù)學(xué)簡單、編程容易和高精度等優(yōu)點，吸引許多國際力學(xué)和數(shù)學(xué)學(xué)者近年來的深入研究。該方法的主要缺點是為了防止控制微分方程根本解的原點奇異性，在物理邊界之外引入了虛假邊界，而虛假邊界的選取隨意性較大，在求解復(fù)雜幾何域問題時容易造成計算不穩(wěn)定，因此根本解法多用于計算常規(guī)幾何形狀問題。為了防止使用虛假邊界，Chen等ADDINNE.Ref.{B080150E-0FA6-41F6-A255-31CA74811DC5}[9,10]提出了邊界節(jié)點法，使用控制微分方程的非奇異徑向基函數(shù)一般解代替奇異根本解。該方法在計算二維和三維復(fù)雜幾何域時，其解的精度和穩(wěn)定性都很高ADDINNE.Ref.{B080150E-0FA6-41F6-A255-31CA74811DC5}[9,10]，但對于一些控制方程，如Laplace方程等，沒有非奇異一般解，從而限制了邊界節(jié)點法的應(yīng)用*圍。為了抑制以上缺點，Young等ADDINNE.Ref.{B0E72E48-CB7D-446F-B8CC-65DD62C0D88A}[11,12]基于雙層勢理論提出了一種奇異無網(wǎng)格法，利用去奇異〔desingularization〕技術(shù)計算插值矩陣中的對角元素，該方法的缺乏是計算精度不高且需要在物理邊界上等間距布點ADDINNE.Ref.{6E5B8502-63EF-4DA6-9AE9-3AAEED9F489E}[14]，難于處理復(fù)雜幾何區(qū)域問題。類似于奇異無網(wǎng)格法，BozidarADDINNE.Ref.{5300B502-A419-4A35-AA10-4B2B0F603E1E}[13]最近提出了修正根本解法，該方法解的精度高于根本解法，并且僅僅需要邊界節(jié)點信息，然而在計算插值矩陣的對角元素時仍需要積分計算?；谝陨涎芯?，ChenADDINNE.Ref.{93B4A145-A7BD-4966-94D5-B6284F28E504}[15]提出了奇異邊界法〔singularboundarymethod:SBM〕，該方法利用控制方程的根本解，通過反插值技術(shù)(inverseinterpolationtechnique)計算插值矩陣的對角元素，編程簡單，是一種真正的無網(wǎng)格法。本文對奇異邊界法中的兩種反插值技術(shù)進(jìn)展了研究，并將得到的結(jié)果與奇異無網(wǎng)格法做了比較。奇異邊界法的近似格式本文以二維Helmholtz方程為例描述奇異邊界法的根本技術(shù)路線：〔1〕〔2〕其中是待求未知量，為波數(shù)，是空間坐標(biāo)，和分別代表計算域及其邊界，和為函數(shù)。二維Helmholtz方程的根本解是〔3〕這里代表和兩點間的歐幾里得距離。根據(jù)根本解法的原理，以根本解為插值基函數(shù)，Helmholtz方程〔1〕的待求函數(shù)可以近似為〔4〕其中是邊界離散點的數(shù)目，為待求插值系數(shù)。當(dāng)配點和源點重疊時，不存在，即產(chǎn)生原點奇異性。為了防止根本解的奇異性，根本解方法的策略是將源點虛擬地布置在物理域以外的虛假邊界上，而將配點*i布置在真實的物理邊界上，即插值源點和配點是兩組完全不同的點。但到目前為止，對于復(fù)雜幾何域或多連通幾何域問題，如何較好地布置虛擬源點以保證計算結(jié)果可靠和穩(wěn)定收斂，仍是根本解方法中一個未能解決的關(guān)鍵問題。奇異邊界法不同于根本解方法的關(guān)鍵之處在于，插值源點和配點是同一組物理邊界的離散點，因而不存在根本解方法中的虛假邊界選取問題。奇異邊界法的插值公式為ADDINNE.Ref.{18FD56D2-72BD-4EA0-8EAE-5B91A989FF77}[15]：〔5〕這里是插值點的總數(shù)，是待定插值系數(shù)。值得注意的是，插值公式(5)和(4)都用根本解為基函數(shù)，但插值公式〔5〕在配點和源點重合處，假設(shè)了一個源點強度因子〔originintensityfactor〕。將插值公式(5)代入方程〔1〕和〔2〕，令，因為根本解滿足控制方程，我們得到下面矩陣形式的邊界條件離散代數(shù)方程：(6)源點強度因子實際上是插值矩陣的對角線元素。由于根本解的原點奇異性，我們不能夠簡單地用根本解插值基函數(shù)來計算。理想的方法是從數(shù)學(xué)理論上導(dǎo)出一個計算源點強度因子的公式，但目前這是一個非常有挑戰(zhàn)性的數(shù)學(xué)物理問題。下面我們通過反插值技術(shù)來求解源點強度因子。反插值技術(shù)〔inverseinterpolationtechnique〕注意到插值公式(5)和離散代數(shù)方程(6)中的待求插值系數(shù)與邊界上的配點分布，邊界條件和右邊項有關(guān)，而源點強度因子僅依賴于邊界條件和邊界上的配點分布，與右邊項無關(guān)。因而，我們可以設(shè)計一個反插值技術(shù)來計算對角線元素。對方程(1)和(2)所描述的Helmholtz有限域問題，在其物理域的邊界上布置個配點，在物理域內(nèi)部布置個計算輔助點。對于有根本解的控制方程，容易發(fā)現(xiàn)它們的一些特解。對Helmholtz方程，有許多的簡單特解，例如。利用插值公式(5)，我們有 (7)這里，是在邊界配點上的影響系數(shù)。因為輔助點和配點完全不重疊，因而沒有奇異性問題。由方程(7)，我們就可以求得影響系數(shù)。這里內(nèi)部輔助點的個數(shù)可以等于或多于邊界配點的數(shù)目。本文中采取兩種方案：方案1——在物理區(qū)域內(nèi)部布置與邊界配點相等數(shù)目的輔助點〔即〕，可以得到插值方陣；方案2——在物理區(qū)域內(nèi)部布置輔助點的數(shù)目多于邊界配點數(shù)〔即〕，可以得到插值矩陣，需用移動最小二乘近似求解。下一步，我們將計算輔助點換成邊界源點，即配點和源點完全重疊在邊界上；我們有 (8)這里插值矩陣A的非對角線元素可由公式得到。因而利用方程(7)中求得的系數(shù)，我們就能用方程(8)計算出關(guān)鍵的未知對角元素，即源點強度因子。利用上面得到的源點強度因子，我們就可以用邊界插值公式(5)，計算具有一樣幾何形狀和控制方程的任意問題。數(shù)值算例和討論為了驗證奇異邊界法的實際有效性并比較兩種反插值技術(shù)的可行性，本節(jié)中我們用該方法數(shù)值檢驗了Laplace方程和Helmholtz方程問題。當(dāng)式〔1〕中的時，可以得到Laplace方程 (9)本文所用的平均相對誤差定義如下ADDINNE.Ref.{A0CAEFCE-FBBF-448F-80A0-61B0EFE58CDC}[16]： (10)其中當(dāng)時，，當(dāng)時，，代表分布在計算域上的檢驗點的數(shù)目。我們在檢驗點上檢驗邊界奇異法的精度和收斂性。注意奇異邊界法是一個邊界離散數(shù)值方法，計算量最終只涉及了下面各圖中所示的邊界點數(shù)。對Laplace方程和Helmholtz方程問題，這里分別用和做為特解，計算邊界源點的源點強度因子。算例1：圓形域上的Laplace算例考慮以原點為中心的單位圓區(qū)域，解析解取為。用兩種反插值技術(shù)以及奇異無網(wǎng)格法〔圖中用雙層勢表示〕計算得到的平均相對誤差收斂曲線圖如圖1所示。方案2所用內(nèi)部輔助點數(shù)為475，檢驗點數(shù)。圖1算例1平均相對誤差〔RMRSE〕的收斂曲線圖從圖1我們可以看出用方案1求解的結(jié)果精度高且誤差收斂曲線收斂性較好；方案2求解的精度較差且誤差收斂曲線有高度震蕩現(xiàn)象；用奇異無網(wǎng)格法求解的結(jié)果精度比方案1精度低，雖然有較好的收斂性，但收斂速度較低。圖2為圖1誤差對應(yīng)的插值矩陣條件數(shù)曲線圖，從中我們可以看出方案1和方案2的條件數(shù)的變化不規(guī)則，隨著邊界節(jié)點數(shù)的增加總體有增加的趨勢。然而奇異無網(wǎng)格法對應(yīng)的條件數(shù)非常小，且隨著邊界點數(shù)的增加逐漸減小。圖2算例1對應(yīng)的條件數(shù)曲線圖算例2：方形域上的二維Laplace算例本算例考慮以為中心的正方形區(qū)域，方案2所用內(nèi)部輔助點數(shù)441，檢驗點數(shù)。用兩種反插值技術(shù)及奇異無網(wǎng)格法求解，圖3和圖4分別給出了對應(yīng)的誤差收斂曲線圖和條件數(shù)曲線圖。圖3說明當(dāng)邊界節(jié)點數(shù)時，方案2和方案1都有較好的收斂性。盡管方案2的精度比方案1的精度略高，但當(dāng)時，方案2得到的誤差收斂曲線產(chǎn)生振蕩現(xiàn)象且精度低于方案1。奇異無網(wǎng)格法比方案1的求解結(jié)果低一個精度。從圖4可以看出，方案1的插值矩陣條件數(shù)與求解精度大致成反比例的關(guān)系，即插值矩陣條件數(shù)增加一個數(shù)量級，對應(yīng)的數(shù)值解提高一個精度。同算例1，奇異邊界法的插值矩陣條件數(shù)非常小且沒有明顯變化。圖3算例2平均相對誤差〔RMRSE〕的收斂曲線圖圖4算例2對應(yīng)的條件數(shù)曲線圖算例3:不規(guī)則區(qū)域上的Helmholtz算例本算例考慮了不規(guī)則區(qū)域〔圖5〕ADDINNE.Ref.{07AA26B6-83E3-4A44-8CC3-4A7324A9A130}[17]上的Helmholtz方程，解析解取為。方案2所用內(nèi)部輔助點數(shù)350，檢驗點。由于奇異邊界法需要對邊界節(jié)點等間距分布，對圖5所示區(qū)域難于處理，因此本算例沒有給出該方法的結(jié)果。圖5算例3所用的不規(guī)則圖形圖6和圖7分別給出隨著邊界節(jié)點數(shù)的增加，平均相對誤差和條件數(shù)的曲線圖。從圖6可以看到方案1求得解的收斂性較好，而方案2求得解的收斂曲線有嚴(yán)重的振蕩現(xiàn)象。方案2的插值矩陣條件數(shù)比方案1的條件數(shù)高且隨著插值節(jié)點數(shù)的增加而增大〔圖7〕。圖6算例3平均相對誤差〔RMRSE〕的收斂曲線圖圖4算例3對應(yīng)的條件數(shù)曲線圖結(jié)論奇異邊界法不需要網(wǎng)格劃分和奇異積分、數(shù)學(xué)簡單、編程容易，直接利用控制方程的根本解為基函數(shù)，通過反插值技術(shù)計算插值矩陣的對角元素，是一種真正的無網(wǎng)格法。本文對奇異邊界法中的兩種反插值技術(shù)進(jìn)展了研究，通過數(shù)值算例比較了兩種反插值技術(shù)的求解精度及收斂性，并與奇異邊界法的結(jié)果做了比較。數(shù)值結(jié)果說明反插值技術(shù)方案1，即物理區(qū)域內(nèi)部輔助點數(shù)與邊界點數(shù)一樣時，求得解的結(jié)果精度較高且收斂性較好。ADDINNE.Bib參考文獻(xiàn)1 BelytschkoT,KrongauzY,OrganD,etal.Meshlessmethods:Anoverviewandrecentdevelopments.puterMethodsinAppliedMechanicsandEngineering,1996,139:3~47.2 NguyenVP,RabczukT,BordasS,etal.Meshlessmethods:Areviewandputerimplementationaspects.MathematicsandputersinSimulation,2008,79:763~813.3*雄，*巖，馬上.無網(wǎng)格法的理論及應(yīng)用.力學(xué)進(jìn)展,2009,39(1):1~36.4 Li*L,ZhuJL.Themethodoffundamentalsolutionsfornonlinearellipticproblems.EnginneringAnalysiswithBoundaryElements,2009,33(3):322~329.5 AntonioJ,TadeuA,GodinhoL.Athree-dimensionalacousticsmodelusingthemethodoffundamentalsolutions.EngineeringAnalysiswithBoundaryElements,2008,32:525~531.6 ChenCS,KarageorghisA,SmyrlisYS.TheMethodofFundamentalSolutions-AMeshlessMethod.2008:DynamicPublishers.7 FairweatherG,KarageorghisA.Themethodoffundamentalsolutionsforellipticboundaryvalueproblems.AdvancesinputationalMathematics,1998,9:69~95.8 PoullikkasA,KarageorghisA,GeorgiouG.Methodsoffundamentalsolutionsforharmonicandbiharmonicboundaryvalueproblems.putationalMechanics,1998,21:416~423.9 ChenW,TanakaM.AMeshless,Integration-Free,andBoundary-OnlyRBFTechnique.putersandMathematicswithApplications,2002,43:379~391.10 ChenW,HonYC.NumericalinvestigationonconvergenceofboundaryknotmethodintheanalysisofhomogeneousHelmholtz,modi?edHelmholtz,andconvection–diffusionproblems.puterMethodsinAppliedMechanicsandEngineering,2003,192:1859~1875.11 YoungDL,ChenKH,LeeCW.Novelmeshlessmethodforsolvingthepotentialproblemswitharbitrarydomain.JournalofputationalPhysics,2005,209:290~321.12 ChenKH,KaoJH,ChenJT,etal.Regularizedmeshlessmethodformultiply-connected-domainLaplaceproblems.EngineeringAnalysiswithBoundaryElements,2006,30:882~896.13 BozidarS.Amodifiedmethodoffundamentalsolutionsforpotentialproblem,inTheMethodofFundamentalSolutions-AMeshlessMethod,C.S.Chen,A.KarageorghisandY.S.Smyrlis,Editors.2008,DynamicPublishers.p.299~321.14 SongRC,ChenW.Aninvestigationontheregularizedmeshlessmethodforirregulardomainproblems.putationalMethodsinEngineeringandSciences,2009(Accepted).15陳文.奇異邊界法：一個新的、簡單、無網(wǎng)格、邊界配點數(shù)值方法.固體力學(xué)學(xué)報,2009(已接收).16 WangFZ,ChenW,Jiang*R.Investigationofregularizedtechniquesforboundaryknotmethod.municationsinNumericalMethodsinEngineering,2009(Accepted).17 ChenCS,ChoHA,GolbergMA.Somementsontheill-conditioningofthemethodoffundamentalsolutions.EngineeringAnalysiswithBoundaryElements,2006,30:405~410.TwoInverseInterpolationSchemesinSingularBoundaryMethodWenChen,FuzhangWang(Departmentofengineeringmechanics,Hohai