红-计算力学框架编程matlab

1、计算力学框架编程（）班级：号：期：计力（春）许 映 红 1 4 2 5 6 42015-06学日目录第一章题 目1题目选择1题目解读1空间框架1编程对象3第二章程 序4程序体系4基本设想4主程序5函数5程序代码6主程序代码6函数代码19第三章算 例28基本信息28运行结果29运行截图29其他截图37结果文件37第一章题 目1.1 题目选择编制框架结构在外荷载作用下结构的应力、位移等响应。1)2)3)4)5)Level 1（ Level 2（ Level 3（Level 4（60%）：平面框架+静力荷载；70%）：平面框架+动力荷载；80%）：空间框架+静力荷载；90%）：空间框架+动力荷载；L

2、evel 5（100%）：空间框架+静力荷载+动力荷载。选择“Level 3（ 80%）：空间框架+静力荷载”进行编制。1.2 题目解读空间框架采用空间三维有限元，一般用考虑了横向剪切变形的 Timoshenko 梁理论进行考虑。这种梁单元既有局部坐标又有总体坐标。1.2.1 空间框架空间杆单元和平面杆单元的区别在于：杆件除了可能承受轴力和弯矩的作用外，还可能承受扭矩的作用；而且弯矩可能同时在两个坐标面内存在。空间杆单元的系数有弹性模量 E 、剪切弹性模量G 、横截面积 A 、惯性矩 I x 和 I y 、极惯性矩 J 及长度 L 。每个空间杆单元有 2 个节点，并且相对于总体坐标系的 X、Y

3、、Z 轴正向逆时针的倾斜角分别为xX 、yY 、zZ 。如图 1.1 所示。图 1.1 空间杆单元1因此，单元刚度矩阵可根据下面的矩阵给出： = 式中，旋转矩阵T 和弹性模量矩阵k定义如下：（1-1） = 1112（1-2）21221 01 00020002001 02 00002 0001 00020002 000 20230000 = （1-2-1）11 002310002001 02 000020002 00 = 0 0 000 0 = 2（1-2-2）122103 0000 3 1 001 02 0 0001000 20000 = （1-2-3）22 230 0000223式中， =

4、， = ，1262126=， =， =，=，=121231232322=。3 = 000 = 式中， C0 0000 0000（1-3）（1-3-1），coxY ，coxZ ，等等。xYxZ很显然，空间杆单元有 12 个度每个节点有 6 个度（3 个线位移和 3 个角位移）。因此，对一个有 n 个节点的空间框架结构而言，其整体刚度矩阵 K 是 6n6n 的（因为每个节点有 6 个度）。一旦得到整体刚度矩阵 K ，就可以列出以下方程组：2（1-4） = 式中，U 是总体节点位移矢量，F 是总体节点载荷矢量。在这一步中，边界条件被手动赋值给矢量U 和 F 。然后用分解和消去法求解矩阵方程（1-4）

5、。一旦求出未知的位移和支反力，就可以用下式求出每个单元的节点力矢量 f ： = （1-5）式中， f 是 61 的节点力矢量、u 是 121 的单元节点位移矢量。单元位移矢量u的第 1 个、第 2 个和第 3 个元素是第一个节点的位移分量，第 4 个、第 5 个和第 6 个元素是第一个节点的转角分量。单元位移矢量u的第7 个、第 8 个和第 9 个元素是第二个节点的位移分量，第 10 个、第 11 个和第 12个元素是第二个节点的转角分量。1.2.2 编程对象本程序采用进行编程，主要针对的满足以下要求的空间框架进行编制：1.2.3.4.5.水平为 X 层，竖向为 Y 层，横向为 Z 层；水平杆

6、长度为 L（m），竖杆长度为 H（m），横杆长度为 W（m）；杆件截面面积 A（m2），截面惯性矩 Iy（m4）、Iz（m4）,扭转惯性矩J（m4）；杆件材料弹性模量 E（MPa），剪切模量 E（MPa）；顶层节点受 Y 向荷载 P（kN），底层节点固结。其中部分函数可适用于任意形式空间框架计算。3第二章程 序2.1 程序体系2.1.1 基本设想本程序需要达到的目标有：1.模型基本参数输入、节点、单元、荷载施加、约束施加以及模型的三维显示；对框架的节点及单元如图 2.13 所示：1)节点以平行于 XY 面为，再以平行于 X 轴为， X 向增大，Y 向增大，Z 向增大。图 2.1 节点示意图2)

7、单元先编平行于 Y 轴的竖杆，以平行于 XY 面为向增大，Y 向增大，Z 向增大。，再编较低层， X图 2.2 竖杆单元示意图3)单元后编垂直于 Y 轴的水平杆和横杆，以平行于 XY 面为，再编较4低层， X 向增大，Y 向增大，Z 向增大。图 2.3 水平杆、横杆单元示意图2. 刚度矩阵计算、位移及力计算以及计算结果的保存与输出。关键在于空间框架的刚度矩阵计算。2.1.2 主程序主程序共分为 7 大模块，分别为建立模型；集成总刚；求解位移；求解内力；保存结果；输出结果；绘制内力图。见“源代码”文件夹中的“主程序”。2.1.3 函数主程序中主要用到 3 个关键计算用函数以及 6 个绘图用函数，

8、分别为：计算单元刚度矩阵；总体刚度矩阵；求解单元节点力；绘制单元轴力图；绘制单元剪力（Y 轴）图；5绘制单元剪力（Z 轴）图；绘制单元扭矩图；绘制单元弯矩（Y 轴）图；绘制单元弯矩（Z 轴）图。见“源代码”文件夹中的“函数”。2.2 程序代码2.1.1 主程序代码clear建立模型%*%1.建立模型%*%模型基本参数fprf(请输入模型参数n);X=input(水平向层数：);Y=input(竖向层数：);Z=input(横向层数：);H=input(竖杆长度（m）：);L=input(水平杆长度（m）：);W=input(横杆长度（m）：);A=input(杆件截面面积（m2）：);Iy=i

9、nput(在 xz 面内截面惯性矩（m4）：);Iz=input(在 xy 面内截面惯性矩（m4）：);J=input(扭转惯性矩（m4）：);E=input(材料弹性模量（MPa）：);G=input(材料剪切模量（MPa）：);P=input(节点力 N_y（kN）：);%定义节点坐标6Node=zeros(X+1)*(Y+1)*(Z+1),3);x=0;y=0;z=0;t1=1;for s=1:(Z+1)z=W*(s-1);for j=1:(Y+1)y=H*(j-1);for i=1:(X+1)x=L*(i-1);Node(t1,:)=xy z;t1=t1+1;endendend%节点的

10、总数量，用 t1n 表示t1n=t1-1;fprf(n“节点建立”，完成！n);%定义单元Element=zeros(3*X*Z+2*(X+Z)+1)*Y),2);%定义平行于 y 轴的竖杆单元t1=1;t2=1;for s=1:(Z+1)for j=1:(Y+1)for i=1:(X+1)if j1&i1&i=(X+1)Element(t2,:)=t1(t1+(X+1)*(Y+1);t2=t2+1;endt1=t1+1;endendendfor j=1:(Y+1)for i=1:(X+1)if j1&iziLambda=0 0 1;0 1 0;-1 0 0;elseLambda=0 0 -1

11、;0 1 0;1 0 0;endelseCxX=(xj-xi)/L;CxY=(yj-yi)/L;20CxZ=(zj-zi)/L;D=sqrt(CxX*CxX+CxY*CxY);CyX=-CxY/D;CyY=CxX/D;CyZ=0;CzX=-CxX*CxZ/D;CzY=-CxY*CxZ/D;CzZ=D;Lambda=CxX CxY CxZ;.CyX CyY CyZ;.CzX CzY CzZ;endT=Lambda zeros(3) zeros(3) zeros(3);.zeros(3)Lambda zeros(3) zeros(3);.zeros(3)zeros(3) Lambda zeros(

12、3);.zeros(3)zeros(3) zeros(3) Lambda;%计算整体坐标系下的单元刚度矩阵y=T*kprime*T;end总体刚度矩阵%集成总体刚度矩阵function y=Assemble(K,k,i,j)%将第 16 行数据加入总刚for s=1:6%将第 s 行数据加入总刚t=6-s;K(6*i-t,6*i-5)=K(6*i-t,6*i-5)+k(s,1);K(6*i-t,6*i-4)=K(6*i-t,6*i-4)+k(s,2);K(6*i-t,6*i-3)=K(6*i-t,6*i-3)+k(s,3);K(6*i-t,6*i-2)=K(6*i-t,6*i-2)+k(s,4

13、);21K(6*i-t,6*i-1)=K(6*i-t,6*i-1)+k(s,5);K(6*i-t,6*i)=K(6*i-t,6*i)+k(s,6);K(6*i-t,6*j-5)=K(6*i-t,6*j-5)+k(s,7);K(6*i-t,6*j-4)=K(6*i-t,6*j-4)+k(s,8);K(6*i-t,6*j-3)=K(6*i-t,6*j-3)+k(s,9);K(6*i-t,6*j-2)=K(6*i-t,6*j-2)+k(s,10);K(6*i-t,6*j-1)=K(6*i-t,6*j-1)+k(s,11);K(6*i-t,6*j)=K(6*i-t,6*j)+k(s,12);end%将

14、第 712 行数据加入总刚fors=7:12%将第 s 行数据加入总刚t=12-s;K(6*j-t,6*i-5)=K(6*j-t,6*i-5)+k(s,1);K(6*j-t,6*i-4)=K(6*j-t,6*i-4)+k(s,2);K(6*j-t,6*i-3)=K(6*j-t,6*i-3)+k(s,3);K(6*j-t,6*i-2)=K(6*j-t,6*i-2)+k(s,4);K(6*j-t,6*i-1)=K(6*j-t,6*i-1)+k(s,5);K(6*j-t,6*i)=K(6*j-t,6*i)+k(s,6);K(6*j-t,6*j-5)=K(6*j-t,6*j-5)+k(s,7);K(6

15、*j-t,6*j-4)=K(6*j-t,6*j-4)+k(s,8);K(6*j-t,6*j-3)=K(6*j-t,6*j-3)+k(s,9);K(6*j-t,6*j-2)=K(6*j-t,6*j-2)+k(s,10);K(6*j-t,6*j-1)=K(6*j-t,6*j-1)+k(s,11);K(6*j-t,6*j)=K(6*j-t,6*j)+k(s,12);endy=K;end求解单元节点力%计算单元节点力矢量22function y=ElementForce(E,G,A,Iy,Iz,J,xi,yi,zi,xj,yj,zj,u)%计算杆件长度，本程序杆件长度均为 L 或 H 或 WL=(xj

16、-xi)2+(yj-yi)2+(zj-zi)2)(1/2);%计算中的非零项a1=E*A/L;a2=G*J/L;b1=12*E*Iz/(L3);b2=6*E*Iz/(L2);b3=2*E*Iz/L;c1=12*E*Iy/(L3);c2=6*E*Iy/(L2);c3=2*E*Iy/L;%计算局部坐标系下的单元刚度矩阵kprime=a1 0 0 0 0 0 -a1 0 0 0 0 0;.0b1 0 00 b2 0 -b1 0 0 0 b2;.00c1 0-c2 000-c1 0 -c2 0;.000 a20 0 000-a2 0 0;.00-c2 0 2*c3000 c2 0 c3 0;.0b2

17、0 0 0 2*b3 0 -b2 0 0 0 b3;.-a1 0 000 0 a1 0 0 0 0 0;.0-b1 000 -b2 0 b1 0 0 0 -b2;.00-c10c2 0 0 0 c1 0 c2 0;.000 -a20 0 0 0 0 a2 0 0;.00-c2 0c3 0 0 0 c2 0 2*c3 0;.0b2 0 0 0 b3 0 -b2 0 0 0 2*b3;%计算由局部坐标系 x、y、z 轴转换到整体坐标系 X、Y、Z 轴的坐标转换矩阵 Tif xi=xj&yi=yjif zjziLambda=0 0 1;0 1 0;-1 0 0;else23Lambda=0 0 -1

18、;0 1 0;100;endelseCxX=(xj-xi)/L;CxY=(yj-yi)/L;CxZ=(zj-zi)/L;D=sqrt(CxX*CxX+CxY*CxY);CyX=-CxY/D;CyY=CxX/D;CyZ=0;CzX=-CxX*CxZ/D;CzY=-CxY*CxZ/D;CzZ=D;Lambda=CxX CxY CxZ;.CyX CyY CyZ;.CzX CzY CzZ;endT=Lambda zeros(3) zeros(3) zeros(3);.zeros(3)Lambda zeros(3) zeros(3);.zeros(3)zeros(3) Lambda zeros(3);.

19、zeros(3)zeros(3) zeros(3) Lambda;%计算整体坐标系下的单元刚度矩阵y=kprime*T*u;end绘制单元轴力图%绘制节点力矢量沿杆长变化的轴向力曲线图function y=ElementAxialDiagram(f,)x=0;L;z=-f(1);f(7);hold on;24title(轴力图-Element ,num2str(ie);xlabel(L);ylabel(Nx);plot(x,z);y1=0;0;plot(x,y1,k);end绘制单元剪力（Y 轴）图%绘制节点力矢量沿杆长变化的剪力曲线图（沿 Z 轴）function y=ElementShea

20、rYDiagram(f,)x=0;L;z=f(2);-f(8);hold on;title(剪力图(Y)-Element ,num2str(ie);xlabel(L);ylabel(Ny);plot(x,z);y1=0;0;plot(x,y1,k);end绘制单元剪力（Z 轴）图%绘制节点力矢量沿杆长变化的剪力曲线图（沿 Z 轴）function y=ElementShearZDiagram(f,)x=0;L;z=f(3);-f(9);hold on;title(剪力图(Z)-Element ,num2str(ie);xlabel(L);ylabel(Nz);plot(x,z);25y1=0;0;plot(x,y1,k);end绘制单元扭矩图%绘制节点力矢量沿杆长变化的扭矩曲线图function y=ElementTorDiagram(f,)x=0;L;z=f(4);-f(10);hold on;title(扭矩图-Element ,num2str(ie);xlabel(L);ylabel(Mx);plot(x,z);y1=0;0;plot(x,y1,k);end绘制单元弯矩（Y 轴）图%绘制节点力矢量沿杆长变化的弯矩曲线图（沿 Y 轴）functi