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 = 3; 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]; % %悬臂梁节点单元信息、荷载边界条件、材料参数 % 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);