有限元法及其在数学物理中的应用.pdf

乐山师专学报 1986 年第 2 期 己 门 司 刁一一 口 州 有限元法及其在数学物理中的应用 何昆弟 摘 要 本文介绍了有限 元方法 的数学理论 以及应用有限 元方 法的一般规律 同 时 讨论了 在 数学上 用有限元方法求解偏微分方程 及在物理学中求解的 问题 前言 有限 元方法是一种离散化的数值计算方法 它是在4 f 年代就已提出来 直到6 0 年代随着 电子计算机的发展 在解决工程力学问题中发 展起来的 目前有限 元方法 已成为求解偏微分 方程数值解的一个重要方法 它的理论基础 牢靠 物理概念清楚 解题效能高 适用性强 不但可以用来解决固体力学中的大部分问题 而且已经 渗透 到流体力学 热传导 电磁场等 学科领域 由于它对于物理 几何条件复杂的 问 题的特 别有效性 冶 使它 广 泛应 用 于工程技术领域 如 航空 航天 造船 建筑 机械等工程部门 J 有限 元方法 与最 优化 方法相结合 为电子计算机辅助设计 C AD 开 辟了广 阔的前景 有限 元方法的 基础 一是变分原理 二是剖分插值 也就是说 它是古典变 分法 Rit z 一G a l er kin方法 与分块多项式 插值 的结合产物 这 种结 合不仅使有限 元方 法保持了原有变 分法的优点 而且还兼有差分方法 的灵活性 是古典变分法的革新和发展 使古典变分法大 大向前推进了一步 有限 元方 法 的发展 借助于两个重要的工具 一是 在理论推导 中采用了矩阵方法 二是 在 实际计算中采用了电子计算机 有限 元 矩阵 电子 计算机是三位 一体的 本文只讨论有限 元方 法的理论基础 解题的 一般规 律以及在数学物理方 面 的应 用 至 于 在工程技术方面的应用 将另文专题讨论 有限元方法的理论 基础 一 变分原理 1 泛函及其极值 128 先考虑弹性地基梁的例子 设有一个放在弹性地 基上的梁 承受分 布横向载荷 q x 的作用 一端 x o 固结 另一端 x 二 l 自由 问梁取 怎样的 挠度w x 能使这个系 统的总势能二取 最 小值 设梁的弯曲刚度为E J 于是梁的弯 曲 应变能 二 是 一合 E 监羚 二 再设地基的刚度系数为k 则地基中贮存的 能量 二f为 1 T t f 二一龟 长W 一U X 乙J o 由于梁的挠度 载荷 q 的势能变为 二一 q wd 二 因此系统的总势能为三 者 之和 阁 一 以欲鲁 加上边界条件 在 x o 处 kw 一q w dl 1 一 2 W 一 xd一d 这就是一个变 分 问题 另外 常见 的最速降线问题 及悬链 线 问题也属这一类 万了d了乎 几11 卜 一 万 丫 dX I 一 dX 了Zg y 最速降线问题 l e s e s 一一T M y 之 斋 悬链线题问 上 述 问题 的一般提法是要找出函数 或曲线 y y x 或w w x 使得由一巳知 函 数F x y y 通过两 定点y y x 及y Z y x Z 所组成的积分 I y I F x y y d x 1 1 X 为极大或极小 这个积分I y 是一个泛函 其中 y 二 y x y y 1 称为强迫边 界 条件 I y 的极大或极小值称为极值或驻值 2 泛 函 的变分 设对某 函数y 二 y x 给以 变分 即增量6y 二 各y 幼 函数从y 卜 y 6y 则 相 应地 I y 脚甘今I y 乙y 不难证明 129 I y 各y I 了 I 晋 么 1 2 且巾 f f a F a F 升甲 各I 飞 二竺 各y 全专sy ld x J y a y 叫做泛函I 黔的一次变分 简称变分 利用分部积分公式可以使 1 3 式变为 1 3 x一 f f里旦 一其乙 里耳 d 二 JI a y dx 刁y l 4 2 五 X X y 6 F 一 列 分 一 乡 由强迫边界条件知 y在两端 x 及x Z 为已知 即乙y在两端不能变化 从而 1 4 式简化 为 卜 l 器 一 告 翻 yd 1 5 3 泛函的极值条件 泛函I y 的极值条件为 各 I 0 1 6 由 1 5 式知 有 口F d 乡F 一二二一 一二竺 二二 O a y dx ay I 称为E ule r 方程 由 1 1 式及5与l的可交换性知 色 F 0 4 多个独立变量 的情形 1 7 1 6 式等价于 1 8 若 F F x 了 y i 1 2 I 了f 由极值原理有 l F X i X 1 2 1 9 a F O r 夕 曰 如 气尸 尸 a yl 1二l 各yi 0 等价于 1 10 l n nU 一 一一 一 F 一y F 一y 乡一 日分 乡 因此 F的极值问题转化为 等价于 解方程 组 l 1 0 的 问题 以上讨论了单重积分的 泛函 但不难推广到二重 面 三重 空 间 积分的情形 二 几何剖分与分片插值 工30 这部分主要讨论平面区域的几何剖分和相应的分片插值方法 对区域O进行剖分时 基 本单元可取为三角形 四边形等 插值函数可取为一次 二次或高次多项 式等 其中以三角 形单元和相应的三顶点线性插值最简单 1 几何剖分 如 图2所示的平面区域9其边界 Q r r 是 曲的 但总可以用适当的折线 来逼 近 从而把Q剖分为一 系列 的三角形 即剖分 为 点元 A A Z A 线元 B B Z B 二 面 元 C C Z CI 其中 n 一节点总数 m 一线元总数 l一面元 总数 点 线 面元都分别 编号 并给出下列信 息 1 节点坐标 x y i i 2 n 2 线元两顶点编号 m i m i i 二 1 2 m 3 面 元三顶点编号 n n n i i 1 2 l 从此 区 域Q的几何剖分 便完全确定 2 三角形 单元上 的线性插值 在有限元的离散化中 特解函数 u x y 在各单元 上 用适 当的插值函数来代替 最简单的插值方法就是线性 插 值 它 也是其他插值方法的 基础 上 设有任意三角形C 二 A A A 3 顶点A 的坐标为 月 x 为 u U 使得 U U y u x 1 一 2 3 设某函数 u x 力在顶点的 值 2 3 要求作线性函数 x y ax by e 1 11 x y ax by l e u x y Z ax by e 二 u 不 1 1 2 尹 U x 3 从 1 12 y 3 ax 3 by e u3 式可以解出 a b c 代回 1 11 式 得到 3 U x y 乙 N u i 1 131 其中 N N i 盆 y 1 2 3 称为三角形单元上线性插值的形 或基 函数 其本身也是线性函数 且满足 在顶点处 1 当i 二 j N x j ys 6 s 火0当i子j 形函数N 是线性的 它们的偏导数是常数 这种插值方法很容易 推广 到线元 上 二次以上 的插 值 3 6 与此相似 其数值计算 的精度 u i 1 2 n 将会比线性插值有所提高 二 有限元方法的一般规律 下面将以椭圆型方程边值问题为例来说明用有限 元方法求解问 题的一般规律 取平面域 Q上的变系数椭圆型方程 Q 一 卫 了 日里竺 里一日里生 二 f 刁X刁X口y乡y l 2 1 这里日 队 x y o f f x y 都是预先给定的 这类方程具有代表性 物理上众 多 的定态 问题都归结为这个典型方程 或其简化 或推广了的形式 如定常态热传导或扩散 静电磁场 不可 压缩无旋流 定常渗流 定常亚声速流等 为简便计 只讨论方程 2 1 的第一类边值问题 其在边界 a g r r 产上 的条件 是 a Q u y x y 2 2 对于方程 2 1 及边界条件 2 2 可以构成所谓 能量积 分 即 泛函 I 含 器 十 日 器 一 d dy 2 3 Q 可以证明 7 一 所有满 足边界条件 2 3 的 u 函数族中使得I u 达 到极小的 那个函数 即 泛函的极值 u二 u x y 必定在Q内满足方程 2 1 反之 满足 2 1 和 2 2 的函数 u u x y 也必定使 泛函I u 达到极小 值 也 就是 说变 分问题 I 卜 杏 一 器 日 器 一 一 d dy 二m 2 u 二甲 二 2 5 等价于 边值 问题 了 一 会 日登 命 日 登卜 飞 u二 二 2 6 2 7 即两者有相同的解 由此可见 求解偏微分方程 的问题就转化成了变 分 问题 因此 1 首先要找出问题的变分表示形式 如 2 3 2 3 式可以表式为下述 的简化形式 I u u d x d m in Q 2 8 132 2 将平面域 剖分为若干个单元 如图2所示 单元 的大小虽然可以有很大的任意性 但 一个好 的 剖分必须兼顾计算机容量 运算速度和计算精度的要求 因此 剖分必须遵循下述 几点 1 每个单元 的顶点只能是相邻单元 的顶点即不能是相邻单元的内点 如图4 今 长 卜 丫 嘴 夕 场奄 夺 执 闪 对 2 尽量避免出现大的钝角和大 的边长 如图5 否 则将会严重 影 响计算精度 3 在u x y 的梯度变化可能剧烈 的地方 剖分 的网格要密 变化小的 地方 网格可 以稀一点 4 单元编号可以任意 但长点编号 应 使相邻两节点编 号之差 的绝对值中最大者愈小愈 好 因其直接决定线性代数方程 组系数矩 阵 亦叫刚度矩阵 的带宽 3 单元 分析 对方程 2 4一5 的定解域9 径几何剖分以后 自然 可以表示为所有单元之和 于是 能量积分I u 可以分解为有关单元 上能量积分之和 I u E x u e Q e 乳l 合 日 器 日 器 一 d 一 2 9 如图 2 在单元 剖分的 基础上 将未知函数u x y 代以由它在节点A A 的值u u u 产生 的分片插值函数 3 u x y 兄 N iu i l 2 10 a U 3 厌牙 乙 i 1 aN a U 3 E u i 二 1 旦旦 aXI 3 兄 i j 1 a N Ns f旦U 刁Jl 3 云 j 1 aN a N 里过止 2 二 A 13 3 所以 1 u I U 什 尝 令 一 fu d d 或 1 1 u u 3 1 一 U 吸 少 2 曰 a 不e 1 J 3 u u 一 习 e u 2 1 2 i j 1 i 1 其中 U 嚼 努 豁 翻 fN dxdy 2 1 3 2 14 户J 产 e e 矛 a D e 容易证明 7 刀 矛 二 嗜 i j 1 2 3 2 14 从而说明 单元刚度矩阵 州 卿 A e 介钾 仁 冷 冷 是对称阵的 而 b b 圣 e b釜 e b l e T 叫单元载荷向量 4 总体合成 在单元分析中得到的 A e 及 b 按节点编号 采取对号入座的方式 累加到总刚 度矩阵 A 及总体载荷向量 b 的相应元素位置上 由于 A 是对称的 因此总刚度矩阵 A 也是对称 的 其正定性是易于证明的 7 这样能量积分即泛函I u 就完全离散化为多元二次函数 I u l u u 万乙 n j 一 乙 b u i 1 2 15 i j 二 1 5 求解线性代数方程组 根据微积分中的极值原理 1 1 5 当二阶导数矩阵 名I a ui刁以i 134 正定时 极小问题 2 4 等价于解方程组 a I 刁Ui 由 2 1 5 可知 I是二次的 它 的一阶导数是一次的 即 2 16 a I 乡 Ui 习 a u 一 b 二o j l i 1 2 n 2 1 7 从 而有 n 乙 a u b 0 1 2 一 n 2 1 8 j二 1 或 A u b卜 2 1 9 其中 A 为n xn 阶的正定阵 谧 u u u u T b b b Z b T 在求解线性代数方程 2 1 9 之前 还要对约束进 行处理 求解 2 1 9 的方法有多种 二 并且已程序化并应用于计算机 通 过求解 2 19 就能得到各节点的 u i 1 2 n 值 6 结果分析 在求解 u 一次单元 分析 uZ 一 后 在实践中常 常需要知道日 哭 一 日 畏牛 的分布 因此尚需再作 产J 产甘 即按原来的插值原则 补插出所需要的量 三 有限元方法在数学物理中的应用 一 常微分方程边值 问题的有限元方法 在区间 幻内考虑其第一边值问题 厂 d 一 下 p 刃 十qy x 七 又 y o y l 0 3 1 3 2 其 中 p 二 p 二 q q x 率 QX 方程 3 1 一 2 的 变分 问题是 在函数族D 中求函数y x 使泛函 Jr 显然 y p y 尹名 qy么一2 y X m 二 3 3 y o y l 0 这是一个一维 的稳态I q题 对区域 0 l 用n z 个点二 x x x 进 135 行剖分 即剖分成 n 个线元 在每个线元 x 一 x i 1 计 只讨论线性插值函数 Y N T y 2 n 上构造y x 的插值函 数 此处 为简明 3 5 产八O nO 声叮 l e s rl J Z n 其中 N 一 N手 y N 一 N I N T y卜 l y i T X一X 泣 X 卜 z一xi N i 1 2 X一X获一1 X i一泣一Xi N 称为i单元上的形函数 又因为 dy d N T dxdx y B T y 其中 代入 B 擎 Q X 塑 卫 日 I d N dx 3 7 3 3 式 得 p 凌B T y 么 q N T y 一 Zf N T y d x XX 了 气 1 一一 J y e y 一2 1 xi f p B B T q N N T d x y少 1一 N T d x y y T k y 一2 F 丁 y 3 8 X s一1 其中 k k k 单元刚度矩阵 l e e 了 l j J k p 褚B B Tdx X 卜 t X 丈 k 二 q N N rr d x 3 9 x 卜 l 又亩 l J 了 f IF T二 I f凌N Tdx 在单元分析的 基础上 的分别求得单元刚度矩阵 K e 及单元载荷向量 F 再按第二部分 所述的 步骤进行总体合成的约束处理 求解线性代数方程组 从而求得y x 在节点 x x x 一 处的近似值 y t y Z y卜 二 椭圆型方程边值问题的有限元方法 这个问题 已在本文第二部分中作为例子 进行了详细的讨论 椭圆型方程 的边值问题 几乎涉及了物理学中的所有静态或稳态 问题 其应用是很广泛的 136 三 抛物型和双曲型方程边值 问题 的有限元方 法 物理学中热传导 渗流等一类问题就属于抛物型方程的边值问题 现以渗流方程及其边 值问题为例来讨论 里二 p三竺 里一 一 里牛 Q 一 e 理些 ax ax ay 一 刁y 1 at 3 10 u x y o u lr 0 竺 r 刁n J p x y x Q 定压 不渗透 3 1 1 口 l l e s l 其中 a o r r 产 是O的全部边界 只要在某一瞬时 将 aU 乡t 看作仅是 x y 的 函数 则瞬态 问题就可化成稳态 问题 来处理 从而 3 1 0一11 的泛函极值问题为 1 卜 瓦 器 刹 一 Q 一 C 器 d 一 3 一 完全仿照本文第二部分的讨论 在每个三角形单元上有 3 U x y t 乙 N u N T U 手 i 1 3 1 3 其中 N U N x y N N Z N T u e u e u l T 各量为 的函数 N T a N 丁 乡X 3 14 l l t 产 J e s 夕 d U dt U e 令 H a u 一 川 一 a x 川 一 3 15 将 3 1 3 3 1 5 代入 3 1 2 整理后可得 I U e P H T H d xdy U 卜 U 了飞 1 一 2 13 7 U N Qd Xdy一 11 N ee d U 蛋 dt 3 1 6 令上式右端中各项的积分 l s e 矛 l e s s e l K H H d二dy E e N N dX y 3 17 F e e N Qd 于是 I U 艺1 U 一 菩 u K u 一 菩 u T 王N Tdx dy F 书 u E 月釜 上 C T k C e F e U 习 e I JJ T 1 石一 IUI 一 u 月 T 万 C e T E C d U dt 3 1 8 C 为选择矩阵 旦卫亚 二 些 dt dt U d u dt u i u T J 心吸t即 旦 坚 生 T dt l K 奢 e K C F 是 e F 艺 L乙J e L七 J L乙 J L七 J 3 19 3 1 2 可得 处令从此再 则 K U 一 F E d U dt O 3 20 这是一个一阶常微分方程组 完全可用 通常常微分方程数值解法求解 物理学中的波动方程 三兰竺 g t 汤 一a 名 V 名u f 3 21 l召8 就是双曲型方程的典型例 子 这类边值问题 的有限元方法 可参照前面的讨论 最后得到一 个二阶常 微分方程组 d 名 U dt 名 E 担 监 生 K U 一 F 汉令 红 嘴黔 一 兴黔 3 22 其中 d 名 裙U圣 dt么 3 22 式 的数值 解亦是容 易求得的 参试 文 献 0 C Zie n kiewie z T heFiniteElem en t M etho d T hird e ditio n M egr aw 一 H all B o o kC o mp an y U K 19 7 7 E Hinto n an dD R J O nen FiniteElem ent Pr og r am ming A ea d em ie Pr e s s L o n d on 1 9 7 7 龙驭球