李志华　喻　军　杨红光

杭州电子科技大学，杭州，310018

偏微分方程与微分代数方程的一致求解方法

李志华喻军杨红光

杭州电子科技大学，杭州，310018

Modelica 语言是一种复杂物理系统多领域统一建模语言，但目前该语言只能解决由微分代数方程(DAE)描述的问题，而不能解决由偏微分方程(PDE)表达的问题。为此，提出一种偏微分方程与微分代数方程的一致求解方法，利用所构建的径向基函数配点无网格法直接将偏微分方程在空间上离散成一系列的微分代数方程，然后采用成熟的微分代数方程求解器进行求解。实例结果表明，该方法在不改变 Modelica 语法的前提下，能较好地实现偏微分方程与微分代数方程的一致求解，且求解精度高、边界条件处理简单，有利于Modelica直接求解复杂工程系统中多领域耦合、时间域与空间域耦合的复杂问题。

多领域统一建模；Modelica；偏微分方程(PDE)；微分代数方程(DAE)

0　引言

现代复杂机电产品(如航空航天器、机器人、汽车等)通常是集机、电、液、控、磁等多学科领域于一体的复杂物理系统，经常表现出时间依赖(对时间导数)与空间依赖(对坐标偏导)共存的行为特征，而且可能呈现出多领域之间及时间域与空间域之间的耦合特性[1]。

物理系统行为规律的描述通常有两种主要的形式。系统在单纯时间域的物理行为往往由常微分方程(ordinary differential equation，ODE)描述，如果涉及代数约束，则形成微分代数方程(differential-algebraic equation，DAE)，DAE是描述时间域物理规律的普遍形式，如机械多体、电子电路等系统规律的描述；若物理行为涉及空间场，出现对空间变量的偏导，则往往由偏微分方程(partial differential equation，PDE)描述，如位势、传热、波动等相关物理系统的规律描述[1]。

物理系统建模经历了从面向过程建模到面向对象建模、连续域与离散域分散建模到连续-离散混合建模、单一领域独立建模到多领域统一建模的发展阶段。Modelica 语言是近年来欧洲仿真界为解决复杂物理系统建模与仿真问题而提出的一种多领域统一建模语言[2-3]，然而，目前Modelica语言只能对时间域的物理系统进行统一建模，还不能对空间域的物理系统进行描述，更无法对其进行仿真优化，这大大限制了 Modelica 语言的应用范围。为此，国外已有学者着手扩展 Modelica 语言以支持 PDE 问题的建模与仿真[4]。周凡利等[1]提出了解决该问题的思路，但没有具体实现。李志华等[5-6]也开展了这方面的研究工作，初步实现了Modelica语言对PDE问题的描述、建模以及仿真求解。然而现有的这些方法都是采用简单的有限差分法或线上法来求解具有规则区域(如矩形)的PDE问题，对于不规则区域的复杂PDE问题则无法求解。

本文在李志华等原有工作[5-6]的基础上，从多领域统一建模与仿真的角度，针对一般性的PDE问题(包括不规则求解区域、复杂边界条件、线性或非线性PDE系统)，提出一种PDE与DAE的一致求解方法，为Modelica直接求解复杂工程系统中多领域耦合、时间域与空间域耦合的复杂问题奠定基础。

1　PDE的已有解法

PDE的解法主要分为解析法和数值法。到目前为止，只有有限形式的PDE能够得到解析解，在工程实际中一般采用数值求解。PDE的数值求解技术比较成熟，可用的方法包括有限差分法、有限元法、有限体积法以及线上法等[7]。但有限差分法和线上法只适用于规则的求解区域，而有限元法和有限体积法是基于网格的计算方法，在某些工程问题(如动态裂纹扩展、高速撞击、冲击破坏、流固耦合等)中存在网格的束缚，使得计算遇到很大的困难，因此出现了无网格方法[8]。

无网格方法只需要节点的信息，不需要节点与节点之间相互联系的信息，这样很容易在复杂计算区域内布置节点。无网格方法的构建主要包括近似函数的构建方法和微分方程的离散方法两个部分，根据近似函数构建方法和微分方程离散方法的不同，可以构建出许多不同的无网格方法[8]。目前比较常用的近似函数的构建方法有：核函数近似方法、再生核近似方法、移动最小二乘近似方法、单位分解近似方法、径向基函数近似方法、点插值近似方法等。微分方程的离散方法包括加权残量法、配点法、Galerkin法以及局部Petrov-Galerkin法。

Khattak等[9]运用无网格配点法成功求解了一类非线性PDE；Yao等[10]应用全局和局部径向基函数来求解三维抛物线型PDE，并比较了这两种无网格方法的性能；Kamruzzaman等[11]利用多项式点插值和配点法来构造无网格方法，较好地求解了椭圆形、抛物线型和双曲线型PDE；吴宗敏[12]介绍了散乱数据拟合研究中的径向基函数方法,及其在散乱线性泛函信息插值、无网格PDE数值解中的应用；吴孝钿[13]应用Sobolev splines径向基函数和紧支柱正定径向基函数，得到求解PDE边值问题的无网格算法,并针对散乱数据的特点,给出了计算整体稠密度h的算法及如何通过加密节点使h值缩小的一个可行方法。然而上述无网格方法都是直接对时间变量和空间变量进行离散，只适合求解单纯的PDE问题，不适合求解复杂的PDE与DAE耦合问题。

本文借鉴传统的求解PDE的无网格方法，选用径向基函数和配点法来构建径向基函数配点无网格法，并对其进行改进，即只对空间变量离散而保持时间变量连续，直接将PDE在空间上(即配点处)离散成一系列的DAE，然后利用成熟的DAE求解器进行统一求解。

2　径向基函数配点无网格法

对于d维实空间中定义在域Ω上的PDE问题：

(1)

式中，L、B为微分算子；u(X,t)为未知量场函数；f(X,t)、g(X,t)为已知函数。

我们所构建的径向基函数配点无网格法的基本思想是：首先采用径向基函数构造近似函数，将未知量场函数的时-空变量分开，然后运用配点法对空间变量进行离散，而保持时间变量连续，这样就将PDE问题在空间上离散成一系列只含时间变量的DAE问题，具体过程如下。

首先在PDE的不规则求解区域Ω内和边界∂Ω上选定N=Nu+Nb个离散的节点(即配点)，然后应用径向基函数构造u(X,t)的近似函数，并将其构成时间与空间分离的形式：

(2)

X∈Ω⊆Rd

其中，N为节点总数；Nu为域内节点数；Nb为边界节点数；αi(t)为待定系数；Xi为实空间上的配点；φi(‖X-Xi‖)为径向基函数，φi(‖X-Xi‖)=φ(ri(X))；ri(X)为Euclidian范数，ri(X)=‖X-Xi‖。

为确保解的唯一性，φ(ri(X))必须无条件正定，这种径向基函数包括Gaussian函数e-c r2、逆MQ函数(r2+c2)β(β<0)和紧支正定径向基函数等。本文采用Gaussian函数。

将式(2)代入式(1)中，并使这N个点满足式(1)的微分方程和边界条件，得到

(3)

当微分算子L、B为空间变量的偏导时，由于空间变量已与时间变量分离，且径向基函数φi(‖X-Xi‖)对空间变量可求导，在每一节点处，空间坐标已知，因此Lφ(ri(Xk))或Bφ(ri(Xk))就是一个已知值，这样，式(3)中就只剩下时间的导数，即每个节点处对应一个只与时间有关的DAE，这样就将PDE转化为一系列的DAE(具体过程可参见实例部分)。

进一步，PDE问题(式(1))可变为求解一个N×N的线性方程组问题，用矩阵表示为

Sa=b

(4)

其中，S为N×N的矩阵，S=(Ski)N×N；a为待求的系数向量，a=(α1(t),α2(t),…,αN(t));b=(b1,b2,…,bN);且

k,i=1,2,…,N

由式(4)求出未知系数αi(t)(i=1,2,…,N)后，通过式(2)就可以获得域内和边界上任意一点的场函数值u(X,t)。

3　实例及求解过程

通过编程，利用径向基函数配点无网格法对PDE进行空间离散，自动将PDE问题转化成DAE问题。本文采用MATLAB编程，首先将带时间域的PDE问题进行时间与空间分离(若不带时间域则不用分离)，然后通过径向基函数配点无网格法对空间变量进行离散，得到一系列离散点处的DAE，并把这些DAE数据用mat格式保存起来，接着将该mat格式文件导入到基于Modelica语言的多领域统一建模与仿真平台MWorks中[14]，并利用其成熟的DAE/ODE求解器进行PDE与DAE的一致求解。整个流程如图1所示。

图1　PDE与DAE的一致求解过程

以下面一个不规则区域的二维热传导问题为例来说明本文所提方法的有效性：

式中，T为某个时间值。

其求解区域∂Ω由如下边界组成：

y=0,0≤x≤1

0≤y≤1,x=0.25cos2πy+0.75

y=1,0≤x≤1

0≤y≤1,x=0

对该不规则求解区域用配点法进行不规则离散，得其节点分布如图2所示，其中域内节点数Nu=61，边界上节点数Nb=42。然后对这些节点从左到右、从下到上进行编号。

图2　不规则求解域的配点

对该二维热传导PDE问题，运用上述PDE与DAE的一致求解方法和过程进行一致求解，令

对于求解域内的Nu个离散节点(xi,yi)(i=1,2,…,Nu)，满足域内的偏微分方程，即

(5)

对于求解域边界上的Nb个离散节点(xi,yi)(i=1,2,…,Nb)，满足边界条件方程，即

(6)

对于域内及边界上的所有离散节点(xi,yi(i=1,2,…,N)，N=Nu+Nb，满足初始条件方程，即

(7)

本文采用Gaussian径向基函数φ(r)=e-c r2，并取c=7。因此，φj(ri)=e-c[(xi-xj)2+(yi-yj)2]，φ″j(ri)=4ce-c[(xi-xj)2+(yi-yj)2][-1+(xi-xj)2+(yi-yj)2]，在各离散节点处均可计算出它们的值。这样，式(5)～式(7)中就只剩下时间的变量，即利用径向基函数配点无网格法已将PDE在离散节点处转化为一系列的DAE，然后就可以在MWorks环境中利用其自带的DAE求解器进行求解，从而方便地得出场函数在各个离散节点处随时间变化的函数值。图3表示的是场函数在编号为15节点处的仿真结果；图4所示为本文的数值解与其精确解u(x,y,t)=e-2tsin(x+y)之间的比较，此处t=0.5 s。

图3　场函数在编号为15节点处的仿真结果

4　求解精度影响因素分析

为了更好地应用径向基函数配点无网格法来求解PDE问题，本文研究了不同离散节点数、平均节点间距和径向基函数参数c的取值对求解精度的影响，如表1所示。

由表1可以看出，一般情况下，当参数不变时，离散节点数越大、平均节点间距越小，则求解精度越高。例如，在c=1不变的情况下，当节点数为14、平均节点间距为0.47时，误差为2.2475×

(a)本文的数值解

(b)精确解图4　数值解与精确解的比较

离散节点数N平均节点间距h参数c误差er140.4712.2475×10-41.51.1082×10-428.1947×10-52.53.4945×10-437.0821×10-441.8×10-3300.280.57.544×10-50.73.0016×10-511.7142×10-51.53.12×10-521.0034×10-42.52.3097×10-4650.182.757.0608×10-531.2497×10-43.16.0971×10-53.28.8271×10-541.1196×10-451.4230×10-41030.142.754.9346×10-56.63.2072×10-477.1×10-57.51.5381×10-481.5819×10-49.52.6601×10-4

10-4;而当节点数为30、平均节点间距为0.28时，误差则为1.7142×10-5。同时还可以看出，随着参数c的不断增大，求解精度并不是呈递增或递减状态，而是有起伏变化，只有当c取适当的值时，误差er才较小。由此可见，恰当确定径向基函数参数c的取值很关键，然而目前还没有规律可循，只能通过多次反复运算来确定一个合适的值。

5　结论

多领域统一建模要求用偏微分方程和微分代数方程来统一描述和统一求解，本文针对一般性的偏微分方程问题，提出了偏微分方程与微分代数方程的一致求解方法，给出了该方法的实现过程，分析了离散节点数和径向基函数参数对求解精度的影响，得到如下结论：

(1)与传统的无网格方法相比，本文采用只对空间变量离散而保持时间变量连续的策略，能方便地将偏微分方程在离散节点处转化为一系列微分代数方程，从而在不改变 Modelica 语法的前提下，较好地实现偏微分方程与微分代数方程的一致求解，大大简化了复杂的偏微分方程与微分代数方程耦合问题的求解难度。

(2)实例结果表明，本文所提出的方法能较好地解决具有不规则求解区域的偏微分方程问题，且求解精度高，这有利于Modelica直接求解复杂工程系统中多领域耦合、时间域与空间域耦合的复杂问题。

[1]周凡利，陈立平，赵建军，等. 时域-空间耦合物理系统多领域统一建模与仿真及偏微分代数混合方程系统的求解[C]//中国力学学会学术大会.北京：中国力学学会，2007：760-760.

[2]Fritzson P.Introduction to Modeling and Simulation of Technical and Physical Systems with Modelica[M]. New York: Wiley-IEEE Press, 2011.

[3]赵建军，丁建完，周凡利，等. Modelica 语言及其多领域统一建模与仿真机理[J]. 系统仿真学报, 2006, 18(增刊2): 570-573.

Zhao Jianjun, Ding Jianwan, Zhou Fanli, et al. Modelica and Its Mechanism of Multi-domain Unified Modeling and Simulation[J]. Journal of System Simulation, 2006, 18(S2): 570-573.

[4]Saldamli L, Bachmann B, Wiesmann H, et al. A Framework for Describing and Solving PDE Models in Modelica[C]//Proceedings of the 4th International Modelica Conference. Hamburg: Hamburg University of Technology, 2005:113-122.

[5]李志华，张慧丽，郑玲. 基于Modelica语言的PDE与DAE问题的一致表示[J]. 系统仿真学报，2009，21(15):4641-4646.

Li Zhihua, Zhang Huili, Zheng Ling. Consistent Representation of PDE and DAE Problems in Modelica[J]. Journal of System Simulation, 2009，21(15):4641-4646.

[6]Li Zhihua, Zheng Ling, Zhang Huili.Modelling and Simulation of PDE Problems in Modelica[J]. International Journal of Materials and Structural Integrity, 2009, 3(4):318-331.

[7]李荣华. 偏微分方程数值解法[M]. 2版. 北京：高等教育出版社，2010.

[8]刘欣. 无网格方法[M]. 北京：科学出版社，2011.

[9]Khattak A J, Tirmizi S, Islam S. Application of Meshfree Collocation Method to a Class of Nonlinear Partial Differential Equations[J]. Engineering Analysis with Boundary Elements, 2009, 33(5):661-667.

[10]Yao Guangming, Siraj U I. Assessment of Global and Local Meshless Methods Based on Collocation with Radial Basis Functions for Parabolic Partial Differential Equations in Three Dimensions[J]. Engineering Analysis with Boundary Elements, 2012, 36(11):1640-1648.

[11]Kamruzzaman M, Sonar T, Lutz T, et al. A New Meshless Collocation Method for Partial Differential Equations[J]. Communications in Numerical Methods in Engineering, 2008, 24(12):1617-1639.

[12]吴宗敏. 径向基函数、散乱数据拟合与无网格偏微分方程数值解[J]. 工程数学学报, 2002, 19(2):1-12.

Wu Zongmin. Radial Basis Function Scattered Data Interpolation and the Meshless Method of Numerical Solution of PDEs[J]. Journal of Engineering Mathematics, 2002, 19(2):1-12.

[13]吴孝钿. 求解偏微分方程的一类无网格算法[J]. 复旦大学学报(自然科学版), 2004, 43(3):292-299.

Wu Xiaotian. Meshless Method of Solving Partial Differential Equations[J]. Journal of Fudan University (Natural Science), 2004, 43(3):292-299.

[14]吴义忠，刘敏，陈立平. 多领域物理系统混合建模平台开发[J]. 计算机辅助设计与图形学学报，2006, 18(1):120-124.

Wu Yizhong, Liu Min, Chen Liping. Development of Hybrid Modeling Platform for Multi-domain Physical System[J]. Journal of Computer-aided Design & Computer Graphics, 2006, 18(1):120-124.

(编辑苏卫国)

Consistent Solving Method of PDE and DAE

Li ZhihuaYu JunYang Hongguang

Hangzhou Dianzi University,Hangzhou,310018

Modelica is a multi-domain unified modeling language for modeling and simulation of large and complex physical systems. However, it dealt only with DAE but not with PDE. A consistent solving method of PDE and DAE was proposed. The PDE was transferred into a series of DAEs with the meshless method of radial basis function collocation, and was solved by the mature DAE solver in MWorks platform based on Modelica. Results show that this consistent solving method realizes the consistent solution of PDE and DAE under the premise of not changing Modelica grammar, and has high accuracy and the convenience of dealing with boundary conditions, which is conducive to solve complex engineering systems with multi-domain coupling and time domain and space domain coupling.

multi-domain unified modeling; Modelica; partial differential equation (PDE); differential-algebraic equation (DAE)

2013-06-19

2014-12-22

国家自然科学基金资助项目(51275141)；浙江省自然科学基金资助项目(Y1100901)

TH122；TP391DOI：10.3969/j.issn.1004-132X.2015.04.003

李志华，男，1966年生。杭州电子科技大学机械工程学院教授、博士。主要研究方向为多领域建模与仿真、CAD/CAE等。喻军，男，1989年生。杭州电子科技大学机械工程学院硕士研究生。杨红光，男，1985年生。杭州电子科技大学机械工程学院硕士研究生。