Global numbering of nodes (italics), elements, and degrees of freedom (numbers on vectors):2543 1 22 1 3 1 6Form element stiffiness matrix from column vectors (see text p. 56):with(linalg); a1:=matrix(4,1,[-c,-s,c,s]);a2:=matrix(1,4,[-c,-s,c,s]); := a1-c-scs := a2-c -scsMultiply these to get element stiffness matrix of Eq. 2.10: k:=evalm( (A*E/L) * a1 &* a2 ); := kAEc2LAEcsL-AEc2L-AEcsLAEcsLAEs2L-AEcsL-AEs2L-AEc2L-AEcsLAEc2LAEcsL-AEcsL-AEs2LAEcsLAEs2LTrigonometric relations c:=cos(theta);s:=sin(theta); theta:=arctan((y[2]-y[1])/(x[2]-x[1])); := c ()cos θ := s ()sin θ := θarctan−y2y1−x2x1Nodal coordinates of element 1 x[1]:=0;y[1]:=0;x[2]:=1.5;y[2]:=.25;Set precision, get length Digits:=4;L:=sqrt((x[2]-x[1])^2 + (y[2]-y[1])^2 );:= L1.521Define area and modulus (unprotecting E this way is dangerous) A:=3.142e-4;unprotect(E);E:=210e9; := A.0003142 := E.210 1012Evaluate stiffness matrix, save as k1 k1:=map(eval,k); := k1.4221 108.7035 107-.4221 108-.7035 107.7035 107.1173 107-.7035 107-.1173 107-.4221 108-.7035 107.4221 108.7035 107-.7035 107-.1173 107.7035 107.1173 107Redefine nodal coordinates, for element 2 x[1]:=1.5;y[1]:=.25;x[2]:=0;y[2]:=.5;Reevaluate stiffness matrix, save as k2 k2:=map(eval,k); := k2.4221 108-.7035 107-.4221 108.7035 107-.7035 107.1173 107.7035 107-.1173 107-.4221 108.7035 107.4221 108-.7035 107.7035 107-.1173 107-.7035 107.1173 107Define global stiffness matrix K:=matrix( 6,6,[ [ k1[1,1], k1[1,2], k1[1,3], k1[1,4], 0, 0 ], [ k1[2,1], k1[2,2], k1[2,3], k1[2,4], 0, 0 ], [ k1[3,1], k1[3,2], k1[3,3]+k2[1,1], k1[3,4]+k2[1,2], k2[1,3], k2[1,4] ], [ k1[4,1], k1[4,2], k1[4,3]+k2[2,1], k1[4,4]+k2[2,2], k2[2,3], k2[2,4] ], [ 0, 0, k2[3,1], k2[3,2], k2[3,3], k2[3,4] ], [ 0, 0, k2[4,1], k2[4,2], k2[4,3], k2[4,4] ] ] ); :=
View Full Document