SimPC博士：几何非线性有限元基本原理及matlab编程（下）

  
（预测阶段代码）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会在求解内力的函数中用到intForcesULfor i = 1:Model.nelElem(i).L_1     = elemLength(Elem(i),U); %L_1  Length from beggining of stepElem(i).angle_1 = Elem(i).angle; %angle_1 Angle with horizontal axis from beggining of stepElem(i).fi_1    = Elem(i).fi; %fi_1 Vector of internal forces in local system from beggining of stepend% 切线刚度矩阵（根据Ui-1更新Ki0)[Kt,Elem] = tangStiffMtxUL(Model,Anl,Elem,U) ;  %%根据i-1增量步的位结束移向量更新节点位置后得到的单元刚度矩阵（单元长度，单元角度），但返回的Elem只更新了Ke用来接下来计算内力% 预测解位移的切线增量delta_Ui1d_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;%第一增量步中位移切线增量的范数的平方值,会在后续分析步计算GSPelse GSP = n2 / (d_Ut'*d_Ut_1);%公式（41）    % Generalized Stiffness Parameter   if (GSP < 0)s = -s; %公式（42） 增量步符号end        % 增量步调整系数J% J = sqrt((Anl.trg_iter + 1) / (iter + 1));%目标迭代次数/上一增量步步的迭代次数（公式（40））必须将其增加 1 以避免被零除，因为预测的解对应于零次迭代J = 1;      %% 预测阶段荷载比增量  采用圆柱弧长法 % Extract free DOF componentsD_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：圆柱弧长法% Apply correct signd_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）预测阶段假定残差为零因此，只有前一项切线增量   % Store predicted resultsd_lr_1 = d_lr; d_Ut_1 = d_Ut; d_U_1  = d_U;       % 前步的载荷比和位移增量（后续会随着迭代不断更新叠加）D_lr = d_lr;  D_U  = d_U;   % 总荷载比和位移lr = lr + d_lr;U  = U  + d_U;   %  

  
（校正迭代阶段代码）   % Start iterative processiter = 0;conv = 0;while (conv == 0 && iter < Anl.max_iter)% External and internal forcesP = lr * Pref;[F,Elem] = intForcesUL(Model,Elem,U,D_U,1);%采用上一次迭代计算得到的Ke,并根据D_U对Ele.angle和Ele.fi进行更新R = P - F;    %残余力矢量   % Store iteration resultsResult.U_iter(:,size(Result.U_iter,2)+1) = U;Result.lr_iter(size(Result.lr_iter,2)+1) = lr;% Check for  convergenceconv = (norm(R(1:Model.neqf))/norm(Pref(1:Model.neqf)) < Anl.tol);if (conv == 1)break;end% Start/keep corrective cycle of iterationsiter = iter + 1;[Kt,Elem] = tangStiffMtxUL(Model,Anl,Elem,U);   %每次迭代更新，标准牛顿拉普森方法，Elem.ke更新
% Tangent and residual increments of displacementsd_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 stepD_lr = D_lr + d_lr;D_U  = D_U  + d_U;
% Total values of load ratio and displacementslr = lr + d_lr;U  = U  + d_U; end%----------------------------------------------------------------------------------  

  
（内力计算）function [F,Elem] = intForcesUL(Model,Elem,U,D_U,update_angle)% Initialize global vector of internal forcesF = zeros(Model.neq,1);for i = 1:Model.nel% Lengths: Beginning of step, current, and step incrementL_1  = Elem(i).L_1;%每个增量步开始时进行更新L_c  = elemLength(Elem(i),U);D_L  = L_c - L_1;
% Rigid body rotation from step beginning and current anglerbr   = elemAngleIncr(Elem(i),U,D_U);%rbr带有符号，增加正值，减少负值angle = Elem(i).angle_1 + rbr;
% Update element angleif (update_angle)Elem(i).angle = angle;end
% Rotation matrix from local to global coordinate systemc = cos(angle);s = sin(angle);rot = [ c -s  0  0  0  0;s  c  0  0  0  0;0  0  1  0  0  0;0  0  0  c -s  0;0  0  0  s  c  0;0  0  0  0  0  1 ];% Relative rotations（变形转角）r1 = D_U(Elem(i).n1.dof(3)) - rbr;r2 = D_U(Elem(i).n2.dof(3)) - rbr;% Vector of local displacementsdl = [0; 0 ; r1; D_L; 0; r2];%没有y向变形，只有轴向和转动变形？？？% Increment of internal forces in local systemD_fl = Elem(i).ke * dl;
% Total internal forces in local systemfl = Elem(i).fi_1 + D_fl;
% Store total internal forces in local systemElem(i).fi = fl;% Transform internal forces from local system to global systemfg = rot * fl;% Assemble element internal forces to global vectorgle    = Elem(i).gle;F(gle) = F(gle) + fg;endend