1、2023-2-121/23有限元法基础及有限差分法有限差分法基础有限元法有限差分法有限差分法2023-2-122/162有限元法基础n有限元发展过程n有限元应用n有限元发展方向2023-2-123/162有限元法的基本思想n基本思想 1）将连续的求解系统离散为一组由节点相互联在一起的单元组合体 2）在每个单元内假设近似函数来分片表示系统的求解场函数 2023-2-124/162有限元法的基本思想2023-2-125/162有限元法的基本思想2023-2-126/162有限元法的基本思想2023-2-127/162有限元法的基本思想2023-2-128/162有限元法的基本思想n离散为单元网格的
2、冲压件仍然要保证是一个连续体，单元与单元之间没有裂缝、不能重叠，所有单元通过单元节点相互关联着n板料无论产生多大的塑性变形，单元与单元之间依然不会产生裂缝、交叉和重叠，关联单元的节点也不能脱开2023-2-129/162有限元法的基本思想n不合格单元 单元裂缝单元裂缝单元重叠单元重叠2023-2-1210/162有限元法的基本思想n变形前后单元之间都是连续的变形前的网格变形前的网格变形后的网格变形后的网格2023-2-1211/162有限元法的基本思想n基本思想n通过在单元内假设不同的插值函数，建立不同的单元模型，适应各种各样的变形模式和受力模式XFXF2023-2-1212/162有限元法的
3、基本思想n有限元法分类 1）位移法：基于最小势能原理或虚功原理 2）力法：基于最小余能原理 3）杂交法：基于修正余能原理 4）混合法：基于Reissner变分原理 2023-2-1213/162有限元法的基本思想n位移法基本过程 1）离散化过程 3）约束处理过程 2）单元平衡方程组装过程 5）应变、应力回代过程 4）方程组求解过程 2023-2-1214/162离散化过程n最小势能原理 P V A G 弹性体 弹性体的势能peipWW 为弹性体变形后所具有的内能 iW为弹性体所受的外力功 eWdVVT21VAdVdAGuPuTT2023-2-1215/162离散化过程 为弹性体的应变 为弹性体
4、的应力 u为弹性体的可容位移 弹性体处于平衡状态时，其势能应为最小 P0TTTVAVdVdAdVGuPu02023-2-1216/162离散化过程n单元插值关系 eNuu n单元几何关系 n单元本构关系 Lu DeN为单元形函数矩阵 L为单元几何微分算子 为单元弹性矩阵 eD0TTTvavPdvdadvGuPu0)()()(TTTTTTvevaeeeedvdadvGNuPNuuBDBu0TTTvvaeedvdadvGNPNuBDBeu单元节点自由度向量2023-2-1217/162离散化过程0)()()(TTTTTTvevaeeeedvdadvGNuPNuuBDBu0TTTvvaeedvdad
5、vGNPNuBDB0TTTvavPdvdadvGuPuB 称为应变矩阵 LNB fku e单元平衡方程或单元刚度方程 k 称为单元刚度矩阵 vedvBDBkT f 称为单元载荷向量 vadvdaGNPNfTT2023-2-1218/162单元刚度矩阵的特性 n对称性 n奇异性 n主元恒正且对角占优 离散化过程2023-2-1219/162线弹性问题几何方程三维问题 Luwvuxzyzxyzyxzuxwywzvxvyuzwyvxuzxyzxyzzyyxx000000000三维问题Lu 2023-2-1220/162Luwvuxzyzxyzyxzuxwywzvxvyuzwyvxuzxyzxyzzy
6、yxx000000000线弹性问题几何方程二维问题 Luvuxyyxyvxuyvxuxyyyxx00二维问题平面应力和平面应变状态 2023-2-1221/162线弹性问题几何方程二维问题 二维问题轴对称状态 Luwurzzrrrwzuzwrururzzzrr0010Luwvuxzyzxyzyxzuxwywzvxvyuzwyvxuzxyzxyzzyyxx0000000002023-2-1222/162线弹性问题几何方程一维问题 一维问题 LuuxxuxxLuvuxyyxyvxuyvxuxyyyxx002023-2-1223/162线弹性问题本构方程三维问题 三维问题De221000000221
7、000000221000000100010001)21)(1(EeDE为弹性模量；为泊松比 2023-2-1224/162线弹性问题本构方程平面应力 二维问题平面应力状态 00yzxz000yzxzzzzxyzxyzzyyxx000 xyyyxxxyyyxxzxyzxyzzyyxxxyzzyyxx2023-2-1225/162线弹性问题本构方程平面应力 zxyzxyzzyyxxezxyzxyzzyyxxD2100010112EeDxyyyxxxyyyxxE2100010112平面应力状态 00000 xyzzyyxxexyyyxxD2023-2-1226/162线弹性问题本构方程平面应变 二维
8、问题平面应变状态 000yzxzzz00yzxzzxyzxyzzyyxx00 xyzzyyxxxyzzyyxxzxyzxyzzyyxxxyyyxx2023-2-1227/16200000 xyyyxxexyzzyyxxD线弹性问题本构方程平面应变 zxyzxyzzyyxxezxyzxyzzyyxxDxyyyxxxyyyxxE221000101)21)(1(平面应变状态 221000101)21)(1(EeD2023-2-1228/162线弹性问题本构方程轴对称 二维问题轴对称状态 00yzxz00yzxzzxyzxyzzyyxx00 xyzzyyxxxyzzyyxxzxyzxyzzyyxxxy
9、zzyyxx2023-2-1229/162线弹性问题本构方程轴对称 二维问题轴对称状态 00zr00zrzrzrzzrrzrzzrr00zrzzrrzrzrzzrrzrzzrr2023-2-1230/162线弹性问题本构方程轴对称 轴对称状态 zxyzxyzzyyxxezxyzxyzzyyxxDzrzrzzrrezrzrzzrrDzrzzrrezrzzrr0000D221000010101)21)(1(EeDzrzzrrzrzzrrE221000010101)21)(1(2023-2-1231/162线弹性问题本构方程一维问题 一维问题xxxxEEeD2023-2-1232/162常用单元模型
10、 n单元模型插值关系一一对应n单元类型一维单元、二维单元、三维单元等参单元、超参单元、次参单元2023-2-1233/162常用单元模型n一维单元 2节点线单元 1 2 1 2 3 1 2 3节点线单元梁单元2023-2-1234/162常用单元模型n二维单元3节点三角形线性单元 1 2 3 6节点三角形二次单元 1 2 3 5 6 4 2023-2-1235/162常用单元模型n二维单元10节点三角形三次单元 10 1 2 3 6 9 4 5 7 8 4节点四边形双线性单元 1 2 3 4 2023-2-1236/162常用单元模型n二维单元8节点四边形二次单元12节点四边形三次单元 1 2
11、 3 4 8 7 6 5 8 7 1 2 3 4 11 9 5 6 10 12 2023-2-1237/162常用单元模型n三维单元4节点四面体线性单元10节点四面体二次单元 1 2 3 4 1 10 9 8 4 7 2 3 6 5 2023-2-1238/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 2023-2-1239/162常用单元模型n准三维空间单元桁架单元一维2节点线单元+单元局部随体坐标系 为什么要建立单元局部随体坐标系？
12、简化分析问题的复杂程度。在局部坐标系中，空间桁架的每根杆每变成了一维2节点线单元2023-2-1240/162常用单元模型n准三维空间单元框架单元三维梁单元+一维2节点线单元+单元局部随体坐标系 两端都是刚性联结 可以要承受拉压、弯曲、扭转3种变形模式 框架单元的特点2023-2-1241/162常用单元模型n准三维空间单元板单元薄板单元中厚板单元弯曲和横向剪切2种变形模式抵抗板的变形如果板很薄，忽略横向剪切抗力，认为抵抗载荷的主要因素是弯矩2023-2-1242/162常用单元模型n准三维空间单元壳单元 抵抗拉压变形的二维单元+板单元+单元局部随体坐标系。适合于薄壳单元和中厚壳单元从几何上分
13、为薄壳单元和中厚壳单元 组合单元2023-2-1243/162常用单元模型n准三维空间单元 壳理论单元 由空间壳理论严格构造的壳单元。适合于薄壳单元和中厚壳单元 退化单元 由三维实体单元退化成的壳单元。只适合于中厚壳单元 2023-2-1244/162单元模型构造 n有限元法的基本思想 通过单元分片近似，在每个单元内假设近似函数来分片表示系统的场函数 n选择近似函数简单、实用的原则在有限元法中，近似函数称为插值函数 2023-2-1245/162单元模型构造n插值函数 一般都采用多项式函数，主要原因是:采用多项式插值函数比较容易推导单元平衡方程，特别是易于进行微分和积分运算。随着多项式函数阶次
14、的增加，可以提高有限元法的计算精度。从理论上说，无限提高多项式的阶数，可以求得系统的精确解。2023-2-1246/162单元模型构造方法 n整体坐标系法n局部坐标系法 nLagrange插值方法nHermite插值方法2023-2-1247/162单元模型构造方法n2节点线单元12 oxu1u2x1x2ux1.假设插值多项式xaaxu10)(2.利用节点值求 a0 和 a1 21021101xaauxaau12121xxuua1212210 xxxuxua2023-2-1248/162单元模型构造方法3.代入a0 和 a1，得插值多项式 u(x)xxxuuxxxuxuxu1212121221
15、)(4.按u1 和 u2合并同类项，设 l=x2-x1212121122112)(uuNNuulxxlxxulxxulxxxu2023-2-1249/162单元模型构造方法n关键 如何构造插值多项式 u?二维问题三维问题，如何构造插值多项式?2023-2-1250/162n收敛性条件 在单元内，场函数必须是连续的；完备性：插值多项式的阶次必须由低到高依次增加，不能出现跳跃现象；协调性：各单元边界必须连续，单元边界不能出现开裂现象。插值多项式收敛性条件 n收敛：当单元逐渐缩小时，如果插值多项式满足收敛性条件，则数值解将收敛于精确解 2023-2-1251/162插值多项式收敛性条件n协调单元 满
16、足插值多项式收敛性条件和的单元 n完备单元 满足插值多项式收敛性条件的单元ncr 阶连续性 插值多项式的第r阶导数是连续的 2023-2-1252/162插值多项式收敛性条件n非协调单元与部分协调单元 对于一般固体力学问题来说，协调性要求单元在变形时，相邻单元之间不应引起开裂、重叠或其它不连续现象。例如，梁、板、壳等单元，在单元边界不但要求位移是连续的，而且其一阶导数也必须是连续的。板、壳单元位移函数沿单元边界的法向导数（转角）的连续性一般比较难实现，因此出现了许多不完全满足协调性要求的“非协调单元”或“部分协调单元”，有时它们的精度也很好。2023-2-1253/162插值多项式选择条件 n
17、插值多项式应该尽可能满足其收敛性条件（收敛性）n由插值多项式所确定的场函数变化应该与局部坐标系的选择无关（各向同性）n假设的插值多项式系数的数量应该等于单元的节点数（解的唯一性）选择条件2023-2-1254/162插值多项式选择条件n深入分析由收敛性条件可知，插值多项式中必须含有常数项（刚体位移项），高阶项的次数必须依次增加，不允许有跳跃332210)(xxxxu3310)(xxxu2023-2-1255/162插值多项式选择条件由选择条件可知，插值多项式函数在所有自由度方向上要满足各向同性性，这样就不会随局部坐标系变化而改变了 n深入分析xyyxy,xu3210)(23210)(xyxy,
18、xu2023-2-1256/162插值多项式选择条件n深入分析选择条件是为了能由单元节点值唯一确定插值多项式 4节点四边形的插值多项式应该是 xyyxy,xu3210)(插值多项式系数i(i=0,1,2,3)也是4个 2023-2-1257/162单元模型构造整体坐标系法n基本思想 针对弹性体有限元网格建立一个统一的坐标系，每个单元的插值多项式都在这个坐标系上建立 y x o 弹性体 2023-2-1258/162单元模型构造整体坐标系法n2节点线单元12 oxu1u2x1x2ux1.假设插值多项式xaaxu10)(2.利用节点值求 a0 和 a1 21021101xaauxaau12121x
19、xuua1212210 xxxuxua2023-2-1259/162单元模型构造整体坐标系法3.代入a0 和 a1，得插值多项式 u(x)xxxuuxxxuxuxu1212121221)(4.按u1 和 u2合并同类项，设 l=x2-x1212121122112)(uuNNuulxxlxxulxxulxxxu2023-2-1260/162单元模型构造整体坐标系法lxxN2111 122122()euu xN uN uNNuNulxxN12N1 和 N2 称为单元的形函数；N 称为单元的形函数矩阵；ue 称为单元节点位移向量。n2节点线的单元形函数2023-2-1261/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 2023-2-1262/162单元模型构造整体坐标系法1.假设插值多项式yxyxu210),(2.首先，利用节点值求 0、1 和 2 n二维3节点三角形单元 323103222102121101yxuyxuyxuyxyxv210),(2023-2-1263/162单元模型构造整体坐标系法)(21)(21)(21332211233221113322110ucucucAubububAuauauaA3322111112
21、1yxyxyxA A为单元面积2023-2-1264/162单元模型构造整体坐标系法23132113311xxcyybyxyxa31213221122xxcyybyxyxa12321332233xxcyybyxyxa3.将 0、1 和 2 代入插值多项式，按u1、u2、u3合并同类项332211),(uNuNuNyxu2023-2-1265/162单元模型构造整体坐标系法)(21)(21)(21333322221111ycxbaANycxbaANycxbaAN4.同理可得332211),(vNvNvNyxv2023-2-1266/162单元模型构造整体坐标系法5.单元插值多项式为3322113
22、32211),(),(vNvNvNyxvuNuNuNyxu)(21)(21)(21333322221111ycxbaANycxbaANycxbaAN2023-2-1267/162单元模型构造整体坐标系法332211321321000000vuvuvuNNNNNNvu6.单元插值多项式写成矩阵形式（常用）2023-2-1268/162单元模型构造整体坐标系法321321321321000000vvvuuuNNNNNNvu7.单元插值多项式的另一种矩阵形式（不常用）2023-2-1269/162单元模型构造整体坐标系法n4节点四面体单元 2 1(u1,v1,w1)x y z 4 3(x1,y1,z
23、1)(x3,y3,z3)(x4,y4,z4)(u3,v3,w3)(u2,v2,w2)(x2,y2,z2)(u4,v4,w4)(x,y,z)(u,v,w)2023-2-1270/162单元模型构造整体坐标系法zyxzyxwzyxzyxvzyxzyxu421042104210),(),(),(1.假设插值多项式443322114433221144332211),(),(),(wNwNwNwNzyxwvNvNvNvNzyxvuNuNuNuNzyxu2.插值多项式为2023-2-1271/162单元模型构造整体坐标系法)(61zdycxbaVNiiiii（i=1,2,3,4）4443332221zyx
24、zyxzyxa 4433221111zyzyzyb4433221111zxzxzxc1114433221yxyxyxd 444333222111111161zyxzyxzyxzyxV 循环轮换脚标1、2、3、4，相应可以得到a2，b2，c2，d2、a3，b3，c3，d3、a4，b4，c4，d4 2023-2-1272/162单元模型构造整体坐标系法444333222111432143214321000000000000000000000000wvuwvuwvuwvuNNNNNNNNNNNNwvu3.单元插值多项式写成矩阵形式（常用）2023-2-1273/162单元模型构造整体坐标系法4.单元
25、插值多项式另一种矩阵形式（不常用）432143214321432143214321000000000000000000000000wwwwvvvvuuuuNNNNNNNNNNNNwvu2023-2-1274/162单元模型构造整体坐标系法n从理论上讲，整体坐标系法可以求任意单元的形函数，但计算过程太复杂n只能求一维2节点线单元、二维3节点三角形单元和三维4节点四面体单元3种简单单元的形函数n复杂的或二次以上的单元必须采用局部坐标系法求n位移场 u 是形函数 Ni 的线性组合，因此形函数Ni同样具有插值多项式的特性2023-2-1275/162单元刚度矩阵2节点线单元n一维2节点线单元n单元插值
26、关系 eNuu n单元几何关系 n单元本构关系 Lu DeN=N1 N2 De=E21uueuxLlxxNlxxN/)(/)(12212023-2-1276/162单元刚度矩阵2节点线单元n单元刚度矩阵vedvBDBkTvedvBDBkT x yzedxdydzBDBTxedxABDBTA为单元截面积；l为单元长度n矩阵BLNB LNB/21xNxN/1/1llBDBeAlT1111lAE/1/1/1/1llEllAl2023-2-1277/162单元刚度矩阵三角形单元n二维3角形单元n单元插值关系 eNuu T332211vuvuvueu321321000000NNNNNNN)(21)(21
27、)(21333322221111ycxbaANycxbaANycxbaAN2023-2-1278/162单元刚度矩阵三角形单元n单元几何关系 Lu Luvuxyyxxyyyxx002023-2-1279/162单元刚度矩阵三角形单元n单元本构关系 DexyyyxxxyyyxxE2100010112平面应力问题2023-2-1280/162单元刚度矩阵三角形单元n矩阵BLNB LNB xNyNxNyNxNyNyNyNyNxNxNxN33221132132100000033221132132100000021bcbcbccccbbbA2023-2-1281/162单元刚度矩阵三角形单元n单元刚度矩
28、阵vedvBDBkTvedvBDBkT x yzedxdydzBDBTx yedxdyhBDBTh为单元厚度BDBehATk为对称的6*6常数矩阵A为单元面积2023-2-1282/162作业n求4节点四面体单元的单元刚度矩阵 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)2023-2-1283/162单元模型构造整体坐标系法n单元形函数的特性正规性：单元形函数之和等于1。正交性：形函数在本节点的值等于1，在其它节点的值等于0。l
29、xxN21lxxN12例如:2节点线单元形函数2023-2-1284/162单元模型等参单元n等参单元 单元内任意一点的位移u与单元节点位移ue之间的关系为 eNuu 一般单元坐标的插值关系也采用与位移插值关系相同的变换关系即单元内任意一点的坐标x与单元节点坐标xe之间的关系为 eNxx 2023-2-1285/162单元模型等参单元n等参单元凡是几何形状和位移场采用同阶同参数插值关系来描述的单元，称为等参单元 前面介绍的所有单元都属于等参单元 在描述单元的几何形状和位移场时，并不一定非采用同阶插值关系 2023-2-1286/162单元模型等参单元332211321321000000vuvu
30、vuNNNNNNvu332211321321000000yxyxyxNNNNNNyxn等参单元3节点三角形等参单元 2023-2-1287/162单元模型等参单元n超参单元如果几何形状插值函数的阶数高于位移场插值函数的阶数，称为超参单元 n次参单元如果几何形状插值函数的阶数低于位移场插值函数的阶数，称为次参单元 2023-2-1288/162单元平衡方程组装过程 n为什么要组装?2 F 1 3 消除内力n组装的原则是什么?单元自由度与结构自由度对应2023-2-1289/162单元平衡方程组装过程 2 F 1 3 U3U4U2U1U5U6654321UUUUUUU结构自由度向量U2023-2-
31、1290/162单元平衡方程组装过程 2 1 U3U4U2U11u1u2u3u4 3 U611 U2U12u1u2U5u3u443214321UUUUuuuueu65214321UUUUuuuueu22023-2-1291/162单元平衡方程组装过程14131211432114414314214113413313213112412312212111411311211143214321ffffuuuukkkkkkkkkkkkkkkk2 1 U3U4U2U11u1u2u3u422023-2-1292/162单元平衡方程组装过程000000000000000000000065432165432114
32、131211654321144143142141134133132131124123122121114113112111ffffUUUUUUkkkkkkkkkkkkkkkk组装单元2023-2-1293/162单元平衡方程组装过程 3 U611 U2U12u1u2U5u3u424232221432124424324224123423323223122422322222121421321221165216521ffffuuuukkkkkkkkkkkkkkkk2023-2-1294/162单元平衡方程组装过程242314132212211165432124424324224123423323223
33、114414314214113413313213122422312412322212222112121421311411321211221111100000000654321654321ffffffffUUUUUUkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk再组装单元FKU 总体刚度方程 K 称为总体刚度矩阵 U 称为位移向量 F 称为载荷向量 2023-2-1295/162总体刚度矩阵K的特性 n对称性 n奇异性 n稀疏性 n非零元素带状分布 2023-2-1296/162约束处理过程 n为什么要约束处理?总体平衡方程组是奇异的消除无限制的刚体运动 使总体平衡方程组存在唯
34、一一组解2023-2-1297/162约束处理过程边界条件n 边界条件分类 力(载荷)边界条件位移边界条件 集中载荷力 表面分布力 自重力热交换引起的温度载荷 固定位移约束 强制位移约束 关联位移约束 002.U 005.U 011.U 514.U CkUU782023-2-1298/162约束处理过程模型简化xyU2023-2-1299/162约束处理过程模型简化yxUxyU2023-2-12100/162约束处理过程约束方程123456789101112yxU0.0321UUU0.012963VVVVUUUU1211102023-2-12101/162约束处理过程约束处理方法n位移约束处理
35、方法 赋0赋1法 乘大数法 2023-2-12102/162约束处理过程赋0赋1法n强制位移约束条件处理 U4=CFKU 654321654321666564636261565554535251464544434241363534333231262524232221161514131211FFFFFFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK2023-2-12103/162约束处理过程赋0赋1法n强制位移约束条件处理 U4=C6543216543216665646362615655545352514645444342413635343332312625
36、24232221161514131211FFFFFFUUUUUUKKCKKKKKKCKKKKKKCKKKKKKCKKKKKKCKKKKKKCKKKKCKFCKFCKFCKFCKFCKFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK6465454443432421416543216665636261565553525146454342413635333231262523222116151312110000002023-2-12104/162约束处理过程赋0赋1法n有6个方程，5个未知数，如果约束方程可以消除有限元平衡方程组的奇异性，则取任意5个方程联立求解，都会得到方程
37、组的唯一一组解。n系数矩阵由原来的对称的变成了非对称的，这对于大规模有限元方程组求解是十分不利的，采用相同的求解方法，在求解时间和矩阵存贮容量方面都增加了一倍。2023-2-12105/162约束处理过程赋0赋1法CKFCKFCKFCKFCKFCKFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK646545444343242141654321666563626156555352514645434241363533323126252322211615131211000000CKFCKFCKFCKFCKFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKK646
38、54534324214165432166656362615655535251363533323126252322211615131211000000000000为了保证系数矩阵的对称性，去掉方程组第4行2023-2-12106/162CKFCKFCKFCKFCKFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKK64654534324214165432166656362615655535251363533323126252322211615131211000000000000约束处理过程赋0赋1法引入强制位移约束方程 U4=C，使方程组求解时直接将自由度U4求出CKFCKFCCKF
39、CKFCKFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKK64654534324214165432166656362615655535251363533323126252322211615131211000010000002023-2-12107/162约束处理过程赋0赋1法n固定位移约束条件处理 U4=06532165432166656362615655535251363533323126252322211615131211000001000000FFFFFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKK2023-2-12108/162约束处理过程赋0赋1法n基
40、本原理 利用初等变换对求解方程组进行相同的行列变换，既保证方程组解不会改变，又可以保持方程组系数矩阵的对称性。在进行初等变换时，只要保证对方程组系数矩阵做相同的行列变换，就可以保持方程组系数矩阵的对称性。2023-2-12109/162约束处理过程乘大数法n乘大数法n基本原理 利用矩阵的初等变换不改变方程组解的思想。2023-2-12110/162654321654321666564636261565554535251464544434241363534333231262524232221161514131211FFFFFFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKKKKKK
41、KKKKKKK约束处理过程乘大数法n强制位移边界条件 CU 42023-2-12111/162约束处理过程乘大数法CU 4强制约束方程 ACAU 4A是一个大数，是系数矩阵中对角线元素K44的1010倍量级以上 为什么要乘以大数A？放大位移约束方程的优势2023-2-12112/162654321654321666564636261565554535251464544434241363534333231262524232221161514131211FFACFFFFUUUUUUKKKKKKKKKKKKKKAKKKKKKKKKKKKKKKKKKKKKK约束处理过程乘大数法n强制位移边界条件 CU
42、 42023-2-12113/162约束处理过程乘大数法n固定位移边界条件 C=0 约束后的方程组简化为 654321654321666564636261565554535251464544434241363534333231262524232221161514131211FFFFFFUUUUUUKKKKKKKKKKKKKKAKKKKKKKKKKKKKKKKKKKKKK654321654321666564636261565554535251464544434241363534333231262524232221161514131211FFACFFFFUUUUUUKKKKKKKKKKKKKKAK
43、KKKKKKKKKKKKKKKKKKKKK2023-2-12114/162约束处理过程两种方法比较n赋0赋1法在约束处理过程中是严格精确的，而乘大数法是一种近似约束处理方法，它的精度取决于所乘大数A值 两种方法都可以消除有限元平衡方程的奇异性，得到符合实际边界条件的唯一一组解。但两种方法还是有很大的区别 2023-2-12115/162约束处理过程两种方法比较n采用乘大数法约束处理后的有限元平衡方程在求解时可能造成解的失真，大数A值越大可能解的偏差会越大，而赋0赋1法就不会出现类似的问题，它在约束过程和求解过程都是精确的n乘大数法相对于赋0赋1法在约束处理过程上简单一些 2023-2-1211
44、6/162约束处理过程两种方法比较n赋0赋1法实际上是将关联位移约束方程代入到有限元平衡方程中的，是代入法。而乘大数是将占绝对优势的关联位移约束方程合并到有限元平衡方程中的，是罚方法，计算误差来自于合并过程，计算精度取决于关联位移约束方程的优势大小n商业软件中，位移边界条件的约束处理都采用赋0赋1法，乘大数很少被采用主要原因是它是一种近似方法，而且大数的大小也不好确定，有时还会造成求解失败 2023-2-12117/162约束处理过程弹簧单元假设柔性弹簧 kOXYU4f f=kU4k2023-2-12118/162弹簧约束方程 f=kU46543216543216665646362615655
45、54535251464544434241363534333231262524232221161514131211FFFFFFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK654321654321666564636261565554535251464544434241363534333231262524232221161514131211FFfFFFFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK654432165432166656463626156555453525146454443424136353433323126
46、2524232221161514131211FFkUFFFFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK654321654321666564636261565554535251464544434241363534333231262524232221161514131211FFFFFFUUUUUUKKKKKKKKKKKKKKkKKKKKKKKKKKKKKKKKKKKKK约束处理过程弹簧单元2023-2-12119/162方程组求解过程特点n方程组求解有限元计算过程中很重要的一部分，在有限元法的发展过程中，有限元方程的求解效率一直是其应用的最大瓶颈之一 n有
47、限元方程组的特点：有限元方程组的系数矩阵具有对称、稀疏、带状分布以及正定、主元占优。有效地利用这些特点，以减少系数矩阵的存贮量，提高方程组求解效率 2023-2-12120/162方程组求解过程分类比较n线性方程组的解法主要分两大类:直接解法：以高斯消去法基础，以等带宽或变带宽方式存贮系数矩阵内元素，对于求解规模比较大的问题，要存贮的元素非常巨大。迭代解法：只需要存贮系数矩阵中非零元素，存贮量很小，一般是变带宽存贮量的20%或更少，有些算法的求解效率也非常高，适合求解大规模线性方程组。但是这种解法对接近病态的方程组很难保证收敛性。2023-2-12121/162方程组求解过程带宽定义n有限元方
48、程组系数矩阵是稀疏的、非零元素呈带状分布，带宽就是它的宽度，带宽的大小是由系统有限元网格的节点号排序决定的，具体求法是 带宽=(单元最大节点号之差+1)*节点自由度数 n带宽是网格节点标注方法直接决定的，不同标注方法带宽可能相关很大 2023-2-12122/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
49、 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 2023-2-12123/162方程组求解过程带宽n所示四边形网格的三种节点号标注方法，每个节点是2个自由度n结构的带宽分别是12，18，56，相差很大，其中12和56之间相差近5倍，这就意味着系数矩阵的存贮量也是相差5倍，因此，对于大规模复杂系统的节点号优化是十分必要的 2023-2-12124/162方程组求解过程系数矩阵存贮 n系数矩阵存贮 如果节点号排序优化的比较好，系数矩阵的存贮量就会减
50、少很多。根据系数矩阵的对称性，一般都是按半带宽存贮。n系数矩阵存贮的方法 二维等带宽存贮 一维变带宽存贮2023-2-12125/162方程组求解过程二维等带宽存贮 n二维等带宽存贮 D D 0 0 K K n 2023-2-12126/162方程组求解过程二维等带宽存贮n二维等带宽存贮消除了最大带宽以外的全部零元素，节省了系数矩阵元素的存贮量。但是由于取最大带宽为存贮范围，因此不能排除在带宽内的大量零元素。当系数矩阵的各行带宽变化不大时，适合采用二维等带宽存贮，方程组求解过程中系数矩阵元素的寻址也比较方便，求解效率较高。n当出现局部带宽特别大的情况时，采用二维等带宽存贮时，将由于局部带宽过大
