数值分析方法报告

《数值分析方法报告》由会员上传分享，免费在线阅读，更多相关内容在工程资料-天天文库。

1、AtaskofVariationalCalculusandFEM变分法与有限元课程大作业学号：SX2501101姓名：刘畅南京航空大学2016年4月14日目录1•问题描述12•有限单元法12.1网格划分22.1丄网格节点编号22.1.2节点处标32.2刚度矩阵32.3约朿处理32.4计算结果33•验证4附录51•问题描述考虑一个二维的拉普拉斯方程。图1中的三角形区域D内，u满足如卞方程:图1区域D=o*$〃(1)/70z(T)=1Xgsi(2)“(jv)=0xes2(3)nt/(%)=0g5*3(4)采用冇限元法对区域D进行网格划分，并求出各节点的位移。2•有限单元法基于MATLAB采用有限

2、元法编写程序研究上述问题，研究思路如下:图2研究流程图2.1网格划分2.1.1网格节点编号节点编号安装按照白上而卜、自左而右的规则进行(如图3)。取直角边被等分的数目n为变量，那么单元总数为图34单元和16单元示意图按上述规则划分后便只冇两种类型的网格：每个网格节点的编号顺序取逆时针。每一行网格屮，(1)型网格从左到右为奇数，(2)型网格为偶数。网格节点编号的求解思路为：(1)先确定网格所处行数；(2)再确定网格在所处行数位于第几个；(3)根据具所处个数的奇他按照确定网格类型并按逆时针得到网格节点编号。上述思路在MATLAB中采用循环语句实现：NE=nA2;(网格单元总数)ME=zeros(N

3、E,3);(预先一个网格节点编号矩阵)forj=l:NEfork=l:(n)(网格行数)ifj-(k-1)A2-l>=0&j-kA2<=0(确定网格行数)xx=mod((j-(k-l)A2),2);(确定奇偶)ifxx==lME(j,:)=[(k+1)*k/2+l+(j-(k-1)A2+l)/2(k-1)*k/2+(j-(k-1)A2+l)/2(k+1)*k/2+(j-(k-1)A2+l)/2];(奇数)elseME(jz:)=[(k-1)*k/2+l+(j-(k-1)A2)/2-1(k+1)*k/2+l+(j-(k-1)A2)/2(k-1)*k/2+l+(j-(k-1)A2)/2];(偶数)

4、endendendend2.1.2节点坐标取直角边被等分的数冃n为变量，那么节点总数为（"+1）3+2）。确定节点坐标的主要2思路为：（1）先确定节点处于笫儿行；（2）在以该节点编号减去该行起始节点编号得到该节点所处列数；（3）根据和似三角形求出节点坐标；（4）形成节点朋标列阵。上述思路在MATLAB屮采用循环语句实现：NP=（n+l）*（n+2）/2;（节点总数）MP=zeros（NPZ2）;（预先一个节点处标矩阵）fori=l:NPfork=l:（n+1）（彳亍数）ifi-（k-1）*k/2-l>=0&i-（k+1）*k/2<=0（确定节点的行数）MP（iz:）=[（i-（k-1）*k/2

5、-l）*l/n（n-（k-1））/n*l];（确定列数同时求出节点坐标）end2.2刚度矩阵单元刚度矩阵求解的主要思路为：（1）提取单元的节点他标；（2）选取插值函数；（3）求出单元刚度矩阵。单元刚度矩阵求出后再组集即得到总体刚度矩阵。2.3约朿处理总体刚度矩阵是由单元刚度矩阵集合而成，是奇异矩阵，不能求逆，木文采用置大数法进行处理。其主要思想是：采用较人数例如：10巴対矩阵相应主元素进行替换，以便在不影响计算结果前提下使得奇异矩阵能够进行求逆运算。2.4计算结果考虑1节点处的位移，由于计算的结果和网格划分的密数ri冇一定的关系，变更单元数得到不同网格数目下1节点位移的计算结果，得到图4。图4

6、1节点位移计算结果随网格数变化从图4小我们能看出：(1)随着网格数目的增加，1节点位移逐渐增加最终趋于稳定；(2)1节点位移在网格数目增加下趋于2.7,说明1节点位移约为2.7。2.验证为验证上述程序计算结果的正确性，采用MATLAB的PDEToolBox模块求解前述问题。主要过程如下：Coloru图5PDEToolBox计算过程及结果从图5(c)中能看出PDEToolBox模块计算结果高于2.5且接近2.7,和前述程序的计算结果吻合，说明了前述所编程序的准确性。附录程序名为“naee.m”,变量所存m迫文件名为："naae.mat"。主要变量如下：n：直角边被等分数目；「直角边长度；q：均布

7、载荷；NF：每个节点白山度数；NP:节点总数；MP:节点坐标矩阵;NE:单元总数；ME：单元节点编号矩阵；NR：位移被约束个数；NRR：被约束位移列阵；NDF:每个单元被约束节点白由度数目；Ke：总体刚度矩阵;Ke:单元刚度矩阵；NL:载荷组数；LL:载荷列阵；P：节点力列阵；ui:节点位移列阵。下面给出程序操作方法：(1)更改n即能改变单元总数r?;(2)点击run即能计算。