用matlab进行有限元法编程Word文档格式.docx
《用matlab进行有限元法编程Word文档格式.docx》由会员分享，可在线阅读，更多相关《用matlab进行有限元法编程Word文档格式.docx（14页珍藏版）》请在冰豆网上搜索。
建议人们当开发完一个matlab程序后，回去看看如何/是否可以消除任何循环。
通过实践，这将成为matlab编程的第二特性。
3一个典型有限元程序的组成部分
一个典型的有限元程序包含以下三个部分
1.预处理部分
2.处理部分
3.后处理部分
在预处理部分定义了问题陈述的数据和结构会被定。
这些措施包括有限元离散，材料的性能，解决参数等等。
在处理部分，有限元素对象即STI内斯矩阵，计算力矢量等被计算，边界条件被实施，系统被解决。
在后处理部分，处理部分的结果将被分析。
这里应力可能被计算并且数据进行可视化。
本文中，我们将主要关注处理部分。
许多前期和后期处理操作已经在matlab中进行编程并且在网上有参考文献；
有兴趣的人可以直接参看MATLAB脚本文件或在MATLAB命令行中键入帮助“functionname'
，以获得有关如何使用这些功能的进一步信息。
4matlab中有限元数据结构
在这里，我们讨论的是在有限元法中使用，尤其在示例代码中实现的数据结构。
人们可以想出多种储存有限元程序数据的方式，这有些随意，而我们尽力使用最灵活而且最有利于matlab的结构。
这些数据结构的设计可能取决于所使用的编程语言，但通常不会与这里列出的结构明显不同。
4.1节点坐标矩阵
由于我们进行有限元法编程，因此我们需要一种代表元素离散域的方式。
这样做，我们定义一组节点和用一定方式连接这些节点的元素。
节点的坐标存储在节点坐标矩阵中。
这是一个简单的节点坐标的矩阵（可以想象一下）。
这个矩阵的维是nn×
sdim，其中nn是节点的数量，sdim是空间维数。
所以，如果我们考虑一个节点坐标矩阵的节点，第n个节点的y坐标是节点（N2）。
图1显示了一个简单的有限元离散矩阵。
对于这个简单的网格节点坐标矩阵如下：
节点=
（1）
4.2元素连接矩阵
元素的定义存储在元素连接矩阵中。
这是一个表示节点数目的矩阵，矩阵的每一行包含一个元素的连性数。
因此，如果我们考虑连接矩阵元素，即描述了一个4节点四边形的网格,那么第36个元素是由连接向量元素（36，:
）定义的，例如其可能是[36421314]或元素连接节点36→42→13→14。
因此，图1中的简单网格的元素连接矩阵是
元素=
（2）
注意：
该元素的连接数都在逆时针排列的；
如果没有这样做，一些雅可比矩阵的元素会是负数，而会导致刚度矩阵奇异（而且显然是错误的！
）。
4.3边界的定义
在有限元法中，边界条件用于组成力矢量（自然或诺伊曼边界条件）或指定边界上（必要或狄利克雷边界条件）的未知领域的值。
在这两种情况下对于边界的定义是必要的。
实现这一点的最通用的方式是保持必要的边界有限元离散。
这个网格的维度将比问题空间维度小1（即二维边界网格对应于三维问题，一维边界网格对应于二维问题等等）。
让我们再次考虑一下图1中的简单网格。
假设我们希望边界条件适用于网格的右边缘，那么边界网格的定义由以下元素的2个节点的行元素的连接矩阵定义。
右边缘=
（3）
请注意，边界连接矩阵中的数字指同一节点坐标矩阵，就如同内部元素的连接矩阵中的数字一样。
如果我们想要一个必要的边界条件适用于一个边缘，我们需要这个边缘上的节点列表。
这可以在matlab中通过其独特的功能而很容易的实现。
nodesOnBoundary=unique（rightEdge）；
这将设置向量nodesOnBoundary等于[246]。
如果我们想要在这个边缘上通过自然边界条件组成一个力向量，我们只需循环遍历该元素并整合该边缘上的力，就像我们将域内部，即刚度矩阵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,i}→n（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中的连续部分，如下：
（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）;
Forq=1:
长度（qWt）
xi=qPt（问）;
%正交点
球形坐标x%得到在正交分11人阵容
%和导数行列式在正交点,jac
…
f_int=f_int+x^3**qWtjac（q）;
结束
5.2算子“散射”
一元算子的计算需要分散到整个球面系。
散射元素力向如一个整体力矢量如图2所示。
散射依赖在元件连接和自由度映射选择。
下面的代码执行散射如图2所示
elemConn=元素（e,:
）;
%元素连接
enn=长（elemConn）;
forI=1:
enn;
%在节点元素循环
2%loop空间维度
Ii=nn*（i-1）+sctr（I）;
%自由度
f（Ii）=f（Ii）+f（（i-1）*enn+I）;
但使用了一个嵌套的循环的（环中的环）。
这是一个更异乎寻常的行为，考虑到这样一个事实时,即它发生一个元素循环,这真的是个能减慢执行时间的程序（由数量级在许多例）。
我们将矩阵有四个嵌套的循环，刚度散射矩阵操作,情况会变得更糟。
幸运的是,Matlab允许一个简单的解决方法,下面的代码执行相同的散射行为即上述代码,但五任何循环,所以其执行时间是大大改善（不用说,这更洁。
）
sctr=element（e,:
元素连接
sctrVct=[sctrsctr+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有限元法已经帮助解决理论和实际的理论的区别，然后坐下来,写一些人自身的有限的元素的代码。
在附录的例子应该看着,而且我建议试着写一个简单的一维或二维有限元模型。
从起步到代码的方法真的固化。
实例可以作为一个参考减少错误。
祝你好运!
AMatlab程序实例的安装
所有的功能,需要运行实例程序以及实例能在http:
//www.tam.northwestern.edu/jfc795/Matlab/找到。
我相信,以下文件是必要的,但是如果一个系统得到一个运行错误的职能没有找到机会,以至于都忘了列出这里却是之一，Matlab的目录在上述网站中可以找到。
MeshGenerationsquare节点的数组。
m:
以一个数组的形式返回节点2
MeshGenerationmake。
以一个数组的形式生成元素节点.
MeshGenerationmsh2mlab。
看在Gmsh领域
MeshGenerationplot。
展示了昨天元素网格
PostProcessingplot。
展示了有限元领域
正交。
返回各种正交规则
拉格朗日的基础。
返回形函数、梯度的形态
在父坐标系统功能为各种各样的元素
有很多额外的文件,一个人可能会发现有用的兴趣个人可以探索这些有自己的。
这些文件也应该被复制目录，包含例子脚本文件或目录在Matlab的搜索路径。
B，例如:
梁弯曲问题
第一个例子程序解决静态线性弹性梁的弯曲。
配置的问题如图3确切的解决这个问题是这样的：
这个问题可以运行三元素类型;
三个节点三角形元素,一个四节点四边形元素和九个节点四边形元素。
同样,一个人可以选择平面应变或平面应力之间的假设。
%beam.m
%解决了一个二维线弹性梁问题（平面应力或应变）
%与几个元素类型
%用边界条件下
%u_x=0at（0,0）,（0,-c）and（0,c）
%u_y=0at（0,0）
%t_x=yalongtheedgex=0
%t_y=P*（x^2-c^2）alongtheedgex=L
%这个文件和matlab文件可以在一些网址中被发现
%http:
//www.tam.northwestern.edu/jfc795/Matlab
%杰克•Chessa
%西北大学
clear
colordefblack
state=0;
tic;
dispdisp（'
***STARTINGRUN***'
disp
disp（[num2str（toc）,'
START'
]）
%材料的性能
E0=10e7;
%杨氏模量
nu0=0.30;
%泊松比
%光速特性
%梁的长度
%t与外层纤维束的距离
%网络性能
elemType='
Q9'
;
%元素类型使用在有限元仿真计算
%节点恒应变三角单元,'
Q4'
isfor
%一个四节点四边形元素
%节点四边形元素
numy=4;
%t的元素个数x-direction（梁长度）
numx=18;
%如果设置1策划了初始网格（来确定一下正确性）
%TIPLOAD
P=-1;
%thepeakmagnitudeofthetractionattherightedge
%STRESSASSUMPTION
stressState='
PLANE_STRESS'
%设置'
PLANE_STRAIN'
或者"
%
PRE-PROCESSING***
I0=2*c^3/3;
%第二个极地惯性矩的梁的截面
%弹性矩阵计算
if（strcmp（stressState,'
））%PlaneStraincase
C=E0/（1-nu0^2）*[1nu00;
nu010;
00（1-nu0）/2];
else%PlaneStraincase
C=E0/（1+nu0）/（1-2*nu0）*[1-nu0nu00;
nu01-nu00;
001/2-nu0];
%生成有限元网格
%在这里,我们的gnerate有限元网格（用话分子）
%“我不想再过多的细节如何使用这些功能
%一个人能有兴趣”类型-帮助'
函数名称”在matlabcomand
%来更多地了解它。
%数据结构的folowing用于描述有限元
%离散化:
%,i是一个矩阵的节点坐标.e.node（I,j）->
x_Ij
%element-isamatrixofelementconnectivities,i.e.theconnectivity
%ofelementeisgivenby>
element（e,:
）->
[n1n2n3...];
%
%应用边界条件的边界描述是必要的
%我们用一个独立完成这有限元离散化每个的边界
%为一个二维离散化问题的边界是一组1D元素.
%rightEdge-aelementconnectivitymatrixfortherightedge
%leftEdge-I'
llgiveyouthreeguesses
%这些连接matricies参考节点编号的定义协调矩阵节点中。
”
GENERATINGMESH'
switchelemType
case'
%我们在这里产生营运的滤网元素
nnx=numx+1;
nny=numy+1;
node=square_node_array（[0-c],[L-c],[Lc],[0c],nnx,nny）;
inc_u=1;
inc_v=nnx;
node_pattern=[12nnx+2nnx+1];
element=make_elem（node_pattern,numx,numy,inc_u,inc_v）;
%在这里我们产生一个mehs九方的元素nnx=2*numx+1;
nny=2*numy+1;
inc_u=2;
inc_v=2*nnx;
node_pattern=[132*nnx+32*nnx+12nnx+32*nnx+2nnx+1nnx+2];
%否则'
T3'
%最后但并非不重要T3的元素
node_pattern1=[12nnx+1];
node_pattern2=[2nnx+2nnx+1];
element=[make_elem（node_pattern1,numx,numy,inc_u,inc_v）;
make_elem（node_pattern2,numx,numy,inc_u,inc_v）];
%界定
%在这里我们定位边界离散
uln=nnx*（nny-1）+1;
%upperleftnodenumber
urn=nnx*nny;
%upperrightnodenumber
lrn=nnx;
%lowerrightnodenumber
lln=1;
%lowerleftnodenumber
cln=nnx*（nny-1）/2+1;
%nodenumberat（0,0）
rightEdge=[lrn:
2*nnx:
（uln-1）;
（lrn+2*nnx）:
urn;
（lrn+nnx）:
urn]'
leftEdge=[uln:
-2*nnx:
（lrn+1）;
（uln-2*nnx）:
1;
（uln-nnx）:
1]'
edgeElemType='
L3'
otherwise%samediscretizationsforQ4andT3meshes
nnx:
-nnx:
L2'
%得到节点位移边界
%在这里,我们得到的结点基本界域
fixedNodeX=[ulnllncln]'
一个向量的节点被固定在数字
fixedNodeY=[clnxdirection]'
%avectorofnodenumberswhicharefixedin
%they-direction
plot_mesh（node,rightEdge,edgeElemType,'
bo-'
plot_mesh（node,leftEdge,edgeElemType,'
plot（node（fixedNodeX,1）,node（fixedNodeX,2）,'
r>
'
plot（node（fixedNodeY,1）,node（fixedNodeY,2）,'
r^'
axisoff
axis（[0L-cc]）
disp（'
（paused）'
pause
%计算元素在应力应变和应力的观点strain=B*U（sctrB）;
stress（e,q,:
）=C*strain;
end%ofelementloop
stressComp=1;
figure（fn）
clf
plot_field（node+scaleFact*[U（xs）U（ys）],element,elemType,stress（:
:
stressComp））;
holdon
plot_mesh（node+scaleFact*[U（xs）U（ys）],element,elemT
升级会员 