Structural-Mechanical
/
FEM-Course-Matlab
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
FEM-Course-Matlab/16.几何非线性有限元matlab编程/几何非线性有限元-Lee_frame/main_rewrite.m

185 lines
6.8 KiB
Mathematica
Raw Permalink Normal View History

init project
2024-01-28 16:46:36 +00:00
clear all;clc;
% Ԥ<EFBFBD>ȶ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>־<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
flags_pos
init_incr = 0.01;% <EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD>ز<EFBFBD>Ԥ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĺ<EFBFBD><EFBFBD>ر<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
max_ratio = 1.00;% <EFBFBD><EFBFBD>ֹ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ر<EFBFBD>
max_steps = 1000;% <EFBFBD><EFBFBD>ֹ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ز<EFBFBD><EFBFBD><EFBFBD>
max_iter = 100; % ÿ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ز<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
trg_iter = 3; % ÿ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ز<EFBFBD><EFBFBD><EFBFBD>Ԥ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ر<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD>Ŀ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
tol = 1e-05; % ƽ<EFBFBD><EFBFBD><EFBFBD>ж<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% form = SOLVE_UL; %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>շ<EFBFBD><EFBFBD><EFBFBD>
% mtx = MTX_SR2T; % <EFBFBD><EFBFBD><EFBFBD>θնȾ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Сת<EFBFBD><EFBFBD>-<EFBFBD><EFBFBD><EFBFBD><EFBFBD>-<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% alg = ALG_ALCM_CYL; %<EFBFBD>̶<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>-Բ<EFBFBD><EFBFBD> Constant Arc-Length Control Method -Cylindrical
% incr = CONSTANT_INCREMENT; % <EFBFBD><EFBFBD><EFBFBD>ز<EFBFBD>Ԥ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĺ<EFBFBD><EFBFBD>ر<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% iter = NR_STANDARD; %<EFBFBD><EFBFBD>׼ţ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɭ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%% Plotting Options
plot_type = PLOT_BOTH;% PLOT_GRAPH -> Plot equilibrium path (load ratio vs. displacements)PLOT_MODEL -> Plot deformed structurePLOT_BOTH -> Plot equilibrium path and deformed structure
res_type = RES_STEPS;% RES_STEPS -> Plot results for every step (converged equilibrium points) RES_ITER -> Plot results for every iteration RES_BOTH -> Plot results for every step and iteration
node_id = 13;
dof = 2;%<EFBFBD><EFBFBD><EFBFBD>Ƶڶ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɶ<EFBFBD>
%% Pre-Processing
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ṹ<EFBFBD><EFBFBD>
Anl = struct('init_incr',init_incr,'max_ratio',max_ratio,'max_steps',max_steps,'max_iter',max_iter,'trg_iter',trg_iter,'tol',tol);
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ģ<EFBFBD>ͻ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ϣ<EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <EFBFBD><EFBFBD>Ԫ<EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD>Ϣ <EFBFBD>߽<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% % Williams frame (1 element)
% Coords = [[0 0];%<EFBFBD>̶<EFBFBD><EFBFBD>ڵ<EFBFBD>
% [25.886 0];%<EFBFBD>̶<EFBFBD><EFBFBD>ڵ<EFBFBD>
% [12.943 0.386]];%<EFBFBD><EFBFBD><EFBFBD>ɽڵ<EFBFBD>
% Ebc = [[1 1 1];
% [1 1 1];
% [0 0 0]];
% Loads = [[0 0 0];
% [0 0 0];
% [0 -100 0]];
% Connect = [[1 3];
% [2 3]];
% E = [1.03e+07, 1.03e+07];
% A = [0.182979, 0.182979];
% I = [0.000900394, 0.000900394];
% %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ڵ㵥Ԫ<EFBFBD><EFBFBD>Ϣ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ر߽<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϲ<EFBFBD><EFBFBD><EFBFBD> plot_node=11
% Coords = [[0.0 0.0];[0.1 0.0];[0.2 0.0];[0.3 0.0];[0.4 0.0];[0.5 0.0];[0.6 0.0];[0.7 0.0];[0.8 0.0];[0.9 0.0];[1.0 0.0]];
% Ebc = [[1 1 1];[0 0 0];[0 0 0];[0 0 0];[0 0 0];[0 0 0];[0 0 0]; [0 0 0];[0 0 0];[0 0 0];[0 0 0]];
% Loads = [[0 0 0];[0 0 0];[0 0 0]; [0 0 0]; [0 0 0];[0 0 0];[0 0 0];[0 0 0]; [0 0 0];[0 0 0];[0 -1000 0]];
% Connect = [[1 2]; [2 3];[3 4];[4 5];[5 6];[6 7];[7 8];[8 9]; [9 10]; [10 11];];
% E(1:10) = 1e+07;
% A(1:10) = 1e-02;
% I(1:10) = 1e-05;
% Lee Frame<EFBFBD>ڵ㵥Ԫ<EFBFBD><EFBFBD>Ϣ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ر߽<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϲ<EFBFBD><EFBFBD><EFBFBD> plot_node=13
Coords = [[0.00 0.00];[0.00 0.12];[0.00 0.24];[0.00 0.36];[0.00 0.48];[0.00 0.60];[0.00 0.72];
[0.00 0.84]; [0.00 0.96];[0.00 1.08];[0.00 1.20];[0.12 1.20];[0.24 1.20];[0.36 1.20];
[0.48 1.20];[0.60 1.20]; [0.72 1.20];[0.84 1.20];[0.96 1.20];[1.08 1.20];[1.20 1.20]];
Ebc = [[1 1 0];[0 0 0];[0 0 0];[0 0 0];[0 0 0];[0 0 0];[0 0 0];[0 0 0];[0 0 0];
[0 0 0];[0 0 0];[0 0 0];[0 0 0];[0 0 0];[0 0 0];[0 0 0];[0 0 0]; [0 0 0]; [0 0 0];[0 0 0]; [1 1 0]];
Loads = [[0 0 0];[0 0 0]; [0 0 0];[0 0 0];[0 0 0];[0 0 0];[0 0 0];[0 0 0];
[0 0 0];[0 0 0];[0 0 0]; [0 0 0];[0 -2 0];[0 0 0]; [0 0 0];[0 0 0];
[0 0 0]; [0 0 0];[0 0 0];[0 0 0];[0 0 0]];
Connect = [[1 2];[2 3]; [3 4];[4 5];[5 6];[6 7]; [7 8];[8 9]; [9 10]; [10 11];[11 12];
[12 13];[13 14];[14 15]; [15 16]; [16 17];[17 18]; [18 19];[19 20];[20 21]];
E(1:20) = 7200000;
A(1:20) = 0.0006;
I(1:20) = 0.00000002;
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ģ<EFBFBD>ͽṹ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
nnp = size(Coords,1);% nnp -> <EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
nel = size(Connect,1);% nel -> <EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
neq = 3 * nnp;% neq -> <EFBFBD><EFBFBD><EFBFBD≯<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɶȸ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% <EFBFBD>նȾ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>װʱ<EFBFBD>Ľڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɶ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% δ<EFBFBD><EFBFBD>Լ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɶ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ȱ<EFBFBD><EFBFBD><EFBFBD>
ID = zeros(nnp,3);
neqfixed = 0;
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѭ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ڶ<EFBFBD><EFBFBD><EFBFBD>IDû<EFBFBD><EFBFBD><EFBFBD>ã<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ϊ<EFBFBD>˵õ<EFBFBD>neqfixed<EFBFBD><EFBFBD><EFBFBD>̶<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɶȸ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
for i = 1:nnp
for j = 1:3
if (Ebc(i,j) == 1)
neqfixed = neqfixed + 1;
ID(i,j) = 1;
end
end
end
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD>i<EFBFBD><EFBFBD>j<EFBFBD><EFBFBD><EFBFBD>ɶ<EFBFBD>ID
neqfree = neq - neqfixed;
countS = neqfree;
countF = 0;
for i = 1:nnp
for j = 1:3
if ID(i,j) == 0
countF = countF + 1;
ID(i,j) = countF;
else
countS = countS + 1;
ID(i,j) = countS;
end
end
end
% <EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɶ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> (DX1, DY1, RZ1, DX2, DY2, RZ2)<EFBFBD><EFBFBD>ȫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɶ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
GLE = zeros(nel, 6);
for i = 1:nel
for j = 1:6
if (j <= 3)
GLE(i,j) = ID(Connect(i,1),j);
else
GLE(i,j) = ID(Connect(i,2),j-3);
end
end
end
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>װ<EFBFBD><EFBFBD>
P = zeros(neq,1);
for i = 1:nnp
for j = 1:3
P(ID(i,j)) = P(ID(i,j)) + Loads(i,j);
end
end
Pref=P;
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>model<EFBFBD>ṹ<EFBFBD><EFBFBD>
Model = struct('nnp',nnp,'nel',nel,'neq',neq,'neqf',neqfree,'neqc',neqfixed);
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>node <EFBFBD>ṹ<EFBFBD><EFBFBD>
Node(Model.nnp,1) = struct('x',[],'y',[],'px',[],'py',[],'mz',[],'dof',[]);
for i = 1:Model.nnp
Node(i).x = Coords(i,1);
Node(i).y = Coords(i,2);
Node(i).px = Loads(i,1);
Node(i).py = Loads(i,2);
Node(i).mz = Loads(i,3);
Node(i).dof = ID(i,:);
end
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>elem<EFBFBD>ṹ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
Elem(Model.nel,1) = struct('n1',[],'n2',[],'E',[],'A',[],'I',[],...
'L',[],'L_1',[],'fi',[],'fi_1',[],'fn',[],'fn_1',[],'angle',[],'angle_1',[],...
'gle',[],'rot',[],'ke',[],'kt',[]);
for i = 1:Model.nel
% <EFBFBD><EFBFBD>Ԫ<EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD>
n1 = Node(Connect(i,1));
n2 = Node(Connect(i,2));
% <EFBFBD>ֲ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>x<EFBFBD>ظ˼<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ֱ<EFBFBD>˼<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>z<EFBFBD><EFBFBD>ֱ<EFBFBD>ڶ<EFBFBD>άƽ<EFBFBD>淽<EFBFBD><EFBFBD>
vl = [n2.x-n1.x n2.y-n1.y 0];
vx = vl/norm(vl);
vz = [0 0 1];
vy = cross(vz, vx);
% Assemble rotation transformation matrix for one node
% T = [vx ;vy;vz];
% Set element properties
Elem(i).n1 = n1;
Elem(i).n2 = n2;
Elem(i).E = E(i);
Elem(i).A = A(i);
Elem(i).I = I(i);
Elem(i).L = norm(vl);
Elem(i).L_1 = norm(vl);%%%Length from beggining of step
Elem(i).fi = zeros(6,1); %Vector of internal forces in local system [fx1,fy1,mz1,fx2,fy2,mz2]
Elem(i).fi_1 = zeros(6,1);%%%Vector of internal forces in local system from beggining of step
Elem(i).fn = zeros(3,1);%Vector of internal forces in natural system [N2,M1,M2]
Elem(i).fn_1 = zeros(3,1);%%%Vector of internal forces in natural system from beggining of step
Elem(i).angle = atan2(n2.y-n1.y,n2.x-n1.x);%Angle with horizontal axis
Elem(i).angle_1 = atan2(n2.y-n1.y,n2.x-n1.x);%%%%Angle with horizontal axis from beggining of step
Elem(i).gle = GLE(i,:);
% Elem(i).rot = blkdiag(T,T);%Rotation matrix from global to local coordinate system
end
%% Solve System
tic
Result = solve(Model,Anl,Elem,Pref);
fprintf(1,'Time:%.6f\n',toc);
%% Post-processing
pos_processing(Model,Node,Elem,Result,plot_type,res_type,node_id,dof);