2.1.2 有限元方程的建立与应用_材料成形过程数值模拟（第二版）-QQ阅读男生都市网

书名：材料成形过程数值模拟（第二版）
作者名：傅建肖兵
本章字数：6982字
更新时间：2020-08-28 06:50:42

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

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

2.1.2.1　预备知识

（1）弹性力学基本方程

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

Lσ+b=0　　在（Ω域内）　　（2-1）

式中　L——微分算子，对于平面问题；

　σ——应力，对于平面问题σ^T={σ_x　σ_y　τ_xy}；

　b——体积力（一般为重力）向量，对于平面问题b^T={b_x　b_y}。

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

ε=Lu　　（在Ω域内）　　（2-2）

式中　ε——应变，对于平面问题ε^T={ε_x　ε_y　γ_xy}；

　u——质点位移矢量，对于平面问题u^T={u_x　u_y}。

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

σ=Dε　　（在Ω域内）　　（2-3）

式中　D——弹性矩阵，对于平面问题；

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

④边界条件

Γ=Γ_F+Γ_u　　（在Ω的边界上）　　（2-4）

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

　Γ_u——位移边界。

（2）变分法

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

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

　　（2-5）

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

　 k——热导率；

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

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

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

　　（2-6）

式中　　u——未知函数；

　F、E——特定算子；

Ω——连续求解域；

Γ——Ω的边界；

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

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

δΠ=0　　（2-7）

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

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

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

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

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

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

　　（2-8）

式中　δε——虚应变（即对应虚位移的任意可能应变），δε^T={δε_xδε_yδε_zδγ_yzδγ_zxδγ_xy}；

　δu——虚位移，δu^T={δu　δv　δw}；

　Ω——弹性体内部；

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

，b——面力和体积力。

式（2-8）表明，外力（包括Ω内的体积力和边界上的面力）使弹性体内质点产生虚位移所做的功等于弹性体内部虚应变产生的能量。

2.1.2.2　平面刚度问题

（1）离散处理

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

图2-4　离散后的分析对象（a）与单元结构（b）

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

（2）单元刚度分析

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

　　（2-9）

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

　　（2-10）

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

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

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

　　（2-11）

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

　　（2-12）

式中　　　　u_k、v_k——节点位移，k=i，j，m（下同）；

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

——三角形单元的面积；

a_k，b_k，c_k——与节点坐标有关的常数，a_k=x_jy_m-x_my_j，b_k=y_j-y_m，c_k=-x_j+x_m。

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

　　（2-13）

式中　I——单位矩阵；

　N——形函数矩阵；

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

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

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

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

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

N_i（x，y）+N_j（x，y）+N_m（x，y）=1

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

　　（2-14）

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

　　（2-15）

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

σ=Dε=DBa^e=sa^e　　（2-16）

式中　s=DB=D［B_i　B_j　B_m］=［s_i　s_j　s_m］——单元应力矩阵，其分块矩阵为：

　　（2-17）

式中　E₀、μ₀——材料常数，对于平面应力问题E₀=E，μ₀=μ，对于平面应变问题。

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

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

③单元刚度矩阵与单元刚度方程　现利用2.1.2.1小节介绍的虚位移原理建立单元刚度方程。假设：在来自相邻单元“外力”q^e的作用下，单元节点i，j，m发生了虚位移δa^e，虚位移引起虚应变δε，由此得到单元的虚位移方程：

　　（2-18）

式中　t——单元厚度。

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

q^e=k^ea^e　　（2-19）

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

　　（2-20）

其中

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

a.对称性，即k_rs=k_sr；

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

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

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

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

　　（2-21）

式中　，分别为体积力b、面力和集中力p_c产生的等效节点力。

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

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

图2-5　体积力、面力和集中力举例

表2-1　对应于图2-5的等效节点力

（3）整体刚度分析

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

　　（2-22）

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

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

　　（2-23）

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

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

　 n_e——比邻单元个数。

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

　　（2-24）

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

　　（2-25）

用矩阵表示，有

Ka=Q　　（2-26）

式中　K——总刚度矩阵，；

　Q——总节点载荷列矩阵，；

　a——总节点位移列矩阵。

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

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

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

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

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

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

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

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

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

　　（2-27）

观察式（2-27）可以发现：总刚度矩阵中与零位移节点a_m对应的主对角元素k_mm为1，其余相关元素（即_足标中含有m的k元素）均为0；载荷列矩阵中与零位移节点对应的元素Q_m为0。

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

　　（2-28）

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

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

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

（4）求解总刚度方程

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

　　（2-29）

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

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

（5）平面刚度问题应用举例

①离散处理　假设：将如图2-6所示平面构件离散成8个三节点三角形单元；已知面载荷P，其等效结点力F=P/3均布在4、6、8三个节点上且垂直于x-z坐标面；各结点在z轴方向上的位移忽略不计（z方向尺寸很小，可固定设置为t），并忽略各节点的转动；构件受力时，节点1和9的x位移分量为0，1、9、2、10的y位移分量也为0。

图2-6　平面刚度问题举例

②单元刚度分析　例如：针对图2-6中的单元1建立节点1、4、2的平衡方程（单元刚度方程）

式中　k₁₁～k₆₆——2×2阶子矩阵，其子矩阵中的各元素与材料的弹性模量E、泊松系数μ、节点坐标（x_i，y_i），以及单元厚度t有关，见式（2-20）。

③整体刚度分析　将所有单元分析结果依次组合，得到图2-6平面构件的总刚度方程（有限元方程）。

{F}=［K］{δ}

式中　{F}={F₁F₂…F₁₀}^T={F₁_xF₁_yF₂_xF₂_y…F₁₀_xF₁₀_y}^T；

{δ}={δ₁δ₂…δ₁₀}^T={δ₁_xδ₁_yδ₂_xδ₂_y…δ₁₀_xδ₁₀_y}^T；

④设置边界条件　将位移约束δ₁_x=δ₉_x=δ₁_y=δ₉_y=δ₂_y=δ₁₀_y=0和等效节点力F₄_y=F₆_y=F₈_y=P/3分别代入总刚度方程的位移列矩阵和载荷列矩阵，整理并化简。

⑤计算求解　求解总刚度方程中剩余节点的位移分量，然后再分别利用几何方程和本构方程求解单元应变与单元应力。

⑥分析求解结果　可将计算获得的单元最大等效应力值与构件的许用强度进行比较，将节点最大位移值与构件的许用扰度（刚度）进行比较，以确定构件在面载荷P的作用下是否失效。取最大等效应力值和最大位移值作为失效判断的基本依据是：单元应力和节点位移在离散区域内不会发生突变，这与实际工程问题相吻合。

2.1.2.3　平面变温问题

变化的温度场会在弹性体内部诱导热应力，这种热应力主要来自于温度场变化所引起的弹性体不均匀热胀冷缩。平面变温问题的本构方程可表示为

σ=D（ε-ε₀）　在Ω域内　　（2-30）

式中　ε₀——温度变化引起的温度应变，ε₀=α（ϕ-ϕ₀）{1　1　0}^T；

　α——材料的线胀系数（假设各向同性）；

　ϕ——温度场；

　ϕ₀——初始温度场；

　ε——力载荷引起的应变；

　σ——应力；

　D——弹性矩阵。

同无温度场变化的本构方程相比，式（2-30）用（ε-ε₀）替代了式（2-3）中的ε。

将单元应变方程ε=Ba^e代入式（2-30），得平面变温问题的单元本构方程

σ^e=DBa^e-Dε₀　　（2-31）

若不考虑外力作用，则变温条件下单元位能的泛函表达式为

　　　　（2-32）

式中　——单元刚度矩阵；

　——等效节点热载荷。

总位能的泛函表达式

　　（2-33）

式中　——总刚度矩阵；

　——总变温载荷列矩阵。

根据能量泛函取驻值条件（最小位能原理）δΠ_p=0，对式（2-33）求一阶变分并化简，得

Ka=Q　　（2-34）

式（2-34）即为求解平面变温引起节点位移的有限元刚度方程。需要注意的是：求解变温问题的平面应力时，应将式（2-34）的计算结果代回变温问题的本构方程式（2-30）。利用泛函变分取驻值法建立限元方程是变分原理在弹性力学中的又一应用。同样，方程（2-34）中的Q也可包括表面、体积、集中等载荷。

2.1.2.4　平面稳态热传导问题

（1）热传导基本方程与边界条件

　在Ω域内　　（2-35）

　　　　　　　　　　　　　　　　　　　在Γ₁边界上　　（2-36）

　在Γ₂边界上　　（2-37）

　在Γ₃边界上　　（2-38）

式中　　　　ρ——密度；

　　c——比热；

　　t——时间；

k_x、k_y、k_z——沿x、y、z方向的传热系数；

　　Q——固体材料内部的热源密度，Q=Q（x，y，z，t）；

n_x、n_y、n_z——边界外法线的方向余弦；

　　——Γ₁边界上给定温度，；

　　q——Γ₂边界上给定热流密度（当q=0时，绝热边界），q=q（Γ，t）；

h（ϕ_a-ϕ）——Γ₃边界上给定界面传热；

　　h——界面传热系数；

　　ϕ——待求温度场（即与空间和时间相关的温度函数）；

　　ϕ_a——环境温度（自然对流条件下）或边界层绝热温度（强制对流条件下），ϕ_a=ϕ_a（Γ，t）；

　　Γ——域Ω的边界，Γ=Γ₁+Γ₂+Γ₃。

热传导基本方程（2-35）表示：当图2-7所示微元体的传热达到平衡时，微元体升温（或降温）所需（或减少）热量（方程左边第1项）应与传入（或传出）微元体的热量（方程左边第2～4项）和微元体内部释放（或吸收）的热量（方程左边第5项）相等。其中，热源密度Q与具体工程问题有关，可以是固态相变产生的相变潜热，也可以是压力加工的能量转换热，或热处理感应加热产生的欧姆热等。式（2-36）～式（2-38）分别称为第一类（强制）边界条件、第二类（自然）边界条件和第三类（自然）边界条件。