1.4整体分析 13 每个节点都可列出如上所述的一组平衡方程。如果有限元计算模型共有 N个节点，则可得到2N阶线性方程组 Ka=P 1.21) 此即整体刚度方程。其中K为整体刚度矩阵：α和P分别为整体节点位移向量 和整体节点荷载向量，即 a11 P 1.4.2虚功原理 也可以采用虚位移原理来建立整体刚度方程。为此，将虚位移原理应用于 整个结构，则整体节点荷载所做虚功等于所有单元的虚应变能之和，即 dn-prdn+r) 将虚位移以及虚应变与虚节点位移之间的关系代入上式得 ☒aBad=yar（NTan+Npar) 将单元节点位移向量a用整体节点位移向量a来表示，即 a°=Aa 其中，A°为单元连接矩阵。于是，得到整体平衡方程 ∑AT,BTodQ=∑AT(+P) (1.22a) 将式1.15)代入上式，可得到式1.21)，且 K=∑ATKA,P=∑A(P+P) (1.22b) 为简便起见，上式通常写成 K=K,P=(P+P) (1.22c 此时的求和号表示集成，而非简单相加。 必须指出，在前面的分析中引入节点力并不意味着近似，因为在应用虚位移 原理进行整体分析时，根本涉及不到单元节点力问题。仅从整体分析考虑，内部
每个节点都可列出如上所述的一组平衡方程。如果有限元计算模型共有 Ｎ 个节点，则可得到２Ｎ 阶线性方程组 Ｋａ＝Ｐ （１２１） 此即整体刚度方程。其中 Ｋ 为整体刚度矩阵；ａ和Ｐ 分别为整体节点位移向量 和整体节点荷载向量，即 ａ＝ ａ１ … ａ 烅 烄 烆 烍 烌 Ｎ烎 ， Ｐ＝ Ｐ１ … Ｐ 烅 烄 烆 烍 烌 Ｎ烎 １４２ 虚功原理 也可以采用虚位移原理来建立整体刚度方程。为此，将虚位移原理应用于 整个结构，则整体节点荷载所做虚功等于所有单元的虚应变能之和，即 ∑ｅ ∫Ω ｅδεＴσｄΩ ＝ ∑ｅ ∫Ω ｅδｕＴ ｆｄΩ＋∫Γ ｅ σ （ ） δｕＴ珔ｐｄΓ 将虚位移以及虚应变与虚节点位移之间的关系代入上式得 ∑ｅ ∫Ω ｅδａｅＴＢＴσｄΩ ＝ ∑ｅ δａｅＴ ∫Ω ｅＮＴ ｆｄΩ＋∫Γ ｅ σ （ ） ＮＴ珔ｐｄΓ 将单元节点位移向量ａｅ 用整体节点位移向量ａ来表示，即 ａｅ＝Ａｅａ 其中，Ａｅ 为单元连接矩阵。于是，得到整体平衡方程 ∑ｅ Ａｅ ∫Ｔ Ω ｅＢＴσｄΩ ＝ ∑ｅ ＡｅＴ（Ｐｅ ｆ＋Ｐｅ ｐ） （１２２ａ） 将式（１１５）代入上式，可得到式（１２１），且 Ｋ ＝ ∑ｅ ＡｅＴＫｅＡｅ， Ｐ ＝ ∑ｅ ＡｅＴ（Ｐｅ ｆ＋Ｐｅ ｐ） （１２２ｂ） 为简便起见，上式通常写成 Ｋ ＝ ∑ｅ Ｋｅ， Ｐ ＝ ∑ｅ （Ｐｅ ｆ＋Ｐｅ ｐ） （１２２ｃ） 此时的求和号表示集成，而非简单相加。 必须指出，在前面的分析中引入节点力并不意味着近似，因为在应用虚位移 原理进行整体分析时，根本涉及不到单元节点力问题。仅从整体分析考虑，内部 １４ 整 体 分 析 ３１
第1章有限单元法基本程式 等效节点力的概念似乎可以放弃。但是，在有限单元法中，定义节点力是有意义 的，因为从功的意义上它们与节点位移共轭。为使节点力与节点位移之积给出 功的正确表达式，规定节点力的正向与节点位移的正向一致。 1.4.3直接集成 在实际有限元分析中，建立整体刚度矩阵最常用的方法是直接集成法或直 接刚度法，即直接由单元刚度矩阵集合而成，其关键是把所有单元刚度矩阵的各 元素安放到K中的适当位置。具体地说，就是将单元刚度矩阵K扩大成单元 的贡献矩阵K:然后将各单元的K:直接相加得出K。 1)分块集成 对于三角形单元，分块集成时需把K中的6个元素搬家，按照整体码的顺 序在扩大后的矩阵中重新排列，并在空白处用零元素填补起来。一般地说，若单 元e的局部码1,2,3分别对应于整体码I,J,M,且设I<J<M,则该单元的贡 献矩阵如式(1.23)所示。 整体码 1…I…J…M…N 1「0… 0 … 0 … 0 … 07 … … I0…k11… k12 …k13 … 01 …… … K=J0…k21…k2 …k23 …02 1.23) … M0…k31 …k2… k3 …03 …… … … NL0…0 … 0 0 …0」 1 2 3 局部码 ②)元素集成 在实际计算中，需要按元素形式集成整体刚度矩阵。如果单元的局部码 1,2,3分别对应于整体码1，J,M,则该单元节点自由度编码如表1.1所示。 表1.1自由度局部码与整体码 节点自由度 41 42 2 自由度局部码 2 3 4 5 6 自由度整体码 21-1 21 2-1 2J 2M-12M
等效节点力的概念似乎可以放弃。但是，在有限单元法中，定义节点力是有意义 的，因为从功的意义上它们与节点位移共轭。为使节点力与节点位移之积给出 功的正确表达式，规定节点力的正向与节点位移的正向一致。 １４３ 直接集成 在实际有限元分析中，建立整体刚度矩阵最常用的方法是直接集成法或直 接刚度法，即直接由单元刚度矩阵集合而成，其关键是把所有单元刚度矩阵的各 元素安放到 Ｋ 中的适当位置。具体地说，就是将单元刚度矩阵 Ｋｅ 扩大成单元 的贡献矩阵Ｋｅ ｃ；然后将各单元的 Ｋｅ ｃ直接相加得出 Ｋ。 （１）分块集成 对于三角形单元，分块集成时需把 Ｋｅ 中的６个元素搬家，按照整体码的顺 序在扩大后的矩阵中重新排列，并在空白处用零元素填补起来。一般地说，若单 元ｅ的局部码１，２，３分别对应于整体码Ｉ，Ｊ，Ｍ，且设Ｉ＜Ｊ＜Ｍ，则该单元的贡 献矩阵如式（１２３）所示。 整体码 １ … Ｉ … Ｊ … Ｍ … Ｎ Ｋｅ ｃ＝ １ … Ｉ … Ｊ … Ｍ … Ｎ ０ … ０ … ０ … ０ … ０ … … … … … ０ … ｋ１１ … ｋ１２ … ｋ１３ … ０ … … … … … ０ … ｋ２１ … ｋ２２ … ｋ２３ … ０ … … … … … ０ … ｋ３１ … ｋ３２ … ｋ３３ … ０ … … … … … ０ … ０ … ０ … ０ … 熿 燀 燄 ０燅 １ ２ ３ （１２３） １ ２ ３ 局部码 （２）元素集成 在实际计算中，需要按元素形式集成整体刚度矩阵。如果单元ｅ的局部码 １，２，３分别对应于整体码Ｉ，Ｊ，Ｍ，则该单元节点自由度编码如表１１所示。 表１１ 自由度局部码与整体码 节点自由度 ｕ１ ｖ１ ｕ２ ｖ２ ｕ３ ｖ３ 自由度局部码 １ ２ ３ ４ ５ ６ 自由度整体码 ２Ｉ－１ ２Ｉ ２Ｊ－１ ２Ｊ ２Ｍ－１ ２Ｍ ４１ 第１章 有限单元法基本程式
1.4整体分析 15 利用上述编码表，不难确定单元刚度系数与整体刚度系数的对应关系。例如对 于单元刚度系数k5,从编码表第2格取出21，从第5格取出2M-1,则对应的 整体刚度系数为K2,2M-1,即 k25→K21.2M-1 同理 k11→K21-1,21-1 把全部单元的刚度系数都按照编码表叠加到相应的整体刚度系数中去，就 可得到整体刚度矩阵。 3)荷载集成 荷载集成的任务是形成整体节点荷载向量P。总的节点荷载包括集中力和 体力、面力移置而成的等效荷载。集中力作用点通常都被取做节点，因此只要将 给定的集中力直接送入P中的适当位置即可。例如在整体码为I的节点上沿x 方向作用集中力Q,则该集中力在P中的位置为2I-1。 体力和面力按照等效的原则化成节点荷载。一般是先按单元移置，以后按 节点叠加。例如，设单元e的体力等效节点荷载向量已经得到，记为 P=[Fir Fiy F2r F2y F3r F3y]T 若该单元节点局部码1,2,3对应的整体码为1，J,M,则F2,在P中的位置为 2J-1。 1.4.4K的性质 )K的意义 整体刚度矩阵K中每一列元素的物理意义是：要迫使结构的某节点位移自 由度发生单位位移，而其他节点位移都保持为零的变形状态，在所有各节点上需 要施加的节点荷载。例如，令1=1，其余的节点位移均为零，则式(1.21)的节 点荷载向量显然等于K的第一列元素的列阵。从物理上说，K之对角线上的 主元素K:总是正的，否则作用力的方向将与它引起的对应位移的方向相反。 2)对称性 根据功的互等定理，可知K是对称矩阵。由于K是对称的，故在实际计算 中只形成并存贮上三角阵或下三角阵即可。 3)奇异性 从物理上讲，当结构的几何约束尚未设置、刚体位移未被排除之前，不可能 有唯一的位移解。这个物理事实在数学上表现为K的奇异性，即其行列式的值
利用上述编码表，不难确定单元刚度系数与整体刚度系数的对应关系。例如对 于单元刚度系数ｋ２５，从编码表第２格取出２Ｉ，从第５格取出２Ｍ－１，则对应的 整体刚度系数为 Ｋ２Ｉ，２Ｍ－１，即 ｋ２５→Ｋ２Ｉ，２Ｍ－１ 同理 ｋ１１→Ｋ２Ｉ－１，２Ｉ－１ … 把全部单元的刚度系数都按照编码表叠加到相应的整体刚度系数中去，就 可得到整体刚度矩阵。 （３）荷载集成 荷载集成的任务是形成整体节点荷载向量Ｐ。总的节点荷载包括集中力和 体力、面力移置而成的等效荷载。集中力作用点通常都被取做节点，因此只要将 给定的集中力直接送入Ｐ 中的适当位置即可。例如在整体码为Ｉ的节点上沿ｘ 方向作用集中力 Ｑ，则该集中力在Ｐ 中的位置为２Ｉ－１。 体力和面力按照等效的原则化成节点荷载。一般是先按单元移置，以后按 节点叠加。例如，设单元ｅ的体力等效节点荷载向量已经得到，记为 Ｐｅ ｆ＝［ ］ Ｆ１ｘ Ｆ１ｙ Ｆ２ｘ Ｆ２ｙ Ｆ３ｘ Ｆ３ｙ Ｔ 若该单元节点局部码１，２，３对应的整体码为Ｉ，Ｊ，Ｍ，则Ｆ２ｘ在Ｐ 中的位置为 ２Ｊ－１。 １４４ Ｋ 的性质 （１）Ｋｉｊ的意义 整体刚度矩阵 Ｋ 中每一列元素的物理意义是：要迫使结构的某节点位移自 由度发生单位位移，而其他节点位移都保持为零的变形状态，在所有各节点上需 要施加的节点荷载。例如，令ｕ１＝１，其余的节点位移均为零，则式（１２１）的节 点荷载向量显然等于 Ｋ 的第一列元素的列阵。从物理上说，Ｋ 之对角线上的 主元素Ｋｉｉ总是正的，否则作用力的方向将与它引起的对应位移的方向相反。 （２）对称性 根据功的互等定理，可知 Ｋ 是对称矩阵。由于 Ｋ 是对称的，故在实际计算 中只形成并存贮上三角阵或下三角阵即可。 （３）奇异性 从物理上讲，当结构的几何约束尚未设置、刚体位移未被排除之前，不可能 有唯一的位移解。这个物理事实在数学上表现为 Ｋ 的奇异性，即其行列式的值 １４ 整 体 分 析 ５１
16 第1章有限单元法基本程式 等于零。 ④)稀疏性 K是稀疏矩阵，且划分的单元越多越稀疏。如果节点编号恰当，那些非零 元素都将集中于刚度矩阵的对角线附近，呈斜带状（图1.7）。显然，带外的零元 素不必存贮，也不参加运算，这是有限元整体刚度矩阵的优良数值特性。 在整体刚度矩阵的各行中，由对角线到带边界所包 含最多的元素数目称为半带宽(Semi-bandwidth),用B 表示。若单元各节点的整体码I,J,M等非常接近，则 该单元的刚度系数便集中在K的对角线附近。不难发 现，半带宽决定于各单元中节点整体码的最大差值D。 如果每个节点的自由度为，则有 0 +.子半带 B=n(D+1) (1.24) 在平面问题中，n=2。 1.5数值求解 1.5.1边界条件的引入 如前所述，整体刚度矩阵K是奇异的，只有引入位移边界条件对刚度矩阵 加以修改，即消除刚体位移后才能求解整体刚度方程。对于平面问题来说，要消 除刚体位移，至少要有三个位移约束条件。 1)零位移的实现 对于具有N个节点的平面结构，其平衡方程为 [KL K1,2K13 K14 K1,2N 11 (XI K21 K2,2 K23 K2,4 K2,2N U Y K3.1 K3.2K3.3 K3.4 … K3,2N 1.25) K41 K42K43 K44 K4.2 V2 Y2 … K2N.I K2N.2 K2N.3 K2N,4 … K2N.2NJ UN) 例如，为实现1=0的条件，在式1,25)中可作如下变化：在整体刚度矩阵K 中，除了保留与1相对应的并在主对角线上的系数K.,外，第一行和第一列的 其余系数均改为零：在荷载列阵中，令u1对应的X1=0,从而有
等于零。 （４）稀疏性 Ｋ 是稀疏矩阵，且划分的单元越多越稀疏。如果节点编号恰当，那些非零 元素都将集中于刚度矩阵的对角线附近，呈斜带状（图１７）。显然，带外的零元 素不必存贮，也不参加运算，这是有限元整体刚度矩阵的优良数值特性。 图１７ 半带宽 在整体刚度矩阵的各行中，由对角线到带边界所包 含最多的元素数目称为半带宽（Ｓｅｍｉ?ｂａｎｄｗｉｄｔｈ），用 Ｂ 表示。若单元各节点的整体码Ｉ，Ｊ，Ｍ 等非常接近，则 该单元的刚度系数便集中在 Ｋ 的对角线附近。不难发 现，半带宽决定于各单元中节点整体码的最大差值 Ｄ。 如果每个节点的自由度为ｎ，则有 Ｂ＝ｎ（Ｄ＋１） （１２４） 在平面问题中，ｎ＝２。 １５ 数 值 求 解 １５１ 边界条件的引入 如前所述，整体刚度矩阵 Ｋ 是奇异的，只有引入位移边界条件对刚度矩阵 加以修改，即消除刚体位移后才能求解整体刚度方程。对于平面问题来说，要消 除刚体位移，至少要有三个位移约束条件。 （１）零位移的实现 对于具有 Ｎ 个节点的平面结构，其平衡方程为 Ｋ１，１ Ｋ１，２ Ｋ１，３ Ｋ１，４ … Ｋ１，２Ｎ Ｋ２，１ Ｋ２，２ Ｋ２，３ Ｋ２，４ … Ｋ２，２Ｎ Ｋ３，１ Ｋ３，２ Ｋ３，３ Ｋ３，４ … Ｋ３，２Ｎ Ｋ４，１ Ｋ４，２ Ｋ４，３ Ｋ４，４ … Ｋ４，２Ｎ … … … … … … Ｋ２Ｎ，１ Ｋ２Ｎ，２ Ｋ２Ｎ，３ Ｋ２Ｎ，４ … Ｋ２Ｎ，２ 熿 燀 燄 Ｎ燅 ｕ１ ｖ１ ｕ２ ｖ２ … ｖ 烅 烄 烆 烍 烌 Ｎ烎 ＝ Ｘ１ Ｙ１ Ｘ２ Ｙ２ … Ｙ 烅 烄 烆 烍 烌 Ｎ烎 （１２５） 例如，为实现ｕ１＝０的条件，在式（１２５）中可作如下变化：在整体刚度矩阵 Ｋ 中，除了保留与ｕ１相对应的并在主对角线上的系数 Ｋ１，１外，第一行和第一列的 其余系数均改为零；在荷载列阵中，令ｕ１对应的Ｘ１＝０，从而有 ６１ 第１章 有限单元法基本程式
1.5数值求解 17 K10 0 0 0 (ur 0 0K2,2K2,3 K2,4 K2.2N 4 0 K32 K33 K34 K3.2N X2 0 K42 1.26) K4,3 K44 … K4.2 … 0 K2N,2K2N,3K2N.4 K2N.2N Y 按式1.26)求解就能实现1=0的条件，其他零位移条件的实现可类推。 ②)有限位移的实现 例如，为实现w2=b,可仿照上述做法将式1.25)改为 「K1 K1.2 0 K14 …K1,2N1u1 X1-K1,3b K21 K22 0 K2,4 … K2,2N 1 Y1-K2.3b 0 0 K33 0 0 u2 K3.3b K41K4,20 K4,4… KA.2N Y2-K4.3b LK2N.I K2N.2 0 K2N.4 K2N.2N YN-K2N.3b) 1.27) 不难发现，按1.27)求解就能实现2=b的条件。 然而，为程序设计方便起见，通常采用如下近似方法：把与42对应的对角 线上的刚度系数K3.3乘以一个大数α（例如100）：把u2对应的节点荷载换成 abK33,其余均保持不变，即把平衡方程改为 K1.2 K1,3 K1,4 K1.2N7fu1 K21 K22 K2,3 K24 … K2.2N Y K3.1K3.2aK3.3 K34 K3,2N abK3.3 1.28) K4 K4,2 K4,3 K4,4… K4.2N Y … LK2N.1 K2N.2 K2N.3 K2N.4 K2N.2NN 显然，由于αK3.3比其他刚度系数大得多，故第三个方程实际上可写成 aK33u2≈abK33 从而近似地实现了规定的边界条件u2=b。 1.5.2未知量的求解 引入边界条件消除奇异性后，便可求解整体刚度方程。线性代数方程组的
Ｋ１，１ ０ ０ ０ … ０ ０ Ｋ２，２ Ｋ２，３ Ｋ２，４ … Ｋ２，２Ｎ ０ Ｋ３，２ Ｋ３，３ Ｋ３，４ … Ｋ３，２Ｎ ０ Ｋ４，２ Ｋ４，３ Ｋ４，４ … Ｋ４，２Ｎ … … … … … … ０ Ｋ２Ｎ，２ Ｋ２Ｎ，３ Ｋ２Ｎ，４ … Ｋ２Ｎ，２ 熿 燀 燄 Ｎ 燅 ｕ１ ｖ１ ｕ２ ｖ２ … ｖ 烅 烄 烆 烍 烌 Ｎ烎 ＝ ０ Ｙ１ Ｘ２ Ｙ２ … Ｙ 烅 烄 烆 烍 烌 Ｎ烎 （１２６） 按式（１２６）求解就能实现ｕ１＝０的条件，其他零位移条件的实现可类推。 （２）有限位移的实现 例如，为实现ｕ２＝ｂ，可仿照上述做法将式（１２５）改为 Ｋ１，１ Ｋ１，２ ０ Ｋ１，４ … Ｋ１，２Ｎ Ｋ２，１ Ｋ２，２ ０ Ｋ２，４ … Ｋ２，２Ｎ ０ ０ Ｋ３，３ ０ … ０ Ｋ４，１ Ｋ４，２ ０ Ｋ４，４ … Ｋ４，２Ｎ … … … … … … Ｋ２Ｎ，１ Ｋ２Ｎ，２ ０ Ｋ２Ｎ，４ … Ｋ２Ｎ，２ 熿 燀 燄 Ｎ燅 ｕ１ ｖ１ ｕ２ ｖ２ … ｖ 烅 烄 烆 烍 烌 Ｎ烎 ＝ Ｘ１－Ｋ１，３ｂ Ｙ１－Ｋ２，３ｂ Ｋ３，３ｂ Ｙ２－Ｋ４，３ｂ … ＹＮ－Ｋ２Ｎ，３ 烅 烄 烆 烍 烌 ｂ烎 （１２７） 不难发现，按（１２７）求解就能实现ｕ２＝ｂ的条件。 然而，为程序设计方便起见，通常采用如下近似方法：把与ｕ２ 对应的对角 线上的刚度系数 Ｋ３，３乘以一个大数α（例如１０１０）；把ｕ２ 对应的节点荷载换成 αｂＫ３，３，其余均保持不变，即把平衡方程改为 Ｋ１，１ Ｋ１，２ Ｋ１，３ Ｋ１，４ … Ｋ１，２Ｎ Ｋ２，１ Ｋ２，２ Ｋ２，３ Ｋ２，４ … Ｋ２，２Ｎ Ｋ３，１ Ｋ３，２ αＫ３，３ Ｋ３，４ … Ｋ３，２Ｎ Ｋ４，１ Ｋ４，２ Ｋ４，３ Ｋ４，４ … Ｋ４，２Ｎ … … … … … … Ｋ２Ｎ，１ Ｋ２Ｎ，２ Ｋ２Ｎ，３ Ｋ２Ｎ，４ … Ｋ２Ｎ，２ 熿 燀 燄 Ｎ燅 ｕ１ ｖ１ ｕ２ ｖ２ … ｖ 烅 烄 烆 烍 烌 Ｎ烎 ＝ Ｘ１ Ｙ１ αｂＫ３，３ Ｙ２ … Ｙ 烅 烄 烆 烍 烌 Ｎ 烎 （１２８） 显然，由于αＫ３，３比其他刚度系数大得多，故第三个方程实际上可写成 αＫ３，３ｕ２≈αｂＫ３，３ 从而近似地实现了规定的边界条件ｕ２＝ｂ。 １５２ 未知量的求解 引入边界条件消除奇异性后，便可求解整体刚度方程。线性代数方程组的 １５ 数 值 求 解 ７１