【轉載】弧長法(Riks method)通用求解程序（Code by myself）

2016-11-16 by:CAE仿真在線來源:互聯網

算例1. 如圖1所示的平面桿系結構,頂點受到豎直向下的力P作用,用本程序(Riks method)進行計算,并將計算結果與精確解進行比較,如圖2所示,通過對比可以說明本程序是正確的。

弧長法(Riks method)通用求解程序(Code by myself) - Yusd - 玉石殿

圖1 計算簡圖

圖2 跨中節(jié)點荷載—位移曲線對比

算例2:圖3是經典的Lee's frame簡圖,一個在端部正交的鉸接約束平面剛架,在距離正交點一定距離處有集中力F作用。之所以稱其為經典算例是因為它的荷載位移曲線同時集中了跳躍(snap-through )和回彈(snap-back)現象,傳統(tǒng)的求解策略根本無法對其進行荷載—位移路徑跟蹤,在此,弧長法展現了很大優(yōu)勢,圖4是運用本程序得到剛架的變形動畫,圖5是加載點的荷載位移曲線,并將其與ABAQUS計算結果進行對比,通過對比表明該程序的是正確的。

圖3 Lee's Frame 簡圖

圖4 Lee's frame變形動畫

圖2 加載節(jié)點荷載—位移曲線對比

程序核心部分:

讀取數據文件(節(jié)點、單元、約束、截面屬性、參考力、控制弧長、最大控制參量)
while 控制參量(如位移、最大荷載因子) < 最大控制參量
計算當前切線剛度矩陣 K_Global
計算參考位移 X_Ref= 參考力\K_Global
計算初始荷載因子 lamda0=Arclength/sqrt(1+X_Ref'*X_Ref);
判定初始荷載因子方向 +/- lamda0
更新節(jié)點坐標,更新外力
計算當前節(jié)點反力
計算節(jié)點不平衡力Val
while norm(Val)>1e-6
計算不平衡力產生的位移X_Val
計算荷載因子修正參數delta=X_Val'*X_Ref/(1+X_Ref'*X_Ref);
修正荷載因子lamda1=lamda0-delta;
更新初始荷載因子lamda0=lamda1;
更新節(jié)點坐標,更新外力
計算當前節(jié)點反力
計算節(jié)點不平衡力Val