4 有限差分法与有限单元法,4.1 有限差分法,4.2 有限单元法,4.1 有限差分法,有限差分方法是数值计算中应用非常广泛的一种方法，是求解微分方程的主要方法之一。其实质就是以有限差分代替无限微分、以差分代数方程代替微分方程、以数值计算代替数学推导的过程，从而将连续函数离散化。以有限的、离散的数值代替连续的函数分布。,4.1.1 概述,4.1.2 有限差分法的主要步骤,1、构成差分格式,差分方程通常是一组数量较多的线性代数方程（即：线性方程组）。其求解方法有下列两种：（1）精确法，又称直接法，即消元法；（2）近似法，又称间接法，即迭代法。,2、求解差分方程,3、对所得到的数值解进行精度与收敛性分析和检验。,首先选择网格布局、差分形式和步长；其次，以有限差分代替无限微分，即以 代替dx以差商 代替微商（导数） ，以差分方程代替微分方程及边界条件。,4.1.3 差分方程的建立,建立差分方程是有限差分法的关键环节。导出差分方程的途径可有两种：,（1）从微分方程出发，以泰勒级数截断，从有限差分的数学含义去建立有限差分和差分方程。,（2）从由网格所划分的单元体的能量平衡分析出发、由积分方 程去建立差分方程，该方法又称单元体平衡法。,两种方法各具特色，但无论采取何种差分方程的推导方法，在建立差分方程前，均需对所论区域进行离散化 。,差分方程的建立过程（之一）,合理选择网格布局及步长,在实施有限差分法中首先在如图4.1所示的求解区域内，将自变量x，y 分别沿x，y轴方向的连续变化，离散为x0，x1 ，x2，xn及y0，y1 ，y2，yn个不连续点形成离散化网格；网格交点称为结点(或节点)，依次将结点编号，与区域自变量离散化相对应，区域内函数也将同时被离散化 。,离散化后各相邻离散点之间的距离，或离散化单元的长度称为步长，步长的大小可以是常量，也可以是变量。网格的粗细与是否均匀，要根据求解区域物理场的实际分布和对结果所要求的精确度而定。 一般说来，对均质、形状简单且规则、物理量变化不剧烈的物体或求解精度要求不高时，可采用等步长、大步长，即采用均匀网格；而对形状复杂、组分不同、物理量变化剧烈的物体，或求解精度要求较高时，则采用小步长、变步长。,图4.1 求解区域离散化,合理选择网格布局及步长（续）,另外，对一些较复杂的问题，在选择网格与步长前，往往要对所论区域的物理场作出粗略估计，然后以较粗的网格、较大的步长计算出参考性物理场，根据这一参考性物理场再选择合理的离散化网格。,离散化网格的布局，要根据所要求解的问题的性质及求解要求确定。一般说来，有两种方法：,（1）物理划分法：这种方法是根据问题的物理特性划分，如建筑物墙壁内外层面砖、普通砖和内灰泥层组成；若拟求各层界面壁温，则离散化时应按不同材料组分划分区域。,将微分方程转化为差分方程,差分方程的建立过程（之二）,微分方程转化为差分方程实际上就是以差分代替微分、以差商代替微商的过程，是以有限小量去代替无限微量的近似化过程。,方法：写出微分方程中各微分与微商所对应的差分与差商形式，代入原微分方程即可。,差 分,某物理量的有限增量,差 商,函数的差分与自变量差分之比,用差分代替微分方程中的微分，用差商代替微分方程中的微商，即可将微分方程转化为差分方程。 差分方程通常是一个线性方程组，利用以前介绍的直接法（消元法）或间接法（迭代法）即可解之，从而得到原微分方程的解。,4.2 有限单元法,有限单元法（又称为有限元素法，简称有限元法），是20世纪50年代初才出现的 一种新的数值分析方法，最早应用于航空航天领域，主要用于力学与结构分析中， 20 世纪 70 年以来被应用到传热学计算中。与有限差分法相比较，有限元法的准确性和稳定性都比较好，且由于其单元的灵活性，使它更适应于数值求解非线性热传导问题以及具有不规则几何形状与边界，特别是要求同时得到热应力场的各种复杂导热问题；有限元法在传热学中的应用正处于开拓与发展阶段，迄今为止，其应用已波及热传导、 对流传热及换热器设计与计算。,4.1.1 概述,从“有限元”的名字出现到今天，经历了几十年的发展，其基本理论已经日趋完善，复杂非线性问题的各种算法得到很大的发展，并且在工程领域（如：结构力学、热传导、电磁场、流体力学等连续域问题）得到广范的应用。,4.2.2 有限元法的基本思想,（1）假想把连续系统（包括杆系，连续体，连续介质）分割成数目有限的单元，单元之间只在数目有限的指定点（称为结点）处相互连接，构成一个单元集合体来代替原来的连续系统。在结点上引进等效载荷（或边界条件），代替实际作用于系统上的外载荷（或 边界条件）。 （2）对每个单元由分块近似的思想，按一定的规则（由力学关系或选择一个简单函数）建立待求未知量与结点相互作用（力）之间的关系（力位移、热量温度、电压电流等）。 （3）把所有单元的这种特性关系按一定的条件（变形协调条件、连续条件或变分原理及能量原理）集合起来，引入边界条件，构成一组以结点变量（位移、温度、电压等）为未知量的代数方程组， 解之就可得到有限个节点处的待求变量 。,可见：有限元法实质上是把具有无限个自由度的连续系统，理想化为只有有限个自由度的单元集合体，使问题转化为适合于数值求解的结构型问题。,实例1,实例2,图4.4,4.2.3 有限元法的基本概念与一般步骤,直接刚度法,在刚提出有限元法的时候采用的是直接刚度法它源于结构分析的刚度法。因刚度法只能处理些比较简单的实际问题，现在已很少使用了，但它对我们理解和明确有限元法的一些物理概念是很有帮助的。所以我们首先通过一个例子来介绍直接刚度法，同时说明有限元法求解的一般步骤。,直接刚度法（例）,考虑一个变截面杆，如图4.5所示，杆的一端固定，另一端承受 P1000N的载荷杆的顶部宽w1=2cm， 杆的底部宽w2=1cm ，杆的厚度t=0.125cm，长度L 10cm、杆的弹性模量E10.4106MPa。试分析该杆沿长度方向不同位置的变形情况，假设杆的质量可以忽略不计。,图4.5 受轴向载荷的变截面杆,有限元基本概念与一般步骤,1 前处理过程,图4.6 将杆划分为单元和结点,先将求解的问题分解为结点和单元。为简单起见，将杆划分成五个结点和四个单元（如图4.6所示）。 给定的变截而杆简化为四个独立的部分，每部分的截面面积恒定（为组成该单元的两个结点处的面积的平均值）。,(1) 求解域离散化,（2）直接刚度法分析（结点分析法）,图4.7 具有均匀截面的固体单元在力F作用下的变形,单元中的平均应力为：,平均正应变为：,在弹性范围内，由Hooke定律E可得：,注意：上式与线性弹簧等式F=kx相似。 因此可用弹簧的变形来模拟固态单元的变形，从而得：（keq称为等价刚度）,将上述结果扩展到整个变截面杆上，则：杆可以用一个由四个弹簧（五个结点）串联组成的模型来表示，每个单元模型的弹性行为可以用等价线性弹簧来表示：,根据静态平衡的要求：作用在每个结点上的力的总和为零，因此得到如下方程：,将这些方程进行变化可得到：,矩阵形式为：,在载荷矩阵中区分施加力与反作用力也是必要的。因此矩阵式又可写为 ：, 反作用力矩阵 刚度矩阵 位移矩阵 载荷矩阵 ,对我们所讨论的问题，因为杆的顶端固定，结点1的位移应力0，即 u10，将此边界条件用于上式，得到如下矩阵：,解此方程组即可得到各结点的位移,（3）有限单元法分析（单元分析法）,讨论 一般单元的单元刚度矩阵，并讨论总体刚度矩阵的集成。,先以一个单元（含两个结点）作为研究对象，其传输力可用图4.9表示。,确定每个单元的方程：,集成单元：,首先写出各单元刚度矩阵及其在总体刚度矩阵中的位置如下：,单元 1,单元 2,单元 3,单元 4,组装（相加），即得总体刚度矩阵：,施加边界条件和载荷,杆的顶端固定，应满足边界条件 u1=0载荷P 施加在结点5上。应用这些条件得到下式 ：, 刚度矩阵 位移矩阵 载荷矩阵 ,有限元分析的一般格式（固体力学问题）,解此方程组即可得到各结点的位移,2 求解阶段,杆的截面沿 y 轴的变化可以描述为：,计算得：,由等价单元刚度（刚度系数）公式：,计算得：,各单元刚度矩阵为：,解此矩阵方程（线性方程组），即可得到各结点的位移值：,将单元刚度矩阵组合，得到总体刚度矩阵：,应用边界条件并施加载荷，得到：,3 后处理阶段,因为结点的位移已知、上式也可直接从应力应变关系得到：,对分析结果进行后处理，我们可以得到其他一些信息。如：每个单元中的平均正应力。这些值可以从下式得到：,对本例，可计算出各单元的平均正应力值如下：,4.2.4 有限元程序的结构和特点,有限元法程序总体可分为三个组成部分：前处理部分、有限元分析本体程 序、后处理部分。,有限元法的实现必须通过计算机，全部有限元法的计算原理和数值方法集中反映在有限元法的程序中，因此有限元法的程序极为重要。它应具有分析准确可靠、计算效率高、使用方便、易于扩充和修改等特点：,（1）特点,（2）结构,对于一个实际的工程问题，离散模型的数据文件十分庞大，靠人工处理和生成一般是不可能的。为了解决这一问题，有限元分析程序必须有前处理程序。前处理程序是根据使用者提供的对计算模型外形及网格要求的简单数据描述，自动或半自动地生成离散模型的数据文件，并要生成网格图供使用者检查和修改。这部分程序的功能很大程度上决定了程序使用的方便性。,有限元分析本体程序是有限元分析程序的核心，它根据离散模型的数据文件进行有限元分析；有限元分析的原理和采用的数值方法集中于此，因此它是有限元分忻准确可靠的关键。选用计算方法的合理与否决定了有限元分析程序的计算效率和结果的精度及可靠性 。,有限元程序的结构和特点 （续1）,同样，有限元分析程序的计算结果也是针对离散模型得到的。例如静力平衡问题可以得到离散模型各结点的位移、各单元的应力等，输出的文本文件量很 大，但却不易得到所分析对象的全貌，例如位移哪里最大、应力集中发生在什么部位以及变化趋势如何等。因此一个使用方便的有限元分析程序木不仅要有可供选择输出内容的文本文件，还需有结果的图形显示，如位移图、等应力线图或截面应力分布图等；这部分程序称后处理程序。与前处理程序相似，后处理程序对程序使用的方便性有举足轻重的作用。,有限元分析程序的三个组成部分对于一个较好的用于实际问题分析的有限元程序来说，前后处理的程序量常常超出有限元分析的本体程序，前后处理功能越强，程序的使用就越方便。有限元分析程序中前后处理程序一般可占到全部程序条数的 2/3 4/5。有的近期发展的通用程序更注重程序的“包装”和使用功能。有限元分析本体程序以外部分的比例更高。,有限元程序的结构和特点 （续2）,ANSYS（Analysis System）世界著名力学分析专家、匹兹堡大学J.Swansan教授创立的SASI（Swansan Analysis System Inc.）的大型通用有限元分析软件世界最权威的有限元产品。 SAP （Structural Analysis Program）美国加州大学伯克利分校M.J.Wilson教授的线性静、动力学结构分析程序。 NASTRAN（NASA Structural Analysis ）美国国家航空宇航局（NASA）的结构分析程序。 ADINA（A Finite Element Program for Automatic Dynamic Incremental Nonlinear Analysis ）美国麻省理工学院机械工程系的自动动力增量非线性分析有限元程序。 IDEAS（Integrate Design Engineering Analysis System）美同SDRC公司的机械通用软件，集成化设计工程分析系统。集设计、分析、数控加工、塑料模具设计和测试数据为一体的工作站用软件。 AGOBJ美国AGOBJ 公司在SAP5和ADINA 有限元分析程序的基础上针对微机平台开发的通用有限元分忻系统。,