3.2.1　基本概念与思路

（1）基本概念

有限元法（或称有限单元法、有限元素法）是复杂工程和产品结构强度、刚度、屈曲稳定性、动力响应、热传导、三维多体接触、弹塑性等力学性能的分析计算以及结构性能的优化设计等问题的一种近似数值分析方法。它的基本概念是将一个复杂的连续体的求解区域分解为有限个形状简单的子区域（单元），即将一个连续体简化为由有限个单元组成的等效组合体；通过将连续体离散化，把求解连续体的场变量（应力、位移、压力和温度等）问题简化为求解有限个单元节点上的场变量值。此时求解的基本方程将是一个代数方程组，而不是原来的描述真实连续体场变量的微分方程，得到近似的数值解。求解的近似程度取决于所用的单元类型、数量以及对单元的差值函数。

（2）基本思路

有限元法的思路是：将一个连续求解域（对象）离散（剖分）成有限个形状简单的子域（单元），利用有限个节点将各子域连接起来，使其分别承受相应的等效节点载荷（如应力载荷、热载荷、流速载荷等），并传递子域间的相互作用；在此基础上，借助子域插值函数和“平衡”条件构建各子域的物理场控制方程；将这些方程按照某种规则组合起来，在给定的初始条件和边界条件下进行综合计算求解，从而获得对复杂工程问题的近似数值解。其中，离散和子域插值是有限元法的技术基础。

图3-2是对离散概念的图解说明。离散求解域的目的是为了将原来具有无限自由度的连续变量微分方程和边界条件转换成只含有限个节点变量的代数方程组，以方便计算机处理。

分片插值的概念可以借助图3-3加以说明。假设真实函数为曲线 c ₁ ，求解域为[ a ， b ]。理论上讲，只要定义在[ a ， b ]上的试探函数（亦称插值函数） c ₂ 具有足够高的阶次就能逼近真实函数 c ₁ ，但实际上 c ₂ 对 c ₁ 局部特征的逼近并不理想。如果将求解域划分成若干长度不等的小区间，则可在每一个小区间内用较低阶（例如一阶或二阶）的试探函数 c ₃ 来逼近 c ₁ ，并且通过适当调整求解域局部小区间的数量或尺寸来提高逼近精度，从而获得真实函数 c ₁ 的近似解。

通常，将构建子域物理场方程的过程称为“单元分析”，将在初边值条件支持下综合求解子域方程组的过程称为“整体分析”，因此，有限元法的中心思想又可简略描述为：离散求解域→单元分析→整体分析。

利用子域（单元体）离散连续求解域（实体模型或对象）的过程又被形象地称为网格划分，由此得到的离散模型被称为网格模型。

（3）技术优势

有限元法的技术优势主要体现在：该方法把连续体简化成由有限个单元组成的等效体（物理上的简化），针对等效体建立的基本方程是一组代数方程，而不是原先用于描述真实连续体的常微分或偏微分方程。由于不存在数学上的近似，故有限元法的物理概念清晰，通用性强，能够灵活处理各种复杂的工程问题。

3.2.2　有限元方程的建立与应用

针对不同工程问题构建的有限元方程（有限元模型）是有限元法应用的基础，而变分法和加权余量法是建立有限元方程的两类常用数学方法。本节将以求解平面弹性力学刚度问题为例，分别介绍这两类方法的实施要点。

3.2.2.1　预备知识

（1）弹性力学基本方程

①平衡方程。当弹性体中任一质点上的应力达到平衡时，有：

（3-1）

式中 L ——微分算子；

σ ——应力；

p ——体积力（一般为重力）向量。

对于平面问题有 ， ， 。

②几何方程。几何方程表征质点位移与应变之间的关系。

（3-2）

式中 ε ——应变；

u ——质点位移矢量。

对于平面问题有 ε ^T ={ ε _x ε _y γ _xy }， u ^T ={ u _x u _y }。

③本构方程。即材料的应力-应变关系。

（3-3）

式中 D ——弹性矩阵；

E ， μ ——材料的弹性模量和泊松比。

对于平面问题有 。

④边界条件。

（3-4）

式中 Γ _F ——面力和集中力边界；

Γ _u ——位移边界。

（2）变分法

变分法的基础是变分原理，而变分原理是求解连续介质问题的常用数学方法之一。

例如：一维稳态热传导的定解问题。

（3-5）

式中 ϕ ——未知温度场函数；

k ——热导率；

Q （ x ）——沿 x 方向分布的热载荷。

式（3-5）中的温度场 ϕ 可以借助傅里叶积分求得。

如果采用变分原理求解这类定解问题，则应首先建立定解问题的积分形式：

（3-6）

式中 u ——未知函数；

F ， E ——特定算子；

Ω ——连续求解域；

Γ —— Ω 的边界；

Π ——泛函数（即未知函数 u 的函数）。

此时，如果连续介质问题有解 u ，则解 u 必定使泛函 Π 对于微小变化δ u 取驻值（极值），即泛函的“一阶变分”等于零。

（3-7）

这就是所谓求解连续介质问题的变分原理。

可以证明：用微分方程加边界条件求解连续介质问题同用约束或非约束变分法求解连续介质问题等价：即一方面满足微分方程及其边界条件的函数将使泛函取得驻值（极值）；另一方面，从变分角度看，使泛函取得驻值（极值）的函数正是满足连续介质问题的微分方程及其边界条件的解。

在建立有限元方程时应用变分法的目的，就是将求微分方程的定解问题转变成求泛函的驻值问题，以方便试探函数的分片插值和分片积分。

（3）虚位移方程（位移变分方程）

虚位移方程是变分原理在求解弹性力学问题中的具体应用，它等价于几何微分方程和应力边界条件。所谓虚位移，是指在约束条件允许的范围内，弹性体内质点可能发生的任意微小位移。虚位移的产生与弹性体所受外力及其时间无关。

弹性体在外力的作用下发生变形，表明外力对弹性体做了功。若不考虑变形过程中的热损失、弹性体的动能和外界阻尼，则外力功将全部转换为储存于弹性体内的位能（应变能）。根据能量守恒定律，有：

（3-8）

式中 Ω ——弹性体内部；

Γ _f —— Ω 的面力边界（假设边界上无离散的集中力）；

， p ——面力和体积力；

δ u ——虚位移；

δ ε ——虚应变（即对应虚位移的任意可能应变）。

式（3-8）表明，外力（包括 Ω 内的体积力和边界 Γ _f 上的面力）使弹性体产生虚位移所做的功等于弹性体因虚应变而提高的能量。

3.2.2.2　有限元方程建立与计算

下以平面刚度问题为例，来讲述有限元方程的建立与求解过程。

（1）离散处理

假设分析对象（构件）的厚度尺寸非常小，可以近似将其处理成厚度为常数的平面问题。用三节点三角形单元离散该对象及其边界条件[图3-4（a）]，并从中任取一子域进行单元分析[图3-4（b）]。

图3-4中，单元节点 i 、 j 、 m 按逆时针排布。

（2）单元刚度分析

①单元位移与单元形函数　一个三节点三角形单元共有12个自由度（6个节点位移分量和6个节点转动分量），在线弹性范围内，6个节点转动分量可以忽略不计，由此可给出单元位移 a ^e 的向量表达式如下：

（3-9）

求解离散后的平面刚度问题存在两种情况：①位移分量是节点坐标的已知函数，此时，可直接利用单元节点位移分量求出单元应变分量（几何方程），再由单元应变分量求出单元应力分量（本构方程），最后综合起来便可得到整个平面刚度问题的解；②只有几个结点的位移分量已知，不能直接求出应变和应力分量。对于后者，为了利用结点位移表示单元应变和应力，必须构造一个位移模式（位移函数）。理论上，定义于某一闭区域内的任意函数总可以被一个多项式逼近，所以，位移函数常常取为多项式。现选用一次多项式作为三节点三角形单元的位移模式（位移插值函数），于是有：

（3-10）

式中 u ， v ——单元内任意一点的位移分量；

a ₁ ～ a ₆ ——待定系数。

将三角形的三个节点坐标代入式（3-10），得：

（3-11）

应用克莱姆法则求解式（3-11），并化简：

（3-12）

式中 u _k ， v _k ——节点位移， k = i ， j ， m （下同）；

N _k （ x ， y ）——三节点三角形单元的形函数（一个与单元类型和单元内部节点坐标有关的连续函数）；

A ——三角形单元的面积；

a _k ， b _k ， c _k ——与节点坐标有关的常数， a _k = x _j y _m - x _m y _j ， b _k = y _j - y _m ， c _k =- x _j + x _m 。

用矩阵表示式（3-12），得：

（3-13）

式中 I ——单位矩阵；

N ——形函数矩阵；

a ^e ——单元节点位移列阵。

式（3-12）和式（3-13）即为采用多项式插值构造的三节点三角形单元位移模式，它表明只要知道了节点位移，就能借助单元形函数求得单元内任意一点的位移。换句话说，节点位移通过形函数控制整个单元的位移分布。

单元形函数具有以下两条基本性质：

a.在任一节点上，形函数满足：

b.针对单元内任一点，有：

②单元应变矩阵与应力矩阵　将式（3-13）代入几何方程式（3-2），可得单元应变表达式：

（3-14）

式中 B ——单元应变矩阵，其分块子矩阵为：

（3-15）

再将单元应变表达式（3-14）代入本构方程（3-3），可得单元应力表达式：

（3-16）

式中 s ——单元应力矩阵。

其分块矩阵为：

（3-17）

式中 E ₀ ， μ ₀ ——材料常数；

E ， μ ——材料的弹性模量和泊松比。

对于平面应力问题， E ₀ = E ， μ ₀ = μ ；对于平面应变问题， ， 。

由于应变矩阵 B 中的每一个非零元素均是由节点坐标决定的常数，而节点坐标为定值，所以，应变矩阵 B 是常量矩阵。同理，应力矩阵 s 也是常量矩阵。这表明由式（3-14）和式（3-16）计算获得的单元内各点的应变值相同，应力值亦相同。可以证明，在弹性力学范围内，由节点位移引发的应变和应力主要通过相邻单元的边界节点进行传递。

③单元刚度矩阵与单元刚度方程　现利用3.2.2.1小节介绍的虚位移原理建立单元刚度方程。假设：在外力 q ^e （相邻单元或体积力、边界载荷）的作用下，单元结点 i ， j ， m 发生了虚位移δ a ^e ，虚位移引起虚应变δ ε ，由此得到单元的虚位移方程：

（3-18）

式中 t ——单元厚度。

因虚位移可以是任意的，由此得单元刚度方程（单元刚度模型）：

（3-19）

式中 k ^e ——单元刚度矩阵，表明单元 e 中各节点产生单位位移而引发（或需要）的节点力。

（3-20）

（3-21）

单元刚度矩阵 k ^e 具有以下特征。

a.对称性，即 k _rs = k _sr 。

b.奇异性，即单元刚度矩阵不存在其逆矩阵。

c.各行、各列元素之和恒为0（因在力的作用下，整个单元处于平衡状态）。

d.主对角元素恒为正，即 k _ii >0，表明节点位移与施加其上的节点力同向。

④单元等效节点力　作用在单元上的外力 q ^e 由等效节点力构成，即：

（3-22）

式中 ， ， ——体积力 p 、面力 和集中力 P _c 产生的等效节点力。

所谓等效节点力，是指由作用在单元上的体积力、面力和集中力按照能量等效原则（即原载荷与等效节点载荷在虚位移上做功相等原则）移置到节点上而形成的节点力。

图3-5是作用在单元上的外力举例，表3-1是对应于图3-5的等效节点力表达式。

（3）整体刚度分析

①建立总刚度矩阵与总刚度方程　将式（3-19）改写成如下形式：

（3-23）

式（3-23）表明，当单元内任一节点发生位移时，都将在该节点处产生节点力（实际上是由节点力引起位移），大小等于单元中各节点位移所引起的节点力之和。

因为在被离散对象的整体结构中，一个节点往往为若干比邻单元所共有，根据线性叠加原理，作用在该节点上的力应等于所有比邻单元的等效节点力之和，即：

（3-24）

式中 Q _i ——作用在节点 i 上的合力；

——比邻单元施加给 i 的等效节点力之和；

n _e ——比邻单元个数。

将式（3-23）代入式（3-24），取 r = i ，得到当节点 i 处于平衡状态时的合力表达式：

（3-25）

假设被离散对象（区域）有 n 个节点，于是该对象中所有节点的平衡方程为：

（3-26）

用矩阵表示，有：

（3-27）

式中 K ——总刚度矩阵；

Q ——总结点载荷列阵；

a ——总节点位移列阵。

式（3-27）为利用虚位移方程和节点力叠加原理建立起来的平面刚度问题的整体平衡方程（整体刚度模型），即用有限元格式表示的平面问题总刚度方程。

总刚度矩阵 K 具有以下性质。

a.对称性。总刚度矩阵由单元刚度矩阵叠加形成，所以，与单元刚度矩阵一样具有对称性，即 K = K ^T 。

b.稀疏性。由于任一节点仅与少数节点和单元比邻，相对于该节点，其无关节点在 K 中的对应元素为零，因此，存在于 K 中的大量零元素使总刚度矩阵具有稀疏性（稀疏矩阵）。

c.带状性。总刚度矩阵中的所有非零元素均集中分布在主对角线附近，从而形成所谓带状矩阵。

d.奇异性。对于无任何节点约束或约束不足的结构件，其总刚度矩阵奇异（ K |=0），物理上表现为在力的作用下，结构件整体作刚性运动，即方程（3-27）的解不唯一。

总刚度矩阵的上述特性，为其在计算机中的存储与处理，以及解的唯一性判别和计算方法的选取提供了良好的数学依据。

②设置位移边界条件　由于总刚度矩阵具有奇异性，因此，在求解总刚度方程（3-27）之前，必须设置一些节点位移约束，以防止结构件产生整体刚性运动。节点的位移约束一般施加在被离散对象的边界上，故称之为位移边界条件，见图3-4（a）。边界上节点的位移约束通常分为两大类。

a.零位移约束。设某一结点沿某一方向的位移为零（例如 a _m =0），则平面构件的总刚度方程变成：

（3-28）

观察式（3-28）可以发现：总刚度矩阵中与零位移节点 a _m 对应的主对角元素 k _nn 为1，其余相关元素（即下角标中含有 n 的 k 元素）均为0；载荷列阵中与零位移结点对应的元素 Q _n 为0。

b.非零已知位移约束。已知某一结点沿某约束方向位移为 ，其平面构件的总刚度方程转变成：

（3-29）

此时，式（3-29）表现为：总刚度矩阵中与已知结点位移 对应的主对角元素 k _mm 乘以一个足够大的正数 α （例如， α =10 ²⁰ ），相关列元素（即第二个下角标是 m 的 k 元素）均为零，并且载荷列阵中的 Q _m 被 所代替。

实际上，两类位移边界条件常常同时存在，所以可根据情况将式（3-28）和式（3-29）组合在一块进行化简。

③设置载荷边界条件　设置载荷边界条件实际上是为总载荷列阵中的各分量元素赋初值。结合位移边界条件设置可以简化赋值过程，即只针对未经式（3-28）和式（3-29）处理的载荷元素赋值。需要注意的是：除零位移和已知位移节点外，载荷中的等效体积力应平均加载到所有剩余节点上，等效面力应加载到面力边界节点上，而等效集中力则应根据情况加载到相应单元的节点上。这样处理后的总载荷列阵中，有的载荷分量可能为零（对应零位移节点），有的载荷分量可能是2～3种等效节点力的叠加（例如，某边界节点既受等效面力作用，又受等效体积力作用）。

（4）求解总刚度方程

代入约束和载荷边界条件的总刚度方程可写成：

（3-30）

式中 ——经约束边界条件和载荷边界条件处理后的总刚度矩阵与总节点载荷列阵。

式（3-30）是一个关于节点位移分量的线性方程组，利用迭代法、高斯法、波前法或带宽法等数值方法可以求出节点位移列阵 a 中的各分量，然后代入几何方程求解应变分量，再代入本构方程求解应力分量，于是便获得已知边界条件下平面刚度问题的全部近似解。

3.2.3　有限元解的收敛性与误差控制

基于单元形函数插值的有限元解是原应用问题的近似解。近似解是否收敛于真实解，近似解收敛速度有多块，近似解是否稳定，近似解误差有多大，这些都是决定有限元法能否成功用于具体工程问题的关键。

3.2.3.1　解的收敛性

有限元解的收敛性取决于所构建的试探函数（插值函数）逼近真实函数的程度，因此，试探函数的选择是关键。试探函数的选择必须遵循两条基本准则。

①完备性准则　如果出现在泛函（3-6）中的场函数最高导数是 m 阶，则有限元解收敛的条件之一是单元内场函数的试探函数至少是 m 次完全多项式。

满足完备性准则的单元被称为完备单元。

②协调性准则　如果出现在泛函（3-6）中的场函数最高导数是 m 阶，则试探函数在单元的交界面上必须具有 C _{m -1} 连续性，即在相邻单元的界面上，试探函数应具有直至 m -1阶连续导数。

满足协调性准则的单元被称为协调单元。

完备性是有限元解收敛的必要条件，而协调性是有限元解收敛的充分条件。显然，对于绝大多数工程问题而言，选择多项式作为单元插值函数一般都能满足完备性和协调性准则。采用既完备又协调的单元离散求解域，所获得的有限元解一定收敛，即当单元尺寸趋近于零时，有限元解趋近于真实解。但是，某些非协调单元在一定的条件下也能使有限元解收敛于真实解。

3.2.3.2　误差来源与控制

有限元解的误差主要来源于有限元建模误差和数值计算误差，其中，有限元建模误差包括求解对象（域）离散误差、边界条件误差和单元形状误差。图3-6是有限元解（或有限元分析）的误差组成。

（1）离散误差

离散误差有物理离散误差与几何离散误差之分。前者是插值函数（试探函数）与真实函数之间的差异，后者是单元组合体与原求解对象在几何形状上的差异。

物理离散误差量级可以定性地用以下公式估计：

（3-31）

式中 h ——单元特征尺寸；

p ——单元插值函数（一般为多项式）的最高阶次；

m ——单元插值函数的最高阶导数。

对于三节点三角形单元，其插值函数是线性的，即 p =1。由于求解位移不涉及插值函数的求导，所以 m =0，于是可推断出位移误差的量级为 O （ h ² ），位移解的收敛速度量级与位移误差相同，也为 O （ h ² ）。反之，节点位移与单元应变的几何方程中存在位移函数的一阶求导（ m =1），故应变误差的量级为 O （ h ），应变解收敛速度的量级也为 O （ h ）。同理，根据本构方程可推断出应力误差与应力解收敛速度的量级均为 O （ h ）。由此可见，在弹性力学刚度分析中，节点位移误差对其后应力应变解的精确度影响极大。

控制物理离散（即用插值函数代替真实函数）误差的方法主要有两种。

①减小单元特征尺寸　换句话说，就是增加网格划分密度。通过减小单元特征尺寸来提高有限元解精确度的方法被称为 h 法。

②提高插值函数的阶次　即采用较高阶的多项式插值单元离散求解域。通过提高插值函数的阶次来提高有限元解精确度的方法被称为 p 法。

此外，根据工程实际问题，混合使用高、低阶插值函数单元离散求解域的不同部分，也是控制物理离散误差的有效方法之一。

几何离散误差主要来自被离散对象的边界。如图3-7所示，中心带孔的平面圆被四边形单元均匀离散，其四边形的直边与内外圆周边之间存在间隙（即几何离散误差）。控制几何离散误差常用的方法有两种。

①网格局部加密（图3-8）。

②选择边和（或）面上带有节点的单元，因为这类单元的边和（或）面可以弯曲（图3-9）。

（2）边界条件误差

边界条件误差主要源于两个方面。

①物理边界的量化误差　物理边界量化误差实际上包括边界数据采集误差和边界数学模型误差，前者与某些边界条件的复杂性和边界数据采集的难度有关，后者与边界条件的理论研究和数学建模有关。例如：在钣金件冲压过程中，板料与模具之间的摩擦是动态变化的，而且依赖于具体的界面工况，很难准确地测定其摩擦数据，并利用基于理论分析与物理实验的界面摩擦模型加以描述。

②边界载荷等效移置误差　这类误差主要影响与边界载荷等效移置相关的局部区域特征，而对工程问题的整体求解影响不大。

控制边界条件误差的途径如下。针对第一种误差，尽量采用各种先进的技术方法与实验手段，准确测定工程问题求解所需的边界数据；同时深入开展理论与实验研究，不断完善描述界面条件的数学模型。针对第二种误差，细分重点关注的边界区域网格（例如：模拟焊接热应力时的热影响区），以减少或消除因载荷等效移置带来的求解精度损失。

（3）单元形状误差

单元形状误差由极度不规则的单元结构产生。例如：当图3-10所示三节点三角形单元的底高比（即 a/b ）非常大时，就会造成单元的严重畸变或退化，从而影响有限元求解精度，极端情况下还将导致求解失败或数值计算无法进行。

单元形状误差的影响一般仅限于畸变单元内部或相邻单元，因此，应有针对性地通过局部细分或单元编辑调整关键区域（例如，应力集中区）的网格。

（4）计算误差

计算误差可能来源于计算方法、程序设计、运算次数、误差累积，以及解题性质与解题规模等多个方面。计算误差又分为舍入误差和截断误差，前者主要与计算机内部用于数据存储和字长处理、求解线性方程组和数值积分所需要的运算次数、数值计算采用的计算方法，以及计算方法在程序实现中的误差控制等因素有关，而后者主要与数值计算采用的计算方法、解题性质和解题规模有关。

对于计算误差的控制，可根据实际计算情况，具体问题具体解决。例如：如果计算误差是由解题规模过大引起，则应采用适当措施降低解题规模，以减少运算次数和由此带来的累积误差；如果属于因计算方法选择不当而导致计算效率降低和计算误差增大，则应重新选择高效高精度的计算方法。

3.2.4　非线性问题的有限元法

前面结合有限元方程建立与应用所举的例子为线性问题，在材料成形领域，线弹性体系的刚度分析可用于各种成形模具、工装夹具和焊接结构等的开发与设计，稳态传热分析可用于工件保温、退火、自然冷却、人工时效和砂型铸模传热等过程模拟。然而，材料成形过程中大量遇到的却是一些非线性工程问题，例如：固态金属冲压和锻造中的冷、热塑性变形，液态金属铸造中的流动、凝固与传热，固体材料熔化焊中的传热、物理冶金与焊缝凝固，热塑性塑料注射成型中的黏性流动与冷却固化，模具零件淬火热处理中的传热与相变等。

有限元法在主要材料成形领域（如冲压、锻压、铸造、焊接、注射等）中应用所涉及的某些专业知识，将根据需要分散到后续章节讨论与补充，本小节仅简要介绍与材料成形相关的非线性问题基本概念和求解非线性方程组的基本方法。

3.2.4.1　材料成形领域中的非线性问题

材料成形领域遇到的非线性问题主要体现在以下三个方面。

①材料非线性　材料非线性是指材料变形时的应力与应变关系（本构关系）非线性。图3-11分别为弹塑性、刚塑性、刚黏塑性和黏弹塑性材料在拉伸或压缩变形时的典型应力应变曲线。

比较图3-11中的各曲线特征可知，当前三类材料的变形进入塑性区后，均存在某种程度的应变硬化现象，即随着材料塑性变形量的增加，维持其变形所需的应力（流动应力）也增加，并且两者之间的关系是非线性的；刚（黏）塑性材料与弹塑性材料的区别在于前者的弹性变形相对于其塑性变形可以忽略不计，即假设材料屈服前为刚性；刚塑性材料与刚黏塑性材料的区别在于后者的应力应变关系还与应变速率有关；最后一类黏弹塑性材料的应力应变关系属于高度非线性，其典型代表为（中等结晶的）热塑性塑料。

实际上，不同的成形加工方法会使同样材料呈现不同的非线性性质。例如：金属材料的冲压成形，其应力应变遵循弹塑性本构关系；锻压成型，遵循刚（黏）塑性本构关系；铸造成型，遵循牛顿黏性流体本构关系。

②几何非线性　几何非线性是指由物体内质点的大位移和大转动引起的非线性，力学上表现为研究对象的几何方程不满足线性关系。例如图3-12所示的材料弯曲，塑变区的材料质点不仅存在大位移，而且还存在某种程度的转动。

③边界非线性　边界非线性是指边界条件呈现非线性变化。例如：模锻时，毛坯与模膛表面的接触和摩擦（图3-13），即使不考虑软硬质点的黏着问题，其接触点3也不会沿模膛表面呈线性滑动。

3.2.4.2　求解非线性方程组的基本方法

无论是材料非线性问题或是几何非线性问题，经有限元法处理后，最终都将被归结为求解非线性方程组。离散化的非线性方程组一般可表示为：

（3-32）

或

式中 a ——未知场函数的近似解；

K （ a ）——非线性方程组的系数矩阵；

Q ——外载荷列阵。

由式（3-32）可知，非线性方程组的系数矩阵是变量矩阵。在工程上，常常借助增量法将载荷或时间离散成若干个增量步，针对每一步载荷或时间增量，“线性化”方程组（3-32），即将非线性问题转化成一系列线性问题进行求解。具体做法概括起来如下：

①离散载荷或时间为 m 个增量步；

②设全局载荷初值或时间初值，利用迭代法计算第一增量步（ i =1）内的“线性”方程组；

③当第一增量步内的迭代计算误差小于规定值后，即将最后一次的迭代结果作为第一增量步的解；

④判断 m 个增量步是否全部计算完毕，即不等式 i > m 是否成立；

⑤如果 i ≤ m ，则 i = i +1，并以当前增量步的迭代解作为初值，进行下一增量步的迭代计算；

⑥循环第④、⑤步工作，直到 i > m 。

3.2.4.3　非线性有限元解的稳定性

当利用增量-迭代混合法求解方程组（3-32）时，增量步长的选取对有限元解的稳定性影响极大。所谓有限元解的稳定性，是指当载荷步或时间步的长度（步长）取不同值时，方程组（3-32）的收敛误差是否趋于恒定或波动最小。如果增量步长取任意值，误差都不会无限增长，则称有限元解为无条件稳定；如果增量步长只有在满足一定条件时，误差才不会无限增长，则称有限元解为条件稳定。图3-14表示计算某瞬态传热过程，当时间步长Δ t 分别取1.5和2.6时所对应的有限元解收敛误差变化轨迹。其中，横坐标表示迭代次数，纵坐标表示迭代计算的收敛误差。增量步长的选取受多种因素影响，具体方法请参阅后续章节的相关内容。

3.2.5　有限元求解应用问题的一般步骤

对于不同物理性质和数学模型的问题，有限元求解法的基本步骤是相同的，只是具体公式推导和运算求解不同。有限元求解问题的基本步骤如下。

第一步：问题及求解域定义。根据实际问题近似确定求解域的物理性质和几何区域。

第二步：求解域离散化。将求解域近似为具有不同有限大小和形状且彼此相连的有限个单元组成的离散域，习惯上称为有限元网络划分。显然，单元越小（网格越细），则离散域的近似程度越好，计算结果也越精确，但计算量将增大，因此求解域的离散化是有限元法的核心技术之一。

第三步：确定状态变量及控制方法。一个具体的物理问题通常可以用一组包含问题状态变量边界条件的微分方程式表示，为适合有限元求解，通常将微分方程化为等价的泛函形式。

第四步：单元推导。对单元构造一个适合的近似解，即推导有限单元的列式，其中包括选择合理的单元坐标系，建立单元试函数，以某种方法给出单元各状态变量的离散关系，从而形成单元矩阵（结构力学中称刚度阵或柔度阵）。

为保证问题求解的收敛性，单元推导有许多原则要遵循。对工程应用而言，重要的是应注意每一种单元的解题性能与约束。例如，单元形状应以规则为好，畸形时不仅精度低，而且有缺秩的危险，将导致无法求解。

第五步：总装求解。将单元总装形成离散域的总矩阵方程（联合方程组），反映对近似求解域的离散域的要求，即单元函数的连续性要满足一定的连续条件。总装是在相邻单元结点进行，状态变量及其导数（可能的话）连续性建立在结点处。

第六步：联立方程组求解和结果解释。有限元法最终导致联立方程组。联立方程组的求解可用直接法、迭代法和随机法。求解结果是单元结点处状态变量的近似值。对于计算结果的质量，将通过与设计准则提供的允许值比较来评价并确定是否需要重复计算。

简而言之，有限元分析可分成三个阶段：前置处理、计算求解和后置处理。前置处理是建立有限元模型，完成单元网格划分；后置处理则是采集处理分析结果，使用户能简便提取信息，了解计算结果。 FqYturIgzaobo5IwwnqlfikOg6PsOi6xJMCgtPSQl3jv4YDL/kIdBb9XuvWOf38i