线性有限元方法

例如， 例如，利用虚功原理
δ {q}e
∫ ([ B ][ N ]) [ D ] ([ B ][ N ]) dV {q} = δ {q} ∫ [ N ] {b} dV + δ {q} ∫ [ N ] {t}
T T V e e T e T T e T V e e Sσ
e
dSe
[ K ]e {q}e = {F}e
上述都是仅以位移为基本变量的原理（弱形式）； 仅以位移为基本变量的原理 注：上述都是仅以位移为基本变量的原理（弱形式）； 变分形式因问题会有所不同，如梁和板问题。 变分形式因问题会有所不同，如梁和板问题。
2011-10-19 8
控制方程（4/6） 控制方程（4/6）
积分（ 积分（弱）形式的控制方程： 形式的控制方程：
1）自然（面力）边界条件： ）自然（面力）边界条件：
σ ij n j − tei = 0
∂ ∂x [ B] = ∂ ∂y ∂ ∂z ∂ ∂y ∂ ∂z ∂ ∂y ∂ ∂x
σ ij = λδ ij ε kk + 2 µε ij
ε ij =
2011-10-19
1 +ν ν σ ij − σ kk δ ij E E
µ = E 2 (1 +ν ) λ = Eν (1 +ν )(1 − 2ν )
拉梅（ 拉梅（Lame）常数与工程参数 ） 杨氏模量、泊松比）的关系。 （杨氏模量、泊松比）的关系。
再定义应力偏量和应变偏量分别为 再定义应力偏量和应变偏量分别为 应力偏量和应变偏量
eij = ε ij − ε mδ ij Sij = σ ij − σ mδ ij
K = E G = E 3 (1 − 2ν ) 2 (1 +ν )
本构关系又可写成
2011-10-19 11
有限元方程的建立（1/3） 有限元方程的建立（1/3）
本构关系、几何关系及控制方程，无论是采用解析方法， 本构关系、几何关系及控制方程，无论是采用解析方法， 解析方法 还是采用半解析方法 半解析方法， 还是采用半解析方法，对求解一个固体力学问题都必须用 有限元方法与这些方法的最主要区别在于 与这些方法的最主要区别在于， 到。有限元方法与这些方法的最主要区别在于，它一开始 就着眼于通过数值途径获得结果，因而， 就着眼于通过数值途径获得结果，因而，在离散的有限单 元上进行插值逼近是有限元方法的主要特征。 元上进行插值逼近是有限元方法的主要特征。
2011-10-19 9
控制方程（5/6） 控制方程（5/6）
积分（ 积分（弱）形式的控制方程： 形式的控制方程：
4）基于混合变量的广义变分原理－两变量变分原理 ）基于混合变量的广义变分原理－ 广义变分原理 简称HR） 海－耐（简称 ）变分原理的变分为
δΠ ( u , σ ) =
∫ δ {σ } ([ B ]{u} − [ D ] {σ })dV
在塑性力学中广泛使用。 在塑性力学中广泛使用。
2011-10-19
4
几何关系
小变形问题的几何关系为： 小变形问题的几何关系为：
应变－位移关系的张量形式： 应变－位移关系的张量形式：
1 ε ij = ( u i , j + u j ,i ) 2
有限元插值
{u} = [ N ]{q}
微分算子矩阵 （作用于形函数 N） ）
有限元方法是基于单元的， 有限元方法是基于单元的，这是与其它类型数值方法的 本质区别，它的优点和缺点都是因之而产生。 本质区别，它的优点和缺点都是因之而产生。 有限元方程通常都是在单元层面上建立的， 有限元方程通常都是在单元层面上建立的，但必须经 过组装，在整个结构层面上施加约束、加以求解。 过组装，在整个结构层面上施加约束、加以求解。
单元刚度矩阵
2011-10-19
[ K ] = ([ B ][ N ])T [ D ] ([ B ][ N ]) dV e ∫V e T T { F }e = ∫V [ N ] {b}e dVe + ∫S [ N ] {t}e dS e σ
单元等效载荷
13
有限元方程的建立（3/3） 有限元方程的建立（3/3）
on Sσ
[T ] {σ } − {te } = 0
T
[T]与[B]相似 [T]与[B]相似，但用法向 相似， 余弦替换其中的偏微分。 余弦替换其中的偏微分。
注：自然边界条件将包含在积分形式的控制方程中。 自然边界条件将包含在积分形式的控制方程中。
∂ ∂x ∂ ∂z
2）本质（位移）边界条件： ）本质（位移）边界条件：
3
线性本构关系（3/3） 线性本构关系（3/3）
若定义平均正应变和平均正应力（或称八面体正应力） 若定义平均正应变和平均正应力（或称八面体正应力）分别为 平均正应变和平均正应力
ε m = (ε x + ε y + ε z ) 3 σ m = (σ x + σ y + σ z ) 3
2011-10-19 1
线性本构关系（1/3） 线性本构关系（1/3）
线弹性本构关系（ 定律） 线弹性本构关系（Hooke定律） 定律
Hooke定律的张量形式： 定律的张量形式： 定律的张量形式
σ ij = Dijkl ε kl
σ11 σ12 σ13 σ σ σ 21 22 23 σ31 σ32 σ33
∂ ∂x [ B] = ∂ ∂y ∂ ∂z ∂ ∂z ∂ ∂y ∂ ∂x
应变－位移关系的矩阵形式： 应变－位移关系的矩阵形式：
{ε } = [ B ]{u }
∂ ∂y
T 2）矩阵形式： [ B ] {σ } + {b} = 0 ）矩阵形式：
几何 关系
2011-10-19
{ε } = [ B ]{u }
与几何关系中相同的微 分算子矩阵（6×3） 分算子矩阵（ × ）
注：只有在笛卡尔坐标系 中才有这样的关系
6
控制方程（2/6） 控制方程（2/6）
边界条件： 边界条件：
2011-10-19
2
线性本构关系（2/3） 线性本构关系（2/3）
[D]的元素个数：原始矩阵为6×6＝36个； 的元素个数：原始矩阵为 × ＝ 个 的元素个数 根据对称性变成21个 根据对称性变成 个； 对于正交各向异性（复合）材料为 个 对于正交各向异性（复合）材料为9个； 对于横观各向同性（复合）材料为 个 对于横观各向同性（复合）材料为5个； 对于各向同性材料为2个 对于各向同性材料为 个。 对各向同性材料， 的显式很容易导得， 对各向同性材料，Dijkl与[D]的显式很容易导得，例如 的显式很容易导得
3）基于混合变量的广义变分原理－三变量变分原理 ）基于混合变量的广义变分原理 广义变分原理－ 简称HW）变分原理的广义泛函为 胡－鹫（简称 ）
1 T T Π ( ε , σ , u ) = ∫ {σ } ([ B ]{u} − {ε } ) + {ε } [ D ]{ε } dV V 2
第一章 线性有限元方法
目的：以固体的弹性静力学为例， 目的：以固体的弹性静力学为例，介绍线性有限元方程 建立过程中的关键环节。 建立过程中的关键环节。 地位： 地位：以后章节的预备知识 内容： 内容：
线性本构关系 几何关系 控制方程 有限元方程的建立 有限单元形状函数的构造 有限单元形状函数插值精度 四边形单元向三角形单元的退化
ε m = σ m ( 3K ) eij = Sij ( 2G )
体积、剪切模量与工程参数（ 体积、剪切模量与工程参数（杨氏 模量、泊松比）的关系。 模量、泊松比）的关系。 注：这些本构方程都可化成轴 对称、平面应力、平面应 对称、平面应力、 变等问题的具体形式。 变等问题的具体形式。
Sσ T
− ∫ δ {u }
V
T
([ B ] {σ } + {b})dV
T
({t} − {t })dS
e
5）余能原理 ）
δΠ c
2011-10-19
δ U c − δ Wc = 0
δ Wc = uiδ bi dV + uiδ tei dS ∫V ∫Sσ δ U c = ∫V ε ijδσ ij dV
三类变分宗量
− ∫ {b} {u}dV − ∫ {te } {u}dS
T T V S
可以证明， 变分原理等价于材料本构关系 可以证明，HW变分原理等价于材料本构关系、几何关 变分原理等价于材料本构关系、 微分控制方程及面力边界条件。 系、微分控制方程及面力边界条件。 注：原我校教授匡震邦老师在这方面做出了重要 贡献， 年获得国家自然科学二等奖。 贡献，于1983年获得国家自然科学二等奖。 年获得国家自然科学二等奖
{u}e = [ N ]e {q}e
单元上的位移场
2011-10-19
单元上的 结点位移
单元的形状插值函数
12
有限元方程的建立（2/3） 有限元方程的建立（2/3）
单元位移 单元应变 单元应力
{u}e = [ N ]e {q}e
{ε }e = [ B ][ N ]e {q}e
{σ }e = [ D ][ B][ N ]e {q}e
注：可化成轴对称，平面应力、应变等问题的具体形式。 可化成轴对称，平面应力、应变等问题的具体形式。
2011-10-19
∂ ∂x ∂ ∂z
5
控制方程（1/6） 控制方程（1/6）
微分（ 微分（强）形式的控制方程：（应用力学常用） 形式的控制方程： 应用力学常用）
分布体积力 1）张量形式： 1）张量形式： σ ij , j + bi = 0
T −1 V
两类变分宗量
）－4） 注：1）－ ）都含有 ）－ 位移做为基本变 可以看出， 变分原理等价于几何关系、 变分原理等价于几何关系 可以看出，HR变分原理等价于几何关系、 可统称为势 量，可统称为势 微分控制方程及面力边界条件。 微分控制方程及面力边界条件。 能类变分原理。 能类变分原理。
+ ∫ δ {u }
ε11 ε12 ε13 ε ε ε 21 22 23 ε31 ε32 ε33
Hooke定律的矩阵形式： 定律的矩阵形式： 定律的矩阵形式
{σ } = [ D ]{ε }
T {ε} = {ε11 ε 22 ε33
2ε12 2ε 23 2ε31}
T
T T {σ} = {σ11 σ22 σ33 σ12 σ23 σ31}
特点： 特点：以（应）力为 基本变量
10
控制方程（6/6） 控制方程（6/6）
积分（ 积分（弱）形式的控制方程： 形式的控制方程：
5）Galerkin加权残值法－数学上更加通用和一般 ） 加权残值法－ 加权残值法 一般步骤： 一般步骤： 对微分控制方程、 对微分控制方程、边界条件等的残差在积 分意义下进行加权平均 加权平均； 分意义下进行加权平均； 识别（确定）权函数； 识别（确定）权函数； 进行求解。 进行求解。 派生的方法包括：最小二乘法、配点法、子域法、 派生的方法包括：最小二乘法、配点法、子域法、迦辽金法和 有限体积法等。 有限体积法等。 适用的问题包括：热传导、流体力学、空气动力学、电磁场、 适用的问题包括：热传导、流体力学、空气动力学、电磁场、 声场等几乎所有用二阶偏微分方程描述的问题。 声场等几乎所有用二阶偏微分方程描述的问题。
ui − ui 0 = 0
on
Su
注：位移边界条件就是通常所说的约束条件，必须予以强制满足。 位移边界条件就是通常所说的约束条件，必须予以强制满足。
2011-1wenku.baidu.com-19 7
控制方程（3/6） 控制方程（3/6）
积分（ 积分（弱）形式的控制方程：（有限元法常用） 形式的控制方程： 有限元法常用）
1）虚功（位移）原理： ）虚功（位移）原理：
∫
2）势能原理： ）势能原理：
V
δ {ε } {σ } dV = ∫ δ {u} {b} dV + ∫ δ {u} {te } dS
T T T V Sσ
δΠ δ U − δ W = 0
表面分布力
δ U = σ ijδε ij dV ∫V δ W = ∫V biδ u i dV + ∫Sσ t eiδ ui dS