有限元法的前后处理

8 有限元法的前后处理8-1 引言有限元法分成以下几步：（1）计算模型的几何剖分、数据生成和数据准备；（2）计算分析；（3）计算结果的分析、整理和图形显示。有限元发展初期，（1）、（3）都是由人工来做。有统计资料表明，（1）的工作量为45%，（3）的工作量为50%，而（2）只占到5%。 以（1）为例，首先按比例画到坐标纸上，然后按一定的顺序编号，再整理出如下信息，以供有限元分析之用：*节点信息节点编号和节点坐标；*单元信息单元编号和单元中节点号的排列顺序；*材料信息计算模型的材料性质（弹性模量，泊松比，比热，导热系数等）；*载荷信息计算模型所受的负载信息（集中力，体积力，表面力，温度，压力等）；*约束信息初始条件和边界条件。 最后还要将这些数据一一输入计算机。以上工作乏味而且容易出错。有限元分析的前处理就是使计算机部分或全部完成计算模型的几何剖分、数据生成和数据输入，有限元的后处理则是将有限元计算结果由计算机整理成易于阅读或分析的数值或图形形式。8-2 有限元分析的前处理技术有限元法的前处理主要有以下内容：（1）计算模型的几何表示；（2）模型网格的自动分划（或剖分）；（3）刚度矩阵的带宽优化；（4）模型网格图的计算机绘制。 一 计算模型的几何表示1 对于设计阶段的零部件进行有限元分析，如果采用计算机辅助设计，则模型的几何表示可采用计算机造型系统中的几何表示。常用的有表面模型和立体模型。2 对于已有的零部件做有限元分析，既可以将其输入到计算机中，采用第一种办法表示其几何形状，也可以根据零件的几何形状来决定其表示法。常用的有整体表示法和分块表示法。（1）整体表示法：用点表示线，用线表示面，用面表示体。（2）分块表示法：把整体看成是由简单个体的组合，而简单个体则可用多边形和多面体表示。 二 有限元网格的自动剖分有限元网格的自动剖分与计算模型的几何表示方法有密切的关系。整体表示的几何模型，适宜于采用整体剖分要用到较多的数学知识。下面讨论分块表示几何模型时常用的一种分块剖分法。以二维问题为例。图示的二维区域可以看成是A,B,C三个部分组成。 A,B,C都是简单的四边形，四边形顶点的坐标可以表示其形状和位置。 网格自动剖分的方法分成四个步骤：1）分块映射；2）网格剖分；3）顺序编码；4）总体合成。1 分块映射首先把计算区域手工粗分成若干个四边形区域的组合，每个四边形称为一个大单元。 B CA 为了既能表示直边大单元，又能表示曲边大单元，采用8个节点来描述一个大单元。右图中，（a）是真实图形在总体坐标系下的样子，经过等参数单元的变换后，得到（b）图中边长为2的正方形。同样的方法可以得到B和C的映射结果。 2 3 467 11 2 3 4 5567 88yO OA )(b)(a x 2 网格剖分网格剖分是在映射后的正方形区域中进行的。剖分时可以有两种方法：等边剖分和不等边剖分。（1）等边剖分将=-1边和=1边剖分成相等的份数。同样，将=-1和=1边也剖分成相等的份数。则最终剖份单元数为 。 NN O （2）不等边剖分某一方向两个边界上的剖分数不一致。比如 =-1分割数为3， =1一侧的分割数为7。需要计算沿另一个方向（等值、 方向）的分割数，使相邻的两条等值线的分割数差1。 O3 4 6 75CBA 对于相邻的两条等值线分割数不相等的情况，应将分割数多的等值线上的最后两点与相邻的、分割数少1的等值线上的最后一点组成一个三角形单元。已经形成的四边形单元，只要连接四边形的两个对角点，即可得到两个三角形单元。或者保持四边形单元，这样，求解区域既有四边形单元，也有三角形单元，相应的单元刚度矩阵也有两种类型，总体合成时，按照总体节点编号累加单元中的元素值即可。完成局部坐标系o中的大单元的分割后，对应的总体坐标系xoy中的大四边形也同时被分割了。从o到xoy的坐标变换关系是 式中 为映射函数，也就是8节点四边形等参数单元位移插值函数的形函数。 是总体坐标系中大四边形各边节点的坐标。 ii ii ii yNy xNx 8181 , 8,2,1,iNi ,84 1121 2,6 1121 7,5,3,1 11141, 22 iiiN ii iiiii 8,2,1,iyx ii 按照坐标变换关系式，将已经分割后的局部坐标系中的正方形中的网格节点坐标(,)代入，即可得到对应的总体坐标系下的节点坐标(x,y)。例如， 方向分割数N=4， 方向的分割数N=3，求图中K点的坐标(xK,yK)。 KO 如果以O为起点，则K点的局部坐标为根据已知的 ， 将(K,K)代入坐标变换式，即可求出对应的(x K,yK)。 方向的分割线序数。方向和点算起从，；起算点的局部坐标 OKK KN KNOO OK OK 1,1, 3113212 2 114212 8,2,1,iN i 8,2,1,iyx ii 3 顺序编码编码就是将所有节点连续无缺地顺序编号，并对每个单元编号和形成单元信息。编码可以分成两步来做。首先对每个大四边形编码。对基本的大四边形单元编码时，可以按行（列）进行。以按行编码为例，每计算一个节点，节点号加1。单元编号也可类似进行（比如按列）。 )2( )3(1k2k )1( )4( )5(5k )1(NS )2(NS )3(NS )4(NS 为了形成与某一个单元相关的节点号，在节点编号时，对每一个新行的第一个节点专门形成一个数组NS(i)。其中的i代表节点编码的行号， NS(i)的值代表该行起点节点号。如果单元编码按列进行，则第二个单元的节点是NS(3)， NS(3)+1， NS(2)+1， NS(2)。如果单元是8节点等参元，还要计入单元每边中间一点的节点号。对每个大四边形顺序编码，这样，可能会有一些节点在大四边形组合后，同一节点拥有两个节点号。这个问题可以在下一步总体合成中予以解决。 4 总体合成将经过剖分、编码后的大四边形单元合在一起，调整具有多个节点号的节点编码，重排节点号，称为总体合成。总体合成的主要任务是使整体节点号的编号连续无缺，并修改分开编码形成的单元信息。对于同一节点具有多个节点编号多重编号的节点，可以采用比较节点坐标的办法将其消除。比如，如果同时成立，说明i和j重合或很接近。 yji xji yy xx 21 21 式中(xi1,yi1)为第一大单元中的节点i的坐标；(xj1,yj1)为第二大单元中的节点j的坐标； x,y是预先给定的距离很小的数值。找出重复编号的节点之后，修改第二大单元中的节点信息和单元信息：（1）保留节点号i，取消节点编号j；（2）将第二大单元中比i大的节点编号全部减1；（3）修改第二大单元中该节点所在的单元信息和比该节点编号要大的那些单元的信息。重复以上过程，直至第二大单元中的所有节点都循环完毕。 对于三维问题，也可以实现有限元网格的自动分割，只是分块区域为六面体或五面体，映射后的区域为正六面体和直三棱柱。 1 2 3 4 56 1016 1511 17 18 19 2021 2224 2627 2930 3233 34 35 1823262932 33 2528大单元连续编号重复编号的剔除222123 3134 三 刚度矩阵的带宽优化为了节省计算机的存储空间，需要优化总刚矩阵的带宽。节点编号对应总刚矩阵中的非零元素的位置。节点编号变动，总刚矩中非零元素的位置也发生变动，反之亦然。带宽优化的原理：通过调整总刚矩阵中非零元素的位置，对应地修改单元信息，从而减少刚阵带宽。 为此引入邻接矩阵B，其阶数与刚度矩阵相同，其中的元素非0即1式中，aij是刚度矩阵中的元素。B也具有带状的样子。带宽优化有许多实用算法，有些要用到图论或较深的数学知识。下面介绍两种易于理解的带宽优化方法，它们都是采用变换邻接矩阵中的行或列的办法来减少刚阵带宽的。 0 0 0 1 ijijij aab 1 罗森算法罗森（Roson）算法的基本思想：找出引起B矩阵最大带宽的一对顶点，取其中一个与其他顶点互换，看能否减小带宽。若可以，则进行交换。交换后，再继续寻找引起当时B矩阵最大带宽的一对顶点，并重复上述过程。为此，首先将B中的“顶点”序号纪录在一个数组中，称为“顶点表”。“顶点”是B中每一行离主对角线元素最远的那一点。然后在顶点表中：（1）确定B矩阵的带宽及两个达到最大带宽的第一对顶点。在这对顶点中，若较高编号顶点能与一低编号顶点互换以减小带宽，则转向（6）； （2）若较低编号顶点能与一高编号顶点互换以减小带宽，则转向（6）；（3）若较高编号顶点能与一低编号顶点互换二保持带宽不变，则转向（6）；（4）若较低编号顶点能与一高编号顶点互换而保持带宽不变，则转向（6）；（5）如果在（3）、（4）步中已连续执行过规定的经验次数，或已交换过的顶点又重新选出来交换，则算法终止；否则，转向（6）；（6）执行所指出的顶点交换，即分别交换顶点表的对应分量以及B矩阵的行和列，并转向（1）。可见，整个交换过程都在顶点表中进行。 2 阿基茨厄特库算法阿基茨(Akyuz)厄特库(Utku)算法也在1968年提出一种减小平均带宽的方法，简称AU算法。平均带宽的定义：AU算法的基本步骤：（1）从B矩阵中取出相邻的两行列进行交换，并计算平均带宽。如果满足下列两个条件之一，则交换有效：行的带宽。矩阵第刚度矩阵的阶数；式中in n ini i 11 平均带宽减少；平均带宽保持不变，但有较多元素的行从矩阵中心向外移。（2）在一个指定的交换循环内，按(1) 及(1) 执行行列交换，其交换顺序规定为(1,2), (n,n-1), (2,3), (n-1,n-2),直至中心行。（3）如果在一个循环内没有发生交换，或经过经验次数 次的循环而平均带宽不减小，交换运算停止；否则，重复执行(1)、(2)。一般来讲，实用的带宽优化算法并不能使刚度矩阵的带宽达到极小化，但并不影响其实用性。大型的有限元分析软件包一般都具有带宽极小化的功能。1003 n 四 模型网格的自动生成目的：检查所生成的几何数据的正确性。方法：根据单元数组信息（单元和单元节点关系）和节点数组信息（节点号和节点坐标关系），将相邻节点连接起来，一个单元一个单元地去形成有限元剖分网格。注意：所得结果还要消除隐藏线，否则看不清结果。这部分内容属于计算机图形学，在此不做讨论。 8-3 有限元分析的后处理有限元数据后处理的原因：1）计算结果是位移和应力；2）数据量大，不直观。有限元数据处理工作包括两个方面；1）数值处理 将有限元结果转化成工程中常用的形式，或设计师熟悉的形式；2 ）图形处理 用图形直观地表示设计结果，使设计结果一目了然。 一 数值处理以应力分析为例。有限元位移法计算得到节点的位移，进而计算单元的应力：单元内任意一点处的应力为 。上述结果不论是数值积分点处的应力，节点处的应力，还是单元中心处的应力，也不论其精度高低，都不好使用，原因是：我们需要的不是应力分量，而是主应力或等效应力。说明：单元内不同点处的应力值的精度是有差别的。一般，数值积分点处的应力值较为准确。 )(, eBDzyx zxyzxyzyx , 1 主应力计算和等效应力计算工程中往往用一点的主应力 来计算或判断结构的安全程度。对于二维问题，有限元分析可得任意一点的 ，则主应力 与它们的关系为主应力方向 321 , xyyx 和, 21, xyyxyx xyyxyx 22 21 22 22 yx xy 2arctan21 最大剪应力和最小剪应力为对于三维应力问题，主应力可以通过求解下述三次方程得到：2 21minmax 2223 22221 322132 0 xyzzxyyzxzxyzxyzyx zxyzxyxzzyyx zyxIII III 其中 上式的解可表示为将 按从大到小的顺序排列则 是最大主应力， 是最小主应力。 23221 32131221 13 12 11 32 2792cos ,332 34cos3 32cos3 3cos3 II IIIIIIR RI RI RI 其中 321 , 321 1 3 从理论上可以证明，三个根都是实根，因此有否则，一定是有限元计算结果有误。主应力的方向可由以下三式中的任意两式三式联立，得出主应力的方向余弦1cos1- ,03 221 II 1 000 222 iii iiziyzizx iyziiyixy izxixyiix nml nml nml nml以及 3,2,1, inml iii 上式中前两式与第四式联立求解的结果为令 ，即可得到 的方向余弦。 21222 212 212 12 12 1 11 xyiyixixyzxyzx iyzxyzxy xyiyixi ixyzxyzxi iyzxyzxyiA An Am Al 式中3,2,1i 321 , 另外，要核算所分析的零件强度是否足够，还要根据一定的强度理论来计算某点的综合应力。各种强度理论的等效应力如下： 213232221xd3 31xd3 321xd2 11xd21 第四强度理论第三强度理论第二强度理论第一强度理论 2 应力修匀有限元保证位移在单元边界的连续性。但是，应力在单元边界并不连续。怎样处理计算结果才能使应力连续和直观呢？（1）简单平均或加权平均以三角形单元为例。三角形为常应力单元，但不同三角形单元的应力值不同。简单平均：加权平均： 有关单元面积所在单元面积节点应力值节点应力值有关单元节点应力值节点应力值 （2）最小二乘法在单元内采用最小二乘法修匀应力；单元节点的应力值，取环绕该节点的不同单元的应力值的平均（简单平均或加权平均）。单元内全应力修匀设应力向量 ，修匀后的应力向量为 ，则使A(e)达到最小。其中 )( 21)( eV Te dVA 用的插值函数。值函数和修匀应力时采有限元求解时的单元插节点应力值；应力值和修匀应力后的有限元求解得到的节点ii ii iiiiNN NN , , 根据上式是以修匀后的节点应力向量 为未知量的线性代数方程组。利用解析法或高斯积分，可以求得 ，进而利用单元应力插值，可以计算单元内任何一点的应力值。 )( )()( ,2,1 0 ,2,1 0 )( eV iT eie nidVN niAe 或 i i 单元内部分应力修匀只修匀单元中我们感兴趣的分量。例如，只修匀分量。这时有x )( 2)( 21eV xxxxe dVA )()( ,2,1 0 )( eV ixxxieii xixi xiixxiix nidVNANN NN e 带入变分后的表达式函数。修匀应力时采用的插值值函数；有限元求解时的单元插值；修匀应力后的节点应力应力值；有限元求解得到的节点将 上式为一线性代数方程组，方程的数目为n(e)，未知量为修匀后的单元节点应力分量 。单元内应力线性外推利用单元内的数值积分点处的应力值，外推单元节点的应力值，然后再平均与节点有关的单元节点应力。对于二维四节点等参数单元，若计算单元刚度矩阵时取了 个积分点，则积分点的坐标为xi22 ,31,31: ,31,31: ,31,31: ,31,31: DC BA 带入公式可以求出未修匀时积分点处的应力值设修匀后的应力在单元内按双线性变化，则修匀应力为 )(, eBDzyx DCBA , i i ii NNNN NNNN N 4321 43214321 41, O1 24 3A BCD 令积分点处的应力值等于修匀后的应力值，则 Pi i DCBADCBAP N NDNDNDNDN CNCNCNCN BNBNBNBN ANANANAN 1 43214321 4321 4321 4321 由此可以解出： 将上述求得的节点应力，再根据不同的单元节点应力值予以平均，就可得到最终的节点应力值。一般情况下，采用局部应力磨平处理，可以得到较好的结果，而计算量是很小的，所以得到广泛的应用。应力全域磨平和子域局部应力磨平可参见清华大学王勖成编著的“有限单元法基本原理和数值方法”。二 图形处理图形处理使计算结果简洁，直观，生动。常用以下几种方法。 1 变形图分为两种：1）用网格表示内部变形（当外部变形不大时）；2 ）用边界表示外部变形（只关注外部变形时）。两种方法，都把原图和变形后的图重叠画出，以作对比。绘图时注意：要放大，否则看不出来。例如，原节点坐标(x,y)，位移为(u,v)，则绘制变形时该节点的坐标为 变形放大倍数。式中QP uQyy uPxx, 2 等值线图等值线可以是等温度线，热流线，等应力线。方法：追踪法+插值法。思路：先找到一个等值点，然后根据离散节点坐标和节点的场变量值，求出节点连线间的等值点位置，并和开始找到的等值点相连。就这样一直寻找下去，直到边界。 jA Bia b 入口 入口 出口c 出口 插值算法：设两个节点i和j的坐标为(xi,yi),(xj,yj)，场变量(Ti,Tj)，待求场变量为Ta，则当成立时，在 边上有一等值点a，其坐标为对于二次插值单元，单元边界上有三个节点i,j,k，可采用二次插值或分段线性插值求出等值点。 0 jaia TTTTij ijij iaia ijij iaia yyTT TTyy xxTT TTxx 追踪算法：从A单元a点到b点，再从b点进入B单元，c为B单元的出口。求解关键：找出单元B。为此，只要搜索除A单元外的其他单元，这些单元的特点是：单元包含节点i,j；或单元节点坐标与i,j一致。使用中的几个问题：1）寻找起点没有规律可循。按单元循环，每个单元按边循环。2）等值线可能不只一条；3）只搜索没有搜索过的单元。 3 浓淡色彩图（彩色云图）对显示的区域不断剖分，直到一个像素的大小，同时用插值算法计算剖分后的小单元（像素）对应的场变量值，然后用当前场变量值对应的色彩，显示该像素。4 节点变量变化图采用曲线拟合或曲线插值技术，将所有离散时间点上的某节点场变量值，按一定的函数形式拟合成一条或几条连续曲线，并按比例画出。