用matlab进行有限元法编程.doc

用matlab进行有限元法编程Jack Chessa2002年10月3日1 引言本文的目的在于给用户在使用matlab编写有限元代码方面遇到的一些问题作一个简短的概述和指引。估计读者对于有限元方法的理论知识已经有了基本的认识，因此，本文主要关注有限元方法的实施方法。本文将使用一个用mablab编写的用于分析静态线性弹性力学问题的有限元代码示例，以给大家说明如何进行有限元法的编程。1.1 符号 为清楚起见，我们在本文中采用以下符号；粗体斜体字体V表示一个空间向量，其等同于在一个点上的位移或速度的空间维度，粗斜体字体d表示一个向量或矩阵的维数是未知数在离散系统即像STI内斯矩阵的系统矩阵，大写标表示节点编号，而小写标一般指沿笛卡尔单位向量的向量组成部分。所以，如果d是节点未知数的系统向量，uI是节点I的一个位移向量，那么 uIi是节点I在方向i上的位移向量即uI*ei。MATLAB的语法通常与数学符号混合使用，本文对此增加了增加了明确的解释。本文用打字机字体，字体指示MATLAB的语法被使用。2 有关编写matlab程序 matlab编程语言在说明如何进行有限元法编程方便十分有用，因为其允许我们进行非常迅速的编写数值计算，并且其具有广阔的prede NED数学库。另外，矩阵（疏与密），矩阵和许多线性代数工具极影呗定义了，开发人员可以完全集中在没有定义这些数据结构的算法的实现上。广泛的数学和图形功能进一步使开发者免于开发这些功能本身或是找到适当的现有的库。一个简单的二维有限元程序，在matlab中只需要几百行代码，而在Fortran或C+程序中可能需要几千行。尽管matlab编程语言在数学函数方面的功能非常完整，其依然有一些有限元的有利于发展为独立的功能的特定任务。这些已经被程序化，可参看前面提到的网站。通常有贸易FF这种易于开发。由于matlab是一种解释性的语言，每一行代码由matlab的命令行解释器解释并在运行时被顺序执行，因此，其运行时间可能比编译的编程语言，如Fortran或C + +，要长得多。应当指出的是内置的MATLAB功能已经编制并且效率极高，因此应尽可能多的使用。记住，是matlab的解说性质使得程序放缓，应不惜一切代价避免程序中的循环，尤其是嵌套的循环，因为这些可以使MATLAB程序的运行时间比可能需要的时间要长得多。通常的for循环可以通过使用matlab的量化处理消除。例如，下面的matlab代码设置矩阵A的行和列至零，并在对角线上放一个元素for i=1:size(A,2) A(n,i)=0;endfor i=1:size(A,1) A(i,n)=0;end A(n,n)=1;这种代码一定不能使用因为下面的代码A(:,n)=0;A(:,n)=0;A(n,n)=0;在三行解释代码中做的事情是一样的，如同nr+nc+1解释行，其中A是一个nrnc维矩阵。人们可以很容易地看到，当处理大型系统时（通常是用有限元代码的情况下），这可以迅速增加巨大的超载。for循环有时候是不可避免的，但令人惊讶的是，有几次是这种情况。建议人们当开发完一个matlab程序后，回去看看如何/是否可以消除任何循环。通过实践，这将成为matlab编程的第二特性。3 一个典型有限元程序的组成部分 一个典型的有限元程序包含以下三个部分1. 预处理部分2. 处理部分3. 后处理部分 在预处理部分定义了问题陈述的数据和结构会被定。这些措施包括有限元离散，材料的性能，解决参数等等。在处理部分，有限元素对象即STI内斯矩阵，计算力矢量等被计算，边界条件被实施，系统被解决。在后处理部分，处理部分的结果将被分析。这里应力可能被计算并且数据进行可视化。本文中，我们将主要关注处理部分。许多前期和后期处理操作已经在matlab中进行编程并且在网上有参考文献；有兴趣的人可以直接参看MATLAB脚本文件或在MATLAB命令行中键入帮助“function name，以获得有关如何使用这些功能的进一步信息。4 matlab中有限元数据结构在这里，我们讨论的是在有限元法中使用，尤其在示例代码中实现的数据结构。人们可以想出多种储存有限元程序数据的方式，这有些随意，而我们尽力使用最灵活而且最有利于matlab的结构。这些数据结构的设计可能取决于所使用的编程语言，但通常不会与这里列出的结构明显不同。4.1 节点坐标矩阵由于我们进行有限元法编程，因此 我们需要一种代表元素离散域的方式。这样做，我们定义一组节点和用一定方式连接这些节点的元素。节点的坐标存储在节点坐标矩阵中。这是一个简单的节点坐标的矩阵（可以想象一下）。这个矩阵的维是nnsdim，其中nn是节点的数量，sdim是空间维数。所以，如果我们考虑一个节点坐标矩阵的节点，第n个节点的y坐标是节点（N2）。图1显示了一个简单的有限元离散矩阵。对于这个简单的网格节点坐标矩阵如下：节点 = (1)4.2 元素连接矩阵 元素的定义存储在元素连接矩阵中。这是一个表示节点数目的矩阵，矩阵的每一行包含一个元素的连性数。因此，如果我们考虑连接矩阵元素，即描述了一个4节点四边形的网格,那么第36个元素是由连接向量元素（36，:)定义的，例如其可能是36 4213 14或元素连接节点36421314。因此，图1中的简单网格的元素连接矩阵是元素 = (2)注意：该元素的连接数都在逆时针排列的；如果没有这样做，一些雅可比矩阵的元素会是负数，而会导致刚度矩阵奇异（而且显然是错误的！）。4.3 边界的定义在有限元法中，边界条件用于组成力矢量（自然或诺伊曼边界条件）或指定边界上（必要或狄利克雷边界条件）的未知领域的值。在这两种情况下对于边界的定义是必要的。实现这一点的最通用的方式是保持必要的边界有限元离散。这个网格的维度将比问题空间维度小1（即二维边界网格对应于三维问题，一维边界网格对应于二维问题等等）。让我们再次考虑一下图1中的简单网格。假设我们希望边界条件适用于网格的右边缘，那么边界网格的定义由以下元素的2个节点的行元素的连接矩阵定义。右边缘 = （3）请注意，边界连接矩阵中的数字指同一节点坐标矩阵，就如同内部元素的连接矩阵中的数字一样。如果我们想要一个必要的边界条件适用于一个边缘，我们需要这个边缘上的节点列表。这可以在matlab中通过其独特的功能而很容易的实现。nodesOnBoundary = unique(rightEdge)；这将设置向量nodesOnBoundary等于24 6。如果我们想要在这个边缘上通过自然边界条件组成一个力向量，我们只需循环遍历该元素并整合该边缘上的力，就像我们将域内部，即刚度矩阵K，的任何有限元操作整合。4.4 自由度映射 最终对于所有的有限元程序，我们解决了线性代数系统向量d的形式。Kd = f (4)向量d包含的节点未知数的有限元逼近uh(x) = (5)这里NI(x)是有限元形函数，dI是节点I的节点未知数，其可能是标量或矢量的数量（若uh(x)是一个标量或矢量），nn是离散矩阵中的节点数量。对于标量场节点未知数在在d中的位置最明显的是如下：dI = d(I) (6)但在向量场中节点未知数dIi的位置这里I是节点数，i是向量节点未知数dI的组成部分，有一些歧义。我们需要定义一个从节点数和向量组成到节点位置向量d索引的映射。映射可以写作f:I,in (7)这里f是映射，I是节点数，i是组成，n是d中的索引。因此未知数uIi的位置表示如下：uIi = df(I,i) (8)有两种常用的映射。第一个是交替在节点未知向量d中的每个空间组成部分。在这样的安排下，节点未知向量d的形式如下：d = (9)这里nn是离散矩阵中的节点数量。映射是n = sdim(I - 1) + i (10)在这个映射下，i组成部分在节点I的位置是在d中确定的，如下：dIi = d(sdim*(I-1) + i) (11)另一种选择是令节点未知数的相像组成部分组成向量d中的连续部分，如下：d = (12)这种情况下的映射是n = (i - 1)nn + I (13)因此，这种结构中，组成部分i的在节点I上的位置由以下d确定，如下：dIi = d(i - 1)*nn + I) (14)当我们讨论元素操作符散射进系统操作符时，我们将采用后面的自由度映射。与这些映射相适应是很重要的因为它是一个定期执行在任何有限元程序中的操作符。当然，无论使用哪一个映射，刚度矩阵和力向量应当有相同的结构。5 有限元计算操作 有限元计算的核心计算式有限元的运算符。例如在一个线性矩阵中,他们将转换为外力的向量。有限元研究者将有限元的元素离散化，然后再将元素集中，然后将元素进入整个单元。这个过程是写入和组件。5.1 正交 一个元素的融合算子和一个合适的执行正交规则取决于元素和功能的整合。一般规则是一个积分如下在f()功能集成,q是正交分Wq正交权重。生成一个向量正交函数正交分和一个向量正交权值正交调制的规则。这个函数的语法如下quadWeights,quadPoints =正交(integrationOrder,elementType,dimensionOfQuadrature);所以举一个例子如下，将正交环f = x3的功能上的一条三角形元素qPt,qWt=正交(3,三角,2);For q = 1:长度(qWt)xi= qPt(问);%正交点球形坐标x %得到在正交分11人阵容%和导数行列式在正交点,jacf_int = f_int + x 3 * * qWt jac(q);结束5.2 算子“散射”一元算子的计算需要分散到整个球面系。散射元素力向如一个整体力矢量如图2所示。散射依赖在元件连接和自由度映射选择。下面的代码执行散射如图2所示elemConn = 元素(e,:); % 元素连接enn =长(elemConn);for I=1:enn; %在节点元素循环for i=1:2 % loop空间维度Ii=nn*(i-1)+sctr(I); %自由度f(Ii) = f(Ii) + f(i-1)*enn+I);endend但使用了一个嵌套的循环的(环中的环)。这是一个更异乎寻常的行为，考虑到这样一个事实时,即它发生一个元素循环,这真的是个能减慢执行时间的程序(由数量级在许多例)。我们将矩阵有四个嵌套的循环，刚度散射矩阵操作,情况会变得更糟。幸运的是,Matlab允许一个简单的解决方法,下面的代码执行相同的散射行为即上述代码,但五任何循环,所以其执行时间是大大改善(不用说,这更洁。)sctr = element(e,:);元素连接sctrVct = sctr sctr+nn ;向量分散f(sctrVct) = f(sctrVct) + fe;把一个单元刚度矩阵散成一个整体刚度矩阵以下线需要诀窍。K(sctrVct,sctrVct) = K(sctrVct,sctrVct) + ke;这个数组索引的Matlab简练的有点让人困惑,但是如果一个人花点时间适应它,它就会变得非常自然并且有用的。5.3 执行本质边界条件在最后的问题求解线性代数系统元素方程的执行本质边界条件。通常这包括修改系统。 Kd = f （19） dn = _ dn（20）通常这包括修改系统,其本质边界条件成就的同时保留了原方程在自由度元素表现。在(20)下标n指的是指数的矢量d而不是一个节点数目。一个简单的方法来执行(20)将会修nth改排K矩阵,这样去,N和K的尺寸设置。这减少nth (19)方程(20)。不幸的是,这种破坏对称性的K, 对于许多有效率线性解决者，是一个非常重要的特性。通过调整Knth栏如下：我们可以使系统对称的。当然,这将改变每个方程(19),除非我们修改力向量f。如果我们写的修正方程 (19)可以看出,我们有相同的线性方程组在(19),只是用内力约束自由度。这个程序在Matlab环境下如下：f = f - K(:,fixedDofs)*fixedDofValues;K(:,fixedDofs) = 0;K(fixedDofs,:) = 0;K(fixedDofs,fixedDofs) = bcwt*speye(length(fixedDofs);f(fixedDofs) = bcwt*fixedDofValues;在固定自由度是一个向量的i在d是固定的,固定的自由度观念是一个向量的价值观,固定自由度被指派来和bcwt称量因素是保持制约的刚度矩阵(通常是bcwt =微量(K)= N)。6、希望这个极端编程简单概述有限元方法，用Matlab有限元法已经帮助解决理论和实际的理论的区别，然后坐下来,写一些人自身的有限的元素的代码。在附录的例子应该看着,而且我建议试着写一个简单的一维或二维有限元模型。从起步到代码的方法真的固化。实例可以作为一个参考减少错误。祝你好运!A Matlab程序实例的安装所有的功能,需要运行实例程序以及实例能在/jfc795/Matlab/找到。我相信,以下文件是必要的,但是如果一个系统得到一个运行错误的职能没有找到机会,以至于都忘了列出这里却是之一，Matlab的目录在上述网站中可以找到。MeshGenerationsquare节点的数组。m:以一个数组的形式返回节点2MeshGenerationmake。m:以一个数组的形式生成元素节点.MeshGenerationmsh2mlab。m:看在Gmsh领域MeshGenerationplot。m:展示了昨天元素网格PostProcessingplot。m:展示了有限元领域正交。m:返回各种正交规则拉格朗日的基础。m:返回形函数、梯度的形态在父坐标系统功能为各种各样的元素有很多额外的文件,一个人可能会发现有用的兴趣个人可以探索这些有自己的。这些文件也应该被复制目录，包含例子脚本文件或目录在Matlab的搜索路径。B，例如:梁弯曲问题第一个例子程序解决静态线性弹性梁的弯曲。配置的问题如图3确切的解决这个问题是这样的：这个问题可以运行三元素类型;三个节点三角形元素,一个四节点四边形元素和九个节点四边形元素。同样,一个人可以选择平面应变或平面应力之间的假设。% beam.m%解决了一个二维线弹性梁问题(平面应力或应变)%与几个元素类型% 用边界条件下% u_x = 0 at (0,0), (0,-c) and (0,c)% u_y = 0 at (0,0)% t_x = y along the edge x=0% t_y = P*(x2-c2) along the edge x=L%这个文件和matlab文件可以在一些网址中被发现% /jfc795/Matlab%杰克Chessa%西北大学clearcolordef blackstate = 0;tic;dispdisp(* S T A R T I N G R U N *)dispdisp(num2str(toc), START)% 材料的性能E0 = 10e7; % 杨氏模量nu0 = 0.30; % 泊松比% 光速特性%梁的长度% t与外层纤维束的距离%网络性能elemType = Q9; %元素类型使用在有限元仿真计算 %节点恒应变三角单元, Q4 is for%一个四节点四边形元素%节点四边形元素numy = 4; % t的元素个数x-direction(梁长度)numx = 18; %如果设置1策划了初始网格(来确定一下正确性)% TIP LOADP = -1; % the peak magnitude of the traction at the right edge% STRESS ASSUMPTIONstressState=PLANE_STRESS; %设置 PLANE_STRAIN 或者PLANE_STRESS% P R E - P R O C E S S I N G *I0=2*c3/3; %第二个极地惯性矩的梁的截面% 弹性矩阵计算if ( strcmp(stressState,PLANE_STRESS) ) % Plane Strain caseC=E0/(1-nu02)* 1 nu0 0;nu0 1 0;0 0 (1-nu0)/2 ;else % Plane Strain caseC=E0/(1+nu0)/(1-2*nu0)* 1-nu0 nu0 0;nu0 1-nu0 0;0 0 1/2-nu0 ;end% 生成有限元网格%在这里,我们的gnerate有限元网格(用话分子)%“我不想再过多的细节如何使用这些功能 %一个人能有兴趣”类型-帮助函数名称”在matlab comand%来更多地了解它。%数据结构的folowing用于描述有限元% 离散化:% , i是一个矩阵的节点坐标.e. node(I,j) - x_Ij% element - is a matrix of element connectivities, i.e. the connectivity% of element e is given by element(e,:) - n1 n2 n3 .;%应用边界条件的边界描述是必要的%我们用一个独立完成这有限元离散化每个的边界%为一个二维离散化问题的边界是一组1 D元素.% rightEdge - a element connectivity matrix for the right edge% leftEdge - Ill give you three guesses%这些连接matricies参考节点编号的定义协调矩阵节点中。”disp(num2str(toc), GENERATING MESH)switch elemTypecase Q4 %我们在这里产生营运的滤网元素nnx=numx+1;nny=numy+1;node=square_node_array(0 -c,L -c,L c,0 c,nnx,nny);inc_u=1;inc_v=nnx;node_pattern= 1 2 nnx+2 nnx+1 ;element=make_elem(node_pattern,numx,numy,inc_u,inc_v);case Q9 %在这里我们产生一个mehs九方的元素nnx=2*numx+1;nny=2*numy+1;node=square_node_array(0 -c,L -c,L c,0 c,nnx,nny);inc_u=2;inc_v=2*nnx;node_pattern= 1 3 2*nnx+3 2*nnx+1 2 nnx+3 2*nnx+2 nnx+1 nnx+2 ;element=make_elem(node_pattern,numx,numy,inc_u,inc_v);%否则T3 %最后但并非不重要T3的元素nnx=numx+1;nny=numy+1;node=square_node_array(0 -c,L -c,L c,0 c,nnx,nny);node_pattern1= 1 2 nnx+1 ;node_pattern2= 2 nnx+2 nnx+1 ;inc_u=1;inc_v=nnx;element=make_elem(node_pattern1,numx,numy,inc_u,inc_v);make_elem(node_pattern2,numx,numy,inc_u,inc_v) ;end% 界定%在这里我们定位边界离散uln=nnx*(nny-1)+1; % upper left node numberurn=nnx*nny; % upper right node numberlrn=nnx; % lower right node numberlln=1; % lower left node numbercln=nnx*(nny-1)/2+1; % node number at (0,0)switch elemTypecase Q9rightEdge= lrn:2*nnx:(uln-1); (lrn+2*nnx):2*nnx:urn; (lrn+nnx):2*nnx:urn ;leftEdge = uln:-2*nnx:(lrn+1); (uln-2*nnx):-2*nnx:1; (uln-nnx):-2*nnx:1 ;edgeElemType=L3;otherwise % same discretizations for Q4 and T3 meshesrightEdge= lrn:nnx:(uln-1); (lrn+nnx):nnx:urn ;leftEdge = uln:-nnx:(lrn+1); (uln-