平面四节点等参单元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。
基于matlab的平面四连杆机构设计以及该机构的运动分析参考模板
基于matlab的平面四连杆机构设计以及该机构的运动仿真分析摘要四连杆机构因其结构方便灵活,能够传递动力并实现多种运动形式而被广泛应用于各个领域,因此对其进行运动分析具有重要的意义。
传统的分析方法主要应用几何综合法和解析综合法,几何综合法简单直观,但是精确度较低;解析法精确度较高,但是计算工作量大。
随着计算机辅助数值解法的发展,特别是MATLAB软件的引入,解析法已经得到了广泛的应用。
对于四连杆的运动分析，若应用MATLAB 则需要大量的编程，因此我们引入proe软件，我们不仅可以在此软件中建立实物图，而且还可以对其进行运动仿真并对其运动分析。
在设计四连杆时，我们利用解析综合法建立数学模型，再根据数学模型在MATLAB中编程可以求得其他杆件的长度。
针对范例中所求得的各连杆的长度，我们在proe软件中画出其三维图（如图4）并在proe软件中进行仿真分析得出CB,的角加速度的变化，从而得到CB,两接触处所受到的力是成周期性变化的，可以看出CB,两点处的疲劳断裂，我们提B,两点处极易疲劳断裂，针对C出了在设计四连杆中的一些建议。
关键字：解析法 MATLAB 软件 proe 软件 运动仿真建立用解析法设计平面四杆机构模型对于问题中所给出的连架杆AB 的三个位置与连架杆CD 的三个位置相对应，即三组对应位置为：332211,,,,,ψϕψϕψϕ，其中他们对应的值分别为： 52,45,82,90,112,135，为了便于写代数式，可作出AB 与CD 对应的关系，其图如下：图—2 AB 与CD 三个位置对应的关系通过上图我们可以通过建立平面直角坐标系并利用解析法来求解，其直角坐标系图如下：φααi θi φi图—3 平面机构直角坐标系通过建立直角坐标系OXY ，如上图所示，其中0α与0φ为AB 杆与CD 杆的初始角，各杆件的长度分别用矢量d c b a ,,,,表示，将各矢量分别在X 轴与Y 轴上投影的方程为⎩⎨⎧=++=+)sin(*)sin(*)sin(*)cos(*)cos(*)cos(*φθαφθαc b a c d b a在上述的方程中我们可以消除θ，从而可以得到α与φ之间的关系如下：)cos(2)cos(2)cos(2)(2222αφαφab ac cd b d c a +-=+-++ (1) 为便于化简以及matlab 编程我们可以令：⎪⎪⎪⎩⎪⎪⎪⎨⎧==-++=c d H a d H ac b d c a H 32222212 (2) 通过将（2）式代入（1）式中则可以化简得到如下等式： )cos()cos()cos(321αφαφH H H +-=+ （3）我们可以通过（3）式将两连架杆对应的位置带入（3）式中，我们可以得到如下方程：⎪⎩⎪⎨⎧+-=++-=++-=+)cos()cos()cos()cos()cos()cos()cos()cos()cos(333332123222211311121ϕψϕψϕψϕψϕψϕψH H H H H H H H H （4） 联立（4）方程组我们可以求得321,,H H H ，再根据（2）中的条件以及所给定的机架d 的长度，我们可以求出其它杆件的长度为：⎪⎪⎪⎩⎪⎪⎪⎨⎧-++===1222322acH d c a b H d c H d a （5）四连杆设计范例：在日常生活中，我们经常看到消防门总能自动关上，其实它是利用四连杆机构与弹簧组成的。
平面四边形四节点等参单元Fortran源程序文件
C ************************************************C * FINITE ELEMENT PROGRAM *C * FOR Two DIMENSIONAL ELASticity PROBLEM *C * WITH 4 NODE *C ************************************************PROGRAM ELASTICITYcharacter*32 dat,cchDIMENSION SK(80000),COOR(2,300),AE(4,11),MEL(5,200),& WG(4),JR(2,300),MA(600),R(600),iew(30),STRE(3,200)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)WRITE(*,*)'PLEASE ENTER INPUT FILE NAME'READ(*,'(A)')DATOPEN(4,FILE=dat,STATUS='OLD')OPEN(7,FILE='OUT',STATUS='UNKNOWN')READ(4,*)NP,NE,NM,NRWRITE(7,'(A,I6)')'NUMBER OF NODE---------------------NP=',np WRITE(7,'(A,I6)')'NUMBER OF ELEMENT------------------NE=',ne WRITE(7,'(A,I6)')'NUMBER OF MATERIAL-----------------NM=',nm WRITE(7,'(A,I6)')'NUMBER OF surporting---------------NC=',Nr CALL INPUT (JR,COOR,AE,MEL)CALL CBAND (MA,JR,MEL)DO I=1,NHSK(I)=0.0enddoCALL SK0(SK,MEL,COOR,JR,MA,AE)do I=1,NR(I)=0.0enddopause 'aaa'stopREAD(4,*)NCP,NBE,izWRITE(*,'(5i8)')NCP,NBE,izWRITE(7,'(5i8)')NCP,NBE,izIF(NCP.GT.0)CALL CONCR(NCP,R,JR)IF(NBE.GT.0) CALL BODYR(NBE,R,MEL,COOR,JR,AE)IF(iz.GT.0)thendo jj=1,izREAD (4,*)Js,nse,(WG(I),I=1,4)read(4,*)(iew(m),m=1,nse)CALL FACER(iew,NSE,R,MEL,COOR,JR,WG)enddoendifCALL DECOP (SK,MA)CALL FOBA (SK,MA,R)CALL OUTDISP(NP,R,JR)CALL STRESS (COOR,MEL,JR,AE,R,STRE)WRITE(7,'(A)')' PROGRAM SAFF HAS BEEN ENDED'WRITE(*,'(A)')' PROGRAM SAFF HAS BEEN ENDED'STOPc RETURNENDC *********************************************SUBROUTINE INPUT (JR,COOR,AE,MEL)DIMENSION JR(2,*),COOR(2,*),AE(4,*),MEL(5,*)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHDO 70 I=1,NPREAD(4,*) IP,X,YCOOR(1,IP)=XCOOR(2,IP)=Y70 CONTINUEDO 11 J=1,NEREAD(4,*)NEE,NME,(MEL(I,NEE),I=1,4)MEL(5,NEE)=NME11 CONTINUEDO 10 I=1,NPDO 10 J=1,210 JR(J,I)=1DO 20 I=1,NRREAD(4,*) IP,IX,IYJR(1,IP)=IXJR(2,IP)=IY20 CONTINUEN=0DO 30 I=1,NPDO 30 J=1,2IF (JR(J,I)) 30,30,2525 N=N+1JR(J,I)=N30 CONTINUEDO 55 J=1,NMREAD (4,*)JJ,(AE(I,JJ),I=1,4)WRITE(*,910) JJ,(AE(I,JJ),I=1,4)55 CONTINUE910 FORMAT (/20X,'MATERIAL PROPERTIES'/(3X,I5,4(1x,E8.3))) RETURNENDC **********************************************SUBROUTINE CBAND (MA,JR,MEL)DIMENSION MA(*),JR(2,*),MEL(5,*),NN(8)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHDO 65 I=1,N65 MA(I)=0DO 90 IE=1,NEDO 75 K=1,4IEK=MEL(K,IE)DO 95 M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)95 CONTINUE75 CONTINUEL=NDO 80 I=1,2*4NNI=NN(I)IF(NNI.EQ.0) GO TO 80IF(NNI.LT.L) L=NNI80 CONTINUEDO 85 M=1,2*4JP=NN(M)IF(JP.EQ.0) GO TO 85JPL=JP-L+1IF(JPL.GT.MA(JP)) MA(JP)=JPL85 CONTINUE90 CONTINUEMX=0MA(1)=1DO 10 I=2,NIF(MA(I).GT.MX) MX=MA(I)MA(I)=MA(I)+MA(I-1)10 CONTINUENH=MA(N)WRITE(7,'(A,I8)')'TOTAL DEGREES OF FREEDOM-----------N= ',N WRITE(7,'(A,I8)')'MAX-SEMI-BANDWIDTH-----------------MX=',MX WRITE(7,'(A,I8)')'TOTAL-STORAGE----------------------NH=',NH 500 FORMAT (/5X,'FREEDOM N='*,I5,3X,'SEMI-BANDWI. MX=',I5,3X,* 'STORAGE NH=',I7)RETURNENDC **********************************************SUBROUTINE SK0(SK,MEL,COOR,JR,MA,AE)DIMENSION SK(*),MEL(5,*),COOR(2,*),JR(2,*),MA(*), * AE(4,*),XYZ(2,4),iven(4)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)COMMON /CMN4/ NEE,NMECOMMON /GAUSS/ RSTG(3),H(3)H(1)=0.55560H(2)=0.88890H(3)=H(1)RSTG(1)=-0.14830RSTG(2)=0.00RSTG(3)=-RSTG(1)DO 10 IE=1,NENEE=IENME=MEL(5,IE)DO 75 K=1,4IEK=MEL(K,IE)iven(k)=IEKDO 95 M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)95 XYZ(M,K)=COOR(M,IEK)75 CONTINUECALL STIF(XYZ,AE,iven)DO 60 I=1,8DO 60 J=1,8II=NN(I)JJ=NN(J)IF ((JJ.EQ.0).OR.(II.LT.JJ)) GO TO 60JN=MA(II)-(II-JJ)SK(JN)=SK(JN)+SKE(I,J)60 CONTINUE70 CONTINUEwrite(7,1111) ((ske(i,j),j=1,8),i=1,8)1111 format(2x,8f12.2)10 CONTINUERETURNENDC *********************************************SUBROUTINE STIF(XYZ,AE,iven)DIMENSION AE(4,*),DNX(2,4),XYZ(2,*),iven(*),* RJAC(2,2)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)COMMON /CMN4/ NEE,NMECOMMON /GAUSS/ RSTG(3),H(3)DO 40 I=1,8RF(I)=0.00DO 30 J=1,8SKE(I,J)=0.0030 CONTINUE40 CONTINUEE=AE(1,NME)U=AE(2,NME)GAMA=AE(3,NME)D1=E*(1.00-U)/((1.00+U)*(1.00-2.00*U))D2=E*U/((1.00+U)*(1.00-2.00*U))D3=E*0.50/(1.00+U)DO 120 I=1,4II=2*(I-1)I1=II+1I2=II+2DO 115 J=1,4JJ=2*(J-1)J1=JJ+1J2=JJ+2DXX=0DXY=0DYX=0DYY=0DO 99 IS=1,3S=RSTG(IS)SH=H(IS)DO 98 IR=1,3R=RSTG(IR)RH=H(IR)CALL FDNX (XYZ,DNX,DET,R,S,RJAC,iven,NEE) DNIX=DNX(1,I)DNIY=DNX(2,I)DNJX=DNX(1,J)DNJY=DNX(2,J)DXX=DXX+DNIX*DNJX*DET*RH*SHDXY=DXY+DNIX*DNJY*DET*RH*SHDYX=DYX+DNIY*DNJX*DET*RH*SHDYY=DYY+DNIY*DNJY*DET*RH*SH98 CONTINUE99 CONTINUESKE(I1,J1)=DXX*D1+DYY*D3SKE(I2,J2)=DYY*D1+DXX*D3SKE(I1,J2)=DXY*D2+DYX*D3SKE(I2,J1)=DYX*D2+DXY*D3115 CONTINUE120 CONTINUERETURNENDC ********************************************* SUBROUTINE CONCR(NCP,R,JR)DIMENSION R(*),JR(2,*),XYZ(2)DO 100 I=1,NCPREAD (4,*) IP,PX,PYXYZ(1)=PXXYZ(2)=PYDO 95 J=1,2L=JR(J,IP)IF(L.EQ.0) GO TO 95R(L)=R(L)+XYZ(J)95 CONTINUE100 CONTINUERETURNENDC ********************************************** SUBROUTINE BODYR(NBE,R,MEL,COOR,JR,AE)DIMENSION R(*),MEL(5,*),COOR(2,*),JR(2,*), & AE(4,*),XYZ(2,4),iven(4)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)COMMON /GAUSS/ RSTG(3),H(3)H(1)=1.0H(2)=1.0RSTG(1)=-0.96260RSTG(2)=-RSTG(1)DO 10 IE=1,NBEDO I=1,8RF(I)=0.00ENDDOc READ(4,*)NEENEE=ieNME=MEL(5,NEE)GAMA=AE(3,NME)DO 75 K=1,4IEK=MEL(K,NEE)iven(k)=iekDO 95 M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)95 XYZ(M,K)=COOR(M,IEK)75 CONTINUEDO 99 IS=1,2S=RSTG(IS)SH=H(IS)DO 98 IR=1,2RR=RSTG(IR)RH=H(IR)CALL FUN8 (XYZ,RR,S,DET)DO 30 I=1,4J=2*IRF(J)=RF(J)-FUN(I)*RH*SH*DET*GAMA30 CONTINUE98 CONTINUE99 CONTINUECALL ASLOAD (R)10 CONTINUERETURNENDC *********************************************SUBROUTINE FACER(iew,NSE,R,MEL,COOR,JR,WG)DIMENSION R(*),MEL(5,*),COOR(2,*),JR(2,*),wg(*) * ,XYZ(2,4),iew(*),PR(2)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)COMMON /CMN4/ NEE,NMECOMMON /GAUSS/ RSTG(3),H(3)H(1)=1.0H(2)=1.0RSTG(1)=-0.96260RSTG(2)=-RSTG(1)nwf=0nnf=0ir=wg(1)+0.1if(ir.eq.1)nwf=1if(ir.eq.2)nnf=1DO 510 IE=1,NSEDO I=1,8RF(I)=0.00ENDDOnee=iew(ie)DO 575 K=1,4IEK=MEL(K,NEE)DO 595 M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)595 XYZ(M,K)=COOR(M,IEK)575 CONTINUEIF(NWF.EQ.1) thenGAMA=WG(2)Z0=WG(3)NSU=WG(4)+0.1CALL SURLOD (NSU,XYZ,PR,Z0,GAMA,1)endifIF(NNF.EQ.1) thenq=WG(2)NSU=WG(4)+0.1do j=1,2PR(J)=qenddoCALL SURLOD (NSU,XYZ,PR,Z0,GAMA,2)endifCALL ASLOAD (R)510 CONTINUERETURNENDC *********************************************SUBROUTINE SURLOD (NSU,XYZ,PR,Z0,GAMA,NSI)DIMENSION XYZ(2,*),RST(3),PR(2),KCRD(4),KFACE(2,4), & FVAL(4),NODES(2),FACT(4)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)COMMON /CMN4/ NEE,NMECOMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)COMMON /GAUSS/ RSTG(3),H(3)DATA KCRD/1,1,2,2/DATA KFACE/1, 4,* 2, 3,* 1, 2,* 4, 3/DATA FVAL/-1.00,1.00,-1.00,1.00/ FACT(1)=1.0FACT(2)=-1.0FACT(3)=-1.0FACT(4)=1.0FACTNUS=FACT(NSU)DO I=1,2J=KFACE(I,NSU)NODES(I)=JENDDOIF (NSI.EQ.1) THENDO I=1,2J=NODES(I)Z=Z0-XYZ(2,J)PR(I)=0.00IF (Z.GT.0.00) PR(I)=Z*GAMAENDDOENDIFML=KCRD(NSU)IF(ML.EQ.1)MM=2IF(ML.EQ.2)MM=1RST(ML)=FVAL(NSU)DO 70 LX=1,2RST(MM)=RSTG(LX)CALL FUN8 (XYZ,RST(1),RST(2),DET) PXYZ=0.00DO 25 I=1,2J=NODES(I)PXYZ=PXYZ+FUN(J)*PR(I)25 CONTINUEA1=XJAC(MM,2)A2=-XJAC(MM,1)30 DO 60 I=1,2J=NODES(I)K2=2*JK1=K2-1Q=PXYZ*FUN(J)*H(LX)*FACTNUSRF(K1)=RF(K1)+Q*A1RF(K2)=RF(K2)+Q*A260 CONTINUE70 CONTINUERETURNENDC ********************************************* SUBROUTINE ASLOAD (R)DIMENSION R(*)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)DO 20 I=1,8L=NN(I)IF (L.EQ.0) GO TO 20R(L)=R(L)+RF(I)20 CONTINUERETURNENDC *********************************************** SUBROUTINE DECOP (SK,MA)DIMENSION SK(*),MA(*)COMMON /CMN2/ N,MX,NHDO 50 I=2,NL=I-MA(I)+MA(I-1)+1K=I-1L1=L+1IF (L1.GT.K) GO TO 30DO 20 J=L1,KIJ=MA(I)-I+JM=J-MA(J)+MA(J-1)+1IF (L.GT.M) M=LMP=J-1IF (M.GT.MP) GO TO 20DO 10 LP=M,MPIP=MA(I)-I+LPJP=MA(J)-J+LPSK(IJ)=SK(IJ)-SK(IP)*SK(JP)10 CONTINUE20 CONTINUE30 IF (L.GT.K) GO TO 50DO 40 LP=L,KIP=MA(I)-I+LPLPP=MA(LP)SK(IP)=SK(IP)/SK(LPP)II=MA(I)SK(II)=SK(II)-SK(IP)*SK(IP)*SK(LPP)40 CONTINUE50 CONTINUERETURNENDC *************************************************************SUBROUTINE FOBA (SK,MA,R)DIMENSION SK(*),MA(*),R(*)COMMON /CMN2/ N,MX,NHDO 10 I=2,NL=I-MA(I)+MA(I-1)+1K=I-1IF (L.GT.K) GO TO 10DO 5 LP=L,KIP=MA(I)-I+LPR(I)=R(I)-SK(IP)*R(LP)5 CONTINUE10 CONTINUEDO 20 I=1,NII=MA(I)45 R(I)=R(I)/SK(II)20 CONTINUEDO 30 J1=2,NI=2+N-J1L=I-MA(I)+MA(I-1)+1K=I-1IF (L.GT.K) GO TO 30DO 25 J=L,KIJ=MA(I)-I+J55 R(J)=R(J)-SK(IJ)*R(I)25 CONTINUE30 CONTINUERETURNENDC ***************************************************************** SUBROUTINE STRESS(COOR,MEL,JR,AE,R,STRE)DIMENSION XYZ(2,4),DNX(2,4),AE(4,*),STRE(3,*),& COOR(2,*),MEL(5,*),JR(2,*),RJAC(2,2),SIG(3),& B(3,8),R(*),iven(4)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)DO 106 IE=1,NENME=MEL(5,IE)DO 300 K=1,4IEK=MEL(K,IE)DO 310 M=1,2310 XYZ(M,K)=COOR(M,IEK)DO 320 M=1,2JRR=2*(K-1)+M320 NN(JRR)=JR(M,IEK)300 CONTINUEE=AE(1,NME)U=AE(2,NME)D1=E*(1.00-U)/((1.00+U)*(1.00-2.00*U))D2=E*U/((1.00+U)*(1.00-2.00*U))D3=0.50*E/(1.00+U)SS=0.0RR=0.0CALL FDNX (XYZ,DNX,DET,RR,SS,RJAC,iven,IE) DO 30 I=1,4II=2*(I-1)J1=II+1J2=II+2BI=DNX(1,I)CI=DNX(2,I)B(1,J1)=BIB(2,J1)=0.B(3,J1)=CIB(1,J2)=0.B(2,J2)=CIB(3,J2)=BI30 CONTINUEDO 55 II=1,3SIG(II)=0.0055 CONTINUEDO 70 K=1,8NA=NN(K)IF (NA.EQ.0) GO TO 70DO 60 L=1,3SIG(L)=SIG(L)+B(L,K)*R(NA)60 CONTINUE70 CONTINUESX=D1*SIG(1)+D2*SIG(2)SY=D2*SIG(1)+D1*SIG(2)SXY=D3*SIG(3)STRE(1,IE)=SXSTRE(2,IE)=SYSTRE(3,IE)=SXY106 CONTINUECALL OUTSTRE(NE,STRE)RETURNENDC *********************************************SUBROUTINE FDNX (XYZ,DNX,DET,R,S,RJAC,iven,NEE) DIMENSION XYZ(2,*),DNX(2,*),RJAC(2,2),iven(*)COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)CALL FUN8 (XYZ,R,S,DET)IF (DET.LT.1.0E-5)THENWRITE(7,600) NEE,R,S,detWRITE(7,*) (iven(m),m=1,4)STOPENDIFREC=1.00/DETRJAC(1,1)=REC*XJAC(2,2)RJAC(2,2)=REC*XJAC(1,1)RJAC(2,1)=-REC*XJAC(2,1)RJAC(1,2)=-REC*XJAC(1,2)DO 30 K=1,4DO 20 I=1,2DNX(I,K)=0.DO 25 M=1,2DNX(I,K)=DNX(I,K)+RJAC(I,M)*PN(M,K)25 CONTINUE20 CONTINUE30 CONTINUE600 FORMAT (1X,'ERR0R*** NEGTIVE OR ZERO '* 'JACOBIAN DETERMINANT FOR '* 'ELEMENT'/'ELE.=',I5,' R=',F10.5,6X,'S=',F10.5, * 'det=',f12.5)RETURNENDC *********************************************SUBROUTINE FUN8 (XYZ,R,S,DET)DIMENSION XYZ(2,*),XI(4),ETA(4)COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)DATA XI/-1.0,1.0,1.0,-1.0/DATA ETA/-1.0,-1.0,1.0,1.0/DO 10 I=1,4G1=(1.0+XI(I)*R)G2=(1.0+ETA(I)*S)FUN(I)=0.25*G1*G2PN(1,I)=0.25*XI(I)*G2PN(2,I)=0.25*ETA(I)*G110 CONTINUEDO 80 I=1,2DO 75 J=1,2DET=0.00DO 70 K=1,4DET=DET+PN(I,K)*XYZ(J,K)70 CONTINUEXJAC(I,J)=DET75 CONTINUE80 CONTINUEDET=XJAC(1,1)*XJAC(2,2)* -XJAC(2,1)*XJAC(1,2)RETURNENDC **********************************************SUBROUTINE OUTDISP(NP,R,JR)DIMENSION R(*),JR(2,*),U(2)WRITE(*,650)WRITE(7,650)DO I=1,NPDO M=1,2L=JR(M,I)IF(L.EQ.0)U(M)=0.0IF(L.GT.0)U(M)=R(L)ENDDOWRITE(*,'(5X,I5,10X,2E14.3)') I,UWRITE(7,'(5X,I5,10X,2E14.3)') I,UENDDO650 FORMAT(/25X,'NODAL DISPLACEMENTS'/8X,* 'NODE',13X,'X-COMP.',8X,'Y-COMP.')RETURN ENDC **********************************************SUBROUTINE OUTSTRE(NE,STRE)DIMENSION STRE(3,*),ST(6)WRITE(*,700)WRITE(7,700)DO IE=1,NESX=STRE(1,IE)SY=STRE(2,IE)SXY=STRE(3,IE)ST(1)=SXST(2)=SYST(3)=SXYH1=SX+SYH2=SQRT((SX-SY)*(SX-SY)+4.0*SXY*SXY)ST(4)=(H1+H2)/2.0ST(5)=(H1-H2)/2.0IF(ABS(SXY).LT.1.0E-4)THENIF (SX.GT.SY) ST(6)=0.0IF (SX.LE.SY) ST(6)=90.0ELSEST(6)=ATAN((ST(4)-SX)/SXY)*57.29578ENDIFWRITE(*,'(6X,I4,3X,6F11.3)') IE,STWRITE(7,'(6X,I4,3X,6F11.3)') IE,STENDDO700 FORMAT(/30X,'ELEMENT STRESSES'/5X,*'ELEMENT',5X,'X-STRESS',3X,'Y-STRESS',*2X,'XY-STRESS',1X,'MAX-STRESS',1X,*'MIN-STRESS',4X,'ANGLE'/)RETURNENDC *********************************************。
基于matlab GUI的平面四杆机构的运动分析
基于matlab GUI的平面四杆机构的运动分析一、目的通过matlab对平面四杆机构进行运动仿真，并以GUI界面方式实现输入输出的参数化，对平面四杆机构进行位置分析、速度分析、加速度分析和静力学分析。
此外，通过动画演示，更加形象直观地观察机构的运动过程。
最后，将程序编译成.exe独立可执行文件，可以在其它没有安装matlab的机器上运行。
二、设计思路通过matlab的GUI功能模块，创建一个图形用户界面，在自动生成的代码框架中对初始化函数和回调函数等进行编辑，建立与控件相关联的程序：控件属性、位置分析、速度分析、加速度分析、静力学分析、动画演示等。
图1是平面四杆机构的示意图，输入角q的运动规律为q=pi/50*t^2+q0，r1、r2是从动角。
对t时刻沿着杆长距离原点A的任意一点进行分析。
注意：输入输出角的单位为度，时间t的取值范围为0：0.05：10，任意点lx的取值范围为0～a1+a2+a3，估算的从动角r1、r2的迭代初始值不能偏离平衡位置太大。
图1、平面四杆机构示意图三、设计流程1、通过GUI模块创建图形用户界面命令方式：在Matlab命令窗口键入>>guide；菜单方式：在Matlab的主窗口中，选择File>New>GUI命令，就会显示GUI的设计模板。
如图1所示。
图2、创建图形界面2、设计图形界面在创建之后的图形界面中插入坐标轴axes，静态文本框static text，编辑文本框edit text，按钮push button等等。
如图所示。
图3、图形界面设计3、编辑回调函数1）位置分析：输入角的函数为：q=pi/50*t^2+q0。
在时间t=0~10s内，每一个时间点估算两个初始从动角，根据牛顿-拉普森迭代得到准确的机构位置。
10s刚好主动角经历了360度，记录每一时刻的位置，便可以动画演示。
2）速度分析：输入角速度为：dq=pi/25*t。
选择杆件上的任意一点（坐标表示为质点沿着杆件到原点A的距离）做分析，正确表达出角速度系数和速度系数，便可以求出质点的速度。
单元刚度矩阵(等参元)MATLAB编程
《有限元法》实验报告专业班级力学（实验）1601 姓名田诗豪学号10提交日期实验一（30分）一、实验内容编写一个计算平面3结点三角形单元的应变矩阵、应力矩阵和单元刚度矩阵的MATLAB 函数文件[B3,S3,K3] = ele_mat_tri3(xy3,mat)，其中：输入变量xy3为结点坐标数组，mat为材料参数矩阵；输出变量B3为应变矩阵，S3为应力矩阵，K3为单元刚度矩阵。
（要求给出3个不同算例进行验证，并绘制出单元形状和结点号）二、程序代码通用函数function [B3,S3,K3] = ele_mat_tri3(xy3,mat)%生成平面3结点三角形单元的应变矩阵、应力矩阵和单元刚度矩阵的功能函数%***变量说明****%xy3------------------结点坐标数组%mat------------------材料参数矩阵（弹性模量，泊松比，壁厚）%B3-------------------应变矩阵%S3-------------------应力矩阵%K3-------------------单元刚度矩阵%***xyh=[1,xy3(1,1),xy3(1,2);1,xy3(2,1),xy3(2,2);1,xy3(3,1),xy3(3,2)];A=*det(xyh);A=abs(A);D=mat(1)/(1-mat(2)^2)*[1,mat(2),0;mat(2),1,0;0,0,(1-mat(2))/2];b=zeros(1,3);c=zeros(1,3);%***for i=1:3if i==1j=2;m=3;elseif i==2j=3;m=1;elsej=1;m=2;endb(i)=xy3(j,2)-xy3(m,2);c(i)=xy3(m,1)-xy3(j,1);end%***B31=1/(2*A)*[b(1),0;0,c(1);c(1),b(1)];B32=1/(2*A)*[b(2),0;0,c(2);c(2),b(2)];B33=1/(2*A)*[b(3),0;0,c(3);c(3),b(3)];B3=[B31,B32,B33];%***S3=D*B3;%***K3=A*mat(3)*B3'*D*B3;主程序clear;clc;%***输入结点坐标数组**xy3=[0,0;5,1;1,4];mat=[3e6,,];%****输入材料参数矩阵(弹性模量，泊松比，壁厚)****[B3,S3,K3]=ele_mat_tri3(xy3,mat)三、算例分析算例1：如图1所示三角形单元，结点坐标为1(0,0)，2(5,2)，3(1,4)，弹性模量为200GPa，泊松比为、厚度为。
平面四节点等参单元matlab实现
计算力学报告平面四节点等参单元学生姓名：**学号：********一、问题描述及分析在无限大平面内有一个小圆孔。
孔内有一集中力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软件计算结果计算结果分别罗列圆周方向的单元和半径方向的单元位移和应力。
matlab四节点四边形等参元的刚度矩阵计算程序
-1.2150 0.1440 2.1968 -0.8933 0.5212 -0.0151 -1.5030 0.7645
0.0616 -0.0419 -0.8933 1.8399 0.0673 -0.5250 0.7645 -1.2729
x2=2;y2=0;
x3=2.25;y3=1.5;
x4=1.25;y4=1;
材料参数：
E=30e12;
NU=0.3;
h=1;
ID=1;
Ai=1
Aj=1;
3、程序（由附件给出）
4、计算结果及讨论（需有图表来说明）
计算结果如下：
k =
1.0e+012 *
1.4619 0.2211 -1.2150 0.0616 -0.3716 -0.4236 0.1248 0.1409
-0.3716 -0.4236 0.5212 0.0673 1.1645 0.2763 -1.3141 0.0800
-0.4236 -0.2582 -0.0151 -0.5250 0.2763 0.9061 0.1624 -0.1229
0.1248 0.0585 -1.5030 0.7645 -1.3141 0.1624 2.6923 -0.9854
%输入平面问题性质参数ID（1为平面应力，2为平面应变）
%输出单元刚度矩阵
%-------------------------------------------------------------------
symsst;
a=[-(1-t)*x1+(1-t)*x2+(1+t)*x3-(1+t)*x4]/4;
平面四边形四节点等参单元Fortran源程序复习过程
平面四边形四节点等参单元F o r t r a n源程序C ************************************************C * FINITE ELEMENT PROGRAM *C * FOR Two DIMENSIONAL ELASticity PROBLEM *C * WITH 4 NODE *C ************************************************PROGRAM ELASTICITYcharacter*32 dat,cchDIMENSION SK(80000),COOR(2,300),AE(4,11),MEL(5,200),& WG(4),JR(2,300),MA(600),R(600),iew(30),STRE(3,200)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)WRITE(*,*)'PLEASE ENTER INPUT FILE NAME'READ(*,'(A)')DATOPEN(4,FILE=dat,STATUS='OLD')OPEN(7,FILE='OUT',STATUS='UNKNOWN')READ(4,*)NP,NE,NM,NRWRITE(7,'(A,I6)')'NUMBER OF NODE---------------------NP=',npWRITE(7,'(A,I6)')'NUMBER OF ELEMENT------------------NE=',ne WRITE(7,'(A,I6)')'NUMBER OF MATERIAL-----------------NM=',nm WRITE(7,'(A,I6)')'NUMBER OF surporting---------------NC=',Nr CALL INPUT (JR,COOR,AE,MEL)CALL CBAND (MA,JR,MEL)DO I=1,NHSK(I)=0.0enddoCALL SK0(SK,MEL,COOR,JR,MA,AE)do I=1,NR(I)=0.0enddopause 'aaa'stopREAD(4,*)NCP,NBE,izWRITE(*,'(5i8)')NCP,NBE,izWRITE(7,'(5i8)')NCP,NBE,izIF(NCP.GT.0)CALL CONCR(NCP,R,JR)IF(NBE.GT.0) CALL BODYR(NBE,R,MEL,COOR,JR,AE)IF(iz.GT.0)thendo jj=1,izREAD (4,*)Js,nse,(WG(I),I=1,4)read(4,*)(iew(m),m=1,nse)CALL FACER(iew,NSE,R,MEL,COOR,JR,WG)enddoendifCALL DECOP (SK,MA)CALL FOBA (SK,MA,R)CALL OUTDISP(NP,R,JR)CALL STRESS (COOR,MEL,JR,AE,R,STRE)WRITE(7,'(A)')' PROGRAM SAFF HAS BEEN ENDED'WRITE(*,'(A)')' PROGRAM SAFF HAS BEEN ENDED'STOPc RETURNENDC *********************************************SUBROUTINE INPUT (JR,COOR,AE,MEL)DIMENSION JR(2,*),COOR(2,*),AE(4,*),MEL(5,*)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHDO 70 I=1,NPREAD(4,*) IP,X,YCOOR(1,IP)=XCOOR(2,IP)=Y70 CONTINUEDO 11 J=1,NEREAD(4,*)NEE,NME,(MEL(I,NEE),I=1,4)MEL(5,NEE)=NME11 CONTINUEDO 10 I=1,NPDO 10 J=1,210 JR(J,I)=1DO 20 I=1,NRREAD(4,*) IP,IX,IYJR(1,IP)=IXJR(2,IP)=IY20 CONTINUEN=0DO 30 I=1,NPDO 30 J=1,2IF (JR(J,I)) 30,30,2525 N=N+1JR(J,I)=N30 CONTINUEDO 55 J=1,NMREAD (4,*)JJ,(AE(I,JJ),I=1,4)WRITE(*,910) JJ,(AE(I,JJ),I=1,4)55 CONTINUE910 FORMAT (/20X,'MATERIAL PROPERTIES'/(3X,I5,4(1x,E8.3))) RETURNENDC **********************************************SUBROUTINE CBAND (MA,JR,MEL)DIMENSION MA(*),JR(2,*),MEL(5,*),NN(8)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHDO 65 I=1,N65 MA(I)=0DO 90 IE=1,NEDO 75 K=1,4IEK=MEL(K,IE)DO 95 M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)95 CONTINUE75 CONTINUEL=NDO 80 I=1,2*4NNI=NN(I)IF(NNI.EQ.0) GO TO 80IF(NNI.LT.L) L=NNI80 CONTINUEDO 85 M=1,2*4JP=NN(M)IF(JP.EQ.0) GO TO 85JPL=JP-L+1IF(JPL.GT.MA(JP)) MA(JP)=JPL85 CONTINUE90 CONTINUEMX=0MA(1)=1DO 10 I=2,NIF(MA(I).GT.MX) MX=MA(I)MA(I)=MA(I)+MA(I-1)10 CONTINUENH=MA(N)WRITE(7,'(A,I8)')'TOTAL DEGREES OF FREEDOM-----------N= ',N WRITE(7,'(A,I8)')'MAX-SEMI-BANDWIDTH-----------------MX=',MX WRITE(7,'(A,I8)')'TOTAL-STORAGE----------------------NH=',NH500 FORMAT (/5X,'FREEDOM N='*,I5,3X,'SEMI-BANDWI. MX=',I5,3X,* 'STORAGE NH=',I7)RETURNENDC **********************************************SUBROUTINE SK0(SK,MEL,COOR,JR,MA,AE)DIMENSION SK(*),MEL(5,*),COOR(2,*),JR(2,*),MA(*), * AE(4,*),XYZ(2,4),iven(4)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)COMMON /CMN4/ NEE,NMECOMMON /GAUSS/ RSTG(3),H(3)H(1)=0.5555555555555560H(2)=0.8888888888888890H(3)=H(1)RSTG(1)=-0.7745966692414830RSTG(2)=0.00RSTG(3)=-RSTG(1)DO 10 IE=1,NENEE=IENME=MEL(5,IE)DO 75 K=1,4IEK=MEL(K,IE)iven(k)=IEKDO 95 M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)95 XYZ(M,K)=COOR(M,IEK)75 CONTINUECALL STIF(XYZ,AE,iven)DO 60 I=1,8DO 60 J=1,8II=NN(I)JJ=NN(J)IF ((JJ.EQ.0).OR.(II.LT.JJ)) GO TO 60JN=MA(II)-(II-JJ)SK(JN)=SK(JN)+SKE(I,J)60 CONTINUE70 CONTINUEwrite(7,1111) ((ske(i,j),j=1,8),i=1,8)1111 format(2x,8f12.2)10 CONTINUERETURNENDC *********************************************SUBROUTINE STIF(XYZ,AE,iven)DIMENSION AE(4,*),DNX(2,4),XYZ(2,*),iven(*),* RJAC(2,2)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)COMMON /CMN4/ NEE,NMECOMMON /GAUSS/ RSTG(3),H(3)DO 40 I=1,8RF(I)=0.00DO 30 J=1,8SKE(I,J)=0.0030 CONTINUE40 CONTINUEE=AE(1,NME)U=AE(2,NME)GAMA=AE(3,NME)D1=E*(1.00-U)/((1.00+U)*(1.00-2.00*U))D2=E*U/((1.00+U)*(1.00-2.00*U))D3=E*0.50/(1.00+U)DO 120 I=1,4II=2*(I-1)I1=II+1I2=II+2DO 115 J=1,4JJ=2*(J-1)J1=JJ+1J2=JJ+2DXX=0DXY=0DYX=0DYY=0DO 99 IS=1,3S=RSTG(IS)SH=H(IS)DO 98 IR=1,3R=RSTG(IR)RH=H(IR)CALL FDNX (XYZ,DNX,DET,R,S,RJAC,iven,NEE) DNIX=DNX(1,I)DNIY=DNX(2,I)DNJX=DNX(1,J)DNJY=DNX(2,J)DXX=DXX+DNIX*DNJX*DET*RH*SHDXY=DXY+DNIX*DNJY*DET*RH*SHDYX=DYX+DNIY*DNJX*DET*RH*SHDYY=DYY+DNIY*DNJY*DET*RH*SH98 CONTINUE99 CONTINUESKE(I1,J1)=DXX*D1+DYY*D3SKE(I2,J2)=DYY*D1+DXX*D3SKE(I1,J2)=DXY*D2+DYX*D3SKE(I2,J1)=DYX*D2+DXY*D3115 CONTINUE120 CONTINUERETURNENDC ********************************************* SUBROUTINE CONCR(NCP,R,JR)DIMENSION R(*),JR(2,*),XYZ(2)DO 100 I=1,NCPREAD (4,*) IP,PX,PYXYZ(1)=PXXYZ(2)=PYDO 95 J=1,2L=JR(J,IP)IF(L.EQ.0) GO TO 95R(L)=R(L)+XYZ(J)95 CONTINUE100 CONTINUERETURNENDC ********************************************** SUBROUTINE BODYR(NBE,R,MEL,COOR,JR,AE) DIMENSION R(*),MEL(5,*),COOR(2,*),JR(2,*),& AE(4,*),XYZ(2,4),iven(4)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)COMMON /GAUSS/ RSTG(3),H(3)H(1)=1.0H(2)=1.0RSTG(1)=-0.5773502691896260RSTG(2)=-RSTG(1)DO 10 IE=1,NBEDO I=1,8RF(I)=0.00ENDDOc READ(4,*)NEENEE=ieNME=MEL(5,NEE)GAMA=AE(3,NME)DO 75 K=1,4IEK=MEL(K,NEE)iven(k)=iekDO 95 M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)95 XYZ(M,K)=COOR(M,IEK)75 CONTINUEDO 99 IS=1,2S=RSTG(IS)SH=H(IS)DO 98 IR=1,2RR=RSTG(IR)RH=H(IR)CALL FUN8 (XYZ,RR,S,DET)DO 30 I=1,4J=2*IRF(J)=RF(J)-FUN(I)*RH*SH*DET*GAMA30 CONTINUE98 CONTINUE99 CONTINUECALL ASLOAD (R)10 CONTINUERETURNENDC *********************************************SUBROUTINE FACER(iew,NSE,R,MEL,COOR,JR,WG) DIMENSION R(*),MEL(5,*),COOR(2,*),JR(2,*),wg(*) * ,XYZ(2,4),iew(*),PR(2)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)COMMON /CMN4/ NEE,NMECOMMON /GAUSS/ RSTG(3),H(3)H(1)=1.0H(2)=1.0RSTG(1)=-0.5773502691896260RSTG(2)=-RSTG(1)nwf=0nnf=0ir=wg(1)+0.1if(ir.eq.1)nwf=1if(ir.eq.2)nnf=1DO 510 IE=1,NSEDO I=1,8RF(I)=0.00ENDDOnee=iew(ie)DO 575 K=1,4IEK=MEL(K,NEE)DO 595 M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)595 XYZ(M,K)=COOR(M,IEK)575 CONTINUEIF(NWF.EQ.1) thenGAMA=WG(2)Z0=WG(3)NSU=WG(4)+0.1CALL SURLOD (NSU,XYZ,PR,Z0,GAMA,1)endifIF(NNF.EQ.1) thenq=WG(2)NSU=WG(4)+0.1do j=1,2PR(J)=qenddoCALL SURLOD (NSU,XYZ,PR,Z0,GAMA,2)endifCALL ASLOAD (R)510 CONTINUERETURNENDC *********************************************SUBROUTINE SURLOD (NSU,XYZ,PR,Z0,GAMA,NSI)DIMENSION XYZ(2,*),RST(3),PR(2),KCRD(4),KFACE(2,4), & FVAL(4),NODES(2),FACT(4)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)COMMON /CMN4/ NEE,NMECOMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)COMMON /GAUSS/ RSTG(3),H(3)DATA KCRD/1,1,2,2/DATA KFACE/1, 4,* 2, 3,* 1, 2,* 4, 3/DATA FVAL/-1.00,1.00,-1.00,1.00/FACT(1)=1.0FACT(2)=-1.0FACT(3)=-1.0FACT(4)=1.0FACTNUS=FACT(NSU)DO I=1,2J=KFACE(I,NSU)NODES(I)=JENDDOIF (NSI.EQ.1) THENDO I=1,2J=NODES(I)Z=Z0-XYZ(2,J)PR(I)=0.00IF (Z.GT.0.00) PR(I)=Z*GAMAENDDOENDIFML=KCRD(NSU)IF(ML.EQ.1)MM=2IF(ML.EQ.2)MM=1RST(ML)=FVAL(NSU)DO 70 LX=1,2RST(MM)=RSTG(LX)CALL FUN8 (XYZ,RST(1),RST(2),DET) PXYZ=0.00DO 25 I=1,2J=NODES(I)PXYZ=PXYZ+FUN(J)*PR(I)25 CONTINUEA1=XJAC(MM,2)A2=-XJAC(MM,1)30 DO 60 I=1,2J=NODES(I)K2=2*JK1=K2-1Q=PXYZ*FUN(J)*H(LX)*FACTNUSRF(K1)=RF(K1)+Q*A1RF(K2)=RF(K2)+Q*A260 CONTINUE70 CONTINUEENDC ********************************************* SUBROUTINE ASLOAD (R)DIMENSION R(*)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)DO 20 I=1,8L=NN(I)IF (L.EQ.0) GO TO 20R(L)=R(L)+RF(I)20 CONTINUERETURNENDC *********************************************** SUBROUTINE DECOP (SK,MA)DIMENSION SK(*),MA(*)COMMON /CMN2/ N,MX,NHDO 50 I=2,NL=I-MA(I)+MA(I-1)+1K=I-1L1=L+1IF (L1.GT.K) GO TO 30DO 20 J=L1,KIJ=MA(I)-I+JM=J-MA(J)+MA(J-1)+1IF (L.GT.M) M=LMP=J-1IF (M.GT.MP) GO TO 20DO 10 LP=M,MPIP=MA(I)-I+LPJP=MA(J)-J+LPSK(IJ)=SK(IJ)-SK(IP)*SK(JP)10 CONTINUE20 CONTINUE30 IF (L.GT.K) GO TO 50DO 40 LP=L,KIP=MA(I)-I+LPLPP=MA(LP)SK(IP)=SK(IP)/SK(LPP)II=MA(I)SK(II)=SK(II)-SK(IP)*SK(IP)*SK(LPP)40 CONTINUE50 CONTINUEENDC *************************************************************SUBROUTINE FOBA (SK,MA,R)DIMENSION SK(*),MA(*),R(*)COMMON /CMN2/ N,MX,NHDO 10 I=2,NL=I-MA(I)+MA(I-1)+1K=I-1IF (L.GT.K) GO TO 10DO 5 LP=L,KIP=MA(I)-I+LPR(I)=R(I)-SK(IP)*R(LP)5 CONTINUE10 CONTINUEDO 20 I=1,NII=MA(I)45 R(I)=R(I)/SK(II)20 CONTINUEDO 30 J1=2,NI=2+N-J1L=I-MA(I)+MA(I-1)+1K=I-1IF (L.GT.K) GO TO 30DO 25 J=L,KIJ=MA(I)-I+J55 R(J)=R(J)-SK(IJ)*R(I)25 CONTINUE30 CONTINUERETURNENDC ***************************************************************** SUBROUTINE STRESS(COOR,MEL,JR,AE,R,STRE)DIMENSION XYZ(2,4),DNX(2,4),AE(4,*),STRE(3,*),& COOR(2,*),MEL(5,*),JR(2,*),RJAC(2,2),SIG(3),& B(3,8),R(*),iven(4)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)DO 106 IE=1,NENME=MEL(5,IE)DO 300 K=1,4IEK=MEL(K,IE)DO 310 M=1,2310 XYZ(M,K)=COOR(M,IEK)DO 320 M=1,2JRR=2*(K-1)+M320 NN(JRR)=JR(M,IEK)300 CONTINUEE=AE(1,NME)U=AE(2,NME)D1=E*(1.00-U)/((1.00+U)*(1.00-2.00*U))D2=E*U/((1.00+U)*(1.00-2.00*U))D3=0.50*E/(1.00+U)SS=0.0RR=0.0CALL FDNX (XYZ,DNX,DET,RR,SS,RJAC,iven,IE) DO 30 I=1,4II=2*(I-1)J1=II+1J2=II+2BI=DNX(1,I)CI=DNX(2,I)B(1,J1)=BIB(2,J1)=0.B(3,J1)=CIB(1,J2)=0.B(2,J2)=CIB(3,J2)=BI30 CONTINUEDO 55 II=1,3SIG(II)=0.0055 CONTINUEDO 70 K=1,8NA=NN(K)IF (NA.EQ.0) GO TO 70DO 60 L=1,3SIG(L)=SIG(L)+B(L,K)*R(NA)60 CONTINUE70 CONTINUESX=D1*SIG(1)+D2*SIG(2)SY=D2*SIG(1)+D1*SIG(2)SXY=D3*SIG(3)STRE(1,IE)=SXSTRE(2,IE)=SYSTRE(3,IE)=SXY106 CONTINUECALL OUTSTRE(NE,STRE)RETURNENDC *********************************************SUBROUTINE FDNX (XYZ,DNX,DET,R,S,RJAC,iven,NEE) DIMENSION XYZ(2,*),DNX(2,*),RJAC(2,2),iven(*)COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)CALL FUN8 (XYZ,R,S,DET)IF (DET.LT.1.0E-5)THENWRITE(7,600) NEE,R,S,detWRITE(7,*) (iven(m),m=1,4)STOPENDIFREC=1.00/DETRJAC(1,1)=REC*XJAC(2,2)RJAC(2,2)=REC*XJAC(1,1)RJAC(2,1)=-REC*XJAC(2,1)RJAC(1,2)=-REC*XJAC(1,2)DO 30 K=1,4DO 20 I=1,2DNX(I,K)=0.DO 25 M=1,2DNX(I,K)=DNX(I,K)+RJAC(I,M)*PN(M,K)25 CONTINUE20 CONTINUE30 CONTINUE600 FORMAT (1X,'ERR0R*** NEGTIVE OR ZERO '* 'JACOBIAN DETERMINANT FOR '* 'ELEMENT'/'ELE.=',I5,' R=',F10.5,6X,'S=',F10.5,* 'det=',f12.5)RETURNENDC *********************************************SUBROUTINE FUN8 (XYZ,R,S,DET)DIMENSION XYZ(2,*),XI(4),ETA(4)COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)DATA XI/-1.0,1.0,1.0,-1.0/DATA ETA/-1.0,-1.0,1.0,1.0/DO 10 I=1,4G1=(1.0+XI(I)*R)G2=(1.0+ETA(I)*S)FUN(I)=0.25*G1*G2PN(1,I)=0.25*XI(I)*G2PN(2,I)=0.25*ETA(I)*G110 CONTINUEDO 80 I=1,2DO 75 J=1,2DET=0.00DO 70 K=1,4DET=DET+PN(I,K)*XYZ(J,K)70 CONTINUEXJAC(I,J)=DET75 CONTINUE80 CONTINUEDET=XJAC(1,1)*XJAC(2,2)* -XJAC(2,1)*XJAC(1,2)RETURNENDC ********************************************** SUBROUTINE OUTDISP(NP,R,JR)DIMENSION R(*),JR(2,*),U(2)WRITE(*,650)WRITE(7,650)DO I=1,NPDO M=1,2L=JR(M,I)IF(L.EQ.0)U(M)=0.0IF(L.GT.0)U(M)=R(L)ENDDOWRITE(*,'(5X,I5,10X,2E14.3)') I,UWRITE(7,'(5X,I5,10X,2E14.3)') I,UENDDO650 FORMAT(/25X,'NODAL DISPLACEMENTS'/8X, * 'NODE',13X,'X-COMP.',8X,'Y-COMP.')RETURNENDC ********************************************** SUBROUTINE OUTSTRE(NE,STRE)DIMENSION STRE(3,*),ST(6)WRITE(*,700)WRITE(7,700)DO IE=1,NESX=STRE(1,IE)SY=STRE(2,IE)SXY=STRE(3,IE)ST(1)=SXST(2)=SYST(3)=SXYH1=SX+SYH2=SQRT((SX-SY)*(SX-SY)+4.0*SXY*SXY)ST(4)=(H1+H2)/2.0ST(5)=(H1-H2)/2.0IF(ABS(SXY).LT.1.0E-4)THENIF (SX.GT.SY) ST(6)=0.0IF (SX.LE.SY) ST(6)=90.0ELSEST(6)=ATAN((ST(4)-SX)/SXY)*57.29578ENDIFWRITE(*,'(6X,I4,3X,6F11.3)') IE,STWRITE(7,'(6X,I4,3X,6F11.3)') IE,STENDDO700 FORMAT(/30X,'ELEMENT STRESSES'/5X,*'ELEMENT',5X,'X-STRESS',3X,'Y-STRESS',*2X,'XY-STRESS',1X,'MAX-STRESS',1X,*'MIN-STRESS',4X,'ANGLE'/)RETURNENDC *********************************************。
平面四节点等参单元分析程序
\变分原理与有限元大作业平面四节点等参单元分析程序.姓名：潘清学号：SQ完成时间：2011-4-26：一、概述通常情况下的有限元分析过程是运用可视化分析软件(如ANSYS、ABAQUS、SAP等)进行前处理和后处理，而中间的计算部分一般采用自己编制的程序来运算。
具有较强数值计算和处理能力的Fortran语言是传统有限元计算的首选语言。
随着有限元技术的逐步成熟,它被应用在越来越复杂的问题处理中，但在实际应用中也暴露出一些问题。
有时网格离散化的区域较大,而又限于研究精度的要求，使得划分的网格数目极其庞大，结点数可多达数万个，从而造成计算中要运算的数据量巨大，程序运行的时间较长的弊端，这就延长了问题解决的时间，使得求解效率降低。
因为运行周期长，不利于程序的调试，特别是对于要计算多种运行工况时的情况；同时大数据量处理对计算机的内存和CPU 提出了更高的要求，而在实际应用中，单靠计算机硬件水平的提高来解决问题的能力是有限的。
因此，必须寻找新的编程语言。
随着有限元前后处理的不断发展和完善，以及大型工程分析软件对有限元接口的要求，有限元分析程序不应只满足解题功能，它还应满足软件工程所要求的结构化程序设计条件，能够对存储进行动态分配，以充分利用计算机资源，它还应很容易地与其它软件如CAD 的实体造型，优化设计等接口。
现在可编写工程应用软件的计算机语言较多，其中C语言是一个较为优秀的语言，很容易满足现在有限元分析程序编程的要求。
C语言最初是为操作系统、编译器以及文字处理等编程而发明的。
随着不断完善，它已应用到其它领域，包括工程应用软件的编程。
近年来，C语言已经成为计算机领域最普及的一个编程语言，几乎世界上所有的计算机都装有C的编译器，从PC机到巨型机到超巨型的并行机，C与所有的硬件和操作系统联系在一起。
用C 编写的程序，可移植性极好，几乎不用作多少修改，就可在任何一台装有ANSI、C编译器的计算机上运行。
C既是高级语言，也是低级语言，也就是说，可用它作数值计算，也可用它对计算机存储进行操作。
单元刚度矩阵MATLAB编程
《有限元法》实验报告专业班级力学（实验）1601 姓名田诗豪学号 10提交日期实验一（30分）一、实验内容编写一个计算平面3结点三角形单元的应变矩阵、应力矩阵和单元刚度矩阵的MATLAB 函数文件[B3,S3,K3] = ele_mat_tri3(xy3,mat)，其中：输入变量xy3为结点坐标数组，mat 为材料参数矩阵；输出变量B3为应变矩阵，S3为应力矩阵，K3为单元刚度矩阵。
（要求给出3个不同算例进行验证，并绘制出单元形状和结点号）二、程序代码通用函数function [B3,S3,K3] = ele_mat_tri3(xy3,mat)%生成平面3结点三角形单元的应变矩阵、应力矩阵和单元刚度矩阵的功能函数%*********变量说明****************%xy3------------------结点坐标数组%mat------------------材料参数矩阵（弹性模量，泊松比，壁厚）%B3-------------------应变矩阵%S3-------------------应力矩阵%K3-------------------单元刚度矩阵%*********************************xyh=[1,xy3(1,1),xy3(1,2);1,xy3(2,1),xy3(2,2);1,xy3(3,1),xy3(3,2)];A=*det(xyh);A=abs(A);D=mat(1)/(1-mat(2)^2)*[1,mat(2),0;mat(2),1,0;0,0,(1-mat(2))/2];b=zeros(1,3);c=zeros(1,3);%*********************************for i=1:3if i==1j=2;m=3;elseif i==2j=3;m=1;elsej=1;m=2;endb(i)=xy3(j,2)-xy3(m,2);c(i)=xy3(m,1)-xy3(j,1);end%*********************************B31=1/(2*A)*[b(1),0;0,c(1);c(1),b(1)];B32=1/(2*A)*[b(2),0;0,c(2);c(2),b(2)];B33=1/(2*A)*[b(3),0;0,c(3);c(3),b(3)];B3=[B31,B32,B33];%*********************************S3=D*B3;%*********************************K3=A*mat(3)*B3'*D*B3;主程序clear;clc;%*********输入结点坐标数组********xy3=[0,0;5,1;1,4];mat=[3e6,,];%****输入材料参数矩阵(弹性模量，泊松比，壁厚)****[B3,S3,K3]=ele_mat_tri3(xy3,mat)三、算例分析算例1：如图1所示三角形单元，结点坐标为1(0,0)，2(5,2)，3(1,4)，弹性模量为200GPa，泊松比为、厚度为。
基于MATLAB的四结点等参元的编程
基于MATLAB的四结点等参元的编程
栾小兵;李自林;朱斌
【期刊名称】《天津城市建设学院学报》
【年(卷),期】2006(012)004
【摘要】基于有限元理论,用MATLAB语言编制了四结点等参元计算机程序.计算了悬臂双向板在均布荷载与集中力共同作用下的位移、应力和应变,并与大型通用有限元软件ANSYS计算结果进行了比较,二者吻合较好,验证了本程序编制的正确性.利用MATLAB编程对结构进行分析计算比FORTRAN、C语言简便,具有精度高的优点.
【总页数】4页(P259-262)
【作者】栾小兵;李自林;朱斌
【作者单位】天津城市建设学院,土木工程系,天津,300384;天津城市建设学院,土木工程系,天津,300384;天津城市建设学院,土木工程系,天津,300384
【正文语种】中文
【中图分类】TP314
【相关文献】
1.四结点等参元XFEM程序设计及在裂纹问题中的应用 [J], 苏毅;王生楠;闫晓中
2.机械密封环应力场的四结点等参元分析 [J], 皇甫晶;周启慧;魏昭成
3.简明有效的四结点多变量协调等参元 [J], 史宝军;鹿晓阳
4.平面四结点多变量非协调等参元 [J], 史宝军;袁明武;等
5.基于“求逆犯规”修正方法的平面四结点非协调等参元 [J], 鹿晓阳
因版权原因，仅展示原文概要，查看原文内容请购买。
MATLAB程序：已知三个位置设计平面四杆机构求解程序(位移矩阵法)
%MATLAB程序：已知三个位置设计平面四杆机构求解程序（位移矩阵法）clear;clc;%凡是变量名前带v的为数值变量，不带的是符号变量vxp1=0; vyp1=0; vsita1=0*pi/180;vxp2=-2; vyp2=6; vsita2=40*pi/180;vxp3=-10; vyp3=8; vsita3=90*pi/180; %精确位置P1，P2，P3及各角度vsita12=vsita2-vsita1;vsita13=vsita3-vsita1;vxa=-10; vya=-2;vxd=-5; vyd=-2; %选定A,D点%所有数值均在此确定，更改此处即可解出不同数值的四杆机构位移矩阵方程syms xp1 yp1 xp2 yp2 xp3 yp3 sita12 sita13;syms xa ya xb1 yb1 xb2 yb2 xb3 yb3;f1='(xb2-xa)^2+(yb2-ya)^2=(xb1-xa)^2+(yb1-ya)^2';f2='(xb3-xa)^2+(yb3-ya)^2=(xb1-xa)^2+(yb1-ya)^2'; %前两个机构方程f3='xb2=cos(sita12)*xb1-sin(sita12)*yb1+xp2-xp1*cos(sita12)+yp1*sin(sita12)';f4='yb2=sin(sita12)*xb1+cos(sita12)*yb1+yp2-xp1*sin(sita12)-yp1*cos(sita12)'; %由第一个位移矩阵方程得出f5='xb3=cos(sita13)*xb1-sin(sita13)*yb1+xp3-xp1*cos(sita13)+yp1*sin(sita13)';f6='yb3=sin(sita13)*xb1+cos(sita13)*yb1+yp3-xp1*sin(sita13)-yp1*cos(sita13)'; %由第二个位移矩阵方程得出f1=subs(f1,{xa,ya},{vxa,vya});f2=subs(f2,{xa,ya},{vxa,vya});f3=subs(f3,{xp1,xp2,yp1,sita12},{vxp1,vxp2,vyp1,vsita12});f4=subs(f4,{xp1,yp1,yp2,sita12},{vxp1,vyp1,vyp2,vsita12});f5=subs(f5,{xp1,xp3,yp1,sita13},{vxp1,vxp3,vyp1,vsita13});f6=subs(f6,{xp1,yp1,yp3,sita13},{vxp1,vyp1,vyp3,vsita13}); %代入具体数值[xb1,xb2,xb3,yb1,yb2,yb3]=solve(f1,f2,f3,f4,f5,f6); %解方程vxb1=vpa(xb1);vyb1=vpa(yb1);vxb2=vpa(xb2);vyb2=vpa(yb2);vxb3=vpa(xb3);vyb3=vpa(yb3);(vxb1-vxa)^2+(vyb1-vya)^2;(vxb2-vxa)^2+(vyb2-vya)^2;(vxb3-vxa)^2+(vyb3-vya)^2; %去掉这三行分号可验证B点三个位置是否距离A点相等syms xd yd xc1 yc1 xc2 yc2 xc3 yc3;f7='(xc2-xd)^2+(yc2-yd)^2=(xc1-xd)^2+(yc1-yd)^2';f8='(xc3-xd)^2+(yc3-yd)^2=(xc1-xd)^2+(yc1-yd)^2'; %前两个机构方程f9='xc2=cos(sita12)*xc1-sin(sita12)*yc1+xp2-xp1*cos(sita12)+yp1*sin(sita12)';f10='yc2=sin(sita12)*xc1+cos(sita12)*yc1+yp2-xp1*sin(sita12)-yp1*cos(sita12)'; %由第一个位移矩阵方程得出f11='xc3=cos(sita13)*xc1-sin(sita13)*yc1+xp3-xp1*cos(sita13)+yp1*sin(sita13)';f12='yc3=sin(sita13)*xc1+cos(sita13)*yc1+yp3-xp1*sin(sita13)-yp1*cos(sita13)'; %由第二个位移矩阵方程得出f7=subs(f7,{xd,yd},{vxd,vyd});f8=subs(f8,{xd,yd},{vxd,vyd});f9=subs(f9,{xp1,xp2,yp1,sita12},{vxp1,vxp2,vyp1,vsita12});f10=subs(f10,{xp1,yp1,yp2,sita12},{vxp1,vyp1,vyp2,vsita12});f11=subs(f11,{xp1,xp3,yp1,sita13},{vxp1,vxp3,vyp1,vsita13});f12=subs(f12,{xp1,yp1,yp3,sita13},{vxp1,vyp1,vyp3,vsita13}); %代入具体数值[xc1,xc2,xc3,yc1,yc2,yc3]=solve(f7,f8,f9,f10,f11,f12); %解方程vxc1=vpa(xc1);vyc1=vpa(yc1);vxc2=vpa(xc2);vyc2=vpa(yc2);vxc3=vpa(xc3);vyc3=vpa(yc3);(vxc1-vxd)^2+(vyc1-vyd)^2;(vxc2-vxd)^2+(vyc2-vyd)^2;(vxc3-vxd)^2+(vyc3-vyd)^2; %去掉这三行分号可验证C点三个位置是否距离D点相等%最终答案xb1,yb1,xc1,yc1Lab=sqrt((vxb1-vxa)^2+(vyb1-vya)^2)Lbc=sqrt((vxb1-vxc1)^2+(vyb1-vyc1)^2)Lcd=sqrt((vxc1-vxd)^2+(vyc1-vyd)^2)Lad=sqrt((vxa-vxd)^2+(vya-vyd)^2) %得到四杆长'曲柄存在条件：'%得出四杆长后计算得到'可靠到位条件：'[vxc1-vxb1,vyc1-vyb1]*[vxc1-vxd,vyc1-vyd]'[vxc2-vxb2,vyc2-vyb2]*[vxc2-vxd,vyc2-vyd]'[vxc3-vxb3,vyc3-vyb3]*[vxc3-vxd,vyc3-vyd]''顺序到位条件：'%未完成输出结果：xb1 =(-7-4*sin(2/9*pi)+4*cos(2/9*pi))/(4*cos(2/9*pi)+4*sin(2/9*pi)-5)yb1 =-1xc1 =-6*(27+24*sin(2/9*pi)-64*cos(2/9*pi))/(-31*cos(2/9*pi)-5+39*sin(2/9*pi))yc1 =-2*(72*cos(2/9*pi)-175+192*sin(2/9*pi))/(-31*cos(2/9*pi)-5+39*sin(2/9*pi))Lab =1.0288436025165976748172169832223Lbc =2.9872531417317691216303250912289Lcd =6.9831476545729886023199865357226Lad =5ans =曲柄存在条件：ans =可靠到位条件：ans =14.605219997928496422368168445525ans =19.799913716084881287517588922012ans =20.814756669957613005391246805307ans =顺序到位条件：。
- 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处理四节点单元节点应力与程序处理方式差异有关。