第一章有限元概念 1.1引言 有限元法已成为数值分析中一种实用而又重要的工具,一切 工程领域和应用数学领域几乎都在使用有限元法.论述有限元法 的文献数量极多,而且还在迅速增加。虽有大量文献目录可供查 阅,5,但这些目录仍未包括全部资料,而且这些文献很快就过时 了,尽管有许多现成的讲述各种有限元理论的教科书,但它们中 的大多数都把程序设计这一问题降低到次要的地位,只有 Hinton 和Owen2的教科书例外。本书所持观点与他们二位的相似,即 目的在于对典型的程序设计问题加以全面的、完整的论述,同时把 理论方面的内容限制在最小的程度。本书列举了一些实施有限元 的方法,选用了一些实例,这两方面的理论基础理所当然地是本书 所要提供的。 本章先介绍有限元法中所用的各种典型的积木式程序。这种 模块式程序可用于许多研究领域。以后几章再讨论应用有限元法 的一些具体例子 2有限元法的基础 从数学观点来看,有限元法是建立在积分表达法的基础上的 相比之下,早先的有限差分法则通常是建立在微分表达法的基础 上的.各种问题的有限元模型是根据简单的物理直觉和数学原理 提出来的、历史上,曾利用物理直觉引出一些早期的实用模型,但 目前却更加强调业已公认的该方法的数学基础“.礁61 最初,有元法在数学上不够严谨,而现在,该领域的研究非
常活跃.近代有限元的积分表达法是通过两种不同途径获得的,即 变分法和加权残数法。以下几节将简要地复习一下建立有限元模 型的一般方法。很凑巧,听有这些方法都使丌同一种簿记运算来 产生最后的代数方程组,必须通过求解这些方程式,才能得出待求 的结点参数 有限元模型最早的数学表达式是建立在变分法的基础上。变 分法在发展单元和在解决实际问题方面仍然非常重要,在结构力 学和应力分析领域里尤其如此,这两个领域里的近代分析方法几 乎全靠有限元法。变分法的模型通常是要找出一组结点参数值, 它使某一特定的称作泛函的积分式具有驻点值(即泛函取极大或 极小值)。在大多数情况下,这个取极值的积分式是有物理意义 的.例如,在固体力学中,该积分式代表位能,而在流体力学中则 可能对应于熵的产生率,许多物理问题的变分式都归结成二次型, 从这些二次型又进而得出一组对称正定的代数方程,以用来解决 该物理问题.变分表达式还有另一个重要的、有实际意义的优越 性,那就是变分法往往有一个与之相关的误差范围理论。翻阅一 下现有的许多论述变分法的教科书,就可以找到许多月变分法求 有限元模型的例子.本书后面几章将举例说明这种类型的应用。 众所周知,使泛函积分具有驻值且满足边界条件的解,等价于 对应的微分方程的解,这个微分方程就是 Euler方程,如果泛函 已知,则可比较容易地求出对应的 Euler方程.但大多数工程问 题和物理闰题最早是用微分方酲来确定的。有限元法需要的却是 种积分表达法。因此,必须找出一个泛函,使它的 Euler方程 恰好对应于已知的微分方程(及边界条件)。很不凑巧,一般说来, 这个工作很难,或者是不可能做到的,所以各种各样的加权残数 法越来越受到重视,因为这种方法可以直接从原来的微分方程组 产生出积分表达式。微分方程和积分形式都用物理坐标,比方说 (x:y,z)来确定 现在举一个简单的一维例子来说明积分表达法。考察一个将 在以后的应用实例中加以讨论的泛函
1-1K(d/dx)2+p1d 使此泛函取极小值就相当于满足微分方程 Kdi/drz- hi a 0 另外,此泛函还满足某端点处的自然边界条件d/dx=0,如果在 此端点处并未规定强制边界条件r=的话 利用加权残数法产生有限元模型是一种比较新的方法然而, 在求解微分方程和其他非结构力学问题方面,这种方法也越来越 重要了。加权残数法从基本微分方程 L(中)〓Q 出发,不去寻求在数学上等价的变分式,囡为变分式往往很难找 到。通常是先假设一个近似解,比如φ*,将其代人微分方程。既 然是假没的近似解,这一代入就在微分方程里产住一个残数误差 项 L(中*)-QR 虽然不可能强迫残差项等于零,但却可令残差项在求解域内 的加权积分等于零,所以得到 RW do ed 用插值函数代换近似解d*和加权函数W,就得到一组代 数方程式,解出后可得近似解中的那些未知系数。选用什么样 加权函数决定于所使用加权残数法的类型。要得到 Galerkin判 据,可选 而对于最小二乘方判据, WR/θ小* 給出所需结果。与此类似,选用 Dirac delta函数所得出来的是 种点配置法,即 W=δ 显然,还有其他W可供选用,从而引出其他的加权残数法,如子
域法,前两种方法似乎是有限元法中最常用的.对 Galerkin法 使用部分积分通常会降低对近似函数的连续性要求。如果存在 个变分法,那么从 Galerkin判据将得出与之相同的代数近似解, 因此常可给出有限元解的最佳误差估计值 目前,变分法和加权残数法两者都普遍采用下列限制条件来 作为使有限元模型随网格密度增加而具有收敛性的一种手段7: 1(一种必要判据)当单元尺寸减小而处于极限状态时,单元 的捕值函数必须能模拟因变量或其导数取任意常数,而导数的阶 次高到定义积分式中出现的阶次, 2.(一种充分判据)选择单元的形状函数应能使在单元的交接 面处,因变量及其导数是连续的,而导数的阶次比定义积分式中出 现的阶次低一阶 13通用的有限元分析法 131引 在有限元法里,利用假想的线(或面),把连续介质或更一般 地说是把求解域)的内部和边界分割成有限个数目的、离散的,有 定大小的子域,这些子域就是有服元(见图1.1).假想的网格将 区域进行分割时就构成了一些离散结点。这些结点的位置可以 在网格分格线上的任何地方,也可以在两格分格线内的任何地方, 但通常是在闋格的相交线(或相交面)上,一殷说来,单元具有直 线边界;如果所研究的区域具有曲线边界,则在理想化时就会引人 几何形状方面的假设解 各结点都被赋予一个起识别作用的整数号码结点号),从1 开始增到某一最大值,如M,与此相似,每一个单元也被赋予 个起识别作用的整数号.这些单元号码也是从1开始扩大到某个 最大值,如NE,以后将要讨论,结点号和单元号的赋值方式可能 会对求解所需时间和存储所需容量产生重大影响,工作人员对每 个结点都要规定具有一定数量的(广义的)自由度(dof),如NG
这些白由度就是(待求的)结点参数,工作人员就是选择这些结点 参数来控制如何提出所研究的问题,常见的结点参数有压强、速 度分量位移分量、位移梯度等。结点参数不一定有物理意义,不 过通常是有的。本书将假设系绕中每个结点的结点参数的个数相 同,即结点参数的数目(NG)是一样的,这是通常情况,但并非一 定如此,如图1.1所示,与 个典型结点相关联的往往不止 个典型结点的影响城 一个单元。图1.1也表示一个 典型结点和一个典型单元的影 响域。一个典型单元有好几个 结点与之相联,如N个,而点的 位置则在边界上或在边界内 此处还假设每个单元的单元结 点数(N)是一样的,通常总是 这样假设,但并非所有情况都 必须如此 这种理想化处理就确定了 一个典型单元的影啊 每个典型结点和典型单元的总 自由度数。显然,整个系统的 图1·1影响域 自由度数(若用 NDFREE表示)就是结点数与每一结点参数数 的乘积即 NDFREE一M*NG,与此类似,每个单元的自由度 数(若用 NELFRE表示),就可按 NELFRE一N*NG来确定 应记住:系统自由度的总数对应于结点参数的总数。一般说 来,系统中第1个结点第J个参数所对应的系统自由度号NDF 可(根据归纳推理)由下式确定: NDF=NG*(I-1)十J (L1) 式中1≤I≤M及1≤J≤NG,所以I≤NDF≤ NDFREE 这个基本方程式就是对程序进行“簿记”法的基础,因而这个方程 式非常重要,应该对其理解清楚。方程式(1.1)可用图1.2中一串 一元加以说明,该图所不系统有4个线段元,5个结点,每个结
点有6个自由度(即M=5,NE冖4,NGm6,N一2)。系统 总共有30个自由度,要确定此系统4号结点(I-4)第5个参数 的自由度号码、利用式(1.1),所得结果为NDF〓6(4-1)+5 23,这就是所对应的系统自由度号.对于第3个单元来说,系统 自由度号23对应的是局部自由度号11而对第4个单元来说,对 应的是局部自由度号5。因此,我们可以看出有两个不同的单元 方程组的部分对第23号系统方程式有贡献。 封点号 单元号 2 2 14 26 站点自由度号3 2了 10 16 22 28 5 12 图12目由度号的计算 除了上述这些常数外,还需确定该课题的空间维数,例如用 SPACE表示.从讨论过程中可以看出,利用这些参量来计算在 分析中所生成的矩阵需要多大的储存量。请读入课题数据的实际程 序将在以后讨论 必须提供数据以便确定每一结点的空间坐标位置。这个数组 可用X表示,其维数为M* N SPACE.通常,给每一个结点都编 一个整数号码。这个代码的目的是指出:如果有边界约束条件的 话,该结点哪一个参数具有指定的边界约束,这个数据向量可用 IBC表示,它包含M个整数系数。为完成这一结点边界条件代码 过程,应记住每一结点有NG个参数。因此,可以确定一个(右侧 对齐的)整数代码IBC,包含的位数是NG.令第i位是一个一位 数的指示码,对应该结点的第i个参数.如果这个指示码等于j 0≤≤9,则规定其意义为第i个参数受到j型边界约束。若该 位数指示码为零,就表示在该参数上没有约束条件作用,以后
还会讨论,这种程序可以考虑到好几种共同类型的结点参数边界 约束.图1.3举例说明一套边界条件代码,对应于一组有代表性的 结点,而每一结点有6个参数(NG一6).这种代码法在73节 的例题中及全部应用例题中都将加以使用 结点参数号3的边界码 结点BC了=1 111111 3220000 102 310 4-100 4 100 210201 00I 51000d 5110000 图1.3结点边界码的解释 草元的关联是一个重要概念,所谓关联就是把属于该单元的 总体结点号码排列成表。单元的关联数据确定阙格的拓扑结构 因此,针对每个单元,都需要(以某种前后一贯的顺序)输入与该单 元有关的N个结点号,这个数据所构成的数组可用 nODES表 示,其维数为NE*N.把与某一特定单元相联的(总体)结点号 列成表,这种表通常称为单无关联结点表。识别这套数据对使用 方程(1.1)很重要 L32单元行为和方程式的近似表达 假设整个求解域上的待求变量(也许)连同其导数一起,都 可以由系统结点的结点参数唯一地确定,那么这些结点参数就是 该问题的待求参数,假设一个结点的各参数只影响到与该结点相 联的诸单元内的待求量值。再假设一个插值函数,以便能用与单 元相联的诸结点的结点参数值来表达该单元内任一点的待求量 值.图1.4举例说明一种常用插值函数及其组成部分,此函数由 结点坐标(x,y)单元的面积∥、空间位置坐标(xy)及结点 参数T,T;T+确定 对单元的行为作出假设后,剩下的基本问题就是建立单元矩 阵S和C,一般说来,就是要把插值函数代人基本积分式中
历史上,这两个矩阵一直被分別称为单元刚度姬阵和单元載荷向 量.虽然这两个矩阵有时可按物理直觉推导出来,但通常却用求 泛函极小值的方法得出。或用加仪残数法得出。这些方法在许多 教材中都可见到,本书后面几章将详细加以说明。 点坐标 〔期 t(r,) b;=y“yk (x;y) 2A=aitait d 假设的多项式 结点参数 a r(,y)=u1xy b; 插值所得参量 单元的几何常数 r(,y)=HR(,y)T+Hj(r,y)r+ H(r,D)T r(x,y)"(,y 图1.4常用插值分量 建立单元矩阵,几乎仝要涉及以某种方式给出该单元的特性 或特性系数。有少数有限元问题要求给出结点处的特性。例如 在应力分析问题中,希望能通过指定材料在各结点处的厚度来确 定变厚度元任一点的厚度)有限元法非常适用于非均质问题,所 以大多数有限元程序也要求分析人员指定每个单元的某些特性 一般来说,对于各单元(或各结点)所共有的数据,最好是把它作为 整个系统的各种数据加以存钴.分析人员必须决定哪些数据是需 要的,采用哪种方式输入和取用这些数据最好。 1.33方程式的组装及求解 且建立单元的方程式后,就要把各个单元的贡献加起来组 成整个系统的方程式。组装过程约详纽编程将在72节讨论,有 限元分析得出的系统方程一般具有以下形式
SD=C (12) 其中方阵S的维数是 NDFREE* NDFREE,而向量D和C各 含 NDFREE个系数。向量D包含未知结点参数,而矩阵S C(如后所述)分别由单元矩阵s和C的组装得出,S的维 数是 NELFRE米 NE LFRE,C的维数是 NELFrE*1,大多数 问题的S是对称的,从而S也是对称的。还有,系统矩阵S 般是围绕对角线成带状的。今后假设系统方程组是带状的、对称 的。因而,包括对角线在内的半带宽的大小,是任一有限元法都应 加以考虑的一个重要参量.考虑到在近似表达单元行为时所用的 方法可知:若令由于某一典型单元“e”而使系统方程组产生的半 带宽为IBW,则IBW可由下式确定 IBW=NG水( NDIFF+1) (1.3) 式中 NDIFF是该单元诸结点的结点号码的最大差值的绝对值, 方程(13)将在以后推导.系统的最大半带宽,也就是系统中所存 在的IBW中的最大值,可用 MAXBAN表示,即 MAXBAN IBw (14) 这个参量 MAXBAN是决定系统方程组的存储容量和求辉机时的 重要参量之一。因此,虽然结点号码的赋值可以任意,但实际上, 分析人员应尽力缩小典型单元上各结点编号之间的最大差值(和 带宽).图15说明组装过程,所举例子是一个4单元络,由3结 点三角形组成,而每一结点只有一个参数.图的上方说明系统的 s矩阵和C矩阵的组装,该图用各种阴影线作标记,这种标记用 来指出影响该项的是哪些单元,而不代表该项的数值大小。矩阵 中有许多不同的阴影线区域,与此对应,每个单元也用阴影线作了 标记(见图的下方),矩阵中每有一块阴影区,就表示该项加进了对 应单元的作用。例如,可以看出,载荷向量第C(6)项是把单元2 和单元3的贡献加在一起得出来的,因为单元2的标记是水平线 (-),而单元3的标记是斜线(/).相比之下,第C(1)项只有单 元2的作用。 最后应注意,某些卓有成效的有限元程序并不利用带状矩阵
(a) 图1·5单元对组装的贡献 ()系统矩阵;(b)有贡献的单元 求解法,取而代之的是,运用波前法或稀疏矩阵求解,这些重要课 题将在后面较详细地讨论 组装成系统方程组后,必须在求解未知结点参数之前先施加 边界约束条件。最常见的结点参数边界约束如下