物理信息依賴核函數配點法的研究進展1)

傅卓佳 李明娟 習強 徐文志 劉慶國?,

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

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

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

引言

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

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

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

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

圖1 配點法采用不同試函數求解齊次微分方程問題式(5)和式(6)的計算結果.PIKF-12 表示在計算區域邊界上均勻布置12個節點并采用物理信息依賴核函數φ=作為試函數;RBF-16/36 表示在計算區域上均勻布置16/36 個節點并采用徑向基函數φ=r3作為試函數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

本文依次討論齊次微分方程、非齊次微分方程、非均質微分方程、非穩態微分方程和隱式微分方程這5 類數學物理方程.根據微分方程各自的特點,構造基本解、調和函數、徑向Trefftz 函數以及T 完備函數等包含微分方程物理信息的基函數[12-13],選取合適的配點技術,建立相應的物理信息依賴核函數配點法.最后,通過幾個典型算例驗證本文所提物理信息依賴核函數配點法的有效性.

1 物理信息依賴核函數配點法

1.1 齊次微分方程問題

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

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

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

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

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

1.2 非齊次微分方程問題

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

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

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

為了確保上述高階齊次微分控制方程(14)的解唯一,需要在計算區域邊界 ? Ω 上施加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 常見高階微分方程算子的調和函數 和徑向Trefftz 函數 [12]Table 6 Harmonic functions and radial Trefftz functionsof common-used high-order differential equation operators[12]

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

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

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

1.3 非均質微分方程問題

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

隨后可以直接采用上兩節中構造的物理信息依賴核函數配點法求解變量變換后的各向同性微分方程問題.

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

圖3 物理信息依賴核函數配點法的節點離散及局部子區域示意圖:(a) 節點離散及節點 的局部子區域 Ξk;(b)無奇異物理信息依賴核函數的局部子區域節點離散;(c)奇異物理信息依賴核函數的局部子區域節點離散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

其相應的矩陣向量表達式為

1.4 非穩態微分方程問題

(1) 時間步進法

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

(2) 變換解法

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

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

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

(3) 直接解法

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

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

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

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

1.5 隱式微分方程問題

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

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

類似地,包含非Fick 擴散過程物理信息的基本解可由其粒子運動的統計分布概率密度函數和Heaviside函數表示

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

2 典型數值算例與討論

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

2.1 無限域波傳播問題

本算例采用基于物理信息依賴核函數的局部-全域混合配點法計算三維結構振動引起的聲輻射響應.由于空氣對球體結構振動影響較小,因此忽略空氣對結構的振動影響,考慮如下微分方程問題

在具體計算中,局部物理信息依賴核函數配點法被用于三維結構振動分析,全域邊界型物理信息依賴核函數配點法被用于結構外聲輻射場計算.其中包含三維結構振動方程物理信息的基本解為

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

圖4 球體結構受不同頻率簡諧激勵力作用在點( 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 簡諧激勵力下球體結構附近區域的聲壓級分布Fig.5 Sound pressure level distribution around the spherical structure under simple harmonic excitation force at frequency 100Hz

2.2 大尺寸比結構問題

考慮如圖6 所示的功能梯度材料變半徑圓環瞬態熱傳導問題

圖6 功能梯度材料變半徑圓環區域示意圖Fig.6 Schematic configuration of variable radius FGM ring

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

圖7 SR=100的變半徑圓環中軸線 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 物理信息依賴核函數配點法求解多個時刻不同尺寸比情況下的計算結果(Merr)Table 10 Numerical results(Merr) obtained by using PIKF collocation method at varied time instants under different SRs

2.3 工程反演問題

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

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

式中u和q可由精確解推導得到,r andn(i)為滿足標準正態分布的一組隨機數,e為含噪聲水平.

在物理信息依賴核函數配點法計算(調和函數及徑向Trefftz 函數)中,采用1.2 節介紹的多重互易法確定式(13)中的微分算子L1=Δ-3,選取高階調和函數及徑向Trefftz 函數作為基函數,其中形狀參數c=0.02 ;同時引入截斷奇異值分解技術[31]處理離散得到的病態矩陣方程問題,其中截斷參數kSVD由廣義交叉校驗法確定.由圖9和表11 可知,物理信息依賴核函數配點法結合截斷奇異值技術僅需部分邊界數據即可有效反演位勢問題的解,同時該計算方法對噪聲水平在5%以下的噪音數據具有較好的魯棒性.

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

表11 物理信息依賴核函數配點法采用可測邊界上含5%人工噪音數據反演得到的結果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 移動邊界問題

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

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

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

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

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

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

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

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

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

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

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

3 結論

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

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