数值分析方法报告.docxVIP

变分法与有限元课程大作业A task ofVariational Calculus and FEM学　号：　SX1501101　　　姓 名：　刘畅　　　南京航空航天大学2016年4月14日目录1.问题描述考虑一个二维的拉普拉斯方程。图1中的三角形区域D内，u满足如下方程：S1 S3 S2 图1 区域D（1）（2）（3）（4）采用有限元法对区域D进行网格划分，并求出各节点的位移。2.有限单元法基于MATLAB采用有限元法编写程序研究上述问题，研究思路如下：图2 研究流程图2.1网格划分2.1.1网格节点编号节点编号安装按照自上而下、自左而右的规则进行(如图3)。取直角边被等分的数目n为变量，那么单元总数为n2。图3 4单元和16单元示意图kjj按上述规则划分后便只有两种类型的网格： （2）（1）iik图3 两种类型网格每个网格节点的编号顺序取逆时针。每一行网格中，（1）型网格从左到右为奇数，（2）型网格为偶数。网格节点编号的求解思路为：（1）先确定网格所处行数；（2）再确定网格在所处行数位于第几个；（3）根据其所处个数的奇偶按照确定网格类型并按逆时针得到网格节点编号。上述思路在MATLAB中采用循环语句实现：NE=n^2;（网格单元总数）ME=zeros(NE,3);（预先一个网格节点编号矩阵）for j=1:NEfor k=1:(n)（网格行数）if j-(k-1)^2-1=0j-k^2=0（确定网格行数） xx=mod((j-(k-1)^2),2);（确定奇偶）if xx==1 ME(j,:)=[(k+1)*k/2+1+(j-(k-1)^2+1)/2 (k-1)*k/2+(j-(k-1)^2+1)/2 (k+1)*k/2+(j-(k-1)^2+1)/2];（奇数）else ME(j,:)=[(k-1)*k/2+1+(j-(k-1)^2)/2-1 (k+1)*k/2+1+(j-(k-1)^2)/2 (k-1)*k/2+1+(j-(k-1)^2)/2];（偶数）endendendend2.1.2 节点坐标取直角边被等分的数目n为变量，那么节点总数为。确定节点坐标的主要思路为：（1）先确定节点处于第几行；（2）在以该节点编号减去该行起始节点编号得到该节点所处列数；（3）根据相似三角形求出节点坐标；（4）形成节点坐标列阵。上述思路在MATLAB中采用循环语句实现：NP=(n+1)*(n+2)/2;（节点总数）MP=zeros(NP,2);（预先一个节点坐标矩阵）fori=1:NPfor k=1:(n+1)（行数）ifi-(k-1)*k/2-1=0i-(k+1)*k/2=0（确定节点的行数）MP(i,:)=[(i-(k-1)*k/2-1)*l/n (n-(k-1))/n*l];（确定列数同时求出节点坐标）endendend2.2 刚度矩阵单元刚度矩阵求解的主要思路为：（1）提取单元的节点坐标；（2）选取插值函数；（3）求出单元刚度矩阵。单元刚度矩阵求出后再组集即得到总体刚度矩阵。2.3约束处理总体刚度矩阵是由单元刚度矩阵集合而成，是奇异矩阵，不能求逆，本文采用置大数法进行处理。其主要思想是：采用较大数例如：1030,对矩阵相应主元素进行替换，以便在不影响计算结果前提下使得奇异矩阵能够进行求逆运算。2.4计算结果考虑1节点处的位移，由于计算的结果和网格划分的密数目有一定的关系，变更单元数得到不同网格数目下1节点位移的计算结果，得到图4。图4 1节点位移计算结果随网格数变化从图4中我们能看出：（1）随着网格数目的增加，1节点位移逐渐增加最终趋于稳定；（2）1节点位移在网格数目增加下趋于2.7，说明1节点位移约为2.7。3.验证为验证上述程序计算结果的正确性，采用MATLAB的PDE ToolBox模块求解前述问题。主要过程如下：a.建模b.网格划分c.计算结果图5 PDE ToolBox计算过程及结果从图5（c）中能看出PDE ToolBox模块计算结果高于2.5且接近2.7，和前述程序的计算结果吻合，说明了前述所编程序的准确性。附录程序名为“naee.m”，变量所存mat文件名为：“naae.mat”。主要变量如下：n：直角边被等分数目；l:直角边长度；q:均布载荷；NF：每个节点自由度数；NP:节点总数；MP：节点坐标矩阵；NE:单元总数；ME：单元节点编号矩阵；NR:位移被约束个数；NRR：被约束位移列阵；NDF:每个单元被约束节点自由度数目;Ke:总体刚度矩阵；Ke:单元刚度矩阵；NL：载荷组数；LL:载荷列阵；P:节点力列阵；ui:节点位移列阵。下面给出程序操