平面四节点等参单元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软件计算结果计算结果分别罗列圆周方向的单元和半径方向的单元位移和应力。
选取半径方向的单元编号为：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处理四节点单元节点应力与程序处理方式差异有关。
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 *********************************************。