内容发布更新时间 : 2024/12/22 10:27:24星期一 下面是文章的全部内容请认真阅读。
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