【正文】 （5）提高非線性的程度，分別用常剛度法和切線剛度法計(jì)算，比較兩者的效率，得出相應(yīng)的結(jié)論。（3）采用常剛度算法的UMAT子程序計(jì)算，得出其Mises等效應(yīng)力值的范圍以及應(yīng)力應(yīng)變關(guān)系，然后同ABAQUS計(jì)算的相應(yīng)結(jié)果做比較，看兩者是否一致，驗(yàn)證常剛度法的正確性。6. 程序驗(yàn)證 任何有關(guān)材料模型的開(kāi)發(fā)都需要經(jīng)過(guò)實(shí)踐的檢驗(yàn)和驗(yàn)證，尤其對(duì)于工程問(wèn)題，通過(guò)驗(yàn)證可以發(fā)現(xiàn)程序中的BUG，并做相應(yīng)的修改，本文的驗(yàn)證主要采用綜合比較法，整個(gè)驗(yàn)證過(guò)程如下：（1）確定一個(gè)具體的模型和它的本構(gòu)關(guān)系。 對(duì)于程序沒(méi)有語(yǔ)法錯(cuò)誤卻有算法錯(cuò)誤，這種情況也是很常見(jiàn)的，對(duì)程序的算法做檢查，本文采用的主要方式是通過(guò)輸出變量到文本的方法，具體就是在程序內(nèi)用OPEN語(yǔ)句打開(kāi)一個(gè)文本文件，然后將程序中關(guān)心的中間變量用WRITE語(yǔ)句輸?shù)皆撐谋?，提交給ABAQUS以后，ABAQUS執(zhí)行到該處就會(huì)將其輸出到文本，這樣，只需要查看最后的文本文件，就可以發(fā)現(xiàn)算法中是否純?cè)阱e(cuò)誤，修正算法。39。39。39。. 常剛度法程序設(shè)計(jì)算法設(shè)計(jì)1：定義程序需要用到的常數(shù)和變量2：讀取ABAQUS定義的材料常數(shù)和狀態(tài)變量(這里只定義了一個(gè)狀態(tài)變量)，材料常數(shù)為，彈性模量E，泊松比Mu，屈服應(yīng)力Yield0，參數(shù)A,B,C,并且計(jì)算出剪切模量G,狀態(tài)變量為等效塑性應(yīng)變EQPLAS3：讀取應(yīng)力分量，計(jì)算平均應(yīng)力，應(yīng)力偏量以及Mises等效應(yīng)力平均應(yīng)力： 應(yīng)力偏量：Mises等效應(yīng)力：4：根據(jù)3計(jì)算的Mises等效應(yīng)力和2讀取的屈服應(yīng)力Yield0比較，如果Mises等效應(yīng)力小于屈服應(yīng)力，表明此時(shí)材料未屈服，那么轉(zhuǎn)到5，否則轉(zhuǎn)到65：雅可比矩陣,初始化為0，計(jì)算彈性矩陣，按照彈性理論更新應(yīng)力6：雅可比矩陣,初始化為01） ：計(jì)算切線模量H39。常剛度法采用的是：而切線剛度法采用的是：把握好了兩者的本質(zhì)上的區(qū)別，對(duì)于兩者的算法設(shè)計(jì)和程序開(kāi)發(fā)問(wèn)題便迎刃而解5. UMAT程序設(shè)計(jì)和編碼 本章將嚴(yán)格按照前一章推導(dǎo)的公式展開(kāi)程序設(shè)計(jì)和編碼，為了便于編程，本文將本構(gòu)關(guān)系做了抽象化處理，即將其描述成一個(gè)含參數(shù)的表達(dá)式，改變參數(shù)即可應(yīng)用于不同的模型，這樣做的好處是能保證程序的復(fù)用性，這也是本文反復(fù)強(qiáng)調(diào)的使用UMAT的原則。如果與相關(guān)，則方程的出現(xiàn)非一次項(xiàng)，變成非線性問(wèn)題，在實(shí)際工程中，特別是塑性成型問(wèn)題，材料的幾何方程，本構(gòu)方程以及邊界條件往往是非線性的也體現(xiàn)在中出現(xiàn)了，所以變?yōu)榱朔蔷€性問(wèn)題，要得到最基本的未知量，就必須求解非線性方程組1：直接迭代法又稱常剛度法[11]，在每次求解前，利用上次的解來(lái)求出這一次的值，然后利用和的倒數(shù)的乘積求出的當(dāng)前值 (46) 表達(dá)為迭代形式 (47) 上式可以看出，這種方法首先需要有一個(gè)初始的值，以便開(kāi)始迭代。. 非線性有限元算法理論 對(duì)于非線性問(wèn)題，在有限元求解該問(wèn)題時(shí)，對(duì)一個(gè)自由度總可以表達(dá)成，式中，為基本未矢量。當(dāng)偏離比例變形條件不多時(shí)，全量理論的計(jì)算結(jié)果和實(shí)險(xiǎn)結(jié)果比較接近。因此，又發(fā)展出塑性全量理論，即采用全量應(yīng)力和全量應(yīng)變表示塑性本構(gòu)關(guān)系的理論。用應(yīng)變?cè)隽啃问奖硎舅苄员緲?gòu)關(guān)系的理論稱為塑性增量理論。隨動(dòng)強(qiáng)化模型可以考慮材料的包興格（Bauschinger）效應(yīng)，在循環(huán)加載或可能出現(xiàn)反向屈服的問(wèn)題中，需要采用這種模型。在以上幾種強(qiáng)化模型中，各向同性強(qiáng)化模型應(yīng)用最為廣泛。(3) 混合強(qiáng)化把各向同性強(qiáng)化模型和隨動(dòng)強(qiáng)化模型加以組合，得到混合強(qiáng)化模型。如果用于卸載情況，它只適合反向屈服應(yīng)力等于應(yīng)力反轉(zhuǎn)點(diǎn)的材料，而通常材料不具備這種性質(zhì)，因此在塑性力學(xué)中還發(fā)展了其它強(qiáng)化準(zhǔn)則。而其形狀、中心及其在應(yīng)力空間的方位均保持不變[10]。 (45) 其中，是強(qiáng)化參數(shù)。Von Mises屈服準(zhǔn)則具體形式是，對(duì)于各項(xiàng)同性材料，應(yīng)力偏量第二不變量等于某一定值時(shí)，材料開(kāi)始進(jìn)入了塑性狀態(tài)。當(dāng)已知主應(yīng)力的大小次序時(shí)，使用Tresca屈服條件較為方便；若不知道主應(yīng)力的大小次序，則使用Von Mises屈服條件較為方便。對(duì)于巖土材料則常用Tresca屈服條件、DruckerPrager屈服條件和MohrCoulomb屈服條件。屈服條件是各應(yīng)力分量組合應(yīng)滿足的條件。（5）冪強(qiáng)化模型，為簡(jiǎn)化計(jì)算中的解析式，可將應(yīng)力應(yīng)變關(guān)系的解析式寫為σ＝σy（ε／εy）n，式中σy為屈服應(yīng)力，εy為與σy相對(duì)應(yīng)的應(yīng)變，n為材料常數(shù)。（3）理想剛塑性模型，用于彈性應(yīng)變比塑性應(yīng)變小得多且強(qiáng)化性質(zhì)不明顯的材料。1：本構(gòu)模型塑性力學(xué)的應(yīng)力應(yīng)變曲線通常有5種簡(jiǎn)化模型[8]：（1）理想彈塑性模型，用于低碳鋼或強(qiáng)化性質(zhì)不明顯的材料。所以必須確定二則之間的本構(gòu)關(guān)系，這種本構(gòu)關(guān)系可以用偏微分方程或者增量形式來(lái)描述。典型的彈塑性應(yīng)變?cè)谛遁d后要保持一個(gè)永久的變形。反之如果應(yīng)力達(dá)到后，再增加變形，應(yīng)力也必須增加，則材料是應(yīng)變硬化的，這時(shí)應(yīng)力是塑性應(yīng)變的函數(shù)，可解析為： (41)本構(gòu)關(guān)系反應(yīng)著應(yīng)力應(yīng)變之間的關(guān)系。如繼續(xù)加載，然后在卸載，材料始終保持永久的塑性變形。對(duì)大多數(shù)材料來(lái)說(shuō)，在單調(diào)加載的情況下，存在一個(gè)明顯的極限應(yīng)力，當(dāng)應(yīng)力低于時(shí)，材料保持線彈性。但是一經(jīng)卸載就可以看出兩者的區(qū)別。這是區(qū)別于非線性彈性材料的基本屬性。. 材料的彈塑性本構(gòu)關(guān)系彈塑性材料進(jìn)入塑性的特征是當(dāng)荷載卸去后存在不可恢復(fù)的永久變形。后來(lái)又提出超彈性(即具有彈性勢(shì)的)有限變形理論，由于理論方程的冗長(zhǎng)而復(fù)雜，且工程應(yīng)用也沒(méi)有提出這方面要求而被擱置。直至整個(gè)外力全部施加并得到滿足31的位移。然后回到32，如此循環(huán)，直至32右端為0，也即滿足31。根據(jù)32得到位移增量，此時(shí)要知道力對(duì)位移的梯度KT，以盡快找到滿足平衡條件的位移，由材料的Jacobian 矩陣和單元結(jié)合起來(lái)組裝得到(此處使用UMAT 提供的Jacobian 矩陣)。 以大變形、非性材料為例，整個(gè)計(jì)算步驟是這樣的： 整個(gè)外力不是一次加上，而是一點(diǎn)點(diǎn)加上的，不然會(huì)發(fā)散得不到結(jié)果的。但收斂速度無(wú)論慢快，31式才是判斷結(jié)果準(zhǔn)確與否的唯一標(biāo)準(zhǔn)。剛度矩陣其實(shí)就是力對(duì)位移的梯度。除了小變形、線彈性問(wèn)題，方程21是非性的，要用迭代的方法解出： （32） (33) i 表示一個(gè)增量步內(nèi)的第i 次迭代，n表示第n個(gè)增量步。 該矩陣只影響收斂速度，不影響計(jì)算結(jié)果的準(zhǔn)確性 （當(dāng)然，不收斂自然得不到結(jié)果）。. UMAT的使用方法 我們知道，有限元計(jì)算（增量方法）的基本問(wèn)題[7]是：已知第n 步的結(jié)果（應(yīng)力，應(yīng)變等），然后給出一個(gè)應(yīng)變?cè)隽?，?jì)算，UMAT要完成這一計(jì)算，并要計(jì)算DDSDDE(I,J)=。由于UMAT子程序在單元的積分點(diǎn)上調(diào)用，增量步開(kāi)始時(shí)，主程序路徑將通過(guò)UMAT的接口進(jìn)入U(xiǎn)MAT，單元當(dāng)前積分點(diǎn)必要變量的初始值將隨之傳遞給UMAT的相應(yīng)變量。NSHR：剪切應(yīng)力分量的個(gè)數(shù)。DTIME：增量步的時(shí)間增量。STRAN(NTENS)：應(yīng)變數(shù)組。SSE，SPD，SCD：分別定義每一增量步的彈性應(yīng)變能，塑性耗散和蠕變耗散。材料常數(shù)的個(gè)數(shù)，等于關(guān)鍵詞“*USER MATERIAL”中“CONSTANTS”常數(shù)設(shè)定的值。狀態(tài)變量矩陣的維數(shù)通過(guò)ABAQUS輸入文件中的關(guān)鍵詞“*DEPVAR”定義，關(guān)鍵詞下面數(shù)據(jù)行的數(shù)值即為狀態(tài)變量矩陣的維數(shù)。在增量步的結(jié)束必須更新?tīng)顟B(tài)變量矩陣中的數(shù)據(jù)。UMAT中應(yīng)力張量的度量為柯西應(yīng)力。在增量步的開(kāi)始，應(yīng)力張量矩陣中的數(shù)值通過(guò)UMAT和主程序之間的接口傳遞到UMAT中，在增量步的結(jié)束UMAT將對(duì)應(yīng)力張量矩陣更新。通常Jacobian矩陣是一個(gè)對(duì)稱矩陣，除非在“*USER MATERIAL”語(yǔ)句中加入了“UNSYMM”參數(shù)。各分量之間的順序根據(jù)單元自由度的不同有一些差異，所以編寫UMAT時(shí)要考慮到所使用單元的類別。由于主程序與UMAT之間存在數(shù)據(jù)傳遞，甚至共享一些變量，因此必須遵守有關(guān)UMAT的書寫格式，UMAT中常用的變量在文件開(kāi)頭予以定義，通常格式SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,1 RPL,DDSDDT,DRPLDE,DRPLDT,2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)INCLUDE‘’CHARACTER*80 CMNAMEDIMENSION STRESS(NTENS),STATEV(NSTATV),1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)user coding to define DDSDDE,STRESS,STATEV,SSE,SPD,SCDand,if necessary,RPL,DDSDDT,DRPLDE,DRPLDT,PNEWDTRETURNENDUMAT中的應(yīng)力矩陣、應(yīng)變矩陣以及矩陣DDSDDE、DDSDDT、DRPLDE等，都是直接分量存儲(chǔ)在前，剪切分量存儲(chǔ)在后。（2）幾乎可以用于力學(xué)行為分析的任何分析過(guò)程，幾乎可以把用戶材料屬性賦予ABAQUS中的任何單元。在輸入文件中，使用關(guān)鍵詞“*USER MATERIAL”表示定義用戶材料屬性。. UMAT接口的原理 用戶材料子程序（Userdefined Material Mechanical Behavior，簡(jiǎn)稱UMAT）是ABAQUS提供給用戶定義自己的材料屬性的Fortran程序接口[7][8]，使用戶能使用ABAQUS材料庫(kù)中沒(méi)有定義的材料模型。當(dāng)ABAQUS調(diào)用用戶子程序時(shí)，都會(huì)把當(dāng)前的STEP和INCREMENT利用用戶子程序的兩個(gè)實(shí)參KSTEP和KINC傳給用戶子程序，用戶可把它們輸出到外部文件中，這樣就可清楚的知道ABAQUS何時(shí)調(diào)用該用戶子程序。(2)當(dāng)用戶在用戶子程序中利用OPEN打開(kāi)外部文件時(shí)，要注意以下兩點(diǎn)：一是設(shè)備號(hào)的選擇是有限制的，只能取15～18和大于100的設(shè)備號(hào)，其余的都已被ABAQUS占用；二是用戶需提供外部文件的絕對(duì)路徑而不是相對(duì)路徑。ABAQUS應(yīng)用程序必須由用戶子程序調(diào)用。 在一個(gè)算例中，用戶可以用到多個(gè)用戶子程序。它的一般結(jié)構(gòu)形式是：SUBROUTINE S(x1,x2,……,xn)INCLUDE‘’（用于ABAQUS/Standard用戶子程序中）OR INCLUDE‘’)（用于ABAQUS/Explicit用戶子程序中）……RETURNENDx1，x2，……，xn是ABAQUS提供的用戶子程序的接口參數(shù)，有些參數(shù)是ABAQUS傳到用戶子程序中的，例如SUBROUTINE DLOAD中的KSTEP、KINC、COORDS，有些是需要用戶自己定義的，例如F，它們含有重要的參數(shù)，幫助ABAQUS主求解程序?qū)τ脩糇映绦蜻M(jìn)行編譯和鏈接。3. 用戶材料子程序UMAT. UMAT開(kāi)發(fā)環(huán)境設(shè)置 由于UMAT是采用FORTRAN語(yǔ)言編寫，那么要運(yùn)行UMAT就需要安裝FORTRAN的開(kāi)發(fā)環(huán)境， 同時(shí)還需要ABAQUS的支持，支持INTEL ，Intel Fortran安裝時(shí)又需要安裝Microsoft Visual Studio的相應(yīng)版本，經(jīng)過(guò)比較，+Intel +Microsoft VisualC++ 2005，相對(duì)于ABAQUS來(lái)說(shuō)，UMAT開(kāi)發(fā)環(huán)境的設(shè)置較為繁瑣，這給子程序的使用帶來(lái)諸多不便，為了解決這一問(wèn)題，我用C語(yǔ)言編制了ABAQUS子程序編譯環(huán)境設(shè)置工具，只需要將安裝文件解壓到ABAQUS的安裝目錄，運(yùn)行安裝程序就可以了，整個(gè)過(guò)程不需要人工干預(yù)，也不需要安裝龐大的VisualC++ 2005，如圖31所示 圖31. UMAT注意事項(xiàng) ABAQUS的用戶子程序是根據(jù)ABAQUS提供的相應(yīng)接口，按照Fortran語(yǔ)法，用戶自己編寫的代碼。以上三種基本結(jié)構(gòu)，是組成結(jié)構(gòu)化程序的基本結(jié)構(gòu)形式。(3)循環(huán)結(jié)構(gòu)循環(huán)結(jié)構(gòu)也稱重復(fù)處理結(jié)構(gòu)，即重復(fù)執(zhí)行某一功能塊，直到滿足（或不滿足）某一條件為止。(2)分支選擇結(jié)構(gòu)在給定的條件下，分支選擇結(jié)構(gòu)判斷選擇哪一條路徑執(zhí)行，不同路徑完成的功能是不同的。它有一個(gè)入口和一個(gè)出口，并在入口和出口之間包含著若干個(gè)功能塊，其中每一個(gè)功能塊可以是一個(gè)非轉(zhuǎn)移語(yǔ)句。 FORTRAN語(yǔ)言問(wèn)世以來(lái)，根據(jù)需要幾經(jīng)發(fā)展，先后推出了不同的版本，主要版本有FORTRAN 77，F(xiàn)ORTRAN 90，F(xiàn)ORTRAN 95，A