弹性力学及有限元基础部分作业

图 2 第10-6题图 解：本题只有三个结点，两个单元，计算结果应该是准确度不 会很高。 单元1的刚度矩阵 单元1的刚度矩阵
总刚度矩阵 载荷列阵为 为了消除总刚度矩阵的奇异性，利用乘大数法引入位移边 界条件，总刚度矩阵可以改写为 因此，位移列阵 于是，结点2和3处的位移和转角分别为 ， ， 10-7.如图 3所示薄板结构用四个三角形单元离散： i结 j结 m结 单元 点 点 点 A B C D 求结构的刚度矩阵。 1 3 2 1 4 4 3 3 3 5 5 2
Ke(2*jj-1:2*jj,2*kk-1:2*kk); end end end % 4.求解P矩阵 P=zeros(1,2*(n1+n2+1)); P(2*(n1+1))=P0; % 5.消除刚度矩阵奇异性 K=Kt; K(1,1)=K(1,1)*1e9; K(2,2)=K(2,2)*1e9; K(2*(n1+n2+1)-1,2*(n1+n2+1)-1)=K(2*(n1+n2+1)-1,2* (n1+n2+1)-1)*1e9; K(2*(n1+n2+1),2*(n1+n2+1))=K(2*(n1+n2+1),2*(n1+n2+1))*1e9; % 6.求解位移 a=(K\P')'; % 7.求解节点力 Pt=(Kt*a')'; % 8.求解节点应力 sigma_element=zeros(1,n1+n2); epsilon_element=zeros(1,n1+n2); B=[-1/2 sqrt(3)/2 1/2 -sqrt(3)/2]/Le1; for ii=1:n1 k=[ii,ii+1]; a_e(1:2)=a(2*k(1)-1:2*k(1)); a_e(3:4)=a(2*k(2)-1:2*k(2)); epsilon_element(ii)=B*a_e'; sigma_element(ii)=E*B*a_e'; end B=[-1/sqrt(2) -1/sqrt(2) 1/sqrt(2) 1/sqrt(2)]/Le2; for ii=n1+1:n1+n2 k=[ii,ii+1]; a_e(1:2)=a(2*k(1)-1:2*k(1)); a_e(3:4)=a(2*k(2)-1:2*k(2)); epsilon_element(ii)=B*a_e'; sigma_element(ii)=E*B*a_e'; end
10-4.图 1中所示钢索BC长600mm,BC和BD与铅垂线的夹角分别 为，在B点处沿竖直方向作用一集中力P=80kN，求B点的水平 方向和竖直方向的位移，点C和D处的支反力，单元节点力和 应力。钢索的弹性模量E=207GPa，截面面积。
图 1 第10-4题图 设BC钢索上某个单元内的两个结点的左边分别为，位移分 别为，那么该单元内的正应变的大小为 即 其中， 单位刚度矩阵为 同理BD钢索上某个单元内正应变的大小为 即 其中， 单位刚度矩阵为 由于BC钢索、BD钢索的内的应变均为常量，故将BC钢索 作为一个单元，BD钢索作为一个单元。可以利用如下的matlab 程序求得位移的大小以及相应的单元节点力，应力。
plot(B(1),B(2),'k.', 'Marker','.','MarkerSize',15) plot(D(1),D(2),'k.', 'Marker','.','MarkerSize',15) text(C(1)-30,C(2)-20,'C'); text(B(1)+20,B(2),'B'); text(D(1)+10,D(2)-20,'D'); text(B(1)+20,B(2)-200,'P','FontWeight','Bold'); text(C(1)-100,C(2)-100,'L_{BC}=600mm'); text(B(1)-100,B(2)+200,'\alpha=30^{。}'); text(B(1)+50,B(2)+180,'\beta=45^{。}'); plot([B(1),B(1)],[B(2),B(2)+200],'k-.'),hold on plot([B(1),B(1)],[B(2),B(2)-200],'k-','LineWidth',1.5),hold on arrow_P= [B(1),B(1)+8,B(1)-11;B(2)-200,B(2)-200+30,B(2)-200+30]; fill(arrow_P(1,:),arrow_P(2,:),'k') axis([-200,1000,-800,100]) theta=90:(120-90)/30:120; theta=theta*pi/180; r=100; x0=B(1); y0=B(2); x=x0+r*cos(theta); y=y0+r*sin(theta); plot(x,y,'k-'),hold on theta=45:(120-90)/45:90; theta=theta*pi/180; r=150; x0=B(1); y0=B(2); x=x0+r*cos(theta); y=y0+r*sin(theta); plot(x,y,'k-'),hold on box off
end BB=BB/D; Ke=conj(BB')*DD*BB*A; for jj=1:3 for kk=1:3 Kt((k(jj)*2-1):k(jj)*2,(k(kk)*2-1):k(kk)*2)=... Kt((k(jj)*2-1):k(jj)*2,(k(kk)*2-1):k(kk)*2)+... Ke((jj*2-1):jj*2,(kk*2-1):kk*2); end end end
B点是第2个结点，其位移大小为 C，D两点分别是第1个和第3个结点，由于没有力的加载 为，。 C，D两点的约束反力分别为， ，方向为沿各自钢索的方 向。 相应的正应力为，。 10-6.如图所示梁的结点2和3处的位移和转角。已知梁的弹性模 量为E，截面惯性矩。提示，梁的刚度方程为 其中，v为y方向位移，V为剪力，M为弯矩。
求得总刚度矩阵为令 图 1的matlab绘图指令
C=[0,0]; B=[600/2,-600*sqrt(3)/2]; plot([C(1),B(1)],[C(2),B(2)],'k-'),hold on axis equal D=[B(1)-B(2),0]; plot([D(1),B(1)],[D(2),B(2)],'k-'),hold on plot([C(1)-50,D(1)+50],[C(2),D(2)],'k-','LineWidth',2),hold on n=20; yinying=zeros(2,n+1); yinying(1,1:n+1)=(C(1)-50):((D(1)+50)-(C(1)-50))/n:(D(1)+50); yinying(2,1:n+1)=zeros(1,n+1); yinying(3,1:n+1)=yinying(1,1:n+1)+20*ones(1,n+1); yinying(4,1:n+1)=yinying(2,1:n+1)+20*ones(1,n+1); for ii=1:n+1 plot([yinying(1,ii),yinying(3,ii)], [yinying(2,ii),yinying(4,ii)],'k-'),hold on end plot(C(1),C(2),'k-', 'Marker','.','MarkerSize',15)
图 3 第10-7题图 解：对于某一单元 A是该单元的面积。 单元刚度矩阵 总刚度矩阵，代入诸 后可以得到 其中，t为厚度，E为杨氏模量，ν是泊松比 取h=40mm，w=100mm，E=200GPa，ν=0.32，t=1mm，可以利
用如下的matlab程序求解总刚度矩阵。
clc clear all E=200; % 杨氏模量 niu=0.32; % 泊松比 DD=[1,niu,0;niu,1,0;0,0,(1-niu)/2]; DD=DD*E/(1-niu^2); h=40; w=100; t=1; coord=conj([... 0 0;... 0 h;... w/2 h/2;... w 0;... w h ]'); nc=length(coord); element=[... 1 4 3;... 3 4 5;... 2 3 5;... 1 3 2 ]'; ne=length(element); for ii=1:ne k=element(1:3,ii)'; % 计算Ke x=coord(1,k(1:3)); y=coord(2,k(1:3)); b=[y(2)-y(3),y(3)-y(1),y(1)-y(2)]; c=[x(3)-x(2),x(1)-x(3),x(2)-x(1)]; D=det([1,1,1;x;y]); A=D/2; for jj=1:3 BB(1:3,(2*jj-1):2*jj)=[... b(jj) 0;... 0 c(jj);... c(jj) b(jj)];
clc clear all % 1.常量定义 E=207; A=120; L1=600; L2=L1*sqrt(3)/2*sqrt(2);