185 lines
6.8 KiB
Matlab
185 lines
6.8 KiB
Matlab
clear all;clc;
|
||
% 预先定义标志变量
|
||
flags_pos
|
||
|
||
init_incr = 0.01;% 初始荷载步预测解的荷载比增量
|
||
max_ratio = 1.00;% 终止计算时最大荷载比
|
||
max_steps = 1000;% 终止计算时的最大荷载步数
|
||
max_iter = 100; % 每个荷载步最大迭代次数
|
||
trg_iter = 3; % 每个荷载步的预测解求解时调整荷载比增量时的目标迭代次数
|
||
tol = 1e-05; % 平衡判断条件
|
||
|
||
% form = SOLVE_UL; %更新拉格朗日方程
|
||
% mtx = MTX_SR2T; % 几何刚度矩阵（小转动-二阶-基于轴力）
|
||
% alg = ALG_ALCM_CYL; %固定弧长法-圆柱 Constant Arc-Length Control Method -Cylindrical
|
||
% incr = CONSTANT_INCREMENT; % 荷载步预测解的荷载比增量不变
|
||
% iter = NR_STANDARD; %标准牛顿拉普森迭代法
|
||
|
||
%% 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 = 11;
|
||
dof = 2;%绘制第二自由度
|
||
%% Pre-Processing
|
||
%定义分析参数结构体
|
||
Anl = struct('init_incr',init_incr,'max_ratio',max_ratio,'max_steps',max_steps,'max_iter',max_iter,'trg_iter',trg_iter,'tol',tol);
|
||
% 定义模型基本信息：节点坐标 单元节点信息 边界条件 荷载条件 梁截面参数
|
||
% % Williams frame (1 element)
|
||
% Coords = [[0 0];%固定节点
|
||
% [25.886 0];%固定节点
|
||
% [12.943 0.386]];%自由节点
|
||
% 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];
|
||
%悬臂梁节点单元信息、荷载边界条件、材料参数 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节点单元信息、荷载边界条件、材料参数 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;
|
||
|
||
% 创建模型结构体变量
|
||
|
||
nnp = size(Coords,1);% nnp -> 节点个数
|
||
nel = size(Connect,1);% nel -> 单元个数
|
||
neq = 3 * nnp;% neq -> 方程个数（自由度个数）
|
||
|
||
% 刚度矩阵组装时的节点自由度索引矩阵
|
||
% 未被约束的自由度最先编号
|
||
ID = zeros(nnp,3);
|
||
neqfixed = 0;
|
||
%下面的循环对于定义ID没有用，纯粹是为了得到neqfixed（固定点的自由度个数）
|
||
for i = 1:nnp
|
||
for j = 1:3
|
||
if (Ebc(i,j) == 1)
|
||
neqfixed = neqfixed + 1;
|
||
ID(i,j) = 1;
|
||
end
|
||
end
|
||
end
|
||
%定义节点i的j自由度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
|
||
|
||
% 单元的自由度索引 (DX1, DY1, RZ1, DX2, DY2, RZ2)的全局自由度索引编号
|
||
|
||
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
|
||
|
||
% 荷载向量装配
|
||
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;
|
||
% 创建model结构体
|
||
Model = struct('nnp',nnp,'nel',nel,'neq',neq,'neqf',neqfree,'neqc',neqfixed);
|
||
|
||
% 创建node 结构体
|
||
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
|
||
|
||
% 创建elem结构体向量
|
||
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
|
||
% 单元节点号
|
||
n1 = Node(Connect(i,1));
|
||
n2 = Node(Connect(i,2));
|
||
|
||
% 局部坐标系，x沿杆件方向，y垂直杆件方向，z垂直于二维平面方向
|
||
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);
|