4、幾何非線性matlab代碼實現(xiàn)

（1）非線性求解算法流程總結(jié)

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



圖19

（2）matlab代碼實現(xiàn)

① 切線剛度矩陣的matlab代碼

根據(jù)公式（24），切線剛度矩陣中的彈性剛度矩陣實現(xiàn)代碼如下：



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



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



② 非線性系統(tǒng)求解的matlab代碼

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

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

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

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

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

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

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

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

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

之后，計算單元剛體旋轉(zhuǎn)角度（單元與全局坐標(biāo)系X軸的角度），即每個增量步開始時的角度angle_1+剛體旋轉(zhuǎn)增量 rbr。當(dāng)前角度在單元結(jié)構(gòu)體中更新。

內(nèi)力增量矢量D_f1彈性剛度矩陣與相對位移矢量（變形矢量）的乘積計算得出的。此增量添加到總內(nèi)力矢量，獲得局部坐標(biāo)系中當(dāng)前構(gòu)型下的總內(nèi)力fl存儲在單元數(shù)據(jù)結(jié)構(gòu)中，用于計算切線剛度矩陣（幾何剛度矩陣）。通過將局部坐標(biāo)系中的內(nèi)力矢量乘以旋轉(zhuǎn)矩陣，獲得全局坐標(biāo)系下的內(nèi)力矢量。最后，使用單元索引矢量將單元內(nèi)力插入到全局矢量中。

（預(yù)測階段代碼）
function Result = solve(Model,Anl,Elem,Pref)
% 初始化節(jié)點位移和荷載比
U = zeros(Model.neq,1);
lr = 0;
% 初始化result結(jié)構(gòu)體，并插入初始平衡點
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;   
% 更新參考構(gòu)型下的UL方程（L_1根據(jù)上一增量步結(jié)束時的位移U進(jìn)行更新；angle_1隨著迭代步更新并將上一增量步結(jié)束的angle作為該步angle_1;內(nèi)力fi_1同angle）
% 此處更新的L_1,angle_1和phi_1會在求解內(nèi)力的函數(shù)中用到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
% 切線剛度矩陣（根據(jù)Ui-1更新Ki0)
[Kt,Elem] = tangStiffMtxUL(Model,Anl,Elem,U) ;  %%根據(jù)i-1增量步的位結(jié)束移向量更新節(jié)點位置后得到的單元剛度矩陣（單元長度，單元角度），但返回的Elem只更新了Ke用來接下來計算內(nèi)力
% 預(yù)測解位移的切線增量delta_Ui1
d_Ut = solveLinearSystem(Model,Kt,Pref);%delta_Ui1預(yù)測階段的位移的切線增量，接下來是計算預(yù)測階段的荷載比增量（增量步是否為第一增量步對應(yīng)的荷載比增量公式不同）    
if (step == 1)
s = 1;       %第一增量步的預(yù)測階段的增量的符號默認(rèn)為正
d_lr = Anl.init_incr;   %delta_lamda_11 人為規(guī)定預(yù)測階段d_lr    
n2 = d_Ut'*d_Ut;%第一增量步中位移切線增量的范數(shù)的平方值,會在后續(xù)分析步計算GSP
else 
GSP = n2 / (d_Ut'*d_Ut_1);%公式（41）    % Generalized Stiffness Parameter   
if (GSP < 0)
s = -s; %公式（42） 增量步符號
end        
% 增量步調(diào)整系數(shù)J
% J = sqrt((Anl.trg_iter + 1) / (iter + 1));%目標(biāo)迭代次數(shù)/上一增量步步的迭代次數(shù)（公式（40））必須將其增加 1 以避免被零除，因為預(yù)測的解對應(yīng)于零次迭代
J = 1;      
%% 預(yù)測階段荷載比增量  采用圓柱弧長法 
% 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    
% 荷載比要小于規(guī)定值（==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;%保證恰好達(dá)到最終的荷載比增量1實現(xiàn)荷載的全部施加
end   
d_U = d_lr * d_Ut;%根據(jù)公式（36）預(yù)測階段假定殘差為零因此，只有前一項切線增量   
% Store predicted results
d_lr_1 = d_lr; d_Ut_1 = d_Ut; d_U_1  = d_U;       
% 前步的載荷比和位移增量（后續(xù)會隨著迭代不斷更新疊加）
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,并根據(jù)D_U對Ele.angle和Ele.fi進(jìn)行更新
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);   %每次迭代更新，標(biāo)準(zhǔn)牛頓拉普森方法，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);        


%如果增量步過大可能會出現(xiàn)復(fù)數(shù)
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
%----------------------------------------------------------------------------------

（內(nèi)力計算）
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;%每個增量步開始時進(jìn)行更新
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帶有符號，增加正值，減少負(fù)值
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（變形轉(zhuǎn)角）
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向變形，只有軸向和轉(zhuǎn)動變形？？？
% 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

審核編輯：劉清