弹性力学课程设计

弹性力学课程设计

已知条件：

E=10*10^4 MPa A=0.5 m2 I=1/24 m4

主程序

节点力

F1=[0 -42 -24 0 -18 16 ];

F2=[0 -22.5 -11.25 0 -22.5 11.25];

F4=[0 -40 -30 0 -40 30];

T=[0 1 0 0 0 0

-1 0 0 0 0 0

0 0 1 0 0 0

0 0 0 0 1 0

0 0 0 -1 0 0

0 0 0 0 0 1];

K=zeros(9,9);

E=10*10^4;

A=0.5;

I=1/24;

未转换的单元刚度矩阵

K11=gdjz(E,A,I,4);

K22=gdjz(E,A,I,3);

K33=gdjz(E,A,I,4);

K44=gdjz(E,A,I,3);

K55=gdjz(E,A,I,4);

整体刚度矩阵的集成

K1=PlaneFrameElementStiffness(E,A,I,4,90); p=[0 0 0 1 2 3];

c=find(p);

cc=p(c);

K(cc,cc)=K(cc,cc)+K1(c,c);

K2=PlaneFrameElementStiffness(E,A,I,3,0); p=[1 2 3 4 5 6];

c=find(p);

cc=p(c);

K(cc,cc)=K(cc,cc)+K2(c,c);

K3=PlaneFrameElementStiffness(E,A,I,4,90);

p=[0 0 0 4 5 6];

c=find(p);

cc=p(c);

K(cc,cc)=K(cc,cc)+K3(c,c);

K4=PlaneFrameElementStiffness(E,A,I,3,0); p=[4 5 6 7 8 9];

c=find(p);

cc=p(c);

K(cc,cc)=K(cc,cc)+K4(c,c);

K5=PlaneFrameElementStiffness(E,A,I,4,90); p=[0 0 0 7 8 9 ];

c=find(p);

cc=p(c);

K(cc,cc)=K(cc,cc)+K5(c,c);

整体坐标系下的荷载

F=ans;

U=K\F'

U1=[0 0 0 U(1) U(2) U(3)];

U2=[U(1) U(2) U(3) U(4) U(5) U(6)];

U4=[ U(4) U(5) U(6) U(7) U(8) U(9)];

U3=[0 0 0 U(4) U(5) U(6)];

U5=[0 0 0 U(7) U(8) U(9)];

修正值

f1=-F1';

f2=-F2';

f4=-F4';

杆件的单元反力

FF1=K11*(T*U1')+f1;

FF2=K22*U2'+f2;

FF3=K33*(T*U3');

FF4=K44*U4'+f4;

FF5=K55*(T*U5');

整体坐标系下的单元内力

FF=[T'*FF1 FF2 T'*FF3 FF4 T'*FF5]

子程序1

function y = PlaneFrameElementStiffness(E,A,I,L,theta)

%PlaneFrameElementStiffness This function returns the element % stiffness matrix for a plane frame

% element with modulus of elasticity E,

% cross-sectional area A, moment of

% inertia I, length L, and angle

% theta (in degrees).

% The size of the element stiffness

% matrix is 6 x 6.

x = theta*pi/180;

C = cos(x);

S = sin(x);

w1 = A*C*C + 12*I*S*S/(L*L);

w2 = A*S*S + 12*I*C*C/(L*L);

w3 = (A-12*I/(L*L))*C*S;

w4 = 6*I*S/L;

w5 = 6*I*C/L;

y = E/L*[w1 w3 -w4 -w1 -w3 -w4 ; w3 w2 w5 -w3 -w2 w5 ;

-w4 w5 4*I w4 -w5 2*I ; -w1 -w3 w4 w1 w3 w4 ;

-w3 -w2 -w5 w3 w2 -w5 ; -w4 w5 2*I w4 -w5 4*I];

子程序2

function y=ans

F1=[0 -42 -24 0 -18 16 ];

F2=[0 -22.5 -11.25 0 -22.5 11.25];

F4=[0 -40 -30 0 -40 30];

T=[0 1 0 0 0 0

-1 0 0 0 0 0

0 0 1 0 0 0

0 0 0 0 1 0

0 0 0 -1 0 0

0 0 0 0 0 1 ];

F11=T'*F1';

F22=F2';

F44=F4';

y=[F11(4)+F22(1) F11(5)+F22(2) F11(6)+F22(3) F22(4)+F44(1) F22(5)+F44(2) F22(6)+F44(3) F44(4) F44(5) F44(6)];

子程序3

function y =gdjz(E,A,I,L)

y=E/L*[A 0 0 -A 0 0

0 12*I/L^2 6*I/L 0 -12*I/L^2 6*I/L

0 6*I/L 4*I 0 6*I/L 2*I

-A 0 0 A 0 0

0 -12*I/L^2 -6*I/L 0 12*I/L^2 -6*I/L

0 6*I/L 2*I 0 -6*I/L 4*I];