资源描述
2024/5/61/23有限元法基础及有限差分法有限差分法基础有限元法有限差分法有限差分法有限元法基础及有限差分法基础有限元法有限差分法2024/5/62/162有限元法基础有限元发展过程有限元应用有限元发展方向有限元法基础有限元发展过程2024/5/63/162有限元法的基本思想基本思想 1）将连续的求解系统离散为一组由节点相互联在一起的单元组合体 2）在每个单元内假设近似函数来分片表示系统的求解场函数 有限元法的基本思想基本思想 1）将连续的求解系统离散为2024/5/64/162有限元法的基本思想有限元法的基本思想2024/5/65/162有限元法的基本思想有限元法的基本思想2024/5/66/162有限元法的基本思想有限元法的基本思想2024/5/67/162有限元法的基本思想有限元法的基本思想2024/5/68/162有限元法的基本思想离散为单元网格的件仍然要保证是一个连续体，单元与单元之间没有裂缝、不能重叠，所有单元通过单元节点相互关联着变形体无论产生多大的塑性变形，单元与单元之间依然不会产生裂缝、交叉和重叠，关联单元的节点也不能脱开有限元法的基本思想离散为单元网格的件仍然要保证是一个连续体，2024/5/69/162有限元法的基本思想不合格单元单元裂缝单元裂缝单元重叠单元重叠有限元法的基本思想不合格单元单元裂缝单元重叠2024/5/610/162有限元法的基本思想变形前后单元之间都是连续的变形前的网格变形前的网格变形后的网格变形后的网格有限元法的基本思想变形前后单元之间都是连续的变形前的网格变形2024/5/611/162有限元法的基本思想基本思想通过在单元内假设不同的插值函数，建立不同的单元模型，适应各种各样的变形模式和受力模式XFXF有限元法的基本思想基本思想XFXF2024/5/612/162有限元法的基本思想有限元法分类 1）位移法：基于最小势能原理或虚功原理 2）力法：基于最小余能原理 3）混合法：基于修正余能原理 有限元法的基本思想有限元法分类 1）位移法：基于最小势能2024/5/613/162有限元法的基本思想基本过程离散化过程 约束处理过程 单元平衡方程组装过程 应变、应力回代过程 方程组求解过程 有限元法的基本思想基本过程离散化过程 约束处理过程 单元平衡2024/5/614/162离散化过程最小势能原理 弹性体的势能为弹性体变形后所具有的内能 为弹性体所受的外力功 离散化过程最小势能原理 弹性体的势能为弹性体变形后所具有的内2024/5/615/162离散化过程 为弹性体的应变 为弹性体的应力 u为弹性体的可容位移 弹性体处于平衡状态时，其势能应为最小 0离散化过程 为弹性体的应变 为弹性体的应力 u为弹性体2024/5/616/162离散化过程单元插值关系 单元几何关系 单元本构关系 N为单元形函数矩阵 L为单元几何微分算子 为单元弹性矩阵 单元节点自由度向量离散化过程单元插值关系 单元几何关系 单元本构关系 N为单元2024/5/617/162离散化过程B 称为应变矩阵 单元平衡方程或单元刚度方程 k 称为单元刚度矩阵 f 称为单元载荷向量 离散化过程B 称为应变矩阵 单元平衡方程或单元刚度方程 k 2024/5/618/162单元刚度矩阵的特性 对称性 奇异性 主元恒正且对角占优 离散化过程单元刚度矩阵的特性 对称性 奇异性 主元恒正且对角占优 离散2024/5/619/162线弹性问题几何方程三维问题 三维问题线弹性问题几何方程三维问题 三维问题2024/5/620/162线弹性问题几何方程二维问题 二维问题平面应力和平面应变状态 线弹性问题几何方程二维问题 二维问题平面应力和平面应变状态2024/5/621/162线弹性问题几何方程二维问题 二维问题轴对称状态 线弹性问题几何方程二维问题 二维问题轴对称状态 2024/5/622/162线弹性问题几何方程一维问题 一维问题线弹性问题几何方程一维问题 一维问题2024/5/623/162线弹性问题本构方程三维问题 三维问题E为弹性模量；为泊松比 线弹性问题本构方程三维问题 三维问题E为弹性模量；为泊松2024/5/624/162线弹性问题本构方程平面应力 二维问题平面应力状态 线弹性问题本构方程平面应力 二维问题平面应力状态 2024/5/625/162线弹性问题本构方程平面应力 平面应力状态 线弹性问题本构方程平面应力 平面应力状态 2024/5/626/162线弹性问题本构方程平面应变 二维问题平面应变状态 线弹性问题本构方程平面应变 二维问题平面应变状态 2024/5/627/162线弹性问题本构方程平面应变 平面应变状态 线弹性问题本构方程平面应变 平面应变状态 2024/5/628/162线弹性问题本构方程轴对称 二维问题轴对称状态 线弹性问题本构方程轴对称 二维问题轴对称状态 2024/5/629/162线弹性问题本构方程轴对称 二维问题轴对称状态 线弹性问题本构方程轴对称 二维问题轴对称状态 2024/5/630/162线弹性问题本构方程轴对称 轴对称状态 线弹性问题本构方程轴对称 轴对称状态 2024/5/631/162线弹性问题本构方程一维问题 一维问题线弹性问题本构方程一维问题 一维问题2024/5/632/162常用单元模型 单元模型插值关系一一对应单元类型一维单元、二维单元、三维单元等参单元、超参单元、次参单元常用单元模型 单元模型插值关系一一对应单元类型一维单元、二维2024/5/633/162常用单元模型一维单元 2节点线单元3节点线单元梁单元常用单元模型一维单元 2节点线单元3节点线单元梁单元2024/5/634/162常用单元模型二维单元3节点三角形线性单元6节点三角形二次单元常用单元模型二维单元3节点三角形线性单元6节点三角形二次单元2024/5/635/162常用单元模型二维单元10节点三角形三次单元4节点四边形双线性单元常用单元模型二维单元10节点三角形三次单元4节点四边形双线性2024/5/636/162常用单元模型二维单元8节点四边形二次单元12节点四边形三次单元常用单元模型二维单元8节点四边形二次单元12节点四边形三次单2024/5/637/162常用单元模型三维单元4节点四面体线性单元10节点四面体二次单元常用单元模型三维单元4节点四面体线性单元10节点四面体二次单2024/5/638/162常用单元模型三维单元8节点六面体线性单元20节点六面体二次单元常用单元模型三维单元8节点六面体线性单元20节点六面体二次单2024/5/639/162常用单元模型准三维空间单元桁架单元一维2节点线单元+单元局部随体坐标系 为什么要建立单元局部随体坐标系？1.简化分析问题的复杂程度。2.在局部坐标系中，空间桁架的每根杆每变成了一维2节点线单元常用单元模型准三维空间单元桁架单元一维2节点线单元+单元局部2024/5/640/162常用单元模型准三维空间单元框架单元三维梁单元+一维2节点线单元+单元局部随体坐标系 两端都是刚性联结 可以要承受拉压、弯曲、扭转3种变形模式 框架单元的特点常用单元模型准三维空间单元框架单元三维梁单元+一维2节点线单2024/5/641/162常用单元模型准三维空间单元板单元薄板单元中厚板单元弯曲和横向剪切2种变形模式抵抗板的变形如果板很薄，忽略横向剪切抗力，认为抵抗载荷的主要因素是弯矩常用单元模型准三维空间单元板单元薄板单元中厚板单元弯曲和横向2024/5/642/162常用单元模型准三维空间单元壳单元 抵抗拉压变形的二维单元+板单元+单元局部随体坐标系。适合于薄壳单元和中厚壳单元从几何上分为薄壳单元和中厚壳单元 组合单元常用单元模型准三维空间单元壳单元 抵抗拉压变形的二维单元+板2024/5/643/162常用单元模型准三维空间单元 壳理论单元 由空间壳理论严格构造的壳单元。适合于薄壳单元和中厚壳单元 退化单元 由三维实体单元退化成的壳单元。只适合于中厚壳单元 常用单元模型准三维空间单元 壳理论单元 由空间壳理论严格构2024/5/644/162单元模型构造 有限元法的基本思想 通过单元分片近似，在每个单元内假设近似函数来分片表示系统的场函数 选择近似函数简单、实用的原则在有限元法中，近似函数称为插值函数 单元模型构造 有限元法的基本思想 通过单元分片近似，在每个单2024/5/645/162单元模型构造插值函数 一般都采用多项式函数，主要原因是:采用多项式插值函数比较容易推导单元平衡方程，特别是易于进行微分和积分运算。随着多项式函数阶次的增加，可以提高有限元法的计算精度。从理论上说，无限提高多项式的阶数，可以求得系统的精确解。单元模型构造插值函数 一般都采用多项式函数，主要原因是:采2024/5/646/162单元模型构造方法 整体坐标系法局部坐标系法 Lagrange插值方法Hermite插值方法单元模型构造方法 整体坐标系法Lagrange插值方法2024/5/647/162单元模型构造方法2节点线单元12 oxu1u2x1x2ux1.假设插值多项式2.利用节点值求 a0 和 a1 单元模型构造方法2节点线单元12 oxu1u2x1x2ux12024/5/648/162单元模型构造方法3.代入a0 和 a1，得插值多项式 u(x)4.按u1 和 u2合并同类项，设 l=x2-x1单元模型构造方法3.代入a0 和 a1，得插值多项式 u(2024/5/649/162单元模型构造方法关键 如何构造插值多项式 u?二维问题三维问题，如何构造插值多项式?单元模型构造方法关键 如何构造插值多项式 u?二维问2024/5/650/162收敛性条件 在单元内，场函数必须是连续的；完备性：插值多项式的阶次必须由低到高依次增加，不能出现跳跃现象；协调性：各单元边界必须连续，单元边界不能出现开裂现象。插值多项式收敛性条件 收敛：当单元逐渐缩小时，如果插值多项式满足收敛性条件，则数值解将收敛于精确解 收敛性条件 插值多项式收敛性条件 收敛：当单元逐渐缩小时，如2024/5/651/162插值多项式收敛性条件协调单元 满足插值多项式收敛性条件和的单元 完备单元 满足插值多项式收敛性条件的单元cr 阶连续性 插值多项式的第r阶导数是连续的 插值多项式收敛性条件协调单元 2024/5/652/162插值多项式收敛性条件非协调单元与部分协调单元 对于一般固体力学问题来说，协调性要求单元在变形时，相邻单元之间不应引起开裂、重叠或其它不连续现象。例如，梁、板、壳等单元，在单元边界不但要求位移是连续的，而且其一阶导数也必须是连续的。板、壳单元位移函数沿单元边界的法向导数（转角）的连续性一般比较难实现，因此出现了许多不完全满足协调性要求的“非协调单元”或“部分协调单元”，有时它们的精度也很好。插值多项式收敛性条件非协调单元与部分协调单元2024/5/653/162插值多项式选择条件 插值多项式应该尽可能满足其收敛性条件（收敛性）由插值多项式所确定的场函数变化应该与局部坐标系的选择无关（各向同性）假设的插值多项式系数的数量应该等于单元的节点数（解的唯一性）选择条件插值多项式选择条件 插值多项式应该尽可能满足其收敛性条件（收2024/5/654/162插值多项式选择条件深入分析由收敛性条件可知，插值多项式中必须含有常数项（刚体位移项），高阶项的次数必须依次增加，不允许有跳跃插值多项式选择条件深入分析由收敛性条件可知，插值多项式中必2024/5/655/162插值多项式选择条件由选择条件可知，插值多项式函数在所有自由度方向上要满足各向同性性，这样就不会随局部坐标系变化而改变了 深入分析插值多项式选择条件由选择条件可知，插值多项式函数在所有自由2024/5/656/162插值多项式选择条件深入分析选择条件是为了能由单元节点值唯一确定插值多项式 4节点四边形的插值多项式应该是 插值多项式系数i(i=0,1,2,3)也是4个 插值多项式选择条件深入分析选择条件是为了能由单元节点值唯一2024/5/657/162单元模型构造整体坐标系法基本思想 针对弹性体有限元网格建立一个统一的坐标系，每个单元的插值多项式都在这个坐标系上建立 单元模型构造整体坐标系法基本思想 针对弹性体有限元网格建立2024/5/658/162单元模型构造整体坐标系法2节点线单元12 oxu1u2x1x2ux1.假设插值多项式2.利用节点值求 a0 和 a1 单元模型构造整体坐标系法2节点线单元12 oxu1u2x12024/5/659/162单元模型构造整体坐标系法3.代入a0 和 a1，得插值多项式 u(x)4.按u1 和 u2合并同类项，设 l=x2-x1单元模型构造整体坐标系法3.代入a0 和 a1，得插值多2024/5/660/162单元模型构造整体坐标系法N1 和 N2 称为单元的形函数；N 称为单元的形函数矩阵；ue 称为单元节点位移向量。2节点线的单元形函数单元模型构造整体坐标系法N1 和 N2 称为单元的形函数；2024/5/661/162单元模型构造整体坐标系法二维3节点三角形单元 建立整体坐标系oxy 单元模型构造整体坐标系法二维3节点三角形单元 建立整体坐标2024/5/662/162单元模型构造整体坐标系法1.假设插值多项式2.首先，利用节点值求 0、1 和 2 二维3节点三角形单元 单元模型构造整体坐标系法1.假设插值多项式2.首先2024/5/663/162单元模型构造整体坐标系法A为单元面积单元模型构造整体坐标系法A为单元面积2024/5/664/162单元模型构造整体坐标系法3.将 0、1 和 2 代入插值多项式，按u1、u2、u3合并同类项单元模型构造整体坐标系法3.将 0、1 和 2024/5/665/162单元模型构造整体坐标系法4.同理可得单元模型构造整体坐标系法4.同理可得2024/5/666/162单元模型构造整体坐标系法5.单元插值多项式为单元模型构造整体坐标系法5.单元插值多项式为2024/5/667/162单元模型构造整体坐标系法6.单元插值多项式写成矩阵形式（常用）单元模型构造整体坐标系法6.单元插值多项式写成矩阵形式2024/5/668/162单元模型构造整体坐标系法7.单元插值多项式的另一种矩阵形式（不常用）单元模型构造整体坐标系法7.单元插值多项式的另一种矩阵2024/5/669/162单元模型构造整体坐标系法4节点四面体单元单元模型构造整体坐标系法4节点四面体单元2024/5/670/162单元模型构造整体坐标系法1.假设插值多项式2.插值多项式为单元模型构造整体坐标系法1.假设插值多项式2.插值2024/5/671/162单元模型构造整体坐标系法（i=1,2,3,4）循环轮换脚标1、2、3、4，相应可以得到a2，b2，c2，d2、a3，b3，c3，d3、a4，b4，c4，d4 单元模型构造整体坐标系法（i=1,2,3,4）循环轮换脚2024/5/672/162单元模型构造整体坐标系法3.单元插值多项式写成矩阵形式（常用）单元模型构造整体坐标系法3.单元插值多项式写成矩阵形式2024/5/673/162单元模型构造整体坐标系法4.单元插值多项式另一种矩阵形式（不常用）单元模型构造整体坐标系法4.单元插值多项式另一种矩阵形2024/5/674/162单元模型构造整体坐标系法从理论上讲，整体坐标系法可以求任意单元的形函数，但计算过程太复杂只能求一维2节点线单元、二维3节点三角形单元和三维4节点四面体单元3种简单单元的形函数复杂的或二次以上的单元必须采用局部坐标系法求位移场 u 是形函数 Ni 的线性组合，因此形函数Ni同样具有插值多项式的特性单元模型构造整体坐标系法从理论上讲，整体坐标系法可以求任意2024/5/675/162单元刚度矩阵2节点线单元一维2节点线单元单元插值关系 单元几何关系 单元本构关系 N=N1 N2 De=E单元刚度矩阵2节点线单元一维2节点线单元单元插值关系 单元2024/5/676/162单元刚度矩阵2节点线单元单元刚度矩阵A为单元截面积；l为单元长度矩阵B单元刚度矩阵2节点线单元单元刚度矩阵A为单元截面积；l为单2024/5/677/162单元刚度矩阵三角形单元二维3角形单元单元插值关系 单元刚度矩阵三角形单元二维3角形单元单元插值关系 2024/5/678/162单元刚度矩阵三角形单元单元几何关系 单元刚度矩阵三角形单元单元几何关系 2024/5/679/162单元刚度矩阵三角形单元单元本构关系 平面应力问题单元刚度矩阵三角形单元单元本构关系 平面应力问题2024/5/680/162单元刚度矩阵三角形单元矩阵B单元刚度矩阵三角形单元矩阵B2024/5/681/162单元刚度矩阵三角形单元单元刚度矩阵h为单元厚度k为对称的6*6常数矩阵A为单元面积单元刚度矩阵三角形单元单元刚度矩阵h为单元厚度k为对称的62024/5/682/162作业求4节点四面体单元的单元刚度矩阵作业求4节点四面体单元的单元刚度矩阵2024/5/683/162单元模型构造整体坐标系法单元形函数的特性正规性：单元形函数之和等于1。正交性：形函数在本节点的值等于1，在其它节点的值等于0。例如:2节点线单元形函数单元模型构造整体坐标系法单元形函数的特性正规性：单元形函数2024/5/684/162单元模型等参单元等参单元 单元内任意一点的位移u与单元节点位移ue之间的关系为 一般单元坐标的插值关系也采用与位移插值关系相同的变换关系即单元内任意一点的坐标x与单元节点坐标xe之间的关系为 单元模型等参单元等参单元 单元内任意一点的位移u与单元节点2024/5/685/162单元模型等参单元等参单元凡是几何形状和位移场采用同阶同参数插值关系来描述的单元，称为等参单元 前面介绍的所有单元都属于等参单元 在描述单元的几何形状和位移场时，并不一定非采用同阶插值关系 单元模型等参单元等参单元凡是几何形状和位移场采用同阶同参数2024/5/686/162单元模型等参单元等参单元3节点三角形等参单元 单元模型等参单元等参单元3节点三角形等参单元 2024/5/687/162单元模型等参单元超参单元如果几何形状插值函数的阶数高于位移场插值函数的阶数，称为超参单元 次参单元如果几何形状插值函数的阶数低于位移场插值函数的阶数，称为次参单元 单元模型等参单元超参单元如果几何形状插值函数的阶数高于位移2024/5/688/162单元平衡方程组装过程 为什么要组装?消除内力组装的原则是什么?单元自由度与结构自由度对应单元平衡方程组装过程 为什么要组装?消除内力组装的原则是2024/5/689/162单元平衡方程组装过程 2 F 1 3 U3U4U2U1U5U6结构自由度向量U单元平衡方程组装过程 2 F 1 3 U3U2024/5/690/162单元平衡方程组装过程 2 1 U3U4U2U11u1u2u3u4 3 U611 U2U12u1u2U5u3u42单元平衡方程组装过程 2 1 U3U4U2U11u12024/5/691/162单元平衡方程组装过程2 1 U3U4U2U11u1u2u3u42单元平衡方程组装过程2 1 U3U4U2U11u1u2024/5/692/162单元平衡方程组装过程组装单元单元平衡方程组装过程组装单元2024/5/693/162单元平衡方程组装过程 3 U611 U2U12u1u2U5u3u4单元平衡方程组装过程 3 U611 U2U12u1u2024/5/694/162单元平衡方程组装过程再组装单元总体刚度方程 K 称为总体刚度矩阵 U 称为位移向量 F 称为载荷向量 单元平衡方程组装过程再组装单元总体刚度方程 K 称为总体刚2024/5/695/162总体刚度矩阵K的特性 对称性 奇异性 稀疏性 非零元素带状分布 总体刚度矩阵K的特性 对称性 奇异性 稀疏性 非零元素带状分2024/5/696/162约束处理过程 为什么要约束处理?总体平衡方程组是奇异的消除无限制的刚体运动 使总体平衡方程组存在唯一一组解约束处理过程 为什么要约束处理?总体平衡方程组是奇异的消除2024/5/697/162约束处理过程边界条件 边界条件分类 力(载荷)边界条件位移边界条件 集中载荷力 表面分布力 自重力热交换引起的温度载荷 固定位移约束 强制位移约束 关联位移约束 约束处理过程边界条件 边界条件分类 力(载荷)边界条件位移2024/5/698/162约束处理过程模型简化xy约束处理过程模型简化xy2024/5/699/162约束处理过程模型简化yx xy约束处理过程模型简化yxxy2024/5/6100/162约束处理过程约束方程123456789101112yx约束处理过程约束方程123456789101112yx2024/5/6101/162约束处理过程约束处理方法位移约束处理方法 赋0赋1法 乘大数法 约束处理过程约束处理方法位移约束处理方法 赋0赋1法 乘大2024/5/6102/162约束处理过程赋0赋1法强制位移约束条件处理 U4=C约束处理过程赋0赋1法强制位移约束条件处理 U4=C2024/5/6103/162约束处理过程赋0赋1法强制位移约束条件处理 U4=C约束处理过程赋0赋1法强制位移约束条件处理 U4=C2024/5/6104/162约束处理过程赋0赋1法有6个方程，5个未知数，如果约束方程可以消除有限元平衡方程组的奇异性，则取任意5个方程联立求解，都会得到方程组的唯一一组解。系数矩阵由原来的对称的变成了非对称的，这对于大规模有限元方程组求解是十分不利的，采用相同的求解方法，在求解时间和矩阵存贮容量方面都增加了一倍。约束处理过程赋0赋1法有6个方程，5个未知数，如果约束方程2024/5/6105/162约束处理过程赋0赋1法为了保证系数矩阵的对称性，去掉方程组第4行约束处理过程赋0赋1法为了保证系数矩阵的对称性，去掉方程组2024/5/6106/162约束处理过程赋0赋1法引入强制位移约束方程 U4=C，使方程组求解时直接将自由度U4求出约束处理过程赋0赋1法引入强制位移约束方程 U4=C，使2024/5/6107/162约束处理过程赋0赋1法固定位移约束条件处理 U4=0约束处理过程赋0赋1法固定位移约束条件处理 U4=02024/5/6108/162约束处理过程赋0赋1法基本原理 利用初等变换对求解方程组进行相同的行列变换，既保证方程组解不会改变，又可以保持方程组系数矩阵的对称性。在进行初等变换时，只要保证对方程组系数矩阵做相同的行列变换，就可以保持方程组系数矩阵的对称性。约束处理过程赋0赋1法基本原理2024/5/6109/162约束处理过程乘大数法乘大数法基本原理 利用矩阵的初等变换不改变方程组解的思想。约束处理过程乘大数法乘大数法基本原理2024/5/6110/162约束处理过程乘大数法强制位移边界条件 约束处理过程乘大数法强制位移边界条件 2024/5/6111/162约束处理过程乘大数法强制约束方程 A是一个大数，是系数矩阵中对角线元素K44的1010倍量级以上 为什么要乘以大数A？放大位移约束方程的优势约束处理过程乘大数法强制约束方程 A是一个大数，是系数矩阵2024/5/6112/162约束处理过程乘大数法强制位移边界条件 约束处理过程乘大数法强制位移边界条件 2024/5/6113/162约束处理过程乘大数法固定位移边界条件 C=0 约束后的方程组简化为 约束处理过程乘大数法固定位移边界条件 C=0 约束后的2024/5/6114/162约束处理过程两种方法比较赋0赋1法在约束处理过程中是严格精确的，而乘大数法是一种近似约束处理方法，它的精度取决于所乘大数A值 两种方法都可以消除有限元平衡方程的奇异性，得到符合实际边界条件的唯一一组解。但两种方法还是有很大的区别 约束处理过程两种方法比较赋0赋1法在约束处理过程中是严格精2024/5/6115/162约束处理过程两种方法比较采用乘大数法约束处理后的有限元平衡方程在求解时可能造成解的失真，大数A值越大可能解的偏差会越大，而赋0赋1法就不会出现类似的问题，它在约束过程和求解过程都是精确的乘大数法相对于赋0赋1法在约束处理过程上简单一些 约束处理过程两种方法比较采用乘大数法约束处理后的有限元平衡2024/5/6116/162约束处理过程两种方法比较赋0赋1法实际上是将关联位移约束方程代入到有限元平衡方程中的，是代入法。而乘大数是将占绝对优势的关联位移约束方程合并到有限元平衡方程中的，是罚方法，计算误差来自于合并过程，计算精度取决于关联位移约束方程的优势大小商业软件中，位移边界条件的约束处理都采用赋0赋1法，乘大数很少被采用主要原因是它是一种近似方法，而且大数的大小也不好确定，有时还会造成求解失败 约束处理过程两种方法比较赋0赋1法实际上是将关联位移约束方2024/5/6117/162约束处理过程弹簧单元假设柔性弹簧 kOXYU4f f=kU4k约束处理过程弹簧单元假设柔性弹簧 kOXYU4f f=2024/5/6118/162弹簧约束方程 f=kU4约束处理过程弹簧单元弹簧约束方程 f=kU4约束处理过程弹簧单元2024/5/6119/162方程组求解过程特点方程组求解有限元计算过程中很重要的一部分，在有限元法的发展过程中，有限元方程的求解效率一直是其应用的最大瓶颈之一 有限元方程组的特点：有限元方程组的系数矩阵具有对称、稀疏、带状分布以及正定、主元占优。有效地利用这些特点，以减少系数矩阵的存贮量，提高方程组求解效率 方程组求解过程特点方程组求解有限元计算过程中很重要的一部分2024/5/6120/162方程组求解过程分类比较线性方程组的解法主要分两大类:直接解法：以高斯消去法基础，以等带宽或变带宽方式存贮系数矩阵内元素，对于求解规模比较大的问题，要存贮的元素非常巨大。迭代解法：只需要存贮系数矩阵中非零元素，存贮量很小，一般是变带宽存贮量的20%或更少，有些算法的求解效率也非常高，适合求解大规模线性方程组。但是这种解法对接近病态的方程组很难保证收敛性。方程组求解过程分类比较线性方程组的解法主要分两大类:2024/5/6121/162方程组求解过程带宽定义有限元方程组系数矩阵是稀疏的、非零元素呈带状分布，带宽就是它的宽度，带宽的大小是由系统有限元网格的节点号排序决定的，具体求法是 带宽=(单元最大节点号之差+1)*节点自由度数 带宽是网格节点标注方法直接决定的，不同标注方法带宽可能相关很大 方程组求解过程带宽定义有限元方程组系数矩阵是稀疏的、非零元2024/5/6122/162方程组求解过程带宽带宽是网格节点标注方法直接决定的，不同标注方法带宽可能相关很大 方程组求解过程带宽带宽是网格节点标注方法直接决定的，不同标2024/5/6123/162方程组求解过程带宽所示四边形网格的三种节点号标注方法，每个节点是2个自由度结构的带宽分别是12，18，56，相差很大，其中12和56之间相差近5倍，这就意味着系数矩阵的存贮量也是相差5倍，因此，对于大规模复杂系统的节点号优化是十分必要的 方程组求解过程带宽所示四边形网格的三种节点号标注方法，每个2024/5/6124/162方程组求解过程系数矩阵存贮 系数矩阵存贮 如果节点号排序优化的比较好，系数矩阵的存贮量就会减少很多。根据系数矩阵的对称性，一般都是按半带宽存贮。系数矩阵存贮的方法 二维等带宽存贮 一维变带宽存贮方程组求解过程系数矩阵存贮 系数矩阵存贮2024/5/6125/162方程组求解过程二维等带宽存贮 二维等带宽存贮 方程组求解过程二维等带宽存贮 二维等带宽存贮 2024/5/6126/162方程组求解过程二维等带宽存贮二维等带宽存贮消除了最大带宽以外的全部零元素，节省了系数矩阵元素的存贮量。但是由于取最大带宽为存贮范围，因此不能排除在带宽内的大量零元素。当系数矩阵的各行带宽变化不大时，适合采用二维等带宽存贮，方程组求解过程中系数矩阵元素的寻址也比较方便，求解效率较高。当出现局部带宽特别大的情况时，采用二维等带宽存贮时，将由于局部带宽过大而使整体系数矩阵的存贮大大增加。方程组求解过程二维等带宽存贮二维等带宽存贮消除了最大带宽以2024/5/6127/162方程组求解过程一维变带宽存贮 一维变带宽存贮 一维变带宽存贮方法就是把变化的带宽内的元素按一定的顺序存贮在一个一维数组中。由于它不按最大带宽存贮，因此比二维等带宽存贮更节省内存。按照解法可分为按行一维变带宽存贮和按列一维变带宽存贮。方程组求解过程一维变带宽存贮 一维变带宽存贮 2024/5/6128/162按行一维变带宽存贮 方程组求解过程一维变带宽存贮 辅助的寻址数组M 按行一维变带宽存贮 方程组求解过程一维变带宽存贮 辅助的寻2024/5/6129/162一维变带宽存贮是最节省内存的一种方法，但是由于要借助于寻址数组寻找系数矩阵元素的位置，相对二维等带宽存贮方法来说要复杂一些，而且在程序实现时也要复杂得多，方程组求解过程中也要消耗一些数组寻址时间。因此，在选用存贮方法时要权衡二者的利弊，统盘考虑。一般当带宽变化不大，计算机内存允许时，采用二维等带宽存贮方法是比较合适的。方程组求解过程一维变带宽存贮 一维变带宽存贮是最节省内存的一种方法，但是由于要借助于寻址数2024/5/6130/162方程组求解过程求解方法方程组求解方法 高斯消去法 三角分解法 雅可比（Jacobi）迭代法 高斯-赛德尔（Gauss-Seidel）迭代法 方程组求解过程求解方法方程组求解方法2024/5/6131/162应变、应力回代过程 单元应变和应力回代求解 通过求解有限元平衡方程得到有限元节点位移后，就可以进行系统的刚度校核。如果所分析问题要进行强度校核，就要回代求解单元的应变和应力。由插值关系和几何关系可得单元应变，再通过本构关系得到单元应力应变、应力回代过程 单元应变和应力回代求解 2024/5/6132/162有限差分法有限差分法有限差分法2024/5/6133/162 从弹性力学的基本方程建立以来，这些方程在各种问题的边界条件下如何求解，一直是很多数学工作者和力学工作者研究的内容。即弹性力学的经典解法存在一定的局限性，当弹性体的边界条件和受载情况复杂一点，往往无法求得偏微分方程的边值问题的解析解，许多工程重要问题，不能够得出函数式的解答。因此，弹性力学问题的各种数值解法便具有重要的实际意义。从弹性力学的基本方程建立以来，这些方程在各种问2024/5/6134/162 工程中常用的数值解法有有限单元法和差分法。有限单元法 是以有限个单元的集合体来代替连续体，属于物理上的近似。差分法 是把弹性力学的基本方程和边界条件（一般均为微分方程）近似地改用差分方程（代数方程）来表示，把求解微分方程的问题改换成为求解代数方程的问题，属于数学上的近似。工程中常用的数值解法有有限单元法和差分法。有2024/5/6135/162 第一节 差分方程 第二节 应力函数的差分解 第三节 深梁应力函数的差分解 第一节 差分方程 第二节 应力函数的差分解 第三节2024/5/6136/162第一节 差分方程 差分法是沿用已久的一种数值解法。随着计算机的普及和相应的软件发展，此法成为解弹性力学问题的一种有效的方法。第一节 差分方程 差分法是沿用已久的一种数值解法。2024/5/6137/162 我们在弹性体上，用相隔等间距h而平行于坐标轴的两组平行线织成正方形网格，x=y=h，如图。设f=f(x,y)为弹性体内的某一个连续函数。该函数在平行于x轴的一根网线上，如在-上,它只随x坐标的改变而变化。在邻近结点处，函数f可展为泰勒级数如下：我们在弹性体上，用相隔等间距h而平行于坐标轴的两组平2024/5/6138/162 我们将只考虑离开结点充分近的那些结点，即（x-x0）充分小。于是可不计（x-x0）的三次及更高次幂的各项，则上式简写为：在结点，x=x0-h,在结点1,x=x0+h，代入(b)得：我们将只考虑离开结点充分近的那些结点，即（x-x0联立(c),(d),解得差分公式：同理，在网线-上可得到差分公式联立(c),(d),解得差分公式：同理，在网线-上 差分公式(-)及(-)是以相隔2h的两结点处的函数值来表示中间结点处的一阶导数值，可称为中点导数公式。以相邻三结点处的函数值来表示一个端点处的一阶导数值，可称为端点导数公式。应当指出：中点导数公式与端点导数公式相比，精度较高。因为前者反映了结点两边的函数变化，而后者却只反映了结点一边的函数变化。因此，我们总是尽可能应用前者，而只有在无法应用前者时才不得不应用后者。差分公式(-)及(-)是以相隔2h的 以上(-)(-)是基本差分公式,从而可导出其它的差分公式如下：以上(-)(-)是基本差分公式2024/5/6142/162第二节 应力函数的差分解 当不计体力时，我们已把弹性力学平面问题归结为在给定边界条件下求解双调和方程的问题。用差分法解平面问题，就应先将双调和方程变换为差分方程，而后求解之。第二节 应力函数的差分解 当不计体力时，我 一旦求得弹性体全部节点的值后，就可按应力分量差分公式（对节点0）算得弹性体各节点的应力。一旦求得弹性体全部节点的值后，就可按应力分 可见，用差分法解平面问题，共有两大任务：一、建立差分方程 将（1-68）代入双调和方程 对于弹性体边界以内的每一结点，都可以建立这样一个差分方程。整理即得 可见，用差分法解平面问题，共有两大任务：2024/5/6145/162 二、联立求解这些线性代数方程，就能求得各内结点处的值。为了求得边界上各结点处的值，须要应用应力边界条件，即：一般建立和求解差分方程，在数学上不会遇到很大困难。但是，当对于边界内一行的（距边界为h的）结点，建立的差分方程还将涉及边界上各结点处的值，并包含边界外一行的虚结点处的值。二、联立求解这些线性代数方程，就能求得各内结代入上式，即得：l1=cos(N,x)=cos=dy/ds,l2=cos(N,y)=sin=-dx/ds,于是，式（a）可改写为：由右图可见，代入上式，即得：l1=cos(N,x)=cos=dy/d关于边界上任一点处由此得：的值，可将（b）式从A点到B点对s积分得到：关于边界上任一点处由此得：的值，可将（b）式从A点到B点对 将此式亦从A点到B点沿s进行积分，就得到边界上任一点B处的值。为此利用分部积分法，得：由高等数学可知，将此式亦从A点到B点沿s进行积分，就得到边界将式(b),(c)代入，整理得：由前知，把应力函数加上一个线性函数，并不影响应力。因此，可设想把应力函数加上a+bx+cy，然后调整a，b，c三个数值，使得由式(d)及式(c)可见，设即可根据面力分量及导数求得为已知，将式(b),(c)代入，整理得：由前知，从图易看出，式（23）右边的积分式表示A与B之间的x方向的面力之和；式（24）右边的积分式表示A与B之间的y方向的面力之和；式（25）右边的积分式表示A与B之间的面力对于B点的矩。于是式(d)，式(c)即简化为：从图易看出，式（23）右边的积分式表示A与 至此，我们解决了怎样计算边界上各结点的值的问题。至于边界外一行虚结点处 的值，则可用边界上结点处的 或 值和边界内一行相应结点处 的值来表示。例如，对于图1中的虚结点14，因为有 至此，我们解决了怎样的值的问题。至2024/5/6152/162所以有 当求出全部结点上的值以后，我们就可按应力分量的差分公式（21）计算应力分量。用差分法解弹性平面问题时，可按下列步骤进行：（1）在边界上任意选定一个结点作为基点A，取所以有 当求出全部结点上的值以后，我们就可 然后由面力的矩及面力之和算出边界上所有各结点处 的值，以及所必需的一些 及 值，即垂直于边界方向的导数值。（2）应用公式（26），将边界外一行虚结点处的值用边界内的相应结点处的值来表示。（3）对边界内的各结点建立差分方程（22），联立求解这些结点处的值。（4）按照公式（26），算出边界外一行的各虚结点处的值。然后由面力的矩及面力之和算出边界上所有各结点处2024/5/6154/162（5）按照公式（21）计算应力的分量。说明：1以上是针对单连体导出的结果。对于多连体，情况就不象这样简单。2.如果一部分边界是曲线的，或是不与坐标轴正交，则边界附近将出现不规则的内结点。对于这样的结点，差分方程（22）必须加以修正。（5）按照公式（21）计算应力的分量。说明：2024/5/6155/162第三节 例 深梁的应力函数差分解 现以如图所示的混凝土深梁为例，应用应力函数的差分解求出应力分量。已知混凝土深梁上边受有均布向下的铅直荷载q，并由下角点处的反力维持平衡。第三节 例 深梁的应力函数差分解 现以如图所示的混凝土(1）计算边界上各结点的、和 值。取A为基点，且由上面公式所得的计算结果见下表。(1）计算边界上各结点的、和 值。由上面公式所得（2）计算边界外以行各虚结点处的值。由式（2-6）及前表可得（2）计算边界外以行各虚结点处的值。(3）边界内各结点的差分方程,由式（2-2）可知(3）边界内各结点的差分方程,由式（2-2）可知联立求解上式，可得（以qh2为单位)（4）计算结点外一行各结点处的值。由(a)、(b)、(c)可得（5）计算应力。对于结点M，由式（1）可知 联立求解上式，可得（以qh2为单位)（4）计算结点外2024/5/6160/162同理可得 沿着梁的中线MA，的变化如下图和右图所示。同理可得 沿着梁的中线MA，的变化如下图和右图所示。差分方法的优缺点优点：1.差分法是解微分方程边值问题和弹性力学问题的有效方法；2.差分法简便易行；3.对于某些结构，为了更精确地分析局部的应力状态，可以用差分法进行分析；缺点：1.对于曲线边界和斜边界等产生的不等间距网格的处理，比较麻烦和易于出错；2.比较比较适用于求解二维问题或平面问题；3.比较适用于等间距网格，对于应力变化较为剧烈时，需要采用二次网格进行计算。差分方法的优缺点优点：2024/5/6162/162课堂作业课堂作业用差分法计算下图中用差分法计算下图中A A和和B B点的应力分量点的应力分量。FFBAaa234165XY7课堂作业用差分法计算下图中A和B点的应力分量。FFBAaa2
展开阅读全文