2019.4.28

前处理程序 clear;clc XY = [1, 0, 0 2, 2, 0

3, 2, -2]; %XY为N行3列，N为节点总数 ELB = [1, 1, 2, 1

2, 2, 3, 1]; %ELB为MB行4列，MB为单元总数 b = 4.5e-2;

h = 2e-1; %截面宽b和高h A1 = b*h;

IZ1 = 3e-5; %惯性矩

EAIZ = [200e9, A1, IZ1]; %弹性模量与截面积和惯性矩 ELPQ = [1, 1, -2e4, 2

2, 1, -1e4, 2]; %ELPQ第一列为受载荷单元号，第二列为长度C（见下表），三列为载荷大小，四列为载荷类型 SU = [1,0 2, 0

3,0]; %SU第一列为已知节点位移对应的方程号，第2列为节点位移数值，约束即等于0

子程序1

function [PO, ii, jj] = dxjd(ELB, XY, ELPQ1) PO = zeros(6,1); k = ELPQ1(1);

ii = ELB(k, 2); jj = ELB(k, 3);

dxy = [XY(ii, 2), XY(ii, 3) XY(jj, 2), XY(jj, 3)]; DY = dxy(2,2) - dxy(1,2); DX = dxy(2,1) - dxy(1,1); L = sqrt(DX^2 + DY^2); C = ELPQ1(2); Q = ELPQ1(3); type = ELPQ1(4); switch type case 1

PO(2) = PO(2) + 0.5*Q*C*(2-2*C^2/L^2 + C^3/L^3); PO(3) = PO(3) + Q*C^2*(6-8*C/L+3*C^2/L^2)/12; PO(5) = PO(5)+Q*C-PO(2);

PO(6) = PO(6) - Q*C^3*(4-3*C/L)/12/L; case 2

D = L - C;

PO(2) = PO(2) + Q*(L+2*C)*D^2/L^3; PO(3) = PO(3) + Q*C*D^2/L^2; PO(5) = PO(5) + Q-PO(2);

PO(6) = PO(6) + Q*D*C^2/L^2; case 3

D = L - C;

PO(1) = PO(1) + Q*D/L; PO(4) = PO(4) + Q*C/L; case 4

PO(2) = PO(2) + 7*Q*L/20; PO(3) = PO(3) + Q*L^2/20; PO(5) = PO(5) + 3*Q*L/20; PO(6) = PO(6) + Q*L^2/30; end

子程序2

function [KE, T] = ketb(dxy, E, A, IZ) DY = dxy(2,2) - dxy(1,2); DX = dxy(2,1) - dxy(1,1); L = sqrt(DX^2 + DY^2); S = DY/L; C = DX/L; a1 = IZ/L; a2 = a1/L;

a3 = E/L;

KE = a3*[A, 0, 0, -A, 0, 0 0, 12*a2, 6*a1, 0, -12*a2, 6*a1 0, 6*a1, 4*IZ, 0, -6*a1, 2*IZ -A, 0, 0, A, 0, 0 0, -12*a2, -6*a1, 0, 12*a2, -6*a1 0, 6*a1, 2*IZ, 0, -6*a1, 4*IZ]; t = [C, S, 0; -S, C, 0; 0, 0, 1]; t1 = zeros(3,3); T = [t, t1; t1, t]; end

子程序

function KZ = kztb(XY, ELB, EAIZ) [N,M] = size(XY); KZ =zeros(9, 9); [MB, m] = size(ELB); for k = 1 :2

ii = ELB(k, 2); jj = ELB(k, 3); LTB = ELB(k, 4);

dxy = [XY(ii,2), XY(ii,3) XY(jj,2), XY(jj,3)]; E = EAIZ(LTB, 1); A = EAIZ(LTB, 2); IZ = EAIZ(LTB, 3);

[KE, T] = ketb(dxy, E, A, IZ);

CN = [3*ii-2, 3*ii-1, 3*ii, 3*jj-2, 3*jj-1, 3*jj]; KE = (T')*KE*T; for i = 1:6

for j = 1:6

KZ(CN(i), CN(j)) = KZ(CN(i), CN(j)) + KE(i, j); end end end

子程序

function U = weiyi(KZ, P, SU) [LR, m] = size(SU); for k =1:LR

i = SU(k, 1);

P(i) = 1e10*SU(k, 2); KZ(i,i) = 1e10*KZ(i, i); end

U = KZ\\P; end

子程序

function P = ydx(XY, ELB, EAIZ, ELPQ) [N, m] = size(XY) P = zeros(3*N, 1) [LPQ, m] = size(ELPQ) for j = 1:LPQ

ELPQ1 = ELPQ(j, :) k = ELPQ(j, 1)

[PO, ii, jj] = dxjd(ELB, XY, ELPQ1) LTB = ELB(k, 4);

dxy = [XY(ii,2), XY(ii,3) XY(jj,2), XY(jj, 3)]; E = EAIZ(LTB, 1); A = EAIZ(LTB, 2); IZ = EAIZ(LTB, 3);

[KE, T] = ketb(dxy, E, A, IZ ); PO = (T')*PO;

CN = [3*ii-2, 3*ii-1, 3*ii, 3*jj-2, 3*jj-1, 3*jj] for i = 1:6

P(CN(i)) = P(CN(i)) + PO(i); end end

后处理程序 clc;

[N, M] = size(XY); DU = zeros(N,4); for i = 1:N

DU(i,1) = i;

DU(i,2) = U(3*i-2); DU(i,3) = U(3*i-1); DU(i,4) = U(3*i); end

[MB, M] = size(ELB); for i = 1:MB P1(i, 1) = i; end

'节点号 水平位移 竖直位移 转角' DU

主程序 >>GJINPUT;

>>KZ = kztb(XY, ELB,EAIZ);

>>P = ydx(XY, ELB, EAIZ, ELPQ); >>U = weiyi(KZ, P, SU); >>GJOUTPUT;

运行结果： ans =

'节点号 水平位移 竖直位移 转角' DU =

1.0000 -0.0000 -0.0000 -0.0000 2.0000 -0.0000 -0.0111 -0.0100 3.0000 -0.0231 -0.0111 -0.0125