物理信息依赖核函数配点法的研究进展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]进行求解,目前这方面的研究工作还较少,有待进一步研究.此外,本文从齐次微分方程物理信息依赖核函数开始,依次介绍了非齐次、非均质、非稳态微分方程构造物理信息依赖核函数的方法,最后根据物理力学问题统计意义下的概率密度函数构建物理信息依赖核函数,无需建立微分方程问题的显式表达式.更进一步地,是否可以直接根据已有数据特点构建物理信息依赖核函数,进而计算得到物理力学问题的解？这也有待进一步研究.