有限元和有限差分法基础超详细版本.ppt
2023/4/51/23关于有限元与有限差分法基础超详细版本第一张，PPT共一百六十二页，创作于2022年6月2023/4/52/162有限元法基础n有限元发展过程n有限元应用n有限元发展方向第二张，PPT共一百六十二页，创作于2022年6月2023/4/53/162有限元法的基本思想n基本思想 1）将连续的求解系统离散为一组由节点相互联在一起的单元组合体 2）在每个单元内假设近似函数来分片表示系统的求解场函数 第三张，PPT共一百六十二页，创作于2022年6月2023/4/54/162有限元法的基本思想第四张，PPT共一百六十二页，创作于2022年6月2023/4/55/162有限元法的基本思想第五张，PPT共一百六十二页，创作于2022年6月2023/4/56/162有限元法的基本思想第六张，PPT共一百六十二页，创作于2022年6月2023/4/57/162有限元法的基本思想第七张，PPT共一百六十二页，创作于2022年6月2023/4/58/162有限元法的基本思想n离散为单元网格的件仍然要保证是一个连续体，单元与单元之间没有裂缝、不能重叠，所有单元通过单元节点相互关联着n变形体无论产生多大的塑性变形，单元与单元之间依然不会产生裂缝、交叉和重叠，关联单元的节点也不能脱开第八张，PPT共一百六十二页，创作于2022年6月2023/4/59/162有限元法的基本思想n不合格单元单元裂缝单元裂缝单元重叠单元重叠第九张，PPT共一百六十二页，创作于2022年6月2023/4/510/162有限元法的基本思想n变形前后单元之间都是连续的变形前的网格变形前的网格变形后的网格变形后的网格第十张，PPT共一百六十二页，创作于2022年6月2023/4/511/162有限元法的基本思想n基本思想n通过在单元内假设不同的插值函数，建立不同的单元模型，适应各种各样的变形模式和受力模式XFXF第十一张，PPT共一百六十二页，创作于2022年6月2023/4/512/162有限元法的基本思想n有限元法分类 1）位移法：基于最小势能原理或虚功原理 2）力法：基于最小余能原理 3）混合法：基于修正余能原理 第十二张，PPT共一百六十二页，创作于2022年6月2023/4/513/162有限元法的基本思想n基本过程离散化过程 约束处理过程 单元平衡方程组装过程 应变、应力回代过程 方程组求解过程 第十三张，PPT共一百六十二页，创作于2022年6月2023/4/514/162离散化过程n最小势能原理 弹性体的势能为弹性体变形后所具有的内能 为弹性体所受的外力功 第十四张，PPT共一百六十二页，创作于2022年6月2023/4/515/162离散化过程 为弹性体的应变 为弹性体的应力 u为弹性体的可容位移 弹性体处于平衡状态时，其势能应为最小 0第十五张，PPT共一百六十二页，创作于2022年6月2023/4/516/162离散化过程n单元插值关系 n单元几何关系 n单元本构关系 N为单元形函数矩阵 L为单元几何微分算子 为单元弹性矩阵 单元节点自由度向量第十六张，PPT共一百六十二页，创作于2022年6月2023/4/517/162离散化过程B 称为应变矩阵 单元平衡方程或单元刚度方程 k 称为单元刚度矩阵 f 称为单元载荷向量 第十七张，PPT共一百六十二页，创作于2022年6月2023/4/518/162单元刚度矩阵的特性 n对称性 n奇异性 n主元恒正且对角占优 离散化过程第十八张，PPT共一百六十二页，创作于2022年6月2023/4/519/162线弹性问题几何方程三维问题 三维问题第十九张，PPT共一百六十二页，创作于2022年6月2023/4/520/162线弹性问题几何方程二维问题 二维问题平面应力和平面应变状态 第二十张，PPT共一百六十二页，创作于2022年6月2023/4/521/162线弹性问题几何方程二维问题 二维问题轴对称状态 第二十一张，PPT共一百六十二页，创作于2022年6月2023/4/522/162线弹性问题几何方程一维问题 一维问题第二十二张，PPT共一百六十二页，创作于2022年6月2023/4/523/162线弹性问题本构方程三维问题 三维问题E为弹性模量；为泊松比 第二十三张，PPT共一百六十二页，创作于2022年6月2023/4/524/162线弹性问题本构方程平面应力 二维问题平面应力状态 第二十四张，PPT共一百六十二页，创作于2022年6月2023/4/525/162线弹性问题本构方程平面应力 平面应力状态 第二十五张，PPT共一百六十二页，创作于2022年6月2023/4/526/162线弹性问题本构方程平面应变 二维问题平面应变状态 第二十六张，PPT共一百六十二页，创作于2022年6月2023/4/527/162线弹性问题本构方程平面应变 平面应变状态 第二十七张，PPT共一百六十二页，创作于2022年6月2023/4/528/162线弹性问题本构方程轴对称 二维问题轴对称状态 第二十八张，PPT共一百六十二页，创作于2022年6月2023/4/529/162线弹性问题本构方程轴对称 二维问题轴对称状态 第二十九张，PPT共一百六十二页，创作于2022年6月2023/4/530/162线弹性问题本构方程轴对称 轴对称状态 第三十张，PPT共一百六十二页，创作于2022年6月2023/4/531/162线弹性问题本构方程一维问题 一维问题第三十一张，PPT共一百六十二页，创作于2022年6月2023/4/532/162常用单元模型 n单元模型插值关系一一对应n单元类型一维单元、二维单元、三维单元等参单元、超参单元、次参单元第三十二张，PPT共一百六十二页，创作于2022年6月2023/4/533/162常用单元模型n一维单元 2节点线单元3节点线单元梁单元第三十三张，PPT共一百六十二页，创作于2022年6月2023/4/534/162常用单元模型n二维单元3节点三角形线性单元6节点三角形二次单元第三十四张，PPT共一百六十二页，创作于2022年6月2023/4/535/162常用单元模型n二维单元10节点三角形三次单元4节点四边形双线性单元第三十五张，PPT共一百六十二页，创作于2022年6月2023/4/536/162常用单元模型n二维单元8节点四边形二次单元12节点四边形三次单元第三十六张，PPT共一百六十二页，创作于2022年6月2023/4/537/162常用单元模型n三维单元4节点四面体线性单元10节点四面体二次单元第三十七张，PPT共一百六十二页，创作于2022年6月2023/4/538/162常用单元模型n三维单元8节点六面体线性单元20节点六面体二次单元第三十八张，PPT共一百六十二页，创作于2022年6月2023/4/539/162常用单元模型n准三维空间单元桁架单元一维2节点线单元+单元局部随体坐标系 为什么要建立单元局部随体坐标系？1.简化分析问题的复杂程度。2.在局部坐标系中，空间桁架的每根杆每变成了一维2节点线单元第三十九张，PPT共一百六十二页，创作于2022年6月2023/4/540/162常用单元模型n准三维空间单元框架单元三维梁单元+一维2节点线单元+单元局部随体坐标系 两端都是刚性联结 可以要承受拉压、弯曲、扭转3种变形模式 框架单元的特点第四十张，PPT共一百六十二页，创作于2022年6月2023/4/541/162常用单元模型n准三维空间单元板单元薄板单元中厚板单元弯曲和横向剪切2种变形模式抵抗板的变形如果板很薄，忽略横向剪切抗力，认为抵抗载荷的主要因素是弯矩第四十一张，PPT共一百六十二页，创作于2022年6月2023/4/542/162常用单元模型n准三维空间单元壳单元 抵抗拉压变形的二维单元+板单元+单元局部随体坐标系。适合于薄壳单元和中厚壳单元从几何上分为薄壳单元和中厚壳单元 组合单元第四十二张，PPT共一百六十二页，创作于2022年6月2023/4/543/162常用单元模型n准三维空间单元 壳理论单元 由空间壳理论严格构造的壳单元。适合于薄壳单元和中厚壳单元 退化单元 由三维实体单元退化成的壳单元。只适合于中厚壳单元 第四十三张，PPT共一百六十二页，创作于2022年6月2023/4/544/162单元模型构造 n有限元法的基本思想 通过单元分片近似，在每个单元内假设近似函数来分片表示系统的场函数 n选择近似函数简单、实用的原则在有限元法中，近似函数称为插值函数 第四十四张，PPT共一百六十二页，创作于2022年6月2023/4/545/162单元模型构造n插值函数 一般都采用多项式函数，主要原因是:采用多项式插值函数比较容易推导单元平衡方程，特别是易于进行微分和积分运算。随着多项式函数阶次的增加，可以提高有限元法的计算精度。从理论上说，无限提高多项式的阶数，可以求得系统的精确解。第四十五张，PPT共一百六十二页，创作于2022年6月2023/4/546/162单元模型构造方法 n整体坐标系法n局部坐标系法 nLagrange插值方法nHermite插值方法第四十六张，PPT共一百六十二页，创作于2022年6月2023/4/547/162单元模型构造方法n2节点线单元12 oxu1u2x1x2ux1.假设插值多项式2.利用节点值求 a0 和 a1 第四十七张，PPT共一百六十二页，创作于2022年6月2023/4/548/162单元模型构造方法3.代入a0 和 a1，得插值多项式 u(x)4.按u1 和 u2合并同类项，设 l=x2-x1第四十八张，PPT共一百六十二页，创作于2022年6月2023/4/549/162单元模型构造方法n关键 如何构造插值多项式 u?二维问题三维问题，如何构造插值多项式?第四十九张，PPT共一百六十二页，创作于2022年6月2023/4/550/162n收敛性条件 在单元内，场函数必须是连续的；完备性：插值多项式的阶次必须由低到高依次增加，不能出现跳跃现象；协调性：各单元边界必须连续，单元边界不能出现开裂现象。插值多项式收敛性条件 n收敛：当单元逐渐缩小时，如果插值多项式满足收敛性条件，则数值解将收敛于精确解 第五十张，PPT共一百六十二页，创作于2022年6月2023/4/551/162插值多项式收敛性条件n协调单元 满足插值多项式收敛性条件和的单元 n完备单元 满足插值多项式收敛性条件的单元ncr 阶连续性 插值多项式的第r阶导数是连续的 第五十一张，PPT共一百六十二页，创作于2022年6月2023/4/552/162插值多项式收敛性条件n非协调单元与部分协调单元 对于一般固体力学问题来说，协调性要求单元在变形时，相邻单元之间不应引起开裂、重叠或其它不连续现象。例如，梁、板、壳等单元，在单元边界不但要求位移是连续的，而且其一阶导数也必须是连续的。板、壳单元位移函数沿单元边界的法向导数（转角）的连续性一般比较难实现，因此出现了许多不完全满足协调性要求的“非协调单元”或“部分协调单元”，有时它们的精度也很好。第五十二张，PPT共一百六十二页，创作于2022年6月2023/4/553/162插值多项式选择条件 n插值多项式应该尽可能满足其收敛性条件（收敛性）n由插值多项式所确定的场函数变化应该与局部坐标系的选择无关（各向同性）n假设的插值多项式系数的数量应该等于单元的节点数（解的唯一性）选择条件第五十三张，PPT共一百六十二页，创作于2022年6月2023/4/554/162插值多项式选择条件n深入分析由收敛性条件可知，插值多项式中必须含有常数项（刚体位移项），高阶项的次数必须依次增加，不允许有跳跃第五十四张，PPT共一百六十二页，创作于2022年6月2023/4/555/162插值多项式选择条件由选择条件可知，插值多项式函数在所有自由度方向上要满足各向同性性，这样就不会随局部坐标系变化而改变了 n深入分析第五十五张，PPT共一百六十二页，创作于2022年6月2023/4/556/162插值多项式选择条件n深入分析选择条件是为了能由单元节点值唯一确定插值多项式 4节点四边形的插值多项式应该是 插值多项式系数i(i=0,1,2,3)也是4个 第五十六张，PPT共一百六十二页，创作于2022年6月2023/4/557/162单元模型构造整体坐标系法n基本思想 针对弹性体有限元网格建立一个统一的坐标系，每个单元的插值多项式都在这个坐标系上建立 第五十七张，PPT共一百六十二页，创作于2022年6月2023/4/558/162单元模型构造整体坐标系法n2节点线单元12 oxu1u2x1x2ux1.假设插值多项式2.利用节点值求 a0 和 a1 第五十八张，PPT共一百六十二页，创作于2022年6月2023/4/559/162单元模型构造整体坐标系法3.代入a0 和 a1，得插值多项式 u(x)4.按u1 和 u2合并同类项，设 l=x2-x1第五十九张，PPT共一百六十二页，创作于2022年6月2023/4/560/162单元模型构造整体坐标系法N1 和 N2 称为单元的形函数；N 称为单元的形函数矩阵；ue 称为单元节点位移向量。n2节点线的单元形函数第六十张，PPT共一百六十二页，创作于2022年6月2023/4/561/162单元模型构造整体坐标系法n二维3节点三角形单元 建立整体坐标系oxy 第六十一张，PPT共一百六十二页，创作于2022年6月2023/4/562/162单元模型构造整体坐标系法1.假设插值多项式2.首先，利用节点值求 0、1 和 2 n二维3节点三角形单元 第六十二张，PPT共一百六十二页，创作于2022年6月2023/4/563/162单元模型构造整体坐标系法A为单元面积第六十三张，PPT共一百六十二页，创作于2022年6月2023/4/564/162单元模型构造整体坐标系法3.将 0、1 和 2 代入插值多项式，按u1、u2、u3合并同类项第六十四张，PPT共一百六十二页，创作于2022年6月2023/4/565/162单元模型构造整体坐标系法4.同理可得第六十五张，PPT共一百六十二页，创作于2022年6月2023/4/566/162单元模型构造整体坐标系法5.单元插值多项式为第六十六张，PPT共一百六十二页，创作于2022年6月2023/4/567/162单元模型构造整体坐标系法6.单元插值多项式写成矩阵形式（常用）第六十七张，PPT共一百六十二页，创作于2022年6月2023/4/568/162单元模型构造整体坐标系法7.单元插值多项式的另一种矩阵形式（不常用）第六十八张，PPT共一百六十二页，创作于2022年6月2023/4/569/162单元模型构造整体坐标系法n4节点四面体单元第六十九张，PPT共一百六十二页，创作于2022年6月2023/4/570/162单元模型构造整体坐标系法1.假设插值多项式2.插值多项式为第七十张，PPT共一百六十二页，创作于2022年6月2023/4/571/162单元模型构造整体坐标系法（i=1,2,3,4）循环轮换脚标1、2、3、4，相应可以得到a2，b2，c2，d2、a3，b3，c3，d3、a4，b4，c4，d4 第七十一张，PPT共一百六十二页，创作于2022年6月2023/4/572/162单元模型构造整体坐标系法3.单元插值多项式写成矩阵形式（常用）第七十二张，PPT共一百六十二页，创作于2022年6月2023/4/573/162单元模型构造整体坐标系法4.单元插值多项式另一种矩阵形式（不常用）第七十三张，PPT共一百六十二页，创作于2022年6月2023/4/574/162单元模型构造整体坐标系法n从理论上讲，整体坐标系法可以求任意单元的形函数，但计算过程太复杂n只能求一维2节点线单元、二维3节点三角形单元和三维4节点四面体单元3种简单单元的形函数n复杂的或二次以上的单元必须采用局部坐标系法求n位移场 u 是形函数 Ni 的线性组合，因此形函数Ni同样具有插值多项式的特性第七十四张，PPT共一百六十二页，创作于2022年6月2023/4/575/162单元刚度矩阵2节点线单元n一维2节点线单元n单元插值关系 n单元几何关系 n单元本构关系 N=N1 N2 De=E第七十五张，PPT共一百六十二页，创作于2022年6月2023/4/576/162单元刚度矩阵2节点线单元n单元刚度矩阵A为单元截面积；l为单元长度n矩阵B第七十六张，PPT共一百六十二页，创作于2022年6月2023/4/577/162单元刚度矩阵三角形单元n二维3角形单元n单元插值关系 第七十七张，PPT共一百六十二页，创作于2022年6月2023/4/578/162单元刚度矩阵三角形单元n单元几何关系 第七十八张，PPT共一百六十二页，创作于2022年6月2023/4/579/162单元刚度矩阵三角形单元n单元本构关系 平面应力问题第七十九张，PPT共一百六十二页，创作于2022年6月2023/4/580/162单元刚度矩阵三角形单元n矩阵B第八十张，PPT共一百六十二页，创作于2022年6月2023/4/581/162单元刚度矩阵三角形单元n单元刚度矩阵h为单元厚度k为对称的6*6常数矩阵A为单元面积第八十一张，PPT共一百六十二页，创作于2022年6月2023/4/582/162作业n求4节点四面体单元的单元刚度矩阵第八十二张，PPT共一百六十二页，创作于2022年6月2023/4/583/162单元模型构造整体坐标系法n单元形函数的特性正规性：单元形函数之和等于1。正交性：形函数在本节点的值等于1，在其它节点的值等于0。例如:2节点线单元形函数第八十三张，PPT共一百六十二页，创作于2022年6月2023/4/584/162单元模型等参单元n等参单元 单元内任意一点的位移u与单元节点位移ue之间的关系为 一般单元坐标的插值关系也采用与位移插值关系相同的变换关系即单元内任意一点的坐标x与单元节点坐标xe之间的关系为 第八十四张，PPT共一百六十二页，创作于2022年6月2023/4/585/162单元模型等参单元n等参单元凡是几何形状和位移场采用同阶同参数插值关系来描述的单元，称为等参单元 前面介绍的所有单元都属于等参单元 在描述单元的几何形状和位移场时，并不一定非采用同阶插值关系 第八十五张，PPT共一百六十二页，创作于2022年6月2023/4/586/162单元模型等参单元n等参单元3节点三角形等参单元 第八十六张，PPT共一百六十二页，创作于2022年6月2023/4/587/162单元模型等参单元n超参单元如果几何形状插值函数的阶数高于位移场插值函数的阶数，称为超参单元 n次参单元如果几何形状插值函数的阶数低于位移场插值函数的阶数，称为次参单元 第八十七张，PPT共一百六十二页，创作于2022年6月2023/4/588/162单元平衡方程组装过程 n为什么要组装?消除内力n组装的原则是什么?单元自由度与结构自由度对应第八十八张，PPT共一百六十二页，创作于2022年6月2023/4/589/162单元平衡方程组装过程 2 F 1 3 U3U4U2U1U5U6结构自由度向量U第八十九张，PPT共一百六十二页，创作于2022年6月2023/4/590/162单元平衡方程组装过程 2 1 U3U4U2U11u1u2u3u4 3 U611 U2U12u1u2U5u3u42第九十张，PPT共一百六十二页，创作于2022年6月2023/4/591/162单元平衡方程组装过程2 1 U3U4U2U11u1u2u3u42第九十一张，PPT共一百六十二页，创作于2022年6月2023/4/592/162单元平衡方程组装过程组装单元第九十二张，PPT共一百六十二页，创作于2022年6月2023/4/593/162单元平衡方程组装过程 3 U611 U2U12u1u2U5u3u4第九十三张，PPT共一百六十二页，创作于2022年6月2023/4/594/162单元平衡方程组装过程再组装单元总体刚度方程 K 称为总体刚度矩阵 U 称为位移向量 F 称为载荷向量 第九十四张，PPT共一百六十二页，创作于2022年6月2023/4/595/162总体刚度矩阵K的特性 n对称性 n奇异性 n稀疏性 n非零元素带状分布 第九十五张，PPT共一百六十二页，创作于2022年6月2023/4/596/162约束处理过程 n为什么要约束处理?总体平衡方程组是奇异的消除无限制的刚体运动 使总体平衡方程组存在唯一一组解第九十六张，PPT共一百六十二页，创作于2022年6月2023/4/597/162约束处理过程边界条件n 边界条件分类 力(载荷)边界条件位移边界条件 集中载荷力 表面分布力 自重力热交换引起的温度载荷 固定位移约束 强制位移约束 关联位移约束 第九十七张，PPT共一百六十二页，创作于2022年6月2023/4/598/162约束处理过程模型简化xy第九十八张，PPT共一百六十二页，创作于2022年6月2023/4/599/162约束处理过程模型简化yx xy第九十九张，PPT共一百六十二页，创作于2022年6月2023/4/5100/162约束处理过程约束方程123456789101112yx第一百张，PPT共一百六十二页，创作于2022年6月2023/4/5101/162约束处理过程约束处理方法n位移约束处理方法 赋0赋1法 乘大数法 第一百零一张，PPT共一百六十二页，创作于2022年6月2023/4/5102/162约束处理过程赋0赋1法n强制位移约束条件处理 U4=C第一百零二张，PPT共一百六十二页，创作于2022年6月2023/4/5103/162约束处理过程赋0赋1法n强制位移约束条件处理 U4=C第一百零三张，PPT共一百六十二页，创作于2022年6月2023/4/5104/162约束处理过程赋0赋1法n有6个方程，5个未知数，如果约束方程可以消除有限元平衡方程组的奇异性，则取任意5个方程联立求解，都会得到方程组的唯一一组解。n系数矩阵由原来的对称的变成了非对称的，这对于大规模有限元方程组求解是十分不利的，采用相同的求解方法，在求解时间和矩阵存贮容量方面都增加了一倍。第一百零四张，PPT共一百六十二页，创作于2022年6月2023/4/5105/162约束处理过程赋0赋1法为了保证系数矩阵的对称性，去掉方程组第4行第一百零五张，PPT共一百六十二页，创作于2022年6月2023/4/5106/162约束处理过程赋0赋1法引入强制位移约束方程 U4=C，使方程组求解时直接将自由度U4求出第一百零六张，PPT共一百六十二页，创作于2022年6月2023/4/5107/162约束处理过程赋0赋1法n固定位移约束条件处理 U4=0第一百零七张，PPT共一百六十二页，创作于2022年6月2023/4/5108/162约束处理过程赋0赋1法n基本原理 利用初等变换对求解方程组进行相同的行列变换，既保证方程组解不会改变，又可以保持方程组系数矩阵的对称性。在进行初等变换时，只要保证对方程组系数矩阵做相同的行列变换，就可以保持方程组系数矩阵的对称性。第一百零八张，PPT共一百六十二页，创作于2022年6月2023/4/5109/162约束处理过程乘大数法n乘大数法n基本原理 利用矩阵的初等变换不改变方程组解的思想。第一百零九张，PPT共一百六十二页，创作于2022年6月2023/4/5110/162约束处理过程乘大数法n强制位移边界条件 第一百一十张，PPT共一百六十二页，创作于2022年6月2023/4/5111/162约束处理过程乘大数法强制约束方程 A是一个大数，是系数矩阵中对角线元素K44的1010倍量级以上 为什么要乘以大数A？放大位移约束方程的优势第一百一十一张，PPT共一百六十二页，创作于2022年6月2023/4/5112/162约束处理过程乘大数法n强制位移边界条件 第一百一十二张，PPT共一百六十二页，创作于2022年6月2023/4/5113/162约束处理过程乘大数法n固定位移边界条件 C=0 约束后的方程组简化为 第一百一十三张，PPT共一百六十二页，创作于2022年6月2023/4/5114/162约束处理过程两种方法比较n赋0赋1法在约束处理过程中是严格精确的，而乘大数法是一种近似约束处理方法，它的精度取决于所乘大数A值 两种方法都可以消除有限元平衡方程的奇异性，得到符合实际边界条件的唯一一组解。但两种方法还是有很大的区别 第一百一十四张，PPT共一百六十二页，创作于2022年6月2023/4/5115/162约束处理过程两种方法比较n采用乘大数法约束处理后的有限元平衡方程在求解时可能造成解的失真，大数A值越大可能解的偏差会越大，而赋0赋1法就不会出现类似的问题，它在约束过程和求解过程都是精确的n乘大数法相对于赋0赋1法在约束处理过程上简单一些 第一百一十五张，PPT共一百六十二页，创作于2022年6月2023/4/5116/162约束处理过程两种方法比较n赋0赋1法实际上是将关联位移约束方程代入到有限元平衡方程中的，是代入法。而乘大数是将占绝对优势的关联位移约束方程合并到有限元平衡方程中的，是罚方法，计算误差来自于合并过程，计算精度取决于关联位移约束方程的优势大小n商业软件中，位移边界条件的约束处理都采用赋0赋1法，乘大数很少被采用主要原因是它是一种近似方法，而且大数的大小也不好确定，有时还会造成求解失败 第一百一十六张，PPT共一百六十二页，创作于2022年6月2023/4/5117/162约束处理过程弹簧单元假设柔性弹簧 kOXYU4f f=kU4k第一百一十七张，PPT共一百六十二页，创作于2022年6月2023/4/5118/162弹簧约束方程 f=kU4约束处理过程弹簧单元第一百一十八张，PPT共一百六十二页，创作于2022年6月2023/4/5119/162方程组求解过程特点n方程组求解有限元计算过程中很重要的一部分，在有限元法的发展过程中，有限元方程的求解效率一直是其应用的最大瓶颈之一 n有限元方程组的特点：有限元方程组的系数矩阵具有对称、稀疏、带状分布以及正定、主元占优。有效地利用这些特点，以减少系数矩阵的存贮量，提高方程组求解效率 第一百一十九张，PPT共一百六十二页，创作于2022年6月2023/4/5120/162方程组求解过程分类比较n线性方程组的解法主要分两大类:直接解法：以高斯消去法基础，以等带宽或变带宽方式存贮系数矩阵内元素，对于求解规模比较大的问题，要存贮的元素非常巨大。迭代解法：只需要存贮系数矩阵中非零元素，存贮量很小，一般是变带宽存贮量的20%或更少，有些算法的求解效率也非常高，适合求解大规模线性方程组。但是这种解法对接近病态的方程组很难保证收敛性。第一百二十张，PPT共一百六十二页，创作于2022年6月2023/4/5121/162方程组求解过程带宽定义n有限元方程组系数矩阵是稀疏的、非零元素呈带状分布，带宽就是它的宽度，带宽的大小是由系统有限元网格的节点号排序决定的，具体求法是 带宽=(单元最大节点号之差+1)*节点自由度数 n带宽是网格节点标注方法直接决定的，不同标注方法带宽可能相关很大 第一百二十一张，PPT共一百六十二页，创作于2022年6月2023/4/5122/162方程组求解过程带宽n带宽是网格节点标注方法直接决定的，不同标注方法带宽可能相关很大 第一百二十二张，PPT共一百六十二页，创作于2022年6月2023/4/5123/162方程组求解过程带宽n所示四边形网格的三种节点号标注方法，每个节点是2个自由度n结构的带宽分别是12，18，56，相差很大，其中12和56之间相差近5倍，这就意味着系数矩阵的存贮量也是相差5倍，因此，对于大规模复杂系统的节点号优化是十分必要的 第一百二十三张，PPT共一百六十二页，创作于2022年6月2023/4/5124/162方程组求解过程系数矩阵存贮 n系数矩阵存贮 如果节点号排序优化的比较好，系数矩阵的存贮量就会减少很多。根据系数矩阵的对称性，一般都是按半带宽存贮。n系数矩阵存贮的方法 二维等带宽存贮 一维变带宽存贮第一百二十四张，PPT共一百六十二页，创作于2022年6月2023/4/5125/162方程组求解过程二维等带宽存贮 n二维等带宽存贮 第一百二十五张，PPT共一百六十二页，创作于2022年6月2023/4/5126/162方程组求解过程二维等带宽存贮n二维等带宽存贮消除了最大带宽以外的全部零元素，节省了系数矩阵元素的存贮量。但是由于取最大带宽为存贮范围，因此不能排除在带宽内的大量零元素。当系数矩阵的各行带宽变化不大时，适合采用二维等带宽存贮，方程组求解过程中系数矩阵元素的寻址也比较方便，求解效率较高。n当出现局部带宽特别大的情况时，采用二维等带宽存贮时，将由于局部带宽过大而使整体系数矩阵的存贮大大增加。第一百二十六张，PPT共一百六十二页，创作于2022年6月2023/4/5127/162方程组求解过程一维变带宽存贮 n一维变带宽存贮 一维变带宽存贮方法就是把变化的带宽内的元素按一定的顺序存贮在一个一维数组中。由于它不按最大带宽存贮，因此比二维等带宽存贮更节省内存。按照解法可分为按行一维变带宽存贮和按列一维变带宽存贮。第一百二十七张，PPT共一百六十二页，创作于2022年6月2023/4/5128/162n按行一维变带宽存贮 方程组求解过程一维变带宽存贮 辅助的寻址数组M 第一百二十八张，PPT共一百六十二页，创作于2022年6月2023/4/5129/162n一维变带宽存贮是最节省内存的一种方法，但是由于要借助于寻址数组寻找系数矩阵元素的位置，相对二维等带宽存贮方法来说要复杂一些，而且在程序实现时也要复杂得多，方程组求解过程中也要消耗一些数组寻址时间。因此，在选用存贮方法时要权衡二者的利弊，统盘考虑。一般当带宽变化不大，计算机内存允许时，采用二维等带宽存贮方法是比较合适的。方程组求解过程一维变带宽存贮 第一百二十九张，PPT共一百六十二页，创作于2022年6月2023/4/5130/162方程组求解过程求解方法n方程组求解方法 高斯消去法 三角分解法 雅可比（Jacobi）迭代法 高斯-赛德尔（Gauss-Seidel）迭代法 第一百三十张，PPT共一百六十二页，创作于2022年6月2023/4/5131/162应变、应力回代过程 n单元应变和应力回代求解 通过求解有限元平衡方程得到有限元节点位移后，就可以进行系统的刚度校核。如果所分析问题要进行强度校核，就要回代求解单元的应变和应力。n由插值关系和几何关系可得单元应变，再通过本构关系得到单元应力第一百三十一张，PPT共一百六十二页，创作于2022年6月2023/4/5132/162有限差分法有限差分法第一百三十二张，PPT共一百六十二页，创作于2022年6月2023/4/5133/162 从弹性力学的基本方程建立以来，这些方程在各种问题的边界条件下如何求解，一直是很多数学工作者和力学工作者研究的内容。即弹性力学的经典解法存在一定的局限性，当弹性体的边界条件和受载情况复杂一点，往往无法求得偏微分方程的边值问题的解析解，许多工程重要问题，不能够得出函数式的解答。因此，弹性力学问题的各种数值解法便具有重要的实际意义。第一百三十三张，PPT共一百六十二页，创作于2022年6月2023/4/5134/162 工程中常用的数值解法有有限单元法和差分法。有限单元法 是以有限个单元的集合体来代替连续体，属于物理上的近似。差分法 是把弹性力学的基本方程和边界条件（一般均为微分方程）近似地改用差分方程（代数方程）来表示，把求解微分方程的问题改换成为求解代数方程的问题，属于数学上的近似。第一百三十四张，PPT共一百六十二页，创作于2022年6月2023/4/5135/162 第一节 差分方程 第二节 应力函数的差分解 第三节 深梁应力函数的差分解第一百三十五张，PPT共一百六十二页，创作于2022年6月2023/4/5136/162第一节 差分方程 差分法是沿用已久的一种数值解法。随着计算机的普及和相应的软件发展，此法成为解弹性力学问题的一种有效的方法。第一百三十六张，PPT共一百六十二页，创作于2022年6月2023/4/5137/162 我们在弹性体上，用相隔等间距h而平行于坐标轴的两组平行线织成正方形网格，x=y=h，如图。设f=f(x,y)为弹性体内的某一个连续函数。该函数在平行于x轴的一根网线上，如在-上,它只随x坐标的改变而变化。在邻近结点处，函数f可展为泰勒级数如下：第一百三十七张，PPT共一百六十二页，创作于2022年6月2023/4/5138/162 我们将只考虑离开结点充分近的那些结点，即（x-x0）充分小。于是可不计（x-x0）的三次及更高次幂的各项，则上式简写为：在结点，x=x0-h,在结点1,x=x0+h，代入(b)得：第一百三十八张，PPT共一百六十二页，创作于2022年6月联立(c),(d),解得差分公式：同理，在网线-上可得到差分公式第一百三十九张，PPT共一百六十二页，创作于2022年6月 差分公式(-)及(-)是以相隔2h的两结点处的函数值来表示中间结点处的一阶导数值，可称为中点导数公式。以相邻三结点处的函数值来表示一个端点处的一阶导数值，可称为端点导数公式。应当指出：中点导数公式与端点导数公式相比，精度较高。因为前者反映了结点两边的函数变化，而后者却只反映了结点一边的函数变化。因此，我们总是尽可能应用前者，而只有在无法应用前者时才不得不应用后者。第一百四十张，PPT共一百六十二页，创作于2022年6月 以上(-)(-)是基本差分公式,从而可导出其它的差分公式如下：第一百四十一张，PPT共一百六十二页，创作于2022年6月2023/4/5142/162第二节 应力函数的差分解 当不计体力时，我们已把弹性力学平面问题归结为在给定边界条件下求解双调和方程的问题。用差分法解平面问题，就应先将双调和方程变换为差分方程，而后求解之。第一百四十二张，PPT共一百六十二页，创作于2022年6月 一旦求得弹性体全部节点的值后，就可按应力分量差分公式（对节点0）算得弹性体各节点的应力。第一百四十三张，PPT共一百六十二页，创作于2022年6月 可见，用差分法解平面问题，共有两大任务：一、建立差分方程 将（1-68）代入双调和方程 对于弹性体边界以内的每一结点，都可以建立这样一个差分方程。整理即得第一百四十四张，PPT共一百六十二页，创作于2022年6月2023/4/5145/162 二、联立求解这些线性代数方程，就能求得各内结点处的值。为了求得边界上各结点处的值，须要应用应力边界条件，即：一般建立和求解差分方程，在数学上不会遇到很大困难。但是，当对于边界内一行的（距边界为h的）结点，建立的差分方程还将涉及边界上各结点处的值，并包含边界外一行的虚结点处的值。第一百四十五张，PPT共一百六十二页，创作于2022年6月代入上式，即得：l1=cos(N,x)=cos=dy/ds,l2=cos(N,y)=sin=-dx/ds,于是，式（a）可改写为：由右图可见，第一百四十六张，PPT共一百六十二页，创作于2022年6月关于边界上任一点处由此得：的值，可将（b）式从A点到B点对s积分得到：第一百四十七张，PPT共一百六十二页，创作于2022年6月 将此式亦从A点到B点沿s进行积分，就得到边界上任一点B处的值。为此利用分部积分法，得：由高等数学可知，第一百四十八张，PPT共一百六十二页，创作于2022年6月将式(b),(c)代入，整理得：由前知，把应力函数加上一个线性函数，并不影响应力。因此，可设想把应力函数加上a+bx+cy，然后调整a，b，c三个数值，使得由式(d)及式(c)可见，设即可根据面力分量及导数求得为已知，第一百四十九张，PPT共一百六十二页，创作于2022年6月 从图易看出，式（23）右边的积分式表示A与B之间的x方向的面力之和；式（24）右边的积分式表示A与B之间的y方向的面力之和；式（2
