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

2016-11-16 by:CAE仿真在線來(lái)源:互聯(lián)網(wǎng)

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

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

圖1 計(jì)算簡(jiǎn)圖

圖2 跨中節(jié)點(diǎn)荷載—位移曲線對(duì)比

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

圖3 Lee's Frame 簡(jiǎn)圖

圖4 Lee's frame變形動(dòng)畫(huà)

圖2 加載節(jié)點(diǎn)荷載—位移曲線對(duì)比

程序核心部分:

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