📄 k_assembly.m
字号:
function [Kgl, M] = K_assembly (in_data)
% assemble global stiffness & mass matrix for CSQ, CST, BCIZ, 3D-BEAM elements
% INPUT:
% ND - nodes matrix
% EL - elements matrix
% rhoX - density of the element: X direction
% rhoY - density of the element: Y direction
% EE - modulus of Elasticity
% hhh - element thickness
% mmiu - Poisson coefficient
% ---------------------------------------------------------------------------
% memory allocation
if in_data.EL(1,2)==0 | in_data.EL(1,2)==1 | in_data.EL(1,2)==2 % 2-D beams
dof_per_node = 3;
dof = size(in_data.ND)*dof_per_node; % # of system dof's
Kgl = zeros(dof(1),dof(1)); % global K
M = zeros(dof(1),dof(1)); % global M
[M] = M_2D_beam_assembly (in_data);
end;
if in_data.EL(1,2)==5 | in_data.EL(1,2)==4 % "4/5" - CST/CSQ element (triangle/quadrilateral)
dof = size(in_data.ND)*2; % dof - # of system dof's
Kgl = zeros(dof(1),dof(1)); % global K: dof = [x y xy]
M = zeros(dof(1),dof(1)); % global M
end;
if in_data.EL(1,2)==9 % "9" - BCIZ element (bending triangle)
dof_per_node = 3;
dof = size(in_data.ND)*dof_per_node;
Kgl = zeros(dof(1),dof(1)); % global K: dof = [z rotX rotY]
M = zeros(dof(1),dof(1)); % global M
end;
if in_data.EL(1,2)==3 % "3" - 3D-beam element
dof_per_node = 6;
dof = size(in_data.ND)*dof_per_node;
Kgl = zeros(dof(1),dof(1)); % global K: dof = [z rotX rotY]
M = zeros(dof(1),dof(1)); % global M
[M] = M_3D_beam_assembly (in_data);
end;
if in_data.EL(1,2)==6 % "6" - BRICK element (8-nodes)
dof_per_node = 3;
dof = size(in_data.ND)*dof_per_node;
Kgl = zeros(dof(1),dof(1)); % global stiffness matrix: dof = x y xy
M = zeros(dof(1),dof(1));
end;
for i=1:size(in_data.EL,1) % iterating through all elements ------------------------
if in_data.EL(i,2)==0 | in_data.EL(i,2)==1 | in_data.EL(i,2)==2 % "0/1/2" 2D-beams
node1 = find(in_data.ND(:,1)==in_data.EL(i,3));
node2 = find(in_data.ND(:,1)==in_data.EL(i,4));
if in_data.mater.E~=0 E1=in_data.mater.E; else E1=in_data.EL(i,5); end;
if in_data.mater.A~=0 A_1=in_data.mater.A; else A_1 = in_data.EL(i,6); end;
if in_data.mater.I~=0 I_1=in_data.mater.I; else I_1 = in_data.EL(i,7); end;
if in_data.EL(i,2)==0
Klc = FF_BEAM (in_data.ND(node1,2),in_data.ND(node1,3),in_data.ND(node2,2),in_data.ND(node2,3),E1,A_1,I_1);
end;
if in_data.EL(i,2)==1
Klc = FP_BEAM (in_data.ND(node1,2),in_data.ND(node1,3),in_data.ND(node2,2),in_data.ND(node2,3),E1,A_1,I_1);
end;
if in_data.EL(i,2)==2
Klc = PF_BEAM (in_data.ND(node1,2),in_data.ND(node1,3),in_data.ND(node2,2),in_data.ND(node2,3),E1,A_1,I_1);
end;
t1 = node1*3; t2 = node2*3;
% global K
Kgl((t1-2):t1,(t1-2):t1) = Kgl( (t1-2):t1,(t1-2):t1) + Klc(1:3,1:3);
Kgl((t1-2):t1,(t2-2):t2) = Kgl( (t1-2):t1,(t2-2):t2) + Klc(1:3,4:6);
Kgl((t2-2):t2,(t1-2):t1) = Kgl( (t2-2):t2,(t1-2):t1) + Klc(4:6,1:3);
Kgl((t2-2):t2,(t2-2):t2) = Kgl( (t2-2):t2,(t2-2):t2) + Klc(4:6,4:6);
end;
if in_data.EL(i,2)==4 % "4" - CST element (triangle) --------------------
node1 = find(in_data.ND(:,1)==in_data.EL(i,3));
node2 = find(in_data.ND(:,1)==in_data.EL(i,4));
node3 = find(in_data.ND(:,1)==in_data.EL(i,5));
if in_data.mater.E~=0 E1=in_data.mater.E; else E1=in_data.EL(i,6); end;
if in_data.mater.h~=0 h_1=in_data.mater.h; else h_1 = in_data.EL(i,7); end;
if in_data.mater.miu~=0 miu_1=in_data.mater.miu; else miu_1 = in_data.EL(i,8); end;
[Bsys,Esys,Klc,Msys] = D2_CST (in_data.ND(node1,2),in_data.ND(node1,3), ...
in_data.ND(node2,2),in_data.ND(node2,3), in_data.ND(node3,2),in_data.ND(node3,3),...
E1,h_1,miu_1);
Msys = Msys.*in_data.mater.rhoX;
t1 = node1*2; t2 = node2*2; t3 = node3*2;
x1 = in_data.ND(node1,2); y1 = in_data.ND(node1,3);
x2 = in_data.ND(node2,2); y2 = in_data.ND(node2,3);
x3 = in_data.ND(node3,2); y3 = in_data.ND(node3,3);
A = (0.5 * det([1 1 1; x1 x2 x3; y1 y2 y3]) ); % triangle area
% global K
Kgl((t1-1):t1,(t1-1):t1) = Kgl( (t1-1):t1,(t1-1):t1) + Klc(1:2,1:2);
Kgl((t1-1):t1,(t2-1):t2) = Kgl( (t1-1):t1,(t2-1):t2) + Klc(1:2,3:4);
Kgl((t1-1):t1,(t3-1):t3) = Kgl( (t1-1):t1,(t3-1):t3) + Klc(1:2,5:6);
Kgl((t2-1):t2,(t1-1):t1) = Kgl( (t2-1):t2,(t1-1):t1) + Klc(3:4,1:2);
Kgl((t2-1):t2,(t2-1):t2) = Kgl( (t2-1):t2,(t2-1):t2) + Klc(3:4,3:4);
Kgl((t2-1):t2,(t3-1):t3) = Kgl( (t2-1):t2,(t3-1):t3) + Klc(3:4,5:6);
Kgl((t3-1):t3,(t1-1):t1) = Kgl( (t3-1):t3,(t1-1):t1) + Klc(5:6,1:2);
Kgl((t3-1):t3,(t2-1):t2) = Kgl( (t3-1):t3,(t2-1):t2) + Klc(5:6,3:4);
Kgl((t3-1):t3,(t3-1):t3) = Kgl( (t3-1):t3,(t3-1):t3) + Klc(5:6,5:6);
% global M
M((t1-1):t1,(t1-1):t1) = M( (t1-1):t1,(t1-1):t1) + Msys(1:2,1:2);
M((t1-1):t1,(t2-1):t2) = M( (t1-1):t1,(t2-1):t2) + Msys(1:2,3:4);
M((t1-1):t1,(t3-1):t3) = M( (t1-1):t1,(t3-1):t3) + Msys(1:2,5:6);
M((t2-1):t2,(t1-1):t1) = M( (t2-1):t2,(t1-1):t1) + Msys(3:4,1:2);
M((t2-1):t2,(t2-1):t2) = M( (t2-1):t2,(t2-1):t2) + Msys(3:4,3:4);
M((t2-1):t2,(t3-1):t3) = M( (t2-1):t2,(t3-1):t3) + Msys(3:4,5:6);
M((t3-1):t3,(t1-1):t1) = M( (t3-1):t3,(t1-1):t1) + Msys(5:6,1:2);
M((t3-1):t3,(t2-1):t2) = M( (t3-1):t3,(t2-1):t2) + Msys(5:6,3:4);
M((t3-1):t3,(t3-1):t3) = M( (t3-1):t3,(t3-1):t3) + Msys(5:6,5:6);
end;
if in_data.EL(i,2)==3 % "3" - 3D-beam element --------------------
node1 = find(in_data.ND(:,1)==in_data.EL(i,3));
node2 = find(in_data.ND(:,1)==in_data.EL(i,4));
if in_data.mater.E~=0 E1=in_data.mater.E; else E1=in_data.EL(i,5); end;
if in_data.mater.G~=0 G1=in_data.mater.G; else G1 = in_data.EL(i,6); end;
if in_data.mater.b~=0 b1=in_data.mater.b; else b1 = in_data.EL(i,7); end;
if in_data.mater.ho~=0 ho1=in_data.mater.ho; else ho1 = in_data.EL(i,8); end;
Klc = D3_BEAM (in_data.ND(node1,2),in_data.ND(node1,3),in_data.ND(node1,4), ...
in_data.ND(node2,2),in_data.ND(node2,3),in_data.ND(node2,4),E1,G1,b1,ho1);
t1 = node1*6; t2 = node2*6;
% global K
Kgl((t1-5):t1,(t1-5):t1) = Kgl((t1-5):t1,(t1-5):t1) + Klc(1:6,1:6);
Kgl((t1-5):t1,(t2-5):t2) = Kgl((t1-5):t1,(t2-5):t2) + Klc(1:6,7:12);
Kgl((t2-5):t2,(t1-5):t1) = Kgl((t2-5):t2,(t1-5):t1) + Klc(7:12,1:6);
Kgl((t2-5):t2,(t2-5):t2) = Kgl((t2-5):t2,(t2-5):t2) + Klc(7:12,7:12);
end;
if in_data.EL(i,2)==6 % "6" - BRICK element (8-nodes) ------------------
% rhoX is considered = rhoY
node1 = find(in_data.ND(:,1)==in_data.EL(i,3));
node2 = find(in_data.ND(:,1)==in_data.EL(i,4));
node3 = find(in_data.ND(:,1)==in_data.EL(i,5));
node4 = find(in_data.ND(:,1)==in_data.EL(i,6));
node5 = find(in_data.ND(:,1)==in_data.EL(i,7));
node6 = find(in_data.ND(:,1)==in_data.EL(i,8));
node7 = find(in_data.ND(:,1)==in_data.EL(i,9));
node8 = find(in_data.ND(:,1)==in_data.EL(i,10));
if in_data.mater.E~=0 Em=in_data.mater.E; else Em=in_data.EL(i,11); end;
if in_data.mater.miu~=0 miu_1=in_data.mater.miu; else miu_1 = in_data.EL(i,12); end;
if in_data.mater.rho~=0 rho_=in_data.mater.rho; else rho_=in_data.EL(i,13); end; % density
[Klc,Bsys,Esys,V] = D3_LB (in_data.ND(node1,2),in_data.ND(node1,3),in_data.ND(node1,4),...
in_data.ND(node2,2),in_data.ND(node2,3),in_data.ND(node2,4),in_data.ND(node3,2),...
in_data.ND(node3,3), in_data.ND(node3,4),in_data.ND(node4,2),in_data.ND(node4,3),...
in_data.ND(node4,4),in_data.ND(node5,2),in_data.ND(node5,3),in_data.ND(node5,4),...
in_data.ND(node6,2),in_data.ND(node6,3),in_data.ND(node6,4),in_data.ND(node7,2),...
in_data.ND(node7,3),in_data.ND(node7,4),in_data.ND(node8,2),in_data.ND(node8,3),...
in_data.ND(node8,4), Em, miu_1);
% Msys = Msys.*in_data.mater.rhoX;
t1 = node1*3; t2 = node2*3; t3 = node3*3; t4 = node4*3;
t5 = node5*3; t6 = node6*3; t7 = node7*3; t8 = node8*3;
x1 = in_data.ND(node1,2); y1 = in_data.ND(node1,3); z1 = in_data.ND(node1,4);
x2 = in_data.ND(node2,2); y2 = in_data.ND(node2,3); z2 = in_data.ND(node2,4);
x3 = in_data.ND(node3,2); y3 = in_data.ND(node3,3); z3 = in_data.ND(node3,4);
x4 = in_data.ND(node4,2); y4 = in_data.ND(node4,3); z4 = in_data.ND(node4,4);
x5 = in_data.ND(node5,2); y5 = in_data.ND(node5,3); z5 = in_data.ND(node5,4);
x6 = in_data.ND(node6,2); y6 = in_data.ND(node6,3); z6 = in_data.ND(node6,4);
x7 = in_data.ND(node7,2); y7 = in_data.ND(node7,3); z7 = in_data.ND(node7,4);
x8 = in_data.ND(node8,2); y8 = in_data.ND(node8,3); z8 = in_data.ND(node8,4);
sf = ones(24,1);
fM = V*rho_/8.*sf'; % equal mass : is very approximate !!!
% global K
Kgl((t1-2):t1,(t1-2):t1) = Kgl( (t1-2):t1,(t1-2):t1) + Klc(1:3,1:3);
Kgl((t1-2):t1,(t2-2):t2) = Kgl( (t1-2):t1,(t2-2):t2) + Klc(1:3,4:6);
Kgl((t1-2):t1,(t3-2):t3) = Kgl( (t1-2):t1,(t3-2):t3) + Klc(1:3,7:9);
Kgl((t1-2):t1,(t4-2):t4) = Kgl( (t1-2):t1,(t4-2):t4) + Klc(1:3,10:12);
Kgl((t1-2):t1,(t5-2):t5) = Kgl( (t1-2):t1,(t5-2):t5) + Klc(1:3,13:15);
Kgl((t1-2):t1,(t6-2):t6) = Kgl( (t1-2):t1,(t6-2):t6) + Klc(1:3,16:18);
Kgl((t1-2):t1,(t7-2):t7) = Kgl( (t1-2):t1,(t7-2):t7) + Klc(1:3,19:21);
Kgl((t1-2):t1,(t8-2):t8) = Kgl( (t1-2):t1,(t8-2):t8) + Klc(1:3,22:24);
Kgl((t2-2):t2,(t1-2):t1) = Kgl( (t2-2):t2,(t1-2):t1) + Klc(4:6,1:3);
Kgl((t2-2):t2,(t2-2):t2) = Kgl( (t2-2):t2,(t2-2):t2) + Klc(4:6,4:6);
Kgl((t2-2):t2,(t3-2):t3) = Kgl( (t2-2):t2,(t3-2):t3) + Klc(4:6,7:9);
Kgl((t2-2):t2,(t4-2):t4) = Kgl( (t2-2):t2,(t4-2):t4) + Klc(4:6,10:12);
Kgl((t2-2):t2,(t5-2):t5) = Kgl( (t2-2):t2,(t5-2):t5) + Klc(4:6,13:15);
Kgl((t2-2):t2,(t6-2):t6) = Kgl( (t2-2):t2,(t6-2):t6) + Klc(4:6,16:18);
Kgl((t2-2):t2,(t7-2):t7) = Kgl( (t2-2):t2,(t7-2):t7) + Klc(4:6,19:21);
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -