PDE工具箱实现偏微分方程的有限元求解

1、 基金项目 :国家 “ 十五 ” 863课题 (2001AA602021 资助 收稿日期 :2005-05-15 第 23卷 第 9期计 算 机 仿 真2006年 9月 文章编号 :1006-9348(2006 09-0100-03P D E 工具箱实现偏微分方程的有限元求解黄作英 , 阙沛文 , 陈亮(上海交通大学电子信息学院 , 上海 200030摘要 :偏微分方程在科学和工程上有着广泛的应用。 有限元法是一种重要的偏微分方程数值解法。 编程实现从偏微分方程 到有限元求解全过程需要很好的理论基础和编程技巧 , 难度较高。该文介绍了偏微分方程有限元求解的基本理论和一般Neumann 条件下椭

2、圆型方程的有限元求解具体过程 , 并通过两个实例 , , 介绍了使用 P DE 工具箱实现偏微分方程的有限元解法。 实验结果表明这一方法具有操作简单明了 , 运算速度快 , 关键词 :偏微分方程 ; 有限元法 ; 数值解法 中图分类号 :TP391. 9 文献标识码 :ASoluti on of Parti a l ti ti by FE M i n PD E Toolboxying, QUE Pei -wen, CHEN L iang(and I nfor mati on, Shanghai J iaot ong University, Shanghai 200030, China ABST

3、RACT:differential equati on is widely used in p r oble m s of science and engineering . Finite ele ment method (FE M is an i m portant nu merical s oluti on method of partial differential equati on . Good theory foundati on and p r ogra mm ing skill are necessary f or p r ogra mm ing t o s olve part

4、ial differential equati on by FE M method . This article in 2tr oduces the basic theory of FE M method f or s olving partial differential equati on and the s oluti on p r ocedure t o elli p tic partial differential equati on under Neu mann boundary conditi on . T wo exa mp les, magnetic p r oble m o

5、f mot or and heat exchange p r oblem, are s olved in P DE T oolbox . The results show that thismethod is easy t o handle and easy t o contr ol the err or with faster operati on s peed .KE YWO RD S:Partial differential equati on; Finite ele ment method (FE M ; Nu merical s oluti on method1 引言偏微分方程定解问

6、题有着广泛的应用背景 , 科学和工程 中的大多数问题在数学上是由偏微分方程描述的 , 很多重要 的物理力学学科的基本方程本身就是偏微分方程 , 人们用偏 微分方程来描述 , 解释或者预见各种自然现象并用于科学和 技术的各个领域 , 随着计算机技术的飞速发展 , 数值方法已 经成为一种越来越重要的工具 。编制高效的程序求解各种 偏微分方程已成为可能 , 偏微分方程数值解法也成为科学或 工程计算的重要分支 。对于理论研究和实际应用中提出的偏微分方程问题 , 由 于其边界和边界条件复杂等原因 , 寻求解析解是非常困难甚 至是不可能的 , 利用计算机研究相应问题的数值解法是十分 必要的 。 编程实现从

7、偏微分方程数值求解全过程需要很好 的理论基础和编程技巧 , 以有限元法为例 , 求解要经过等价 解的导出 , 区域划分及多项式选取 , 单元刚度矩阵 , 单元载荷 向量的形成 , 总刚度矩阵和总载荷向量的组装 , 求解等多个步骤 。 本文介绍了利用 P DE 工具箱对偏微分方程的有限元 求解 , 借助于这个工具 , 我们可以从繁琐 , 共性的求解步骤中 解脱出来而专注于问题的核心即问题的描述 , 定义及简化 , 边界条件的确定 , 求解方法和精度控制的选择等 , 从而大大 提高了求解效率 。2 偏微分方程的有限元解法2. 1 问题偏微分方程可根据它们的数学特征分为三大类型即椭 圆型方程 , 抛

8、物线型方程和双曲型方程 。 它们表示如下 :椭圆型方程 :- (c u +au =f (在 域中 (1a 是平面有界区域 , c, a, f 以及未知函数 u 是定义在 域上的实或复函数 ;抛物型方程 :d5t- (c u +au =f (在 域中 (1b 双曲型方程 :001d 25t2- (c u +au =f (在 域中 (1c 其中 d 是定义在 上的复函数 , 是待求的特征值 , 在抛 物型和双曲型方程中 , 系数 c, a, f 和 d 可以依赖于时间 t 。边界条件 :D irichlet 条件 :在边界上 hu =r (2a Neumann 条件 :在边界上n (c u +qu

9、 =g(2b 混合边界条件 :是前述两种条件的组合2. 2数值解法对于在理论研究和实际应用问题中提出的许多偏微分 方程 , 由于其边界和边界条件复杂等原因 , 寻求解析解是非 常困难甚至是不可能的 , 所以必须利用计算机来研究偏微分 方程的数值解 , 常用的偏微分方程数值分析方法包括有限差 分法和有限元法 。 有限差分法是以差分原理为基础的一种数 值计算方法 , 即用差分方程代替偏微分方程 , 值问题转化为一组相应的差分问题 , , 方法简单 , , 网格划 分缺乏灵活性 。 似解的方法 。 先从偏微分方程边值问题出发 , 找出一个能量 泛函的积分式 , 并令其在满足第一类边界条件的前提下取极

10、 值 , 即构成条件变分问题 。 这个条件变分问题是和偏微分方 程边值问题等价的 。 之后以条件变分问题为对象来求解偏微 分方程问题 。 在求解过程中 , 将场的求解区域剖分成有限个 单元 , 因此在单元中构造出插值函数 , 将插值函数代入能量 泛函的积分式 , 再把泛函离散化成多元函数 。 通过多元函数 极值求极值方法得到一个代数方程组 。 最后由此方程组求解 得到数值解 。 有限元法由于单元定义灵活 , 处理边界条件容 易 , 具有正定对称系数矩阵而占据主导地位 。下面以一般 Neumann 条件下椭圆型方程为例说明有限 元求解的具体过程如下 :1 导出定解问题对应的弱形式取任意试验函数

11、V, 乘 (1a 式两边 , 并在 上积分-( (c u +audx =f dx (3应用格林公式并加 Neumann 边界条件可得(c u +au -f dx -5(-qu +g ds =0(4式 (4 也被称为式 (1a 的广义解或弱解 , 在一定限制条 件下问题 (1a 和它的弱解形式 (4 是等价的 。2 区域划分和插值多项式的选取区域划分是将连续场域离散为有限数目元素的过程 , 目 前这一过程一般由计算机程序自动生成并能实现自适应网 格划分 , 多种基本单元已被用于有限元分析 , 采用三角形单 元进行划分在几何上具有很大的灵活性对边界 5逼近较 好 , 考虑三角形单元的线性插值函数u

12、 (x, y =ax +by +c (5设三角形三个顶点的 u 的值分别为 u i , u j , u m 可得方程 组ax i +by i +c =u i ax j +by j +c =u j ax m +by m +c =u m(6 求解 (5 (6 式可得u (x, y =N i (x, y u i +N j (x, y u j +N m (x, y u m (7 其中 N i (x, y , N j (x, y , N m (x, y 为顶点坐标的函数(7 式也可以表示为 u =N(8 u 的梯度向量u =(5u /5x , 5y T =Be(9其中 e =Ni , N j , N m

13、 T , B是2。, i (i =1, 2,. . . , N p 为 V 的 N p 维子空间的基函 数 , 按所剖分的单元将 (4 式改写为 N En =1en(c u +au -f dxdy= N En =1rn(g -qu ds (10这里 N E 是单元数 , r n =5e n 5(10 式的积分都只是在某个单元上求积 , 下面讨论每个单元上的积分计算 ,ec u +audxdy=3Te k e e(11ef dxdy =e(N3eTfdxdy=3Te Fe =F i(12 rqu ds =rq (N3eT(N3e ds=3Te K e e(13 rgds = r( T(g ds

14、=r(N3eT gds=3Te F e(14如果 e 单元是一边界单元 , 则单元刚度矩阵应为 K e + K e , 单元载荷向量应为 Fe + F e4 总刚度矩阵和总载荷向量的组装将单元刚度矩阵和单元载荷向量表达式代入 (10 , 为了 便于叠加 , 先对 e , Fe , Ke 进行扩充 , 扩充为 N p 维向 量和 N p ×N p 阶方阵 , 然后对号入座进行叠加 。5 约束处理 , 求解方程组对于 Neu mann 条件 , 由于是自然边界条件 , 5上不需要 满足任何约束条件 , 因此 , 可立即得到线性方程组 :K U =F(15K 为总刚度矩阵 , U 为位移向

15、量 , F 为总载荷向量 , 如果是 D irichlet 条件 , 还需对边界节点进行约束处理 , 然后再解101 方程组 。 3 计算方法和实例3. 1 实现方法 实现上节说明的有限元方法求解偏微分方程主要有两 种方法 , 一种是从偏微分方程出发 , 将有限元方法的各个求 解步骤分别编写计算机程序 , 然后代入具体问题求解并得出 结果 , 这个方法编程工作量大但方法灵活能解决一些特殊问 题 , 对试验一些算法是也是必要的 。 对于大多数偏微分方程 应用问题 , 可以选用已有的偏微分方程工具 , matlab 的偏微 分方程工具箱提供了研究和求解空间二维偏微分方程的一 个强大而又灵活实用的环

16、境 , 它的功能包括 设置 P DE 定解问题 , 即设置求解区域 、 边界条件以及 方程的形式和系数 ; 用有限元法求解 P DE, 包括网格划分 、 方程离散化以 及求出数值解 ; 结果的可视化及数据输出 说明工具箱的用法 。 3. 2 图 1 电机结构图1 电机结构图假设电机较长 使得端面影响可以 忽略 不 计 , 从 而 可 以简化为二维模型 处理 。 如图 1所示 。2 计算参数和流程静磁场的双旋 度控制方程可以转 化为标量椭圆型偏 微分方程 。在 mat 2lab 中输入 pdet ool 进入工具箱 。首先对电机建模 , 这里对整个电机建模以获得更好的视觉效果 , 在 D ra

17、w 模式可以用图 形交互模式或命令方式建模 , 下面语句说明了建模过程 , 为 节省篇幅 , 这里没有对所有线圈建模 。pdecirc (0, 0, 0. 08, C1 pdecirc (0, 0, 0. 086, C2 pdecirc (0, 0, 0. 16, C3 pdecirc (0, 0, 0. 23, C4pderect (-0. 040. 040. 050. 2, R11 pderect (0. 040. 0680. 080. 135, R12 pderect (0. 080. 1350. 040. 068, R13 pderect (0. 050. 2-0. 040. 04,

18、R21之后进入 Boundary 模式删除多余线条 , 并加入边界条 件 , 这里忽略电机外漏磁场的因素 , 故将模型的整个外圆边 界设置 D irichlet 边界条件 , 令 h =1r =0, 即 A =0。然后进入 P DE 模式 , 显示子区域编号如图 2a 所示 。设置 P DE 参数 , 在这个应用中需要设置各个子区域的磁导率 和电流密度 J 。 在线圈区域 4, 6, 10, 11中设置为 =1, J =2, 线圈区域 3, 5, 7, 8中 =1, J =-2, 定子和转子区域 1, 2中 J=0, 导磁率设置为 5000. /(1+0. 053(ux . 2+uy . 2

19、+200以模拟材料的非线性磁特性 , 其中 ux . 2+uy . 2等于 | A |2, 在空气区域 9设为 =1, J =0。2Mesh 模式 , 在这里可以设置网格的划分方 式如最大边长度 , 网格细化比例等 。对电机模型采用初始划 分方式 , 将整个区域划分为 6474个三角形单元如图 2b 所 示 。 下一步进入 Solve 模式 , 由于材料的非线性特性选择非 线性求解器 , 非线性误差设为 1e -4。 然后求解 。3 计算结果及显示求解后对结果如图 3显示 , 这里用箭头显示磁通密度 , 灰度轮廓线显示静磁矢势的等势线 , 试验结果与预期结果相 同。图 3 计算结果3. 3 金

20、属板的热传导问题1 问题描述一个 带 圆 孔 的 金 属 板 热 传 导 问 题 。板 的 左 边 保 持100 , 右边热量从板向环境空气定常流动 , 其他边及内孔边界保持绝缘 。 初始 t =t0时板的温度为 0 , 可概括为抛物型 方程定解问题 :d 5u /5t -u =0(15边界条件 :u =100, 左边界上 ; 5u /5n =-1, 右边界上 ; 5u /5n =0, 其它边界上 ; u |t =t 0=02 计算参数和流程进入工具箱后首先对金属板建模 ,pderect (-0. 50. 5-0. 80. 8, R1 ; pdeelli p (0, 0, 0. 2, 0. 2

21、, 0, E1 ; set (findobj (get (pdefig, C hildren, (下转第 115页 201 多 ; 方位向上 , 对于主瓣宽度各方法相差不大 , 区别在最大旁 瓣级上 , 最邻域近似法的旁瓣级虽然较未矫正时低了十几分 贝 , 但和其他两种方法相比 , 旁瓣级还远不够低 。可见拉格 朗日法和香农法在改善方位向旁瓣级上有优势 , 香农法更好 一点 , 但耗时太长 , 实际应用中可与最邻域近似法结合应用 。5 结论本文分析了 S AS 产生距离徙动的原因 , 详细介绍了用于 徙动矫正的几种插值方法 , 并针对距离 -多普勒算法对点目 标进行仿真成像 , 比较了这几种插

22、值方法的效果及区别 。三 种方法在改善距离向性能时相近 ; 作用于方位向时 , 最邻域 近似法用时最短 , 但方位向旁瓣级较高 , 而拉格朗日插值和 香农插值均能得到较低的方位向旁瓣级 , 性能优于最邻域近 似法 。 但香农法用时较长 , 实际中需酌情采用 。 参考文献 :1 蒋小奎 . 北工业大学博士论文 , 2004.2 封建湖 , 车刚明 , 聂玉峰 . 数值分析原理 M.北京 :科学出版社 , 2002.3 刘纪元 , 周荫清 . 星载 S AR 距离徙动分析与矫正 J .沈阳航空工业学院学报 , 1994(27 :54-62.4 R ichard Ba m ler . A compa

23、ris on of range -dopp ler and wavenumberdomain S AR f ocusing algorithm sJ .I EEE Trans . Geosci . Re mote Sensing, 1992, 30(4 :706-713.5 Mehrdad Soumekh . Synthetic Aperture Radar signal p r ocessing withMAT LAB algorithm sM.John W iley &Sons Publicati on, 1999.作者简介 丁迎迎 (- , (汉族 , 河南偃师人 , 硕士, 水

24、下目标成像 等 ;(1965. 10- , 女 (汉族 , 河南许昌人 , 博士 , 教授 , 博士生导师 , 研究方向为阵列信号处理、 自。(上接第 102页 Tag , P DEEval , String , R 1-E1之后进入 Boundary 模式加入边界条件 , 根据问题将左边 界设为 D irichlet 边界条件 , 令 h =1, r =100; 右边界设为Neu mann 边界条件 , 令 g =-1, q =0; 其它边界设为 Neu mann边界条件 , 令 g =0, q =0。在 P DE 模式中设置为 Parabolic 型方程 , 系数设置为 :d=1, c =1, a =1, f =0。在 MES H 模式中在初始划分基础上进一步细分 , 将求解 区域分为 768个三角形单元 。在 S OLVE 模 式 中 将 求 解 参 数 设 置 为 :Ti m e:0:10,u (t 0 =0, Relativ