d3_beam.m

结构力学中的有限元例子,包含了7个分类文件夹

function K = D3_BEAM (x1, y1, z1, x2, y2, z2, E, G, b, h) 

% returns a 3D fixed-fixed beam element stiffness matrix in global coordinates
% INPUT:
%       x1, y1, z1 - coordinates of joint 1 of the truss elements
%       x2, y2, z2 - coordinates of joint 2 of the truss elements
%       E      - elastic (Young's) modulus
%       G      - shear modulus
%       I      - moment of inertia
L = sqrt((x2-x1)^2 + (y2-y1)^2 + (z2-z1)^2);
A = b*h;
Ax= b*L;      % shear area normal to x
Ay= b*h;
Az= h*b;
Ix= b*L^3/12; % moment of inertia normal to x
Iy= h*b^3/12;
Iz= b*h^3/12;
Jx= Iy+Iz;
if Ix== 0 J=Jx; else J=Ix; end;
phiy = 12*E*Iz/(G*Az*L^2);
phiz = 12*E*Iy/(G*Ay*L^2);
az = 12*E*Iz/(L^3*(1+phiy));
ay = 12*E*Iy/(L^3*(1+phiz));
bz = -12*E*Iz/(L^3*(1+phiy));
by = -12*E*Iy/(L^3*(1+phiz));
cz = 6*E*Iz/(L^2*(1+phiy));
cy = 6*E*Iy/(L^2*(1+phiz));
dz = -6*E*Iz/(L^2*(1+phiy));
dy = -6*E*Iy/(L^2*(1+phiz));
ez = (4+phiy)*E*Iz/(L*(1+phiy));
ey = (4+phiz)*E*Iy/(L*(1+phiz));
fz = (2-phiy)*E*Iz/(L*(1+phiy));
fy = (2-phiz)*E*Iy/(L*(1+phiz));


Ko= [A*E/L  0   0    0      0   0 -A*E/L  0    0    0     0    0;
     0      az  0    0      0   cz  0     bz   0    0     0   cz;
     0      0   ay   0     dy   0   0     0    by   0    dy    0;
     0      0   0  G*J/L    0   0   0     0    0  -G*J/L  0    0;
     0      0   dy   0     ey   0   0     0    cy   0    fy    0;
     0      cz  0    0      0   ez  0    dz    0    0     0   fz;
    -A*E/L  0   0    0      0   0  A*E/L  0    0    0     0    0;
     0      bz  0    0      0   dz  0    az    0    0     0   dy;
     0      0   by   0     cy   0   0     0    ay   0    cz    0;
     0      0   0  -G*J/L   0   0   0     0    0  G*J/L   0    0;
     0      0   dy   0     fy   0   0     0    cz   0    ey    0;
     0      cz  0    0      0   fz  0    dy    0    0     0   ez];
  
  
Lxy = sqrt((x2-x1)^2 + (y2-y1)^2);
d = 0.0001*L;
thau = 0;  % adjustment angle
if Lxy>d S1=(y2-y1)/Lxy; end;
if Lxy<d S1=0; end;
S2 = (z2-z1)/L;
S3 = sin(thau);
if Lxy>d C1=(x2-x1)/Lxy; end;
if Lxy<d C1=1; end;
C2 = Lxy/L;
C3 = cos(thau);

T = [C1*C2                S1*C2                S2;
   (-C1*S2*S3-S1*C3)   (-S1*S2*S3+C1*C3)    S3*C2;
   (-C1*S2*C3+S1*S3)   (-S1*S2*C3-C1*S3)    C3*C2];  % transformation matrix
T0 = zeros(3,3);
Tr = [T   T0   T0   T0;
      T0  T    T0   T0;
      T0  T0   T    T0;
      T0  T0   T0   T];
      
K = Tr'*Ko*Tr;  % K in global coordinates
% ------------------end