摘 要:关键词:本文结合有限元计算结果的特点,发展了一套用于整理计算结果的理论算法及实现步骤。从而可以方便地生成计算结果的等值线图和块图。文中算法以有限元位移插值函数为基础从而和有限元理论形成完美统一,适合有限元结果可视化的需要。
有限元等值线图块图二维三维可视化
Abstract Inthis paper we deal with the problems in the finite element method (FEM)output result visualizations.W e can show the result of FEM by equal value lines graphs,filling contour graphs.These can help us analyzing the output data ofFEM.The methods of showing FEMmeshes,merging equal value lines and filling contours is introduced.The post-process characters of three dimensions solidelement and spaced plate and shell element are also discussed.
SAP ADINA MARC NASTRAN PAFEC NFAP NONSAP ANSYS ASKA SAP ADINA SAP84DDJW X=N X Y=N Y N i ,有限元自从本世纪三、四十年代发明以来在理论,工程实践及计算机程序标准化等许多方面取得了令人瞩目的进步。广大科技工作者,工程技术人员结合自己的专业工作编写了大量的有限元程序。国外在六,七十年代出现了一批以大型机为基本环境功能强大的通用有限元程序包,如,,,,,,,,等。八十年代随着工程工作站的迅猛发展,这些大型程序相继出现了工作站版本并增加了基于图形的前后处理系统。有限元程序前后处理能力已经成为和单元库及材料库同等重要的衡量
有限元程序包功能的重要标准。我国由于在有限元理论和应用方面起步都落后于国外,在大型有限元程序包方面基本上是以引进为主,现在广泛应用的系统有,,ASKA,NASTRAN等国外大型程序的早期版本,我国自主开发的有限元软件包有北大的,大连理工大学的等。这些软件包多数以批处理方式运行,其图形能力仅限于被动地在显示器,笔式绘图仪上产生线框图,图形功能较弱,远达不到国外运
行于工作站上的有限元软件包的图形水平,不能满足交互式检查计算结果的要求。对我国现有“落伍”程序包进行升级改造,使之拥有先进的图形处理能力这对提高我国有限元应用水平,节约大量引进先进有限元系统的费用具有深远意义。本文结合计算机图形学发展对有限元结果可视化过程中的算法进行了研究,并对一些问题进行了讨论。
有限元网格可以用线框表示单元边界,也可以用实体表面显示单元。为加快显示速度我们一般用线框图。由于显示器无论分辨率多高其总是离散的有限个象素点,画出各单元的边界就是画出表示单元边界的象素点,既进行光栅扫描的过程。我们知道二维单元内部任意一点的坐标:Σ,Σ,其中(ξ,η)是有限元的位移插值函数。这是用ξ,η表示的参数方程,单元表示见图1。画出单元边界,既是令=±1,取一系列ξ∈[-1,1]和令=±1,取一系列[-1,1]分别计算出相应的N值,进而计算出象素点坐标的过程。我们无法直接对单元边界进行光栅扫描转换。因于我们事先无法知道多大的ξ,η增量对应X,Y坐标下的一个象素点的间距(这和单元的形状,大小及显示比例有关),实际上我们无法用参
数方程直接得到这些边界上全部象素点。
这里我们采取“以直代曲”的办法。首先我们根据节点的坐标来判断将单元边界分为几个直线段,然后在取值范围内分几等份,求出,进而计算出各插值点相应的N值和坐标,最后以直线连接各等份插值点从而画出单元边界。
有限元计算和经典方法不同:经典方法结果是以满足微分方程及边界条件的解析式形式给出的;有限元计算结果是给出各离散点上的位移,应力,应变,温度,热流,电压,电流等物理量。我们要丛这些离散结果中总结规律其工作量很大。据文献
一、概述
二、有限元网格显示
三、计算结果等值线显示
i i i i i i ξηξ∈ξ,ηξ,η有限元结果处理过程的可视化
郭宇恒YuHengGuo
RESULTOFFINITEELEMENTMETHODVISUALIZATION
介绍有限元输出结果的判别和分析占分析工作总工作量的一半。一个理想的办法是利用计算机图形生成各种物理量的等值线图和以颜表示的块图等物理量图,通过这些图来分析其中的规律。
普通等值线图绘制其离散点数据采样是在矩形网格上进行的,有限元数据由于网格是实在的且一般并非矩形,要重新构造网格并将有限元结果插值到新网格上去必然降低结果的精度。本文作者结合有限元特点发展了一种直接从有限元计算结果形成等值线图的新方法。
1三节点三角形单元等值线算法。三节点三角形单元等值线算法是其他有限元单元等值线算法的基础。其算法虽然简单但具有典型性,我们首先介绍该种方法。其算法如下:
首先遍历所有离散点上要生成等值线图的物理量,出最大,最小值H,H。若以等分方式形成等值线图则按H=H+i*(H-H)/n,(i=0,n);若以不等分方式形成等值线图则输入,以确定各等值线值。
判定三角形三个顶点物理量,,是否有相差极小的情况出现,出现则说明在单元内物理量变化极小不可能有等值线通过,则对下一个单元执行本步操作,否则将三角形三个顶点排列为P,P,P并使其节点物理量满足H≥H≥H关系。
判断所有等值线和当前单元是否相交。若不成立,对下一个单元执行上步操作;若成立,则H通过当前单元,再判断是否,成立则情况见图-2中线;否则情况见图-2中线。
利用线形插值求,线和三角形边的交点。,点计算公式:
()()(),()()()。点计算公式:
()()(),()()()。点计算公式:
()()(),()()()。连接或既是当前单元的等值线
。连接所有成对的交点既得当前单元所有等值线。
对所有单元循环执行至即可得到等值线图。
这里我们是在直角坐标系插值,对三角形也可采用三角形面积坐标L1,L2,L3进行插值,原理相同,区别仅为用L1,L2,L3代替X,Y。
2-节点三角形单元等值线算法。采用三角形面积坐标,,表示其形函数,为单元的节点数,具体表达式见参考文献。首先用,,三条参数曲线分割大三角形为四个小三角形区域见图-3,其次将四个小三
角形区域分别按1中步骤进行处理。为不失一般性4-6节点三角形应按曲边三角形考虑,这时不应再在X—Y坐标系中进行线形插值,而应在L1,L2,L3参数域内进行插值,故在求等值线与单元边界点时不同于三节点三角形,现以3—6—5区域为例说明如下:
1L1,L2,L3参数域内坐标确定。等值线与边界5—3的交点坐标:L1=0,
L2=1/2+(0-1/2)(Hi-H5)/(H3-H5),
L3=1/2+(1-1/2)(Hi-H5)/(H3-H5)。
等值线与边界3—6的交点坐标:L2=0,
L1=1/2+(0-1/2)(Hi-H6)/(H3-H6),
L3=1/2+(1-1/2)(Hi-H6)/(H3-H6)。
等值线与边界5—6的交点坐标:
L3=1/2,
L2=1/2+(0-1/2)(Hi-H5)/(H6-H5),
L1=0+(1/2-0)(Hi-H5)/(H3-H5)。
2计算形函数的值。Ni(L1,L2,L3)。
3计算节点坐标。
X=Ni(L,L,L)Xi,Y=Ni(L,L,L)Yi。
.(1)(2)H H H (3)i K M (4)K M K M X=X2+X1-X2*Hi-H2/H1-H2Y=Y2+Y1-Y2*Hi-H2/H1-H2K2X=X3+X1-X3*Hi-H3/H1-H3Y=Y3+Y1-Y3*Hi-H3/H1-H3M2X=X3+X2-X3*Hi-H3/H2-H3Y=Y3+Y2-Y3*Hi-H3/H2-H3(5)K1K2M1M2Hi (6)(2)(5).46L1L2L3Ni i L1=1/2L2=1/2L3=1/2()()()123123MAXMINi12312312311MINMAXMINiiiHHH≥H≥HH≥H≥HH≥H≥HΣΣ12312312
如果单元节点数是4个或5个,则可利用线性
插值补插1个或2个节点,从而形成6个节点的三
角形,再按上述步骤进行处理。这里假定节点3,
6,5的H值是按降幂排列的,若非如此只要改动
节点顺序即可。
34—9节点四边形单元等值线图算法。对
4节点四边形单元可连接任意两个对顶点,以形成
两个三节点三角形,由于连接的方法不同而生成的
等值线也不同,具体见图-4a和图-4b。在网
格较大情况下会出现对称计算结果等值线明显不对
称情况,图-4b。按图-5形成四个三角形可克
服上述问题,但三角形数目增加一倍且须补插中心
点。对5-9节点四边形单元可分为8个三角形,
再求交点坐标。同样计算也必须先在参数域上对
进行线形插值。
有限元计算结果不但可以用等值线方法表示,
同时还可以用块图表示。这种方法往往更加直
观,形象。块图表示物理量场既是把处在某一值
范围内的物理量用特定的颜块来表示,也就是
在相邻等值线间填充特定颜。在许多计算机编程
语言中图形填充函数都要求指定填充区域边界和填
充区域内一点,即“填充种子”。一种填充方法是
先形成整个区域的等值线图,然后在相邻等值线间
调用计算机图形学的填充算法完成。该种方法具有
一次填充相邻等值线间所有单元,填充速度快,不
易溢出等优点,但给定“填充种子”有一定的难
度。第二种方法是以单元为基础的,既先形成单元
等值线,然后再填充单元,逐个进行直至完成所有
单元。本文采用了第二种方法。下面以图-2的三
角形为例进行说明。
1先用形成等值线的算法计算出当前单元内所
有等值线。假设单元内有两条等值线,其M线为H
i,K线为Hi+1。
2连接P1和P2P3的中点,并在该线上接
近P1点的附近选取一点作为当前单元的“填充种
子”,一般该点距P1点大约3-5个象素点即可
保证该点在P1P2和P1P3之间。
3首先以上步选择的象素点为种子在整个三角
形中用界于Hi,Hi-1间的颜填充形成图-
6a。再用和单元边界同种颜画出等值线Hi。
4用同一种子在区域P1M1M2P3内,以
界于HiHi+1的颜填充,丛而形成图-
6b。
5用同一种子在区域P1K1K2内,以界
于Hi+1Hi+2的颜填充,丛而形成图-
6c。
这里由于对当前三角形的多区域填充采用了先
填充“低值区”后填充“高值区”逐次覆盖方法从
而仅指定一次“填充种子”,降低了算法的难度。
对于以上步骤中较小三角形或填充区域较小的情
况,这时若不能保证种子在域内则该用按光栅扫描
线顺序逐点填充该区域。shell程序的编写流程
上面我们结合二维问题讨论了有限元计算结果
的可视化。对于三维实体问题同样可以解决。其问
题的关键是首先要解决遮挡问题,也即图形学中的
消隐问题。三维有限元模型中,要处理的表面数量
十分惊人。例如图-7所示,一个由个“砖块
单元”构成的中小型模型的四边形表面数既达到
个至多。何况每个表面上又有多种颜或多根
等值线,直接处理计算
量太大。结合有限元网
格的特点我们提出以下
方法
1将各单元表面分
类为物体表面的面,物
体内部的面。分类的准
则是单元间的表面均为
物体内部的面,非单元
.
.
.
.
.
.
1000
6000
.
四、计算结果块图显示
五、三维问题的表示方法
ξ,η
图-4a     图-4b    图-5 
图-6a   图-6b   图-6c
图-7图-1     图-2   
     图-3
间的表面均为物体表面的面。具体判断可按构成单元表面的节点号集合进行,凡两个和两个以上表面的节点编号集合相同则为物体内部的面。否则为物体表面的面。经过这样处理后要进行消隐处理的表面降为个。
2对物体表面的面进行进行遮挡判断,具体方法可计算机系统的图形能力而定。对象SGI等具有固化Z-BUFFER的机器,可采用此算法以加快消隐速度。
3在可见表面上按和上面二维问题类似的方法进行显示处理。
对空间膜元和壳元也可采用类似办法,但必须增加Z坐标。同样在有遮挡情况时应先进行消隐。
依据上述理论,我们在中文版WINDOWS32环境用BORLAND C++语言进行了实现。取得了预想的效果。
例一是一个平面应力问题的小孔应力集中问题,一个长宽均为10CM,厚为1CM的A3钢板上钻了一个1CM的圆孔。在左右两边上均匀作用集度为100MP载荷。整个网格共有310个节点,552个三节点三角形单元,图-8a表示的是其X方向位移图。图-8b表示的是其x图。图中应力集中因子等于3十分清楚。可见即使在这样网格大小变化很大,应力场极不均匀的情况下仍能得到理想的结果,这说明上述算法是成功的。
二是一个平面热传导问题,图中突出部分是墙体,其下扩大部分代表基础。墙体的右侧为室内,其墙体和地面表面温度均为20摄氏度,左侧为室外,表面温度是零下30摄氏度。整个网格共有1181个节点,364个8节点四边形单元。图-9是墙体和基础内的温度场。
600 .
. .?
六、应用举例
(作者单位:510620广州智海建筑工程技术有限公司;
2006.5)
[1]孙家广陈玉健黄汉文计算机辅助设计技术基
础清华大学出版社
[2]任仲贵张关康CAD/CAM原理清华大学出
版社
[3]谢贻全何福保弹性和塑性力学中的有限单元
法机械工业出版社
[4]王来生计算机图形处理技术西安交通大学出版
[5]王勖成邵敏有限单元法基本原理与数值方法清
华大学出版社
[6]赵兴华徐福娣线形与非线形有限元分析(AD
INA程序理论部分)
参考文献
...
.,1990.9
...
,1991.9
.
.
.,1981.8
..
,1988.6
...
,1988.9
..
,1982.3
图-9
图-8a
图-8b