%% Function to solve nonlinear system of equations % % INPUT ARGUMENTS: % * Model -> Structure with model parameters: % - nnp: Number of nodes % - nel: Number of elements % - neq: Number of equation % - neqf: Number of free DOF equations % - neqc: Number of constrained DOF equations % * Anl -> Structure with analysis parameters: % - form: Nonlinear beam formulation % - alg: Solution algorithm % - incr_type: Increment size update strategy in predictor step % - iter_type: Iteration strategy % - geom_mtx: Geometric stiffness matrix % - init_incr: Initial load ratio increment % - max_ratio: Maximum load ratio % - max_steps: Maximum number of steps % - max_iter: Maximum number of iterations % - trg_iter: Target number of iterations % - tol: Tolerance for iteration convergence % * Node -> Vector of structures with nodes information: % - x: X-axis coordinate % - y: Y-axis coordinate % - px: Nodal load in global X-axis direction % - py: Nodal load in global Y-axis direction % - mz: Concentrated moment about global Z-axis % - dof: Vector with global DOF numbers % * Elem -> Vector of structures with elements information: % - n1: Handle to initial node % - n2: Handle to final node % - E: Elasticity modulus % - A: Cross-section area % - I: Cross-section moment of inertia % - L: Length % - L_1: Length from beggining of step % - fi: Vector of internal forces in local system [fx1,fy1,mz1,fx2,fy2,mz2] % - fi_1: Vector of internal forces in local system from beggining of step % - fn: Vector of internal forces in natural system [N2,M1,M2] % - fn_1: Vector of internal forces in natural system from beggining of step % - angle: Angle with horizontal axis % - angle_1: Angle with horizontal axis from beggining of step % - gle: Gather vector with global DOF numbers % - rot: Rotation matrix from global to local coordinate system % - ke: Elastic stiffness matrix % - kt: Tangent stiffness matrix % * Pref -> Reference load vector (total applied nodal load) % % OUTPUT: % * Result -> Vector of structures with results to be plotted: % - U_step: Matrix of nodal displacements (each column is the displacements vector in each step of analysis) % - U_iter: Matrix of nodal displacements (each column is the displacements vector in each iteration of analysis) % - lr_step: Vector of load ratio in all steps % - lr_iter: Vector of load ratio in all iterations % % NOTATION: % * lr -> load ratio % * U -> Vector of nodal displacements % * D_ -> Indicate accumulated step increment of a variable % * d_ -> Indicate iterative increment of a variable % function Result = solve(Model,Anl,Elem,Pref) % 初始化节点位移和荷载比 U = zeros(Model.neq,1); lr = 0; % 初始化result结构体，并插入初始平衡点 Result = struct('U_step',[],'U_iter',[],'lr_step',[],'lr_iter',[]); Result.U_step(:,1) = U; Result.U_iter(:,1) = U; Result.lr_step(1) = lr; Result.lr_iter(1) = lr; %==================================================================================================================================== % 开始增量步求解过程 step = 0; while (step < Anl.max_steps) step = step + 1; % 更新参考构型下的UL方程（L_1根据上一增量步结束时的位移U进行更新；angle_1随着迭代步更新并将上一增量步结束的angle作为该步angle_1;内力fi_1同angle） % 此处更新的L_1,angle_1和phi_1会在求解内力的函数中用到intForcesUL for i = 1:Model.nel Elem(i).L_1 = elemLength(Elem(i),U); %L_1 Length from beggining of step Elem(i).angle_1 = Elem(i).angle; %angle_1 Angle with horizontal axis from beggining of step Elem(i).fi_1 = Elem(i).fi; %fi_1 Vector of internal forces in local system from beggining of step end % 切线刚度矩阵（根据Ui-1更新Ki0) [Kt,Elem] = tangStiffMtxUL(Model,Anl,Elem,U) ; %%根据i-1增量步的位结束移向量更新节点位置后得到的单元刚度矩阵（单元长度，单元角度），但返回的Elem只更新了Ke用来接下来计算内力 % 预测解位移的切线增量delta_Ui1 d_Ut = solveLinearSystem(Model,Kt,Pref);%delta_Ui1（预测阶段的位移的切线增量，接下来是计算预测阶段的荷载比增量（增量步是否为第一增量步对应的荷载比增量公式不同） if (step == 1) s = 1; %第一增量步的预测阶段的增量的符号默认为正 d_lr = Anl.init_incr; %delta_lamda_11 人为规定预测阶段d_lr n2 = d_Ut'*d_Ut;%第一增量步中位移切线增量的范数的平方值,会在后续分析步计算GSP else % Generalized Stiffness Parameter GSP = n2 / (d_Ut'*d_Ut_1);%公式（41） d_Ut_1 上一个增量步 d_Ut 当前增量步 % 增量步符号 if (GSP < 0) s = -s; %公式（42） end % 增量步调整系数J % J = sqrt((Anl.trg_iter + 1) / (iter + 1));%目标迭代次数/上一增量步步的迭代次数（公式（40））必须将其增加 1 以避免被零除，因为预测的解对应于零次迭代 J = 1; %% 预测阶段荷载比增量 采用ALG_ALCM_CYL 固定弧长法-圆柱 % Extract free DOF components D_U_temp = D_U(1:Model.neqf); d_Ut_temp = d_Ut(1:Model.neqf); d_lr = J * sqrt((D_U_temp'*D_U_temp) /(d_Ut_temp'*d_Ut_temp));%公式38：圆柱弧长法 % d_lr = J * abs(d_lr_1);%% ALG_LCM (Load Increment) % d_lr = J * sqrt(abs((D_lr*Pref'*D_U) /(Pref'*d_Ut))); %ALG_WCM (Work Increment) % Apply correct sign d_lr = s * d_lr; end % 荷载比要小于规定值（==1） if ((Anl.max_ratio > 0.0 && lr + d_lr > Anl.max_ratio) ||... (Anl.max_ratio < 0.0 && lr + d_lr < Anl.max_ratio)) d_lr = Anl.max_ratio - lr;%保证恰好达到最终的荷载比增量1实现荷载的全部施加 end % 预测阶段位移增量 d_U = d_lr * d_Ut;%根据公式（36）预测阶段假定残差为零因此，只有前一项切线增量 % % 存储预测解 d_lr_1 = d_lr; d_Ut_1 = d_Ut;%判断GSP符号，公式（41） d_U_1 = d_U; % 前步的载荷比和位移增量（后续会随着迭代不断更新叠加）大写D代表增量步对应的位移和荷载比增量，小写d代表迭代步的位移和荷载比增量 D_lr = d_lr; D_U = d_U; % 总荷载比和位移 lr = lr + d_lr; U = U + d_U; %------------------------------------------------------------------------------------------------------------------------------- % Start iterative process iter = 0; conv = 0; while (conv == 0 && iter < Anl.max_iter) % External and internal forces P = lr * Pref; [F,Elem] = intForcesUL(Model,Elem,U,D_U,1);%采用上一次迭代计算得到的Ke,并根据D_U对Ele.angle和Ele.fi进行更新 %内力是基于位移计算，与位移存在非线性特性，因而产生残差 % Unbalanced forces R = P - F; % Store iteration results Result.U_iter(:,size(Result.U_iter,2)+1) = U; Result.lr_iter(size(Result.lr_iter,2)+1) = lr; % Check for convergence conv = (norm(R(1:Model.neqf))/norm(Pref(1:Model.neqf)) < Anl.tol); if (conv == 1) break; end % Start/keep corrective cycle of iterations iter = iter + 1; [Kt,Elem] = tangStiffMtxUL(Model,Anl,Elem,U); %每次迭代更新，标准牛顿拉普森方法，Elem.ke更新 % Tangent and residual increments of displacements d_Ut = solveLinearSystem(Model,Kt,Pref); d_Ur = solveLinearSystem(Model,Kt,R); % 荷载比修正（公式（45））恒定弧长法-圆柱法 a = d_Ut'*d_Ut; b = d_Ut'*(d_Ur + D_U); c = d_Ur'*(d_Ur + 2*D_U);%D_U为上一迭代步的值，d_Ur为本步更新后的值 s_iter = sign(D_U'*d_Ut); d_lr = -b/a + s_iter*sqrt((b/a)^2 - c/a); %为何会出现复数 if (~isreal(d_lr)) conv = -1; break; end % 公式（36） d_U = d_lr * d_Ut + d_Ur; % Increments of load ratio and displacements for current step D_lr = D_lr + d_lr; D_U = D_U + d_U; % Total values of load ratio and displacements lr = lr + d_lr; U = U + d_U; end %-------------------------------------------------------------------------------------------------------------------------------- % Check for convergence fail reached or complex value of increment if (conv == 0 && Anl.method == METH_INCR_ITER) fprintf(1,'Convergence not reached!\nStep = %d\nIter = %d\n',step,iter); return; elseif (conv == -1) fprintf(1,'Unable to compute load increment!\nStep = %d\nIter = %d\n',step,iter); return; end % Store equilibrium configuration Result.U_step(:,step+1) = U; Result.lr_step(step+1) = lr; fprintf(1,'Step:%d | ratio:%.2f | Iter:%d\n',step,lr,iter); % Check if maximum load ratio was reached if (lr >= 0.999*Anl.max_ratio) break; end end end