4、幾何非線性matlab代碼實現

（1）非線性求解算法流程總結

在前面的部分中，描述了通過增量迭代數值方法求解線性化平衡方程組的過程。通用求解算法的流程圖如圖19所示。



圖19

（2）matlab代碼實現

① 切線剛度矩陣的matlab代碼

根據公式（24），切線剛度矩陣中的彈性剛度矩陣實現代碼如下：



根據公式（25），切線剛度矩陣中的幾何剛度矩陣實現代碼如下，注釋掉的部分為應變張量的高階項，可根據需要決定計算與否：



在上述彈性剛度矩陣和幾何剛度矩陣得到后，單元的切線剛度矩陣實現代碼為：



② 非線性系統求解的matlab代碼

主函數首先初始化節點位移向量 U、載荷比值 lr。然后，增量迭代過程以一個while循環開始，該循環一直運行直到達到最大步數。增量步數計數器在增量循環開始時立即遞增。預測階段從更新UL公式的參考構型開始。

在當前的平衡構型下，即在上一步結束時獲得的單元內力和節點位移基礎上，計算預測階段的切線剛度矩陣 Kt。因此，tangStiffMtxUL()函數接收單元結構體數據的向量和節點位移的向量，兩者都使用上一步的結果進行了更新。然后，切線剛度矩陣用于計算位移的切線增量d_Ut，通過使用參考載荷矢量求解線性系統，求解線性系統函數solveLinearSystem()。

下一步是計算預測階段荷載比增量。當增量步是第一步時，預測增量的初始符號s設置為正數，且負載率增量 d_lr認為指定，另外存儲第一步中位移切線增量d_Ut的范數的平方值，n2。該值用來在后續步驟中計算GSP。對于剩余的增量步（增量步>1），GSP根據公式(41)計算，增量符號在根據公式（42）進行調整，即每次 GSP 值為負時必須反轉預測階段荷載比增量的符號。

之后，根據式(40)獲得增量大小的調整因子J。需要指出表達式中的變量 iter存儲上一步執行的迭代次數，必須將其增加1以避免被零除，因為預測的解對應于零次迭代。如果用戶選擇以恒定增量執行分析，則調整因子J的值設置為一。最后，根據公式（38）計算預測的荷載比增量。得到荷載率預測增量后，根據式(36)計算位移預測增量d_U（僅使用位移的切線增量），并對當前增量步的載荷比和位移增量（D_lr 和 D_U）進行更新。

迭代校正階段，首先迭代次數的計數器iter被初始化為零，因為假設第一次校正迭代僅在預測階段的收斂之后才開始。然后迭代循環開始于一個while循環中，該循環運行直到達到最大迭代次數或者不收斂時停止。

在迭代循環開始時，載荷比和節點位移的總值用上次迭代獲得的修正增量更新，或者如果剛進入循環則用預測階段的增量更新。外部節點力矢量 P 很容易從總載荷比lr中獲得，而內部節點力矢量 F 需要在通過intForcesUL()函數針對當前節點位移進行計算。然后通過外力和內力矢量之間的差獲得殘余力矢量 R。

基于殘余力范數的標準檢查當前迭代的收斂性，如公式(43) 所示。如果與參考荷載向量的比足夠小，則迭代收斂，算法退出循環進入下一個增量步。否則，算法繼續進行下一次迭代校正并立即遞增迭代次數inter。 荷載和位移增量校正迭代前先要選擇迭代方案。如果選擇標準的 Newton-Raphson 迭代，將使用更新后的構型來計算新的切線剛度矩陣，即使用在上一次迭代結束時獲得的單元內力和節點位移。否則，如果選擇修正Newton-Raphson 迭代法，則繼續使用在預測階段的得到的切線剛度矩陣。

位移校正的切線增量d_Ut和殘余增量d_Ur分別使用參考載荷Pref和殘余力R用線性系統計算。荷載比比的迭代校正選用恒定圓柱弧長法的表達式，由方程式(45)給出。最終根據方程(36)得到節點位移的迭代修正。一旦獲得載荷和位移的修正值，運用其增量值來更新前增量步。最后，更新總的位移和荷載。

對于給定的位移計算內力基于UL公式通過函數intForcesUL()實現，首先初始化內力的全局向量，然后，對所有單元執行循環，計算局部坐標系下每個單元的內力矢量，將其轉換為全局坐標系下的內力矢量，并將其組裝至全局矢量中。在循環開始時，計算單元伸長率D_L，是通過當前單元長度 L_c 與參考配置（每個增量步開始時）的長度 L_1 之間的差值獲得的。

之后，計算單元剛體旋轉角度（單元與全局坐標系X軸的角度），即每個增量步開始時的角度angle_1+剛體旋轉增量 rbr。當前角度在單元結構體中更新。

內力增量矢量D_f1彈性剛度矩陣與相對位移矢量（變形矢量）的乘積計算得出的。此增量添加到總內力矢量，獲得局部坐標系中當前構型下的總內力fl存儲在單元數據結構中，用于計算切線剛度矩陣（幾何剛度矩陣）。通過將局部坐標系中的內力矢量乘以旋轉矩陣，獲得全局坐標系下的內力矢量。最后，使用單元索引矢量將單元內力插入到全局矢量中。

（預測階段代碼）
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 
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 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：圓柱弧長法
% 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）預測階段假定殘差為零因此，只有前一項切線增量   
% Store predicted results
d_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 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進行更新
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
%----------------------------------------------------------------------------------

（內力計算）
function [F,Elem] = intForcesUL(Model,Elem,U,D_U,update_angle)
% Initialize global vector of internal forces
F = zeros(Model.neq,1);
for i = 1:Model.nel
% Lengths: Beginning of step, current, and step increment
L_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 angle
rbr   = elemAngleIncr(Elem(i),U,D_U);%rbr帶有符號，增加正值，減少負值
angle = Elem(i).angle_1 + rbr;


% Update element angle
if (update_angle)
Elem(i).angle = angle;
end


% Rotation matrix from local to global coordinate system
c = 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 displacements
dl = [0; 0 ; r1; D_L; 0; r2];%沒有y向變形，只有軸向和轉動變形？？？
% Increment of internal forces in local system
D_fl = Elem(i).ke * dl;


% Total internal forces in local system
fl = Elem(i).fi_1 + D_fl;


% Store total internal forces in local system
Elem(i).fi = fl;
% Transform internal forces from local system to global system
fg = rot * fl;
% Assemble element internal forces to global vector
gle    = Elem(i).gle;
F(gle) = F(gle) + fg;
end
end

審核編輯：劉清