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
審核編輯:劉清
-
matlab
+關注
關注
185文章
2974瀏覽量
230424 -
GSPM
+關注
關注
0文章
2瀏覽量
6120 -
MATLAB編程
+關注
關注
1文章
13瀏覽量
8428
原文標題:SimPC博士:幾何非線性有限元基本原理及matlab編程(下)
文章出處:【微信號:sim_ol,微信公眾號:模擬在線】歡迎添加關注!文章轉載請注明出處。
發布評論請先 登錄
相關推薦
評論