椭圆型偏微分方程的三角形单元有限元的数值解法

俞 辉， 杨臻豪， 奉继明， 瞿少成

(1.三峡大学 理学院， 湖北 宜昌 443000; 2.华中师范大学 电子信息工程系， 武汉 430079)



椭圆型偏微分方程的三角形单元有限元的数值解法

俞 辉1*， 杨臻豪1， 奉继明1， 瞿少成2

(1.三峡大学 理学院， 湖北 宜昌 443000; 2.华中师范大学 电子信息工程系， 武汉 430079)

本文以椭圆型偏微分方程问题为背景研究三角形单元的有限元方法.随着计算机图形学的发展，三角形单元划分取得了巨大的成就，可以得到质量非常好的三角形单元，进而提高偏微分方程的有限元方法的数值解的精度.本文采用节点增量算法，对问题区域进行三角形单元划分，得到的三角形单元满足Delaunay条件，再对三角形单元的所有节点采用自适应编号，最后运用三角形单元的有限元方法得到椭圆型偏微分方程的数值解.通过数值实验，得出相比传统的三角形单元的有限元方法，本文的三角形单元的有限元方法减小了舍入误差，提高了计算精度.

Delaunay三角形单元； 自适应编号； 舍入误差

随着经济的发展，工程的各个领域都呈现蓬勃发展的趋势.但是在发展的同时遇到了很多的工程问题.很大一部分的工程问题通过数学建模[1]，可以转化成偏微分方程的问题.通过对流体力学的分析，应用数学知识，可以将流体力学中的问题转化成偏微分方程，即Navier-Stokes[2]方程，通过方程的分析和求解其数值解，足以描述边界层、涡流等多种现象.在基于股票的衍生证券市场上，对买入期权、卖出期权、股票的买入价格、卖出价格的分析，布莱尔(Black[3])和休尔斯(Scholes[3])得出了欧式期权的无套利价格的倒向偏微分方程.可以看出偏微分方程在很多领域都涉及到了.

有限元方法[4]是求解偏微分方程数值解的重要方法.有限元求解偏微分方程时，一个重要的步骤是[5]，运用数学知识将偏微分方程转化为积分方程.对偏微分方程的原问题区域进行单元划分[6-8]，如四边形单元划分和三角单元划分.20世纪80年代三角网格剖分在众多领域引起了广泛的关注.英国Bath大学数学分校的P.J Green和R. Sibson等从数学的角度对三角网格剖分进行了研究[9].Bath大学数学分校的A. Bowyer和澳大利亚悉尼大学Geology and Geophyics系的D.F.Watson于1981年发表了文章[10]，提出各不相同的三角网格构造方法，但最终得到的结果都是以Dirichlet/Voronoi图为理论基础的Delaunay三角网格.

本文将先进的三角形单元网格划分运用到有限元方法中.对问题区域采用成熟的节点增量算法[11]进行单元剖分，得到的三角形单元都是满足Delaunay条件的三角形.在运用有限元方法时，单元的“剖分”、“设置”、“编号”将影响有限元方程系数矩阵的形状，最后影响结果的计算精度.本文对问题区域进行三角形单元剖分后，对单元中的所有节点进行的编号，计算单元上系数矩阵和组装单元矩阵时，对单元的节点采用自适应编号排序，得到单元节点关系.最后从数值仿真例子可以得出，本文的有限元方法，比传统的有限元方法，计算精度更高.

1 有限元构思

考虑以下经典的二阶椭圆型方程的边值问题[12]：

(1)

其中，Ω是平面有界区域， Γ是区域的边界.

.

引理[5]设F是一个微分算子，则

(2)

通过上述引理(2)，可以将经典的椭圆型偏微分方程转化成积分方程了，即在边值问题中的偏微分方程两边分别同时乘以v，再对两边进行积分，从而将偏微分方程问题，转化如下成积分方程：

(3)

其中

从而将式(3)化简得到下面式子

(4)

记

于是得到如下式子

(5)

考虑到二阶椭圆型偏微分方程(1)的边界条件，即

其中，Γ是问题域的边界.所以可得到

于是椭圆型偏微分方程边值问题可以转化为

(6)

上式(6)称为椭圆型偏微分方程边值问题的积分形式.

从问题的积分形式可以看出，二阶椭圆型偏微分方程的条件是u的二阶偏导存在，通过运用Green格林公式，将u的二阶偏导问题转化为一阶偏导问题，这就使得原问题的解u的光滑性的要求大为下降，只要u和偏导u′x，u′y存在并且可积就可以了.将有限元方法应用于二阶椭圆型偏微分方程的边值问题，其本质是将有限元方法应用于积分形式的方程，故本文以后将放弃二阶椭圆型偏微分方程的边值问题的原形式，而直接从积分形式出发来研究.

设Pm为次数不超过m的多项式空间，则本文所要用到的多项式样条函数空间，可以写为

(7)

(8)

2 三角形单元的有限元方法实现

有限元方法中一个重要的思想是，将问题区域进行划分，得到一个一个单元区域，对单元区域应用有限元方法得到单元矩阵，通过单元区域之间的关系，组装单元矩阵得到线性方程组的系数矩阵.故三角形单元划分[14]包括2大部分，一是对问题区域进行三角形划分；二是建立单元与单元之间的关系即单元关系矩阵.

2.1三角形单元划分

本文仿真例子的问题区域是Ω={(x，y)|x2+y2=1}.取步长hx=hy=0.5，对圆形区域进行三角形单元划分.在三角形单元划分时，要避免畸形单元，即三角单元中存在角度非常小的单元.为了避免畸形单元，本文采用成熟的节点增量算法，使三角形单元满足Delaunay三角化特性[11].得图1.

图1 圆形区域三角形单元划分Fig.1 Circular area triangle unit division

2.2单元变换

图2 单元变换Fig.2 Unit transform

设

2.3单元插值函数

设

运用Matlab编程计算，可得：c1=2，c2=0，c3=0，c4=-1，c5=0，c6=0

Ni=2ξ2-ξ.

同理可得

Nj=2η2-η，

Nk=2ξ2+2η2+4ξη-3ξ-3η+1，

Nr=4ξη，

Nq=4ξ-4ξη-4ξ2， Np=4η-4ξη-4η2.

至此就建立了标准三角单元上Lagrange二次型的6个基函数了，故Lagrange二次型插值函数如下

Npup+Nquq+Nrur.

2.4有限元方法的实现

在建立有限元方程[12]时，用的是Lagrange二次型插值函数，直接从微分方程的积分方程出发，积分方程的形式如下

因为Lagrange二次型插值函数为

Npup+Nquq+Nrur,

Npvp+Nqvq+Nrvr,

将插值函数代入得

aijuivj+aikuivk+aipuivp+…+

ariurvi+arjurvj+arkurvk+arpurvp+

μijuivj+μiruivr+…+μriurvi+

djvj+dkvk+dpvp+dqvq+drvr],

其中,γij是边界Γh的一条边，长度为lij，并且在γij上，ξ+η=1.

q·NmNl]Jdξdη,m，l=i，j，k，p，q，r,

令

故

上式即为三角单元的有限元方程，通过解线性方程就可以求解出椭圆型偏微分方程有限元的数值解了.

2.5单元自适应编号排序

下图3和图4是2种单元变换.

图3 单元变换1Fig.3 1 transform unit

图4 单元变换2Fig.4 2 transform unit

故Delaunay三角单元有限元求解椭圆型偏微分方程的流程图如图5.

图5 有限元方法求解流程图Fig.5 Finite element method for solving the flow chart

3 数值例子

本文有限元方法采用的是Lagrange二次元插值型，在计算数值积分时可用到了高斯积分技术[1]，所有的程序是在Matlab 2012编程实现的.本文单元划分步长取的是0.1.通过Matlab程序运行，得到真实值和半有限元法解如图6.

图6 有限元求解对比图Fig.6 Finite element solution comparison chart

上图5中，红色线表示真实解，绿色线表示有限元数值解.左边图，在求有限元数值解时，没有进行Delaunay三角划分和单元自适应编号处理.右边图，在求有限元数值解时，单元划分时，对三角单元进行了Delaunay三角化处理，并且对单元节点采用了自适应编号排序处理.从上面两图比较可以看出，通过对三角单元进行了Delaunay三角化处理和自适应编号处理，可以减少舍入误差，计算结果的精度更高.从右图可以得出，红色线和绿色线十分吻合，可以说明本文有限元方法得到的数值解，逼近真实解，即本文的有限元方法是可行的.

4 结论

有限元方法作用于微分方程，其本质是通过数学方法，将微分方程转化成积分方程，对问题区域进行单元划分，建立单元插值函数，通过在各个单元上运用有限元方法，应用高斯数值积分求取单元上的系数矩阵，通过单元关系矩阵将单元上的系数矩阵，组装成总系数矩阵，求解有限元方程，得到偏微分方程的有限元方法的数值解.本文将Delaunay三角化运用到三角单元划分，并对单元节点采用自适应编号排序，减少了舍入误差，从而提高了数值解的精度.

[1] 倪 兴. 常微分方程数值解法及其应用[D]. 合肥：中国科学技术大学, 2010．

[2] GIRAULT V, RAVIART P A. Finite element method for Navier-Stokes equations: theory and algorithms[M]. Berlin: Heidelber Springer-Verlag, 1981．

[3] BLACK F, SCHOLES M. The pricing of options and corporate liabilities[ J] . Journal of Political Economy , 1973 , 81(4): 633-654．

[4] Ho-Le K. Finite element mesh generation methods: a review and classification[J]. Computer Aided Design, 1988, 20(1):27-38．

[5] 傅永华.有限元分析基础[M]. 武汉:武汉大学出版社, 2003．

[6] Zienkeiwicz O C. The finite element method[M]. 3rd edition, McGraw-Hill, 1977．

[7] COGGON J H. Electromagnetic and electrical modeling by the finite element method[J]. Geophysics, 1971, 36(2):132-155．

[8] MILLER G L, TALMOR D, TENG S H, et al. Control valume meshes using sphere packing: generation, refinement and coarsening[J]. Proceeding of 5th International Meshing Roundtable, 1996, 47-61．

[9] Green P J, Sibson R. Computing Dirichlet tessellation in the plane[J]. The Computer Journal, 1987, 2(2):168-173．

[10] Bowyer A. Computing Dirichlet Tessellations[J]. The Computer Journal, 1981, 24(2):162-166．

[11] 杨 钦. 限定Delaunay三角网络剖分技术[M].北京: 电子工业出版社， 2005.

[12] 林 群. 微分方程数值解法基础教程[M].第二版.北京:科学出版社， 2003．

[13] 华东师范大学数学系. 数学分析[M]. 北京: 高等教育出版社, 2001．

[14] 罗智中. 一种任意二维区域图形的三角剖分方法[J].机床与液压, 2002(6):145-146．

Triangle finite element numerical solution of elliptic partial differential equations

YU Hui1, YANG Zhenhao1, FENG Jiming1, QU Shaocheng2

(1.College of Science，China Three Gorges University，Yichang, Hubei 443000;2.Department of Electronics and Information Engineering，Central China Normal University，Wuhan 430079)

In the present work, the finite element method for triangular element is studied based on the elliptic partial differential equation. With the development of computer graphics, the triangle element division has made great achievements. High-quality triangular elements are able to be generated, which improved the accuracy of numerical solutions for partial differential equations. Here, the triangular element satisfying the Delaunay conditions is calculated through unit division of the region with the node incremental algorithm. All nodes of the triangular element are adopted with adaptive numbers and the numerical solution of elliptic partial differential equations is obtained utilizing the finite element method of triangular elements. Finally, through numerical experiments, the finite element method is compared with the traditional one. Results indicate that the finite element method of triangular element reduces the rounding error and elevates the calculation precision.

delaunay triangular element; adaptive number; rounding error

2016-03-16.

国家自然科学基金项目(61273183, 61374028, 61304162).

1000-1190(2016)04-0489-07

O175.25

A

*E-mail: yuhui@ctgu.edu.cn.