材料成形过程有限元法和有限差分法

1、1 材料成形过程 数值模拟 2 课程简介 n课时数：40课时 n学分数：2.5 n先修课程先修课程 n线性代数、工程力学、传质传热学、材料科学基础、材料 成型原理、热处理原理及工艺、铸造工艺及设备、焊接工 艺及设备等 3 n目的： n系统掌握材料成形过程数值模拟的基本原理和一般步骤 n了解材料成形时各加工方式的特点、基本规律、数值分析 原理与技术要点 n培养学习者分析解决实际工程问题的能力。 n要求： n掌握材料成形数值模拟的基本知识和理论 n掌握各种材料成形过程的相关理论、数值方法、实现过程 和实际应用 n了解材料加工过程中建模与数值分析技术现状和发展趋势 课程目的和要求 4 第一章 绪论

2、n基本概念 nCAE n将一个成形过程（或过程的某一方面）定义为由一组控制将一个成形过程（或过程的某一方面）定义为由一组控制 方程加上边界条件构成的定解问题，从而获得对成形过程方程加上边界条件构成的定解问题，从而获得对成形过程 的定量认识。的定量认识。 n材料成形包括: n金属金属铸造成形、塑性成形塑性成形、连接成形 n高分子材料黏流态注射成形 5 1.2 材料成形数值模拟的工程意义及应用现状 n工程意义工程意义 n制定和优化材料成形方案与模具设计方案 n解决工模具调试或产品试成形过程中的技术问题 n解决成形制品批量生产中的质量控制问题 n应用现状应用现状 n铸造成形、塑性成形、焊接成形、黏流

3、态成形 6 1.3 材料成形数值模拟的发展趋势 n（1）模拟分析由宏观进入微观 n（2）加大多物理场的耦合分析 n（3）拓宽数值模拟在特种成形中的应用 n（4）强化基础性研究 n（5）关注反向模拟技术应用 7 第二章 有限元与有限差分法基础 nCAE的工具： n有限元法（FEM）、有限差分法（FDM）、边界元法 （BEM）、有限体积法（FVM）、无网格法等等 n在材料成形的CAE中主要使用的是有限元法和有限差分法 8 “ ” 的早在20世纪40年代初期就 有人提出，但真正用于工程中则是电子计算机出现以后。 “ ” 这一名称是1960年美国的克拉夫 （Clough，R.W.）在一篇题为 “平面应

4、力分析的有限元 法” 论文中首先使用。此后，有限元法的应用得到蓬勃 发展。 到20世纪80年代初期国际上较大型的结构分析有限 元通用程序多达几百种，从而为工程应用提供了方便条件。 由于有限元通用程序使用方便，计算精度高，其计算结果 已成为各类工业产品设计和性能分析的可靠依据。 9 最初用于飞机结构的，由于它在 理论上的通用性，因而它可用于解决工程中的许多问题。 目前，它可以解决几乎所有的连续介质和场的问题， 包括热传导、电磁场、流体动力学、地质力学、原子工程 和生物医学等方面的问题。 中，从齿轮、轴、轴承等通用零部件到机 床、汽车、飞机等复杂结构的应力和变形分析（包括热应 力和热变形分析）。

5、不仅可以解决工程中的线性问题、非线性问 题，而且对于各种不同性质的固体材料，如各向同性和各 向异性材料，粘弹性和粘塑性材料以及流体均能求解； 对于工程中最有普遍意义的非稳态问题也能求解。 10 2.1 有限元法基础 基本思想：基本思想： n将一个连续求解域（对象）离散（剖分）成有限个形状 简单的子域（单元） n利用有限个节点将各子域连接起来 n在给定的初始条件和边界条件下进行综合计算求解，从 而获得对复杂工程问题的近似数值解 11 为什么要离散？为什么要离散？ n1.无法得到复杂实际问题的解析解 n2.将域划分成一些微小而形状规则的单元后，便于在一个 单元内得到近似解 n3.域中所有单元的解可

6、视为该复杂问题的近似解 12 有限元分析的过程 n 1.连续体离散化 n 2.单元分析 n 3.整体分析 n 4.确定约束条件 n 5.方程求解 n 6.结果分析与讨论 13 1.连续体离散化 n 连续体：是指所求解的对象（如物体或结构）。 n离散化（划分网格或网络化）：是将所求解的对象划 分为有限 n个具有规则形状的微小块体，把每个微小块体称为单元， 相邻两个 n单元之间只通过若干点互相连接，每个连接点称为节点。 n相邻单元只在节点处连接，载荷也只通过节点在各单 元之间传 n递，这些有限个单元的集合体，即原来的连续体。 n *单元划分后，给每个单元及节点进行编号； n *选定坐标系，计算各个

7、节点坐标； n *确定各个单元的形态和性态参数以及边界条件等。 14 基本上是任意的，一个结构体可以有 。但应遵循以下划分原则： (1) 分析清楚所讨论对象的性质，例如，是桁架结构 还是结构物，是平面问题还是空间问题等等。 (2) 单元的几何形状取决于结构特点和受力情况，单 元的几何尺寸(大小)要按照要求确定。一般来说，单元几 何形体各边的长度比不能相差太大。 (3) 有限元模型的网格划分越密，其计算结果越精确， 但计算工作量就越大。因此，在保证计算精度的前提下， 单元网格数量应尽量少。 (4) 在进行网格疏密布局时，应力集中或变形较大的 部位，单元网格应取小一些，网格应划分得密一些，而其 他

8、部分则可疏一些。 15 (5) 在设计对象的厚度或者弹性系数有突变的情况下， 应该取相应的突变线作为网格的边界线； (6) 相邻单元的边界必须相容，不能从一单元的边或者 面的内部产生另一个单元的顶点。 (7) 网格划分后，要将全部单元和节点按顺序编号，不 允许有错漏或者重复。 (8) 划分的单元集合成整体后，应精确逼近原设计对象。 原设计对象的各个顶点都应该取成单元的顶点。 所有网格的表面顶点都应该在原设计对象的表面上。所 有原设计对象的边和面都应被单元的边和面所逼近。 16 图例 n将悬臂梁划分为许多三角形单元 n三角形单元的三个顶点都是节点 n载荷直接施加在节点上 悬臂梁及其有限元模型 1

9、7 2.单元分析 连续体离散化后，即可对单元体进行特性分析，简称 为单元分析。 单元分析工作主要有两项： (1)选择单元位移模式(位移函数) 用节点位移来表示单元体内任一点的位移、应变和 应力，就需 搞清各单元中的位移分布。 一般是假定单元位移是坐标的某种简单函数，用其 模拟内位移的分布规律，这种函数就称为位移模式或位移 函数。通常采用的函数形式多为多项式。 根据所选定的位移模式，就可以导出用节点位移来 表示单元体内任一点位移的关系式。 18 2.单元分析(2) n(2) 分析单元的特性，建立单元刚度矩阵 n 进行单元力学特性分析，将作用在单元上的所有力 （表面 n力、体积力、集中力）等效地移

10、置为节点载荷； n 采用有关的力学原理建立单元的平衡方程，求得单 元内节 n点位移与节点力之间的关系矩阵单元刚度矩阵。 19 3. 整体分析 n n把各个单元的刚度矩阵集成为总体刚度矩阵，以及将 各单元的节点力向量集成总的力向量，求得整体平衡 方程。 n集成过程所依据的原理是节点变形协调条件和平衡条 件。 20 4. 确定约束条件 n由上述所形成的是一组线性代数方程， 在求解之前，必修根据具体情况分析，确定 的边界约束条件，并对这些方程进行适当修正。 21 5. 有限元方程求解 n应用求解机械结构应力类问题时，根据未知 量和分析 n有三种基本解法： 位移法 力法 混合法 22 (1)位移法 以

11、节点位移作为基本未知量，通过选择适当的位移函 数，进行单元的力学特性分析。在节点处建立单元刚度方 程，再组合成整体刚度矩阵，求解出节点位移后，进而由 节点位移求解出应力。 位移法优点是比较简单，规律性强，易于编写计算机 程序。所以得到广泛应用，其缺点是精度稍低。 (2)力法 以节点力作为基本未知量，在节点处建立位移连续方 程，求解出节点力后，再求解节点位移和单元应力。 力法的特点是计算精度高。 (3)混合法 取一部分节点位移和一部分节点力作为基本未知量， 建立平衡方程进行求解。 23 单元特性的推导方法 的推导是的基本步骤之 一。目前，建立单元刚度矩阵的方法主要有以下四种： 直接刚度法 虚功原

12、理法 能量变分法 加权残数法 24 1. 直接刚度法直接刚度法 是直接应用物理概念来建立单元的有限 元方程和分析单元特性的一种方法。这一方法仅能适用 于简单形状的单元，如梁单元。但它可以帮助理解有限 元法的物理概念。 图1所示是xoy平面中的一，现以它为例， 来说明用建立单元刚度矩阵的思想和过程。 图1平面简支梁元及其计算模型 25 梁在横向外载荷（可以是集中力或分布力或力矩等）作用下产 生弯曲变形，在水平载荷作用下产生线位移。 对于该平面简支梁问题： 梁上任一点受有三个力的作用： 水平力Fx， 剪切力Fy ， 弯矩Mz。 相应的位移为： 水平线位移u， 挠度v ， 转角 z 。 由上图可见：

13、 水平线位移和水平力向右为正， 挠度和剪切力向上为正， 转角和弯矩逆时针方向为正。 通常规定： 26 为使问题简化，可把图示的梁看作是一个梁单元。 如图1所示，当令左支承点为节点 i ，右支承点为节点 j 时， 则该单元的节点位移和节点力可以分别表示为： 称为单元的节点位移列阵。 称为；若 F 为，则称为载荷列阵。 （1-1） （1-2） 写成矩阵形式为 q(e)=ui ,vi , zi ,vj ,uj , zjT ui ,vi , zi ,vj ,uj , zj F(e)=Fxi ,Fyi ,Mzi ,Fxj ,Fyj ,MzjT Fxi ,Fyi ,Mzi ,Fxj ,Fyj ,Mzj 2

14、7 显然，梁的节点力和节点位移是有联系的。在弹性 小变形范围 内，这种关系是线性的，可用下式表示 111213141516 212223242526 313233343536 414243444546 515253545556 6162636 xi yi zi xj yj zj F kkkkkk F kkkkkk Mkkkkkk Fkkkkkk Fkkkkkk kkkkM 46566 i i zi j j zj u v u v kk ( )( )( )eee FKq或 （1-3b） （1-3a） 28 上式(1-3b)称为，或称为单元刚度方 程，它代表了单元的载荷与位移之间（或力与变形之间）

15、的联系； 式中，K(e)称为单元刚度矩阵，它是单元的特性矩阵。 对于图1所示的平面梁单元问题，利用材料力学中的杆 件受力与变形间的关系及叠加原理，可以直接计算出单元 刚度矩阵K(e)中的各系数 kst( s, t = i, j ) 的数值 29 2. 虚功原理法虚功原理法 下面以平面问题中的三角形单元为例，说明利用虚功原理法来建 立单元刚度矩阵的步骤。 如前所述，将一个连续的弹性体分割为一定形状和数量的单元， 从而使连续体转换为有限个单元组成的组合体。单元与单元之间仅通 过节点连结，除此之外再无其他连结。也就是说，一个单元上的只能 通过节点传递到相邻单元。 从分析对象的组合体中任取 一个三角形

16、单元： 设其编号为 e ， 三个节点的编号为i、j、m， 在定义的坐标系 xoy 中，节点坐 标分别为(x j , y j)、(xi , y i)、(xm, ym)，如图2所示。 图2三节点三角形单元 30 由弹性力学平面问题的特点可知，有两个位移分 量，即每个单元有6个自由度，相应有6个节点载荷，写成矩阵形式， 即 单元节点载荷矩阵： F(e)=Fxi ,Fyi ,Fxj ,Fyj ,Fxm ,FymT 单元节点位移矩阵： q(e)=ui ,vi ,uj ,vj ,um ,vmT 图2三节点三角形单元 31 (1)设定位移函数 按照有限元法的基本思想：首先需设定一种函数来近似表达单元 内部的

17、实际位移分布，称为位移函数，或位移模式。 三节点三角形单元有6个自由度，可以确定 6个待定系数，故三 角形单元的位移函数为 （1-4） 式(1-4)为线性多项式，称为线性位移函数，相应的单元称为线性单元。 u=u(x,y)= 1+ 2x+ 3y v=v(x,y)= 4+ 5x+ 6y 32 1 2 3 4 5 6 1 0 0 0 0 0 0 1 ux y ds vx y 上式(5-5)也可用矩阵形式表示，即 式中，d为单元内任意点的位移列阵。 （1-5） 33 由于节点 i、j、m 在单元上，它们的位移自然也就满足位移函 数式(1-4)。 设三个节点的位移值分别为( ui, vi)、( uj,

18、 vj )、( um, vm )，将节 点位移和节点坐标代入式(1-4)，得 123 123 123 iii jjj mmm uxy uxy uxy 456 456 456 iii jjj mmm vxy vxy vxy 34 ( ) 1 2 3 4 5 6 000 000 000 1 0002 000 000 e ijmi ijmi ijnj ijmj ijmm ijmm aaau bbbv cccu aaav bbbu cccv （1-6） 1 111 1 222 1 ii jjijjmmijimjim mm xy xyx yx yx yx yx yx y xy 式中 （1-7） 由上可知

19、，共有6个方程，可以求出6个待定系数。解方程，求得 各待定系数和节点位移之间的表达式为 为三角形单元的面积。其中： 35 , , , , , , ijmmjjmiimmijji ijmjmimij imjjimmji ax yx yax yx yax yx y byybyybyy cxxcxxcxx （1-8） 将式(1-7)及式(1-8)、式(1-9)代入式(1-6)中，得到 （1-9） （1-10） 36 式中，矩阵N 称为单元的形函数矩阵； 为单元节点位移列阵。 其中， 为单元的形函数，它们反映单元内位移的 分布形态，是x, y 坐标的连续函数，且有 ( ) e q , , ijm NN

20、N 2 2 2 iiii jjjj mmmm Nab xc y Nab xc y Nab xc y （1-11） 式(1-10)又可以写成 , , , , iijjmmii i i j m iijjmmii i i j m uN uN uN uN u vN vN vN vN v （1-12） 上式清楚地表示了单元内任意点位移可由节点位移插值求出。 37 (2) 利用几何方程由位移函数求应变 根据弹性力学的几何方程 ， 线应变 剪切应变 则应变列阵可以写成 /, / , xy uxuy / , xy uyux ( )( ) 000 1 000 2 i i xijm ee j yijm j xyi

21、ijjmm m m u u v x bbb u u cccB q vy cbcbcb uuu yxv 式中，B称为单元应变矩阵，它是仅与单元几何尺寸有关的常量矩阵， 即 （1-13） 38 000 1 000 2 ijm ijm iijjmm bbb Bccc cbcbcb （1-14） 上述方程(1-13)称为单元应变方程，它的意义在于： 单元内任意点的应变分量亦可用基本未知量即节点位移分量来表 示。 39 (3)利用广义虎克定律求出单元应力方程 根据广义虎克定律，对于平面应力问题 1 () 1 () 12(1) xxy yyx xyxyxy E E GE 上式（1-15）也可写成 ( )(

22、 )ee D （1-15） （1-16） 式中， 为应力列阵； D 称为弹性力学平面问题的弹性矩阵，并有 , T xyxy 40 2 10 10 1 1 00 2 E D ( ) ( )( )( ) e x eee y xy DDBq 则有如下单元应力方程 由式(1-18)可求单元内任意点的应力分量，它也可用基本未知量即 节点位移分量来表示。 （1-17） （1-18） 41 (4)由虚功原理求单元刚度矩阵 根据虚功原理，当弹性结构受到外载荷作用处于平衡状态时， 在任意给出的微小的虚位移上，外力在虚位移上所做的虚功 AF等 于结构内应力在虚应变上所存储的虚变形势能 A ，即 F AA 设处于平

23、衡状态的弹性结构内任一单元发生一个微小的虚位移， 则单元各节点的虚位移 为 ( ) e q ( ) , , , , , eT iijjmm quvuvuv （1-20） （1-21） （1-19） 则单元内部必定产生相应的虚应变，故单元内任一点的虚应变 为 ( ) e ( ) , , eT xyxy 42 显然，虚应变和虚位移之间关系为 ( )( ) ee Bq 设节点力为 ( ) , , , , , T e xiyixjyjxmym FFFFFFF 则外力虚功为 ( ) ( ) T e e F AqF （1-24） （1-22） （1-23） 单元内的虚变形势能为 ( ) T e V Adv

24、 43 根据虚功原理 ( )( ) ( ) TT ee e V qFdv ( )( )ee DDBq ( )( )( ) TTT eee T BqBq 因为 （1-26） （1-25） 代人式(1-25)，则有 ( )( ) ( )( ) TT ee Tee V qFBqDBqdv 式中， ，均与坐标 x, y 无关，故可以从积分符号中 提出，可得： ( ) T e q ( )e q 44 ( )( )( )( )Teeee V FBDB dvqKq 其中，单元刚度矩阵 ( )eT V KBDB dv （1-27） 式(1-27)称为单元有限元方程，或称单元刚度方程，其中 是单元刚度矩阵。 (

25、 )e K （1-28） 因为三角形单元是常应变单元，其应变矩阵B 、弹性矩阵D 均为常量，而 ，所以式(1-28)可以写成 V dvtdxdyt ( )eT KtBDB （1-29） 45 式中，t 为三角形单元的厚度； 为三角形单元的面积。 ( ) ( ) e iiijim e jijjjm mimjmm KKK KKKK KKK 对于图2所示的三角形单元，将D 及B代入式(1-28)，可以得 到单元刚度为 （1-30） 式中: K为66阶矩阵，其中每个子矩阵为22阶矩阵，由下式给出 2 11 22 114(1) c 22 ( , ,) rsrsrsrs rs rsrsrsrs b bc

26、cb cc b Et K c bb ccb b r si j m （1-31） 46 按照力学的一般说法，任何一个实际状态的是 这个系统从实际状态运动到某一参考状态(通常取弹性体外载荷为 零时状态为参考状态)时。弹性体的总位 能 是一个函数的函数，即泛函，位移是泛函的容许函数。 从考虑，变形弹性体受外力作用处于平衡状态时，在 很多可能的变形状态中，使总位能最小的就是弹性体的真正变形， 这就是最小位能原理。用求能量泛函的极值方法就是能量变 分原理。能量变分原理除了可解机械结构位移场问题以外，还扩展 到求解热传导、电磁场、流体力学等连续性问题。 3. 能量变分原理法 47 该方法是将假设的场变量的

27、函数(称为试函数)引入问题的控制方 程式及边界条件，利用最小二乘法等方法使残差最小，便得到近似的 场变量函数形式。 该方法的优点是不需要建立要解决问题的泛函式，所以，即使没 有泛函表达式也能解题。 4. 加权残数法 48 n有限元解的收敛性有限元解的收敛性 n有限元解是近似解 n近似解是否收敛于真实解、近似解收敛速度、近似解的稳定性 n近似解的收敛条件: n1.完备性准则（必要条件） n试探函数（插值函数）的次数（m）不小于场函数的最高可导阶数 n2.协调性准则（充分条件） n试探函数在m-1次连续可导。 49 有限元分析的误差 有限元 分析误差 建模误差 计算误差 离散误差 物理离散误差 几

28、何离散误差 边界条件误差 单元形状误差 舍入误差 截断误差 插值函数与真实函数之间的差异 1.减小单元特征尺寸，称为减小单元特征尺寸，称为h法法 2.提高插值函数的阶次，称为提高插值函数的阶次，称为p法法 单元组合体与求解对象几何形状的差异 1.网格局部加密网格局部加密 2.选用边或面上带有节点的单元选用边或面上带有节点的单元 边界条件的复杂性 1.准确测定，完善模型准确测定，完善模型 2.细分边界网格细分边界网格 单元严重畸变而退化 细分局部网格或者控制调整关键区域的网格细分局部网格或者控制调整关键区域的网格 数据储存 计算方法、解题性质、解题规模 注意网格的划分注意网格的划分 选择合适的解

29、算方法选择合适的解算方法 控制解题的规模控制解题的规模 减少运算次数，降低解题规模减少运算次数，降低解题规模 选择合适的解算方法，控制解题规模选择合适的解算方法，控制解题规模 50 材料成形中的非线性问题 n1.材料非线性 n材料本构方程非线性 弹塑性、刚弹性、刚黏塑性、黏弹塑性 n2.几何非线性 n3.边界非线性 51 2.2 有限差分法基础 n一种直接将微分问题转变成代数问题的近似数值解法。 n基本思想基本思想 n数值微分法数值微分法 n是把连续的定解区域划分成差分网格，用有限个节点代替 原连续求解域。 n把原方程和定解条件中的微商用差商来近似 n把原微分方程和定解条件近似地用代数方程组代

30、替，即有 限差分方程 52 差分网格通常为矩形 在边界不规则或者形状复杂时精度降低 有限元网格有限差分网格 53 差分概念 n自变量x的解析函数 y =f (x)，则有： ndx,dy自变量和函数的微分 n 函数对自变量的一阶导数 n 函数对自变量的一阶差商 00 ()( ) limlim xx dyyf xxf x dxxx y x dy dx 差商差商 54 差分方向 n向前差分 n向后差分 n中心差分 ()( )yf xxf x ( )()yf xf xx 11 ()() 22 yf xxf xx 55 差商的截断误差 n将函数f (x+x)按Taylor级数展开 n向前 n向后 n中心

31、 23 4 ()() ()( )( )( )( )0() ) 1!2!3! xxx f xxf xfxfxfxx 2 3 ()( )() ( )( )( )0() ) 2!3! f xxf xxx fxfxfxx x 2 3 ( )()() ( )( )( )0() ) 2!3! f xf xxxx fxfxfxx x 2 3 ()() () 22 ( )( )0() ) 2! 3! xx f xf x x fxfxx x 56 n二阶中心差商 n通常采用向前差商的向后差商 n截断误差与(x)2 同一数量级 n一阶向前差商 一阶向后差商 一阶计算精度 n一阶中心差商 二阶中心差商 二阶计算精度

32、 2 22 22 3 2 ()2 ( )() () 2() ( )( )0() ) 4! yf xxf xf xx xx yx fxfxx x 57 我们在弹性体上，用相隔等 间距h而平行于坐标轴的两组平行 线织成正方形网格，x=y=h， 如图。 设f=f(x,y)为弹性体内的某一 个连续函数。该函数在平行于x 轴的一根网线上，如在- 上,它只随x坐标的改变而变化。 在邻近结点处，函数f可展为 泰勒级数如下： .)( ! 3 1 )( ! 2 1 )( 3 0 0 3 3 2 0 0 2 2 0 0 0 xx x f xx x f xx x f ff 58 我们将只考虑离开结点充分近的那些结

33、点，即（x-x0）充分小。于是可不计（x-x0）的 三次及更高次幂的各项，则上式简写为： 在结点，x=x0-h, 在结点1, x=x0+h，代入(b) 得： )(.)( ! 2 1 )( 2 0 0 2 2 0 0 0 bxx x f xx x f ff )(. 2 0 2 22 0 03 c x fh x f hff )(. 2 0 2 22 0 01 d x fh x f hff 59 联立(c),(d),解得差分公式： )11( 2 31 0 h ff x f )21 ( 2 2 031 0 2 2 h fff x f 同理，在网线-上可得到差分公式 )41 ( 2 2 )31 ( 2

34、042 0 2 2 42 0 h fff y f h ff y f 60 从而可导出其它的差分公式如下： )()(46 1 )()(24 1 )()(46 1 )()( 4 1 1210420 4 0 4 4 876543210 4 0 22 4 119310 4 0 4 4 7586 0 2 2 fffff hy f fffffffff hyx f fffff hx f ffff h yx f 61 相隔2h的两结点处的函数值来表示中间结点处的一 阶导数值，可称为中点导数公式。 以相邻三结点处的函数值来表示一个端点处的一阶 导数值，可称为端点导数公式。 中点导数公式与端点导数公式相比，精度较

35、高。 因为前者反映了结点两边的函数变化，而后者却只反 映了结点一边的函数变化。 62 边界元法简介 n边界元法（boundary element method） n一种结合有限元法和边界积分法发展起来的一种新数值方 法 n只在定义域的边界上划分单元，用满足控制方程的函数去 逼近边界条件。 n适用于应力（薄板）、流体力学、声场、电磁场等问题 63 边界元法基本思想 n以微分控制方程的基本解为权函数，利用加权余量法加权余量法将区 域积分转化为边界积分，并结合求解域边界的离散，构建 基于边界单元的代数方程组，然后进行计算求解 n以定义在边界上的边界积分方程为控制方程，通过对边界元插值离散， 化为代数

36、方程组求解 n降低了问题的维数，从而显著降低了自由度数 n边界的离散比区域的离散方便得多，可用较简单的单元准确地模拟边 界形状，最终得到阶数较低的线性代数方程组 64 加权余量法简介 n一种直接从所需求解的微分方程及边界条件出发，寻求边值问题近似 解的数学方法。 n从静力发展到了动力、稳定、材料非线性和几何非线性等各方面。 n在求解域上建立一个试函数 试函数由完备函数集的子集构成。已被采用过的试函数有幂级数、三角级数、 样条函数、贝赛尔函数、切比雪夫和勒让德多项式等等。 n试函数与真实解之间的偏差，即余量（内部和边界） n引入权函数，定义消除余量的条件 n加权余量法就是一种定义近似解与真解之间

37、余量，并设法使其最小的 方法。 65 设问题的控制微分方程为: 在V域内 在S边界上 式中 : L、B分别为微分方程和边界条件中的微分算子； f、g 为与未知函数u无关的已知函数域值； u为问题待求的未知函数。 ( )0B ug ( )0L uf 66 当利用加权余量法求近似解时，首先在求解域上建立一个试函数 ， 一般具有如下形式： u 1 n ii i uC NNC 式中： 待定系数，也可称为广义坐标； i C 取自完备函数集的线性无关的基函数。 i N 由于 一 般只是待求函数u的近似解，因此代入后将得 不到满足，若记： u ( ) ( ) I B RL uf RB ug 在V域内 在S边

38、界上 显然 反映了试函数与真实解之间的偏差，它们分别称 做内部和边界余量余量。 BI RR、 67 若在域V内引入内部权函数 ，在边界S上引入边界权函数 则可建立n个消除余量的条件，一般可表示为： I W B W 0(1,2, ) IiIBiB VS W R dVW R dSin 不同的权函数 和 反映了不同的消除余量的准则。从上 式可以得到求解待定系数矩阵C的代数方程组。一经解得待定 系数，由式(5.1.3)即可得所需求解边值问题的近似解。 Bi W Ii W 68 Ii W 由于试函数 的不同，余量 和 可有如下三种情况， 依此加权余量法可分为： 1内部法 试函数满足边界条件，也即 此时消

39、除余量的条件成为: 2边界法 试函数满足控制方程，也即 此时消除余量的条件为: ( )0 B RB ug ( )0 I RL uf 0(1,2, ) IiI V W R dVin 0(1,2, ) BiB S W R dSin u I R B R 69 n3混合法 n 试函数不满足控制方程和边界条件，此时消除余量的条件为: n显然，混合法对于试函数的选取最方便，但在相同精度条件下，工作 量最大。对内部法和边界法必须使基函数事先满足一定条件，这对复 杂结构分析往往有一定困难，但试函数一经建立，其工作量较小。 0(1,2, ) IiIBiB VS W R dVW R dSin 70 边界元法特点

40、n 1 前处理工作量小 n 2 解算规模小 n 3 求解奇异性问题时计算精度高 n 4 在载荷集中和半无限域等问题上有优势 n相对于有限元法，边界元法发展较慢， 71 有限元法解决应力集中问题 n在应力分析中对于应力集中区域必须划分很多的 单元，从而增加了求解方程的阶数，计算费用也 就随之增加 n用位移型有限元法求解出的应力的精度低于位移 的精度，对于一个比较复杂的问题必须划分很多 单元，相应的数据输人量就很大，同时，在输出 的大量信息中，又有许多并不是人们所需要的。 72 n1.使问题的维数降低一维 n原为三维空间的可降为二维空间，原为二维空间的问题可原为三维空间的可降为二维空间，原为二维空

41、间的问题可 降为一维。降为一维。 n2.只需将边界离散而需将区域离散化 n所划分的单元数目远小于有限元，这样它减少了方程组的所划分的单元数目远小于有限元，这样它减少了方程组的 方程个数和求解问题所需的数据，不但减少了准备工作，方程个数和求解问题所需的数据，不但减少了准备工作， 而且节约了计算时间。而且节约了计算时间。 边界元法与有限元法比较 73 n3.由于直接建立在问题控制微分方程和边界条件上的，不 需要事先寻找任何泛函。有限元法是以变分问题为基础， 如果泛函不存在就难于使用。 n可以求解经典区域法无法求解的无限域类问题。 n4.由于边界元法引入基本解，具有解析与离散相结合的特 点，因而具有

42、较高的精度。 74 边界元法的缺点 n对于非均值和非线性问题求解比较困难 n用边界元法解非线性问题时，遇到同非线性项相对应的区 域积分，这种积分在奇异点附近有强烈的奇异性，使求解 遇到困难。 75 有限体积法 n基本思想 n将计算区域划分为一系列不重复的控制体积，并使每个网 格点周围有一个控制体积。将待解的微分方程对每一个控 制体积积分积分，便得出一组离散方程。其中的未知数是网格 点上的因变量的数值。 n n有限体积法得出的离散方程，要求因变量的积分守恒对任 意一组控制体积都得到满足，对整个计算区域，自然也得 到满足 n n属于加权余量法中的子域法属于加权余量法中的子域法 76 加权余量法子域

43、法(Subdomain Method) 此法首先将求解域V划分成n个子域 ，在每个子 域 内令权函数等于1，而在子域之外取权函数为零，也 即: i V 1() 0() i Ii i V W V 内 外 如果在各个子域里分别选取试函数，那么它的求解 在形式上将类似于有限元法。 77 有限体积法特点 n即使在粗网格情况下，也显示出准确的积分守恒 n就离散方法而言，有限体积法可视作有限元法和 有限差分法的中间物 n有限体积法在寻求控制体积的积分时，必须假定 值在网格点之间的分布；有限元法必须假定值在 网格点之间的变化规律（既插值函数），并将其 作为近似解。 n有限体积法只寻求结点值；有限差分法只考虑

44、网 格点上的数值而不考虑值在网格点之间如何变化。 78 2.4 应用数值方法模拟的若干注意事项 n1 简化模型 n简化细节（如倒角等等）；利用对称结构简化细节（如倒角等等）；利用对称结构 n2 选择单元 n节点尽量少；曲（线）面边界用带（边）面节点的单元；轴对称结构节点尽量少；曲（线）面边界用带（边）面节点的单元；轴对称结构 采用轴对称单元采用轴对称单元 n3 划分网格 n细分网格；网格自适应（细分网格；网格自适应（remeshing） n4 建立初始条件和边界条件 n5定义材料参数 79 有限元常用软件 n从20世纪70年代开始，基于有限元法在结构线性 分析方面已经成熟，并被工程界广泛采用，

45、一批 由专业软件公司研制的大型通用商业软件（如 NASTRAN，ASKA，SAP，ANSYS，MARC，ABAQUS， JIFEX等）公开发行和被应用。它包含众多的单元 型式、材料模型及分析功能，并具有网格自动划 分、结果分析和显示等先后处理功能。 80 编号编号软件名称软件名称开发单位开发单位编程语言编程语言( (程序规模程序规模) ) 1ANSYSSwanson Analysis System(美国美国)FORTRAN77(150000FORTRAN77(150000行行) ) 2ADINAADINA工程公司工程公司(美国美国)FORTRAN(150000FORTRAN(150000行行)

46、 ) 3MARCMARC公司公司(美国美国)FORTRAN4FOR66FOR77(100000FORTRAN4FOR66FOR77(100000行行) ) 4NASTRAN NASA(美美)主持，主持，MSC公司公司(美美)开开 发发 FORTRAN4Assembler(600000FORTRAN4Assembler(600000行行) ) 5ABAQUS Hilbitt, Karlson, and Sorensen公公 司司(美美) FORTRAN77(140000FORTRAN77(140000行行) ) 6FENRISNTH, SINTEF(挪威挪威)FORTRAN77(160000FO

47、RTRAN77(160000行行) ) 7PAFECPAFEC公司公司(美美)FORTRAN(400000FORTRAN(400000行行) ) 8ASKA 斯图加特大学静动力学研究所斯图加特大学静动力学研究所(西西 德德) FORTRAN4(600000FORTRAN4(600000行行) ) 9EALEISE公司公司(美美) ANSI FOR66Assembler(150000ANSI FOR66Assembler(150000 行行) ) 10SAMCEFL.T.A.S(比利时比利时)FORTRAN4(300000FORTRAN4(300000行行) ) 11 LARSTRAN8 0 斯

48、图加特大学静动力学研究所斯图加特大学静动力学研究所(西西 德德) FORTRAN4FORTRAN77(200000FORTRAN4FORTRAN77(200000行行) ) 12HAJIF系列系列航空部航空部(中国中国)FORTRAN4(280000FORTRAN4(280000行行) ) 81 几种有限元程序及其应用范围 程序名程序名 应用应用 范围范围 ADINAANSYSASKAMARCNASTRANSAP 非线性分析非线性分析 塑性分析塑性分析 断裂力学断裂力学 热应力与蠕变热应力与蠕变 厚板厚壳厚板厚壳 管道系统管道系统 船舶结构船舶结构 焊接接头焊接接头 粘弹性材料粘弹性材料 82

49、 几种有限元程序及其应用范围（续） 程序名程序名 应用应用 范围范围 ADINAANSYSASKAMARCNASTRANSAP 热传导热传导 薄板薄壳薄板薄壳 复合材料复合材料 结构稳定性结构稳定性 流体力学流体力学 气动弹性力学气动弹性力学 电场电场 83 ANSYS 简介 n1970 年成立的美国ANSYS公司是世界CAE行业 最著名的公司之一,长期以来一直致力于设计分析 软件的开发、研制,其先进的技术及高质量的产品 赢得了业界的广泛认可。 n在我国,ANSYS用户也越来越多,三峡工程、二滩 电站、黄河下游特大型公路斜拉桥、国家大剧院、 浦东国际机场等在结构设计时都采用了ANSYS作 为分

50、析工具。 84 ANSYS的分析方法 1. 建立有限元模型建立有限元模型 3. 查看结果查看结果 2. 施加载荷求解施加载荷求解 主菜单主菜单 分析的三个主要步骤可在主菜单中得到明确体现分析的三个主要步骤可在主菜单中得到明确体现. Objective 2-2. ANSYS分析步骤在分析步骤在GUI中的体现中的体现. 85 ANSYS的分析方法(续) ANSYS GUI中的功能排列中的功能排列 按照一种动宾结构，以动按照一种动宾结构，以动 词开始（如词开始（如Create), 随后随后 是一个名词是一个名词 (如如Circle). 菜单的排列，按照由前到后菜单的排列，按照由前到后 、由简单到复杂

51、的顺序，与、由简单到复杂的顺序，与 典型分析的顺序相同典型分析的顺序相同. 86 ANSYS文件及工作文件名 n一些特殊的文件 n数据库文件jobname.db二进制 nLog 文件jobname.log文本 n结果文件jobname.rxx二进制 n图形文件jobname.grph二进制 ANSYS的数据库，是指在前处理、求解及后处理过 程中，ANSYS保存在内存中的数据。数据库既存储输入的数 据，也存储结果数据: n输入数据 - 必须输入的信息 (模型尺寸、材料属性、载荷等). n结果数据 - ANSYS计算的数值 (位移、应力、应变、温度 等). 87 启动ANSYS Objective

52、 1-1. 启动启动ANSYS软件软件. 注意注意: 以下过程默认为在以下过程默认为在 Windows NT环境下环境下. 要启动要启动ANSYS:Windows NT 屏幕屏幕: Start Programs ANSYS X.x 1. 选择选择 Interactive. Procedure 1. . 2. . 3. . 88 启动ANSYS(续) 3. 选择选择ANSYS的工作目录，的工作目录，ANSYS所有生成的所有生成的 文件都将写在此目录下。缺省为上次运行定义文件都将写在此目录下。缺省为上次运行定义 的目录。的目录。 4. 如果配置了如果配置了3D显卡则选择显卡则选择3D. 5. 设定

53、初始工作文件名，缺省为上次运行定义的设定初始工作文件名，缺省为上次运行定义的 工作文件名，第一次运行缺省为工作文件名，第一次运行缺省为 file. 7. 选择选择 Run to 运行运行ANSYS. Note: 还可以设定其他的交互选项，但通常以上几项通常还可以设定其他的交互选项，但通常以上几项通常 要设置好要设置好. 如果在第如果在第1步中选择步中选择 Run Interactive Now，将读取，将读取 上一次的设置，跳过此窗口，直接运行上一次的设置，跳过此窗口，直接运行ANSYS, 2. 选择选择ANSYS产品产品. 6. 设定设定ANSYS工作空间及数据库大小工作空间及数据库大小 (

54、参考参考 ANSYS安装及配置手册安装及配置手册). 89 ANSYS窗口 Objective 1-2. ANSYS GUI中六个窗口的总体功能中六个窗口的总体功能 输入输入 显示提示信息，输入显示提示信息，输入 ANSYS命令，所有输入命令，所有输入 的命令将在此窗口显示的命令将在此窗口显示 。 主菜单主菜单 包含包含ANSYS的主要功能的主要功能 ，分为前处理、求解、，分为前处理、求解、 后处理等。后处理等。 输出输出 显示软件的文本输出。显示软件的文本输出。 通常在其他窗口后面，通常在其他窗口后面， 需要查看时可提到前面需要查看时可提到前面 。 应用菜单应用菜单 包含例如文件管理、选包含

55、例如文件管理、选 择、显示控制、参数设择、显示控制、参数设 置等功能置等功能. 工具条工具条 将常用的命令制成工具将常用的命令制成工具 条，方便调用条，方便调用. 图形图形 显示由显示由ANSYS创建或传创建或传 递到递到ANSYS的图形的图形. 90 一、主要功能简介 n1. 结构分析 n 1） 静力分析 - 用于静态载荷. 可以考虑结构的线性及 非线性行为。 n 线性结构静力分析 n 非线性结构静力分析 n 几何非线性：大变形、大应变、应力强化、旋转软 化 n 材料非线性：塑性、粘弹性、粘塑性、超弹性、多 线性弹性、蠕变、肿胀等 n 接触非线性：面面/点面/点点接触、柔体/柔体刚体 接触、

56、热接触 n 单元非线性：死/活单元、钢筋混凝土单元、非线性 阻尼/弹簧元、预紧力单元等 91 n 2）模态分析 - 计算线性结构的自振频率及振形. 谱 分析 是模态分析的扩展，用于计算由于随机振动引起 的结构应力和应变 (也叫作 响应谱或 PSD). n 3）谐响应分析 - 确定线性结构对随时间按正弦曲 线变化的载荷的响应. n 4）瞬态动力学分析 - 确定结构对随时间任意变化 的载荷的响应. 可以考虑与静力分析相同的结构非线性 行为. n 5）谱分析 n 6）随机振动分析等 n 7）特征屈曲分析 - 用于计算线性屈曲载荷并确定 屈曲模态形状. (结合瞬态动力学分析可以实现非线性屈 曲分析.)

57、 n 8）专项分析: 断裂分析, 复合材料分析，疲劳分析 92 n2. 高度非线性瞬态动力分析（ANSYS/LS-DYNA） n 全自动接触分析，四十多种接触类型 n 任意拉格郎日欧拉（ALE）分析 n 多物质欧拉、单物质欧拉 n 适应网格、网格重划分、重启动 n 100多种非线性材料模式 n 多物理场耦合分析：结构、热、流体、声学 n 爆炸模拟，起爆效果及应力波的传播分析 n 侵彻穿甲仿真，鸟撞及叶片包容性分析，跌落分析 n 失效分析，裂纹扩展分析 n 刚体运动、刚体柔体运动分析 n 实时声场分析 n BEM边界元方法，边界元、有限元耦合分析 n 光顺质点流体动力（SPH）算法 93 n3.

58、 热分析 n 稳态、瞬态温度场分析 n 热传导、热对流、热辐射分析 n 相变分析 n 材料性质、边界条件随温度变化 n4. 电磁分析 n 静磁场分析计算直流电（DC)或永磁体产生的磁 场 n 交变磁场分析 计算由于交流电(AC)产生的磁场 n 瞬态磁场分析计算随时间随机变化的电流或外 界引起的磁场 n 电场分析用于计算电阻或电容系统的电场. 典 型的物理量有电流密度、电荷密度、电场及电阻热等。 n 高频电磁场分析用于微波及RF无源组件，波导、 雷达系统、同轴连接器等分析。 94 n5. 流体动力学分析 n 定常/非定常分析 n 层流/湍流分析 n 自由对流/强迫对流/混合对流分析 n 可压缩流/不可压缩流分析 n 亚音速/跨音速/超音速流动分析 n 任意拉格郎日欧拉分析（ALE） n 多组份流动分析（多达6组份） n 牛顿流与非牛顿流体分析 n 内流和外流分析 n 共轭传热及热辐射边界 n 分布阻尼和风扇模型 n 移动壁面及自由界面分析 95 n6. 声学分析 n 定常分析 n 模态分析 n 动力响应分析 n7. 压电分析 n 稳态、瞬态分析 n 模态分析 n 谐响应分析 n8.