qzb.m

四面体有限元程序,单元刚度的计算及总刚度矩征集成

axis([0 1 0 1 0 1])
xlabel('x')
ylabel('y')
zlabel('z')
hold on
plot3([0 1],[0 0],[0 0],'k',...
      [0 0],[0 1],[0 0],'k',...
      [0 0],[0 0],[0 1],'k')
view(110,15)
   
r=[0,1]; a=pi/4; b=pi/4;  %画半径
[x1,y1,z1]=sph2cart([0,a],[0,b],r);
[x2,y2,z2]=sph2cart([0,a+pi*0.05],[0,b],r);
[x3,y3,z3]=sph2cart([0,a],[0,b+pi*0.05],r);
[x4,y4,z4]=sph2cart([0,a+pi*0.05],[0,b+pi*0.05],r);
plot3(x1,y1,z1,'b',x2,y2,z2,'b',...
      x3,y3,z3,'b',x4,y4,z4,'m :')

r1=ones(1,11);      %画小体元外表面
theta1=a*ones(1,11);
phi1=b:0.005*pi:b+pi*0.05;
[x1,y1,z1]=sph2cart(theta1,phi1,r1);
theta2=(a+pi*0.05)*ones(1,11);
phi2=b:0.005*pi:b+pi*0.05;
[x2,y2,z2]=sph2cart(theta2,phi2,r1);
theta3=a:0.005*pi:a+pi*0.05;
phi3=b*ones(1,11);
[x3,y3,z3]=sph2cart(theta3,phi3,r1);
theta4=a:0.005*pi:a+pi*0.05;
phi4=(b+pi*0.05)*ones(1,11);
[x4,y4,z4]=sph2cart(theta4,phi4,r1);
plot3(x1,y1,z1,'b',x2,y2,z2,'b',...
   x3,y3,z3,'b',x4,y4,z4,'b')

r2=0.85*r1;   %画小体元内表面
[x5,y5,z5]=sph2cart(theta1,phi1,r2);
[x6,y6,z6]=sph2cart(theta2,phi2,r2);
[x7,y7,z7]=sph2cart(theta3,phi3,r2);
[x8,y8,z8]=sph2cart(theta4,phi4,r2);
plot3(x5,y5,z5,'b',x6,y6,z6,'c :',...
      x7,y7,z7,'b',x8,y8,z8,'c :')