基于C語言的鉆柱有限元單元剛度矩陣計算_圖文

1、2010年第39卷 石油礦場機(jī)械第3期第36貞 OIL FIEI。D EQUIPMENT 2010.39(3:3639文章編號:100卜3482(201003003604基于C語言的鉆柱有限元單元剛度矩陣計算秦 偉(中國石油大學(xué)(華東機(jī)電工程學(xué)院,山東東營257061摘要:單元剛度矩陣的建立是有限元分析中的一個重要步驟。就實際鉆柱受力與變形分析中出現(xiàn) 的20節(jié)點(diǎn)六面體單元進(jìn)行研究,給出了計算單元剛度矩陣的方法和步驟,并運(yùn)用C語言編寫了相 應(yīng)程序,為解決大型三維有限元問題提供了較好的思路。關(guān)鍵詞:鉆柱;單元剛度矩陣;有限元;20節(jié)點(diǎn);C語言中圖分類號:TE921.2文獻(xiàn)標(biāo)識碼:ACalculat

2、ion of Eiement Stiffness Matrix about FEA of DrillString Using C Programming LanguageQIN Wei(College of Mechanical and Electronic Engineering,China University of Petroleum,Dongying 257061,ChinaAbstract:The calculation of element stiffness matrix is an important step of Finite Element analy sis.A stu

3、dy was made on 20-Node and 6一Side element which appeared on actual analysis of de formation and stress of drill string and put on the detailed methods and procedures of the caleula tion of element stiffness matrix.Moreover,the corresponding procedures using C Programming Language was prepared.All pr

4、ovide a better idea for the settlement of largescale threedimen sionaI finite element.Key words:drill string;element stiffness matrix;finite element analysis;20一node;C programming language有限元法是一種求解復(fù)雜工程結(jié)構(gòu)非常有效的 數(shù)值方法,是將所研究的工程系統(tǒng)轉(zhuǎn)化成一個結(jié)構(gòu) 近似的有限元系統(tǒng),該有限元系統(tǒng)由節(jié)點(diǎn)及單元組 合而成,以取代原有的工程系統(tǒng)1。有限元系統(tǒng)可 以轉(zhuǎn)化成一個數(shù)學(xué)模式,并根據(jù)該數(shù)學(xué)模式,進(jìn)而

5、得 到該有限元系統(tǒng)的解答,并通過節(jié)點(diǎn)、單元表現(xiàn)出 來。單元剛度矩陣是有限元分析中的第1步,也是 后面形成總體剛度矩陣和有限元方程組求解的基 礎(chǔ),但是目前尚無20節(jié)點(diǎn)六面體單元剛度矩陣的 相關(guān)文獻(xiàn)。問題的提出鉆柱是石油鉆井設(shè)備的主要組成部分,其主要 功能是將井口處的動力傳遞到井底處的鉆頭以破碎 巖石,并通過鉆井液將巖屑帶到地面。鉆柱受到井 壁的限制,承受自重和井口處的大鉤載荷、扭矩及井 底鉆頭處巖石的反作用力等。采用20節(jié)點(diǎn)六面體 單元對鉆柱的受力和變形進(jìn)行有限元分析3罐。 1.1鉆柱受力復(fù)雜在鉆進(jìn)過程中,鉆柱除了承受上述作用力外,還收稿日期:20090908基金項目:山東省自然科學(xué)基金項目“細(xì)

6、長柔性體動力學(xué)及其計算機(jī)仿真研究”(Q2006A08作者簡介:秦偉(1985一,男,山東鄒城人,碩士研究生主要從事鉆柱有限元分析和C語言編程實現(xiàn)方面的研究,E mail:qwl20060307163.tom。第39卷第3期秦 偉:基于C語青的鉆柱有限元單元剛度矩陣計算=要承受鉆柱繞井眼軸線公轉(zhuǎn)的離心力及環(huán)空鉆井液 柱的外擠壓力,另外還會發(fā)生各種振動,例如縱向振 動、扭轉(zhuǎn)振動和橫向擺振等3。在有限元計算中,受 力越復(fù)雜,單元選擇也越復(fù)雜,通過上述分析可知 一、二維單元無法解決這一復(fù)雜問題,需要選擇三維 單元進(jìn)行分析,20節(jié)點(diǎn)單元陽1叩是其中之一。 1.2鉆柱變形復(fù)雜正是由于鉆柱的受力過于復(fù)雜,鉆

7、柱由原來的 關(guān)于軸線對稱的空心圓柱變?yōu)殡S井眼發(fā)生彎曲并且 關(guān)于井眼軸線不對稱的復(fù)雜形狀,因此,把它看成是 軸對稱結(jié)構(gòu)或者空間梁結(jié)構(gòu)都是不合適的,而只能 看成三維實體結(jié)構(gòu)。由于鉆柱自身的曲面結(jié)構(gòu)與 20節(jié)點(diǎn)六面體正好相符合,因此應(yīng)選擇20節(jié)點(diǎn)六面體單元1¨。但是文獻(xiàn)1和文獻(xiàn)3都沒有給出詳細(xì)的20節(jié)點(diǎn)六面體單元剛度矩陣的求解方法和步驟,部分 學(xué)者都繞過這一難點(diǎn),將鉆柱1171簡化為空間梁或 者把它當(dāng)作軸對稱結(jié)構(gòu)來解決,結(jié)果產(chǎn)生了較大誤 差,因此研究該單元剛度矩陣的建立方法是必要的。2單元剛度矩陣的求解20節(jié)點(diǎn)六面體單元采用等參變換,即單元幾何 形狀的變換和單元內(nèi)的場函數(shù)采用相同數(shù)目的節(jié)點(diǎn)

8、B6×60=OOa_一azOa3y a3x。2矸EP(1(-卜ak2.“_rNl00N10LNl1笆1一 笆一 1一戶OO O 笆 1“ 1皂 1岸參數(shù)及相同的插值函數(shù)進(jìn)行變換。其位移模式及坐 標(biāo)變換公式為“=N。口l,。=N,u。,叫=Nf硼f(1z=N。Xi,y=NiYi,z=Nizf(2其中,形函數(shù)為Nt5言(1+島(1+珈(1+如(島+咖+如一2,(i一1,3,5,7,13,15,17,19Nz5寺(1一孑(1+伽(1+知,(i=2,6,14,18 Nr2÷(1一礦(1+島(1+島,(i一4,8,16,20 Nt一寺(1一P(1-+-毒o(1+伽,(i一9,10,1

9、1,12島2車,伽一孕17,函3§善式中,、孕、繭為等參六面體單元中對應(yīng)點(diǎn)的坐標(biāo)值。單元剛度矩陣總的求解公式為K。=J.二_。-r:。:。B7DBJded7d箏式中,B為幾何矩陣;D為彈性矩陣(常數(shù)陣;J為雅可比矩陣。笆1一“笆1一肛 1OO 0O OO OOOOO熱02(1一“。 蠆1(卜-2盧豈_。蠆1(1-一2笆盧_、l oo o ooo一 o oMo oMooo 旦如。旦知旦如oa一如ooa一毋oa一弛石油礦場機(jī)械 2010年3月式中,E為彈性模量,MPa;/。為泊松比。J2麓努一§8Nf 皇蕊疊§aNf 蘭而置雖aN, 毫蕊3i§aN。魚瓦y

10、z§aNi暑而孓§aNi直瓦yt§aNi 三蕊§aNf 喜百翻§aNf 宣可將上述形狀函數(shù)和對應(yīng)的坐標(biāo)值代入求解,而 后進(jìn)行三重積分就可求解出對應(yīng)的單元剛度矩陣。3編程求解的難點(diǎn)3.1雅可比矩陣的求解該工作計算量龐大,首先要求解各個形函數(shù)Nl 對手、'7、善求偏導(dǎo),而后與20點(diǎn)坐標(biāo)矩陣相乘(60X 3,得到3×3的矩陣,這一矩陣是隨手、叩、改變 的。除此之外,還要求解求J、I J l,因為對應(yīng)的曰 矩陣要求形函數(shù)對z、Y、z的偏導(dǎo),這需要一個關(guān) 系,即盟8x,坐8y,塑8z=J-1等,嚳,等1 L .J J L aea'

11、;7a善J 、u 7 3.2數(shù)值積分的求解高斯積分法是解決這類上下限為±1的積分的 有效方法。在此積分方案中,積分點(diǎn)不是等間距分 布,積分點(diǎn)的位置由下述方法確定。首先定義療次 多項式P(亭,即P(搴=(e一6(e一邑(亭一己=II拿一8(4 而后根據(jù)I亭P(d￥-0(i=0,1,行一1來確定 ,1個積分點(diǎn)的位置。對于三維數(shù)值積分,有I 2J-1j-_1J一。F(8吩fdSdr/d#2HiH,H。F(8,孕,蚤 (5 式中,H。、HJ、H。為積分權(quán)系數(shù);毫、孕、蚤為積分點(diǎn) 坐標(biāo)。高斯積分法的積分階次的選擇按照公式竹=夕一,l+1來確定,P是插值函數(shù)中完全多項式的方 次(此處為4,rn為

12、微分算子L中導(dǎo)數(shù)的階次(此處 為1,因此20節(jié)點(diǎn)單元應(yīng)采取4×4×4高斯積分 來進(jìn)行單元剛度矩陣的計算。對于四階高斯積 分m1,所用積分點(diǎn)坐標(biāo)和權(quán)函數(shù)如表1所示。表1四階高斯積分參數(shù)用計算機(jī)編程時,需要采用三重循環(huán),才能將 H,、H,、H。求解完整,即它對應(yīng)的是i、.、m在O4變化的乘積。4算例運(yùn)用上述求法,進(jìn)一步對結(jié)構(gòu)位移進(jìn)行求解,研 究了東營某化工廠一受內(nèi)壓的旋轉(zhuǎn)厚壁筒(與鉆柱 類似。該厚壁筒長為12cm,內(nèi)徑10cm,外徑20 cm,兩端自由,內(nèi)徑受壓p=1.2×108Pa,并以角速 度叫=209rad/s繞中心軸轉(zhuǎn)動。材料的彈性模量 E一2X1011N/m

13、2,泊松比=0.3,容重y一0.7×105N/m3。取中央橫截面為xy坐標(biāo)平面,由對稱 性取1/8采用20節(jié)點(diǎn)六面體等參單元分析,單元劃 分如圖1所示,共4個單元51個節(jié)點(diǎn),采用4×4×4高斯積分,形成的單元剛度矩陣為60×60,采用C 語言編程求解。由有限元法計算出的節(jié)點(diǎn)1、2、3的 徑向位移分別為“1=0.03934ITlm、u2=0.04436 mm、U3=0.05810mm,而精確值為Ul=0.04000 mm、u2=0.04517mm、u30.05810mm,相對誤 差均在2%以下。圖1厚壁圓筒的單元和節(jié)點(diǎn)第39卷第3期 秦 偉:基于(、語胄的

14、鉆枉有限元單元剛度矩陣計算 39皇 掛論 分析J石油礦場機(jī)械,2008,37(10:5658.J ;嗣叱 16趙國山,管志川,王以法,等.鉆柱結(jié)構(gòu)聲傳輸特性試 1 目前國內(nèi)對三維有限元的研究已經(jīng)取得了 驗研究J.石濉礦場機(jī)械,2009,38(11:4549 :定竺量:璺麓登莖苧笆毫燃簍 口73呈數(shù)霽囂蓍勰j莩。囂鬻茄裟鬟要羹 糊。本文提出的方法簡單、可行,具有很強(qiáng)的可操 i_:?！癏1”“”作性。E18李春明.巨大節(jié)點(diǎn)數(shù)有限元的數(shù)值計算要點(diǎn)J/OL.2 20節(jié)點(diǎn)三維塊單元剮度矩陣的建立不僅中國科技論文在線,2009一lo13.對鉆柱的受力分析有用。還可運(yùn)用于其他三維有限元問題當(dāng)中。3 算例表明

15、,該求解方法簡單實用,誤差較小。參考文獻(xiàn):1署橫木,全興華.工程有限單元法M.東營:石油大 學(xué)出版社,2003.2趙國珍,龔偉安.鉆井力學(xué)基礎(chǔ)M.北京;石油工業(yè) 出版社。1988:45103.3鄧學(xué)峰.水平井套管受力變形分析的有限元方法J. 石油礦場機(jī)械,2009,38(5:15-18.43王素玲,葉鑫銳,秦慶忠.注水井反沖洗套管承載能力 有限元分析J.石油礦場機(jī)械。2009,38(1:3134. 5張瑩,吳澤兵。王嘉容.三牙輪鉆頭滑動軸承副接觸 有限元分析J.石油礦場機(jī)械。2009,38(2:5962. 6戴相富.鉆修機(jī)絞車滾筒的有限元分析33.石油礦場 機(jī)械,2009,38(3:4648.

16、7張水澤,梁政,蔣發(fā)光,等.復(fù)雜結(jié)構(gòu)有限元分析強(qiáng)度 判定方法J.石油礦場機(jī)械,2009,38(5:5-8.8李鵬,牟新明。高明,等.頂驅(qū)中心管強(qiáng)度有限元分 析J.石油礦場機(jī)械,2009,38(9:5962.9趙先瓊,楊曉紅.ANSYS有限元分析與20節(jié)點(diǎn)塊單 元J.岳陽師范學(xué)院學(xué)報:自然科學(xué)版,2001,14(3: 2931.10方鵬。江進(jìn)國.有限元法及CAE技術(shù)在鉆柱研究 中的應(yīng)用J.裝備制造技術(shù),2006(5:3940.11楊東全.鉆柱靜態(tài)變形的有限元分析J.暨南大學(xué) 學(xué)報,2005,26(1:137140.12尹曉麗,牛文杰,張中慧.鉆桿自動傳送系統(tǒng)及設(shè)計 方案J.石油礦場機(jī)械。2009

17、,38(7:4246.133湯云霞,王立軍,王文云.鉆鋌連接螺紋斷裂失效分 析及結(jié)構(gòu)優(yōu)化J.石油礦場機(jī)械。2009,38(6:25 28.14崔曉華,周學(xué)芹,李玉海.井下鉆柱縱向振動特性分析 J.石油礦場機(jī)械。2009,38(11:1215.151朱煥剛,張 寶.基于ANSYS軟件的鉆柱縱向振動 附:部分源程序(C語言char Input(char*instr;讀人前處理信息 void Calculate(; 有限元法求解部分主程序 float Asek(int n,int iask;求解單元剛度矩陣 void Astk(;組裝總剛度矩陣void Astp(;組裝總載矩陣void Cons(;

18、處理約束支承條件void Oned(;一維變帶寬存儲函數(shù)void Solve(;求解方程組void Result(;輸出結(jié)果constdouble H4=0.3478548451, 0.6521451549,0.6521451549,0.347854845 1;const double z4=0.8611363116, 0.3399810436,一0.3399810436,一0.86l 136 3116;for(int i=O;i<4;i+for(int v=0;v<4;v+for(int w=0;w<4;w+F(z(i,z(v,z(w,n;for(int i5=0;i5&l

19、t;60;i5+for(int i6=0;i6<60;i6+EKi5i6+=H(i*H(v*H(w *KEi5i6;/采用高斯積分 基于C語言的鉆柱有限元單元剛度矩陣計算作者:秦偉 , QIN Wei作者單位:中國石油大學(xué)(華東機(jī)電工程學(xué)院,山東,東營,257061刊名:石油礦場機(jī)械英文刊名:OIL FIELD EQUIPMENT年,卷(期:2010,39(3參考文獻(xiàn)(18條1. 趙先瓊;楊曉紅 ANSYS有限元分析與20節(jié)點(diǎn)塊單元 期刊論文-岳陽師范學(xué)院學(xué)報(自然科學(xué)版 2001(032. 李鵬;牟新明;高明 頂驅(qū)中心管強(qiáng)度有限元分析 期刊論文-石油礦場機(jī)械 2009(093. 張永澤;梁政;蔣發(fā)光 復(fù)雜結(jié)構(gòu)有限元分析強(qiáng)度判定方法 期刊論文-石油礦場機(jī)械 2009(054. 李春明 巨大