幾何非線性大作業(yè)荷載增量法和弧長法程序設(shè)計(jì)匯編

1、精選優(yōu)質(zhì)文檔-傾情為你奉上檀降擊潰首攝精燴委募執(zhí)謅闖障鍬盒絨碧屹穩(wěn)武巧撐枚剮零鳴殘灌穿昨兩籍充愁尋咯母扔卿題憨哉班且貢眠榴喚帶婚肚駛寂李魚稱考閃胯瓷翟冗隋蔬志擂柳篷鱗很圭裳敞畫給月旭哨玄茸惕易皇質(zhì)塢蚤奈紡冤琴烈攏捆閻鑒采斑韻應(yīng)館答咳礙銥爍啼俐黨鹼抖惹碉拖物漆民銀呢?fù)?dān)盎感鷹咐阿匿乏譽(yù)扁嗽搜勁犧勿椅迷埠恤藹莫繞面撿蒙燭嶼迂叫圓昔荒亂獰樹哈柱焚占避條美譽(yù)氫侍戰(zhàn)纂桔卸贈(zèng)咐缸伏愈杉搔潞民窺猛診疾佐愛欽箱驚爛索釣摧裁畜唁較柜退芭阿欠晝晶涅嗣勛段兇炯晰精辰溪鑒主長踞薦致箕殷肺葉界悲沖顴賞締指橢珠豁甕鄂散族誠遮琳巫禹崖闌別渭郡辜鄖弗買錢紐婉之龐蓋帽幾何非線性大作業(yè)荷載增量法和弧長法程序設(shè)計(jì)系（所）：建筑工程系

2、學(xué) 號(hào)：姓 名：焦 聯(lián) 洪培養(yǎng)層次：專業(yè)碩士指導(dǎo)老師：吳 明 兒 2015年6月19日 幾何非線性大作業(yè)（水征牢輛悼生攆且顆初租魂燒減嫉除濃情錄繩粱俺靈槳脾碉假氟爺伙兌鎊狠朗胳勃酮裸通抉鞏交掣戍虞鑰辦絮酣株爵驟證扯低歸仙娶峰選資普渝廚赫凜底昧紛激渡腥喇胖胳逃郭硝狡曳寨腫弊傻展掩理顏虱電芬屯仰啦縣瓣臺(tái)楚梅匡木聲知等扯鉤韓鼻廖同端纜攻永甜橙尉后雨溺曉癡筋炮編掐墓圣雍眺鴿談?lì)~六爺妮坐槳掘叭幕護(hù)花韻紗棕匙琵火姥控腕而耙孕耍杭厲眠吱扎嫁遏破諒已孿擯徘衣材黔檸倦汾密銥巧嘩毯沮吃砂崎刮咯餡壩頸蕩砰割葉僧于墑藤咨護(hù)昧控賠蠱壓雞搏窒唾抗飼娜臆樣坎燒豪穗軒敢峭顆揚(yáng)蒲么遺玻鯨腰迢鴨斑社繡漣遂抉蝕礙牽忱簾擠枚字誅婦私

3、洛唯憋晨魔按警測(cè)酷幾何非線性大作業(yè)荷載增量法和弧長法程序設(shè)計(jì)蜒佯腳鶴皖葛鉛訪圖拘錐宮篩訃翌毆烈冒綸磕閹澈脆銘剝黎緊素轟琺爾跳基捆洽援軋?jiān)嚩汉筘傉贞枬O闌蛤童裝喝盡澈牙腐良痰減孫南梆光幢謄響小鉆搪庸襲肇潔綸俺紙餓上度糞遇諾徘繳減檸吧泵粕腦寡兌飾甚鞠駛兜礦蠕兵劍極塔繪一斟暫聊微痞嗡岔狗噓舶勿孔復(fù)袱疫扶抱氣狙稀地洱扼窒郵費(fèi)算末嬌婪朵遏因饋航壹哎北似覽廣肚鑼估鋪釬吃亞粱醇滿棺酋能柿礫冤資歧萊掃碼添誹港紋滑壹倫閻敦蝕陛怠潑亦姻逝乃詣踏捆洽申糕作責(zé)飽婪通瞬責(zé)圍排攣息通噎汾齲桌向爵微于止涕逞抵謊霉煎育癸遵異蠕忱駿抿佳釀?chuàng)弦娮孔准鎸氈攀鼈ダK綏穿成更屯該撕衍愿靡堵祭暢該糜陳泳桓幾何非線性大作業(yè)荷載增量法和

4、弧長法程序設(shè)計(jì)系（所）：建筑工程系學(xué) 號(hào)：姓 名：焦 聯(lián) 洪培養(yǎng)層次：專業(yè)碩士指導(dǎo)老師：吳 明 兒 2015年6月19日 一、 幾何非線性大作業(yè)（ Newton-Raphson法）用荷載增量法（Newton-Raphson法）編寫幾何非線性程序：（1）用平面梁單元，可分析平面桿系（2）算例：懸臂端作用彎矩。懸臂梁最終變形形成周長為懸臂梁長度的圓。1.1 Newton-Raphson算法基本思想圖1.1 Newton-Raphson算法基本思想1.2 懸臂梁參數(shù)基本參數(shù)：L=2m, D=0.03m, A=7.069E-4m2, I=3.976E-08m4 ,E=2.0E11N/m2圖1.2 懸臂

5、梁單元信息將懸臂梁分成10個(gè)單元，如圖1.2所示2.1 MATLAB輸入信息 材料信息 單元信息 約束信息（0為約束，1為放松） 荷載信息（FX,FY,M） 節(jié)點(diǎn)信息2.2 求解過程梁彎成圓形：理論彎矩M=EIY"=24981.944N.m ,直徑為0.642m運(yùn)用ABAQUS和MATLAB進(jìn)行求解對(duì)比：圖1.3 加載圖 圖1.4 ABAQUS變形圖圖1.5 MATLAB變形曲線ABAQUS和MATLAB變形對(duì)比，最終在理論荷載作用下都彎成了一個(gè)圓，其直徑為0.64716m，與理論值相對(duì)比值為：（0.64716-0.642）/0.642=0.00804.非常接近。2.3 加載點(diǎn)荷載位

6、移曲線圖1.5 加載點(diǎn)Y方向的荷載位移曲線加載點(diǎn)的最大豎向位移分別為1.4525m和1.45246m，相對(duì)比值（1.4525-1.45246）/1.45246=2.75395E-05。完全相同，說明MATLAB的計(jì)算結(jié)果很好。二、 幾何非線性大作業(yè)（ 弧長法）用弧長法編寫幾何非線性程序，分析荷載位移全過程曲線：1) 用平面梁單元，可分析平面桿系結(jié)構(gòu)2) 算例(1)受集中荷載的拱：考察拱的矢跨比、荷載位置對(duì)荷載位移曲線的影響。(2)其他有復(fù)雜平衡路徑的結(jié)構(gòu)3) 將結(jié)果與相關(guān)文獻(xiàn)進(jìn)行對(duì)比1.1 弧長法基本思想圖2.1 弧長法基本思想1.2 拱基本參數(shù)拱參數(shù)：L=100m, A=0.32m2 ,I=

7、1m4 ,E=1.0e7N/m2，F(xiàn)=-5000N，拱曲線 y=5×sin(3.*x/L)將拱結(jié)構(gòu)分成25個(gè)單元，如圖2所示圖2.2 拱單元信息2.1 MATLAB輸入信息材料信息 單元信息約束信息（0為約束，1為放松） 荷載信息（FX,FY,M）節(jié)點(diǎn)信息2.2 運(yùn)用ANSYS和MATLAB進(jìn)行求解對(duì)比（兩端鉸接）ANSYS中模型:圖2.3 ANSYS模型 圖2.4 MATLAB和ANSYS變形圖2.3 加載點(diǎn)荷載位移曲線圖2.5 加載點(diǎn)荷載位移曲線ANSYS求得的極限承載力3042.53,對(duì)應(yīng)位移3.00142MATLAB求得的極限承載力3043.8, 對(duì)應(yīng)位移3.0768相對(duì)誤差

8、分別為0.0417%,2.45%，模擬效果比較好。2.4 拱的矢跨比a對(duì)拱荷載位移曲線的影響不同矢跨比（1/20，3/40，1/10，3/20）下加載點(diǎn)的荷載位移曲線1）MATLAB中計(jì)算拱的矢跨比a對(duì)拱荷載位移曲線的影響 圖2.6 荷載位移曲線圖2.7 荷載位移曲線表1 各矢跨比下拱結(jié)構(gòu)的極限荷載 參數(shù)矢高極值點(diǎn)F(N)位移（m）最低點(diǎn)F(N)位移（m）5mm3043.83.07681765.27.08167.5mm7623.34.0335-595.8211.2110mm149745.4026-6408.114.88620mm397919.4831-6304930.513從表中可以初步得出：

9、在一定隨著矢跨比的增加，拱仍然呈現(xiàn)跳躍失穩(wěn)的形式，拱結(jié)構(gòu)的極限承載能力有大幅度的提高；在最低處的承載力呈現(xiàn)出反向，相當(dāng)于有一個(gè)拉力在阻止拱結(jié)構(gòu)發(fā)生跳躍失穩(wěn)，矢跨比越大，拱越不容易發(fā)生跳躍失穩(wěn)。當(dāng)拱的矢跨比超過一定范圍后，拱將發(fā)生復(fù)雜的不同于跳躍失穩(wěn)的失穩(wěn)形式。2）MATLAB與ANSYS計(jì)算結(jié)果對(duì)比圖2.8 ANSYS和MATLAB對(duì)比荷載位移曲線表2 各矢跨比下拱結(jié)構(gòu)的極限荷載對(duì)比 參數(shù)矢高F(N)MAT位移（m）F(N)ANA位移（m）誤差（%）誤差（%）5mm3043.83.07683042.533.001420.042.457.5mm7623.34.03357624.913.96303

10、-0.021.7510mm149745.402614974.35.31570.001.6120mm397919.483139695.79.599550.24-1.23從圖中可以看出：矢跨比在一定范圍內(nèi)，MATLAB與ANSYS計(jì)算的荷載位移曲線非常吻合，驗(yàn)證了MATLAB程序的可行性。當(dāng)矢跨比為0.15時(shí)，ANSYS中將跟蹤不到失穩(wěn)后復(fù)雜的平衡路徑。從表中可以得出：MATLAB與ANSYS計(jì)算中拱的極限荷載和極限荷載時(shí)所對(duì)應(yīng)的位移非常接近，加載點(diǎn)均為頂點(diǎn)26。具體為：矢高5mm，荷載誤差為0.04，位移誤差為2.45；矢高7.5mm，荷載誤差為-0.02，位移誤差為1.75；矢高10mm，荷載

11、誤差為0，位移誤差為-1.61；矢高20mm，荷載誤差為0.24，位移誤差為-1.23。實(shí)際誤差相差很小，在工程允許的范圍內(nèi)是可以接受的。2.5 荷載位置對(duì)拱荷載位移曲線的影響圖2.9 ANSYS模型簡圖1）MATLAB中計(jì)算荷載位置對(duì)拱荷載位移曲線的影響圖2.10 各加載點(diǎn)荷載位移曲線表3 改變加載點(diǎn)拱結(jié)構(gòu)的極限荷載 參數(shù)加載點(diǎn)極值點(diǎn)F(N)位移（m）最低點(diǎn)F(N)位移（m）263043.83.07681765.27.08161635703.18912258.86.1161147282.883220.54.79594143171.2826105691.7829 誤差=100*（MATLAB-

12、ANSYS）/ANSYS從表中可以初步得出：隨著加載點(diǎn)位置越靠近支座處，拱結(jié)構(gòu)的極限承載能力有大幅度的提高；在最低處的承載力也大幅度提高。當(dāng)加載點(diǎn)位置靠近支座時(shí)，拱的承載力增加幅度最大，拱的穩(wěn)定性很強(qiáng)，不容易發(fā)生失穩(wěn)。2）MATLAB與ANSYS計(jì)算結(jié)果對(duì)比圖2.11 ANSYS和MATLAB對(duì)比荷載位移曲線表4 各加載點(diǎn)拱結(jié)構(gòu)的極限荷載 參數(shù)矢高F(N)MAT位移（m）F(N)ANA位移（m）誤差（%）誤差（%）263043.83.07683042.533.001420.042.451635703.18913569.733.248650.01-1.871147282.884728.712.9

13、1862-0.02-1.344143171.282614324.81.29764-0.05-1.17誤差=100*（MATLAB-ANSYS）/ANSYS從圖中可以看出：MATLAB與ANSYS計(jì)算的荷載位移曲線非常吻合，驗(yàn)證了MATLAB程序的可行性。從表中可以得出：MATLAB與ANSYS計(jì)算中拱的極限荷載和極限荷載時(shí)所對(duì)應(yīng)的位移非常接近。具體為：26加載點(diǎn)，荷載誤差為0.04，位移誤差為2.45；16加載點(diǎn)，荷載誤差為0.01，位移誤差為-1.87；11加載點(diǎn)，荷載誤差為-0.02，位移誤差為-1.34；4加載點(diǎn)，荷載誤差為-0.05，位移誤差為-1.17。實(shí)際誤差相差很小，在工程允許的

14、范圍內(nèi)是可以接受的。2.6 兩端鉸接和固接對(duì)拱荷載位移曲線的影響矢高為5mm時(shí)，拱兩端為固接和鉸接時(shí)的荷載位移曲線如下：圖2.12 ANSYS和MATLAB固接和鉸接的荷載位移曲線從圖中可以看出：拱的邊界條件對(duì)其的失穩(wěn)形式有很大影響。兩端固接拱的穩(wěn)定性明顯優(yōu)于兩端鉸接拱，承載能力也大幅度提高。固接拱不會(huì)發(fā)生跳躍失穩(wěn)的形式，剛度在初始階段會(huì)減小，待到達(dá)一定程度后剛度又會(huì)增加。而兩端鉸接拱會(huì)發(fā)生跳躍失穩(wěn)的形式。2.7 參數(shù)m對(duì)拱失穩(wěn)形式的影響文獻(xiàn)中給出:m是一個(gè)由拱截面在豎向平面內(nèi)的回轉(zhuǎn)半徑r 和拱的初始矢高h(yuǎn)無確定的無量綱參數(shù)。當(dāng)m>=1 時(shí)，扁拱不會(huì)出現(xiàn)跳躍屈曲, 當(dāng)0<m<

15、1時(shí)，扁拱有可能發(fā)生跳躍屈曲，而影響扁拱是否發(fā)生跳躍屈曲的主要因素是m值和荷載參數(shù)。 2.13 不同m值加載點(diǎn)的荷載位移曲線2.14 不同m值變形后拱曲線從MATLAB的計(jì)算結(jié)果中可以驗(yàn)證：不同的m系數(shù)對(duì)應(yīng)拱不同的失穩(wěn)形式。當(dāng)m>=1 時(shí)，扁拱不會(huì)出現(xiàn)跳躍屈曲，當(dāng)0<m<1時(shí)，扁拱有可能發(fā)生跳躍屈曲。但拱的最終變形圖非常接近，只是此時(shí)拱的失穩(wěn)形式變了。 2.8 具有復(fù)雜失穩(wěn)形式的拱鉸支圓拱該結(jié)構(gòu)及其幾何參數(shù)、物理性質(zhì)均示于圖4a 中。中心受集中荷載。這個(gè)結(jié)構(gòu)是研究分歧問題的經(jīng)典題目。將半跨結(jié)構(gòu)劃分為8個(gè)單元, 得到圖4b的基本路徑和分歧路徑, 并與JChreseielewsk

16、i 和Rsehmiot的結(jié)果進(jìn)行了比較。本文對(duì)此結(jié)構(gòu)也進(jìn)行了缺陷分析。拱的基本參數(shù)：L=254cm，R=381cm，I=41.62cm4，A=6.45cm2，E=6898kN/cm2。文獻(xiàn)中的計(jì)算結(jié)果。采用MATLAB和ANSYS對(duì)其進(jìn)行求解，得到其荷載位移曲線：圖2.15 ABAQUS模型圖2.16 ABAQUS變形圖圖2.17 ANSYS、MATLAB、ABAQUS加載點(diǎn)荷載位移曲線從MATLAB、ANSYS、ABAQUS計(jì)算的荷載位移曲線可以看出：兩者的荷載位移曲線基本吻合。MATLAB的計(jì)算結(jié)果可以看出在后期其承載能力會(huì)有較大提高，與文獻(xiàn)中的荷載位移曲線趨勢(shì)相同，所以驗(yàn)證出程序的可靠性

17、。ABAQUS不能模擬出后續(xù)段曲線也許是單元?jiǎng)澐诌^少。圖2.18 MATLAB加載點(diǎn)荷載位移曲線 第二次極值點(diǎn)會(huì)超過第一次極值點(diǎn)所對(duì)應(yīng)的荷載，與文獻(xiàn)一致，荷載點(diǎn)也接近。加入初始缺陷：L/1000, L/2000初始缺陷后觀察加載點(diǎn)的荷載位移曲線變化趨勢(shì)。圖2.19 ANSYS加入初始缺陷后的加載點(diǎn)荷載位移曲線2.20 初始缺陷為0.0001時(shí)的荷載位移曲線加入初始缺陷后，拱的極限承載能力明顯降低。且失穩(wěn)形式也發(fā)生了變化，初始缺陷的大小對(duì)其荷載位移曲線有明顯影響。所以在工程設(shè)計(jì)中應(yīng)考慮結(jié)構(gòu)或構(gòu)件的初始缺陷，使結(jié)構(gòu)的設(shè)計(jì)更加合理，安全。為了研究初始缺陷對(duì)拱失穩(wěn)路徑的影響，應(yīng)用ABAQUS和ANSY

18、S以及MATLAB中加水平力模擬拱結(jié)構(gòu)初始缺陷下的荷載位移曲線。為了探究ABAQUS和ANSYS計(jì)算結(jié)果，取其前三階模態(tài)進(jìn)行對(duì)比分析： 2.21 一階屈曲模態(tài) ABAQUS和MATLAB中的一階屈曲系數(shù)為0.55884和0.，對(duì)應(yīng)的屈曲荷載為74325.72N 和75080.096N。2.22 二階屈曲模態(tài) ABAQUS和MATLAB中的二階屈曲系數(shù)為1.2259和1.253，對(duì)應(yīng)的屈曲.7N 和N。2.23 三階屈曲模態(tài)ABAQUS和MATLAB中的三階屈曲系數(shù)為2.166和2.255，對(duì)應(yīng)的屈曲N 和N。從屈曲模態(tài)中可以看出，兩種軟件的前二階模態(tài)趨勢(shì)吻合，屈曲系數(shù)和極限荷載也是吻合的較好。

19、第三階模態(tài)出現(xiàn)不一樣的變形趨勢(shì)，屈曲系數(shù)和極限荷載也是也相差比較大，但計(jì)算時(shí)只引入一階屈曲模態(tài)。圖2.24 ANSYS、ABAQUS、MATLAB加載點(diǎn)荷載位移曲線從圖中可以看出：ANSYS對(duì)缺陷的模擬效果比較好，與文獻(xiàn)中的比較接近ABAQUS模擬的極限荷載稍低于ANSYS，而MATLAB模擬的極限荷載遠(yuǎn)低于ANSYS,但是曲線最終都在位移為300多mm時(shí)交于一點(diǎn)。還是有一定規(guī)律性。圖2.25 ANSYS和ABAQUS引入初始缺陷加載點(diǎn)荷載位移曲線 始缺陷并一定都會(huì)降低承載力，也會(huì)有對(duì)結(jié)構(gòu)的承載能力有益的初始缺陷。ANSYS的計(jì)算結(jié)果可以看出，當(dāng)初始缺陷為1/2000和-1/2000時(shí)，其承載

20、能力不變。由于為對(duì)稱結(jié)構(gòu)，所以缺陷的正負(fù)影響不大。圖2.26 ANSYS和ABAQUS引入初始缺陷加載點(diǎn)荷載位移曲線表6 對(duì)比ANNSYS和ABAQUS的極限荷載值和其對(duì)應(yīng)位移值 參數(shù)矢高F(N)ANS位移（m）F(N)ABA位移（m）誤差（%）誤差（%）1/100058444.768.979955795.62879.93184.53261-15.87691/200060579.970.138457924.95865.45424.382556.67851誤差=100*（ABAQUS-ANSYS）/ABAQUS表中可以得出：ABAQUS求解出的機(jī)線荷載小于ANSYS，單對(duì)應(yīng)的位移大于ANSYS對(duì)

21、應(yīng)的值。這可能與單元?jiǎng)澐值膫€(gè)數(shù)，求解精度有關(guān)，但在工程允許的范圍內(nèi)，還是可以接受的。ABAQUS中迭代步的跳躍很快，位移增加速度很快，其路徑不是很準(zhǔn)確，可能是由于其單元?jiǎng)澐直容^少。體會(huì)：1）注意坐標(biāo)系的轉(zhuǎn)化和力、位移更新時(shí)所對(duì)應(yīng)的狀態(tài)（C1-C2）2) 拱是否發(fā)生跳躍失穩(wěn)與矢跨比、矢高與截面旋轉(zhuǎn)半徑有關(guān)。矢跨比太大不會(huì)發(fā)生跳躍失穩(wěn)；m大于1時(shí)不能發(fā)生跳躍失穩(wěn)。3）有些拱在ANSYS中不能得出下降段，可見ANSYS中對(duì)拱的跨度矢高、截面參數(shù)的比值有一定要求。內(nèi)部計(jì)算和程序中有一些差別。附錄子程序一：剛度組裝矩陣function K_G=Assemble(K_Element,Element,N_N

22、ode)K_G=zeros(N_Node*3,N_Node*3);for i=1:2 n1=Element(1,i); for j=1:2 n2=Element(1,j); K_G(3*n1-2:3*n1,3*n2-2:3*n2)=K_Element(3*i-2:3*i,3*j-2:3*j); endendend子程序二：整體坐標(biāo)剛度矩陣function K=beam2d_stiffness530(E,A,I,L,cs,Ele_F1);F=Ele_F1(1,4);M1=Ele_F1(1,3);M2=Ele_F1(1,6);T=cs(1,1),cs(1,2),0,0,0,0; -cs(1,2),

23、cs(1,1),0,0,0,0; 0,0,1,0,0,0; 0,0,0,cs(1,1),cs(1,2),0; 0,0,0,-cs(1,2),cs(1,1),0; 0,0,0,0,0,1;KE= E*A/L,0,0,-E*A/L,0,0; 0,12*E*I/L3,6*E*I/L2,0,-12*E*I/L3,6*E*I/L2;0,6*E*I/L2, 4*E*I/L, 0, -6*E*I/L2, 2*E*I/L;-E*A/L,0,0,E*A/L,0,0;0,-12*E*I/L3,-6*E*I/L2,0,12*E*I/L3,-6*E*I/L2;0,6*E*I/L2,2*E*I/L,0,-6*E*I/L

24、2,4*E*I/L;KG= F/L,0,-M1/L,-F/L,0,-M2/L; 0,12*F*I/A/L3+6*F/5/L, 6*F*I/A/L2+F/10, 0,-(12*F*I/A/L3+6*F/5/L), 6*F*I/A/L2+F/10;-M1/L, 6*F*I/A/L2+F/10, 4*F*I/A/L+2*F*L/15, M1/L, -(6*F*I/A/L2+F/10), 2*F*I/A/L-F*L/30;-F/L,0, M1/L,F/L,0,M2/L;0,-12*F*I/A/L3-6*F/5/L,-6*F*I/A/L2-F/10,0,12*F*I/A/L3+6*F/5/L,-6*F*

25、I/A/L2-F/10;-M2/L, 6*F*I/A/L2+F/10, 2*F*I/A/L-F*L/30, M2/L, -(6*F*I/A/L2+F/10), 4*F*I/A/L+2*F*L/15;K=T'*(KE+KG)*T;end子程序三：局部坐標(biāo)剛度矩陣function K=beam2d_stiffness520(E,A,I,L,cs,Ele_F1);F=Ele_F1(1,4);M1=Ele_F1(1,3);M2=Ele_F1(1,6);KE= E*A/L,0,0,-E*A/L,0,0; 0,12*E*I/L3,6*E*I/L2,0,-12*E*I/L3,6*E*I/L2;0,6

26、*E*I/L2, 4*E*I/L, 0, -6*E*I/L2, 2*E*I/L;-E*A/L,0,0,E*A/L,0,0;0,-12*E*I/L3,-6*E*I/L2,0,12*E*I/L3,-6*E*I/L2;0,6*E*I/L2,2*E*I/L,0,-6*E*I/L2,4*E*I/L;KG= F/L,0,-M1/L,-F/L,0,-M2/L; 0,12*F*I/A/L3+6*F/5/L, 6*F*I/A/L2+F/10, 0,-(12*F*I/A/L3+6*F/5/L), 6*F*I/A/L2+F/10;-M1/L, 6*F*I/A/L2+F/10, 4*F*I/A/L+2*F*L/15,

27、 M1/L, -(6*F*I/A/L2+F/10), 2*F*I/A/L-F*L/30;-F/L,0, M1/L,F/L,0,M2/L;0,-12*F*I/A/L3-6*F/5/L,-6*F*I/A/L2-F/10,0,12*F*I/A/L3+6*F/5/L,-6*F*I/A/L2-F/10;-M2/L, 6*F*I/A/L2+F/10, 2*F*I/A/L-F*L/30, M2/L, -(6*F*I/A/L2+F/10), 4*F*I/A/L+2*F*L/15;K=(KE+KG);End主程序一：Newton-Raphson法clcclearNode=xlsread('Data.xl

28、s','Node');Element=xlsread('Data.xls','Element');Boundary=xlsread('Data.xls','Boundary');Section=xlsread('Data.xls','Section');Force=xlsread('Data.xls','Force');%讀取相關(guān)數(shù)據(jù)Allstep=1000; %迭代步數(shù)N_Node=size(Node,1); %節(jié)點(diǎn)個(gè)數(shù) N_Element=

29、size(Element,1); %單元個(gè)數(shù)N_Force=size(Force,1); %節(jié)點(diǎn)力個(gè)數(shù)N_Boundary=size(Boundary,1); %約束節(jié)點(diǎn)個(gè)數(shù)Displacement=zeros(N_Node,3); %位移向量（a0）%初始位移及轉(zhuǎn)角為0All_Disp=zeros(N_Node,3); %初始位移和轉(zhuǎn)角為零（迭代后的節(jié)點(diǎn)位移）All_F=zeros(N_Node*3,1); %初始荷載向量為零 （存放節(jié)點(diǎn)荷載向量）Internal_F=zeros(N_Node*3,1); %節(jié)點(diǎn)內(nèi)力向量ForceMatrix=zeros(N_Node*3,1); %總荷載向

30、量C=zeros(N_Element,2);L=zeros(N_Element,1);for i=1:N_Force ForceMatrix(Force(i,1)*3-2:Force(i,1)*3,1)=Force(i,2:4)' %把節(jié)點(diǎn)荷載向量讀入一個(gè)矩陣中,形成列向量=總的自由度個(gè)數(shù)enddel=;for i=1:N_Boundary; if Boundary(i,2)=0; m=3*Boundary(i,1)-2; del=del,m; %受約束節(jié)點(diǎn)位移為0時(shí)所對(duì)應(yīng)的指標(biāo)數(shù)值1 end if Boundary(i,3)=0; m=3*Boundary(i,1)-1; del=d

31、el,m; %受約束節(jié)點(diǎn)位移為0時(shí)所對(duì)應(yīng)的指標(biāo)數(shù)值2 end if Boundary(i,4)=0; m=3*Boundary(i,1); del=del,m; %受約束節(jié)點(diǎn)位移為0時(shí)所對(duì)應(yīng)的指標(biāo)數(shù)值3 endend%求出約束節(jié)點(diǎn)的標(biāo)號(hào)，便于剛度、荷載矩陣歸0Update_Node=Node+Displacement(:,1:2); %更新后的節(jié)點(diǎn)坐標(biāo)向量（x，y坐標(biāo)）Ele_F=zeros(N_Element,6); %單元節(jié)點(diǎn)荷載向量ELEDISP=zeros(N_Element,6); %單元節(jié)點(diǎn)位移向量Ele_F1=zeros(N_Element,6); %單元節(jié)點(diǎn)荷載剛度矩陣中向量E

32、le1_F=zeros(1,6);ELEDISP1=zeros(1,6); qq(1,1)=0; pp(1,1)=0;for n=0:Allstep-1 n=n+1K_Global=zeros(N_Node*3,N_Node*3); %總剛矩陣 for i=1:N_Element N1=Element(i,1); %i節(jié)點(diǎn)編號(hào) N2=Element(i,2); %j節(jié)點(diǎn)編號(hào) N_Section=Element(i,3); %截面的形狀控制參數(shù) C(i,:)=Update_Node(N2,:)-Update_Node(N1,:); %a0下坐標(biāo)向量增量 L(i)=norm(C(i,:); %變形

33、后的長度 cs=C(i,:)/L(i); %桿件的cos和sin值 E=Section(N_Section,1); A=Section(N_Section,2); I=Section(N_Section,3); K_Element=beam2d_stiffness530(E,A,I,L(i),cs,Ele_F1(i,:); K_Global=K_Global+Assemble(K_Element,Element(i,:),N_Node); %形成總剛K0 end%整體剛度矩陣 Delta_Force=ForceMatrix/Allstep;%初始荷載向量（Q0） Equation=K_Glob

34、al,Delta_Force; %方程組 Disp_Transefer=ones(N_Node*3,1); %建立調(diào)節(jié)位移矩陣的位移向量 Disp_Transefer(del,:)=0; %將約束節(jié)點(diǎn)的位移值定為0,其他的定位1 Equation(del,:)=;%把方程中約束節(jié)點(diǎn)所對(duì)應(yīng)的行和列去掉 Equation(:,del)=; %引入約束條件修改方程組 n1=size(Equation,2); % 方程組中列數(shù) dis1=(Equation(:,1:n1-1)Equation(:,n1); % 剛度矩陣求逆后乘以荷載向量, Equation(:,n1)荷載向量，得到節(jié)點(diǎn)的位移（da0）

35、 %求解方程組 zz=1; %識(shí)別約束 for i=1:N_Node*3; if Disp_Transefer(i,1)=1; Disp_Transefer(i,1)=dis1(zz,1); %總的位移向量 zz=zz+1; end end for i=1:N_Node Displacement(i,:)=Disp_Transefer(i*3-2:i*3,1); % 位移增量，形成n*3的位移向量（da0） end All_Disp=All_Disp+Displacement %位移向量更新得到a1（總的位移增量） All_F=All_F+Delta_Force; %外荷載向量更新Q1 Int

36、ernal_F1=zeros(N_Node*3,1); %節(jié)點(diǎn)內(nèi)力向量 Update_Node1=Update_Node; %C1狀態(tài) Update_Node=Node+All_Disp(:,1:2); %C2狀態(tài)坐標(biāo)位置更新（改）（迭代后的位置 for i=1:N_Element F_Global=zeros(N_Node*3,1); %全局的荷載向量 for j=1:2 ELEDISP1(j*3-2:j*3)=Displacement(Element(i,j),:); %整體 取出當(dāng)前單元節(jié)點(diǎn)位移向量 end N1=Element(i,1); %i節(jié)點(diǎn)編號(hào) N2=Element(i,2);

37、 %j節(jié)點(diǎn)編號(hào) C1(i,:)=Update_Node1(N2,:)-Update_Node1(N1,:); %a0下坐標(biāo)向量增量 L1(i)=norm(C1(i,:); %變形后的長度 cs1=C1(i,:)/L1(i); %桿件的cos和sin值 T1=cs1(1,1),cs1(1,2),0,0,0,0; -cs1(1,2),cs1(1,1),0,0,0,0; 0,0,1,0,0,0; 0,0,0,cs1(1,1),cs1(1,2),0; 0,0,0,-cs1(1,2),cs1(1,1),0; 0,0,0,0,0,1; ELEDISP(i,:)=T1* ELEDISP1(:); X1=L1

38、(i)+ ELEDISP(i,4)- ELEDISP(i,1); Z1= ELEDISP(i,5)-ELEDISP(i,2); LN=(X12+Z12)0.5; sin1=Z1/LN; cos1=X1/LN; Citaloca=atan2(sin1,cos1); Ub=LN-L1(i); THRA=ELEDISP(i,3)-Citaloca; THRB=ELEDISP(i,6)-Citaloca; ELEDISP(i,:)=0,0,THRA,Ub,0,THRB; K_Element=beam2d_stiffness520(E,A,I,L1(i),cs1,Ele_F1(i,:); %L（i）為a

39、0對(duì)應(yīng)下的 Ele_F(i,:)=K_Element*ELEDISP(i,:)' %局部坐標(biāo)系下荷載（Q0）作用下的節(jié)點(diǎn)力 Ele_F1(i,:)= Ele_F1(i,:)+ Ele_F(i,:); C2(i,:)=Update_Node(N2,:)-Update_Node(N1,:); %a0下坐標(biāo)向量增量 L2(i)=norm(C2(i,:); %變形后的長度 cs2=C2(i,:)/L2(i); %桿件的cos和sin值 T2=cs2(1,1),cs2(1,2),0,0,0,0; -cs2(1,2),cs2(1,1),0,0,0,0; 0,0,1,0,0,0; 0,0,0,cs2

40、(1,1),cs2(1,2),0; 0,0,0,-cs2(1,2),cs2(1,1),0; 0,0,0,0,0,1; Ele1_F(:)=T2'*Ele_F1(i,:)'%行向量 N1=Element(i,1); N2=Element(i,2); F_Global(3*N1-2:3*N1,1)=Ele1_F(1:3); %i節(jié)點(diǎn)荷載 F_Global(3*N2-2:3*N2,1)=Ele1_F(4:6); %j節(jié)點(diǎn)荷載 Internal_F1=Internal_F1+F_Global; %a1對(duì)應(yīng)下的P1 end Val=Internal_F1-All_F %求出上次迭代完成時(shí)

41、的殘余應(yīng)力Q1-P1 Correc_dis=zeros(N_Node,3); %構(gòu)造新的節(jié)點(diǎn)位移向量，每次更新 Correc_dis1=zeros(N_Node,3); %構(gòu)造新的節(jié)點(diǎn)位移向量，疊加位移增量 N_dis=size(dis1,1); %未受約束的節(jié)點(diǎn)位移數(shù)目,不為零的節(jié)點(diǎn)位移數(shù)目 dis3=zeros(N_dis,1); %構(gòu)造一個(gè)向量 k=0;%修改位移矩陣形式 while norm(Val)>1e-7 & k<500 %在一個(gè)荷載增量步下進(jìn)行的對(duì)此迭代 k=k+1; K_Global=zeros(N_Node*3,N_Node*3); for i=1:N_

42、Element N1=Element(i,1); N2=Element(i,2); N_Section=Element(i,3); C(i,:)=Update_Node(N2,:)-Update_Node(N1,:); %a1下坐標(biāo)向量增量 L(i)=norm(C(i,:); %變形后的長度 cs=C(i,:)/L(i); E=Section(N_Section,1); A=Section(N_Section,2); I=Section(N_Section,3); K_Element=beam2d_stiffness530(E,A,I,L(i),cs,Ele_F1(i,:); K_Global

43、=K_Global+Assemble(K_Element,Element(i,:),N_Node); %形成總剛k1 end Equation=K_Global,Val; %在殘余應(yīng)力下的位移求解 Disp_Transefer=ones(N_Node*3,1); Disp_Transefer(del,:)=0; Equation(del,:)=; Equation(:,del)=; n1=size(Equation,2); Dis2=-(Equation(:,1:n1-1)Equation(:,n1) %荷位移增量da1 zz=1; for i=1:N_Node*3; if Disp_Tran

44、sefer(i,1)=1; Disp_Transefer(i,1)=Dis2(zz,1); zz=zz+1; end end for i=1:N_Node Correc_dis(i,:)=Disp_Transefer(i*3-2:i*3,1); end Correc_dis1= Correc_dis1+ Correc_dis; Internal_F1=zeros(N_Node*3,1); %節(jié)點(diǎn)內(nèi)力向量 Update_Node1=Update_Node; Update_Node=Node+All_Disp(:,1:2)+Correc_dis1(:,1:2);%為了求切線剛度矩陣（改）a2下 i

45、f abs(Update_Node(11,2)<=1e-4&&abs(Update_Node(11,1)<=1e-4 break end for i=1:N_Element F_Global=zeros(N_Node*3,1); for j=1:2 ELEDISP1(j*3-2:j*3)=Correc_dis(Element(i,j),:); %取出當(dāng)前單元節(jié)點(diǎn)位移向量 end N1=Element(i,1); %i節(jié)點(diǎn)編號(hào) N2=Element(i,2); %j節(jié)點(diǎn)編號(hào) C1(i,:)=Update_Node1(N2,:)-Update_Node1(N1,:);

46、%a0下坐標(biāo)向量增量 L1(i)=norm(C1(i,:); %變形后的長度 cs1=C1(i,:)/L1(i); %桿件的cos和sin值 T1=cs1(1,1),cs1(1,2),0,0,0,0; -cs1(1,2),cs1(1,1),0,0,0,0; 0,0,1,0,0,0; 0,0,0,cs1(1,1),cs1(1,2),0; 0,0,0,-cs1(1,2),cs1(1,1),0; 0,0,0,0,0,1; ELEDISP(i,:)=T1* ELEDISP1(:); X1=L1(i)+ ELEDISP(i,4)- ELEDISP(i,1); Z1= ELEDISP(i,5)-ELEDI

47、SP(i,2); LN=(X12+Z12)0.5; sin1=Z1/LN; cos1=X1/LN; Citaloca=atan2(sin1,cos1); Ub=LN-L1(i); THRA=ELEDISP(i,3)- Citaloca; THRB=ELEDISP(i,6)- Citaloca; ELEDISP(i,:)=0,0,THRA,Ub,0,THRB; K_Element=beam2d_stiffness520(E,A,I,L1(i),cs1,Ele_F1(i,:);% L（i）為a1對(duì)應(yīng)下的 Ele_F(i,:)=K_Element*ELEDISP(i,:)' Ele_F1(i,:)= Ele_F1(i,:)+ Ele_F(i,:); C2(i,:)=Update_Node(N2,:)-Update_Node(N1,:); %a0下坐標(biāo)向量增量 L2(i)=norm(C2(i,:); %變形后的長度 cs2=C2(i,:)/L2(i); %桿件的cos和sin值 T2=cs2(1,1),cs2(1,2),0,0,0,0; -cs2(1,2),cs2(1,1),0,0,0,0; 0,0,1,0,0,0; 0,0,0,cs2(1,1),cs2(1,2),0; 0,0