三角形常应变单元程序的编制与使用
常应变三角形单元
宜
不宜
5.相邻单元的尺寸尽可能接近。 6.结点所连接的单元个数尽可能一致。 宜 二.结点编码 尽可能使相关结点的结点编码差值最小. 1 2 3 4 5 6 7 1 3 5 7 11 13 不宜
8
9
10
11 12
13
14
2
4
6
8
10
12
14
总刚半带宽=(相关结点最大差值+1)*结点位移数 总刚半带宽=(7+1)*2=16 总刚半带宽=(2+1)*2=6 总刚需占用的存贮空间为: 6 *14*2=168
宜
不宜
5.相邻单元的尺寸尽可能接近。 6.结点所连接的单元个数尽可能一致。 宜 不宜
§3.3 有限元分析应注意的问题和结果整理
一.结点的选择和单元划分 1.集中力作用点、分布力突变点、支承点应选作结点。 2.不同厚度、不同材料的部分不应划在同一个单元。 3.应力变化大处单元应密集一些。结点的多少与疏密要考虑计算 机的容量和计算精度。 4.单元边界的边长之比应尽可能靠近1。
§3 平面问题的有限元分析
§3.1 常应变三角形单元 §3.2 矩形双线性单元 §3.3 有限元分析应注意的问题和结果整理
一.结点的选择和单元划分 1.集中力作用点、分布力突变点、支承点应选作结点。 2.不同厚度、不同材料的部分不应划在同一个单元。 3.应力变化大处单元应密集一些。结点的多少与疏密要考虑计算 机的容量和计算精度。 4.单元边界的边长之比应尽可能靠近1。
总刚需占用的存贮空间为: 16 *14*2=448
三.充分利用结构的对称性 P P P P P
四.应力结果的整理
位移的计算结果一般比应力、内力结果精度高。位移达到满意结果， 由几何方程求应变，再由物理方程求应力，结果的精度较差。上述三角形 单元为常应力，矩形单元应力线性变化，而工程问题的应力是比较复杂的。 为更好地反应实际应力情况，需要对计算结果进行整理。常用处理方法有 两种:绕结点平均法和两单元平均法。
2021新有限元用三角形单元分析专业资料
平面三角形单由元整图体5分-2析2，四个单元的局部码与总码的对应关系为 平面三角形单单元整元体①分：析1，2,，4——1，2，4； 第4 二平步 面，三把角各形单单单元元元的分②贡析：献1矩，阵2叠,，加，4—即—得出2，整体4，矩5阵； 平面三角形单单元整元体③分：析1，2,，4——5，4，2； 平面三角形单单元整元体④分：析1，2,，4——4，5，6；
0 2
3 1
1 2 1 3 0 1
0 0 2 0 2 0
1 0 1 1 0 1
（t——板厚）
(4-26)
4.4 平面三角形单元分析
弹性力学平面问题三角形单元的刚度矩阵 k e 具有对称性和奇异性等特点。 三角形单元刚度矩阵 k e 可写成分块形式：
FF12e
k11 k21
有限元用三角形单元分析
4.4 平面三角形单元分析
4.4.1 平面三角形单元分析步骤
对弹性力学平面问题三角形单元进行单元分析，建立单元刚度矩阵。
图 5-16 为一个三角形单元。单元三个结点按逆时针方向排序编码 1、2、4。结点坐标(x1，y1)、(x2， y2)和(x4，y4)为已知。
在弹性力学平面问题中，每个结点有 2 个位移分量，因此，三角形单元共有 6 个自由度(u1，v1)、(u2， v2)和(u4，v4)。设结点 1 的位移向量为：
4 平面三角形单元分析
4 平面三角形各单元小分单析元 e 刚度矩阵 k e 是 6×6 阶矩阵，分块形式为式(4-29)，其中 9 个子块按局部码排列。
整体刚度矩阵 K 是 12×12 阶矩阵，分块形式为式(4-40)，其中 36 个子块按总码排列的．
平面三角形单元整体分析
三结点三角形有限单元程序设计C 语言
fclose(fp) ;
int NN2=NN*2; int NBW=NN2;
//半带宽 NBW
float GK[50][50];
float DU[50];
int *count=new int [NN];
float (*strain)[3]=new float [NN][3];
float (*stress)[3]=new float [NN][3];
int ngn[3]; for(i=0;i<3;i++)ngn[i]=NGN[NOE-1][i];
float cn[3][2]; for(i=0;i<3;i++) for(j=0;j<2;j++) cn[i][j]=CN[ngn[i]-1][j];
float a[3]={0,0,0}; float b[3]={0,0,0}; float c[3]={0,0,0};
ｖｏｉｄ ｅｌｅｓｔｒａｉｎ＿ｓｔｒｅｓｓ（
ｉｎｔＮＯＥ，ｉｎｔＮＧＮ［］［３］，ｆｌｏａｔＣＮ［］［２］，
ｆｌｏａｔＤＵ［５０］，ｆｌｏａｔＥ０，ｆｌｏａｔＵ０，ｆｌｏａｔｔ，
ｉｎｔＮＥ，ｉｎｔＮＮ，ｆｌｏａｔｅｌｅｓｔｒａｉｎ［３］，ｆｌｏａｔｅｌｅｓｔｒｅｓｓ［３］）；
＃ｅｎｄｉｆ
Ⅴ．节点应力－应变头文件 ＮｏｄｅＳｔｒａｉｎＳｔｒｅｓｓ．ｈ
fscanf(fp,"%s",&newlin); for(i=0;i<NN;i++)
for(j=0;j<2;j++) fscanf(fp, "%f", &CN[s",&newlin); int NRE=0; fscanf(fp, "%d", &NRE) ; int *restrict=new int [NRE]; for(i=0;i<NRE;i++)
常应变三角形单元
四.应力结果的整理
位移的计算结果一般比应力、内力结果精度高。位移达到满意结果， 由几何方程求应变，再由物理方程求应力，结果的精度较差。上述三角形 单元为常应力，矩形单元应力线性变化，而工程问题的应力是比较复杂的。 为更好地反应实际应力情况，需要对计算结果进行整理。常用处理方法有 两种:绕结点平均法和两单元平均法。
四.应力结果的整理
位移的计算结果一般比应力、内力结果精度高。位移达到满意结果， 由几何方程求应变，再由物理方程求应力，结果的精度较差。上述三角形 单元为常应力，矩形单元应力线性变化，而工程问题的应力是比较复杂的。 为更好地反应实际应力情况，需要对计算结果进行整理。常用处理方法有 两种:绕结点平均法和两单元平均法。
宜
不宜
5.相邻单元的尺寸尽可能接近。 6.结点所连接的单元个数尽可能一致。 宜 二.结点编码 尽可能使相关结点的结点编码差值最小. 1 2 3 4 5 6 7 1 3 5 7 9 11 13 不宜
8
9
10
11 12
13
14
2
4
6
8
10
12
14
总刚半带宽=(相关结点最大差值+1)*结点位移数 总刚半带宽=(7+1)*2=16 总刚半带宽=(2+1)*2=6 总刚需占用的存贮空间为: 6 *14*2=168
宜
不宜
5.相邻单元的尺寸尽可能接近。 6.结点所连接的单元个数尽可能一致。 宜 不宜
§3.3 有限元分析应注意的问题和结果整理
一.结点的选择和单元划分 1.集中力作用点、分布力突变点、支承点应选作结点。 2.不同厚度、不同材料的部分不应划在同一个单元。 3.应力变化大处单元应密集一些。结点的多少与疏密要考虑计算 机的容量和计算精度。 4.单元边界的边长之比应尽可能靠近1。
第6章——常应变三角形单元
[S1
S 2 S3 ]
ci 1-µ bi 2
bi E µb = Si DB = i i 2 A(1 − µ 2 ) 1-µ ci 2
µ ci
平面应变：用平面应变弹性矩阵代入得到类似结果。
22:42
有限单元法
崔向阳
16
单元应变和应力矩阵
由于同一单元中的D、B矩阵都是常数矩阵，所以S矩阵也是常 数矩阵。也就是说，三角形三节点单元内的应力分量也是常 量。 当然，相邻单元的E， µ， A和bi、ci(i，j，m)一般不完全相同， 因而具有不同的应力，这就造成在相邻单元的公共边上存在 着应力突变现象。但是随着网格的细分，这种突变将会迅速 减小。
m j
h
1 i U = ∫∫ (σ xε x + σ yε y + τ xyγ xy )hdxdy x 2 A 1 T T T T = ∫∫ σ T εhdxdy σ = ( D ε ) = ε D 2 A
−1
u1 u2 u 3
u ( x, y ) = {1 x
1 x1 y} 1 x2 1 x3
y1 y2 y3
−1
u1 u2 u 3
4
22:42
有限单元法
崔向阳
平面三角形单元
假设
{ N1
N2
N 3 } = {1 x
m 相邻单元的位移在公共边上是连续的 形函数在单元上的面积分和边界上的线积分公式为 i
j p
A ∫ ∫A Ni dxdy = 3
式中 lij 为 ij 边的长度。
Ni =1 i m j
22:42
1 ∫ij Ni dl = 2 lij
三角形单元有限元程序设计
三角形单元有限元程序设计一、引言
⑴文档背景
⑵文档目的
⑶读者对象
⑷名词解释
二、程序结构设计
⑴程序流程图
⑵程序组成模块描述
⑶主要数据结构设计
⑷代码逻辑设计
三、数据预处理
⑴数据输入格式与读取
⑵数据清洗与去噪
⑶数据预处理方法
⑷数据分割与划分
四、网格
⑴网格算法介绍
⑵网格质量评估与改善
⑶网格的实现方法
五、有限元方法
⑴有限元离散的原理
⑵有限元单元的选择
⑶有限元离散的网格节点选取
⑷有限元插值函数的计算六、模型求解
⑴线性方程组的求解方法
⑵模型参数的设置与调整
⑶迭代算法的选择与实现七、模型评估与验证
⑴结果评估指标的选择
⑵模型验证方法
⑶结果可视化与分析
八、性能优化
⑴程序性能分析与评估
⑵程序高效化的方法与技巧
⑶并行计算与优化
九、实例与案例
⑴实例介绍与问题描述
⑵实例数据处理过程
⑶实例模型求解与结果分析十、总结与展望
⑴工作总结
⑵程序改进与优化展望
⑶研究方向与发展趋势
十一、附件
●附件1、程序源代码
●附件2、输入数据样例文件
●附件3、网格结果数据文件附录：法律名词及注释
●法律名词1、注释1
●法律名词2、注释2
●法律名词3、注释3
请注意，附件的具体文件名和数据内容将根据你的实际情况进行调整。
法律名词及注释部分需要根据实际需要进行扩充和修改。
请合理对文档结构和章节内容进行调整，以符合你的具体要求。
三角形常应变单元
第五章P77，求解三角形常应变单元刚度矩阵****************************************************************************** function y=linear_triangle_element_stiffness(E,NU,t,xi,yi,xj,yj,xm,ym,p)%A=(xi*(yj-ym)+xj*(ym-yi)+xm*(yi-yj))/2;betai=yj-ym;betaj=ym-yi;betam=yi-yj;gammai=xm-xj;gammaj=xi-xm;gammam=xj-xi;B=[betai 0 betaj 0 betam 0;0 gammai 0 gammaj 0 gammam;gammai betai gammaj betaj gammam betam]/(2*A);if p==1D=(E/(1-NU*NU))*[1 NU 0;NU 1 0;0 0 (1-NU)/2];else if p==2D=(E/(1+NU)/(1-2*NU))*[1-NU NU 0;NU 1-NU 0;0 0 (1-2*NU)/2];endendBDDB=D*By=t*A*B'*D*B;******************************************************************************* function y=linear_triangle_assemble(K,k,i,j,m)%K(2*i-1,2*i-1)=K(2*i-1,2*i-1)+k(1,1);K(2*i-1,2*i)=K(2*i-1,2*i)+k(1,2);K(2*i-1,2*j-1)=K(2*i-1,2*j-1)+k(1,3);K(2*i-1,2*j)=K(2*i-1,2*j)+k(1,4);K(2*i-1,2*m-1)=K(2*i-1,2*m-1)+k(1,5);K(2*i-1,2*m)=K(2*i-1,2*m)+k(1,6);K(2*i,2*i-1)=K(2*i,2*i-1)+k(2,1);K(2*i,2*i)=K(2*i,2*i)+k(2,2);K(2*i,2*j-1)=K(2*i,2*j-1)+k(2,3);K(2*i,2*j)=K(2*i,2*j)+k(2,4);K(2*i,2*m-1)=K(2*i,2*m-1)+k(2,5);K(2*i,2*m)=K(2*i,2*m)+k(2,6);K(2*j-1,2*i-1)=K(2*j-1,2*i-1)+k(3,1);K(2*j-1,2*i)=K(2*j-1,2*i)+k(3,2);K(2*j-1,2*j-1)=K(2*j-1,2*j-1)+k(3,3);K(2*j-1,2*j)=K(2*j-1,2*j)+k(3,4);K(2*j-1,2*m-1)=K(2*j-1,2*m-1)+k(3,5);K(2*j-1,2*m)=K(2*j-1,2*m)+k(3,6);K(2*j,2*i-1)=K(2*j,2*i-1)+k(4,1);K(2*j,2*i)=K(2*j,2*i)+k(4,2);K(2*j,2*j-1)=K(2*j,2*j-1)+k(4,3);K(2*j,2*j)=K(2*j,2*j)+k(4,4);K(2*j,2*m-1)=K(2*j,2*m-1)+k(4,5);K(2*j,2*m)=K(2*j,2*m)+k(4,6);K(2*m-1,2*i-1)=K(2*m-1,2*i-1)+k(5,1);K(2*m-1,2*i)=K(2*m-1,2*i)+k(5,2);K(2*m-1,2*j-1)=K(2*m-1,2*j-1)+k(5,3);K(2*m-1,2*j)=K(2*m-1,2*j)+k(5,4);K(2*m-1,2*m-1)=K(2*m-1,2*m-1)+k(5,5);K(2*m-1,2*m)=K(2*m-1,2*m)+k(5,6);K(2*m,2*i-1)=K(2*m,2*i-1)+k(6,1);K(2*m,2*i)=K(2*m,2*i)+k(6,2);K(2*m,2*j-1)=K(2*m,2*j-1)+k(6,3);K(2*m,2*j)=K(2*m,2*j)+k(6,4);K(2*m,2*m-1)=K(2*m,2*m-1)+k(6,5);K(2*m,2*m)=K(2*m,2*m)+k(6,6);y=K;******************************************************************************* function y=linear_tiangle_element_stress(E,NU,xi,yi,xj,yj,xm,ym,p,u)%A=(xi*(yj-ym)+xj*(ym-yi)+xm*(yi-yj))/2;betai=yj-ym;betaj=ym-yi;betam=yi-yj;gammai=xm-xj;gammaj=xi-xm;gammam=xj-xi;B=[betai 0 betaj 0 betam 0;0 gammai 0 gammaj 0 gammam;gammai betai gammaj betaj gammam betam]/(2*A);if p==1D=(E/(1-NU*NU))*[1 NU 0;NU 1 0;0 0 (1-NU)/2];else if p==2D=(E/(1+NU)/(1-2*NU))*[1-NU NU 0;NU 1-NU 0;0 0 (1-2*NU)/2];endendy=D*B*u;******************************************************************************* ******************************************************************************* *************************************************************************clc;E= ;v= ;t= ;k1=linear_triangle_element_stiffness(E,v,t, , , , , , ,1) %(E,v,t,x1,y1,x2,y2, x3,y3,1)k2=linear_triangle_element_stiffness(E,v,t, , , , , , ,1) %(E,v,t,x1,y1,x2,y2, x3,y3,1)k3=linear_triangle_element_stiffness(E,v,t, , , , , , ,1) %(E,v,t,x1,y1,x2,y2, x3,y3,1)k4=linear_triangle_element_stiffness(E,v,t, , , , , , ,1) %(E,v,t,x1,y1,x2,y2, x3,y3,1)K=zeros();%节点坐标个数，节点个数乘以2K=linear_triangle_assemble(K,k1, , , );%节点编号K=linear_triangle_assemble(K,k2, , , );%节点编号K=linear_triangle_assemble(K,k3, , , );%节点编号K=linear_triangle_assemble(K,k4, , , );%节点编号KB=K([ , , , , , , , ,],:) ; %取特定的行5 6 7 8 9 10 11 12k=B(:,[ , , , , , , , ]) %取特定的列5 6 7 8 9 10 11 12f=[ ]’u=k\f %求解位移列向量U=[u(1);0;u(2);u(3);0;0;0;0]%总的位移列向量u1=[U( );U( );U( );U( );U( );U( )];%三角形单元一的位移列向量u2=[U( );U( );U( );U( );U( );U( )]; %三角形单元二的位移列向量u3=[U( );U( );U( );U( );U( );U( )];%三角形单元三的位移列向量u4=[U( );U( );U( );U( );U( );U( )]; %三角形单元四的位移列向量sigma1=linear_tiangle_element_stress(E,v, , , , , , ，1,u1)%单元一的应力大小，(E,v,x1,y1,x2,y2,x3,y3,1,u)sigma1=linear_tiangle_element_stress(E,v, , , , , , ,1,u2) %单元二的应力大小(E,v,x1,y1,x2,y2,x3,y3,1,u)sigma1=linear_tiangle_element_stress(E,v, , , , , , ,1,u3)%单元三的应力大小，(E,v,x1,y1,x2,y2,x3,y3,1,u)sigma1=linear_tiangle_element_stress(E,v, , , , , , ,1,u4) %单元四的应力大小(E,v,x1,y1,x2,y2,x3,y3,1,u)。
abaqus三角形单元创建
abaqus三角形单元创建
带有3个节点的三角形是最基础的二维有限元单元。
相比起其他单元，3节点三角形最大的特点是应变和应力在单元上的分布是相同的，所以也被称为Constant Strain Triangular (CST) 单元。
本文介绍如何通过用户自定义单元子程序UEL在ABAQUS环境中开发CST单元。
单元为平面应变单元，对标ABAQUS中的CPE3。
由于ABAQUS不支持显示用户自己开发的单元，所以需要通过UMAT辅助显示UEL计算得到的结果。
大致的流程为：
（1）在input中定义两种单元，一种是用户自己开发的“真”单元（UEL），另一种是“假”单元（abaqus内置的单元，但是材料用UMAT），用于后处理显示。
这两种单元享用相同的节点标号，但是有不同的单元标号；
（2）在UEL和UMAT中都定义COMMON block，比如: common/KUSER/ UVAR；
（3）UEL计算得到结果，将需要显示的结果赋值给SVARS （UEL中的场变量）；
（4）在UEL中把SVARS 的数值赋予UVAR；
（5）UMAT中的SDV接收UVAR，另外需要写一个弹性的本构，弹性模量 E 可以设置得很小.。
计算力学三角形单元划线法
计算力学三角形单元划线法力学三角形单元划线法（Method of Sections in Mechanics）是一种用于计算结构静力学的常用方法之一、这种方法通过在结构内部划线，将结构分为两个或更多部分，然后根据平衡条件和受力分析来求解结构的未知力和反力。
力学三角形单元划线法的基本原理是根据力的平衡条件，将结构以任意截面划分为两部分，并在划分的截面处应用力的平衡条件来求解未知力和反力。
这种方法适用于简支梁、连续梁、悬臂梁等结构的受力分析。
具体的计算步骤如下：1.选择合适的切截面：根据问题的要求和结构的特点，选择一个或多个合适的截面作为切割面。
2.划分结构：在选择的切割面上绘制一个或多个划线，将结构分为两个或更多部分。
在划割面上的划线应该是从结构的受力部分出发，沿着力的路径进行划线。
3.受力分析：对每个划分的部分应用力的平衡条件，即根据力的平衡方程（ΣF=0）来求解未知力和反力。
这些未知力和反力的数量与划割面上的未知力和反力相等。
4.解方程：根据所得到的受力分析结果，建立相应的方程组。
利用方程求解方法，如代入或联立方程等方法，求解出未知力和反力。
5.验证结果：将求解得到的未知力和反力代入原结构的其他部分，检验是否符合平衡条件。
如果结构的未知力和反力满足平衡条件，说明计算结果正确。
需要注意的是，在力学三角形单元划线法中，选择合适的切割面是十分关键的。
选择不当的切割面可能导致解方程困难或解得的方程不足以求解所需的未知力和反力。
因此，在选择切割面时，应根据结构的几何形状、受力条件和计算的目的来确定。
力学三角形单元划线法是一种简便有效的结构受力分析方法，可以应用于各种结构的计算。
通过使用这种方法，可以快速准确地求解结构的未知力和反力，对于解决工程设计和分析中的实际问题具有重要的意义。
三角形常应变单元程序的编制与使用共18页文档
三角形常应变单元程序的编制与使用有限元法是求解微分方程边值问题的一种通用数值方法，该方法是一种基于变分法（或变分里兹法）而发展起来的求解微分方程的数值计算方法，以计算机为手段，采用分片近似，进而逼近整体的研究思想求解物理问题。
有限元分析的基本步骤可归纳为三大步：结构离散、单元分析和整体分析。
对于平面问题，结构离散常用的网格形状有三角形、矩形、任意四边形，以三个顶点为节点的三角形单元是最简单的平面单元，它较矩形或四边形对曲边边界有更好Array的适应性，而矩形或四边形单元较三节点三角形有更高的计算精度。
Matlab语言是进行矩阵运算的强大工具，因此，用Matlab语言编写有限元中平面问题的程序有优越性。
本章将详细介绍如何利用Matlab语言编制三角形常应变单元的计算程序，程序流程图见图1。
有限元法中三节点三角形分析结构的步骤如下：1）整理原始数据，如材料性质、荷载条件、约束条件等，离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。
2）单元分析，建立单元刚度矩阵。
3）整体分析，建立总刚矩阵。
4）建立整体结构的等效节点荷载和总荷载矩阵5）边界条件处理。
6）解方程，求出节点位移。
7）求出各单元的单元应力。
8）计算结果整理。
计算结果整理包括位移和应力两个方面；位移计算结果一般不需要特别的处理，利用计算出的节点位移分量，就可画出结构任意方向的位移云图；而应力解的误差表现在单元内部不满足平衡方程，图1 程序流程图单元与单元边界处应力一般不连续，在边界上应力解一般与力的边界条件不相符合。
1.1 程序说明% 三角形常应变单元求解结构主程序●功能：运用有限元法中三角形常应变单元解平面问题的计算主程序。
●基本思想：单元结点按右手法则顺序编号。
●荷载类型：可计算结点荷载。
●说明：主程序的作用是通过赋值语句、读取和写入文件、函数调用等完成算法的全过程，即实现程序流程图的程序表达。
1 程序准备format short e %设定输出类型clear all %清除所有已定义变量clc %清屏●说明：format short e －设定计算过程中显示在屏幕上的数字类型为短格式、科学计数法；clear all －清除所有已定义变量，目的是在本程序的运行过程中，不会发生变量名相同等可能使计算出错的情况；clc －清屏，使屏幕在本程序运行开始时2 全局变量定义global NNODE NPION NELEM NVFIX NFORCE COORD LNODS YOUNG POISS THICKglobal FORCE FIXEDglobal BMATX DMATX SMATX AREAglobal ASTIF ASLOD ASDISPglobal FP1说明：NNODE—单元结点数，NPION—总结点数， NELEM—单元数，NVFIX—受约束边界点数，NFORCE—结点力数，COORD—结构结点坐标数组，LNODS —单元定义数组，YOUNG—弹性模量，POISS—泊松比，THICK—厚度FORCE —节点力数组(n,3) n:受力节点数目,(n,1):作用点,(n,2):x方向,(n,3):y方向； FIXED—约束信息数组(n,3) n:受约束节点数目, (n,1):约束点 (n,2)与(n,3)分别为约束点x方向和y方向的约束情况,受约束为1否则为0BMATX—单元应变矩阵(3*6)， DMATX—单元弹性矩阵(3*3)，SMATX—单元应力矩阵(3*6)，AREA—单元面积ASTIF—总体刚度矩阵，ASLOD—总体荷载向量，ASDISP—结点位移向量FP1—数据文件指针3 打开文件FP1=fopen('input.txt','rt'); %打开输入数据文件存放初始数据●说明：FP1=fopen('input.txt','rt'); －打开已存在的输入数据文件input.txt，且设置其为只读格式，使程序在执行过程中不能改变输入文件中的数值，并用文件句柄FP1来执行FP2=fopen('output.txt','wt'); －打开输出数据文件，该文件不存在时，通过此命令创建新文件，该文件存在时则将原有内容全部删除。
平面三角形单元有限元程序设计
一、题目如图1所示，一个厚度均匀的三角形薄板，在顶点作用沿板厚方向均匀分布的竖向载荷。
已知：P=150N/m，E=200GPa，=0.25，t=0.1m，忽略自重。
试计算薄板的位移及应力分布。
要求：1.编写有限元计算机程序，计算节点位移及单元应力。
（划分三角形单元，单元数不得少于30个）；2.采用有限元软件分析该问题（有限元软件网格与程序设计网格必须一致），详细给出有限元软件每一步的操作过程，并将结果与程序计算结果进行对比（任选取三个点，对比位移值）；3.提交程序编写过程的详细报告及计算机程序；4.所有同学参加答辩，并演示有限元计算程序。
有限元法中三节点三角形分析结构的步骤如下：1）整理原始数据，如材料性质、荷载条件、约束条件等，离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。
2）单元分析，建立单元刚度矩阵。
3）整体分析，建立总刚矩阵。
4）建立整体结构的等效节点荷载和总荷载矩阵5）边界条件处理。
6）解方程，求出节点位移。
7）求出各单元的单元应力。
8）计算结果整理。
一、程序设计网格划分如图，将薄板如图划分为6行，并建立坐标系，则X Y P X Y P刚度矩阵的集成建立与总刚度矩阵等维数的空矩阵，已变单元刚度矩阵的集成。
由单元分析已知节点、单元的排布规律，继而通过循环计算求得每个单元对应的节点序号。
通过循环逐个计算：（1）每个单元对应2种单元刚度矩阵中的哪一种；（2）该单元对应总刚度矩阵的那几行哪几列（3）将该单元的单元刚度矩阵加入总刚度矩阵的对应行列循环又分为3层循环：（1）最外层：逐行计算（2）中间层：该行逐个计算（3）最里层：区分为第奇/偶数个计算单元刚度的集成：[][][][][][]''''''215656665656266256561661eZeeeZeZeeeekkkKkkkkkk+⋯++=⇓=⇒==⇒==⇒=⨯⨯⨯⨯⨯⨯边界约束的处理：划0置1法适用：这种方法适用于边界节点位移分量为已知(含为0)的各种约束。
第一节三角形常应变单元
第三章平面问题的有限元法本章通过三角形常应变单元，介绍有限元法应用于弹性体应力分析的基本原理和方法：包括弹性体的离散化，单元特性的分析，刚度矩阵的建立，等效节点力的计算，解答的收敛性以及实施步骤和注意事项，同时对形函数的性质也作简要的叙述。
第一节三角形常应变单元一、结构离散化用有限元法分析弹性力学平面问题，第一步就是把原来的连续的弹性体离散化。
(a) (b)图3.1 弹性体和有限元模型将整个结构（平板）划分成有限个三角形。
这样的三角形称为单元（三角形单元）。
三角形单元的顶点取为节点。
3节点三角形单元用边界节点间的直线段来近似板的曲线边界。
这些三角形在其节点处相互连接，组成一个单元集合体，以代替原来的弹性体。
注：1. 全部节点和全部单元一般由1开始按自然顺序编号。
2. 节点编码：总码-----------用于整体分析，如1,2，…，按自然顺序编号局部码--------用于单元分析，i，j，m 要求按逆时针方向的顺序进行编码每个单元的节点局部码i，j，m和节点总码有一一对应的关系3. 单元间不能有重叠4. 一个单元的任一顶点不许为另一单元任一边的内点5. 所有作用在单元上的载荷，包括集中载荷、表面载荷和体积力，都按虚功等效的原则移置到节点上，成为等效节点载荷。
二、 位移模式1. 单元节点位移列阵iu图 3.2 平面三角形单元设单元e 的节点号码为i ，j ，m 。
由弹性力学平面可知，单元内任意一点有两个位移分量u ，v ，记为{}Tf u v ⎡⎤⎣⎦=故每个节点也有两个位移分量，因此称节点自由度为2。
3个节点得位移分量分别是 ,,,,,m m i i j j u v u v u v ，用列阵表示为{}i ei i e j j j m m m u v u v u v δδδδ⎧⎫⎪⎪⎪⎪⎧⎫⎪⎪⎪⎪⎪⎪⎪⎪⎨⎬⎨⎬⎪⎪⎪⎪⎪⎪⎪⎪⎩⎭⎪⎪⎪⎪⎩⎭== （3-1）称单元自由度是6。
其中任一子矩阵为{}Ti ii u v δ⎡⎤⎢⎥⎣⎦= （i ，j ，m 轮换）2. 位移模式结构承受载荷以后发生位移，但位移分布事先并不知道。
用常应变三角形单元解弹性力学平面问题的程序
用常应变三角形单元解弹性力学平面问题的程序******************************************************************** ANALYSIS PROGTAM OF FINITE ELEMENT METHOD ** FOR PLANE STRESS/STRAIN OF TRIANGULAR ELEMENT ** ----- FEMT3.FOR ----- **------------------------------------------------------------- ** Subroutines: 1-SDATA, 2-STE, 3-ATE, 4-DTE, 5-BTE, 6-STIFF ** 7-EQUPE, 8-INSCD, 9-BGSMT, 10-SIGME ********************************************************************DIMENSION LND(50,3),X(100),Y(100),JR(20,3),PJ(20,3),P(200)REAL KS(200,100)OPEN(5,FILE='FEMT3.DAT')OPEN(6,FILE='FEMT3.OUT',STATUS='NEW')READ(5,*) NJ,NE,NS,NPJ,IPS(结点、单元、支承、荷载、类型)WRITE(6,*)' FINITE ELEMENT ANALYSIS IN PLANE PROBLEM'WRITE(6,*)' SOURCE DATA OUTPUT'WRITE(6,20) NJ,NE,NS,NPJ,IPS20 FORMAT(4X,'NJ',3X,'NE',3X,'NS',3X,'NPJ',2X,'IPS'/1X,5I5)IF(IPS.EQ.0) WRITE(6,*)' PLANE STRESS PROBLEM'IF(IPS.EQ.1) WRITE(6,*)' PLANE STRAIN PROBLEM'CALL SDATA(NJ,NE,NS,NW,NPJ,IPS,E,PR,T,V,LND,X,Y,JR,PJ)NJ2=2*NJWRITE(6,50) NJ250 FORMAT(/1X,'DEGREES OF FREEDOM=',I5)WRITE(6,60) NW60 FORMAT(1X,'BAND WIDTH=',I5)CALL STIFF(NJ,NE,NJ2,NW,LND,X,Y,E,PR,T,KS)(总刚6)CALL EQUPE(NJ,NE,NPJ,NJ2,T,V,LND,X,Y,PJ,P)({P}7)CALL INSCD(NS,NW,NJ2,JR,KS,P)(引入支承条件8)CALL BGSMT(NJ,NJ2,NW,KS,P)（解方程9）CALL SIGME(NE,NJ,NJ2,E,PR,LND,X,Y,P)（求应力10）CLOSE(5)CLOSE(6)END*--------------------------------------------------------C SUBPROGRAM-1C INPUT STRUCTURAL DATASUBROUTINE SDATA(NJ,NE,NS,NW,NPJ,IPS,E,PR,* T,V,LND,X,Y,JR,PJ)DIMENSION LND(NE,3),X(NJ),Y(NJ),JR(NS,3),PJ(NPJ,3)READ(5,*) E,PR,T,V(弹性模量、泊松比、厚度、容重)WRITE(6,10) E,PR,T,V10 FORMAT(/6X,'E',10X,'PR',9X,'T',9X,'V'/,4F10.2)READ(5,*)((LND(I,J),J=1,3),I=1,NE)（结点编码）WRITE(6,20)20 FORMAT(/1X,'ELEMENT INFORMATION'/3X,'ELEM',3X,* 'I J K'/)WRITE(6,30)(I,(LND(I,J),J=1,3),I=1,NE)30 FORMAT(1X,4I5)READ(5,*)(X(I),Y(I),I=1,NJ)（结点坐标）WRITE(6,40)40 FORMAT(/1X,'COORDINATES OF NODES'/3X,'NODES',* 8X,'X',13X,'Y')WRITE(6,50)(I,X(I),Y(I),I=1,NJ)50 FORMAT(1X,I5,2E15.6)READ(5,*)((JR(I,J),J=1,3),I=1,NS)（约束信息）WRITE(6,60)60 FORMAT(/1X,'CONSTRAINED NODES'/3X,'NODE',3X,'X',4X,'Y') WRITE(6,70)((JR(I,J),J=1,3),I=1,NS)70 FORMAT(1X,3I5)READ(5,*)((PJ(I,J),J=1,3),I=1,NPJ)（荷载信息）WRITE(6,80)80 FORMAT(/1X,'LOAD CASES'/3X,'NODE',8X,'X',13X,'Y')WRITE(6,90)((PJ(I,J),J=1,3),I=1,NPJ)90 FORMAT(1X,F5.0,2E15.6)100 NW=0（半带宽）DO 110 IE=1,NEDO 110 I=1,3DO 110 J=1,3IW=IABS(LND(IE,I)-LND(IE,J))IF(NW.LT.IW) THENNW=IWENDIF110 CONTINUENW=(NW+1)*2IF(IPS.NE.0) THENE=E/(1.0-PR*PR)PR=PR/(1.0-PR)ENDIFEND*---------------------------------------------------------C SUBPROGRAM-2C CALCULATE ELEMENT STIFFNESS MATRIXSUBROUTINE STE(IE,NJ,NE,LND,X,Y,E,PR,T,KE)DIMENSION LND(NE,3),X(NJ),Y(NJ),B(3,6),D(3,3)REAL KE(6,6)CALL ATE(IE,NJ,NE,LND,X,Y,AE)CALL DTE(E,PR,D)CALL BTE(IE,NJ,NE,LND,X,Y,AE,B)DO 10 I=1,6DO 10 J=1,6KE(I,J)=0.DO 10 K=1,3DO 10 K1=1,310 KE(I,J)=KE(I,J)+B(K,I)*D(K,K1)*B(K1,J)C=AE*TDO 30 I=1,6DO 30 J=1,630 KE(I,J)=KE(I,J)*CEND*------------------------------------------------ C SUBPROGRAM-3C CALCULATE ELEMENT AREASUBROUTINE ATE(IE,NJ,NE,LND,X,Y,AE)DIMENSION LND(NE,3),X(NJ),Y(NJ)I=LND(IE,1)J=LND(IE,2)K=LND(IE,3)XIJ=X(J)-X(I)YIJ=Y(J)-Y(I)XIK=X(K)-X(I)YIK=Y(K)-Y(I)AE=.5*(XIJ*YIK-XIK*YIJ)END*---------------------------------------------- C SUBPROGRAM-4C CALCULATE ELASTICITY MATRIXSUBROUTINE DTE(E,PR,D)DIMENSION D(3,3)DO 10 I=1,3DO 10 J=1,310 D(I,J)=0.D(1,1)=E/(1.-PR*PR)D(1,2)=E*PR/(1.-PR*PR)D(2,1)=D(1,2)D(2,2)=D(1,1)D(3,3)=.5*E/(1.+PR)END*------------------------------------------------ C SUBPROGRAM-5C CALCULATE MATRIX [B]SUBROUTINE BTE(IE,NJ,NE,LND,X,Y,AE,B)DIMENSION LND(NE,3),X(NJ),Y(NJ),B(3,6)I=LND(IE,1)J=LND(IE,2)K=LND(IE,3)DO 10 II=1,3DO 10 JJ=1,610 B(II,JJ)=0.B(1,1)=Y(J)-Y(K)B(1,3)=Y(K)-Y(I)B(1,5)=Y(I)-Y(J)B(2,2)=X(K)-X(J)B(2,4)=X(I)-X(K)B(2,6)=X(J)-X(I)B(3,1)=B(2,2)B(3,2)=B(1,1)B(3,3)=B(2,4)B(3,4)=B(1,3)B(3,5)=B(2,6)B(3,6)=B(1,5)DO 20 I1=1,3DO 20 J1=1,620 B(I1,J1)=.5/AE*B(I1,J1)END*------------------------------------------------------- C SUBPROGRAM-6C CALCULATE GLOBAL STIFFNESS MATRIXSUBROUTINE STIFF(NJ,NE,NJ2,NW,LND,X,Y,E,PR,T,KS) DIMENSION LND(NE,3),X(NJ),Y(NJ)REAL KS(NJ2,NW),KE(6,6)DO 5 I=1,NJ2DO 5 J=1,NW5 KS(I,J)=0.DO 10 IE=1,NECALL STE(IE,NJ,NE,LND,X,Y,E,PR,T,KE)DO 10 I=1,3IZ=LND(IE,I)DO 10 II=1,2IH =2*(I -1)+IIIDH=2*(IZ-1)+IIDO 10 J=1,3JZ=LND(IE,J)DO 10 JJ=1,2L =2*(J -1)+JJIL=2*(JZ-1)+JJIF(IL.GE.IDH) THENIDL=IL-IDH+1KS(IDH,IDL)=KS(IDH,IDL)+KE(IH,L)ENDIF10 CONTINUEEND*-------------------------------------------------------- C SUBPROGRAM-7C CALCULATE NODAL LOAD VECTORSUBROUTINE EQUPE(NJ,NE,NPJ,NJ2,T,V,LND,X,Y,PJ,P) DIMENSION LND(NE,3),X(NJ),Y(NJ),PJ(NPJ,3),P(NJ2) DO 10 I=1,NJ210 P(I)=0.DO 20 I=1,NPJII=PJ(I,1)P(2*II-1)=PJ(I,2)20 P(2*II)=PJ(I,3)30 IF(V.EQ.0.) GOTO 50DO 40 IE=1,NECALL ATE(IE,NJ,NE,LND,X,Y,AE)PE=-V*AE*T/3.DO 40 I=1,3II=LND(IE,I)40 P(2*II)=P(2*II)+PE50 RETURNEND*---------------------------------------------C SUBPROGRAM-8C INTRODUCE BOUNDARY CONDITIONSUBROUTINE INSCD(NS,NW,NJ2,JR,KS,P)DIMENSION P(NJ2),JR(NS,3)REAL KS(NJ2,NW)DO 30 I=1,NSIR=JR(I,1)DO 30 J=2,3IF(JR(I,J).EQ.0) GOTO 30II=2*IR+J-3KS(II,1)=1.DO 10 JJ=2,NW10 KS(II,JJ)=0.IF(II.GT.NW) JO=NWIF(II.LE.NW) JO=IIDO 20 JJ=2,JO20 KS(II-JJ+1,JJ)=0.P(II)=0.30 CONTINUEEND*-------------------------------------------C SUBPROGRAM-9C SOLVE EQUATIONS BY GS METHODSUBROUTINE BGSMT(NJ,NJ2,NW,KS,P)DIMENSION P(NJ2)REAL KS(NJ2,NW)NJ1=NJ2-1DO 20 K=1,NJ1IF(NJ2.GT.K+NW-1) IM=K+NW-1IF(NJ2.LE.K+NW-1) IM=NJ2K1=K+1DO 20 I=K1,IML=I-K+1C=KS(K,L)/KS(K,1)IW=NW-L+1DO 10 J=1,IWM=J+I-K10 KS(I,J)=KS(I,J)-C*KS(K,M)20 P(I)=P(I)-C*P(K)P(NJ2)=P(NJ2)/KS(NJ2,1)DO 40 I1=1,NJ1I=NJ2-I1IF(NW.GT.NJ2-I+1) JO=NJ2-I+1IF(NW.LE.NJ2-I+1) JO=NWDO 30 J=2,JOK=J+I-130 P(I)=P(I)-KS(I,J)*P(K)40 P(I)=P(I)/KS(I,1)WRITE(6,50)50 FORMAT(/1X,'NODAL DISPLACEMENTS'/3X,* 'NODE',5X,'X-DISP.',8X,'Y-DISP.')DO 60 I=1,NJ60 WRITE(6,70) I,P(2*I-1),P(2*I)70 FORMAT(1X,I5,2E15.6)END*--------------------------------------------------- C SUBPROGRAM-10C CALCULATE ELEMENT STRESS MATRIXSUBROUTINE SIGME(NE,NJ,NJ2,E,PR,LND,X,Y,P)DIMENSION LND(NE,3),X(NJ),Y(NJ),D(3,3),B(3,6), * S(3,6),ST(3),P(NJ2),DE(6)WRITE(6,*)WRITE(6,*)' ELEMENT STRESSES'CALL DTE(E,PR,D)DO 50 IE=1,NECALL ATE(IE,NJ,NE,LND,X,Y,AE)CALL BTE(IE,NJ,NE,LND,X,Y,AE,B)DO 10 I=1,3DO 10 J=1,6S(I,J)=0.DO 10 K=1,310 S(I,J)=S(I,J)+D(I,K)*B(K,J)DO 20 I=1,3DO 20 J=1,2IH=2*(I-1)+JIW=2*(LND(IE,I)-1)+J20 DE(IH)=P(IW)DO 30 I=1,3ST(I)=0.DO 30 J=1,630 ST(I)=ST(I)+S(I,J)*DE(J)SGX=ST(1)SGY=ST(2)TXY=ST(3)ASG=(SGX+SGY)*.5RSG=SQRT(.25*(SGX-SGY)**2+TXY*TXY)SGMA=ASG+RSGSGMI=ASG-RSGIF(SGY.EQ.SGMI) CETA=0.IF(SGY.NE.SGMI) CETA=90.-57.29578*ATAN* (TXY/(SGY-SGMI))50 WRITE(6,60) IE,SGX,SGY,TXY,SGMA,SGMI,CETA60 FORMAT(1X,'ELEMENT NO.=',I4/2X,'SIGX=',E10.4, * 2X,'SIGY=',E10.4,2X,'TXY =',E10.4/2X,'SGMA=', * E10.4,2X,'SGMI=',E10.4,2X,'CETA=',E10.4)END。
8 第四章 用常应变三角形单元解力学平面问题 (2)解析
um xm ym
1 um ym
1 xm um
其中
1 xi yi
2 1 x j y j
1 xm ym
(c) (d) (1)
从解析几何可知，式中的 就是三角形i、j、m的面积。
为保证求得的面积为正值，节点i、j、m的编排次序必须是逆 时针方向，如图1所示。
7. 由单元的节点位移列阵计算单元应力
解出整体结构的节点位移列阵 后，再根据单元节点的 编号找出对应于单元的位移列阵 e，将 e代入（3-3）式就
可求出各单元的应力分量值。
8. 计算结果输出
求解出整体结构的位移和应力后，可有选择 地整理输出某些关键点的位移值和应力值，特别 要输出结构的 变形图、应力图、应变图、结构仿 真变形过程动画图及整体结构的弯矩、剪力图等 等。
平面问题可用三角元，四边元等。
例如：
3. 选择单元的位移模式
结构离散化后，要用单元内结点的位移通过插值来获得 单元内各点的位移。在有限元法中，通常都是假定单元的位 移模式是多项式，一般来说，单元位移多项式的项数应与单 元的自由度数相等。它的阶数至少包含常数项和一次项。至 于高次项要选取多少项，则应视单元的类型而定。
有限元法的实质是：把有无限个自由度的连续体， 理想化为只有有限个自由度的单元集合体，使问题简化 为适合于数值解法的结构型问题。
二、经典解与有限元解的区别：
微分 经 典 解 法 —— （解析法）
数目增到∞ 大小趋于 0
建立一个描述连续体 性质的偏微分方程
有限单元 离散化 集合
总体分析解
有限元法——连续体——单元——代替原连续体
式中：
Re ke e
（3-4）
——单元刚度矩阵
ke BT DBdxdydz
第一节三角形常应变单元
第三章平面问题的有限元法本章通过三角形常应变单元，介绍有限元法应用于弹性体应力分析的基本原理和方法：包括弹性体的离散化，单元特性的分析，刚度矩阵的建立，等效节点力的计算，解答的收敛性以及实施步骤和注意事项，同时对形函数的性质也作简要的叙述。
第一节三角形常应变单元一、结构离散化用有限元法分析弹性力学平面问题，第一步就是把原来的连续的弹性体离散化。
(a) (b)将整个结构（平板）划分成有限个三角形。
这样的三角形称为单元（三角形单元）。
三角形单元的顶点取为节点。
3节点三角形单元用边界节点间的直线段来近似板的曲线边界。
这些三角形在其节点处相互连接，组成一个单元集合体，以代替原来的弹性体。
注：1. 全部节点和全部单元一般由1开始按自然顺序编号。
2. 节点编码：总码-----------用于整体分析，如1,2，…，按自然顺序编号局部码--------用于单元分析，i，j，m 要求按逆时针方向的顺序进行编码每个单元的节点局部码i，j，m和节点总码有一一对应的关系3. 单元间不能有重叠4. 一个单元的任一顶点不许为另一单元任一边的内点5. 所有作用在单元上的载荷，包括集中载荷、表面载荷和体积力，都按虚功等效的原则移置到节点上，成为等效节点载荷。
二、位移模式1. 单元节点位移列阵iu设单元e 的节点号码为i ，j ，m 。
由弹性力学平面问题可知，单元内任意一点有两个位移分量u ，v ，记为{}Tf uv ⎡⎤⎣⎦=故每个节点也有两个位移分量，因此称节点自由度为2。
节点i 位移分量记为{}Ti ii u v δ⎡⎤⎣⎦= （i ，j ，m 轮换）则3个节点的位移分量用列阵表示为{}i ei i e j j j m m m u v u v u v δδδδ⎧⎫⎪⎪⎪⎪⎧⎫⎪⎪⎪⎪⎪⎪⎪⎪⎨⎬⎨⎬⎪⎪⎪⎪⎪⎪⎪⎪⎩⎭⎪⎪⎪⎪⎩⎭==（3-1）称为单元节点位移列阵（向量）。
单元自由度是6。
2. 位移模式结构承受载荷以后发生位移，但位移分布事先并不知道。
三角形常应变单元与二维问题求解
5.1 引言 二维单元适用于在平面内承受载荷和 变形的结构问题。
虽然实际结构不可能是真正二维的，而有经验的分 析人员通常可以把许多实际问题抽象为而为问题并采用二 维模型来进行分析，从而得到满意的结果，这比采用三维 分析更有效和经济。
平面应力：物体的厚度与其他几何尺寸相比小
nd 1 i j , j 1, 2,..., 0 i j ,i , j 1, 2,...,ttention!
N ( x) 1
i 1 i
n
（2） 推导应变矩阵 B （3） 利用矩阵N 和 B 计算单元刚度矩阵、 单元质量矩阵、单元力向量。
5.2 有限元模型 对于复杂几何适应性而言，三角形单元非常 实用，但精度不高。
5.3 常应变三角形单元
根据形函数单位分解的性质：
Ni ( x j )
得多，所有的力都作用在 x y平面内，因此位移 只是 x 和 y 的函数。
平面应变:物体的厚度与结构的其他几何尺寸相
比大很多，并且载荷（压力）沿长度方向是均匀 分布，因此，沿其在整个厚度上的变形几乎是相 同的。
求解步骤： （1） 构造满足下式的形函数矩阵 N
Ni ( x j )
nd 1 i j , j 1, 2,..., 0 i j ,i , j 1, 2,..., nd
第7章——有限元三角形单元程序设计
C ( 4, 1)
B(0, 2)
12:45
有限单元法
崔向阳
三角形单元程序
clear all first_time= first_time =cputime cputime; ; format long g %-------------------------------------------%input data for control parameters %-------------------------------------------lengthx=4; lengthx =4; %length of xx-axis side of problem lengthy=2; %length of yy-axis side of problem emodule=1.0; emodule=1.0; ; poisson=0.0; poisson =0.0; fload= fload =-1; %elastic modulus %Poisson's ratio % the total load
B(0, 2)
A(0, 0) D (4, 0)
F 1
C ( 4, 4 1)
12:45
有限单元法
崔向阳
三角形单元程序
lx=16; % number of element in x x-axis ly=8; ly =8; % number of element in yy-axis nel=2*lx* nel =2*lx*ly ly; ; % number of element nnel=3; nnel =3; %number of nodes per element ndof=2; ndof =2; %number of dofs per node nnode=(lx+1)*(ly+1); nnode =(lx+1)*(ly+1); (lx 1) (ly 1); %total number of nodes in system sdof= sdof =nnode* nnode*ndof; ndof; %total system dofs edof= edof =nnel* nnel*ndof; ndof; %degrees of freedom per element %------------------------------------------------------%input data for nodal coordinate values %------------------------------------------------------x0=[]; for i=1:lx+1 for j j=1:ly+1 1:ly+1 x0=[x0; (i(i-1)* 1)*lengthx lengthx/lx /lx end end 12:45
平面三角形单元常应变单元matlab程序的编制
———————————————————————————————— 作者：
———————————————————————————————— 日期：
ﻩ
三角形常应变单元程序的编制与使用
有限元法是求解微分方程边值问题的一种通用数值方法,该方法是一种基于变分法(或变分里兹法）而发展起来的求解微分方程的数值计算方法，以计算机为手段，采用分片近似，进而逼近整体的研究思想求解物理问题。
有限元法中三节点三角形分析结构的步骤如下:
1）整理原始数据,如材料性质、荷载条件、约束条件等,离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。
2）单元分析,建立单元刚度矩阵。
3）整体分析，建立总刚矩阵。
4）建立整体结构的等效。
POISＳ=fsｃanf（FＰ1，'%f'，1）%泊松比
TＨICＫ=ｆscａnf（ＦＰ1,'%d',1)％厚度
LＮODS=fscanf（ＦＰ1,'%ｄ',［3,NＥLEM］)'%单元定义数组(单元结点号）
功能:运用有限元法中三角形常应变单元解平面问题的计算主程序。
基本思想：单元结点按右手法则顺序编号。
荷载类型：可计算结点荷载。
说明:主程序的作用是通过赋值语句、读取和写入文件、函数调用等完成算法的全过程，即实现程序流程图的程序表达。
%---－--－-----－---－-----－--－---－－---－--－-----－-－-－－－－-－-－－-－-----------------－----－----－----－----------
glｏｂａｌFＯRCＥFIXＥＤ
globaｌＢMATＸＤMＡTＸSＭATＸAＲＥA
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
三角形常应变单元程序的编制与使用有限元法是求解微分方程边值问题的一种通用数值方法，该方法是一种基于变分法（或变分里兹法）而发展起来的求解微分方程的数值计算方法，以计算机为手段，采用分片近似，进而逼近整体的研究思想求解物理问题。
有限元分析的基本步骤可归纳为三大步：结构离散、单元分析和整体分析。
对于平面问题，结构离散常用的网格形状有三角形、矩形、任意四边形，以三个顶点为节点的三角形单元是最简单的平面单元，它较矩形或四边形对曲边边界有更好的适应性，而矩形或四边形单元较三节点三角形有更高的计算精度。
Matlab语言是进行矩阵运算的强大工具，因此，用Matlab语言编写有限元中平面问题的程序有优越性。
本章将详细介绍如何利用Matlab语言编制三角形常应变单元的计算程序，程序流程图见图1。
有限元法中三节点三角形分析结构的步骤如下：1）整理原始数据，如材料性质、荷载条件、约束条件等，离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。
2）单元分析，建立单元刚度矩阵。
3）整体分析，建立总刚矩阵。
4）建立整体结构的等效节点荷载和总荷载矩阵5）边界条件处理。
6）解方程，求出节点位移。
7）求出各单元的单元应力。
8）计算结果整理。
计算结果整理包括位移和应力两个方面；位移计算结果一般不需要特别的处理，利用计算出的节点位移分量，就可画出结构任意方向的位移云图；而应力解的误差表现在单元内部不满足平衡方程，单元与单元边界处应力一般不连续，在边界上应力解一般与力的边界条件不相符合。
图1 程序流程图1.1 程序说明%******************************************************************* % 三角形常应变单元求解结构主程序%******************************************************************* ●功能：运用有限元法中三角形常应变单元解平面问题的计算主程序。
●基本思想：单元结点按右手法则顺序编号。
●荷载类型：可计算结点荷载。
●说明：主程序的作用是通过赋值语句、读取和写入文件、函数调用等完成算法的全过程，即实现程序流程图的程序表达。
%----------------------------------------------------------------------------------------------------- 1 程序准备format short e %设定输出类型clear all %清除所有已定义变量clc %清屏●说明：format short e －设定计算过程中显示在屏幕上的数字类型为短格式、科学计数法；clear all －清除所有已定义变量，目的是在本程序的运行过程中，不会发生变量名相同等可能使计算出错的情况；clc －清屏，使屏幕在本程序运行开始时%----------------------------------------------------------------------------------------------------- 2 全局变量定义global NNODE NPION NELEM NVFIX NFORCE COORD LNODS YOUNG POISS THICKglobal FORCE FIXEDglobal BMATX DMATX SMATX AREAglobal ASTIF ASLOD ASDISPglobal FP1●说明：NNODE—单元结点数，NPION—总结点数，NELEM—单元数，NVFIX—受约束边界点数，NFORCE—结点力数，COORD—结构结点坐标数组，LNODS —单元定义数组，YOUNG—弹性模量，POISS—泊松比，THICK—厚度FORCE —节点力数组(n,3) n:受力节点数目,(n,1):作用点,(n,2):x方向,(n,3):y 方向；FIXED—约束信息数组(n,3) n:受约束节点数目, (n,1):约束点(n,2)与(n,3)分别为约束点x方向和y方向的约束情况,受约束为1否则为0BMATX—单元应变矩阵(3*6)，DMATX—单元弹性矩阵(3*3)，SMATX—单元应力矩阵(3*6)，AREA—单元面积ASTIF—总体刚度矩阵，ASLOD—总体荷载向量，ASDISP—结点位移向量FP1—数据文件指针3 打开文件FP1=fopen('input.txt','rt'); %打开输入数据文件存放初始数据●说明：FP1=fopen('input.txt','rt'); －打开已存在的输入数据文件input.txt，且设置其为只读格式，使程序在执行过程中不能改变输入文件中的数值，并用文件句柄FP1来执行FP2=fopen('output.txt','wt'); －打开输出数据文件，该文件不存在时，通过此命令创建新文件，该文件存在时则将原有内容全部删除。
该文件设置为可写格式，可在程序执行过程中向输出文件写入数据。
%----------------------------------------------------------------------------------------------------- 4 读入程序控制信息NPION=fscanf(FP1,'%d',1) %结点个数（结点编码总数）NELEM=fscanf(FP1,'%d',1) %单元个数（单元编码总数）NFORCE=fscanf(FP1,'%d',1) ％结点荷载个数NVFIX=fscanf(FP1,'%d',1) ％受约束边界点数YOUNG=fscanf(FP1,'%e',1) ％弹性模量POISS=fscanf(FP1,'%f',1) ％泊松比THICK=fscanf(FP1,'%d',1) %厚度LNODS=fscanf(FP1,'%d',[3,NELEM])' %单元定义数组（单元结点号）●说明：建立LNODS矩阵，该矩阵指出了每一单元的连接信息。
矩阵的每一行针对每一单元，共计NELEM；每一列相应为单元结点号（编码）、按逆时针顺序输入。
命令中，[3,NELEM]’表示读取NELEM行3列数据赋值给LNODS矩阵。
显然，LNODS(i,1:3)依次表示i单元的i，j，k结点号。
COORD=fscanf(FP1,'%f',[2,NPION])' %结点坐标数组●说明：建立COORD矩阵，该矩阵用来存储各结点x，y方向的坐标值。
从FP1文件中读取全部结点个数NPOIN的坐标值，从1开始按顺序读取。
COORD(i，1：2)表示第i个结点的x，y坐标。
FORCE=fscanf(FP1,'%f',[3,NFORCE])' %结点力数组●说明：(n,3) n:受力结点数目,(n,1):作用点,(n,2):x方向,(n,3):y方向FIXED=fscanf(FP1,'%d',[3,inf])' %约束信息数组●说明：(n,3) n:受约束节点数目, (n,1):约束点(n,2)与(n,3)分别为约束点x方向和y 方向的约束情况,受约束为1否则为0●总体说明：从输入文件FP1中读入结点个数，单元个数，结点荷载个数，受约束边界点数，弹性模量，泊松比，厚度，单元定义数组，结点坐标数组，结点力数组，约束信息数组；程序中弹性模量仅输入了一个值，表明本程序仅能求解一种材料构成的结构，如：钢筋混凝土结构、钢结构，不能求解钢筋混凝土－钢组合结构。
采用了命令fscanf，其中：’%d’表示读入整数格式，’%f'’表示读入浮点数；1表示读取1个数，[A,B]形式表示读A行B列数组，[A,B]’表示将[A,B]转置，inf 表示正无穷。
%----------------------------------------------------------------------------------------------------- 5 调用子程生成单刚，组成总刚并加入约束信息function ASSEMBLE()%----------------------------------------------------------------------------------------------------- 6 调用子程生成荷载向量function FORMLOAD()%----------------------------------------------------------------------------------------------------- 7 计算结点位移向量ASDISP=ASTIF\ASLOD'%----------------------------------------------------------------------------------------------------- 8调用子程计算单元应力function WRITESTRESS()%******************************************************************* 9 关闭输出数据文件fclose(FP2);%*******************************************************************读取ASSEMBLE子程%*******************************************************************function ASSEMBLE()% 所引用的全局变量：global NPION NELEM NVFIX LNODS ASTIF THICKglobal BMATX SMATX AREA FIXED%------------------------------------------------------------------------------------------------ ---- % 计算单刚并生成总刚ASTIF(1:2*NPION,1:2*NPION)=0; %张成特定大小总刚矩阵并置0 ●说明：建立单元刚度矩阵ASTIF，该矩阵的行列数均为2*NPION ，NPION表示结点数，每个结点有两个方向的力和位移。