弹性力学及有限元基础部分作业
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
clc clear all % 1.常量定义 E=207; A=120; L1=600; L2=L1*sqrt(3)/2*sqrt(2);
P0=-80; % 2.网格划分 n1=1; n2=1; Le1=L1/n1; Le2=L2/n2; % 3.求解Kt矩阵 Ke=zeros(4); Kt=zeros(2*(n1+n2+1)); Ke=[... 0.2500 -0.4330 -0.2500 0.4330;... -0.4330 0.7500 0.4330 -0.7500;... -0.2500 0.4330 0.2500 -0.4330;... 0.4330 -0.7500 -0.4330 0.7500 ]; Ke=Ke*E*A/Le1; for ii=1:n1 k=[ii,ii+1]; for jj=1:2 for kk=1:2 Kt(2*k(jj)-1:2*k(jj),2*k(kk)-1:2*k(kk))=... Kt(2*k(jj)-1:2*k(jj),2*k(kk)-1:2*k(kk))+... Ke(2*jj-1:2*jj,2*kk-1:2*kk); end end end Ke=[... 0.5000 0.5000 -0.5000 -0.5000;... 0.5000 0.5000 -0.5000 -0.5000;... -0.5000 -0.5000 0.5000 0.5000;... -0.5000 -0.5000 0.5000 0.5000]; Ke=Ke*E*A/Le2; for ii=n1+1:n1+n2 k=[ii,ii+1]; for jj=1:2 for kk=1:2 Kt(2*k(jj)-1:2*k(jj),2*k(kk)-1:2*k(kk))=... Kt(2*k(jj)-1:2*k(jj),2*k(kk)-1:2*k(kk))+...
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
图 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)];
ux=a(2*(n1+1)-1) uy=a(2*(n1+1)) PC=sqrt(P(1)^2+P(2)^2) PD=sqrt(P(2*(n1+n2+1)-1)^2+P(2*(n1+n2+1))^2) F_C=E*A*epsilon_element(1) F_D=E*A*epsilon_element(n1+n2) sigma_C=sigma_element(1) sigma_D=sigma_element(n1+n2)
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
图 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
求得总刚度矩阵为百度文库单位是
附录 三幅图的matlab绘图指令 图 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)
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 程序求得位移的大小以及相应的单元节点力，应力。
图 2的matlab绘图指令为
% cantilever_beam 绘图 clear all A=[0,-10]; B=[0,10];
C=[50,-10]; D=[50,10]; plot([A(1),C(1)],[A(2),C(2)],'k-','LineWidth',1.5),hold on plot([B(1),D(1)],[B(2),D(2)],'k-','LineWidth',1.5),hold on plot([D(1),C(1)],[D(2),C(2)],'k-','LineWidth',1.5),hold on plot([A(1),B(1)],[A(2)-10,B(2)+10],'k-','LineWidth',1.5),hold on axis equal E=[0,0]; F=[100,0]; plot([E(1),F(1)],[E(2),F(2)],'k-.'),hold on G=[50,-5]; H=[50,5]; K=[100,-5]; L=[100,5]; plot([G(1),K(1)],[G(2),K(2)],'k-','LineWidth',1.5),hold on plot([H(1),L(1)],[H(2),L(2)],'k-','LineWidth',1.5),hold on plot([K(1),L(1)],[K(2),L(2)],'k-','LineWidth',1.5),hold on plot([C(1),C(1)],[C(2),C(2)-10],'k-','LineWidth',1),hold on plot([K(1),K(1)],[K(2),K(2)-15],'k-','LineWidth',1),hold on plot([0,22],[-15,-15],'k-','LineWidth',1),hold on plot([28,72],[-15,-15],'k-','LineWidth',1),hold on plot([78,100],[-15,-15],'k-','LineWidth',1),hold on arrow_P=[0+4,0,0+4;-15-1,-15,-15+1]; fill(arrow_P(1,:),arrow_P(2,:),'k') arrow_P=[50-4,50,50-4;-15-1,-15,-15+1]; fill(arrow_P(1,:),arrow_P(2,:),'k') arrow_P=[50+4,50,50+4;-15-1,-15,-15+1]; fill(arrow_P(1,:),arrow_P(2,:),'k') arrow_P=[100-4,100,100-4;-15-1,-15,-15+1]; fill(arrow_P(1,:),arrow_P(2,:),'k') plot([F(1),F(1)+15],[F(2),F(2)],'k-','LineWidth',1),hold on arrow_P=[F(1)+11,F(1)+15,F(1)+11;F(2)-1,F(2),F(2)+1]; fill(arrow_P(1,:),arrow_P(2,:),'k') n=30; yinying=zeros(4,n+1); yinying(1,1:n+1)=zeros(1,n+1); yinying(2,1:n+1)=(A(2)-10):(((B(2)+10))-(A(2)-10))/n:(B(2)+10);
B点是第2个结点，其位移大小为 C，D两点分别是第1个和第3个结点，由于没有力的加载 为，。 C，D两点的约束反力分别为， ，方向为沿各自钢索的方 向。 相应的正应力为，。 10-6.如图所示梁的结点2和3处的位移和转角。已知梁的弹性模 量为E，截面惯性矩。提示，梁的刚度方程为 其中，v为y方向位移，V为剪力，M为弯矩。
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
P0=-80; % 2.网格划分 n1=1; n2=1; Le1=L1/n1; Le2=L2/n2; % 3.求解Kt矩阵 Ke=zeros(4); Kt=zeros(2*(n1+n2+1)); Ke=[... 0.2500 -0.4330 -0.2500 0.4330;... -0.4330 0.7500 0.4330 -0.7500;... -0.2500 0.4330 0.2500 -0.4330;... 0.4330 -0.7500 -0.4330 0.7500 ]; Ke=Ke*E*A/Le1; for ii=1:n1 k=[ii,ii+1]; for jj=1:2 for kk=1:2 Kt(2*k(jj)-1:2*k(jj),2*k(kk)-1:2*k(kk))=... Kt(2*k(jj)-1:2*k(jj),2*k(kk)-1:2*k(kk))+... Ke(2*jj-1:2*jj,2*kk-1:2*kk); end end end Ke=[... 0.5000 0.5000 -0.5000 -0.5000;... 0.5000 0.5000 -0.5000 -0.5000;... -0.5000 -0.5000 0.5000 0.5000;... -0.5000 -0.5000 0.5000 0.5000]; Ke=Ke*E*A/Le2; for ii=n1+1:n1+n2 k=[ii,ii+1]; for jj=1:2 for kk=1:2 Kt(2*k(jj)-1:2*k(jj),2*k(kk)-1:2*k(kk))=... Kt(2*k(jj)-1:2*k(jj),2*k(kk)-1:2*k(kk))+...
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
图 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)];
ux=a(2*(n1+1)-1) uy=a(2*(n1+1)) PC=sqrt(P(1)^2+P(2)^2) PD=sqrt(P(2*(n1+n2+1)-1)^2+P(2*(n1+n2+1))^2) F_C=E*A*epsilon_element(1) F_D=E*A*epsilon_element(n1+n2) sigma_C=sigma_element(1) sigma_D=sigma_element(n1+n2)
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
图 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
求得总刚度矩阵为百度文库单位是
附录 三幅图的matlab绘图指令 图 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)
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 程序求得位移的大小以及相应的单元节点力，应力。
图 2的matlab绘图指令为
% cantilever_beam 绘图 clear all A=[0,-10]; B=[0,10];
C=[50,-10]; D=[50,10]; plot([A(1),C(1)],[A(2),C(2)],'k-','LineWidth',1.5),hold on plot([B(1),D(1)],[B(2),D(2)],'k-','LineWidth',1.5),hold on plot([D(1),C(1)],[D(2),C(2)],'k-','LineWidth',1.5),hold on plot([A(1),B(1)],[A(2)-10,B(2)+10],'k-','LineWidth',1.5),hold on axis equal E=[0,0]; F=[100,0]; plot([E(1),F(1)],[E(2),F(2)],'k-.'),hold on G=[50,-5]; H=[50,5]; K=[100,-5]; L=[100,5]; plot([G(1),K(1)],[G(2),K(2)],'k-','LineWidth',1.5),hold on plot([H(1),L(1)],[H(2),L(2)],'k-','LineWidth',1.5),hold on plot([K(1),L(1)],[K(2),L(2)],'k-','LineWidth',1.5),hold on plot([C(1),C(1)],[C(2),C(2)-10],'k-','LineWidth',1),hold on plot([K(1),K(1)],[K(2),K(2)-15],'k-','LineWidth',1),hold on plot([0,22],[-15,-15],'k-','LineWidth',1),hold on plot([28,72],[-15,-15],'k-','LineWidth',1),hold on plot([78,100],[-15,-15],'k-','LineWidth',1),hold on arrow_P=[0+4,0,0+4;-15-1,-15,-15+1]; fill(arrow_P(1,:),arrow_P(2,:),'k') arrow_P=[50-4,50,50-4;-15-1,-15,-15+1]; fill(arrow_P(1,:),arrow_P(2,:),'k') arrow_P=[50+4,50,50+4;-15-1,-15,-15+1]; fill(arrow_P(1,:),arrow_P(2,:),'k') arrow_P=[100-4,100,100-4;-15-1,-15,-15+1]; fill(arrow_P(1,:),arrow_P(2,:),'k') plot([F(1),F(1)+15],[F(2),F(2)],'k-','LineWidth',1),hold on arrow_P=[F(1)+11,F(1)+15,F(1)+11;F(2)-1,F(2),F(2)+1]; fill(arrow_P(1,:),arrow_P(2,:),'k') n=30; yinying=zeros(4,n+1); yinying(1,1:n+1)=zeros(1,n+1); yinying(2,1:n+1)=(A(2)-10):(((B(2)+10))-(A(2)-10))/n:(B(2)+10);
B点是第2个结点，其位移大小为 C，D两点分别是第1个和第3个结点，由于没有力的加载 为，。 C，D两点的约束反力分别为， ，方向为沿各自钢索的方 向。 相应的正应力为，。 10-6.如图所示梁的结点2和3处的位移和转角。已知梁的弹性模 量为E，截面惯性矩。提示，梁的刚度方程为 其中，v为y方向位移，V为剪力，M为弯矩。
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