物理信息依賴核函數(shù)配點(diǎn)法的研究進(jìn)展1)

傅卓佳 李明娟 習(xí)強(qiáng) 徐文志 劉慶國?,

*(河海大學(xué)力學(xué)與材料學(xué)院,工程與科學(xué)數(shù)值模擬軟件中心,南京 211100)

?(盧布爾雅那大學(xué)機(jī)械工程學(xué)院,斯洛文尼亞盧布爾雅那1000)

**(金屬與技術(shù)研究所,斯洛文尼亞盧布爾雅那1000)

引言

近幾十年來,計(jì)算機(jī)數(shù)值模擬逐漸成為解決科學(xué)和工程問題的主流方法,其與理論分析、實(shí)驗(yàn)研究和數(shù)據(jù)探索一起成為當(dāng)前科學(xué)研究的4 類研究手段之一[1].如表1 所示,有限元法和有限差分法等傳統(tǒng)網(wǎng)格類數(shù)值離散方法作為目前主流的計(jì)算機(jī)數(shù)值模擬方法已被廣泛用于眾多科學(xué)與工程領(lǐng)域.而無網(wǎng)格與粒子方法、邊界元法和離散元法等數(shù)值離散方法各具特色,在傳統(tǒng)數(shù)值計(jì)算方法求解不理想的特大變形、無限域、顆粒介質(zhì)等物理力學(xué)問題中占據(jù)了一席之地.

表1 Web of Science 數(shù)據(jù)庫中各計(jì)算方法的論文發(fā)表數(shù)(搜索日期:2022-9-5)Table 1 Bibliographic database search of computational methods based on the Web of Science(search date:September 5,2022)

本文主要研究無網(wǎng)格與粒子方法中一類最古老且簡單的配點(diǎn)法[2-10].以下述微分方程為例

其中 ?和B 分別表示微分控制方程算子和邊界條件方程算子,u(x)為計(jì)算域上x點(diǎn)處的待求物理量,f(x)和g(x) 分別為計(jì)算域內(nèi)的已知源項(xiàng)激勵和計(jì)算域邊界上的已知邊界條件.對于上述微分方程問題式(1)和式(2),引入試函數(shù)矩陣 Φ={φ}和待求系數(shù)向量 α,構(gòu)建微分方程數(shù)值解=Φα,采用加權(quán)殘量法[11]進(jìn)行計(jì)算得到

圖1 配點(diǎn)法采用不同試函數(shù)求解齊次微分方程問題式(5)和式(6)的計(jì)算結(jié)果.PIKF-12 表示在計(jì)算區(qū)域邊界上均勻布置12個節(jié)點(diǎn)并采用物理信息依賴核函數(shù)φ=作為試函數(shù);RBF-16/36 表示在計(jì)算區(qū)域上均勻布置16/36 個節(jié)點(diǎn)并采用徑向基函數(shù)φ=r3作為試函數(shù)Fig.1 Numerical results by using collocation method with different trial functions for solving homogeneous differential equation problem Eqs.(5) and(6).PIKF-12 representsis the trial function and 12 nodes are placed at the boundary of computational domain;RBF-16/36 representsφ=r3 is the trial function and 16/36 nodes are uniformly placed in computational domain

本文依次討論齊次微分方程、非齊次微分方程、非均質(zhì)微分方程、非穩(wěn)態(tài)微分方程和隱式微分方程這5 類數(shù)學(xué)物理方程.根據(jù)微分方程各自的特點(diǎn),構(gòu)造基本解、調(diào)和函數(shù)、徑向Trefftz 函數(shù)以及T 完備函數(shù)等包含微分方程物理信息的基函數(shù)[12-13],選取合適的配點(diǎn)技術(shù),建立相應(yīng)的物理信息依賴核函數(shù)配點(diǎn)法.最后,通過幾個典型算例驗(yàn)證本文所提物理信息依賴核函數(shù)配點(diǎn)法的有效性.

1 物理信息依賴核函數(shù)配點(diǎn)法

1.1 齊次微分方程問題

表2 常見微分方程算子的基本解 GF[12]Table 2 Fundamental solutions GFof common-used differential equation operators[12]

表3 常見微分方程算子的調(diào)和函數(shù) GH和徑向Trefftz 函數(shù) GRT[12]Table 3 Harmonic functions GHand radial Trefftz functions GRTof common-used differential equation operators[12]

表4 常見微分方程算子的T 完備函數(shù) GT[13]Table 4 T-complete functions GTof common-used differential equation operators[13]

圖2 物理信息依賴核函數(shù)配點(diǎn)法的節(jié)點(diǎn)離散示意圖:(a) 奇異基本解 GF和(b)無奇異物理信息依賴核函數(shù)(GH,GRT,GT)Fig.2 Sketch of node discretization of PIKF collocation methods:(a) singular fundamental solutions GFand(b) nonsingular PIKFs(GH,GRT,GT)

由于虛假邊界系數(shù)d的選取對基于基本解的配點(diǎn)法(基本解法[14])的計(jì)算精度和穩(wěn)定性有很大的影響,因此國內(nèi)外學(xué)者對如何選取合適的虛假邊界系數(shù)進(jìn)行了大量的研究.同時也有一些學(xué)者通過邊界積分方程和數(shù)值反插值等技術(shù)構(gòu)造去奇異基本解來避免虛假邊界系數(shù)d的選取.有興趣的讀者可以查閱參考文獻(xiàn)[15-18].

1.2 非齊次微分方程問題

上一節(jié)中討論了f(x)=0的情況,那么對于f(x)≠0的微分方程問題(1)和(2),是否也能構(gòu)造得到自動滿足非齊次微分控制方程(1)的物理信息依賴核函數(shù)？答案是肯定的.

一般來說,非齊次微分方程問題式(1)和式(2)的數(shù)值解可分解為齊次解和特解兩部分,即

則可將方程(10)轉(zhuǎn)化成高階齊次微分控制方程

為了確保上述高階齊次微分控制方程(14)的解唯一,需要在計(jì)算區(qū)域邊界 ? Ω 上施加N組附加約束條件

表5 常見高階微分方程算子的基本解 [12]Table 5 Fundamental solutions of common-used highorder differential equation operators[12]

表5 常見高階微分方程算子的基本解 [12]Table 5 Fundamental solutions of common-used highorder differential equation operators[12]

表6 常見高階微分方程算子的調(diào)和函數(shù) 和徑向Trefftz 函數(shù) [12]Table 6 Harmonic functions and radial Trefftz functionsof common-used high-order differential equation operators[12]

表6 常見高階微分方程算子的調(diào)和函數(shù) 和徑向Trefftz 函數(shù) [12]Table 6 Harmonic functions and radial Trefftz functionsof common-used high-order differential equation operators[12]

表7 常見高階微分方程算子的T 完備函數(shù) [23]Table 7 T-complete functions of common-used high-order differential equation operators[23]

表7 常見高階微分方程算子的T 完備函數(shù) [23]Table 7 T-complete functions of common-used high-order differential equation operators[23]

1.3 非均質(zhì)微分方程問題

上兩節(jié)中討論的微分方程算子 ? 都是均勻各向同性的.對于非均質(zhì)微分方程問題,本節(jié)將以熱傳導(dǎo)方程為例進(jìn)行介紹.對于一些特殊的非均質(zhì)微分方程問題,可以通過變量變換技術(shù)將其轉(zhuǎn)換成均勻各向同性微分方程問題進(jìn)行計(jì)算.

隨后可以直接采用上兩節(jié)中構(gòu)造的物理信息依賴核函數(shù)配點(diǎn)法求解變量變換后的各向同性微分方程問題.

若熱傳導(dǎo)系數(shù)矩陣Kij(x) 是一般形式的函數(shù)矩陣,則考慮引入局部子區(qū)域思想[24],使得Kij(x) 在局部子區(qū)域 Ξk內(nèi)具有上述幾種特殊形式的函數(shù)矩陣抑或可以近似簡化成常數(shù)矩陣.進(jìn)而構(gòu)造在局部子區(qū)域 Ξk上自動滿足微分控制方程(1)的物理信息依賴核函數(shù)對問題進(jìn)行求解.

圖3 物理信息依賴核函數(shù)配點(diǎn)法的節(jié)點(diǎn)離散及局部子區(qū)域示意圖:(a) 節(jié)點(diǎn)離散及節(jié)點(diǎn) 的局部子區(qū)域 Ξk;(b)無奇異物理信息依賴核函數(shù)的局部子區(qū)域節(jié)點(diǎn)離散;(c)奇異物理信息依賴核函數(shù)的局部子區(qū)域節(jié)點(diǎn)離散Fig.3 Sketch of node discretization and local subdomain of PIKF collocation method:(a) node discretization and local subregion Ξkfor node;(b) node discretization of local subdomain Ξkfor nonsingular PIKFs;(c) node discretization of local subdomain Ξkfor singular PIKFs

其相應(yīng)的矩陣向量表達(dá)式為

1.4 非穩(wěn)態(tài)微分方程問題

(1) 時間步進(jìn)法

以向前差分格式(m=1)為例,可離散得

(2) 變換解法

另一類常用的時間離散技術(shù)為變換解法,其主要思想是采用Laplace 變換、Fourier 變換等變換技術(shù)將原時間依賴微分方程問題轉(zhuǎn)化成變換空間下的時間無關(guān)微分方程問題,然后采用相應(yīng)的數(shù)值逆變換將變換空間下的微分方程問題解轉(zhuǎn)化為原時間依賴微分方程問題解.本文以Laplace 變換為例進(jìn)行介紹,采用如下Laplace 變換定義

可將微分方程問題式(28)和式(30)轉(zhuǎn)化為

其中p為Laplace 變換系數(shù),物理量上標(biāo)(L)表示為Laplace 變換后的變量.可采用前文中介紹的方法根據(jù)微分算子和右端源項(xiàng)構(gòu)造相應(yīng)的物理信息依賴核函數(shù),并建立對應(yīng)的配點(diǎn)法計(jì)算得到u(L)(x,p).隨后選取合適的數(shù)值Laplace逆變換得到時間依賴微分方程問題的解u(x,t).需要注意的是,數(shù)值Laplace 逆變換本質(zhì)上屬于曲線擬合,使用時需要謹(jǐn)慎選取變換參數(shù)p的數(shù)目.由于篇幅所限,更詳細(xì)的數(shù)值Laplace 逆變換介紹可參考文獻(xiàn)[27].

(3) 直接解法

另一種時間離散技術(shù)為直接解法[28].采用該方法的關(guān)鍵在于構(gòu)造自動滿足時間依賴微分方程(28)的物理信息依賴核函數(shù),如表8和表9 所示.隨后,類似式(7),可用一組時間依賴物理信息核函數(shù)的線性組合表示時間依賴微分方程問題式(28)～ 式(30)的解

表8 時間依賴微分方程算子的基本解GFTable 8 Fundamental solutions GFof time-dependent differential equation operators

表9 時間依賴微分方程算子的徑向Trefftz 函數(shù)GRTTable 9 Radial Trefftz functions GRTof time-dependent differential equation operators

其中試函數(shù)φ(x,t,sj,τn) 可選為基本解和徑向Trefftz函數(shù)等物理信息依賴核函數(shù).此時,只需對邊界條件方程(29)和初始條件方程(30)進(jìn)行離散,有效降低了計(jì)算成本.

1.5 隱式微分方程問題

在本文之前討論的經(jīng)典微分方程問題中,其微分控制方程的表達(dá)式均可被顯式表示,此類顯式的經(jīng)典微分方程已被廣泛用于描述眾多物理力學(xué)現(xiàn)象.然而,最近的一些研究發(fā)現(xiàn)反常擴(kuò)散和頻率冪律依賴聲波衰減等現(xiàn)象難以用上述顯式經(jīng)典微分方程來描述.本節(jié)將介紹能描述這些反常物理力學(xué)現(xiàn)象的隱式微分方程建模方法[29],其基本思路是直接從物理力學(xué)問題出發(fā)構(gòu)造相應(yīng)的物理信息依賴核函數(shù),跳過建立相應(yīng)微分方程模型的過程.這里“隱式”是指微分控制方程的顯式表達(dá)式難以推導(dǎo)得到甚至不需要.

也就是說,經(jīng)典Fick 擴(kuò)散方程的物理信息依賴基本解可由其粒子運(yùn)動的高斯分布概率密度函數(shù)和Heaviside 函數(shù)表示.

類似地,包含非Fick 擴(kuò)散過程物理信息的基本解可由其粒子運(yùn)動的統(tǒng)計(jì)分布概率密度函數(shù)和Heaviside函數(shù)表示

此時,采用直接解法將物理信息依賴基本解(39)代入計(jì)算表達(dá)式(36),并只需對邊界條件方程(29)和初始條件方程(30)進(jìn)行離散確定待定系數(shù) αj,n,即可數(shù)值模擬此類反常擴(kuò)散過程.

2 典型數(shù)值算例與討論

本節(jié)將通過四個典型算例驗(yàn)證物理信息依賴核函數(shù)配點(diǎn)法的有效性.為了量化數(shù)值計(jì)算結(jié)果,分別定義平均相對誤差Rerr和最大相對誤差Merr為

2.1 無限域波傳播問題

本算例采用基于物理信息依賴核函數(shù)的局部-全域混合配點(diǎn)法計(jì)算三維結(jié)構(gòu)振動引起的聲輻射響應(yīng).由于空氣對球體結(jié)構(gòu)振動影響較小,因此忽略空氣對結(jié)構(gòu)的振動影響,考慮如下微分方程問題

在具體計(jì)算中,局部物理信息依賴核函數(shù)配點(diǎn)法被用于三維結(jié)構(gòu)振動分析,全域邊界型物理信息依賴核函數(shù)配點(diǎn)法被用于結(jié)構(gòu)外聲輻射場計(jì)算.其中包含三維結(jié)構(gòu)振動方程物理信息的基本解為

考慮圖4 所示的球體結(jié)構(gòu)受沿外法向方向大小為 200N/m2的均布面力荷載作用情況,其中球體結(jié)構(gòu)半徑為0./5m,彈性模量為E=210GPa,密度為ρE=7800kgm3,泊松比為vE=0.3.為了驗(yàn)證物理信息依賴核函數(shù)的局部-全域混合配點(diǎn)法的有效性,此算例采用COMSOL 軟件(FEM)的仿真結(jié)果作為參考解.其中物理信息依賴核函數(shù)的局部-全域混合配點(diǎn)法計(jì)算時在球體結(jié)構(gòu)區(qū)域布置了3553 個節(jié)點(diǎn)用于結(jié)構(gòu)振動分析,并在球體結(jié)構(gòu)表面布置了1000個節(jié)點(diǎn)用于結(jié)構(gòu)外聲輻射分析;而COMSOL 軟件計(jì)算時對球體結(jié)構(gòu)附近的截?cái)鄥^(qū)域劃分了178 497 個四面體網(wǎng)格和9298 個三角形網(wǎng)格.圖4 給出了球體結(jié)構(gòu)在頻率為 50～200Hz 的簡諧激勵力作用下振動引起的點(diǎn)( 2,0,0) 處的聲壓級變化結(jié)果.圖5 給出了受100Hz 簡諧激勵力下球體結(jié)構(gòu)附近區(qū)域的聲壓級分布.由圖4和圖5 中可以看出物理信息依賴核函數(shù)配點(diǎn)法計(jì)算得到的聲壓級結(jié)果與COMSOL 軟件的仿真結(jié)果吻合較好.

圖4 球體結(jié)構(gòu)受不同頻率簡諧激勵力作用在點(diǎn)( 2,0,0) 處的聲壓級變化Fig.4 Sound pressure level results of spherical structure under simple harmonic excitation force on(2,0,0) with different frequencies

圖5 受100Hz 簡諧激勵力下球體結(jié)構(gòu)附近區(qū)域的聲壓級分布Fig.5 Sound pressure level distribution around the spherical structure under simple harmonic excitation force at frequency 100Hz

2.2 大尺寸比結(jié)構(gòu)問題

考慮如圖6 所示的功能梯度材料變半徑圓環(huán)瞬態(tài)熱傳導(dǎo)問題

圖6 功能梯度材料變半徑圓環(huán)區(qū)域示意圖Fig.6 Schematic configuration of variable radius FGM ring

隨后采用1.2 節(jié)介紹的多重互易法確定式(13)中的微分算子L2=L1=Δ-1.最后采用Fixed Talbot數(shù)值Laplace 逆變換[27]得到原問題的解.表10列出了物理信息依賴核函數(shù)配點(diǎn)法求解多個時刻不同尺寸比情況下的計(jì)算結(jié)果.圖7 給出SR=100的變半徑圓環(huán)中軸線rc=R-0.475+0.45θ/π 上不同時刻(t=0.4,2,10s)的溫度分布.由表10和圖7 可知,物理信息依賴核函數(shù)配點(diǎn)法采用T-完備函數(shù)作為基函數(shù),無需布置非常稠密的節(jié)點(diǎn),僅在變半徑圓環(huán)區(qū)域邊界布置91 個離散節(jié)點(diǎn)即可得到高精度的結(jié)果.需要指出的是,表10中同一時刻不同尺寸比情況下本文配點(diǎn)法得到相似的計(jì)算結(jié)果誤差,這是由于物理信息依賴核函數(shù)配點(diǎn)法計(jì)算得到了非常精確的頻域解,而計(jì)算結(jié)果誤差主要產(chǎn)生于Fixed Talbot 數(shù)值Laplace 逆變換過程中.

圖7 SR=100的變半徑圓環(huán)中軸線 rc=R-0.475+0.45θ/π 上不同時刻的溫度分布Fig.7 Temperature distributions along with the curve rc=R-0.475+0.45θ/πinside the FGM under different time instants

表10 物理信息依賴核函數(shù)配點(diǎn)法求解多個時刻不同尺寸比情況下的計(jì)算結(jié)果(Merr)Table 10 Numerical results(Merr) obtained by using PIKF collocation method at varied time instants under different SRs

2.3 工程反演問題

考慮如圖8 所示類腫瘤區(qū)域的各向異性位勢柯西反問題

圖8 三維類腫瘤區(qū)域及可測邊界測量點(diǎn)示意圖Fig.8 Schematic configuration of 3D tumor-like domain with measurement points(○).

式中u和q可由精確解推導(dǎo)得到,r andn(i)為滿足標(biāo)準(zhǔn)正態(tài)分布的一組隨機(jī)數(shù),e為含噪聲水平.

在物理信息依賴核函數(shù)配點(diǎn)法計(jì)算(調(diào)和函數(shù)及徑向Trefftz 函數(shù))中,采用1.2 節(jié)介紹的多重互易法確定式(13)中的微分算子L1=Δ-3,選取高階調(diào)和函數(shù)及徑向Trefftz 函數(shù)作為基函數(shù),其中形狀參數(shù)c=0.02 ;同時引入截?cái)嗥娈愔捣纸饧夹g(shù)[31]處理離散得到的病態(tài)矩陣方程問題,其中截?cái)鄥?shù)kSVD由廣義交叉校驗(yàn)法確定.由圖9和表11 可知,物理信息依賴核函數(shù)配點(diǎn)法結(jié)合截?cái)嗥娈愔导夹g(shù)僅需部分邊界數(shù)據(jù)即可有效反演位勢問題的解,同時該計(jì)算方法對噪聲水平在5%以下的噪音數(shù)據(jù)具有較好的魯棒性.

圖9 物理信息依賴核函數(shù)配點(diǎn)法采用可測邊界 Γ1 上不同人工噪音水平的100個測量點(diǎn)數(shù)據(jù)反演得到位勢問題的結(jié)果Fig.9 Numerical results obtained by using PIKF collocation method with 100measurement points under different noise levels on accessible boundary Γ1

表11 物理信息依賴核函數(shù)配點(diǎn)法采用可測邊界上含5%人工噪音數(shù)據(jù)反演得到的結(jié)果Table 11 Numerical results obtained by using PIKF collocation method with different boundary measurement points on the accessible boundary under 5% noise level

2.4 移動邊界問題

三維數(shù)值波浪水槽(圖10)長度b=30m,寬度w=0.1 m,靜水深h=0.4 m.水槽邊界主要分為自由表面邊界條件 Γ1、入射邊界 Γ2、底部邊界 Γ3、吸收邊界 Γ4以及側(cè)面邊界 Γ5.考慮無黏無旋不可壓縮理想流體,水槽區(qū)域 ΩW的控制方程為

圖10 三維數(shù)值波浪水槽和梯形潛堤橫向剖面示意圖Fig.10 Schematic diagram of 3D numerical wave flume and transverse profile of trapezoidal submerged breakwater

其中φW為速度勢,ΩW為水槽區(qū)域.在自由表面邊界Γ1上,采用半拉格朗日觀點(diǎn)得到考慮海綿層的動力和運(yùn)動邊界條件

其中 νW(x1)為阻尼系數(shù),用于防止入射波反射

式中 ωW為入射波圓頻率,這里采用二階顯式龍格庫塔法進(jìn)行離散.在邊界 Γ4施加輻射邊界條件

其中,入射波相速度CW=λW/T,λW為入射波波長,T為入射波周期.在入射邊界 Γ2施加邊界條件

式中U(x3,t)為沿x方向水平速度U,斜坡函數(shù)

不透水側(cè)邊界 Γ3和底部邊界 Γ5滿足無通量條件

其中n為沿不透水邊界的向外單位法向量.

在上游邊界 Γ2采用二階斯托克斯波作為入射波,其中周期T=2 s,波長λW=3.693 m,波高H=0.02 m.數(shù)值水槽計(jì)算域內(nèi)布點(diǎn)間距dx1=dx2=dx3=0.025 m,時間步長 dt=0.01 s,鄰近點(diǎn)數(shù)nx=ns=40,總時長為12T,海綿層厚度等于波長 λW.由圖11 可知,基于局部物理信息依賴核函數(shù)配點(diǎn)法(Laplace方程基本解)建立的數(shù)值波浪水槽能有效模擬自由水面高程變化,與實(shí)驗(yàn)數(shù)據(jù)[32]吻合較好.

圖11 最后兩個周期內(nèi)三維數(shù)值波浪水槽自由表面兩個觀測點(diǎn)處(G3和G5)的高程演化Fig.11 Elevation evolutions at two gauges(G3 and G5) of 3D numerical wave flume in the last 2 periods

3 結(jié)論

本文介紹了一類基于物理信息依賴核函數(shù)的無網(wǎng)格配點(diǎn)法.該類方法的核心思想在于構(gòu)建包含問題微分控制方程物理信息的基函數(shù),隨后無需/僅需少量配點(diǎn)對所求微分控制方程進(jìn)行離散,能有效提高計(jì)算效率.通過無限域波傳播、大尺寸比結(jié)構(gòu)、工程反演和移動邊界4 個典型算例,驗(yàn)證了本文所提物理信息依賴核函數(shù)配點(diǎn)法的有效性.

需要指出的是,本文主要討論了線性物理力學(xué)問題解連續(xù)/分段連續(xù)的情況,對于解不連續(xù)及非線性物理力學(xué)問題,所提物理信息依賴核函數(shù)配點(diǎn)法需要附加解不連續(xù)處理技術(shù)[33]和非線性處理技術(shù)[12]進(jìn)行求解,目前這方面的研究工作還較少,有待進(jìn)一步研究.此外,本文從齊次微分方程物理信息依賴核函數(shù)開始,依次介紹了非齊次、非均質(zhì)、非穩(wěn)態(tài)微分方程構(gòu)造物理信息依賴核函數(shù)的方法,最后根據(jù)物理力學(xué)問題統(tǒng)計(jì)意義下的概率密度函數(shù)構(gòu)建物理信息依賴核函數(shù),無需建立微分方程問題的顯式表達(dá)式.更進(jìn)一步地,是否可以直接根據(jù)已有數(shù)據(jù)特點(diǎn)構(gòu)建物理信息依賴核函數(shù),進(jìn)而計(jì)算得到物理力學(xué)問題的解？這也有待進(jìn)一步研究.