有限元方法ppt课件.ppt

1 有限元方法 2 有限元法是求解偏微分方程问题的一种重要数值方法 它的基础分两个方面 一是变分原理 二是剖分插值 从第一方面看 有限元法是Ritz Galerkin方法的一种变形 它提供了一种选取 局部基函数 的新技巧 从而克服了Ritz Galerkin方法选取基函数的固有困难 从第二方面看 它是差分方法的一种变形 差分法是点近似 它只考虑在有限个离散点上函数值 而不考虑在点的邻域函数值如何变化 有限元方法考虑的是分段 块 的近似 因此有限元方法是这两类方法相结合 取长补短而进一步发展了的结果 在几何和物理条件比较复杂的问题中 有限元方法比差分方法有更广泛的适应性 3 7 两点边值问题的有限元方法 本节以两点边值问题为例 并从Ritz法和Galerkin法两种观点出发来叙述有限元法的基本思想及解题过程 7 1基于Ritz法的有限元方程考虑两点边值问题其中 4 1 写出Ritz形式的变分问题 与边值问题 7 1 7 2 等价的变分问题是 求 使其中 7 3 式 7 3 是应用有限元法求解边值问题 7 1 7 2 的出发点 5 2 区域剖分 剖分原则与差分法相同 即将求解区域剖分成若干个互相连接 且不重叠的子区域 这些子区域称为单元 单元的几何形状可以人为选取 一般是规则的 但形状与大小可以不同 对于一维情形最为简单 将求解区间分成若干个子区间 其节点为每个单元的长度为 单元在区间中分布的疏密程度或单元尺寸的大小 可根据问题的物理性质来决定 一般来说 在物理量变化剧烈的地方 单元尺寸要相对小一些 排列要密一些 6 设为的有限维子空间 它的元素为 要构造 只需构造单元基函数 构造单元基函数所遵循的原则是 其中 是单元节点序号为的节点 1 每个单元中的基函数的个数和单元中的节点数相同 每个节点对应一个基函数 本例中 单元有两个节点 因此基函数有两个 2 基函数应具有性质 7 3 确定单元基函数 有限元法与Ritz Galerkin方法的主要区别之一 就在于有限元方法中的基函数是在单元中选取的 由于各个单元具有规则的几何形状 而且可以不必考虑边界条件的影响 因此在单元中选取基函数可遵循一定的法则 8 9 10 11 7 6 令 4 形成有限元方程 便得到确定的线性代数方程组 称式 7 5 为有限元方程 12 7 8 7 7 值得注意的是 在实际计算中 并不是按照上述步骤形成有限元方程的 而是先进行单元分析 即在单元上建立有限元特征式 然后再进行总体合成 即将各单元的有限元特征式进行累加 合成为有限元方程 具体过程如下 第一步 单元分析 注意到 作变换 13 7 10 并引入记号 其中 于是 或写成 7 9 其中 从而有 7 11 这里 7 12 称为单元刚度矩阵 其中 7 13 7 16 7 14 式中 7 15 将式 7 11 7 14 代入式 7 7 便有 对式 7 7 右端第二项积分 有 这样 我们就得到了单元有限元特征式的一般表示形式 于是有 第二步 总体合成 总体合成就是将单元上的有限元特征式进行累加 合成为总体有限元方程 这一过程实际上是将单元有限元特征式中的系数矩阵 称为单元刚度矩阵 逐个累加 合成为总体系数矩阵 称为总刚度矩阵 同时将右端单元荷载向量逐个累加 合成为总荷载向量 从而得到关于的线性代数方程组 为此 记 从而式 7 16 右端第一个和式为 7 17 其中 未标明的元素均为0 这就是总刚度矩阵 对式 7 16 右端第二个和式 有 其中 这就是总荷载向量 从总刚度矩阵和总荷载向量的形成过程可以看出 的计算 实际上是把中四个元素在适当的位置上 对号入座 地叠加 的计算也是如此 我们引入 只是为了叙述方便 实际上 在编制程序时并不需要 显然 方程组 7 18 的系数矩阵是对称正定的三对角矩阵 因此可采用追赶法求出在节点上的近似值 7 18 其这样 就可将式 7 16 写成 因此 有限元方程为 19 7 两点边值问题的有限元方法 本节以两点边值问题为例 并从Ritz法和Galerkin法两种观点出发来叙述有限元法的基本思想及解题过程 7 1基于Ritz法的有限元方程7 2基于Galerkin法的有限元方程从Galerkin法出发形成有限元方程的过程与前面完全一样 针对边值问题 7 1 7 2 所得到的结果也是一致的 但是从Galerkin法出发形成的有限元方程更具一般性 它不仅适用于对称正定的算子方程 而且也适用于非对称正定的算子方程 所以我们今后主要是依据这一观点建立有限元方程 20 与边值问题 7 1 7 2 等价的Galerkin变分问题是 求 使得 7 19 其中仍用分段线性函数构成的试探函数空间替代 将代入 7 19 则得到所满足的线性代数方程组 7 20 这和方程组 7 6 是完全一样的 与容易看出 方程组 7 20 的系数矩阵就是总刚度矩阵 在总刚度矩阵形成的过程中 注意到 7 21 而从而有即故有这就是有限元方程 7 18 由上述看出 按Galerkin法推导有限元方程更加直接方便 尤其重要的是 按这一观点推导的有限元方程 不仅适用于定常的微分方程定解问题 而且也适用于不定常的微分方程定解问题 因此具有广泛的适应性 例7 1用有限元方法解边值问题将区间 0 1 等分成4个单元 解利用上述分析结果 我们只需构造出单元刚度矩阵和单元荷载向量 然后合成为总刚度矩阵和总荷载向量 注意到 7 13 和 7 15 并将形成单元上的中点值则不难得到其中 单元的中点为于是有 如果把单元刚度矩阵和单元荷载向量 扩大 便得到和为类似地 可写出和 然后进行叠加 便得到总刚度矩阵和总荷载向量 依边界条件即在中划去首末两行和首末两列 在中划去首末两行 便得到如下线性代数方程组 解之 得 8 二维椭圆边值问题的有限元方法 用有限元方法求解二维椭圆边值问题的过程与两点边值问题的有限元方法大体相同 只是在具体处理时比一维情形更复杂些 考虑如下椭圆型方程的第一边值问题 8 1 8 2 其中 是平面的一个有界域 其边值是分段光滑的简单闭曲线 以下我们从Galerkin法出发 叙述有限元求解问题 8 1 8 2 的全过程 与边值问题 8 1 8 2 等价的Galerkin变分问题是 求 使得 8 3 其中8 1区域剖分正如前章所言 对高维区域的剖分与对一维区域的剖分有很大不同 对一维区域无论作哪一种剖分 其单元仍然是一个区间 对不同的剖分只是区间长度不同而已 对高维区域而言 不同的剖分其单元的形状各异 如对二维区域 剖分后的子区域可以是三角形 矩形或四边形 限于篇幅 本书只讨论剖分后所得的子区域是三角形的情况 这种剖分称为三角形剖分 将区域划分成有限个三角形单元 剖分方法见前章 那里曾假定剖分的单元应是锐角三角形 现在我们去掉这一限制 只假定不同的单元是无重叠的内部 且单元的顶点不是其它单元边的内点 当然还要尽量避免出现大钝角的三角形 在物理量变化剧烈的地方 单元要划分得细密一些 变化缓和的地方 划分得稀一些 划分好单元之后 要对单元和节点进行编号 设是区域中的单元总数 将全区域中的单元统一编号 单元号记为 全区域中的节点也要按一定的顺序统一编号 记全区域中共有个节点 节点号记为节点编号的一般原则是尽可能使同一单元内的节点号比较接近 以后可以看到 单元内节点序号的差值决定了总体系数矩阵的带宽 8 2确定单元基函数与一维情形一样 为了构造试探函数空间我们只需在每个单元上构造插值基函数 这里 我们仅考虑三角形单元上的线性插值函数 为了便于后面积分的计算 我们先将直角坐标转换为面积坐标 1 面积坐标及有关公式 1 面积坐标的定义设是以为顶点的任意三角形单元 面积为 我们规定的次序按逆时针方向排列 在中 图8 1 任意一点的位置 可用它在直角坐标系中的两个坐标值来确定 图 8 1 如果我们过点作与三个顶点的连线 形成三个小三角形那么一旦的值确定后 这三个三角形的面积也就有了确定的值 反之 这三个小三角形的面积确定之后 点也就有了确定的位置 由此可见 三角形单元中任意一点的位置 除了可用直角坐标来确定外 还可以用连接点与与三个节点所形成的小三角形的面积来确定 用分别表示这三个小三角形的面积 显然 令 8 4 则称这三个比值为点的面积坐标 由定义可知 所以 并不是互相独立的 其中任意一个面积坐标都可以用另外两个面积坐标来表示 而且它们与直角坐标系的选取方法是无关的 这也是采用面积坐标的一个优点 显然 三个节点的面积坐标是节点节点节点 2 面积坐标与直角坐标的关系我们知道于是 有其中 由此可得到面积坐标与直角坐标之间的如下转换关系 8 6 8 7 从上述关系中可以看出 面积坐标和直角坐标之间是线性变换的关系 它实际上是将平面上的任意形状的三角形变换到平面上的直角三角形单元 经过这种变换 使得在任意三角形区域上的积分问题转化为在直角边为1的直角三角形区域上的积分问题 所以在计算上会带来很大的方便 3 面积坐标函数对直角坐标的偏导数设面积坐标函数为是的函数 由复合函数的求导法则 有注意到式 8 6 可得 所以面积坐标函数对直角坐标的偏导数为 8 8 4 面积坐标的积分单元分析中的积分 由于基函数几乎无例外地均采用多项式函数 被积函数一般都是以幂函数形式出现的 因此在单元分析中经常考虑的是如下典型形式的积分其中是任意非负整数 以面积坐标替换直角坐标 并利用重积分变量替换公式 不难算出 8 9 为证明此式 只需注意到变换的Jacobi行列式以及积分关系式 8 10 容易看出 在一维有限元分析中 由式 7 8 给出的变换正与上面讨论的面积坐标相当 如果说变换 8 6 将平面上的任意形状的三角形单元变换到平面上直角边为1的直角三角形单元 那么变换 7 8 则把轴上的线段单元变换到轴上的参考单元 0 1 由此可见 在一维的情况下 与公式 8 9 相应的积分公式为 8 11 其中 为线段单元的长度 2 构造单元基函数任一个三角形单元上 可唯一确定一个线性插值函数其中 为三角形单元顶点处给定的函数值 为处相应的单元节点基函数 它们都是线性函数 且满足条件根据前面的讨论容易验证 即面积坐标正好是三角形单元上线性插值函数的基函数 于是在任一三角形单元上 8 12 称式 8 12 为单元形状函数 将每一个单元上构造的函数合并恰来就得到在整个区域上的分块近似函数 由于每个节点对应一个基函数 所以整个区域上共有个基函数 容易证明 由所生成的试探函数空间是的维子空间 中的任一函数均可表为其中 是在节点处的值 8 3单元分析令 8 13 下面利用式 8 6 8 8 对式 8 13 右端的积分逐项进行分析 设单元为在三个顶点上 的值分别为和记为的面积 则在单元上于是式 8 13 右端的第一个积分可化为 其中 8 14 称为单元刚度矩阵 同理 可将式 8 13 右端的第二个积分化为 其中 8 15 称为单元载荷向量 综合以上分析 我们得到8 4总体合成令 8 16 我们知道 单元刚度矩阵或单元载荷向量中的元素下标值对应于单元节点序号 所谓总体合成 就是将这些序号转 换为总体节点序号 然后把这个元素加到总体系数矩阵某个位置上去 这个位置的行列序号 正是相应的总体节点序号 必须指出 对于每个单元中的节点号统一按逆 或顺 时针方向排列 但是它们在总的节点编号中就不一点按原来的次序排列了 设已经给出了单元节点编号与总的节点编号的对应关系 令对单元有其中 是一个的矩阵 若单元的节点在总节点编号中的序号为 则的第一行第个元素为1 其余元素均为0 的第二行 第三行分别也只有一个非零元素 其值为1 其位置由单元的节点在总编号中的序号来决定 同理 可以写出于是式 8 16 右端的第一项成为其中 8 17 就是总刚度矩阵 表示对个单元求和 显然 只是单元刚度矩的九个元素在总节点编号下重新排列和 扩展 的结果 而总刚度矩阵则是将各个单元的 贡献 叠加起来 它是一个的对称 正定矩阵 对于式 8 16 右端的第二项 同样得其中 8 18 就是总载荷向量 它是每个单元上单元载荷向量贡献的叠加 这样由式 8 3 可知 对 有 故得到所满足的线性代数方程组 8 19 方程组 8 19 的系数矩阵对称 正定 故方程组 8 19 有唯一解求得它们后 就有8 5边界条件的处理1 第一边值条件若第一边值条件为非齐次的 8 20 则应像内点一样 在界点也引进基函数 引进的方法及计算公式同内点完全一样 假定内节点和边界点的总个数为 则与上面的推导完全一样 可得到所满足的线性代数方程组或写成 8 21 要得到问题 8 1 8 20 的解在内部节点上的近似值 则可按如下方法对方程组 8 21 进行处理 假设有个边界节点 个内节点 为叙述方便 假定在总体节点编号时 把边界上的节点排在最前面 于是为了得到内点所满足的有限元方程 可先从方程组 8 21 中去掉前个方程 即 8 22 然后用边值代入左端相应项 并移至右端 便得到有限元方程 它是一个具个未知数 个方程的线性代数方程组若记其中 是的方阵 是的方阵 都是的 都是的 则方程 8 23 可写成如下形式 8 24 其中是从中划去头行列元素而得到的 显然 若边界条件是齐次的 则此方程组就是式 8 19 若边界上的节点不是排在前个 则在中划去相应的行列 在中划去相应的行 然后用边值代入左端相应项 并移至右端 同样得到内点所满足的线性代数方程组 以上的约束处理方法 将改为时 要重新存储矩阵的元素 从而给编制程序带来麻烦 因此在实际计算时 应把方程组 8 24 改写为 8 25 其中 是阶单位矩阵 是由边界上节点的值所组成 当边界条件为齐次时 这一方程组的系数矩阵保留了的阶数 它与方程组 8 24 等价 2 第二和第三边值条件假定问题中给定的是下