有限元与有限差分法基础(超详细版本)

1、2022-4-281/23有限元法基础及有限差分法有限差分法基础有限元法有限差分法有限差分法2022-4-282/162有限元法基础n有限元发展过程n有限元应用n有限元发展方向2022-4-283/162有限元法的基本思想n基本思想 1）将连续的求解系统离散为一组由节点相互联在一起的单元组合体 2）在每个单元内假设近似函数来分片表示系统的求解场函数 2022-4-284/162有限元法的基本思想2022-4-285/162有限元法的基本思想2022-4-286/162有限元法的基本思想2022-4-287/162有限元法的基本思想2022-4-288/162有限元法的基本思想n离散为单元网格的

2、件仍然要保证是一个连续体，单元与单元之间没有裂缝、不能重叠，所有单元通过单元节点相互关联着n变形体无论产生多大的塑性变形，单元与单元之间依然不会产生裂缝、交叉和重叠，关联单元的节点也不能脱开2022-4-289/162有限元法的基本思想n不合格单元 单元裂缝单元裂缝单元重叠单元重叠2022-4-2810/162有限元法的基本思想n变形前后单元之间都是连续的变形前的网格变形前的网格变形后的网格变形后的网格2022-4-2811/162有限元法的基本思想n基本思想n通过在单元内假设不同的插值函数，建立不同的单元模型，适应各种各样的变形模式和受力模式XFXF2022-4-2812/162有限元法的基

3、本思想n有限元法分类 1）位移法：基于最小势能原理或虚功原理 2）力法： 基于最小余能原理 3）混合法：基于修正余能原理 2022-4-2813/162有限元法的基本思想n基本过程离散化过程 约束处理过程 单元平衡方程组装过程 应变、应力回代过程 方程组求解过程 2022-4-2814/162离散化过程n最小势能原理 P V A G 弹性体 弹性体的势能peipWW 为弹性体变形后所具有的内能 iW为弹性体所受的外力功 eWdVVT21VAdVdAGuPuTT2022-4-2815/162离散化过程 为弹性体的应变 为弹性体的应力 u为弹性体的可容位移 弹性体处于平衡状态时，其势能应为最小 P

4、0TTTVAVdVdAdVGuPu02022-4-2816/162离散化过程n单元插值关系 eNuu n单元几何关系 n单元本构关系 Lu DeN为单元形函数矩阵 L为单元几何微分算子 为单元弹性矩阵 eD0TTTvavPdvdadvGuPu0)()()(TTTTTTvevaeeeedvdadvGNuPNuuBDBu0TTTvvaeedvdadvGNPNuBDBeu单元节点自由度向量2022-4-2817/162离散化过程0)()()(TTTTTTvevaeeeedvdadvGNuPNuuBDBu0TTTvvaeedvdadvGNPNuBDB0TTTvavPdvdadvGuPuB 称为应变矩阵

5、 LNB fku e单元平衡方程或单元刚度方程 k 称为单元刚度矩阵 vedvBDBkT f 称为单元载荷向量 vadvdaGNPNfTT2022-4-2818/162单元刚度矩阵的特性 n对称性 n奇异性 n主元恒正且对角占优 离散化过程2022-4-2819/162线弹性问题几何方程三维问题 Luwvuxzyzxyzyxzuxwywzvxvyuzwyvxuzxyzxyzzyyxx000000000三维问题Lu 2022-4-2820/162Luwvuxzyzxyzyxzuxwywzvxvyuzwyvxuzxyzxyzzyyxx000000000线弹性问题几何方程二维问题 Luvuxyyxy

6、vxuyvxuxyyyxx00二维问题平面应力和平面应变状态 2022-4-2821/162线弹性问题几何方程二维问题 二维问题轴对称状态 Luwurzzrrrwzuzwrururzzzrr0010Luwvuxzyzxyzyxzuxwywzvxvyuzwyvxuzxyzxyzzyyxx0000000002022-4-2822/162线弹性问题几何方程一维问题 一维问题 LuuxxuxxLuvuxyyxyvxuyvxuxyyyxx002022-4-2823/162线弹性问题本构方程三维问题 三维问题De221000000221000000221000000100010001)21)(1 (EeD

7、E为弹性模量；为泊松比 2022-4-2824/162线弹性问题本构方程平面应力 二维问题平面应力状态 00yzxz000yzxzzzzxyzxyzzyyxx000 xyyyxxxyyyxxzxyzxyzzyyxxxyzzyyxx2022-4-2825/162线弹性问题本构方程平面应力 zxyzxyzzyyxxezxyzxyzzyyxxD2100010112EeDxyyyxxxyyyxxE2100010112平面应力状态 00000 xyzzyyxxexyyyxxD2022-4-2826/162线弹性问题本构方程平面应变 二维问题平面应变状态 000yzxzzz00yzxzzxyzxyzzyy

8、xx00 xyzzyyxxxyzzyyxxzxyzxyzzyyxxxyyyxx2022-4-2827/16200000 xyyyxxexyzzyyxxD线弹性问题本构方程平面应变 zxyzxyzzyyxxezxyzxyzzyyxxDxyyyxxxyyyxxE221000101)21)(1 (平面应变状态 221000101)21)(1 (EeD2022-4-2828/162线弹性问题本构方程轴对称 二维问题轴对称状态 00yzxz00yzxzzxyzxyzzyyxx00 xyzzyyxxxyzzyyxxzxyzxyzzyyxxxyzzyyxx2022-4-2829/162线弹性问题本构方程轴对

9、称 二维问题轴对称状态 00zr00zrzrzrzzrrzrzzrr00zrzzrrzrzrzzrrzrzzrr2022-4-2830/162线弹性问题本构方程轴对称 轴对称状态 zxyzxyzzyyxxezxyzxyzzyyxxDzrzrzzrrezrzrzzrrDzrzzrrezrzzrr0000D221000010101)21)(1 (EeDzrzzrrzrzzrrE221000010101)21)(1 (2022-4-2831/162线弹性问题本构方程一维问题 一维问题xxxxEEeD2022-4-2832/162常用单元模型 n单元模型插值关系一一对应n单元类型一维单元、二维单元、三

10、维单元等参单元、超参单元、次参单元2022-4-2833/162常用单元模型n一维单元 2节点线单元 1 2 1 2 3 1 2 3节点线单元梁单元2022-4-2834/162常用单元模型n二维单元3节点三角形线性单元 1 2 3 6节点三角形二次单元 1 2 3 5 6 4 2022-4-2835/162常用单元模型n二维单元10节点三角形三次单元 10 1 2 3 6 9 4 5 7 8 4节点四边形双线性单元 1 2 3 4 2022-4-2836/162常用单元模型n二维单元8节点四边形二次单元12节点四边形三次单元 1 2 3 4 8 7 6 5 8 7 1 2 3 4 11 9

11、5 6 10 12 2022-4-2837/162常用单元模型n三维单元4节点四面体线性单元10节点四面体二次单元 1 2 3 4 1 10 9 8 4 7 2 3 6 5 2022-4-2838/162常用单元模型n三维单元8节点六面体线性单元20节点六面体二次单元 8 6 1 2 3 4 5 7 8 16 17 10 3 20 19 18 15 14 6 13 12 11 9 2 4 1 5 7 2022-4-2839/162常用单元模型n准三维空间单元桁架单元一维2节点线单元+单元局部随体坐标系 为什么要建立单元局部随体坐标系 ？1.简化分析问题的复杂程度。2.在局部坐标系中，空间桁架的

12、每根杆每变成了一维2节点线单元2022-4-2840/162常用单元模型n准三维空间单元框架单元三维梁单元+一维2节点线单元+单元局部随体坐标系 两端都是刚性联结 可以要承受拉压、弯曲、扭转3种变形模式 框架单元的特点2022-4-2841/162常用单元模型n准三维空间单元板单元薄板单元中厚板单元弯曲和横向剪切2种变形模式抵抗板的变形如果板很薄，忽略横向剪切抗力，认为抵抗载荷的主要因素是弯矩2022-4-2842/162常用单元模型n准三维空间单元壳单元 抵抗拉压变形的二维单元+板单元+单元局部随体坐标系。适合于薄壳单元和中厚壳单元从几何上分为薄壳单元和中厚壳单元 组合单元2022-4-28

13、43/162常用单元模型n准三维空间单元 壳理论单元 由空间壳理论严格构造的壳单元。适合于薄壳单元和中厚壳单元 退化单元 由三维实体单元退化成的壳单元。只适合于中厚壳单元 2022-4-2844/162单元模型构造 n有限元法的基本思想 通过单元分片近似，在每个单元内假设近似函数来分片表示系统的场函数 n选择近似函数简单、实用的原则在有限元法中，近似函数称为插值函数 2022-4-2845/162单元模型构造n插值函数 一般都采用多项式函数，主要原因是: 采用多项式插值函数比较容易推导单元平衡方程，特别是易于进行微分和积分运算。随着多项式函数阶次的增加，可以提高有限元法的计算精度。从理论上说，

14、无限提高多项式的阶数，可以求得系统的精确解。2022-4-2846/162单元模型构造方法 n整体坐标系法n局部坐标系法 nLagrange插值方法nHermite插值方法2022-4-2847/162单元模型构造方法n2节点线单元12 oxu1u2x1x2ux1. 假设插值多项式xaaxu10)(2. 利用节点值求 a0 和 a1 21021101xaauxaau12121xxuua1212210 xxxuxua2022-4-2848/162单元模型构造方法3. 代入a0 和 a1，得插值多项式 u(x)xxxuuxxxuxuxu1212121221)(4. 按u1 和 u2合并同类项，设

15、l = x2- x1212121122112)(uuNNuulxxlxxulxxulxxxu2022-4-2849/162单元模型构造方法n关键 如何构造插值多项式 u ?二维问题三维问题，如何构造插值多项式?2022-4-2850/162n收敛性条件 在单元内，场函数必须是连续的； 完备性：插值多项式的阶次必须由低到高依次增加，不能出现跳跃现象； 协调性：各单元边界必须连续，单元边界不能出现开裂现象。 插值多项式收敛性条件 n收敛：当单元逐渐缩小时，如果插值多项式满足收敛性条件，则数值解将收敛于精确解 2022-4-2851/162插值多项式收敛性条件n协调单元 满足插值多项式收敛性条件和的

16、单元 n完备单元 满足插值多项式收敛性条件的单元ncr 阶连续性 插值多项式的第r阶导数是连续的 2022-4-2852/162插值多项式收敛性条件n非协调单元与部分协调单元 对于一般固体力学问题来说，协调性要求单元在变形时，相邻单元之间不应引起开裂、重叠或其它不连续现象。例如，梁、板、壳等单元，在单元边界不但要求位移是连续的，而且其一阶导数也必须是连续的。板、壳单元位移函数沿单元边界的法向导数（转角）的连续性一般比较难实现，因此出现了许多不完全满足协调性要求的“非协调单元”或“部分协调单元”，有时它们的精度也很好。 2022-4-2853/162插值多项式选择条件 n插值多项式应该尽可能满足

17、其收敛性条件（收敛性）n由插值多项式所确定的场函数变化应该与局部坐标系的选择无关（各向同性） n假设的插值多项式系数的数量应该等于单元的节点数（解的唯一性） 选择条件2022-4-2854/162插值多项式选择条件n深入分析由收敛性条件可知，插值多项式中必须含有常数项（刚体位移项），高阶项的次数必须依次增加，不允许有跳跃332210)(xxxxu3310)(xxxu2022-4-2855/162插值多项式选择条件由选择条件可知，插值多项式函数在所有自由度方向上要满足各向同性性，这样就不会随局部坐标系变化而改变了 n深入分析xyyxy, xu3210)(23210)(xyxy, xu2022-4

18、-2856/162插值多项式选择条件n深入分析选择条件是为了能由单元节点值唯一确定插值多项式 4节点四边形的插值多项式应该是 xyyxy, xu3210)(插值多项式系数i (i = 0,1,2,3) 也是4个 2022-4-2857/162单元模型构造整体坐标系法n基本思想 针对弹性体有限元网格建立一个统一的坐标系，每个单元的插值多项式都在这个坐标系上建立 y x o 弹性体 2022-4-2858/162单元模型构造整体坐标系法n2节点线单元12 oxu1u2x1x2ux1. 假设插值多项式xaaxu10)(2. 利用节点值求 a0 和 a1 21021101xaauxaau12121xx

19、uua1212210 xxxuxua2022-4-2859/162单元模型构造整体坐标系法3. 代入a0 和 a1，得插值多项式 u(x)xxxuuxxxuxuxu1212121221)(4. 按u1 和 u2合并同类项，设 l = x2- x1212121122112)(uuNNuulxxlxxulxxulxxxu2022-4-2860/162单元模型构造整体坐标系法lxxN2111 122122( )euu xN uN uNNuNulxxN12N1 和 N2 称为单元的形函数；N 称为单元的形函数矩阵；ue 称为单元节点位移向量。 n2节点线的单元形函数2022-4-2861/162单元模

20、型构造整体坐标系法n二维3节点三角形单元 (x1 , y1) (u1 , v1) 1 y x o v 3 2 (x2 , y2) (x3 , y3) u (u3 , v3) (u2 , v2) (u , v) (x , y) 建立整体坐标系oxy 2022-4-2862/162单元模型构造整体坐标系法1. 假设插值多项式yxyxu210),(2. 首先，利用节点值求 0 、 1 和 2 n二维3节点三角形单元 323103222102121101yxuyxuyxuyxyxv210),(2022-4-2863/162单元模型构造整体坐标系法)(21)(21)(213322112332211133

21、22110ucucucAubububAuauauaA33221111121yxyxyxA A为单元面积2022-4-2864/162单元模型构造整体坐标系法23132113311xxcyybyxyxa31213221122xxcyybyxyxa12321332233xxcyybyxyxa3. 将 0 、 1 和 2 代入插值多项式，按u1、u2、u3合并同类项332211),(uNuNuNyxu2022-4-2865/162单元模型构造整体坐标系法)(21)(21)(21333322221111ycxbaANycxbaANycxbaAN4. 同理可得332211),(vNvNvNyxv2022

22、-4-2866/162单元模型构造整体坐标系法5. 单元插值多项式为332211332211),(),(vNvNvNyxvuNuNuNyxu)(21)(21)(21333322221111ycxbaANycxbaANycxbaAN2022-4-2867/162单元模型构造整体坐标系法332211321321000000vuvuvuNNNNNNvu6. 单元插值多项式写成矩阵形式（常用）2022-4-2868/162单元模型构造整体坐标系法321321321321000000vvvuuuNNNNNNvu7. 单元插值多项式的另一种矩阵形式（不常用）2022-4-2869/162单元模型构造整体坐

23、标系法n4节点四面体单元 2 1 (u1 , v1 , w1) x y z 4 3 (x1 , y1 , z1) (x3 , y3 , z3) (x4 , y4 , z4) (u3 , v3 , w3) (u2 , v2 , w2) (x2 , y2 , z2) (u4 , v4 , w4) (x , y , z) (u , v , w) 2022-4-2870/162单元模型构造整体坐标系法zyxzyxwzyxzyxvzyxzyxu421042104210),(),(),(1. 假设插值多项式443322114433221144332211),(),(),(wNwNwNwNzyxwvNvNv

24、NvNzyxvuNuNuNuNzyxu2. 插值多项式为2022-4-2871/162单元模型构造整体坐标系法)(61zdycxbaVNiiiii（i=1,2,3,4） 4443332221zyxzyxzyxa 4433221111zyzyzyb4433221111zxzxzxc1114433221yxyxyxd 444333222111111161zyxzyxzyxzyxV 循环轮换脚标1、2、3、4，相应可以得到a2，b2 ， c2 ， d2 、 a3 ， b3 ， c3 ， d3 、 a4 ， b4 ， c4 ， d4 2022-4-2872/162单元模型构造整体坐标系法4443332

25、22111432143214321000000000000000000000000wvuwvuwvuwvuNNNNNNNNNNNNwvu3. 单元插值多项式写成矩阵形式（常用）2022-4-2873/162单元模型构造整体坐标系法4. 单元插值多项式另一种矩阵形式（不常用）432143214321432143214321000000000000000000000000wwwwvvvvuuuuNNNNNNNNNNNNwvu2022-4-2874/162单元模型构造整体坐标系法n从理论上讲，整体坐标系法可以求任意单元的形函数，但计算过程太复杂n只能求一维2节点线单元、二维3节点三角形单元和三维4节

26、点四面体单元3种简单单元的形函数n复杂的或二次以上的单元必须采用局部坐标系法求n位移场 u 是形函数 Ni 的线性组合，因此形函数Ni同样具有插值多项式的特性2022-4-2875/162单元刚度矩阵2节点线单元n一维2节点线单元n单元插值关系 eNuu n单元几何关系 n单元本构关系 Lu DeN=N1 N2 De=E21uueuxLlxxNlxxN/ )(/ )(12212022-4-2876/162单元刚度矩阵2节点线单元n单元刚度矩阵vedvBDBkTvedvBDBkT x yzedxdydzBDBTxedxABDBTA为单元截面积；l为单元长度n矩阵BLNB LNB /21xNxN/

27、1/1llBDBeAlT1111lAE/1/1/1/1llEllAl2022-4-2877/162单元刚度矩阵三角形单元n二维3角形单元n单元插值关系 eNuu T332211vuvuvueu321321000000NNNNNNN)(21)(21)(21333322221111ycxbaANycxbaANycxbaAN2022-4-2878/162单元刚度矩阵三角形单元n单元几何关系 Lu Luvuxyyxxyyyxx002022-4-2879/162单元刚度矩阵三角形单元n单元本构关系 DexyyyxxxyyyxxE2100010112平面应力问题2022-4-2880/162单元刚度矩阵三

28、角形单元n矩阵BLNB LNB xNyNxNyNxNyNyNyNyNxNxNxN33221132132100000033221132132100000021bcbcbccccbbbA2022-4-2881/162单元刚度矩阵三角形单元n单元刚度矩阵vedvBDBkTvedvBDBkT x yzedxdydzBDBTx yedxdyhBDBTh为单元厚度BDBehATk为对称的6*6常数矩阵A为单元面积2022-4-2882/162作业n求4节点四面体单元的单元刚度矩阵 2 1 (u1 , v1 , w1) x y z 4 3 (x1 , y1 , z1) (x3 , y3 , z3) (x4

29、, y4 , z4) (u3 , v3 , w3) (u2 , v2 , w2) (x2 , y2 , z2) (u4 , v4 , w4) (x , y , z) (u , v , w) 2022-4-2883/162单元模型构造整体坐标系法n单元形函数的特性正规性：单元形函数之和等于1。 正交性：形函数在本节点的值等于1，在其它节点的值等于0。 lxxN21lxxN12例如: 2节点线单元形函数2022-4-2884/162单元模型等参单元n等参单元 单元内任意一点的位移u与单元节点位移ue之间的关系为 eNuu 一般单元坐标的插值关系也采用与位移插值关系相同的变换关系即单元内任意一点的坐

30、标x与单元节点坐标xe之间的关系为 eNxx 2022-4-2885/162单元模型等参单元n等参单元凡是几何形状和位移场采用同阶同参数插值关系来描述的单元，称为等参单元 前面介绍的所有单元都属于等参单元 在描述单元的几何形状和位移场时，并不一定非采用同阶插值关系 2022-4-2886/162单元模型等参单元332211321321000000vuvuvuNNNNNNvu332211321321000000yxyxyxNNNNNNyxn等参单元3节点三角形等参单元 2022-4-2887/162单元模型等参单元n超参单元如果几何形状插值函数的阶数高于位移场插值函数的阶数，称为超参单元 n次参

31、单元如果几何形状插值函数的阶数低于位移场插值函数的阶数，称为次参单元 2022-4-2888/162单元平衡方程组装过程 n为什么要组装 ? 2 F 1 3 消除内力n组装的原则是什么 ? 单元自由度与结构自由度对应2022-4-2889/162单元平衡方程组装过程 2 F 1 3 U3U4U2U1U5U6654321UUUUUUU结构自由度向量U2022-4-2890/162单元平衡方程组装过程 2 1 U3U4U2U11u1u2u3u4 3 U611 U2U12u1u2U5u3u443214321UUUUuuuueu65214321UUUUuuuueu22022-4-2891/162单元平

32、衡方程组装过程14131211432114414314214113413313213112412312212111411311211143214321ffffuuuukkkkkkkkkkkkkkkk2 1 U3U4U2U11u1u2u3u422022-4-2892/162单元平衡方程组装过程000000000000000000000065432165432114131211654321144143142141134133132131124123122121114113112111ffffUUUUUUkkkkkkkkkkkkkkkk组装单元2022-4-2893/162单元平衡方程组装过程 3

33、U611 U2U12u1u2U5u3u424232221432124424324224123423323223122422322222121421321221165216521ffffuuuukkkkkkkkkkkkkkkk2022-4-2894/162单元平衡方程组装过程242314132212211165432124424324224123423323223114414314214113413313213122422312412322212222112121421311411321211221111100000000654321654321ffffffffUUUUUUkkkkkkkkkkk

34、kkkkkkkkkkkkkkkkkkkkk再组装单元FKU 总体刚度方程 K 称为总体刚度矩阵 U 称为位移向量 F 称为载荷向量 2022-4-2895/162总体刚度矩阵K的特性 n对称性 n奇异性 n稀疏性 n非零元素带状分布 2022-4-2896/162约束处理过程 n为什么要约束处理 ?总体平衡方程组是奇异的消除无限制的刚体运动 使总体平衡方程组存在唯一一组解2022-4-2897/162约束处理过程边界条件n 边界条件分类 力(载荷)边界条件位移边界条件 集中载荷力 表面分布力 自重力热交换引起的温度载荷 固定位移约束 强制位移约束 关联位移约束 002.U 005.U 011.

35、U 514.U CkUU782022-4-2898/162约束处理过程模型简化xyU2022-4-2899/162约束处理过程模型简化yxUxyU2022-4-28100/162约束处理过程约束方程123456789101112yxU0 . 0321UUU0 . 012963VVVVUUUU1211102022-4-28101/162约束处理过程约束处理方法n位移约束处理方法 赋0赋1法 乘大数法 2022-4-28102/162约束处理过程赋0赋1法n强制位移约束条件处理 U4=CFKU 6543216543216665646362615655545352514645444342413635

36、34333231262524232221161514131211FFFFFFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK2022-4-28103/162约束处理过程赋0赋1法n强制位移约束条件处理 U4=C654321654321666564636261565554535251464544434241363534333231262524232221161514131211FFFFFFUUUUUUKKCKKKKKKCKKKKKKCKKKKKKCKKKKKKCKKKKKKCKKKKCKFCKFCKFCKFCKFCKFUUUUUUKKKKKKKKKKKKKKK

37、KKKKKKKKKKKKKKK6465454443432421416543216665636261565553525146454342413635333231262523222116151312110000002022-4-28104/162约束处理过程赋0赋1法n有6个方程，5个未知数，如果约束方程可以消除有限元平衡方程组的奇异性，则取任意5个方程联立求解，都会得到方程组的唯一一组解。 n系数矩阵由原来的对称的变成了非对称的，这对于大规模有限元方程组求解是十分不利的，采用相同的求解方法，在求解时间和矩阵存贮容量方面都增加了一倍。 2022-4-28105/162约束处理过程赋0赋1法CKFC

38、KFCKFCKFCKFCKFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK646545444343242141654321666563626156555352514645434241363533323126252322211615131211000000CKFCKFCKFCKFCKFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKK64654534324214165432166656362615655535251363533323126252322211615131211000000000000为了保证系数矩阵的对称性，去掉方程组第4行2022-4-281

39、06/162CKFCKFCKFCKFCKFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKK64654534324214165432166656362615655535251363533323126252322211615131211000000000000约束处理过程赋0赋1法引入强制位移约束方程 U4=C，使方程组求解时直接将自由度U4求出CKFCKFCCKFCKFCKFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKK64654534324214165432166656362615655535251363533323126252322211615131211000

40、010000002022-4-28107/162约束处理过程赋0赋1法n固定位移约束条件处理 U4=06532165432166656362615655535251363533323126252322211615131211000001000000FFFFFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKK2022-4-28108/162约束处理过程赋0赋1法n基本原理 利用初等变换对求解方程组进行相同的行列变换，既保证方程组解不会改变，又可以保持方程组系数矩阵的对称性。 在进行初等变换时，只要保证对方程组系数矩阵做相同的行列变换，就可以保持方程组系数矩阵的对称性。 2022-4-

41、28109/162约束处理过程乘大数法n乘大数法n基本原理 利用矩阵的初等变换不改变方程组解的思想。 2022-4-28110/162654321654321666564636261565554535251464544434241363534333231262524232221161514131211FFFFFFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK约束处理过程乘大数法n强制位移边界条件 CU 42022-4-28111/162约束处理过程乘大数法CU 4强制约束方程 ACAU 4A是一个大数，是系数矩阵中对角线元素K44的1010倍量级以上 为什

42、么要乘以大数A ？放大位移约束方程的优势2022-4-28112/162654321654321666564636261565554535251464544434241363534333231262524232221161514131211FFACFFFFUUUUUUKKKKKKKKKKKKKKAKKKKKKKKKKKKKKKKKKKKKK约束处理过程乘大数法n强制位移边界条件 CU 42022-4-28113/162约束处理过程乘大数法n固定位移边界条件 C = 0 约束后的方程组简化为 6543216543216665646362615655545352514645444342413635

43、34333231262524232221161514131211FFFFFFUUUUUUKKKKKKKKKKKKKKAKKKKKKKKKKKKKKKKKKKKKK654321654321666564636261565554535251464544434241363534333231262524232221161514131211FFACFFFFUUUUUUKKKKKKKKKKKKKKAKKKKKKKKKKKKKKKKKKKKKK2022-4-28114/162约束处理过程两种方法比较n赋0赋1法在约束处理过程中是严格精确的，而乘大数法是一种近似约束处理方法，它的精度取决于所乘大数A值 两种方法

44、都可以消除有限元平衡方程的奇异性，得到符合实际边界条件的唯一一组解。但两种方法还是有很大的区别 2022-4-28115/162约束处理过程两种方法比较n采用乘大数法约束处理后的有限元平衡方程在求解时可能造成解的失真，大数A值越大可能解的偏差会越大，而赋0赋1法就不会出现类似的问题，它在约束过程和求解过程都是精确的n乘大数法相对于赋0赋1法在约束处理过程上简单一些 2022-4-28116/162约束处理过程两种方法比较n赋0赋1法实际上是将关联位移约束方程代入到有限元平衡方程中的，是代入法。而乘大数是将占绝对优势的关联位移约束方程合并到有限元平衡方程中的，是罚方法，计算误差来自于合并过程，计

45、算精度取决于关联位移约束方程的优势大小n商业软件中，位移边界条件的约束处理都采用赋0赋1法，乘大数很少被采用主要原因是它是一种近似方法，而且大数的大小也不好确定，有时还会造成求解失败 2022-4-28117/162约束处理过程弹簧单元假设柔性弹簧 kOXYU4f f = kU4k2022-4-28118/162弹簧约束方程 f = kU4654321654321666564636261565554535251464544434241363534333231262524232221161514131211FFFFFFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK

46、KKK654321654321666564636261565554535251464544434241363534333231262524232221161514131211FFfFFFFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK6544321654321666564636261565554535251464544434241363534333231262524232221161514131211FFkUFFFFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK65432165432166656463626156555

47、4535251464544434241363534333231262524232221161514131211FFFFFFUUUUUUKKKKKKKKKKKKKKkKKKKKKKKKKKKKKKKKKKKKK约束处理过程弹簧单元2022-4-28119/162方程组求解过程特点n方程组求解有限元计算过程中很重要的一部分，在有限元法的发展过程中，有限元方程的求解效率一直是其应用的最大瓶颈之一 n有限元方程组的特点： 有限元方程组的系数矩阵具有对称、稀疏、带状分布以及正定、主元占优。有效地利用这些特点，以减少系数矩阵的存贮量，提高方程组求解效率 2022-4-28120/162方程组求解过程分类比

48、较n线性方程组的解法主要分两大类: 直接解法：以高斯消去法基础，以等带宽或变带宽方式存贮系数矩阵内元素，对于求解规模比较大的问题，要存贮的元素非常巨大。 迭代解法：只需要存贮系数矩阵中非零元素，存贮量很小，一般是变带宽存贮量的20%或更少，有些算法的求解效率也非常高，适合求解大规模线性方程组。但是这种解法对接近病态的方程组很难保证收敛性。 2022-4-28121/162方程组求解过程带宽定义n有限元方程组系数矩阵是稀疏的、非零元素呈带状分布，带宽就是它的宽度，带宽的大小是由系统有限元网格的节点号排序决定的，具体求法是 带宽=(单元最大节点号之差+1)*节点自由度数 n带宽是网格节点标注方法直

49、接决定的，不同标注方法带宽可能相关很大 2022-4-28122/162方程组求解过程带宽n带宽是网格节点标注方法直接决定的，不同标注方法带宽可能相关很大 1 2 3 4 5 6 9 10 7 8 11 12 14 13 16 15 20 19 18 17 21 22 23 24 28 26 27 25 1 2 3 4 5 6 9 10 7 8 11 12 14 13 16 15 20 19 18 17 21 22 23 24 28 26 27 25 1 2 3 4 5 6 9 10 7 8 11 12 14 13 16 15 20 19 18 17 21 22 23 24 28 26 27

50、25 2022-4-28123/162方程组求解过程带宽n所示四边形网格的三种节点号标注方法，每个节点是2个自由度n结构的带宽分别是12，18，56，相差很大，其中12和56之间相差近5倍，这就意味着系数矩阵的存贮量也是相差5倍，因此，对于大规模复杂系统的节点号优化是十分必要的 2022-4-28124/162方程组求解过程系数矩阵存贮 n系数矩阵存贮 如果节点号排序优化的比较好，系数矩阵的存贮量就会减少很多。根据系数矩阵的对称性，一般都是按半带宽存贮。n系数矩阵存贮的方法 二维等带宽存贮 一维变带宽存贮2022-4-28125/162方程组求解过程二维等带宽存贮 n二维等带宽存贮 D D 0

51、 0 K K n 2022-4-28126/162方程组求解过程二维等带宽存贮n二维等带宽存贮消除了最大带宽以外的全部零元素，节省了系数矩阵元素的存贮量。但是由于取最大带宽为存贮范围，因此不能排除在带宽内的大量零元素。当系数矩阵的各行带宽变化不大时，适合采用二维等带宽存贮，方程组求解过程中系数矩阵元素的寻址也比较方便，求解效率较高。n当出现局部带宽特别大的情况时，采用二维等带宽存贮时，将由于局部带宽过大而使整体系数矩阵的存贮大大增加。2022-4-28127/162方程组求解过程一维变带宽存贮 n一维变带宽存贮 一维变带宽存贮方法就是把变化的带宽内的元素按一定的顺序存贮在一个一维数组中。由于它

52、不按最大带宽存贮，因此比二维等带宽存贮更节省内存。按照解法可分为按行一维变带宽存贮和按列一维变带宽存贮。 2022-4-28128/162n按行一维变带宽存贮 对称 20181296431M方程组求解过程一维变带宽存贮 辅助的寻址数组M 2022-4-28129/162n一维变带宽存贮是最节省内存的一种方法，但是由于要借助于寻址数组寻找系数矩阵元素的位置，相对二维等带宽存贮方法来说要复杂一些，而且在程序实现时也要复杂得多，方程组求解过程中也要消耗一些数组寻址时间。因此，在选用存贮方法时要权衡二者的利弊，统盘考虑。一般当带宽变化不大，计算机内存允许时，采用二维等带宽存贮方法是比较合适的。 方程组

53、求解过程一维变带宽存贮 2022-4-28130/162方程组求解过程求解方法n方程组求解方法 高斯消去法 三角分解法 雅可比（Jacobi）迭代法 高斯-赛德尔（Gauss-Seidel）迭代法 2022-4-28131/162应变、应力回代过程 n单元应变和应力回代求解 通过求解有限元平衡方程得到有限元节点位移后，就可以进行系统的刚度校核。如果所分析问题要进行强度校核，就要回代求解单元的应变和应力。n由插值关系和几何关系可得单元应变，再通过本构关系得到单元应力eeBuLNuDe2022-4-28132/162有限差分法有限差分法2022-4-28133/162 从弹性力学的基本方程建立以来

54、，这些方程在各种问题的边界条件下如何求解，一直是很多数学工作者和力学工作者研究的内容。即弹性力学的经典解法存在一定的局限性，当弹性体的边界条件和受载情况复杂一点，往往无法求得偏微分方程的边值问题的解析解，许多工程重要问题，不能够得出函数式的解答。 因此，弹性力学问题的各种数值解法便具有重要的实际意义。2022-4-28134/162 工程中常用的数值解法有有限单元法和差分法。 有限单元法 是以有限个单元的集合体来代替连续体，属于物理上的近似。 差分法 是把弹性力学的基本方程和边界条件（一般均为微分方程）近似地改用差分方程（代数方程）来表示，把求解微分方程的问题改换成为求解代数方程的问题，属于数

55、学上的近似。2022-4-28135/162 第一节 差分方程 第二节 应力函数的差分解 第三节 深梁应力函数的差分解2022-4-28136/162第一节 差分方程 差分法是沿用已久的一种数值解法。随着计算机的普及和相应的软件发展，此法成为解弹性力学问题的一种有效的方法。2022-4-28137/162 我们在弹性体上，用相隔等间距h而平行于坐标轴的两组平行线织成正方形网格，x=y=h，如图。.)(! 31)(! 21)(3003320022000 xxxfxxxfxxxfff 设f=f(x,y)为弹性体内的某一个连续函数。该函数在平行于x轴的一根网线上，如在-上,它只随x坐标的改变而变化。

56、在邻近结点处，函数f可展为泰勒级数如下：2022-4-28138/162 我们将只考虑离开结点充分近的那些结点，即（x-x0）充分小。于是可不计（x-x0）的三次及更高次幂的各项，则上式简写为：在结点，x=x0-h, 在结点1, x=x0+h，代入(b) 得：)()(! 21)(20022000bxxxfxxxfff2230200( )2fhfffhcxx2210200( )2fhfffhdxx联立(c),(d),解得差分公式： )11(2310hffxf)21 (22031022hfffxf同理，在网线-上可得到差分公式2402240220(13)22(14)2fffyhffffyh 差分公

57、式(-)及(-)是以相隔2h的两结点处的函数值来表示中间结点处的一阶导数值，可称为中点导数公式。 以相邻三结点处的函数值来表示一个端点处的一阶导数值，可称为端点导数公式。 应当指出：中点导数公式与端点导数公式相比，精度较高。因为前者反映了结点两边的函数变化，而后者却只反映了结点一边的函数变化。因此，我们总是尽可能应用前者，而只有在无法应用前者时才不得不应用后者。 以上(-)(-)是基本差分公式,从而可导出其它的差分公式如下：2268570401391144040123456782240402410124401()()(1 5)4164()()(1 6)142()()(1 7)164()()(1

58、 8)fffffx yhffffffxhffffffffffxyhffffffyh 2022-4-28142/162第二节 应力函数的差分解 当不计体力时，我们已把弹性力学平面问题归结为在给定边界条件下求解双调和方程的问题。用差分法解平面问题，就应先将双调和方程变换为差分方程，而后求解之。2240220021302200257682001()21()2(21)1()()4xyxyyhxhx yh 一旦求得弹性体全部节点的值后，就可按应力分量差分公式（对节点0）算得弹性体各节点的应力。 可见，用差分法解平面问题，共有两大任务：一、建立差分方程 将（1-68）代入双调和方程012345678910

59、1112208() 2() () 0(2 2) 对于弹性体边界以内的每一结点，都可以建立这样一个差分方程。024422444yyxx整理即得2022-4-28145/162 二、联立求解这些线性代数方程，就能求得各内结点处的值。 为了求得边界上各结点处的值，须要应用应力边界条件，即： 一般建立和求解差分方程，在数学上不会遇到很大困难。但是，当对于边界内一行的（距边界为h的）结点，建立的差分方程还将涉及边界上各结点处的值，并包含边界外一行的虚结点处的值。yyxyxxyxpllpll2121代入上式，即得： yxxyxyyx22222,l1=cos(N,x)=cos=dy/ds, l2=cos(N

60、,y)=sin=-dx/ds,yxpxsxyxsypyxsxysy222222dddddddd于是，式（a）可改写为：由右图可见，2212222122axyllpyx yllpx yx 关于边界上任一点处BAyBABAxBAspxspyddyxpxsxyxsypyxsxysy222222dddddddddbdddxypsypsx（ ）由此得： xy的值，可将（b）式从A点到B点对s积分得到： 将此式亦从A点到B点沿s进行积分，就得到边界上任一点B处的值。为此利用分部积分法，得： 由高等数学可知，syysxxsdd.dd.dd BABABABABAsysyyysxsxxx,dddddddcdBx