计算力学(有限元)matlab编程大作业.docx

1、 计算力学编程大作业（网架结构）一、 编程算法： 图11.总刚矩阵形成：根据图1，对于任意的i,j节点，不论原来的连接状况如何，只要I,j两节点有了联系，则在总刚矩阵中对应的位置需要加上有ij杆的产生带来的刚度贡献，贡献具体可表现为在Kii处需要加上有ij杆件贡献的K11（经过坐标变换后），Kjj处需要加上ij杆件贡献的K22（经过坐标变换后），Kij处需要加上有ij杆件贡献的K12（经过坐标变换后），Kji处需要加上有ij杆件贡献的K21（经过坐标变换后）,由此可以直接形成总刚矩阵。2.非线性实现（分段力，存在系统误差）假定杆件一端发生了总位移超过了其自身杆长的1/1000（此参数可调）为大

2、位移，需要考虑由于其端点发生位移而导致结构刚度变化的影响，将由原来不考虑非线性的结构各节点位移算出，求出各节点相对于自身杆件发生的位移，取其中最大值，将之设定为1/1000，由此可得出需要考虑非线性的临界荷载，将原荷载分成n份（n取整）临界荷载，逐级加载在结构上，每一次加载都考虑结构节点位移对结构刚度的影响。由此计算出结构的最终位移。3.牛顿拉普森迭代法实现非线性利用牛顿拉普森迭代法，在原来线性的基础上，计算出在加上变形的基础上求出的结构新的刚度矩阵K1，由F=F-K1*p得到差量F，再将结构在F作用下的位移求出，并求出发生该位移后结构新的刚度，重复以上步骤，直到前后两次节点变形最大值小于1e

3、-8，循环结束，输出结果。二、命令窗口输入命令（不考虑非线性）：n=input(请输入节点个数:);E=input(请输入E的大小:);A=input(请输入A的大小:);N=dianzuobiao(n); %3*n的矩阵，用于存储n个节点的坐标M=jiedianlianjie(n); %n*n的矩阵，用于存储n个节点的连接关系，连接则对应元素为1，否则为0L=l(M,N,n); %n*n矩阵，用于存储I,j节点连接所得的杆件的长度Y=ganjiangeshu(L,n); %数量值，表示杆件的个数t=T(L,N,n); %3n*3n矩阵，用于存储以上各杆件的坐标转换矩阵S=dangangjih

4、e(E,A,n,Y,L,t); %6*6y矩阵，用于存储y根杆件的各自的单刚矩阵K=zonggang(E,A,L,t,n); %3n*3n矩阵，描述总刚P=zhiling(K); %总刚矩阵主元根据约束条件置零Q=jia1(P,n); %总刚矩阵置零主元加1z=lixiangliang(n); %由输入力组成的列向量p=inv(Q)*z; %求解位移F=neili(n,Y,S,L,p); %由位移及单刚矩阵求得各杆件内力三、命令执行输入参数及计算结果：（5节点九单元，实际模型见下图）请输入节点个数:5请输入E的大小:2e8请输入A的大小:0.01从左往右从下往上编号输入第i节点x坐标：0从左往

5、右从下往上编号输入第i节点y坐标：0从左往右从下往上编号输入第i节点z坐标：0从左往右从下往上编号输入第i节点x坐标：0从左往右从下往上编号输入第i节点y坐标：0从左往右从下往上编号输入第i节点z坐标：1从左往右从下往上编号输入第i节点x坐标：0从左往右从下往上编号输入第i节点y坐标：1从左往右从下往上编号输入第i节点z坐标：0从左往右从下往上编号输入第i节点x坐标：1从左往右从下往上编号输入第i节点y坐标：0从左往右从下往上编号输入第i节点z坐标：0从左往右从下往上编号输入第i节点x坐标：1从左往右从下往上编号输入第i节点y坐标：1从左往右从下往上编号输入第i节点z坐标：0逐一输入第i个节点

6、和其他节点的连接情况,若连接,输入1,否则输入0,令和自身连接情况为0：0逐一输入第i个节点和其他节点的连接情况,若连接,输入1,否则输入0,令和自身连接情况为0：1逐一输入第i个节点和其他节点的连接情况,若连接,输入1,否则输入0,令和自身连接情况为0：1逐一输入第i个节点和其他节点的连接情况,若连接,输入1,否则输入0,令和自身连接情况为0：1逐一输入第i个节点和其他节点的连接情况,若连接,输入1,否则输入0,令和自身连接情况为0：0逐一输入第i个节点和其他节点的连接情况,若连接,输入1,否则输入0,令和自身连接情况为0：1逐一输入第i个节点和其他节点的连接情况,若连接,输入1,否则输入0

7、,令和自身连接情况为0：0逐一输入第i个节点和其他节点的连接情况,若连接,输入1,否则输入0,令和自身连接情况为0：1逐一输入第i个节点和其他节点的连接情况,若连接,输入1,否则输入0,令和自身连接情况为0：1逐一输入第i个节点和其他节点的连接情况,若连接,输入1,否则输入0,令和自身连接情况为0：1逐一输入第i个节点和其他节点的连接情况,若连接,输入1,否则输入0,令和自身连接情况为0：1逐一输入第i个节点和其他节点的连接情况,若连接,输入1,否则输入0,令和自身连接情况为0：1逐一输入第i个节点和其他节点的连接情况,若连接,输入1,否则输入0,令和自身连接情况为0：0逐一输入第i个节点和其

8、他节点的连接情况,若连接,输入1,否则输入0,令和自身连接情况为0：1逐一输入第i个节点和其他节点的连接情况,若连接,输入1,否则输入0,令和自身连接情况为0：1逐一输入第i个节点和其他节点的连接情况,若连接,输入1,否则输入0,令和自身连接情况为0：1逐一输入第i个节点和其他节点的连接情况,若连接,输入1,否则输入0,令和自身连接情况为0：1逐一输入第i个节点和其他节点的连接情况,若连接,输入1,否则输入0,令和自身连接情况为0：1逐一输入第i个节点和其他节点的连接情况,若连接,输入1,否则输入0,令和自身连接情况为0：0逐一输入第i个节点和其他节点的连接情况,若连接,输入1,否则输入0,令

9、和自身连接情况为0：1逐一输入第i个节点和其他节点的连接情况,若连接,输入1,否则输入0,令和自身连接情况为0：0逐一输入第i个节点和其他节点的连接情况,若连接,输入1,否则输入0,令和自身连接情况为0：1逐一输入第i个节点和其他节点的连接情况,若连接,输入1,否则输入0,令和自身连接情况为0：1逐一输入第i个节点和其他节点的连接情况,若连接,输入1,否则输入0,令和自身连接情况为0：1逐一输入第i个节点和其他节点的连接情况,若连接,输入1,否则输入0,令和自身连接情况为0：0输入结构约束总数：11按照编号逐一输入结构约束编号（先输入节点号）：1输入约束节点对应的位移号（输入范围：13）：1按

10、照编号逐一输入结构约束编号（先输入节点号）：1输入约束节点对应的位移号（输入范围：13）：2按照编号逐一输入结构约束编号（先输入节点号）：1输入约束节点对应的位移号（输入范围：13）：3按照编号逐一输入结构约束编号（先输入节点号）：3输入约束节点对应的位移号（输入范围：13）：1按照编号逐一输入结构约束编号（先输入节点号）：3输入约束节点对应的位移号（输入范围：13）：2按照编号逐一输入结构约束编号（先输入节点号）：3输入约束节点对应的位移号（输入范围：13）：3按照编号逐一输入结构约束编号（先输入节点号）：4输入约束节点对应的位移号（输入范围：13）：1按照编号逐一输入结构约束编号（先输入节

11、点号）：4输入约束节点对应的位移号（输入范围：13）：2按照编号逐一输入结构约束编号（先输入节点号）：5输入约束节点对应的位移号（输入范围：13）：1按照编号逐一输入结构约束编号（先输入节点号）：5输入约束节点对应的位移号（输入范围：13）：2按照编号逐一输入结构约束编号（先输入节点号）：5输入约束节点对应的位移号（输入范围：13）：3输入施加在结构上的总的力的个数：3按照编号逐一输入力作用节点编号（先输入节点号）：2逐一输入力对应的向量中的坐标值：0逐一输入力对应的向量中的坐标值：0逐一输入力对应的向量中的坐标值：-1输入力的大小：1000按照编号逐一输入力作用节点编号（先输入节点号）：2逐

12、一输入力对应的向量中的坐标值：0逐一输入力对应的向量中的坐标值：1逐一输入力对应的向量中的坐标值：0输入力的大小：1230按照编号逐一输入力作用节点编号（先输入节点号）：2逐一输入力对应的向量中的坐标值：1逐一输入力对应的向量中的坐标值：0逐一输入力对应的向量中的坐标值：0输入力的大小：1000四、计算结果：p = (m)0.000000000000000 0.000000000000000 0.000000000000000 0.002272807092008 0.000440269119346 0.000115000000000 0.000000000000000 0.0000000000

13、00000 0.000000000000000 0.000000000000000 0.000000000000000 -0.002157807092008 0.000000000000000 0.000000000000000 0.000000000000000 F = (KN)00006.82E-13100000000023001000000-23000-230-6.82E-13-10000000000-6.82E-13-1000000000-2300-1000000230002306.82E-131000000五、有限元软件验算模型建立： 模型计算结果：软件内力输出：（此处内力是x,y,

14、z方向的合力）1节点荷载0.0000000.0000002节点荷载0.0000000.0000003节点荷载0.0000000.0000004节点荷载-325.269119-325.2691195节点荷载230.000000230.0000006节点荷载0.0000000.0000007节点荷载0.0000000.0000008节点荷载0.0000000.0000009节点荷载-1732.050808-1732.050808软件位移输出：1节点荷载0.0000000.0000000.0000000.0000000.0000000.0000002节点荷载0.0022730.0004400.000

15、1150.0000000.0000000.0000003节点荷载0.0000000.0000000.0000000.0000000.0000000.0000004节点荷载0.0000000.000000-0.002158 0.000000 0.000000 0.0000005节点荷载0.0000000.0000000.0000000.0000000.0000000.000000完全一致。六、窗口命令输入（分段力考虑非线性，在上文窗口命令的基础上输入）a=maxbianxing(L,p,n,Y); %找出所有杆件中（位移/杆长）最大值f,nmax=fenduanli(z,a); %求出单元加载力

16、向量以及需要加载的次数（nmax）d,G,LL=feixianxing(E,A,f,Q,N,n,M,P,nmax); %求出考虑非线性下的位移d，并输出最终变形后的坐标值，杆长tt=T(LL,G,n); %求出最终变形后的各杆件的坐标转换矩阵KK=dangangjihe(E,A,n,Y,LL,tt); %求出最终变形后的新的当刚矩阵集合FF=neili(n,Y,KK,LL,d); %求出最终变形后的杆件的内力FF七、计算结果：nmax = 3d = (m)0.000000000000000 0.000000000000000 0.000000000000000 0.00227397482896

17、3 0.000438694002248 0.000111640092579 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 -0.002155294821124 0.000000000000000 0.000000000000000 0.000000000000000 FF = (KN)-0.72010-12.391-0.68606-4.91152997.9415-4.3808900-0.1390600226.17130.0028841000.394.380888

18、0-12.391-237.56800.035619-226.3384.941355-1001.130.0125930-0.035620.720097012.390980.6860574.91152-997.9414.380888000.13906300-226.171-0.00288-1000.39-4.38089012.39098237.56780-0.03562226.3377-4.941361001.126-0.0125900.035619八、窗口命令输入（牛顿拉普森迭代法考虑非线性，在不考虑非线性的窗口命令基础上输入）GG,LLL,s,m=newlp(E,A,N,n,M,P,p,z);

19、 %牛顿拉普森迭代法，输出最终变形坐标，杆长，位移，迭代次数tt=T(LLL,GG,n); %最终坐标下的转换矩阵KKK=dangangjihe(E,A,n,Y,LLL,tt); %最终坐标下的单刚集合FFF=neili(n,Y,KKK,LLL,s); %最终坐标下的各杠内力九、计算结果：m = 3S=(m)0.000000000000000 0.000000000000000 0.000000000000000 0.002274516177708 0.000437904574547 0.000109951291515 0.000000000000000 0.000000000000000 0

20、.000000000000000 0.000000000000000 0.000000000000000 -0.002154025184182 0.000000000000000 0.000000000000000 0.000000000000000 FFF=(KN)-0.524870-9.26307-0.519450.04666467999.9479-3.27500-0.1004400228.1069-2.04E-051001.7933.2750010-9.26307-230.61100.019917-228.23-0.0468766-1002.330.0070420-0.019920.52

21、486809.2630690.519449-0.0466647-999.9483.275001000.10043600-228.1072.04E-05-1001.79-3.27509.263069230.61150-0.01992228.23020.046876611002.334-0.0070400.019917十、附件（函数文件）（1）%输入节点坐标function M=dianzuobiao(n)for ii=1:n M(1,ii)=input(从左往右从下往上编号输入第i节点x坐标：);%输入ii节点的x轴坐标 M(2,ii)=input(从左往右从下往上编号输入第i节点y坐标：);%

22、输入ii节点的y轴坐标 M(3,ii)=input(从左往右从下往上编号输入第i节点z坐标：);%输入ii节点的z轴坐标endend（2）%形成节点连接function M=jiedianlianjie(n)for ii=1:n for jj=1:n M(ii,jj)=input(逐一输入第i个节点和其他节点的连接情况,若连接,输入1,否则输入0,令和自身连接情况为0：);%逐一输入第i个节点和其他节点的连接情况,若连接,输入1,否则输入0,令和自身连接情况为0：endend（3）%根据坐标求杆长function l=l(M,N,n) l=zeros(n,n);for ii=1:n for j

23、j=1:n if(M(ii,jj)=1) l(ii,jj)=sqrt(N(1,ii)-N(1,jj)2+(N(2,ii)-N(2,jj)2+(N(3,ii)-N(3,jj)2); %根据ij节点的坐标求出杆件的长度 else l(ii,jj)=0; end endendend（4）function y=ganjiangeshu(L,n)z=0;for ii=1:n; for jj=1:n; if(L(ii,jj)=0) z=z+1; %根据杆长矩阵元素，求得结构中共有的杆件个数 end endendy=z/2;end（5）%转换矩阵T分块function T=T(L,N,n) T=zeros(

24、3*n,3*n);for ii=1:n for jj=1:n if (L(ii,jj)=0&N(1,jj)=N(1,ii)|N(3,jj)=N(3,ii) %非竖直杆件的坐标转换矩阵 t11=(N(1,jj)-N(1,ii)/L(ii,jj); t12=(N(2,jj)-N(2,ii)/L(ii,jj); t13=(N(3,jj)-N(3,ii)/L(ii,jj); T(3*ii-2:3*ii,3*jj-2:3*jj)=t11,t12,t13;-t11*t12/sqrt(t112+t132),sqrt(t112+t132),-t12*t13/sqrt(t112+t132);-t13/sqrt(

25、t112+t132),0,t11/sqrt(t112+t132); end if (L(ii,jj)=0&N(1,jj)=N(1,ii)&N(3,jj)=N(3,ii) %竖直杆件的坐标转换矩阵 t22=(N(2,jj)-N(2,ii)/L(ii,jj); T(3*ii-2:3*ii,3*jj-2:3*jj)=0,t22,0;-t22,0,0;0,0,1; end if L(ii,jj)=0 %不存在的杆件的坐标转换矩阵，置为eye矩阵 T(3*ii-2:3*ii,3*jj-2:3*jj)=eye(3,3); end endendend（6）%按照杆件编号写出各杆件的单刚矩阵，并集合在一矩阵中

26、function K=dangangjihe(E,A,n,y,L,T)z=0;M=1 0 0 -1 0 0;0 0 0 0 0 0;0 0 0 0 0 0;-1 0 0 1 0 0;0 0 0 0 0 0;0 0 0 0 0 0;K=zeros(6,6*y);Q=zeros(6,6);for tt=1:n; for pp=tt:n if(L(tt,pp)=0) Q(1:3,1:3)=T(3*tt-2:3*tt,3*pp-2:3*pp); Q(4:6,4:6)=T(3*tt-2:3*tt,3*pp-2:3*pp); Q(1:3,4:6)=0; Q(4:6,1:3)=0; K(:,6*z+1:6*

27、z+6)=Q*E*A/L(tt,pp)*M*Q; %把每根杆的单刚矩阵集成到K矩阵中，这个矩阵区别于总刚矩阵 z=z+1; end endendend（7）%集成总刚矩阵function K=zonggang(E,A,l,T,n)K=zeros(3*n,3*n);for ii=1:n for jj=ii:n if l(ii,jj)=0 K(3*ii-2:3*ii,3*ii-2:3*ii)=K(3*ii-2:3*ii,3*ii-2:3*ii)+T(3*ii-2:3*ii,3*jj-2:3*jj)*E*A/l(ii,jj)*1,0,0;0,0,0;0,0,0*T(3*ii-2:3*ii,3*jj-

28、2:3*jj);%加上ij杆件的单刚K11对总刚的贡献 K(3*ii-2:3*ii,3*jj-2:3*jj)=K(3*ii-2:3*ii,3*jj-2:3*jj)+T(3*ii-2:3*ii,3*jj-2:3*jj)*E*A/l(ii,jj)*-1,0,0;0,0,0;0,0,0*T(3*ii-2:3*ii,3*jj-2:3*jj);%加上ij杆件的单刚K12对总刚的贡献 K(3*jj-2:3*jj,3*ii-2:3*ii)=K(3*jj-2:3*jj,3*ii-2:3*ii)+T(3*ii-2:3*ii,3*jj-2:3*jj)*E*A/l(ii,jj)*-1,0,0;0,0,0;0,0,0

29、*T(3*ii-2:3*ii,3*jj-2:3*jj);%加上ij杆件的单刚K21对总刚的贡献 K(3*jj-2:3*jj,3*jj-2:3*jj)=K(3*jj-2:3*jj,3*jj-2:3*jj)+T(3*ii-2:3*ii,3*jj-2:3*jj)*E*A/l(ii,jj)*1,0,0;0,0,0;0,0,0*T(3*ii-2:3*ii,3*jj-2:3*jj);%加上ij杆件的单刚K22对总刚的贡献 end endendend（8）%矩阵置零，主元置1%function K=zhiling(M)r=input(输入结构约束总数：);for ii=1:rjj=input(按照编号逐一输

30、入结构约束编号（先输入节点号）：);mm= input(输入约束节点对应的位移号（输入范围：13）：);M(3*jj-3+mm,:)=0;M(:,3*jj-3+mm)=0;endK=M;End（9）%主元加1%function K=jia1(M,n)for kk=1:3*nif M(kk,kk)=0 M(kk,kk)= M(kk,kk)+1; %主元加1endendK=M;End（10）%力矩阵%function z=lixiangliang(n)r=input(输入施加在结构上的总的力的个数：);z=zeros(3*n,1);for ii=1:rjj=input(按照编号逐一输入力作用节点编

31、号（先输入节点号）：); for kk=1:3 mm(kk,1)= input(逐一输入力对应的向量中的坐标值：); endf=input(输入力的大小：);z(3*jj-2:3*jj,1)=z(3*jj-2:3*jj,1)+f*mm(1,1)/sqrt(mm(1,1)2+(mm(2,1)2+(mm(3,1)2);mm(2,1)/sqrt(mm(1,1)2+(mm(2,1)2+(mm(3,1)2);mm(3,1)/sqrt(mm(1,1)2+(mm(2,1)2+(mm(3,1)2);%根据力向量将力分解成三个方向，分别加到对应的荷载值中endend（11）function M=neili(n,

32、y,Q,L,p)M=zeros(6,y);z=1;for ii=1:n for jj=ii:n if(L(ii,jj)=0) c=zeros(6,1); c(1:3,1)=p(3*ii-2:3*ii,1);c(4:6,1)=p(3*jj-2:3*jj,1); M(1:6,z)=Q(:,6*z-5:6*z)*c; %单刚与对应位移相乘，得到对应杆件的内力 z=z+1; end endendend（12）function a=maxbianxing(L,p,n,y)M=zeros(2*y,1);z=0;for ii=1:n for jj=ii:n if(L(ii,jj)=0) M(z+1,1)=s

33、qrt(p(3*ii-2,1)2+(p(3*ii-1,1)2+(p(3*ii,1)2)/L(ii,jj); %求出i节点的总位移 M(z+2,1)=sqrt(p(3*jj-2,1)2+(p(3*jj-1,1)2+(p(3*jj,1)2)/L(ii,jj); %求出j节点的总位移 z=z+2; end endenda=max(M(:); %求出所有节点的总位移最大值end（13）%分段力大小以及所分整数段function f,nmax=fenduanli(z,a)b=1/1000; %设定大位移界限值，该参数可调整nmax=ceil(a/b); %求得需要循环加载的次数，取整f=z/nmax;

34、%求得每次加载的力的大小end（14）function d,G,LL=feixianxing(E,A,f,Q,N,n,M,P,ttg)p=inv(Q)*f;d=zeros(3*n,1);pp,tt=size(p);dd=pp/3;B=zeros(3,dd);for jj=1:ttg for ii=1:dd B(1:3,ii)=p(3*ii-2:3*ii,1); end N=N+B; %用原来的坐标加上一级荷载施加在结构上产生的位移，得到新的节点坐标 L=l(M,N,n); %求出新的坐标下各杆件的长度 t=T(L,N,n); %求出新的坐标下各杆件的转换矩阵 K=zonggang(E,A,L,

35、t,n); %新的坐标下的总刚矩阵 for gg=1:3*n if(P(gg,gg)=0) K(gg,:)=0;K(:,gg)=0; %对新的总刚矩阵进行置零 end end Q=jia1(K,n); %对新的总刚矩阵加1 p=inv(Q)*f; d=d+p; %在原来的各点的位移值上加上新的变形值endfor ii=1:dd B(1:3,ii)=p(3*ii-2:3*ii,1);endN=N+B;G=N;LL=L;End（15）%牛顿拉普森迭代法求解非线性function GG,LLL,s,m=newlp(E,A,N,n,M,P,p,z)d=zeros(3*n,1);m=0;pp,tt=si

36、ze(p);dd=pp/3;B=zeros(3,dd); for ii=1:dd B(1:3,ii)=p(3*ii-2:3*ii,1); end N=N+B; %用原来的坐标加上一级荷载施加在结构上产生的位移，得到新的节点坐标 L=l(M,N,n); %求出新的坐标下各杆件的长度 t=T(L,N,n); %求出新的坐标下各杆件的转换矩阵 K=zonggang(E,A,L,t,n); %新的坐标下的总刚矩阵 for gg=1:3*n if(P(gg,gg)=0) K(gg,:)=0;K(:,gg)=0; %对新的总刚矩阵进行置零 end end Q=jia1(K,n); %对新的总刚矩阵加1 d

37、eltF=z-Q*p; d=inv(Q)*deltF; while (norm(d,inf)1e-8) for ii=1:dd B(1:3,ii)=d(3*ii-2:3*ii,1); end N=N+B; %用原来的坐标加上一级荷载施加在结构上产生的位移，得到新的节点坐标 L=l(M,N,n); %求出新的坐标下各杆件的长度 t=T(L,N,n); %求出新的坐标下各杆件的转换矩阵 K=zonggang(E,A,L,t,n); %新的坐标下的总刚矩阵 for gg=1:3*n if(P(gg,gg)=0) K(gg,:)=0;K(:,gg)=0; %对新的总刚矩阵进行置零 end end Q=jia1(K,n); %对新的总刚矩阵加1 deltF=z-Q*p; d=inv(Q)*deltF; p=p+d;m=m+1; end s=p; for ii=1:dd B(1:3,ii)=d(3*ii-2:3*ii,1);endN=N+B;GG=N;LLL=L;end