abaqus2用戶單元子程序abaqus2用戶單元子程序abaqus2用戶單元子程序xxx公司abaqus2用戶單元子程序文件編號：文件日期：修訂次數(shù)：第1.0次更改批準審核制定方案設(shè)計，管理制度20ABAQUS用戶單元子程序(UEL)在這一章中將列舉兩個在這些年里發(fā)展過的ABAQUS/Standard用戶單元子程序（UEL）。第一個例子是一個非線性的索單元，我們的目的是通過這個比較簡單的例子讓讀者了解用戶單元子程序的基本開發(fā)過程；第二個例子是一個用于計算應(yīng)變梯度理論的單元，應(yīng)變梯度是當今比較熱點的一個科研前沿問題，有各種理論，我們?yōu)榱蓑炞C新的理論，需要數(shù)值結(jié)果與實驗對照來進行評價，整個例子的目的是通過它說明用戶子單元可以求解的問題范圍很廣，但是由于內(nèi)容比較艱深，程序也很長，所以這個例子我們并沒有給出最后的全部程序。另外，到目前為止，ABAQUS還只有隱式求解器ABAQUS/Standard支持用戶自定義單元，而顯式求解器ABAQUS/Explicit中還不支持這一功能。非線性索單元背景鋼索斜拉橋和斜拉索結(jié)構(gòu)廣泛應(yīng)用于土木工程建筑上。索力的計算分析是設(shè)計和施工的關(guān)鍵環(huán)節(jié)。清華大學工程力學系在采用ABAQUS進行荊沙長江斜拉橋的計算機仿真分析（這個項目我們已在第15章“ABAQUS在土木工程中的應(yīng)用（一）——荊州長江大橋南汊斜拉橋結(jié)構(gòu)三維仿真分析”中討論過）時，也曾進行了自行建立索單元的嘗試。本節(jié)介紹的就是這方面的工作。香港理工大學土木與結(jié)構(gòu)工程系采用ABAQUS有限元軟件進行計算，完成了香港TingKau斜拉橋和TsingMa懸索橋的結(jié)構(gòu)計算和分析。對于鋼索計算，他們采用梁單元進行模擬。由于梁單元含有彎曲剛度，計算的高階頻率值偏高，周期較低。一般假設(shè)索是單向受拉力的構(gòu)件。隨著應(yīng)變的非線性增加，索力呈非線性增加。盡管ABAQUS單元庫中有500個以上的單元類型，但是，還沒有索單元。本文發(fā)展了三維非線性索單元模型，形成ABAQUS的用戶單元子程序，可以利用ABAQUS輸入文件調(diào)入到具體的分析中。通過靜態(tài)和動態(tài)例題的計算比較，索單元工作良好?；竟皆谌S索單元計算中，如圖20-1所示，坐標x和位移u的變量表達式為： (x,y,z)(u,v,w) (20-1)應(yīng)變的公式為： (20-2)公式（20-2）中，L為索的長度，索的張力為： (20-3)在公式（20-3）中，A為截面面積，E為彈性模量，N0為初始張力。圖20圖20-1索單元在總體坐標系下，單元剛度矩陣為： (20-4)單元剛度矩陣中的子陣K分別由線性和非線性矩陣項組成： (20-5)在公式（20-5）中的KL和KNL均是3×3的對稱矩陣，分別為： 索單元的節(jié)點質(zhì)量為： (20-6)在公式（20-6）中，為密度。索單元的質(zhì)量矩陣為： (20-7)結(jié)構(gòu)的運動方程為： (20-8)公式（21-8）中Fext為作用在結(jié)構(gòu)上的外力。在不斷變化的索的變形中，求解運動方程，得到節(jié)點的位移值。應(yīng)用舉例圖19-2由五個單元組成的兩端鉸接的索桿結(jié)構(gòu)圖19-2由五個單元組成的兩端鉸接的索桿結(jié)構(gòu)由5個單元組成的兩端鉸接的索桿結(jié)構(gòu)，高5m長10m，6個節(jié)點號碼依次為101~106，如圖19-2所示。計算自由振動的頻率和周期。輸入文件中的用戶單元界面ABAQUS輸入文件（.inp）中的用戶單元界面如下：*HEADINGTwodimensionaloverheadhoistframeusing2nodesself-developedtrusselement,InitialforceNisdefinedinproperty(5)andreferencedbyuserelementSIUnits1-axishorizontal,2-axisvertical……*USERELEMENT,NODES=2,TYPE=U1,PROPERTIES=5,COORDINATES=3,VARIABLES=121,2,3*UELPROPERTY,ELSET=UTRUSS,,,7800,*ELEMENT,TYPE=U1,ELSET=UTRUSS……計算結(jié)果和比較表20-1列出了由用戶索單元計算的圖20-2所示結(jié)構(gòu)的固有周期，并與應(yīng)用ABAQUS梁單元B31的計算結(jié)果進行了比較。索單元與梁單元前4階模態(tài)的周期基本一致；索單元的第6~9階模態(tài)與梁單元第7~10階模態(tài)的周期基本一致。從第11階模態(tài)開始，隨著梁單元彎曲變形的增加，梁的彎曲剛度逐漸發(fā)揮作用并和軸向剛度耦合，與同階模態(tài)的索單元相比，梁單元的振動周期顯著降低，而頻率高于索單元。表20-1ABAQUS用戶索單元和梁單元B31計算的頻率比較振動模態(tài)索單元固有周期(Cycles/time)梁單元固有周期(Cycles/time)123456789101112用戶開發(fā)單元的缺點是不能采用ABAQUS的后處理進行顯示，只能從數(shù)據(jù)文件（.dat）中讀取結(jié)果。另外，ABAQUS的接觸算法等某些功能也無法應(yīng)用。非線性索單元用戶子程序subroutineuel(rhs,amatrx,svars,energy,ndofel,nrhs,nsvars,*props,nprops,coords,mcrd,nnode,u,du,v,a,jtype,time,dtime,*kstep,kinc,jelem,params,ndload,jdltyp,adlmag,predef,npredf,*lflags,mlvarx,ddlmag,mdload,pnewdt,jprops,njprop,period)CInclude''CAllcoordinatesinglobalCdimensionrhs(mlvarx,*),amatrx(ndofel,ndofel),*svars(12),energy(8),props(5),coords(mcrd,nnode),*u(ndofel),du(mlvarx,*),v(ndofel),a(ndofel),time(2),*params(3),jdltyp(mdload,*),adlmag(mdload,*),*ddlmag(mdload,*),predef(2,npredf,nnode),lflags(*),*jprops(*)Cdimensionsresid(6),uji(3),xji(3),smatrx(3,3)CCMaterialpropertiesarea=props(1)e=props(2)anu=props(3)rho=props(4)CInitialtensionforceinuserelementfn0=props(5)CCGeometry,stiffnessandmassparametersdx=coords(1,2)-coords(1,1)dy=coords(2,2)-coords(2,1)dz=coords(3,2)-coords(3,1)alen=sqrt(dx*dx+dy*dy+dz*dz)ang=atan(dy/dx)ak=area*e/alenam=*area*rho*alenCdoi=1,3uji(i)=u(i+3)-u(i)xji(i)=coords(i,2)-coords(i,1)enddostrain=(xji(1)*uji(1)+xji(2)*uji(2)+xji(3)*uji(3)*+*(uji(1)**2+uji(2)**2+uji(3)**2))/alentforce=e*area*strain+fn0eal=e*area/alen**3CCStiffnessmatrixparameterssmatrx(1,1)=eal*xji(1)**2+tforce/alensmatrx(1,2)=eal*xji(1)*xji(2)smatrx(1,3)=eal*xji(1)*xji(3)smatrx(2,1)=smatrx(1,2)smatrx(2,2)=eal*xji(2)**2+tforce/alensmatrx(2,3)=eal*xji(2)*xji(3)smatrx(3,1)=smatrx(1,3)smatrx(3,2)=smatrx(2,3)smatrx(3,3)=eal*xji(3)**2+tforce/alenCdo6k1=1,ndofelsresid(k1)=do2krhs=1,nrhsrhs(k1,krhs)=2continuedo4k2=1,ndofelamatrx(k2,k1)=4continue6continueCif(lflags(3).thenCNormalincrementationif(lflags(1).thenC*StaticCElementstiffnessmatrixdoi=1,3doj=1,3amatrx(i,j)=smatrx(i,j)amatrx(i,j+3)=-smatrx(i,j)amatrx(i+3,j)=-smatrx(i,j)amatrx(i+3,j+3)=smatrx(i,j)enddoenddoCCReactionforceif(lflags(4).thendoi=1,3force=ak*(u(i+3)-u(i))dforce=ak*(du(i+3,1)-du(i,1))sresid(i)=-dforcesresid(i+3)=dforcerhs(i,1)=rhs(i,1)-sresid(i)rhs(i+3,1)=rhs(i+3,1)-sresid(i+3)enddoelsedok=1,3force=ak*(u(k+3)-u(k))sresid(k)=-forcesresid(k+3)=forcerhs(k,1)=rhs(k,1)-sresid(k)rhs(k+3,1)=rhs(k+3,1)-sresid(k+3)enddoendifCelseif(lflags(1).thenC*Dynamicalpha=params(1)beta=params(2)gamma=params(3)Cdadu=(beta*dtime**2)dvdu=gamma/(beta*dtime)Cdo14k1=1,ndofelamatrx(k1,k1)=am*dadurhs(k1,1)=rhs(k1,1)-am*a(k1)14continuedoi=1,3doj=1,3amatrx(i,j)=amatrx(i,i)++alpha)*smatrx(i,j)amatrx(i+3,j+3)=amatrx(i+3,j+3)++alpha)*smatrx(i,j)amatrx(i,j+3)=amatrx(i,j+3)-+alpha)*smatrx(i,j)amatrx(i+3,j)=amatrx(i+3,j)-+alpha)*smatrx(i,j)enddoenddoCdoi=1,3force=ak*(u(i+3)-u(i))sresid(i)=-forcesresid(i+3)=forcerhs(i,1)=rhs(i,1)-(+alpha)*sresid(i)*-alpha*svars(i))rhs(i+3,1)=rhs(i+3,1)-(+alpha)*sresid(i+3)*-alpha*svars(i+3))enddoCdo16k1=1,ndofelsvars(k1+6)=svars(k1)svars(k1)=sresid(k1)16continueendifCelseif(lflags(3).thenCStiffnessmatrixdoi=1,3doj=1,3amatrx(i,j)=smatrx(i,j)amatrx(i,j+3)=-smatrx(i,j)amatrx(i+3,j)=-smatrx(i,j)amatrx(i+3,j+3)=smatrx(i,j)enddoenddoCelseif(lflags(3).thenCMassmatrixdo40k1=1,ndofelamatrx(k1,k1)=am40continueelseif(lflags(3).thenCCHalf-stepresidualcalculationalpha=params(1)doi=1,3force=ak*(u(i+3)-u(i))sresid(i)=-forcesresid(i+3)=forcerhs(i,1)=rhs(i,1)-am*a(i)-+alpha)*sresid(i)*+*alpha*(svars(i)+svars(i+6))rhs(i+3,1)=rhs(i+3,1)-am*a(i+3)-+alpha)*sresid(i+3)*+*alpha*(svars(i+3)+svars(i+9))enddoCelseif(lflags(3).thenCInitialaccelerationcalculationdo60k1=1,ndofelamatrx(k1,k1)=am60continuedoi=1,3force=ak*(u(i+3)-u(i))sresid(i)=-forcesresid(i+3)=forcerhs(i,1)=rhs(i,1)-sresid(i)rhs(i+3,1)=rhs(i+3,1)-sresid(i+3)enddoCdo62k1=1,ndofelsvars(k1)=sresid(k1)62continueCelseif(lflags(3).thenCOutputforperturbationsif(lflags(1).thenC*staticdoi=1,3force=ak*(u(i+3)-u(i))dforce=ak*(du(i+3,1)-du(i,1))sresid(i)=-dforcesresid(i+3)=dforcerhs(i,1)=rhs(i,1)-sresid(i)rhs(i+3,1)=rhs(i+3,1)-sresid(i+3)enddoCdokvar=1,nsvarssvars(kvar)=enddodoj=1,6svars(j)=rhs(j,1)enddoelseif(lflags(1).thenC*Frequencydo90krhs=1,nrhsdoi=1,3dforce=ak*(du(i+3,krhs)-du(i,krhs))sresid(i)=-dforcesresid(i+3)=dforcerhs(i,krhs)=rhs(i,krhs)-sresid(i)rhs(i+3,krhs)=rhs(i+3,krhs)-sresid(i+3)enddo90continuedokvar=1,nsvarssvars(kvar)=enddodoj=1,6svars(j)=rhs(j,1)enddoCendifendifCreturnend利用ABAQUS用戶單元計算應(yīng)變梯度塑性問題我們在本節(jié)內(nèi)容中主要介紹兩種應(yīng)變梯度理論，并在最后給出用這兩種應(yīng)變梯度理論編寫的用戶單元子程序的數(shù)值計算結(jié)果與實驗的對比。本節(jié)的目的在于讓讀者了解ABAQUS即使在面對如此復雜的理論問題時，也可以勝任。引言研究應(yīng)變梯度理論的意義很多試驗表明，當非均勻塑性變形特征長度在微米量級時，材料具有很強的尺度效應(yīng)。例如Fleck等在細銅絲的扭轉(zhuǎn)試驗中觀察到，當銅絲的直徑為12m時，無量綱的扭轉(zhuǎn)硬化將增加至170m直徑時的3倍；Stolken和Evans在薄梁彎曲試驗中也觀察到，當梁的厚度從100m減至m時，無量綱的彎曲硬化也顯著增加；而在單軸拉伸情況這種尺度效應(yīng)并不存在。在微米量級的尺度下，微觀硬度試驗與顆粒增強金屬基復合材料中也觀察到尺度效應(yīng)，當壓痕深度從10m減至1由于在傳統(tǒng)的塑性理論中，本構(gòu)模型不包含任何尺度，所以它不能預測尺度效應(yīng)。然而，在工程實踐中迫切需要處理微米量級的設(shè)計和制造問題，例如，厚度在1m或者更小尺寸下的薄膜；整個系統(tǒng)尺寸不超過10m的傳感器、執(zhí)行器和微電力系統(tǒng)（MEMS）；零部件尺寸小于10m的微電子封裝；顆粒或者纖維的尺寸在微米量級的先進復合材料及微加工等等?，F(xiàn)在的設(shè)計方法，如有限元方法（FEM）促使建立細觀尺寸下連續(xù)介質(zhì)理論的另一個目的是在韌性材料的宏觀斷裂行為和原子斷裂過程之間建立聯(lián)系。在一系列值得注意的試驗中，Elssner等測量了單晶鈮/藍寶石界面的宏觀斷裂韌度和原子分離功，使原子點陣或強界面分離所需要的力約為或者10Y（E為彈性模量，Y為拉伸屈服應(yīng)力）。而按照經(jīng)典的塑性理論，Hutchinson指出裂紋前方最大應(yīng)力水平只能達到4至5倍Y。很明顯這遠遠小于Elssner等在試驗中觀察到的結(jié)果，不足使原子分離?？紤]應(yīng)變梯度的影響有望解釋這一現(xiàn)象。應(yīng)變梯度理論簡介目前發(fā)展的應(yīng)變梯度理論有很多種，包括CS理論（偶應(yīng)力理論）、SG理論（拉伸和旋轉(zhuǎn)梯度理論）、MSG理論（基于細觀機制的應(yīng)變梯度塑性理論）以及TNT理論（基于Taylor關(guān)系的非局部應(yīng)變梯度理論）等。我們利用ABAQUS用戶單元主要進行了MSG和TNT兩種理論應(yīng)變梯度塑性的有限元分析，對于MSG塑性和TNT塑性都包括形變理論和流動理論，TNT塑性還包括了有限變形問題的形變和流動理論的分析。現(xiàn)在對MSG理論和TNT理論加以簡單的介紹。兩種應(yīng)變梯度理論基于細觀機制的MSG應(yīng)變梯度塑性理論基于位錯機制的MSG應(yīng)變梯度塑性理論是由位錯理論出發(fā)的，它通過一個多尺度、分層次的框架由微觀位錯機制推導出了宏細觀的應(yīng)變梯度塑性理論。這個理論相比于其它理論，物理機制更明確，構(gòu)造方法系統(tǒng)，而且第一次提出了材料長度的表達式。圖20-3是MSG應(yīng)變梯度塑性理論的原理圖，在微觀層次上塑性是由位錯運動產(chǎn)生的，在細觀層次上引入應(yīng)變的梯度與微觀的幾何必須位錯密度相關(guān)聯(lián)，通過細觀和微觀的功等效由微觀塑性推導出細觀的本構(gòu)理論。為了在細觀尺度下的應(yīng)變梯度塑性和微尺度下的Taylor硬化關(guān)系之間建立聯(lián)系，在MSG理論框架中采用如下的基本假設(shè)：圖圖20-3MSG理論中采用的多尺度框架 1）假設(shè)微尺度的流動應(yīng)力由位錯運動控制，并且遵守應(yīng)變梯度律給出的Taylor硬化關(guān)系 (20-9) 2）微觀尺度和細觀尺度的聯(lián)系是塑性功相等 (20-10) 3）在微尺度胞元中假設(shè)經(jīng)典塑性的基本結(jié)構(gòu)成立，其J2形變理論可以表示為 (20-11)其中，。微尺度的屈服條件為 (20-12) 基于以上的理論假設(shè)，應(yīng)變梯度塑性MSG形變理論本構(gòu)關(guān)系可以建立如下： (20-13) (20-14)其中 (20-15)式中： (20-16)為材料特征長度。 (20-17) (20-18)其中應(yīng)變梯度表示為 (20-19) (20-20)基于Taylor關(guān)系的非局部應(yīng)變梯度塑性理論（TNT理論）MSG應(yīng)變梯度理論物理機制明確、構(gòu)造方法系統(tǒng)，那么為什么還要發(fā)展TNT理論呢因為前面提到的應(yīng)變梯度塑性理論，無論CS、SG還是MSG，都是高階理論，引入了高階應(yīng)力和附加的邊界條件，而這些高階應(yīng)力和附加的邊界條件都難以測量，難以想象，因此無法得到工業(yè)界的認可，難以走向?qū)嵱谩＝?jīng)典的塑性問題控制方程是2階，而在MSG理論中由于高階應(yīng)力的影響，其控制方程是4階，這就增加了解決問題的復雜性。非局部連續(xù)介質(zhì)力學理論給我們以啟發(fā)，采用應(yīng)變的非局部加權(quán)積分來確定應(yīng)變的梯度，這樣就有了TNT——基于TAYLOR關(guān)系的非局部塑性理論。這種理論既保持MSG理論的所有優(yōu)點，同時與經(jīng)典的塑性理論相比又不增加方程的階數(shù)。但也正是由于TNT理論的非局部性質(zhì)，使其在求取解析解方面比較困難，也正因為如此，有限元解對TNT理論尤其重要。TNT作為非局部塑性理論，有三個基本特點：（1）流動應(yīng)力遵從TAYLOR硬化關(guān)系，這是TNT塑性理論的出發(fā)點。（2）應(yīng)變梯度和幾何必須位錯密度是非局部量，表示為應(yīng)變的加權(quán)平均。這是TNT塑性理論的核心概念。（3）TNT塑性理論保持經(jīng)典塑性理論的基本結(jié)構(gòu)。這個特點使得TNT塑性理論具有很強的實用潛力。和經(jīng)典的塑性理論相比，TNT塑性理論的特別之處在于屈服條件的不同：流動應(yīng)力不僅依賴于應(yīng)變，還同時依賴于應(yīng)變的梯度。這里應(yīng)變的梯度是非局部量。確定應(yīng)變梯度的非局部積分如下：將應(yīng)變在一點附近Taylor展開： (20-21)式中為以x為坐標原點的局部坐標。在包含x的表示體元內(nèi)積分上式： (20-22)假定特征尺寸足夠小，可忽略的高階小量，于是梯度可以表示為應(yīng)變的積分形式： (20-23)從而由關(guān)系式（20-19）、（21-20）可以得出應(yīng)變梯度的值。TNT形變理論的本構(gòu)關(guān)系： (20-24) (20-25)式中的為由（20-16）式給出的材料長度，為非局部變量，由（20-20）式給出。TNT流動理論的本構(gòu)關(guān)系： (20-26) (20-27)加載時，卸載時。ABAQUS用戶單元的使用如上文所述，由于MSG理論和TNT理論本身的復雜性，用它們做解析解比較困難。只是對于很少幾個簡單問題才有解析解。為比較理論與實驗以及解析解的符合程度，必須用有限元計算加以驗證。因為考慮應(yīng)變梯度的本構(gòu)關(guān)系與經(jīng)典理論完全不同，所以在有限元實現(xiàn)中不能利用現(xiàn)有的程序。但是由于ABAQUS具有的用戶子程序功能，為我們實現(xiàn)應(yīng)變梯度的有限元計算提供了很方便的條件。為了保證用戶子程序能夠完成計算應(yīng)變梯度問題的功能，編寫時必須遵守ABAQUS規(guī)定的法則。具體包括INCLUDE聲明、命名約定、重新定義變量、編譯和鏈接、測試和調(diào)試、用戶可用及不可用的通道號、中斷分析等。我們編寫的關(guān)于應(yīng)變梯度的用戶子程序中，用的是9節(jié)點矩形單元，由于節(jié)點變量包括應(yīng)變梯度和高階應(yīng)力的分量，故在編程過程中不能利用原有的幾何關(guān)系，也不能利用ABAQUS提供的前后處理程序，但可以方便地利用它的求解器。我們在用戶子程序中計算出有限元的剛度矩陣和右端項，再利用ABAQUS的求解器解線形代數(shù)方程，在用戶子程序UEL中需要保存的變量，如計算出的廣義應(yīng)變、廣義應(yīng)力等，都必須存在ABAQUS指定的變量SVARS中，下一次調(diào)用UEL時再從變量SVARS中讀取出需要的變量值。關(guān)于用戶子程序的調(diào)試可通過讀寫輸出文件來獲得信息，子程序中可以將調(diào)試信息寫在ABAQUS的消息文件（.msg）（通道號為7）或數(shù)據(jù)文件（.dat）（通道號為6）中，讀這些文件可得到調(diào)試信息，從而驗證程序是否已滿足了計算的要求。有限元計算的結(jié)果MSG理論的有限元計算結(jié)果我們編寫了應(yīng)變梯度本構(gòu)的平面問題和軸對稱問題的用戶子程序。利用軸對稱程序計算了微小尺度下的淺壓頭壓痕問題，并對微尺度壓痕的試驗結(jié)果進行了參數(shù)擬合；對于平面問題的程序，進行了靜態(tài)裂紋的分析，分析過程主要針對I型裂紋。圖20-4MSG理論擬合微壓痕試驗結(jié)果圖20-4是利用MSG理論計算微壓痕的曲線：坐標橫軸是壓痕深度的倒數(shù)，坐標縱軸是無量綱化的硬度的平方，其中H是微壓痕的硬度，H0是h取值很大時的壓痕硬度，它與壓痕的深度h無關(guān)。圖中的試驗點是McElhaney等人在1998年對多晶銅所作的結(jié)果?？梢?，利用MSG理論計算的結(jié)果曲線在從1/10個微米到幾個微米很大的范圍內(nèi)與試驗結(jié)果符合的非常好，計算擬合的材料特征長度l=也符合多晶銅由（20-16）式的計算結(jié)果以及試驗對銅的估算值，這表明MSG理論能夠比較準確地反映微米到亞微米量級材料的塑性行為。下面討論利用MSG本構(gòu)的平面問題程序計算I型靜止裂紋的結(jié)果。整個計算區(qū)域為1000l，遠場施加的是彈性位移邊界條件，外加K場強度為20。圖20-5給出了雙對數(shù)坐標下裂尖的奇異性曲線，即裂紋延長線上（實際上取的是最靠近裂紋面的一列高斯點，）的等效應(yīng)力對應(yīng)距裂尖距離的曲線。橫坐標是無量綱化的點到裂尖的距離，