張尚元， 聶玉峰， 李義強

（西北工業(yè)大學 數(shù)學與統(tǒng)計學院，西安 710129）

引 言

經(jīng)典的擴散現(xiàn)象和彈性力學問題是通過偏微分方程(partial differential equation，PDE)來刻畫的，由于模型在建立之初就默認解具有一定的連續(xù)性，因此這類經(jīng)典模型不適用于解具有奇異性的情形.為此，研究者們提出了非局部擴散模型[1-2]和近場動力學 (peridynamics，PD)模型[3-5]，這兩類非局部模型方程都是以有限作用范圍的積分形式呈現(xiàn)，并且不涉及未知量的空間導數(shù)，物體的內(nèi)力采取非接觸的相互作用，這些特點使得它們在表征材料的損傷以及不連續(xù)現(xiàn)象時具有明顯的優(yōu)勢，已經(jīng)被廣泛應(yīng)用于含裂紋的熱擴散問題[1,6]和固體損傷破壞[7]領(lǐng)域中.

為推動這類有限作用范圍的非局部模型在工程實際問題中的應(yīng)用，國內(nèi)外許多專家嘗試了不同數(shù)值方法來求解這類方程.2005 年，Silling 等通過將積分直接近似為有限和，提出了基于求積的無網(wǎng)格配置點法[8].Chen 等在2011 年使用連續(xù)有限元和不連續(xù)有限元方法求解了一維非局部擴散問題和近場動力學問題，指出在解不連續(xù)的情形下，使用不連續(xù)有限元可以獲得更高的精度[9].針對非局部方程計算量大的問題，Wang 和Tian 等基于剛度矩陣的Toeplitz 結(jié)構(gòu)，給出了一維非局部擴散方程的快速配置點法[10-11]，并在后續(xù)工作中推廣到了二維問題[12-13].2013 年，Tian 和Du 通過將非局部算子看成差分算子的加權(quán)平均，提出了基于求積的有限差分方法，用于一維非局部擴散問題的求解[14].2016 年，Lehoucq 等使用徑向基函數(shù)Galerkin 方法求解了二維非局部擴散問題[15].2017 年，Zhao 等使用局部徑向基函數(shù)偽譜方法求解了二維非局部擴散問題[16]，并在之后推廣到了近場動力學問題[17].2018 年，Pasetto 等將再生核方法應(yīng)用到二維近場動力學問題的數(shù)值求解當中[18].2019 年，Du 和Yin 給出了非局部擴散問題的一種協(xié)調(diào)間斷Galerkin 方法[19].2020 年，Zhang 和Nie 將POD 降階方法用于非穩(wěn)態(tài)非局部擴散問題，加速了非局部離散系統(tǒng)的求解過程[20].2021 年，Liang 等發(fā)展了近場動力學問題的邊界元方法[21].同年，Lu 和Nie 發(fā)展了具有再生性的局部化徑向基函數(shù)方法來求解非局部擴散問題，解決了再生核增強徑向基函數(shù)方法在求解非局部擴散方程時不收斂的問題[22].

由于這兩類模型方程的非局部特性，經(jīng)典的數(shù)值方法均面臨著算法復雜、難以實現(xiàn)和計算量大的問題.對于基于網(wǎng)格的方法，包括分片多項式函數(shù)近似的配置點法[10]，以及連續(xù)/不連續(xù)有限元方法[9]，由于非局部相互作用范圍是一個有限范圍的球形鄰域，因此這些方法在方程離散時，需要處理網(wǎng)格和球形區(qū)域的求交問題，否則將會降低數(shù)值精度.對于基于方程弱形式的方法(包括有限元方法[9]、無網(wǎng)格Galerkin 方法[15]等)，由于模型方程積分的特性，這就導致方程的弱形式會額外嵌套一層積分，從而在剛度矩陣的形成過程中會出現(xiàn)累次重積分的計算，給問題的求解帶來了較大的計算負擔.然而，基于徑向基函數(shù)的無網(wǎng)格方法[23-24]能夠避免處理上述提及的區(qū)域求交的問題.同時，基于方程強形式的配置點法，可以避免額外一層積分的計算，從而可以節(jié)省剛度矩陣生成時的計算量.基于上述分析，無網(wǎng)格配置點法更適合于這兩類非局部模型的數(shù)值求解.值得注意的是，全局RBF 方法雖然避免了網(wǎng)格生成，但是在計算RBF 插值時，由于插值矩陣規(guī)模較大，會出現(xiàn)計算量大、插值矩陣條件數(shù)過大的問題，而局部化的RBF 則會降低插值矩陣的條件數(shù)，故選擇局部化的徑向基函數(shù)插值更為合理.因此，本文考慮了一種局部化的徑向基函數(shù)方法-單位分解徑向基函數(shù)方法來求解這類方程：一方面它不同于文獻[15]中基于方程弱形式的無網(wǎng)格Galerkin 方法，是基于方程強形式的無網(wǎng)格配置點法，節(jié)省了額外一層的積分計算；另一方面，本文方法不同于有限元方法[9]，基函數(shù)無網(wǎng)格依賴，免去了區(qū)域求交問題.故本文方法在求解非局部問題時具備實施簡便、計算量小的優(yōu)勢.

針對二維非局部擴散問題和近場動力學問題，本文采用單位分解徑向基函數(shù)方法來數(shù)值求解.首先概述了RBF 插值和RBF-PU 插值方法；然后給出了RBF-PU 方法應(yīng)用于非局部擴散問題和近場動力學問題求解時離散系統(tǒng)形成的基本過程；最后用數(shù)值實驗驗證了RBF-PU 方法在求解這兩類非局部問題時的有效性.

1 徑向基函數(shù)插值

1.1 全局徑向基函數(shù)插值

1.2 單位分解徑向基函數(shù)插值

圖1 無網(wǎng)格點和PU 覆蓋示意圖Fig. 1 A set of regular meshless points and a set of circular PU patches

定義指標集I(x)={l:x∈?l}，則式(6)中的求和索引可以用k∈I(x)簡化.用Wendland 函數(shù)φ (r)來構(gòu)造函數(shù)φj(x)：

2 非局部問題的RBF-PU 方法

2.1 圓盤上的數(shù)值積分

在PD 的數(shù)值計算中，會大量涉及到圓盤上的數(shù)值積分.為了降低計算復雜度，通常的做法是將圓盤的積分轉(zhuǎn)化為矩形上的積分（如圖2），矩形上的數(shù)值積分直接由一維積分的張量積來獲得.對于二維問題，利用極坐標變換：

圖2 極坐標變換Fig. 2 Polar transformation

有了這些基本工具，下面給出RBF-PU 方法形成離散系統(tǒng)的過程.

2.2 非局部擴散方程的RBF-PU 方法

二維非局部擴散方程為

2.3 近場動力學方程的RBF-PU 方法

考慮向量形式的二維近場動力學運動方程：

3 數(shù) 值 實 驗

3.1 非局部擴散方程

算例1 取問題(16)的精確解為

核函數(shù)為

對區(qū)域采用兩種離散方式：一種是在每個方向上均勻離散為N份，均勻離散下離散點的示意圖和生成的矩陣結(jié)構(gòu)如圖3 所示；另一種采用相同數(shù)量的Halton點進行離散，離散點的分布和形成的剛度矩陣的結(jié)構(gòu)如圖4 所示.分別對兩種離散形式沿每個方向分別取10～60 個離散點(以1 0為增量)，進行數(shù)值計算，均勻離散的數(shù)值結(jié)果如表1 所示.從矩陣條件數(shù)看，RBF-PU 方法的條件數(shù)在離散點數(shù)量為60×60時為1010～1011量級，而全局RBF 方法參見文獻[22]表7 中的結(jié)果為1014～1015量級，由此可見，RBF-PU 方法在一定程度上降低了剛度矩陣的條件數(shù).從求解時間上看，RBF-PU 方法在離散點數(shù)量為60×60時的計算時間約為25 min，對同等規(guī)模的問題，全局RBF 方法約耗時10 h 4 min(參見文獻[22]表6 中的結(jié)果)，基于分片線性基函數(shù)的配置點法約耗時14 h 21 min(參見文獻[13]表1 中的結(jié)果)，可以看出RBF-PU 方法求解時間更短，效率更高.在表1 最后一行60×60離散下的誤差反而比50×50離散下的誤差大，原因有兩點：一是影響RBF-PU 方法精度的參數(shù)包括點的分布，PU 塊內(nèi)點的數(shù)量，由于PU 塊的大小沒有隨著離散點的增加而改變，因此導致插值矩陣進一步惡化，導致誤差增大.二是由于方程本身的核函數(shù)奇異，影響數(shù)值積分的精度，由于數(shù)值積分點數(shù)量固定，在離散點數(shù)目較小的時候，數(shù)值離散的誤差占主導地位，但隨著離散點數(shù)目增加，數(shù)值積分的誤差開始占據(jù)主導地位，最終導致誤差增大，這樣的現(xiàn)象在文獻[16]表2、3 中也可以觀察到.圖5 給出了在50×50的情況下數(shù)值結(jié)果和誤差的分布圖.Halton點離散下的數(shù)值結(jié)果如表2 所示.圖6 給出了50×50的情況下的數(shù)值結(jié)果和誤差的分布圖.數(shù)值實驗結(jié)果表明，RBF-PU 方法在兩種離散方式下在求解非局部擴散方程時都是有效的.

圖3 均勻離散無網(wǎng)格點的分布和非局部擴散方程矩陣結(jié)構(gòu)：(a) 無網(wǎng)格點分布和單位分解劃分； (b) 矩陣結(jié)構(gòu)Fig. 3 The distribution of points and the matrix structure of nonlocal diffusion under uniform discretization: (a) the distribution of meshless points and PU patches;(b) the matrix structure

圖4 Halton 離散下點的分布和非局部擴散方程矩陣結(jié)構(gòu)：(a) Halton 點的分布和單位分解劃分；(b) 矩陣結(jié)構(gòu)Fig. 4 The distribution of points and the matrix structure of nonlocal diffusion under Halton discretization: (a) the distribution of Halton points and circular PU patches； (b) the matrix structure

圖5 均勻離散下非局部擴散方程數(shù)值解和誤差分布：(a) 數(shù)值解； (b) 誤差分布Fig. 5 The numerical solution and the error distribution for the nonlocal diffusion equation under uniform discretization: (a) the numerical solution;(b) the error distribution

圖6 Halton 離散下非局部擴散方程的數(shù)值結(jié)果和誤差分布圖：(a) 數(shù)值解； (b) 誤差分布Fig. 6 The numerical solution and the error distribution for the nonlocal diffusion equation under Halton discretization: (a) the numerical solution;(b) the error distribution

表1 均勻離散非局部擴散方程數(shù)值結(jié)果(δ=0.2)Table 1 Numerical results of the nonlocal diffusion equation under uniform discretization (δ=0.2)

表2 離散非局部擴散方程的數(shù)值結(jié)果(δ=0.2)Table 2 Numerical results of the nonlocal diffusion equation under Halton discretization (δ=0.2)

3.2 近場動力學方程

算例2 取二維近場動力學問題(20)的精確解為

核函數(shù)為

在均勻離散的情況下對該問題進行了數(shù)值實驗，圖7 給出了均勻離散下點的分布和最終形成的矩陣結(jié)構(gòu)，計算給出了兩個方向的位移場在不同數(shù)目的離散點下位移的數(shù)值結(jié)果，位移場在x和y方向的數(shù)值結(jié)果分別如表3 和表4 所示.表中數(shù)據(jù)顯示，RBF-PU 方法在求解近場動力學方程的時間比相同離散的非局部擴散方程更低，原因在于在此算例中采取的單位分解劃分更小，導致局部插值矩陣的規(guī)模更小，節(jié)省了剛度矩陣生成的時間.與此同時，在50×50的離散點下位移和誤差的分布如圖8 和9 所示.計算結(jié)果表明，RBF-PU 方法在求解向量形式的近場動力學方程時也是有效的.

圖7 均勻離散下無網(wǎng)格點的分布和近場動力學方程矩陣結(jié)構(gòu)：(a) 無網(wǎng)格點分布和單位分解劃分； (b) 矩陣結(jié)構(gòu)Fig. 7 The distribution of points and the matrix structure for the peridynamic equation under uniform distribution: (a) the distribution of meshless points and circular PU patches; (b) the matrix structure

圖8 近場動力學方程x方向數(shù)值解和誤差分布：(a) 數(shù)值解； (b) 誤差分布Fig. 8 The numerical solution and the error distribution of the peridynamic equation in the x direction: (a) the numerical solution; (b) the error distribution

表3 模型的x 方向位移的數(shù)值解u 1(δ=0.2)Table 3 Numerical results of displacement in the x direction for the peridynamic model u 1 (δ=0.2)

表4 模型的y方向位移數(shù)值解u 2 (δ=0.2)Table 4 Numerical results of displacement in the y direction for the peridynamic model u 2 (δ=0.2)

圖9 近場動力學方程y方向數(shù)值解和誤差分布：(a) 數(shù)值解； (b) 誤差分布Fig. 9 The numerical solution and the error distribution of the peridynamic equation in the y direction: (a) the numerical solution; (b) the error distribution

4 結(jié) 論

本文發(fā)展了求解二維非局部擴散和近場動力學問題的RBF-PU 方法，該方法無網(wǎng)格依賴性、精度高、易于編程實現(xiàn).通過數(shù)值算例驗證了該方法對此類非局部方程求解的有效性.基于本文的工作，可以考慮將此類方法用于更加復雜的非線性、高維近場動力學模型的數(shù)值計算當中.

參考文獻（ References ） ：

［1］B OBARU F, MONCHAI D. The peridynamic formulation for transient heat conduction[J].International Journal of Heat and Mass Transfer, 2010, 53（19/20）: 4047-4059.

［2］D U Q, GUNZBURGER M, LEHOUCQ R B, et al. Analysis and approximation of nonlocal diffusion problems with volume constraints[J].SIAM Review, 2012, 54（4）: 667-696.

［3］S ILLING S A. Reformulation of elasticity theory for discontinuities and long-range forces[J].Journal of the Mechanics and Physics of Solids, 2000, 48（1）: 175-209.

［4］S ILLING S A, EPTON M, WECKER O, et al. Peridynamic states and constitutive modeling[J].Journal of Elasticity, 2007, 88（2）: 151-184.

［5］黃 丹, 章青, 喬丕忠, 等. 近場動力學方法及其應(yīng)用[J]. 力學進展, 2010, 40（4）: 61-70. （HUANG Dan, ZHANG Qin,QIAO Pizhong, et al. A review on peridynamics(PD) method and its application[J].Advances in Mechanics,2010, 40（4）: 61-70.（in Chinese））

［6］G U X, ZHANG Q, MADENCI E. Refined bond-based peridynamics for thermal diffusion[J].Engineering Computations, 2019, 36（8）: 2557-2587.

［7］李 天一, 章青, 夏曉舟, 等. 考慮混凝土材料非均質(zhì)特性的近場動力學模型[J]. 應(yīng)用數(shù)學和力學, 2018, 39（8）: 913-924. （ LI Tianyi, ZHANG Qing, XIA Xiaozhou, et al. A peridynamic model for heterogeneous concrete materials[J].Applied Mathematics and Mechanics, 2018, 39（8）: 913-924.（in Chinese））

［8］S ILLING S A, ASKARI E. A meshfree method based on the peridynamic model of solid mechanics[J].Computers and Structures, 2005, 83（17/18）: 1526-1535.

［9］C HEN X, GUNZBURGER M. Continuous and discontinuous finite element methods for a peridynamics model of mechanics[J].Computer Methods in Applied Mechanics and Engineering, 2011, 200（9/12）: 1237-1250.

［10］W ANG H, TIAN H. A fast Galerkin method with efficient matrix assembly and storage for a peridynamic model[J].Journal of Computational Physics, 2012, 231（23）: 7730-7738.

［11］T IAN H, WANG H, WANG W. An efficient collocation method for a nonlocal diffusion model[J].International Journal of Numerical Analysis and Modeling, 2013, 10（4）: 815-825.

［12］W ANG C, WANG H. A fast collocation method for a variable-coefficient nonlocal diffusion model[J].Journal of Computational Physics, 2017, 330: 114-126.

［13］W ANG H, TIAN H. A fast and faithful collocation method with efficient matrix assembly for a two-dimensional nonlocal diffusion model[J].Computer Methods in Applied Mechanics and Engineering, 2014, 273: 19-36.

［14］T IAN X C, DU Q. Analysis and comparison of different approximations to nonlocal diffusion and linear peridynamic equations[J].SIAM Journal on Numerical Analysis, 2013, 51（6）: 3458-3482.

［15］L EHOUCQ R B, ROWE S T. A radial basis function Galerkin method for inhomogeneous nonlocal diffusion[J].Computer Methods in Applied Mechanics and Engineering, 2016, 299: 366-380.

［16］Z HAO W, HON Y C, STOLL M. Localized radial basis functions-based pseudo-spectral method (LRBF-PSM) for nonlocal diffusion problems[J].Computers and Mathematics With Applications, 2017, 75（5）: 1685-1704.

［17］Z HAO W, HON Y C. An accurate and efficient numerical method for solving linear peridynamic models[J].Applied Mathematical Modelling, 2019, 74: 113-131.

［18］P ASETTO M, LENG Y, CHEN J, et al. A reproducing kernel enhanced approach for peridynamic solutions[J].Computer Methods in Applied Mechanics and Engineering, 2018, 340: 1044-1078.

［19］D U Q, YIN X B. A conforming DG method for linear nonlocal models with integrable kernels[J].Journal of Scientific Computing, 2019, 80（3）: 1913-1935.

［20］Z HANG S Y, NIE Y F. A POD-based fast algorithm for the nonlocal unsteady problems[J].International Journal of Numerical Analysis and Modeling, 2020, 17（6）: 858-871.

［21］L IANG X, WANG L J, XU J F, et al. The boundary element method of peridynamics[R/OL]. 2021（2021-06-16）[2021-10-29].https://arxiv.org/abs/2009.08008.

［22］L U J X, NIE Y F. A collocation method based on localized radial basis functions with reproducibility for nonlocal diffusion models[J].Computational and Applied Mathematics, 2021, 40（8）: 271-294.

［23］洪 文強, 徐績青, 許錫賓, 等. 求解Bratu型方程的徑向基函數(shù)逼近法[J]. 應(yīng)用數(shù)學和力學, 2016, 37（6）: 617-625.（HONG Wenqiang, XU Jiqing, XU Xibin, et al. The radial basis function approximation method for solving Bratu-type equations[J].Applied Mathematics and Mechanics, 2016, 37（6）: 617-625.（in Chinese））

［24］李 怡, 吳林鍵, 舒丹, 等. 基于Gauss全局徑向基函數(shù)的近岸淺水變形波高數(shù)值計算新方法[J]. 應(yīng)用數(shù)學和力學, 2014,35（8）: 903-912. （LI Yi, WU Linjian, SHU Dan, et al. A new method to calculate the wave height of deformed shallow water based on the gauss global radial basis function[J].Applied Mathematics and Mechanics, 2014,35（8）: 903-912.（in Chinese））

［25］W ENDLAND H.Scattered Data Approximation[M]. Cambridge: Cambridge University Press, 2005.

［26］S HEPARD D. A two-dimensional interpolation function for irregularly-spaced data[C]//Proceedings of the 1968 23rd ACM National Conference（ACM 68）. New York: Association for Computing Machinery, 1968: 517-524.