第五章 偏微分方程的有限元法

HarbinInstituteofTechnologyYangkunkyang 1 43 第五章偏微分方程的有限元法 5 1泛函与变分原理5 2基于变分原理的有限元法5 3matlab有限元法工具箱 HarbinInstituteofTechnologyYangkunkyang 2 43 第五章偏微分方程的有限元法 有限元法 FEA FiniteElementAnalysis FEM 有限元法的基本思想是用较简单的问题代替复杂问题 然后再对简单问题进行求解的数值计算方法 有限元法将求解域看成是由许多被称为有限元的小的互连子域组成 对每一单元假定一个较简单的近似解 然后推导求解这个域总的满足条件 从而得到问题的解 这个解不是准确解 而是近似解 有限元不仅计算精度高 而且能适应各种复杂形状 因而成为行之有效的数值计算方法 有限元法于上世纪50年代首先在力学领域 飞机结构的静 动态特性分析中得到应用 随后很快广泛的应用于求解热传导 电磁场 流体力学等连续性问题 有限元法主要用于求解拉普拉斯方程和泊松方程所描述的各类物理场中 HarbinInstituteofTechnologyYangkunkyang 3 43 第五章偏微分方程的有限元法 有限元法 变分原理 基于变分原理的有限元法是逼近论 偏微分方程 变分与泛函分析的巧妙结合 基于变分原理的有限元法以变分原理为基础 把所要求解的微分方程定解问题 首先转化为相应的变分问题 即泛函求极值问题 它将求解域看成是由许多称为有限元的小的互连子域组成 然后利用剖分插值 对每一单元假定一个合适的 较简单的 近似解 把离散化的变分问题转化为普通多元函数的极值问题 然后推导求解这个域总的满足条件 边界条件 即最终归结为一组多元的代数方程组 求解代数方程组 就得到待求边值问题的数值解 HarbinInstituteofTechnologyYangkunkyang 4 43 第五章偏微分方程的有限元法 有限元法 加权余数法 自从1969年以来 某些学者在流体力学中应用加权余数法中的迦辽金法或最小二乘法等同样获得了有限元方程 因而有限元法可应用于以任何微分方程所描述的各类物理场中 而不再要求这类物理场和泛函的极值问题有所联系 加权余数法的核心思想是 近似解与解析解相比会存在误差R 但是可以通过一个准则使R尽量小 求解这个等式 就可以得到待定常数的值 也就得到了近似解 HarbinInstituteofTechnologyYangkunkyang 5 43 第五章偏微分方程的有限元法 有限元法特点有限元法的物理意义直观明确 理论完整可靠 因为变分原理描述了支配物理现象的物理学中的最小作用原理 如力学中的最小势能原理 优异的解题能力 有限元法对边界几何形状复杂以及媒质物理性质变异等复杂物理问题求解上 有突出优点 不受几何形状和媒质分布的复杂程度限制 不必单独处理第二 三类边界条件 离散点配置比较随意 通过控制有限单元剖分密度和单元插值函数的选取 可以充分保证所需的数值计算精度 HarbinInstituteofTechnologyYangkunkyang 6 43 5 1泛函与变分原理 数学上 通常自变量与因变量间的关系称为函数 而泛函则是函数集合的函数 也就是函数的函数 即自变量为函数 而不是变量 5 1 1泛函的定义泛函通常是指一种定义域为函数 而值域为实数的 函数 设C是函数的集合 B是实数集合 如果对C中的任一元素y x 在B中都有一个元素J与之对应 则称J为y x 的泛函 记为J y x HarbinInstituteofTechnologyYangkunkyang 7 43 5 1泛函与变分原理 例5 1 1质点在重力作用下 沿一条光滑的从A点到B点的曲线运动 如图所示 求下落时间最短的曲线 曲线上任一小段线元长度为 捷线问题 HarbinInstituteofTechnologyYangkunkyang 8 43 5 1泛函与变分原理 线元处的质点速度为 ds线元下落时间为 从A点到B点的下落时间为 HarbinInstituteofTechnologyYangkunkyang 9 43 5 1泛函与变分原理 5 1 2函数的变分 设y x 是泛函J定义域内任一函数 如果y x 变化为新函数Y x 且Y x 属于泛函J的定义域 则Y x 与y x 之差为函数y x 的变分 变分 y是x的函数 它不同于函数的增量 y 性质 函数求导与求变分可以交换次序 HarbinInstituteofTechnologyYangkunkyang 10 43 5 1泛函与变分原理 5 1 3泛函的变分 定义最简泛函 F x y y 称为泛函的 核函数 泛函的变分 最简泛函 核函数只包含自变量x 未知函数y x 以及导数y x HarbinInstituteofTechnologyYangkunkyang 11 43 5 1泛函与变分原理 利用二元函数的泰勒展开 HarbinInstituteofTechnologyYangkunkyang 12 43 5 1泛函与变分原理 其中 分别称为泛函的一阶变分和二阶变分 HarbinInstituteofTechnologyYangkunkyang 13 43 5 1泛函与变分原理 泛函取极值的必要条件 一阶变分为零 性质 对于最简泛函 变分运算可以与积分 微分运算交换次序 HarbinInstituteofTechnologyYangkunkyang 14 43 5 1泛函与变分原理 5 1 4泛函的极值问题 泛函的一阶变分 利用 1泛函的极值问题的间接解法转化为微分方程 欧拉方程 HarbinInstituteofTechnologyYangkunkyang 15 43 5 1泛函与变分原理 对于驻定问题 两边界固定 这就是最简泛函的欧拉方程 等价于泛函取极值的必要条件 把变分问题转化微分方程的定解问题 边值问题 来求解 HarbinInstituteofTechnologyYangkunkyang 16 43 5 1泛函与变分原理 对于例5 1 1求下落时间最短的轨迹 利用最简泛函的欧拉方程 HarbinInstituteofTechnologyYangkunkyang 17 43 5 1泛函与变分原理 代入欧拉方程 HarbinInstituteofTechnologyYangkunkyang 18 43 5 1泛函与变分原理 变换得到 进一步化简得到 积分 HarbinInstituteofTechnologyYangkunkyang 19 43 5 1泛函与变分原理 做变量替换 得 而 HarbinInstituteofTechnologyYangkunkyang 20 43 5 1泛函与变分原理 对上式积分得到 这样就得到了下落时间最短曲线的参数方程 式中常数c1和c2由始末两点位置确定 练习 画出经过 0 0 和 1 1 的下落时间最短曲线 连接两个点上凹的唯一一段旋轮线 HarbinInstituteofTechnologyYangkunkyang 21 43 5 1泛函与变分原理 2泛函的极值问题的直接解法 基本做法 瑞利 里兹 Rayleigh Ritz 法 1 选定一组具有相对完备性的基函数 构造一个线性组合的近似函数 2 将含有n个待定系数的构造函数作为近似的极值函数 代入泛函 HarbinInstituteofTechnologyYangkunkyang 22 43 5 1泛函与变分原理 3 为了求泛函的极值 按照多元函数取极值的必要条件 4 求解以上方程组 求出就可以得到极值函数的近似解 5 再将含有n 1个待定系数的函数作为近似极值函数 重复 2 4 就可以得到极值函数新的近似解 如果连续两次所得到的结果接近 就认为最后得到的函数就是极值函数的近似解 HarbinInstituteofTechnologyYangkunkyang 23 43 5 1泛函与变分原理 例5 1 2求下列泛函的极值函数 解 为了满足边界条件 取基函数为 近似函数为 HarbinInstituteofTechnologyYangkunkyang 24 43 5 1泛函与变分原理 当n 1时 代入泛函 取极值 HarbinInstituteofTechnologyYangkunkyang 25 43 5 1泛函与变分原理 计算得到 近似函数 同理n 2时 利用欧拉方程 得到的精确解 HarbinInstituteofTechnologyYangkunkyang 26 43 5 1泛函与变分原理 HarbinInstituteofTechnologyYangkunkyang 27 43 5 1泛函与变分原理 泛函的极值问题可以通过变分运算产生一个微分方程和相应的边界条件 即欧拉方程 其解对应于最简泛函的极值函数 也就是泛函的极值问题可以等价为在一定边界条件下求解微分方程问题 变分原理通过求解一个相应的泛函的极小函数而得到偏微分方程边值问题的解 有限元法正是里兹法与有限差分法相结合的成果 它取长补短地在理论上以变分为基础 在具体方法构造上又利用了有限差分法网格离散化处理的思想 HarbinInstituteofTechnologyYangkunkyang 28 43 5 1泛函与变分原理 20世纪60年代初首次提出结构力学计算有限元概念的克拉夫 Clough 教授形象地将其描绘为 有限元法 RayleighRitz法 分片函数 有限元法是RayleighRitz法的一种局部化情况 不同于求解满足整个定义域边界条件的允许函数的RayleighRitz法 往往是困难的 有限元法将函数定义在简单几何形状 如二维问题中的三角形或任意四边形 的单元域上 分片函数 且不考虑整个定义域的复杂边界条件 这是有限元法优于其它近似方法的原因之一 HarbinInstituteofTechnologyYangkunkyang 29 43 5 2基于变分原理的有限元法 对于具有不同物理性质和数学模型的问题 有限元法的基本做法是相同的 只是具体公式推导和运算求解不同 有限元法基本做法 首先把待求的偏微分方程边值问题转化为等价的变分问题 然后通过有限单元剖分的离散处理 构造一个分片解析的有限元子空间 通过构造近似函数 把变分问题近似地转化为有限元子空间中的多元函数极值问题 由此直接利用RayleighRitz法探求变分问题的近似解 极值函数解 以此作为所求边值问题的近似解 HarbinInstituteofTechnologyYangkunkyang 30 43 5 2基于变分原理的有限元法 有限元法具体求解步骤 建立积分方程根据变分原理或方程余量与权函数正交化原理 建立与微分方程初边值问题等价的积分表达式 这是有限元法的出发点 区域单元剖分根据求解区域的形状及实际问题的物理特点 将区域剖分为若干相互连接 不重叠的单元 区域单元划分是采用有限元方法的前期准备工作 这部分工作量比较大 除了给计算单元和节点进行编号和确定相互之间的关系之外 还要表示节点的位置坐标 同时还需要列出自然边界和本质边界的节点序号和相应的边界值 HarbinInstituteofTechnologyYangkunkyang 31 43 5 2基于变分原理的有限元法 确定单元基函数根据单元中节点数目及对近似解精度的要求 选择满足一定插值条件的插值函数作为单元基函数 有限元方法中的基函数是在单元中选取的 由于各单元具有规则的几何形状 在选取基函数时可遵循一定的法则 单元分析将各个单元中的求解函数用单元基函数的线性组合表达式进行逼近 再将近似函数代入积分方程 并对单元区域进行积分 可获得含有待定系数 即单元中各节点的参数值 的代数方程组 称为单元有限元方程 总体合成在得出单元有限元方程之后 将区域中所有单元有限元方程按一定法则进行累加 形成总体有限元方程 HarbinInstituteofTechnologyYangkunkyang 32 43 5 2基于变分原理的有限元法 边界条件的处理一般边界条件有三种形式 对于第二类边界条件 一般在积分表达式中可自动得到满足 对于第一类边界条件和第三类边界条件 需按一定法则对总体有限元方程进行修正满足 解有限元方程根据边界条件修正的总体有限元方程组 是含所有待定未知量的封闭方程组 采用适当的数值计算方法求解 可求得各节点的函数值 HarbinInstituteofTechnologyYangkunkyang 33 43 5 2基于变分原理的有限元法 有限元分析可分成三个阶段 前置处理 计算求解和后置处理 前置处理是建立有限元模型 完成单元网格划分 后置处理则是采集处理分析结果 使用户能简便提取信息 了解计算结果 HarbinInstituteofTechnologyYangkunkyang 34 43 5 2基于变分原理的有限元法 1 求解区域离散 离散单元基本要求 各单元只能在顶点处相交 不同单元在边界处相连 既不能相互分离又不能相互重叠 各单元节点编号循序应一致 一律按逆时针方向 从最小节点号开始 同一单元节点编号相差不能太悬殊 对多区域的编号 按区域连续编号 把求解区域分割成有限个单元体的集合 单元体形状原则上是任意的 一般取有规则形体 有限元法计算步骤 HarbinInstituteofTechnologyYangkunkyang 35 43 5 2基于变分原理的有限元法 三角单元是经常使用的单元剖分方法 剖分时应注意几下几点 三角形不能重叠 不能把一个三角形的顶点取为相邻三角形的边上 剖分的三角形应该避免钝角 三角形不可过于狭长 最长边一般不大于最短边的3倍 三角形三边之比尽量接近1 不能把一个三角形跨越不同的介质 每个三角形最多只有一个边在边界上 三角形单元面积越小 计算精度越高 HarbinInstituteofTechnologyYangkunkyang 36 43 5 2基于变分原理的有限元法 把求解区域划分m个三角形有限单元 共有n个节点 在有限单元e j k l 上进行分片线性插值 插值函数为 2 选择近似函数 HarbinInstituteofTechnologyYangkunkyang 37 43 5 2基于变分原理的有限元法 在单元节点上 求解以上方程组可以得到 3 求解单元形函数 HarbinInstituteofTechnologyYangkunkyang 38 43 5 2基于变分原理的有限元法 同理可以求出 HarbinInstituteofTechnologyYangkunkyang 39 43 5 2基于变分原理的有限元法 HarbinInstituteofTechnologyYangkunkyang 40 43 5 2基于变分原理的有限元法 则插值函数可以写为 单元形函数 基函数 HarbinInstituteofTechnologyYangkunkyang 41 43 5 2基于变分原理的有限元法 三角元e插值函数可以改写为矩阵形式 HarbinInstituteofTechnologyYangkunkyang 42 43 5 2基于变分原理的有限元法 下面以泊松方程为例讨论有限元解法 所对应的泛函为 4 建立单元特征式 难点 寻找与微分方程对应的泛函 HarbinInstituteofTechnologyYangkunkyang 43 43 5 2基于变分原理的有限元法 在第e个三角元的泛函 由于 HarbinInstituteofTechnologyYangkunkyang 44 43 5 2基于变分原理的有限元法 改写为矩阵形式 HarbinInstituteofTechnologyYangkunkyang 45 43 5 2基于变分原理的有限元法 其中 同理 HarbinInstituteofTechnologyYangkunkyang 46 43 5 2基于变分原理的有限元法 三角元e的泛函 其中 HarbinInstituteofTechnologyYangkunkyang 47 43 5 2基于变分原理的有限元法 改写Ke到所有n个节点 即把扩充部分添零 以方便总体矩阵的处理 其中 HarbinInstituteofTechnologyYangkunkyang 48 43 5 2基于变分原理的有限元法 求解区域上的总体泛函 其中 变分问题被离散化的多元二次函数的极值问题 5 建立系统有限元方程 HarbinInstituteofTechnologyYangkunkyang 49 43 5 2基于变分原理的有限元法 根据多元函数极值理论 得到第i点有限元方程 即 求解上述有限元方程 线性代数方程组 就可以得到节点上的函数值 HarbinInstituteofTechnologyYangkunkyang 50 43 5 2基于变分原理的有限元法 获得有限元方程之后 就可以选择各种方法求解相应的代数方程组 常用方法有高斯消去法 列元素消去法 迭代法等等 在变分问题中第二类 第三类边界条件已经自然包含在泛函达到极值的要求中 不必单独处理 称为自然满足的边界条件 只需考虑第一类强加边界条件 强加边界条件的处理方法因代数方程组的解法而异 6 有限元方程求解与边界条件处理 HarbinInstituteofTechnologyYangkunkyang 51 43 5 2基于变分原理的有限元法 迭代法求解 凡是遇到边界节点所对应的方程均不迭代 节点值始终保持给定值 不必单独处理边界 直接法求解 节点m为边界 函数值um u0 处理方法为 把对角元素的特征元素设置为1 即kmm 1 然后把m行与m列的其它元素全部设置为0 方程的等式右边改为给定的函数值u0 其它元素则要减去该节点处理前对应的m列的特征系数kim与u0的乘积 HarbinInstituteofTechnologyYangkunkyang 52 43 5 2基于变分原理的有限元法 例5 2 1一个边长为1的二维正方形静电场域 电位函数为 x y 边界条件如图所示 试用有限元法确定二维静电场域的电位分布 解 该二维静电场域的电位函数 x y 可以用下列第一类边界条件的偏微分方程描述 HarbinInstituteofTechnologyYangkunkyang 53 43 5 2基于变分原理的有限元法 按照右图进行三角形单元剖分 单元编号按照从左到右 从下到上的顺序编号 节点编号 1 0 0 2 0 1 3 0 5 0 5 4 1 0 5 1 1 三角形单元编号 e j k l 单元内顶点按逆时针编号1 1 3 2 2 1 4 3 3 2 3 5 4 3 4 5 第一类边界条件 HarbinInstituteofTechnologyYangkunkyang 54 43 5 2基于变分原理的有限元法 对于三角元1 1 3 2 HarbinInstituteofTechnologyYangkunkyang 55 43 5 2基于变分原理的有限元法 扩展到全部节点 HarbinInstituteofTechnologyYangkunkyang 56 43 5 2基于变分原理的有限元法 同样 对于三角元2 1 4 3 HarbinInstituteofTechnologyYangkunkyang 57 43 5 2基于变分原理的有限元法 扩展到全部节点 HarbinInstituteofTechnologyYangkunkyang 58 43 5 2基于变分原理的有限元法 同样 对于三角元3 2 3 5 HarbinInstituteofTechnologyYangkunkyang 59 43 5 2基于变分原理的有限元法 扩展到全部节点 HarbinInstituteofTechnologyYangkunkyang 60 43 5 2基于变分原理的有限元法 同样 对于三角元4 3 4 5 HarbinInstituteofTechnologyYangkunkyang 61 43 5 2基于变分原理的有限元法 扩展到全部节点 HarbinInstituteofTechnologyYangkunkyang 62 43 5 2基于变分原理的有限元法 全部节点K HarbinInstituteofTechnologyYangkunkyang 63 43 5 2基于变分原理的有限元法 有限元方程 可以采用迭代法和直接解法 求解此线性代数方程组 HarbinInstituteofTechnologyYangkunkyang 64 43 5 2基于变分原理的有限元法 迭代法 迭代公式为 代入初值 HarbinInstituteofTechnologyYangkunkyang 65 43 5 2基于变分原理的有限元法 直接解法 需要处理第一类边界条件 HarbinInstituteofTechnologyYangkunkyang 66 43 5 2基于变分原理的有限元法 差分法 得到差分递推公式 HarbinInstituteofTechnologyYangkunkyang 67 43 5 2基于变分原理的有限元法 差分网格 利用边界条件 HarbinInstituteofTechnologyYangkunkyang 68 43 5 3matlab有限元法工具箱 大型通用有限元商业软件 国外软件ANSYS ADINA ABAQUS MSC等国内软件FEPG JFEX KMAS等 Matlab偏微分方程工具箱 PDEToolbox 提供了利用有限元法 图形界面求解偏微分方程的计算环境 PDEtool有较大的局限性 可以求解特殊PDE问题 比如只能求解二阶PDE问题 并且不能解决偏微分方程组 但是它提供了GUI界面 从繁杂的编程中解脱出来了 同时还可以通过File SaveAs直接生成M代码 HarbinInstituteofTechnologyYangkunkyang 69 43 5 3matlab有限元法工具箱 PDEToolbox求解偏微分方程类型 1椭圆型方程 Elliptic 2抛物线型方程 Parabolic 3双曲型方程 Hyperbolic HarbinInstituteofTechnologyYangkunkyang 70 43 5 3matlab有限元法工具箱 4特征值方程 Eigenmodes 上述微分方程中 c a d f在椭圆型方程中可以为函数 但在其它方程中必须为常数 HarbinInstituteofTechnologyYangkunkyang 71 43 5 3matlab有限元法工具箱 PDEToolbox边界条件 1狄里赫利条件 Didchlet 2诺依曼条件 Neumann n为边界上的单位外法线矢量 h r q g可以为函数 HarbinInstituteofTechnologyYangkunkyang 72 43 5 3matlab有限元法工具箱 PDEToolbox启动 1启动 2界面 HarbinInstituteofTechnologyYangkunkyang 73 43 5 3matlab有限元法工具箱 PDEToolbox菜单 Options 打开或关闭栅格调整栅格大小打开或关闭捕捉栅格功能绘图轴的坐标范围打开或关闭绘图方轴关闭帮助信息图形缩放选择应用模式重新显示图形 HarbinInstituteofTechnologyYangkunkyang 74 43 5 3matlab有限元法工具箱 PDEToolbox菜单 Draw 进入绘图模式对角点绘矩形固定中心绘矩形矩形对角点绘椭圆固定中心绘椭圆绘多边形旋转已选图形将几何描述矩阵输出到主工作空间 HarbinInstituteofTechnologyYangkunkyang 75 43 5 3matlab有限元法工具箱 PDEToolbox菜单 Boundary 进入边界模式对已选边界输入条件显示边界区域标识开关显示子区域标识开关删除已选的子域边界删除所有的子域边界将分解几何矩阵 边界条件矩阵输出到主工作空间 HarbinInstituteofTechnologyYangkunkyang 76 43 5 3matlab有限元法工具箱 PDEToolbox菜单 PDE 进入偏微分方程模式显示子区域标识开关调整PDE参数和类型将PDE参数输出到主工作空间 HarbinInstituteofTechnologyYangkunkyang 77 43 5 3matlab有限元法工具箱 PDEToolbox菜单 Mesh 输入网格模式初始化三角形网格加密当前三角形网格优化网格退回上一步用数字化的颜色显示网格质量 大于0 6可接受显示网格节点标识显示三角形网格标识修改网格生成参数输出网格矩阵到主工作空间 HarbinInstituteofTechnologyYangkunkyang 78 43 5 3matlab有限元法工具箱 PDEToolbox菜单 Solve 对已经定义的偏微分方程求解调整解PDE的参数输出解到主工作空间 Plot 显示图形解绘图参数设置输出动画 HarbinInstituteofTechnologyYangkunkyang 79 43 5 3matlab有限元法工具箱 PDEToolbox求解步骤 求解区域设置应用模式设置输入边界条件微分方程参数设定网格剖分初值和误差设置解方程图形解显示参数设置File SaveAs直接生成M代码 HarbinInstituteofTechnologyYangkunkyang 80 43 5 3matlab有限元法工具箱 例5 3 1如图带有矩形孔 0 1 0 8 的金属板 1 1 6 金属板左侧保持在100 右侧热量可以向环境定常流动 上下侧及内孔保持绝热 初始温度为0 求t 0 1 0 3 0 5 1 5s时金属板温度分布 解 此问题可以表示为如下定解问题 HarbinInstituteofTechnologyYangkunkyang 81 43 5 3matlab有限元法工具箱 求解区域设置提示符输入pdetool 1选择画矩形 2画矩形 3双击矩形 弹出对话框 输入准确矩形参数 HarbinInstituteofTechnologyYangkunkyang 82 43 5 3matlab有限元法工具箱 同样画出矩形孔 利用两个矩形运算得到求解区域 图形运算 HarbinInstituteofTechnologyYangkunkyang 83 43 5 3matlab有限元法工具箱 2 应用模式设置 HarbinInstituteofTechnologyYangkunkyang 84 43 5 3matlab有限元法工具箱 3输入边界条件 1点击 显示边界 2双击边界 弹出边界条件窗口 h 1 r 100 HarbinInstituteofTechnologyYangkunkyang 85 43 5 3matlab有限元法工具箱 输入每个边的边界条件 红色 Dirichlet蓝色 Neumann HarbinInstituteofTechnologyYangkunkyang 86 43 5 3matlab有限元法工具箱 4 微分方程参数设定 1点击 设置方程 弹出窗口 2抛物线型 d 1c 1a 0f 0 HarbinInstituteofTechnologyYangkunkyang 87 43 5 3matlab有限元法工具箱 5 网格剖分 点击 网格剖分 点击 加密网格 HarbinInstituteofTechnologyYangkunkyang 88 43 5 3matlab有限元法工具箱 6 初值和误差设置 单击Solve菜单中的Paramenters HarbinInstituteofTechnologyYangkunkyang 89 43 5 3matlab有限元法工具箱 7 解方程 点击 解方程 HarbinInstituteofTechnologyYangkunkyang 90 43 5 3matlab有限元法工具箱 8 图形解显示参数设置 HarbinInstituteofTechnologyYangkunkyang 91 43 5 3matlab有限元法工具箱 不同时刻温度分布 HarbinInstituteofTechnologyYangkunkyang 92 43 5 3matlab有限元法工具箱 不同时刻温度分布 HarbinInstituteofTechnologyYangkunkyang 93 43 5 3matlab有限元法工具箱 MATLAB除了提供有限元工具箱求解二阶偏微分方程之外 还提供了pdepe函数 可以直接求解一维椭圆型和抛物