平面四节点等参单元matlab实现
matlab4节点杆单元计算
篇下载MATLAB4节点杆单元计算在工程结构分析领域中，节点杆单元是一种常用的有限元分析方法。
它通过将结构划分成多个小单元，然后对每个小单元进行力学分析，最终得出整个结构的受力情况。
MATLAB作为一种强大的工程计算工具，被广泛应用于结构分析中。
本文将介绍如何利用MATLAB进行4节点杆单元计算，并提供相应的代码实例。
1. 理论背景在进行4节点杆单元计算之前，首先需要了解节点杆单元的基本理论。
节点杆单元是将结构划分为多个杆件，并在每个节点处考虑位移和受力。
通过分析每个杆件的受力平衡和位移关系，可以得出整个结构的受力和变形情况。
4节点杆单元是其中的一种常用的单元类型，它由4个节点和2个杆件组成，可以用来模拟各种不同形状和受力情况的结构。
2. MATLAB实现在MATLAB中，可以利用有限元分析工具箱进行4节点杆单元计算。
首先需要定义结构的几何形状和材料性质，并将其转化为有限元模型。
然后可以利用有限元分析工具箱提供的函数进行网格划分和边界条件设置。
接下来可以利用求解器进行结构的力学分析，并得出节点的位移和受力情况。
最后可以利用MATLAB的绘图工具对结果进行可视化展示。
3. 代码实例下面是一个简单的MATLAB代码实例，演示了如何利用有限元分析工具箱进行4节点杆单元计算：```matlab定义结构的几何形状和材料性质L = 1; 结构的长度A = 1; 结构的横截面积E = 1; 结构的弹性模量定义节点坐标node = [0, 0; 0, L; L, L; L, 0];定义单元节点关系element = [1, 2; 2, 3; 3, 4; 4, 1];网格划分和边界条件设置model = createpde();geometryFromEdges(model,(p)struct('p',p','e',[]),(p)ones(size(p,2 ),1));generateMesh(model);结构的力学分析structuralProperties(model,'YoungsModulus',E,'PoissonsRatio',0); structuralBC(model,'Edge',1,'Constraint','fixed');节点的位移和受力情况result = solve(model);可视化展示pdeplot(model,'XYData',result.displacement,'Deformation','on'); ```4. 结论通过以上代码实例，可以看到利用MATLAB进行4节点杆单元计算是非常简单和高效的。
四节点四边形单元悬臂梁的matlab有限元编程-概述说明以及解释
四节点四边形单元悬臂梁的matlab有限元编程-概述说明以及解释1.引言1.1 概述悬臂梁是一种常见的结构形式，在工程领域中被广泛应用。
四节点四边形单元是有限元分析中常用的元素类型，能够准确地模拟悬臂梁的受力情况。
Matlab是一种强大的数学工具，可以用来编程实现有限元分析。
本文旨在介绍如何利用Matlab进行四节点四边形单元悬臂梁的有限元编程，并对其进行分析和展望。
通过本文的研究，我们希望能够为工程实践提供一定的参考和指导，同时也为进一步的研究提供基础。
1.2 文章结构本文主要分为三个部分：引言、正文和结论。
引言部分将介绍文章的背景和目的，明确文章研究的问题和意义。
正文部分包括理论基础、Matlab有限元编程介绍和四节点四边形单元悬臂梁建模三个小节。
其中，理论基础将介绍与悬臂梁有关的理论知识，Matlab有限元编程介绍将详细介绍如何使用Matlab进行有限元分析，最后，四节点四边形单元悬臂梁建模将展示具体的悬臂梁建模过程。
结论部分将对实验结果进行分析与总结，探讨本研究的意义和潜在研究方向。
1.3 目的本文旨在利用Matlab编程实现四节点四边形单元悬臂梁的有限元分析，通过建立合适的数学模型，探索悬臂梁在受力状态下的力学特性。
具体目的包括：1. 建立悬臂梁的有限元数学模型，包括节点、单元和材料参数的设置；2. 实现悬臂梁在不同受力情况下的应力、应变、位移等力学性能的计算；3. 分析悬臂梁受力情况下的应力分布情况，探讨悬臂梁的破坏模式和极限承载能力；4. 验证Matlab编程方法的有效性和准确性，为工程实际中悬臂梁等复杂结构的有限元分析提供参考和借鉴。
通过本文的研究，旨在为工程实践提供可靠的数值计算工具和理论分析方法，为解决工程结构强度和稳定性问题提供一定的指导和参考价值。
2.正文2.1 理论基础在介绍四节点四边形单元悬臂梁的Matlab有限元编程之前，我们首先需要了解一些基本的理论知识。
悬臂梁是一种常见的结构形式，在工程领域中广泛应用于桥梁、建筑物等领域。
matlab坐标转换四参数法
matlab坐标转换四参数法1.引言1.1 概述在地理信息系统和测绘学中，坐标转换是一项重要的任务。
由于不同的坐标系统具有不同的基准和投影方式，因此需要进行坐标转换才能将一个点的坐标从一个坐标系统转换到另一个坐标系统。
本文将介绍一种常用的坐标转换方法——四参数法。
四参数法是一种简单而有效的坐标转换方法，通过使用四个参数进行坐标的平移和旋转，实现坐标的转换。
本文的目的是为读者介绍四参数法的原理、应用和优势。
通过深入理解四参数法的原理，读者将能够准确地将坐标在不同的坐标系统之间进行转换。
本文的结构如下：首先，将介绍坐标转换的背景，包括不同坐标系统的特点和应用领域。
其次，将详细介绍四参数法的原理，包括参数的意义和计算方法。
最后，将探讨四参数法在坐标转换中的应用，并对整个文章的内容进行总结。
通过阅读本文，读者将能够全面了解四参数法在坐标转换中的作用，掌握使用四参数法进行坐标转换的基本技巧和要点。
希望本文能够对地理信息系统和测绘学领域的专业人士和学生提供有益的参考和借鉴。
1.2文章结构文章结构部分的内容如下：1.2 文章结构本文分为引言、正文和结论三部分。
每个部分都包含了多个章节，以便清晰地呈现出Matlab坐标转换四参数法的相关内容。
在正文部分，我们将首先介绍坐标转换的背景，包括为什么需要进行坐标转换以及坐标转换的重要性。
然后，我们将详细解释四参数法的原理，包括如何使用四个参数来进行坐标转换，并且说明其适用性和局限性。
在结论部分，我们将探讨四参数法在坐标转换中的实际应用，包括它在地理信息系统和测量等领域中的重要性和实用性。
最后，我们将对整篇文章进行总结，并提出一些展望和未来的研究方向。
通过这种结构，读者将能够系统地了解Matlab坐标转换四参数法的相关知识和应用，同时也可以深入研究并拓展该方法的更多可能性。
1.3 目的本文的目的是介绍和讨论在Matlab中使用四参数法进行坐标转换的方法。
坐标转换是在地理信息系统（GIS）和测量工程中常用的技术，用于在不同的坐标系统或参考框架之间转换地理位置信息。
《有限元基础教程》_【MATLAB算例】4.8.1(1) 基于4节点四面体单元的空间块体分析(Tetrahedron3D4Node)
【MATLAB 算例】4.8.1(1) 基于4节点四面体单元的空间块体分析(Tetrahedron3D4Node)如图4-22所示的一个块体，在右端面上端点受集中力F 作用。
基于MATLAB 平台，计算各个节点位移、支反力以及单元的应力。
取相关参数为：10110Pa,=0.25E μ=⨯，5=110N F ⨯。
图4-22 一个空间块体的分析解答：对该问题进行有限元分析的过程如下。
（1）结构的离散化与编号将结构离散为5个4节点四面体单元，单元编号及节点编号和坐标如图4-22所示，连接关系见表4-8，节点的坐标见表4-9。
表4-8 单元连接关系单元号 节点号 1 2 3 4 51 42 6 1 43 7 6 7 5 1 6 7 84 1 4 6 7表4-9 节点的坐标节点节点坐标/mxyz 1 2 3 4 5 6 7 80 0 0 0.2 0 0 0 0.8 0 0.2 0.8 0 0 0 0.6 0.2 0 0.6 0 0.8 0.6 0.20.80.6节点位移列阵[]111222888 Tu v w u v w u v w =q (4-190)节点外载列阵34780 0 0 0 TT T T T⎡⎤=⎣⎦F F F F F(4-191)其中34785000 00110N ⎡⎤⎡⎤⎢⎥⎢⎥====⎢⎥⎢⎥⎢⎥⎢⎥-⨯⎣⎦⎣⎦F F F F约束的支反力列阵12560000TTT T T ⎡⎤=⎣⎦R R R R R（4-192其中1256112255661256 x x x x y y y y z z z z R R R R R R R R R R R R ⎡⎤⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥====⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎣⎦R R R R总的节点载荷列阵12345678 TT T T T T T T T⎡⎤=+==⎣⎦P F R R R R F F R R F F (4-193)（2）计算各单元的刚度矩阵(以国际标准单位)首先在MA TLAB 环境下，输入弹性模量E 、泊松比NU ，然后针对单元1和单元2，分别5次调用函数Tetrahedron3D4Node_Stiffness ，就可以得到单元的刚度矩阵k1(6×6) ～ k5(6×6)。
《有限元基础教程》_【MATLAB算例】4.8.1(1) 基于4节点四面体单元的空间块体分析(Tetrahedron3D4Node)
【MATLAB 算例】4.8.1(1) 基于4节点四面体单元的空间块体分析(Tetrahedron3D4Node)如图4-22所示的一个块体，在右端面上端点受集中力F 作用。
基于MATLAB 平台，计算各个节点位移、支反力以及单元的应力。
取相关参数为：10110Pa,=0.25E μ=⨯，5=110N F ⨯。
图4-22 一个空间块体的分析解答：对该问题进行有限元分析的过程如下。
（1）结构的离散化与编号将结构离散为5个4节点四面体单元，单元编号及节点编号和坐标如图4-22所示，连接关系见表4-8，节点的坐标见表4-9。
表4-8 单元连接关系单元号 节点号 1 2 3 4 51 42 6 1 43 7 6 7 5 1 6 7 84 1 4 6 7表4-9 节点的坐标节点节点坐标/mxyz 1 2 3 4 5 6 7 80 0 0 0.2 0 0 0 0.8 0 0.2 0.8 0 0 0 0.6 0.2 0 0.6 0 0.8 0.6 0.20.80.6节点位移列阵[]111222888 Tu v w u v w u v w =q (4-190)节点外载列阵34780 0 0 0 TT T T T⎡⎤=⎣⎦F F F F F(4-191)其中34785000 00110N ⎡⎤⎡⎤⎢⎥⎢⎥====⎢⎥⎢⎥⎢⎥⎢⎥-⨯⎣⎦⎣⎦F F F F约束的支反力列阵12560000TTT T T ⎡⎤=⎣⎦R R R R R（4-192其中1256112255661256 x x x x y y y y z z z z R R R R R R R R R R R R ⎡⎤⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥====⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎣⎦R R R R总的节点载荷列阵12345678 TT T T T T T T T⎡⎤=+==⎣⎦P F R R R R F F R R F F (4-193)（2）计算各单元的刚度矩阵(以国际标准单位)首先在MA TLAB 环境下，输入弹性模量E 、泊松比NU ，然后针对单元1和单元2，分别5次调用函数Tetrahedron3D4Node_Stiffness ，就可以得到单元的刚度矩阵k1(6×6) ～ k5(6×6)。
基于MATLAB的平面四连杆机构运动仿真.
图2~4分别为θ4角速度、点C的速度变化曲线。
4结论
本文在复数向量坐标系中推导了四连杆机构运动方程,并应用MATLAB软件进行了连杆机构运动数值仿真。从计算结果可以看出,该方法可以方便快捷地得到连杆运动参数,能够有效提高分析效率和计算精度,可进一步推广到多连杆机构设计及优化计算中。
文章编号:1009-9492(201104-0051-02
引言
四连杆机构因其结构灵活、能够传递动力并有效地实现预定动作,在很多领域得到了广泛应用
[1]
。进行连杆机
构运动分析,传统方法主要是图解法或分析法[2]
,无论设
计精度还是设计效率都相对低下,无法满足现代机械高速高精度的要求。随着计算机技术的飞速发展,特别是以
面四杆机构[J ].机械制造, 2002,
(3:26-28.
[3]周进雄,张陵.机构动态仿真[M ].西安:西安交通大学出
版社, 2002.
[4]李娟玲,张建峰.基于C语言的平面连杆机构的运动分析
[J ].机械研究与应用, 2006, 19(5:117-120.
[5]宋兆基. MATLAB6.5在科学计算中的应用[M ].北京:清
华大学出版社, 2005.
[6]王正林.精通MATLAB科学计算[M ].北京:电子工业出
版社, 2009.
[7]曹惟庆.机构设计[M ].北京:机械工业出版社, 2004. [8]李洪涛,徐巍华.基于MATLAB软件对抽油机连杆运动规律
的仿真研究[J ].机械工程师, 2009(5:99-101.
参考文献:
[1]孙桓,陈作模.机械原理[M ].北京:高等教育出版社,
2006.
四边形四节点等参元matlab程序
悬臂钢梁，尺寸如图一所示；v=0.3。
h=1，E=2.1e11.图一悬臂钢梁图二单元划分与结点编号附录：%---------------四边形四节点等参元matlab计算程序----------------------------% 2013年% 13级建筑与土木工程Bruce% B15-405% 本程序只输出了节点位移、单元中心点的应力%******************************************************************* % 变量说明% E v h% 弹性模量泊松比厚度% NPOIN NELEM NVFIX NNODE NFPOIN% 总结点数, 单元数, 约束结点个数, 单元节点数,受力结点数% COORD LNODS% 结构节点整体坐标数组, 单元定义数组,% FPOIN FORCE FIXED% 结点力数组，总体荷载向量, 约束信息数组% HK DISP% 总体刚度矩阵,结点位移向量%******************************clear allformat short eFP1=fopen(' sjd.txt','rt'); %打开数据文件%%读入控制数据E=fscanf(FP1,'%f',1); %弹性模量v=fscanf(FP1,'%f',1); % 泊松比h=fscanf(FP1,'%f',1); %厚度NELEM=fscanf(FP1,'%d',1); %单元数NPOIN=fscanf(FP1,'%d',1); % 总结点数NNODE=fscanf(FP1,'%d',1); %单元节点数NFPOIN=fscanf(FP1,'%d',1); %受力结点数NVFIX=fscanf(FP1,'%d',1); %约束结点个数LNODS=fscanf(FP1,'%f',[NNODE,NELEM])'; % 单元定义：单元结点号（逆时针）COORD=fscanf(FP1,'%f',[2,NPOIN])'; % 结点号x,y坐标(整体坐标下) FPOIN=fscanf(FP1,'%f',[3,NFPOIN])';% 节点力：结点号、X方向力(向右正)，Y方向力(向上正)FIXED=fscanf(FP1,'%d',[3,NVFIX])';%约束信息数组(n,3) n:受约束节点数目, (n,1):约束点号%(n,2)与(n,3)分别为约束点x方向和y方向的约束情况,受约束为1否则为0% 刚度矩阵的生成%计算刚度矩阵，并对约束条件进行处理Ke=zeros(2*NNODE,2*NNODE); % 单元刚度矩阵并清零HK=zeros(2*NPOIN,2*NPOIN); % 张成总刚矩阵并清零%调用子程序生成单元刚度矩阵for m=1:NELEM %m为单元号Ke=K(E,v,h,...COORD(LNODS(m,1),1),COORD(LNODS(m,1),2),...COORD(LNODS(m,2),1),COORD(LNODS(m,2),2),...COORD(LNODS(m,3),1),COORD(LNODS(m,3),2),...COORD(LNODS(m,4),1),COORD(LNODS(m,4),2)); %调用单元刚度矩阵a=LNODS(m,:); %临时向量,用来记录当前单元的节点编号%对总刚度矩阵的处理for j=1:4for k=1:4HK((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)=HK((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+...Ke(j*2-1:j*2,k*2-1:k*2);endendend %—————————————————————————————————%对荷载向量进行处理FORCE=zeros(2*NPOIN,1); % 张成总荷载向量并清零for i=1:NFPOINb1=FPOIN(i,1)*2-1;b2=FPOIN(i,1)*2; %FPION(i,1)为作用点FORCE(b1)=FPOIN(i,2); %FPION(i,2)为x方向的节点力FORCE(b2)=FPOIN(i,3); %FPION(i,3)为y方向的节点力end %—————————————————————————————————%将约束信息加入总刚,总荷载for i=1:NVFIXif FIXED(i,2)==1c1=2*FIXED(i,1)-1;HK(c1,:)=0; %将一约束序号处的总刚列向量清0HK(:,c1)=0; %将一约束序号处的总刚行向量清0HK(c1,c1)=1; %将行列交叉处的元素置为1FORCE(c1)=0;endif FIXED(i,3)==1c2=2*FIXED(i,1);HK(c2,:)=0;HK(:,c2)=0;HK(c2,c2)=1;FORCE(c2)=0;endendDISP=HK\FORCE %计算节点位移向量%———————————求解单元应力————————————————stress=zeros(3,NELEM);for m=1:NELEMu(1:8)=0;d=LNODS(m,:); %临时向量,用来记录当前单元的节点编号for i=1:NNODEu(i*2-1:i*2)=DISP(d(i)*2-1:d(i)*2);%从总位移向量中取出当前单元的节点位移endD=(E/(1-v*v))*[1 v 0;v 1 0;0 0 (1-v)/2];%弹性矩阵%形成应变矩阵BMBM=zeros(3,8);for i=1:NNODEJ=Jacobi(COORD(LNODS(m,1),1),COORD(LNODS(m,1),2),...COORD(LNODS(m,2),1),COORD(LNODS(m,2),2),...COORD(LNODS(m,3),1),COORD(LNODS(m,3),2),...COORD(LNODS(m,4),1),COORD(LNODS(m,4),2),0,0);[N_s,N_t]=DHS(0,0);B1i=J(2,2)*N_s(i)-J(1,2)*N_t(i);B2i=-J(2,1)*N_s(i)+J(1,1)*N_t(i);BM(1:3,2*i-1:2*i)=[B1i 0;0 B2i;B2i B1i]/det(J);endstressm=D*BM*u';stress(:,m)=stressm;endstress %输出应力function Ke=K(E,v,h,x1,y1,x2,y2,x3,y3,x4,y4)%=========单元刚度矩阵=============== % E 弹性模量% v 泊松比% h 厚度% x1,y1,x2,y2,x3,y3,x4,y4 为4个角结点的坐标%矩阵尺寸为8 x 8Ke=zeros(8,8);D=(E/(1-v*v))*[1 v 0;v 1 0;0 0 (1-v)/2];%弹性矩阵%高斯积分采用3 x 3 个积分点W1=5/9;W2=8/9;W3=5/9; %加权系数W=[W1 W2 W3];r=15^(1/2)/5;x=[-r 0 r];%积分点for i=1:3for j=1:3B=eleB(x1,y1,x2,y2,x3,y3,x4,y4,x(i),x(j));J=Jacobi(x1,y1,x2,y2,x3,y3,x4,y4,x(i),x(j));Ke=Ke+W(i)*W(j)*B'*D*B*det(J)*h;endendfunction B=eleB(x1,y1,x2,y2,x3,y3,x4,y4,s,t)%调用导函数[N_s,N_t]=DHS(s,t);%求Jacobi矩阵J=Jacobi(x1,y1,x2,y2,x3,y3,x4,y4,s,t);%求应变矩阵BB=zeros(3,8);for i=1:4B1i=J(2,2)*N_s(i)-J(1,2)*N_t(i);B2i=-J(2,1)*N_s(i)+J(1,1)*N_t(i);B(1:3,2*i-1:2*i)=[B1i 0;0 B2i;B2i B1i];endB=B/det(J);function J=Jacobi(x1,y1,x2,y2,x3,y3,x4,y4,s,t) %-------Jacobi-----------%单元坐标x=[x1 x2 x3 x4];y=[y1 y2 y3 y4];%%调用形函数对局部坐标的导数[N_s,N_t]=DHS(s,t);%求Jacobi矩阵的行列式的值x_s=0;y_s=0;x_t=0;y_t=0;for i=1:4x_s=x_s+N_s(i)*x(i);y_s=y_s+N_s(i)*y(i);x_t=x_t+N_t(i)*x(i);y_t=y_t+N_t(i)*y(i);endJ=[x_s y_s;x_t y_t];function N=shape(s,t)%ξ,ηN(1)=(1-s)*(1-t)/4;N(2)=(1+s)*(1-t)/4;N(3)=(1+s)*(1+t)/4;N(4)=(1-s)*(1+t)/4;function [N_s,N_t]=DHS(s,t)%形函数求导%ξ,ηN_s(1)=-1/4*(1-t);N_s(2)=1/4*(1-t);N_s(3)=1/4*(1+t);N_s(4)=-1/4*(1+t);N_t(1)=-1/4*(1-s);N_t(2)=-1/4*(1+s);N_t(3)=1/4*(1+s);N_t(4)=1/4*(1-s);sjd.txt 文件数据2.1E11 0.3 1 5 12 4 1 21 2 8 72 3 9 83 4 10 94 5 11 105 6 12 110.0 0.01.0 0.02.0 0.03.0 0.04.0 0.05.0 0.00.0 1.01.0 1.02.0 1.03.0 1.04.0 1.05.0 1.06 0 -100001 1 17 1 1。
最新平面四边形4结点等参有限单元法
有限元程序设计平面四边形4结点等参有限单元法程序设计1、程序功能及特点a.该程序采用四边形4节点等参单元，能解决弹性力学的平面应力应变问题。
b.前处理采用网格自动划分技术，自动生成单元及结点信息。
b.能计算受集中力、自重体力、分布面力和静水压力的作用。
c.计算结点的位移和单元中心点的应力分量及其主应力。
d.后处理采取整体应力磨平求得各个结点的应力分量。
e.算例计算结果与ANSYS计算结果比较，并给出误差分析。
f.程序采用Visual Fortran 5.0编制而成。
2、程序流程及图框图2-１程序流程图图2-２子程序框图其中，各子程序的主要功能为：INPUT――输入原始数据HUAFEN――自动网格划分，形成COOR(2,NP),X,Y的坐标值与单元信息CBAND――形成主元素序号指示矩阵MA（*）SKO――形成整体刚度矩阵[K]CONCR――计算集中力引起的等效结点荷载{R}eBODYR――计算自重体力引起的等效结点荷载{R}eFACER――计算分布面力引起的等效结点荷载{R}eDECOP――支配方程LU三角分解FOBA――LU分解直接解法中的回代过程OUTDISP――输出结点位移分量STRESS――计算单元应力分量OUTSTRE――输出单元应力分量STIF――计算单元刚度矩阵FDNX――计算形函数对整体坐标的导数TiiyNxN⎥⎦⎤⎢⎣⎡∂∂∂∂，=i1,2,3,4。
FUN8――计算形函数及雅可比矩阵[J]SFUN ――应力磨平-单元下的‘K’=NCN‘SCN――应力磨平-单元下的右端项系数‘CN‘SUMSKN――应力磨平-单元下的右端项集成到总体的‘P‘SUMSTRS――应力磨平-单元下的集成到总体的‘K‘GAUSTRSS――高斯消元求磨平后的应力3、输入数据及变量说明当程序开始运行时，按屏幕提示，键入数据文件的名字。
在运行程序之前，根据程序中INPUT需要的数据输入建立一个存放原始数据的文件，这个文件的名字为INDAT.DAT。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
计算力学报告平面四节点等参单元学生姓名：**学号：********一、问题描述及分析在无限大平面内有一个小圆孔。
孔内有一集中力p，试求用有限元法编程和用ANSYS软件求出各点应力分量和位移分量，并比较二者结果。
根据圣维南原理建立半径为10mm的大圆，设小圆孔的半径a=0.5mm，在远离大圆边界的地方模型是比较精确的。
由于作用在小圆孔上的力引起的位移随距离的衰减非常快，所以可以把大圆边界条件设为位移为零。
二、有限元划分描述在划分单元时，单元数量比较多，于是我采取了使用ansys软件建模自动划分单元网格的方法。
具体操作如下：打开ansys，在单元类型中选择solid->Quad 4 node 182单元；建立类半径为0.5外半径为10的圆环；使用mashtool中的智能划分和将单元退化成三角形单元；使用工具栏中List中的Nodes和Elements 选项将节点和单元数据导出并导入Excle中，总共得到了207个单元和229个节点。
如下图：图1三、有限元程序及求解程序求解使用了matlab语言。
具体如下：程序：clcclearE=2e11; %弹性模量NU=0.3; %泊松比t=0.1; %厚度X=xlsread('D:\data','nodes'); %读取节点坐标elem=xlsread('D:\data','elements'); %读取单元编号w=[1,2,3,4,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28]; %有位移约束的节点n=size(X,1); %节点数m=size(elem,1); %单元数K=zeros(2*n); %初始总体刚度矩阵for i=1:msyms Ks Et x y I1 I2 a b B; %定义可能存在的变量e=[1,1;-1,1;-1,-1;1,-1];for j=1:4 %形函数N(j)=0.25*(1+e(j,1)*Ks)*(1+e(j,2)*Et);endx=0;y=0;for j=1:4 %标准母单元映射到真实单元x=x+N(j)*X(elem(i,j),1);y=y+N(j)*X(elem(i,j),2);endJ1=jacobian([x;y],[Ks;Et]); %雅克比矩阵及其转置J=J1';for j=1:4I1=diff(N(j),Ks); %形函数分别对Ks和Et的偏导数I2=diff(N(j),Et);C=(J^-1)*[I1;I2];a=C(1); %形函数对x，y的偏导数b=C(2);B(1,2*j-1)=a; %组成B阵B(1,2*j)=0;B(2,2*j-1)=0;B(2,2*j)=b;B(3,2*j-1)=b;B(3,2*j)=a;endD=(E/(1-NU*NU))*[1,NU,0;NU,1,0;0,0,(1-NU)/2]; %D阵k=zeros(8,8);Kss=[-0.906179,-0.538469,0,0.538469,0.906179]; %5*5高斯积分点Ett=[-0.906179,-0.538469,0,0.538469,0.906179];H=[0.236926,0.478628,0.568888,0.478628,0.236926];%高斯积分权系数for j=1:5 %高斯积分求单元刚度阵for l=1:5Ks=Kss(j);Et=Ett(l);B=subs(B);J=subs(J);k=k+H(j)*H(l)*B'*D*B*det(J);endendG=zeros(8,2*n); %初始总刚变换矩阵G(1,2*elem(i,1)-1)=1; %总刚变换矩阵G(2,2*elem(i,1))=1;G(3,2*elem(i,2)-1)=1;G(4,2*elem(i,2))=1;G(5,2*elem(i,3)-1)=1;G(6,2*elem(i,3))=1;G(7,2*elem(i,4)-1)=1;G(8,2*elem(i,4))=1;K=K+G'*k*G; %总体刚度矩阵合成endKK=K;b=size(w,1);for i=1:bK(2*w(i)-1,2*w(i)-1)=1e20;K(2*w(i),2*w(i))=1e20;end %置大数法f=zeros(2*n,1); %初始载荷矩阵f(10)=-10e3; %加载荷10kNU=K\f; %节点位移for i=1:m %将每个单元各个节点位移集合u(:,i)=[U(2*elem(i,1)-1);U(2*elem(i,1));U(2*elem(i,2)-1);U(2*elem(i,2));U(2*ele m(i,3)-1);U(2*elem(i,3));U(2*elem(i,4)-1);U(2*elem(i,4))];endfor i=1:m %求单元应力syms Ks Et x y I1 I2 a b B;e=[1,1;-1,1;-1,-1;1,-1];for j=1:4N(j)=0.25*(1+e(j,1)*Ks)*(1+e(j,2)*Et);endx=0;y=0;for j=1:4x=x+N(j)*X(elem(i,j),1);y=y+N(j)*X(elem(i,j),2);endJ1=jacobian([x;y],[Ks;Et]);J=J1';for j=1:4I1=diff(N(j),Ks);I2=diff(N(j),Et);C=(J^-1)*[I1;I2];a=C(1);b=C(2);B(1,2*j-1)=a;B(1,2*j)=0;B(2,2*j-1)=0;B(2,2*j)=b;B(3,2*j-1)=b;B(3,2*j)=a; %以上同前面部分为得到B阵endD=(E/(1-NU*NU))*[1,NU,0;NU,1,0;0,0,(1-NU)/2];w=D*B*u(:,i);w1=subs(w,{Ks,Et},{1,1}); %求单元上各节点的应力sigma1(:,i)=double(w1);w2=subs(w,{Ks,Et},{-1,1});sigma2(:,i)=double(w2);w3=subs(w,{Ks,Et},{-1,-1});sigma3(:,i)=double(w3);w4=subs(w,{Ks,Et},{1,-1});sigma4(:,i)=double(w4);endc=[24,29,47,58,78,79,137,149,160,166,186]'; %如截图选取圆半径方向的节点号d=[166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184 ,185]';%圆周方向选择的节点号e=size(c,1);f=size(d,1);sigmar=zeros(3,e); %r相同角度不同节点应力矩阵sigmat=zeros(3,f); %角度不同r不同节点应力矩阵msigmar=zeros(1,e); %半径方向节点处的mises应力msigmat=zeros(1,f); %圆周方向节点处的mises应力for i=1:e %绕节点平均g=0;for j=1:m %判断节点在单元的那个位置并加上相应的应力值if elem(j,1)-c(i)==0g=g+1;sigmar(:,i)=sigmar(:,i)+sigma1(:,j);endif elem(j,2)-c(i)==0g=g+1;sigmar(:,i)=sigmar(:,i)+sigma2(:,j);endif elem(j,3)-c(i)==0g=g+1;sigmar(:,i)=sigmar(:,i)+sigma3(:,j);endif elem(j,4)-c(i)==0g=g+1;sigmar(:,i)=sigmar(:,i)+sigma4(:,j);endendsigmar(:,i)=sigmar(:,i)/g; %求应力绕节点平均msigmar(:,i)=(0.5*((sigmar(1,i)-sigmar(2,i))^2+sigmar(1,i)^2+sigmar(2,i)^2+6*(si gmar(3,i))^2))^0.5; %求节点处的mises应力endmsigmar %mises应力for i=1:f %同上g=0;for j=1:mif elem(j,1)-d(i)==0g=g+1;sigmat(:,i)=sigmat(:,i)+sigma1(:,j);endif elem(j,2)-d(i)==0g=g+1;sigmat(:,i)=sigmat(:,i)+sigma2(:,j);endif elem(j,3)-d(i)==0g=g+1;sigmat(:,i)=sigmat(:,i)+sigma3(:,j);endif elem(j,4)-d(i)==0g=g+1;sigmat(:,i)=sigmat(:,i)+sigma4(:,j);endendsigmat(:,i)=sigmat(:,i)/g;msigmat(:,i)=(0.5*((sigmat(1,i)-sigmat(2,i))^2+sigmat(1,i)^2+sigmat(2,i)^2+6*(si gmat(3,i))^2))^0.5;endmsigmat四、计算结果：1.ANSYS软件计算结果计算结果分别罗列圆周方向的单元和半径方向的单元位移和应力。
选取半径方向的单元编号为：c=[24,29,47,58,78,79,137,149,160,166,186]';选取圆周方向上的单元编号为：d=[166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181, 182,183,184,185]';其在图中的位置为图1中红线标注：图2在处理ANSYS计算结果时，我导出了节点x，y方向的正应力和剪应力，将其导入excle中并去掉无用项如图2所示，并编写matlab 程序将选取的点的mises应力求出来。
图2，求节点mises应力的程序以及计算结果如下：图3求选取点的mises应力程序：clcA=xlsread('D:\data','ansys');c=[24,29,47,58,78,79,137,149,160,166,186]';d=[166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,18 3,184,185]';e=size(c,1)f=size(d,1);for i=1:eansysr(i)=(0.5*((A(c(i),1)-A(c(i),2))^2+A(c(i),1)^2+A(c(i),2)^2+6*(A(c(i),3))^2))^0.5;endfor i=1:fansyst(i)=(0.5*((A(d(i),1)-A(d(i),2))^2+A(d(i),1)^2+A(d(i),2)^2+6*(A(d(i),3))^2))^0.5;endansysransyst计算结果如下：ansysr =1.0e+004 *Columns 1 through 90.0388 2.2423 0.0969 0.0439 0.0521 0.0702 0.1100 0.1444 0.2330Columns 10 through 110.4766 0.5914ansyst =1.0e+003 *Columns 1 through 94.7655 2.3504 0.9826 0.8454 0.6197 2.0322 0.7001 1.0187 1.2914Columns 10 through 181.3306 1.2516 1.0667 1.0775 0.6165 5.4094 4.4263 1.2137 1.3410Columns 19 through 200.5689 0.77002.等参单元编程计算结果msigmar =1.0e+004 *Columns 1 through 90.0365 3.4201 0.0458 0.0415 0.0448 0.0673 0.1396 0.1076 0.1859Columns 10 through 110.4403 3.2183msigmat =1.0e+004 *Columns 1 through 90.4403 0.3125 0.2060 0.1686 0.2256 0.6079 0.2483 0.2184 0.1692Columns 10 through 180.1140 0.2197 0.1035 0.1069 0.0854 1.5036 0.4233 0.1903 0.1553Columns 19 through 200.2950 0.2027五、讨论比较对所选取的节点的mises应力列表做对比：1.沿半径方向（排列顺序为节点号从小到大）表1节点号24 29 47 58 78 79 137 Ansys解0.0388 2.2423 0.0969 0.0439 0.0521 0.0702 0.1100 程序解0.0365 3.4201 0.0458 0.0415 0.0448 0.0673 0.1396 节点号149 160 166 186Ansys解0.1444 0.2330 0.4766 0.5914程序解0.1076 0.1859 0.4403 3.2183（注：单位为1.0e4Pa）2沿圆周方向（排列顺序按单元从大到小）表2节点号166 167 168 169 170 171 Ansys解0.4765 0.2350 0.0983 0.0885 0.0620 0.2032程序解0.4403 0.3125 0.2060 0.1686 0.2256 0.6079节点号172 173 174 175 176 177 Ansys解0.0700 0.1019 0.1291 0.1331 0.1252 0.1067 程序解0.2483 0.2184 0.1692 0.1140 0.2197 0.1035 节点号178 179 180 181 182 183 Ansys解0.1078 0.0617 0.5409 0.4426 0.1214 0.1341 程序解0.1069 0.0854 1.5036 0.4233 0.1903 0.1553 节点号184 185Ansys解0.0569 0.0770程序解0.2950 0.2027（注：单位为1.0e4Pa）对比分析：在沿半径方向数据吻合的很好，误差基本上是很小的；在圆周方向数据基本吻合，但是在一些应力值很小的点存在较大的差别，这可能与ansys处理四节点单元节点应力与程序处理方式差异有关。