李志華　喻　軍　楊紅光

杭州電子科技大學，杭州，310018

偏微分方程與微分代數方程的一致求解方法

李志華喻軍楊紅光

杭州電子科技大學，杭州，310018

Modelica 語言是一種復雜物理系統(tǒng)多領域統(tǒng)一建模語言，但目前該語言只能解決由微分代數方程(DAE)描述的問題，而不能解決由偏微分方程(PDE)表達的問題。為此，提出一種偏微分方程與微分代數方程的一致求解方法，利用所構建的徑向基函數配點無網格法直接將偏微分方程在空間上離散成一系列的微分代數方程，然后采用成熟的微分代數方程求解器進行求解。實例結果表明，該方法在不改變 Modelica 語法的前提下，能較好地實現(xiàn)偏微分方程與微分代數方程的一致求解，且求解精度高、邊界條件處理簡單，有利于Modelica直接求解復雜工程系統(tǒng)中多領域耦合、時間域與空間域耦合的復雜問題。

多領域統(tǒng)一建模；Modelica；偏微分方程(PDE)；微分代數方程(DAE)

0　引言

現(xiàn)代復雜機電產品(如航空航天器、機器人、汽車等)通常是集機、電、液、控、磁等多學科領域于一體的復雜物理系統(tǒng)，經常表現(xiàn)出時間依賴(對時間導數)與空間依賴(對坐標偏導)共存的行為特征，而且可能呈現(xiàn)出多領域之間及時間域與空間域之間的耦合特性[1]。

物理系統(tǒng)行為規(guī)律的描述通常有兩種主要的形式。系統(tǒng)在單純時間域的物理行為往往由常微分方程(ordinary differential equation，ODE)描述，如果涉及代數約束，則形成微分代數方程(differential-algebraic equation，DAE)，DAE是描述時間域物理規(guī)律的普遍形式，如機械多體、電子電路等系統(tǒng)規(guī)律的描述；若物理行為涉及空間場，出現(xiàn)對空間變量的偏導，則往往由偏微分方程(partial differential equation，PDE)描述，如位勢、傳熱、波動等相關物理系統(tǒng)的規(guī)律描述[1]。

物理系統(tǒng)建模經歷了從面向過程建模到面向對象建模、連續(xù)域與離散域分散建模到連續(xù)-離散混合建模、單一領域獨立建模到多領域統(tǒng)一建模的發(fā)展階段。Modelica 語言是近年來歐洲仿真界為解決復雜物理系統(tǒng)建模與仿真問題而提出的一種多領域統(tǒng)一建模語言[2-3]，然而，目前Modelica語言只能對時間域的物理系統(tǒng)進行統(tǒng)一建模，還不能對空間域的物理系統(tǒng)進行描述，更無法對其進行仿真優(yōu)化，這大大限制了 Modelica 語言的應用范圍。為此，國外已有學者著手擴展 Modelica 語言以支持 PDE 問題的建模與仿真[4]。周凡利等[1]提出了解決該問題的思路，但沒有具體實現(xiàn)。李志華等[5-6]也開展了這方面的研究工作，初步實現(xiàn)了Modelica語言對PDE問題的描述、建模以及仿真求解。然而現(xiàn)有的這些方法都是采用簡單的有限差分法或線上法來求解具有規(guī)則區(qū)域(如矩形)的PDE問題，對于不規(guī)則區(qū)域的復雜PDE問題則無法求解。

本文在李志華等原有工作[5-6]的基礎上，從多領域統(tǒng)一建模與仿真的角度，針對一般性的PDE問題(包括不規(guī)則求解區(qū)域、復雜邊界條件、線性或非線性PDE系統(tǒng))，提出一種PDE與DAE的一致求解方法，為Modelica直接求解復雜工程系統(tǒng)中多領域耦合、時間域與空間域耦合的復雜問題奠定基礎。

1　PDE的已有解法

PDE的解法主要分為解析法和數值法。到目前為止，只有有限形式的PDE能夠得到解析解，在工程實際中一般采用數值求解。PDE的數值求解技術比較成熟，可用的方法包括有限差分法、有限元法、有限體積法以及線上法等[7]。但有限差分法和線上法只適用于規(guī)則的求解區(qū)域，而有限元法和有限體積法是基于網格的計算方法，在某些工程問題(如動態(tài)裂紋擴展、高速撞擊、沖擊破壞、流固耦合等)中存在網格的束縛，使得計算遇到很大的困難，因此出現(xiàn)了無網格方法[8]。

無網格方法只需要節(jié)點的信息，不需要節(jié)點與節(jié)點之間相互聯(lián)系的信息，這樣很容易在復雜計算區(qū)域內布置節(jié)點。無網格方法的構建主要包括近似函數的構建方法和微分方程的離散方法兩個部分，根據近似函數構建方法和微分方程離散方法的不同，可以構建出許多不同的無網格方法[8]。目前比較常用的近似函數的構建方法有：核函數近似方法、再生核近似方法、移動最小二乘近似方法、單位分解近似方法、徑向基函數近似方法、點插值近似方法等。微分方程的離散方法包括加權殘量法、配點法、Galerkin法以及局部Petrov-Galerkin法。

Khattak等[9]運用無網格配點法成功求解了一類非線性PDE；Yao等[10]應用全局和局部徑向基函數來求解三維拋物線型PDE，并比較了這兩種無網格方法的性能；Kamruzzaman等[11]利用多項式點插值和配點法來構造無網格方法，較好地求解了橢圓形、拋物線型和雙曲線型PDE；吳宗敏[12]介紹了散亂數據擬合研究中的徑向基函數方法,及其在散亂線性泛函信息插值、無網格PDE數值解中的應用；吳孝鈿[13]應用Sobolev splines徑向基函數和緊支柱正定徑向基函數，得到求解PDE邊值問題的無網格算法,并針對散亂數據的特點,給出了計算整體稠密度h的算法及如何通過加密節(jié)點使h值縮小的一個可行方法。然而上述無網格方法都是直接對時間變量和空間變量進行離散，只適合求解單純的PDE問題，不適合求解復雜的PDE與DAE耦合問題。

本文借鑒傳統(tǒng)的求解PDE的無網格方法，選用徑向基函數和配點法來構建徑向基函數配點無網格法，并對其進行改進，即只對空間變量離散而保持時間變量連續(xù)，直接將PDE在空間上(即配點處)離散成一系列的DAE，然后利用成熟的DAE求解器進行統(tǒng)一求解。

2　徑向基函數配點無網格法

對于d維實空間中定義在域Ω上的PDE問題：

(1)

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

我們所構建的徑向基函數配點無網格法的基本思想是：首先采用徑向基函數構造近似函數，將未知量場函數的時-空變量分開，然后運用配點法對空間變量進行離散，而保持時間變量連續(xù)，這樣就將PDE問題在空間上離散成一系列只含時間變量的DAE問題，具體過程如下。

首先在PDE的不規(guī)則求解區(qū)域Ω內和邊界?Ω上選定N=Nu+Nb個離散的節(jié)點(即配點)，然后應用徑向基函數構造u(X,t)的近似函數，并將其構成時間與空間分離的形式：

(2)

X∈Ω?Rd

其中，N為節(jié)點總數；Nu為域內節(jié)點數；Nb為邊界節(jié)點數；α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‖)對空間變量可求導，在每一節(jié)點處，空間坐標已知，因此Lφ(ri(Xk))或Bφ(ri(Xk))就是一個已知值，這樣，式(3)中就只剩下時間的導數，即每個節(jié)點處對應一個只與時間有關的DAE，這樣就將PDE轉化為一系列的DAE(具體過程可參見實例部分)。

進一步，PDE問題(式(1))可變?yōu)榍蠼庖粋€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語言的多領域統(tǒng)一建模與仿真平臺MWorks中[14]，并利用其成熟的DAE/ODE求解器進行PDE與DAE的一致求解。整個流程如圖1所示。

圖1　PDE與DAE的一致求解過程

以下面一個不規(guī)則區(qū)域的二維熱傳導問題為例來說明本文所提方法的有效性：

式中，T為某個時間值。

其求解區(qū)域?Ω由如下邊界組成：

y=0,0≤x≤1

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

y=1,0≤x≤1

0≤y≤1,x=0

對該不規(guī)則求解區(qū)域用配點法進行不規(guī)則離散，得其節(jié)點分布如圖2所示，其中域內節(jié)點數Nu=61，邊界上節(jié)點數Nb=42。然后對這些節(jié)點從左到右、從下到上進行編號。

圖2　不規(guī)則求解域的配點

對該二維熱傳導PDE問題，運用上述PDE與DAE的一致求解方法和過程進行一致求解，令

對于求解域內的Nu個離散節(jié)點(xi,yi)(i=1,2,…,Nu)，滿足域內的偏微分方程，即

(5)

對于求解域邊界上的Nb個離散節(jié)點(xi,yi)(i=1,2,…,Nb)，滿足邊界條件方程，即

(6)

對于域內及邊界上的所有離散節(jié)點(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]，在各離散節(jié)點處均可計算出它們的值。這樣，式(5)～式(7)中就只剩下時間的變量，即利用徑向基函數配點無網格法已將PDE在離散節(jié)點處轉化為一系列的DAE，然后就可以在MWorks環(huán)境中利用其自帶的DAE求解器進行求解，從而方便地得出場函數在各個離散節(jié)點處隨時間變化的函數值。圖3表示的是場函數在編號為15節(jié)點處的仿真結果；圖4所示為本文的數值解與其精確解u(x,y,t)=e-2tsin(x+y)之間的比較，此處t=0.5 s。

圖3　場函數在編號為15節(jié)點處的仿真結果

4　求解精度影響因素分析

為了更好地應用徑向基函數配點無網格法來求解PDE問題，本文研究了不同離散節(jié)點數、平均節(jié)點間距和徑向基函數參數c的取值對求解精度的影響，如表1所示。

由表1可以看出，一般情況下，當參數不變時，離散節(jié)點數越大、平均節(jié)點間距越小，則求解精度越高。例如，在c=1不變的情況下，當節(jié)點數為14、平均節(jié)點間距為0.47時，誤差為2.2475×

(a)本文的數值解

(b)精確解圖4　數值解與精確解的比較

離散節(jié)點數N平均節(jié)點間距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;而當節(jié)點數為30、平均節(jié)點間距為0.28時，誤差則為1.7142×10-5。同時還可以看出，隨著參數c的不斷增大，求解精度并不是呈遞增或遞減狀態(tài)，而是有起伏變化，只有當c取適當的值時，誤差er才較小。由此可見，恰當確定徑向基函數參數c的取值很關鍵，然而目前還沒有規(guī)律可循，只能通過多次反復運算來確定一個合適的值。

5　結論

多領域統(tǒng)一建模要求用偏微分方程和微分代數方程來統(tǒng)一描述和統(tǒng)一求解，本文針對一般性的偏微分方程問題，提出了偏微分方程與微分代數方程的一致求解方法，給出了該方法的實現(xiàn)過程，分析了離散節(jié)點數和徑向基函數參數對求解精度的影響，得到如下結論：

(1)與傳統(tǒng)的無網格方法相比，本文采用只對空間變量離散而保持時間變量連續(xù)的策略，能方便地將偏微分方程在離散節(jié)點處轉化為一系列微分代數方程，從而在不改變 Modelica 語法的前提下，較好地實現(xiàn)偏微分方程與微分代數方程的一致求解，大大簡化了復雜的偏微分方程與微分代數方程耦合問題的求解難度。

(2)實例結果表明，本文所提出的方法能較好地解決具有不規(guī)則求解區(qū)域的偏微分方程問題，且求解精度高，這有利于Modelica直接求解復雜工程系統(tǒng)中多領域耦合、時間域與空間域耦合的復雜問題。

[1]周凡利，陳立平，趙建軍，等. 時域-空間耦合物理系統(tǒng)多領域統(tǒng)一建模與仿真及偏微分代數混合方程系統(tǒng)的求解[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 語言及其多領域統(tǒng)一建模與仿真機理[J]. 系統(tǒng)仿真學報, 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]. 系統(tǒng)仿真學報，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]吳義忠，劉敏，陳立平. 多領域物理系統(tǒng)混合建模平臺開發(fā)[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.

(編輯蘇衛(wèi)國)

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年生。杭州電子科技大學機械工程學院碩士研究生。