ABAQUS-Fortran二次開發(fā).doc

《ABAQUS-Fortran二次開發(fā).doc》由會員分享，可在線閱讀，更多相關《ABAQUS-Fortran二次開發(fā).doc（73頁珍藏版）》請在裝配圖網(wǎng)上搜索。

目 錄 摘 要 I ABSTRACT II 1. 緒論 1 1.1. 課題的研究背景 1 1.2. 本文的研究內(nèi)容和方法 2 2. 基于ABAQUS軟件的二次開發(fā) 3 2.1. ABAQUS介紹 3 2.2. ABAQUS各模塊簡介 3 2.3. ABAQUS的二次開發(fā)平臺 5 2.4. ABAQUS的二次開發(fā)語言 6 3. 用戶材料子程序UMAT 8 3.1. UMAT開發(fā)環(huán)境設置 8 3.2. UMAT注意事項 9 3.3. UMAT接口的原理 10 3.4. UMAT的使用方法 12 4. 材料非線性問題 14 4.1. 材料的彈塑性本構關系 14 4.2. 非線性有限元算法理論 17 4.3. 增量理論常剛度法公式推導 20 4.4. 增量理論切線剛度法公式推導 21 5. UMAT程序設計和編碼 25 5.1. 本構關系描述 25 5.2. 常剛度法程序設計 27 5.3. 常剛度法程序編碼 29 5.4. 切線剛度法程序設計 32 5.5. 切線剛度法程序編碼 35 5.6. 程序的調(diào)試 38 6. 程序驗證 40 6.1. 問題描述 40 6.2. 本構關系 41 6.3. ABAQUS自帶材料模型計算 41 6.4. 常剛度法的UMAT驗證 43 6.5. 切線剛度法的UMAT驗證 45 6.6. 兩種算法的比較分析 47 7. 結論與展望 51 7.1. 結論 51 7.2. 展望 51 致 謝 53 參考文獻 54 附1：ABAQUS自帶彈塑性材料驗證的INP文件 55 附2：用于算法驗證的INP文件 61 摘 要 ABAQUS軟件功能強大，特別是能夠模擬復雜的非線性問題，它包括了多種材料本構關系及失效準則模型，并具有良好的開放性，提供了若干個用戶子程序接口，允許用戶以代碼的形式來擴展主程序的功能。 本文主要研究了ABAQUS用戶子程序UMAT的開發(fā)方法，采用FORTRAN語言編制了各向同性硬化材料模型的接口程序，研究該類材料的彈塑性本構關系極其實現(xiàn)方法。 本文緊緊圍繞UMAT的二次開發(fā)技術，首先對其接口原理做了詳細介紹，然后針對非線性有限元增量理論中的常剛度法和切線剛度法的算法理論做了深入的剖析，推導出了常剛度法和切線剛度法的算法理論的具體表達式，然后分別編制了兩種算法的UMAT程序，最后建立了一個具體的驗算模型，通過與ABAQUS自帶彈塑性本構關系的計算結果相比較，驗證兩者的正確性。 本文還對常剛度法和切線剛度法得算法效率做了對比，得出了在非線性程度較高時切線剛度法效率高于常剛度法的結論。 關鍵字： ABAQUS、UMAT、有限元、材料非線性、FORTRAN、切線剛度 ABSTRACT ABAQUS software powerful, especially to simulate complex non-linear problem, which includes a wide range of material constitutive model and failure criteria, and has a good open, providing a number of user subroutine interface that allows users to code form to expand the functions of the main program. This paper studies the user subroutine UMAT of ABAQUS development methods, the use of FORTRAN language isotropic hardening material model of the interface program, studied the effects of such material is extremely elastic-plastic constitutive relation method. This article UMAT tightly around the secondary development of technology, the first principle of its interface detail, and then for the theory of nonlinear finite element incremental stiffness of the regular tangent stiffness method and the theory of algorithms to do an in-depth analysis of deduced a regular tangent stiffness and rigidity of the law of the specific expression of algorithm theory, and then the preparation of the two algorithms, respectively, of the UMAT program, and finally the establishment of a specific model checking, bringing with ABAQUS elasto-plastic constitutive relation of the calculated results compared to verify the correctness of the two. This article also often stiffness and tangent stiffness method was to do a comparison of algorithm efficiency is obtained when a higher degree in the non-linear tangent stiffness method more efficient than the conclusions of law often stiffness. KEY WORDS：ABAQUS、UMAT、Finite element、Material nonlinearity、FORTRAN、Tangent stiffness 1. 緒論 1.1. 課題的研究背景 有限單元法基本思想的提出，可以追溯到克勞夫(R.W.Clough)在1943年的工作[1]，他第一次嘗試應用定義在三角形區(qū)域上的分片連續(xù)函數(shù)和最小位能原理相結合，來求解St. Venant扭轉(zhuǎn)問題。1960年克勞夫進一步處理了平面彈性問題，并第一次提出了“有限單元法”的名稱，使人們開始認識了有限單元法的功效。 四十多年來，隨著電子計算機的廣泛應用和發(fā)展，有限單元法的理論和應用都得到迅速的，持續(xù)不斷的發(fā)展，其應用己由彈性力學平面問題擴展到空間問題、板殼問題，由靜力學問題擴展到穩(wěn)定問題、動力問題和波動問題。分析的對象從彈性材料擴展到塑性、粘彈性、粘塑性和復合材料等，從固體力學擴展到流體力學、傳熱學等連續(xù)介質(zhì)力學領域。在工程分析中的作用已從分析和校核擴展到優(yōu)化設計并和計算機輔助設計。 利用有限元軟件解決工程和科學問題，是有限元理論應用于工程設計和科學研究實踐的主要形式。由于工程設計的巨大市場需要，有限元軟件的發(fā)展是很迅速的，目前常用的大型有限元軟件常見的有Sap2000，ADINA，MSC/NASTRAN，MSC Marc，ANSYS，ABAQUS等，這些軟件的共同特點是具有豐富的單元庫和求解器，強大而可靠的分析功能，人們利用這些軟件解決了很多工程建設和工業(yè)產(chǎn)品設計中遇到的問題，取得了巨大的經(jīng)濟技術效益。 由于工程問題的千差萬別，不同的用戶有不同的專業(yè)背景和發(fā)展方向，通用軟件不免在具體的專業(yè)方面有所欠缺，針對這些不足，大部分的通用軟件都提供了二次開發(fā)功能，以幫助用戶減少重復性的編程工作、提高開發(fā)起點、縮短研發(fā)周期、降低開發(fā)成本，并能簡化后期維護工作，給用戶帶來很多方便。基于通用軟件平臺進行開發(fā)，是目前研究的一個重要發(fā)展方向。 ABAQUS也提供了若干用戶子程序（User Subroutines）接口，它是一個功能非 常強大且適用的分析工具，與命令行的程序格式相比，用戶子程序的限制少得多， 從而使用更加靈活方便。針對ABAQUS所提供的本構關系模型種類有限，無法滿足 工程應用需要的問題，用戶子程序中的用戶材料子程序（User-defined Materia Mechanical Behavior，簡稱UMAT）接口可以幫助用戶定義自己的材料本構模型和算 法，這是ABAQUS的獨到之處。由于其操作方便，能被靈活地應用于各個領域中， 尤其受到用戶的青睞。 1.2. 本文的研究內(nèi)容和方法 ABAQUS中用戶材料子程序UMAT的開發(fā)主要解決兩方面的問題：本構模型的建立和積分算法的選擇。 本文主要研究非線性材料的UMAT實現(xiàn)方法，并重點研究其迭代算法部分，目前，用戶材料子程序UMAT的迭代算法主要是常剛度法，常剛度法的優(yōu)點在于算法原理較簡單，程序編寫較方便，缺點是當遇到復雜非線性材料時，其迭代次數(shù)較多，收斂速度也較慢，在這個情況下，本文采取的是一種迭代次數(shù)較少且收斂速度較快的切線剛度法，具體就是采用FORTRAN語言編制了基于Von-Mises模型的接口程序,并采用切線剛度算法，通過與ABAQUS自帶本構關系計算的結果相比較，驗證其正確性。 本文的研究工作緊緊圍繞UMAT的二次開發(fā)技術，首先根據(jù)有限元方法推導材料非線性問題算法的公式，然后參考UMAT接口規(guī)范設計程序的算法流程，繼而編寫出該程序，最后建立一個具體的本構和具體的模型做測試，驗證程序的正確性，在這一過程中，調(diào)試是一個非常重要的過程，占用了大量的時間，在調(diào)試程序時采用了將中間變量輸出到文本的方式，這樣能明確跟進迭代過程，發(fā)現(xiàn)算法或程序的缺陷。 本文采用的本構關系是經(jīng)過歸納和抽象的，也就是說本文的程序并不僅僅是只針對某個具體模型和問題，而是針對所有符合抽象出的各向同性硬化材料，這樣做的好處是能保證程序的通用性和復用性，避免以后的重復勞動，當然，這也是符合ABAQUS軟件設計UMAT接口的宗旨的。 2. 基于ABAQUS軟件的二次開發(fā) 2.1. ABAQUS介紹 ABAQUS是一套功能強大的基于有限元法的工程模擬軟件[2]，其解決問題的 范圍從相對簡單的線性分析到最富有挑戰(zhàn)性的非線性模擬問題。ABAQUS具備十分豐富的、可模擬任意實際形狀的單元庫。并與之對應擁有各種類型的材料模型庫，可以模擬大多數(shù)典型工程材料的性能，其中包括金屬、橡膠、高分子材料、復合材料、鋼筋混凝土、可壓縮彈性的泡沫材料以及巖石和土這樣的地質(zhì)材料。 作為通用的模擬分析工具，ABAQUS 不僅能解決結構分析中的問題，還能模擬和研究各種領域中的問題，如熱傳導、質(zhì)量擴散、電子元器件的熱控制（熱一電耦合分析）、聲學分析、土壤力學分析（滲流——應力耦合分析）和壓電介質(zhì)力學分析。 ABAQUS為用戶提供了廣泛的功能，且使用起來又十分簡明。最復雜的問題也可以很容易地建立模型[3]。例如復雜的多部件問題可以通過對每個部件定義材料模型和幾何形狀，然后再把它們組裝起來而構成。在大部分模擬分析問題中，甚至在高度非線性問題中，用戶也只需要提供結構的幾何形狀、材料性能、邊界條件和荷載工況這樣的工程數(shù)據(jù)就可以進行分析。在非線性分析中，ABAQUS能自動選擇合適的荷載增量和收斂精度。不僅能選擇這些參數(shù)值，而且能在分析過程中不斷地調(diào)整參數(shù)來保證有效地得到高精度的解，很少需用戶去定義這些參數(shù)。 2.2. ABAQUS各模塊簡介 ABAQUS 有兩個主要的分析模塊：ABAQUS/Standard 和ABAQUS/Explicit 。ABAQUS/Standard還有兩個特殊用途的附加分析模塊：ABAQUS/Aqua和ABAQUS/Design。另外，還有ABAQUS 分別與ADAMS/Flex，C-MOLD和Mold flow的接口模塊：ABAQUS/ADAMS，ABAQUS/C-MOLD和ABAQUS/ MOLDFLOW。ABAQUS/CAE是完全的ABAQUS工作環(huán)境模塊，它包括了ABAQUS模型的構造，交互式提交作業(yè)、監(jiān)控作業(yè)過程以及評價結果的能力。ABAQUS/Viewer是ABAQUS/CAE的子集，它具有后處理功能，這些模塊之間的關系見圖2- 1 圖2-1 ABAQUS/Standard ABAQUS/Standard是一個通用分析模塊，在數(shù)值方法上采用有限元方法常用的 隱式積分。它能夠求解廣泛的線性和非線性問題，包括結構的靜態(tài)、動態(tài)問題、熱 力學場和電磁場問題等。對于通常同時發(fā)生作用的幾何、材料和接觸非線性可以采 用自動控制技術處理，也可以由用戶自己控制。 ABAQUS/Explicit ABAQUS/Explicit是一個在數(shù)值方法上采用有限元顯式積分的特殊模塊，它利用對時間的顯式積分求解動態(tài)有限元方程。它適合于分析諸如沖擊和爆炸這樣短暫、瞬時的動態(tài)問題，同時對高度非線性問題如模擬加工成型過程中接觸條件的改變等也非常有效。 ABAQUS/CAE ABAQUS/CAE是ABAQUS進行有限元分析的前后處理模塊，也是建模、分析和后處理的人機交互平臺。該模塊根據(jù)結構的幾何圖形生成網(wǎng)格，將材料和截面的特性分配到網(wǎng)格上，并施加載荷和邊界條件。該模塊可以進一步將生成的模型投入到分析模塊中進行高效率的后臺運行，并對運行情況進行監(jiān)測，對計算結果進行后處理。ABAQUS/CAE的后處理支持ABAQUS分析模塊的所有功能，并且對計算結果的描述和解釋提供了范圍很廣的選擇，除了通常的云圖，等值線和動畫顯示之外，還可以用列表，曲線(包括部分常用運算)等其他常用工具來完成對結果數(shù)據(jù)的處理。該模塊的許多獨特功能與特點，例如CAD特征化建模、參數(shù)化建模、適應設計者要求的數(shù)據(jù)管理系統(tǒng)等極大的方便了ABAQUS的使用者。 ABAQUS/Aqua ABAQUS/Aqua的一系列功能可以附加在ABAQUS/Standard中應用。它偏向于模擬海上結構，如海洋石油平臺。它的功能包括模擬波浪，風載荷及浮力的 影響。在本指南中不討論ABAQUS/Aqua。 ABAQUS/ADAMS ABAQUS/ADAMS允許ABAQUS有限元模型作為柔性部件進入到MDIADAMS產(chǎn)品族中去進行分析。 ABAQUS/C-MOLD ABAQUS/C-MOLD把注模分析軟件C-MOLD中有限元網(wǎng)格、材料性質(zhì)和初始應力數(shù)據(jù)轉(zhuǎn)換成為ABAQUS 輸入文件。 ABAQUS/Design ABAQUS/Design 的一系列功能可附加在ABAQUS/Standard 中進行設計敏度計算。 ABAQUS/MOLDFLOW ABAQUS/MOLDFLOW 模塊把MOLDFLOW 分析軟件中的有限元模型信息 轉(zhuǎn)換成ABAQUVS 輸入文件的一部分。 2.3. ABAQUS的二次開發(fā)平臺 ABAQUS的腳本語言接口非常友好，其自嵌的腳本語言是Python[4]，系國際上廣泛使用、功能強大、具有良好開放性的一種面向?qū)ο蟪绦蛟O計語言。所以，應用Python在ABAQUS中進行二次開發(fā)也比較方便，且可移植性強。ABAQUS以基于Python的語法規(guī)則向二次開發(fā)者提供了許多庫函數(shù)，這些庫函數(shù)主要是用來增強ABAQUS的交互式(GUI)操作功能。用戶可以通過ABAQUS的交互式(GUI)界面實現(xiàn)分析對象的特征造型、指定材料屬性、完成網(wǎng)格剖分和控制、提交并監(jiān)控分析作業(yè)，也可以使用ABAQUS腳本語言越過ABAQUS的交互式(GUI)界面直接高效地向ABAQUS內(nèi)核提交任務。使用Python可以進行參數(shù)化建模，修改交互式建立的模型，還可以一次提交多個作業(yè)。 出了腳本語言接口，ABAQUS還為用戶提供了功能強大的用戶子程序接口（Abaqus User Subroutines ），以幫助用戶開發(fā)基于ABAQUS內(nèi)核的程序，常用的用戶子程序包括UEL(User subroutine to define an element ，用戶單元子程序)，UMAT(User subroutine to define a materials mechanical behavior，用戶材料子程序 )[5]，其中UMAT的使用最為廣泛，它主要用于用戶開發(fā)自己的材料模型，以彌補ABAQUS自帶材料模型的不足，幫助用戶完成各種材料分析，功能極為強大。 在國外，眾多的有限元分析和研究者熱衷于使用ABAQUS，一個很重要的原因就在于ABAQUS給用戶提供了功能強大，使用方便的二次開發(fā)工具和接口，使得用戶可以方便的進行富含個性化的有限元建模、分析和后處理，滿足特定工程問題的需要。通過用戶材料子程序接口，用戶可定義任何補充的材料模型，不但任意數(shù)量的材料常數(shù)都可以作為資料被讀取，而且ABAQUS對于任何數(shù)量的與解相關的狀態(tài)變量在每一材料計數(shù)點都提供了存儲功能，以便在這些子程序中應用。 2.4. ABAQUS的二次開發(fā)語言 ABAQUS的二次開發(fā)語言主要有3種：Python，F(xiàn)ORTRAN，C++ Python語言主要用于GUI開發(fā)，F(xiàn)ORTRAN語言主要用于用戶子程序開發(fā)，而c++語言主要專注于其他高級開發(fā)部分。 本文主要是針對用戶子程序的開發(fā)，所以采用FORTRAN語言，下面簡要介紹一下該語言極其特點： FORTRAN語言是世界上第一個被正式推廣使用的高級語言[6]。它是1954年被提出來的，1956年開始正式使用，至今已有三十多年的歷史，但仍歷久不衰，它始終是數(shù)值計算領域所使用的主要語言。 FORTRAN語言是Formula Translation的縮寫，意為“公式翻譯”。它是為科學、工程問題或企事業(yè)管理中的那些能夠用數(shù)學公式表達的問題而設計的，其數(shù)值計算的功能較強。 FORTRAN語言問世以來，根據(jù)需要幾經(jīng)發(fā)展，先后推出了不同的版本，主要版本有FORTRAN 77，F(xiàn)ORTRAN 90，F(xiàn)ORTRAN 95，ABAQUS采用FORTRAN 77，通常用固定格式編寫代碼。 FORTRAN77語言同C語言一樣，是一種結構化編程語言 結構化程序設計方法規(guī)定，在結構化的程序中，只能有三種基本結構： (1)順序結構 這是一種最簡單的基本結構形式，它的特點是，在這個結構內(nèi)的各個功能模塊或語句序列，是按其出現(xiàn)的先后順序執(zhí)行的，如賦值語句、輸入/輸出語句等。它有一個入口和一個出口，并在入口和出口之間包含著若干個功能塊，其中每一個功能塊可以是一個非轉(zhuǎn)移語句。因此，順序基本結構塊是由一系列的順序執(zhí)行語句組成的。 (2)分支選擇結構 在給定的條件下，分支選擇結構判斷選擇哪一條路徑執(zhí)行，不同路徑完成的功能是不同的。實現(xiàn)分支選擇結構主要由塊IF語句、ELSE語句、END IF語句以及ELSE IF語句組成的IF-THEN-ELSE結構。 (3)循環(huán)結構 循環(huán)結構也稱重復處理結構，即重復執(zhí)行某一功能塊，直到滿足（或不滿足）某一條件為止。實現(xiàn)循環(huán)結構的FORTRAN90語句主要是DO語句、塊IF語句和邏輯IF語句的結合。 以上三種基本結構，是組成結構化程序的基本結構形式。這里有兩層意思：一是結構化的程序中，各個模塊均由這三種基本結構組成；二是結構化程序本身，從宏觀上也是這三種基本結構形式之一。 3. 用戶材料子程序UMAT 3.1. UMAT開發(fā)環(huán)境設置 由于UMAT是采用FORTRAN語言編寫，那么要運行UMAT就需要安裝FORTRAN的開發(fā)環(huán)境， 同時還需要ABAQUS的支持，本文采用的ABAQUS版本為6.81，支持INTEL Fortran9.1-10.1，Intel Fortran安裝時又需要安裝Microsoft Visual Studio的相應版本，經(jīng)過比較，本文選用ABAQUS6.81+Intel Fortran10.1+Microsoft VisualC++ 2005，相對于ABAQUS來說，UMAT開發(fā)環(huán)境的設置較為繁瑣，這給子程序的使用帶來諸多不便，為了解決這一問題，我用C#語言編制了ABAQUS子程序編譯環(huán)境設置工具，只需要將安裝文件解壓到ABAQUS的安裝目錄，運行安裝程序就可以了，整個過程不需要人工干預，也不需要安裝龐大的VisualC++ 2005，如圖3-1所示 圖3-1 3.2. UMAT注意事項 ABAQUS的用戶子程序是根據(jù)ABAQUS提供的相應接口，按照Fortran語法， 用戶自己編寫的代碼。它是一個獨立的程序單元，可以獨立的被存儲和編譯，也能 被其它程序單元引用，因此，利用它可帶回大量數(shù)據(jù)供引用程序使用，也可以用它 來完成各種特殊的功能。它的一般結構形式是： SUBROUTINE S(x1,x2,……,xn) INCLUDE‘ABA_PARAM.INC’（用于ABAQUS/Standard用戶子程序中） OR INCLUDE‘VABA_PARAM.INC’)（用于ABAQUS/Explicit用戶子程序中） …… RETURN END x1，x2，……，xn是ABAQUS提供的用戶子程序的接口參數(shù)，有些參數(shù)是ABAQUS傳到用戶子程序中的，例如SUBROUTINE DLOAD中的KSTEP、KINC、COORDS，有些是需要用戶自己定義的，例如F，文件aba_param.inc和vaba_param.inc隨著ABAQUS軟件的安裝而包含在操作系統(tǒng)中，它們含有重要的參數(shù)，幫助ABAQUS主求解程序?qū)τ脩糇映绦蜻M行編譯和鏈接。當控制遇到RETURN語句時便返回到引用程序單元中去，END語句是用戶子程序結束的標志。 在一個算例中，用戶可以用到多個用戶子程序，但必須把它們放在一個以.for為擴展名的文件中。運行帶有用戶子程序的算例同時有兩種方法：一是在CAE中運行，在EDIT JOB菜單中的GENERAL子菜單的USER SUBROUTINE FILE對話框中選擇用戶子程序所在的文件即可；另外是在ABAQUS.COMMAND中運行，語法如下： abaqus job=job-name user={source-file|object-file} 編制用戶子程序時應注意 (1)用戶子程序相互之間不能調(diào)用，但可以調(diào)用用戶自己編寫的Fortran子程序和ABAQUS應用程序。ABAQUS應用程序必須由用戶子程序調(diào)用。當用戶編寫Fortran子程序時，建議子程序名以K開頭，以免和ABAQUS內(nèi)部程序沖突。 (2)當用戶在用戶子程序中利用OPEN打開外部文件時，要注意以下兩點：一是設備號的選擇是有限制的，只能取15～18和大于100的設備號，其余的都已被ABAQUS占用；二是用戶需提供外部文件的絕對路徑而不是相對路徑。 (3)對于不同的用戶子程序ABAQUS調(diào)用的時間是不同的，有的是在每個STEP 的開始，有的是STEP的結尾，有的是在每個INCREMENT的開始等等。當ABAQUS 調(diào)用用戶子程序時，都會把當前的STEP和INCREMENT利用用戶子程序的兩個實 參KSTEP和KINC傳給用戶子程序，用戶可把它們輸出到外部文件中，這樣就可清 楚的知道ABAQUS何時調(diào)用該用戶子程序。 為保證用戶子程序的正確執(zhí)行，子程序的書寫必須遵循ABAQUS的相關規(guī)定， 下面以用戶材料子程序為例詳細說明。 3.3. UMAT接口的原理 用戶材料子程序（User-defined Material Mechanical Behavior，簡稱UMAT）是 ABAQUS提供給用戶定義自己的材料屬性的Fortran程序接口[7][8]，使用戶能使用 ABAQUS材料庫中沒有定義的材料模型。用戶材料子程序UMAT通過與ABAQUS主求解程序的接口實現(xiàn)與ABAQUS的資料交流。在輸入文件中，使用關鍵詞“*USER MATERIAL”表示定義用戶材料屬性。 UMAT子程序具有強大的功能，使用UMAT子程序： （1）可以定義材料的本構關系，使用ABAQUS材料庫中沒有包含的材料進行 計算，擴充程序功能。 （2）幾乎可以用于力學行為分析的任何分析過程，幾乎可以把用戶材料屬性賦 予ABAQUS中的任何單元。 （3）必須在UMAT中提供材料本構的雅可比（Jacobian）矩陣，即應力增量對 應變增量的變化率。 由于主程序與UMAT之間存在數(shù)據(jù)傳遞，甚至共享一些變量，因此必須遵守有 關UMAT的書寫格式，UMAT中常用的變量在文件開頭予以定義，通常格式 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‘ABA_PARAM.INC’ CHARACTER*80 CMNAME DIMENSION 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,SCD and,if necessary,RPL,DDSDDT,DRPLDE,DRPLDT,PNEWDT RETURN END UMAT中的應力矩陣、應變矩陣以及矩陣DDSDDE、DDSDDT、DRPLDE等，都是直接分量存儲在前，剪切分量存儲在后。直接分量有NDI個，剪切分量有NSHR個。各分量之間的順序根據(jù)單元自由度的不同有一些差異，所以編寫UMAT時要考慮到所使用單元的類別。下面對UMAT中用到的一些變量進行說明： DDSDDE(NTENS NTENS)：一個NTENSNTENS的矩陣，稱作Jacobian矩陣 是應力的增量，是應變的增量，DDSDDE(i，j)表示增量步結束時第j個應變分量的改變引起的第i個應力分量的變化。通常Jacobian矩陣是一個對稱矩陣，除非在“*USER MATERIAL”語句中加入了“UNSYMM”參數(shù)。 STRESS(NTENS)：應力張量數(shù)組，對應NDI個直接分量和NSHR個剪切分量。在增量步的開始，應力張量矩陣中的數(shù)值通過UMAT和主程序之間的接口傳遞到UMAT中，在增量步的結束UMAT將對應力張量矩陣更新。對于包含剛體轉(zhuǎn)動的有限應變問題，一個增量步調(diào)用UMAT之前就已經(jīng)對應力張量進行了剛體轉(zhuǎn)動，因此UMAT中只需處理應力張量的共旋部分。UMAT中應力張量的度量為柯西應力。 STATEV(NSTATEV)：用于存儲與解有關的狀態(tài)變量的數(shù)組，在增量步開始時將數(shù)值傳遞到UMAT中，也可在子程序USDFLD或UEXPAN中先更新數(shù)據(jù)，然后增量步開始時將更新后的資料傳遞到UMAT中。在增量步的結束必須更新狀態(tài)變量矩陣中的數(shù)據(jù)。和應力張量矩陣不同的是：對于有限應變問題，除了材料本構行為引起的資料更新以外，與解有關的狀態(tài)變量矩陣中的任何向量或者張量都必須通過旋轉(zhuǎn)來考慮材料的剛體運動。狀態(tài)變量矩陣的維數(shù)通過ABAQUS輸入文件中的關鍵詞“*DEPVAR”定義，關鍵詞下面數(shù)據(jù)行的數(shù)值即為狀態(tài)變量矩陣的維數(shù)。 PROPS(NPROPS)：材料常數(shù)數(shù)組。材料常數(shù)的個數(shù)，等于關鍵詞“*USER MATERIAL”中“CONSTANTS”常數(shù)設定的值。矩陣中元素的數(shù)值對應于關鍵詞“USER MATERIAL”下面的數(shù)據(jù)行。 SSE，SPD，SCD：分別定義每一增量步的彈性應變能，塑性耗散和蠕變耗散。 它們對計算結果沒有影響，僅僅作為能量輸出。 STRAN(NTENS)：應變數(shù)組。 DSTRAN(NTENS)：應變增量數(shù)組。 DTIME：增量步的時間增量。 NDI：直接應力分量的個數(shù)。 NSHR：剪切應力分量的個數(shù)。 NTENS：總應力分量的個數(shù)，NTENS＝NDI＋NSHR。 由于UMAT子程序在單元的積分點上調(diào)用，增量步開始時，主程序路徑將通過 UMAT的接口進入UMAT，單元當前積分點必要變量的初始值將隨之傳遞給UMAT 的相應變量。在UMAT結束時，變量的更新值將通過接口返回主程序。 3.4. UMAT的使用方法 我們知道，有限元計算（增量方法）的基本問題[7]是：已知第n 步的結果（應力，應變等），，然后給出一個應變增量，計算，UMAT要完成這一計算，并要計算DDSDDE(I,J)=。是應力增量矩陣，是應變增量矩陣，DDSDDE(I,J) 定義了第J 個應變分量的微小變化對第 I 個應力分量帶來的變化。 該矩陣只影響收斂速度，不影響計算結果的準確性 （當然，不收斂自然得不到結果）。 有限元計算的中心問題就是求得節(jié)點的位移 （進而應變、應力），以使內(nèi)力和 外力達到平衡： （3-1） d 是節(jié)點位移矩陣，黑體字表示矩陣或矢量。除了小變形、線彈性問題，方程2-1是非性的，要用迭代的方法解出： （3-2） (3-3) i 表示一個增量步內(nèi)的第i 次迭代，n表示第n個增量步。KT是切線剛度，由材料的Jacobian 矩陣結合單元計算組裝而得。剛度矩陣其實就是力對位移的梯度。要想快速收斂，位移增量應沿該梯度方向變化，也就是說，如果Jacobian 矩陣不是那么準確，自然KT 也不怎么準確，那么滿足3-1式的位移被找到的速度也就變慢，甚至發(fā)散，根本找不到。但收斂速度無論慢快，3-1式才是判斷結果準確與否的唯一標準。所以Jacobian 矩陣不影響結果的準確性，只影響收斂速度的快慢。 以大變形、非性材料為例，整個計算步驟是這樣的： 整個外力不是一次加上，而是一點點加上的，不然會發(fā)散得不到結果的。所以，每一個增量步開始時就是在原來的外力上加上一點點，得到。根據(jù)3-2得到位移增量，此時要知道力對位移的梯度KT，以盡快找到滿足平衡條件的位移，由材料的Jacobian 矩陣和單元結合起來組裝得到(此處使用UMAT 提供的Jacobian 矩陣)。然后可計算應變增量，調(diào)用UMAT，得到新的應力，進而得到新的內(nèi)力，所以，程序不在乎新的應力是由增量方法得到，還是全量方法得到，而只在乎新應力是否準確。然后回到3-2，如此循環(huán)，直至3-2右端為0，也即滿足3-1。這樣第n+1 步就完成了，然后開始第n+2 步，即 外力加上一點點，按同樣的方法求解新的位移。直至整個外力全部施加并得到滿足3-1的位移。 4. 材料非線性問題 彈性力學作為精確理論，從本質(zhì)上都是非線性的，早期Cauchy, Green,Kirchhoff和Kelvin在這些方面都作出了重要貢獻。后來又提出超彈性(即具有彈性勢的)有限變形理論，由于理論方程的冗長而復雜，且工程應用也沒有提出這方面要求而被擱置。20世紀40年代以后，由于橡膠材料、高分子合成材料的迅速發(fā)展和工業(yè)領域的大量應用，非線性彈性與超彈性的研究再次引起科學和工程界的重視，除了一般理論研究有了新的發(fā)展以外，工程應用計算方法也得到長足的發(fā)展。 4.1. 材料的彈塑性本構關系 彈塑性材料進入塑性的特征是當荷載卸去后存在不可恢復的永久變形。所以，在卸載情況下，應力應變之間不再是唯一的對應關系。這是區(qū)別于非線性彈性材料的基本屬性。只以加載時應力應變關系成非線性，還不足以判斷材料是非線性彈性還是彈塑性。但是一經(jīng)卸載就可以看出兩者的區(qū)別。非線性彈性材料沿原路徑返回，而彈塑性材料將依據(jù)不同的加載歷史卸載后產(chǎn)生不同的永久變形。 對大多數(shù)材料來說，在單調(diào)加載的情況下，存在一個明顯的極限應力，當應力低于時，材料保持線彈性。而當應力達到以后，則材料開始進入彈塑性狀態(tài)。如繼續(xù)加載，然后在卸載，材料始終保持永久的塑性變形。如果應力達到后，應力不再增加，而材料變形可以繼續(xù)增加，及變形處于不定的流動狀態(tài)，則稱材料為理想彈塑性的。反之如果應力達到后，再增加變形，應力也必須增加，則材料是應變硬化的，這時應力是塑性應變的函數(shù)，可解析為： (4-1) 本構關系反應著應力應變之間的關系。對于彈性材料變形是可以恢復的；而塑性材料變形是不可以恢復的。典型的彈塑性應變在卸載后要保持一個永久的變形。如圖3-2 圖4-1 塑性應變有下列特性： （1）總應變分為彈性和塑性兩部分，即 (4-2) 或者： (4-3) （2）塑性變形取決于加載路徑，而應力應變之間沒有一一對應的關系。所以必須確定二則之間的本構關系，這種本構關系可以用偏微分方程或者增量形式來描述。 總之，彈塑性理論主要包括以下幾個方面： （1）應變張量的分解； （2）應力空間的屈服條件； （3）流動法則； （4）強化法則； （5）協(xié)調(diào)性條件。 1：本構模型 塑性力學的應力-應變曲線通常有5種簡化模型[8]： （1）理想彈塑性模型，用于低碳鋼或強化性質(zhì)不明顯的材料。 （2）線性強化彈塑性模型，用于有顯著強化性質(zhì)的材料。 （3）理想剛塑性模型，用于彈性應變比塑性應變小得多且強化性質(zhì)不明顯的材料。 （4）線性強化剛塑性模型，用于彈性應變比塑性應變小得多且強化性質(zhì)明顯的材料。 （5）冪強化模型，為簡化計算中的解析式，可將應力-應變關系的解析式寫為 σ＝σy（ε／εy）n，式中σy為屈服應力，εy為與σy相對應的應變，n為材料常數(shù)。 圖4-2 2：屈服條件 在復雜應力狀態(tài)下，判斷物體屈服狀態(tài)的準則稱為屈服條件[9]。屈服條件是各應力分量組合應滿足的條件。對于金屬材料，最常用的屈服條件為最大剪應力屈服條件（又稱Tresca屈服條件）和彈性形變比能屈服條件（又稱Von Mises條件）。對于巖土材料則常用Tresca屈服條件、Drucker-Prager屈服條件和Mohr-Coulomb屈服條件。對于強化或軟化材料，屈服條件將隨塑性變形的增長而變化，改變后的屈服條件稱為后繼屈服條件。當已知主應力的大小次序時，使用Tresca屈服條件較為方便；若不知道主應力的大小次序，則使用Von Mises屈服條件較為方便。對于韌性較好的材料，Von Mises屈服條件與試驗數(shù)據(jù)符合較好。 Von Mises屈服準則具體形式是，對于各項同性材料，應力偏量第二不變量等于某一定值時，材料開始進入了塑性狀態(tài)。 (4-4) 3：強化法則 對理想的彈塑性材料而言，因無強化作用，所以，整個塑性變形過程中，屈服函數(shù)值保持一個常量，強化定義了屈服面在應力空間的演化準則。 (4-5) 其中，是強化參數(shù)。 通常采用的強化法則有以下幾種： (1) 各向同性強化 此法則規(guī)定材料進入塑性變形以后，加載曲面在各方向均勻的向外擴張，沒有畸變。而其形狀、中心及其在應力空間的方位均保持不變[10]。需要指出的是：各向同性強化法則主要適用于單調(diào)加載情況。如果用于卸載情況，它只適合反向屈服應力等于應力反轉(zhuǎn)點的材料，而通常材料不具備這種性質(zhì)，因此在塑性力學中還發(fā)展了其它強化準則。 (2) 隨動強化 此法則規(guī)定材料進入塑性狀態(tài)以后，加載曲面在應力空間作剛體移動而沒有轉(zhuǎn)動，因此初始屈服面的形狀、大小和方向仍然保持不變。 (3) 混合強化 把各向同性強化模型和隨動強化模型加以組合，得到混合強化模型。它假定在塑性變形過程中，加載曲面不但作剛性平移，還同時在各個方向作均勻擴大。 在以上幾種強化模型中，各向同性強化模型應用最為廣泛。本文也是采用該硬化法則，這一方面是由于它便于進行數(shù)學處理；另一方面，如果在加載過程中應力方向（或各個應力分量的比值）變化不大，采用各向同性強化模型的計算結果與實際情況也比要符合。隨動強化模型可以考慮材料的包興格（Bauschinger）效應，在循環(huán)加載或可能出現(xiàn)反向屈服的問題中，需要采用這種模型。 由于塑性變形與變形歷史有關， 因此反映塑性應力-應變關系的本構關系用應變增量形式給出比較方便。用應變增量形式表示塑性本構關系的理論稱為塑性增量理論。增量理論的本構關系在理論上是合理的，但應用比較麻煩，因為要積分整個變形路徑才能得到最后結果。因此，又發(fā)展出塑性全量理論，即采用全量應力和全量應變表示塑性本構關系的理論。在比例變形的條件下，可通過積分增量理論的本構關系獲得全量理論的本構關系。當偏離比例變形條件不多時，全量理論的計算結果和實險結果比較接近。本文的程序都是基于增量理論。 4.2. 非線性有限元算法理論 對于非線性問題，在有限元求解該問題時，對一個自由度總可以表達成，式中，為基本未矢量。如果是線性問題，與無關，而是一次項，顯然這是一個線性方程。如果與相關，則方程的出現(xiàn)非一次項，變成非線性問題，在實際工程中，特別是塑性成型問題，材料的幾何方程，本構方程以及邊界條件往往是非線性的也體現(xiàn)在中出現(xiàn)了，所以變?yōu)榱朔蔷€性問題，要得到最基本的未知量，就必須求解非線性方程組 1：直接迭代法 又稱常剛度法[11]，這是種最簡.單的求解方法，在每次求解前，利用上次的解來求出這一次的值，然后利用和的倒數(shù)的乘積求出的當前值 (4-6) 表達為迭代形式 (4-7) 上式可以看出，這種方法首先需要有一個初始的值，以便開始迭代。另外，每一次求解都需要對求倒數(shù)，如果求解方程組，就是對剛度矩陣求逆，這種方法在求解中控制兩次求解之差，當其值很小時，就認為接近真實值了，迭代結束 圖4-3 2：Newton-Raphson方法 Newton-Raphson方法的算法與常剛度法不同[12]，如果得近似表達式是不成立的，存在著殘余值，即，此式也可以作為近似值與真實值的差值量度，實際上在具體計算時，也可以控制其值，當極小時，就認為接近真實值了，當?shù)诖蔚闹? 是真實解，則可以按照Taylor級數(shù)展開得到 圖4-4 3：切線剛度法 在復雜非線性問題求解中，剛度與的大小是有一定關系的，在用增量法來求解這種問題時，就等于結構任一點處力與位移的曲線的局部梯度，稱為切線剛度[13]，剛度矩陣的倒數(shù)很難用自變量顯示表達，通過增量方式求解，在每一步荷載增量范圍內(nèi)把問題線性化，求解方法與Newton-Raphson方法相同 總結以上可以得到： 以上幾種算法中，通過比較，不難發(fā)現(xiàn)，直接迭代法采用了固定的剛度，適合解決非線性程度不高的本構關系，而切線剛度法采用了變化的剛度，在每一步上都做了實時的修正，對非線性程度較高本構關系任然有效，在效率和迭代精度方面，切線剛度法采用的修正更符合非線性材料的應力應變關系，具有較大的優(yōu)勢，這也是本文采用切線剛度法計算的原因，當然，非線性有限元算法還有很多，切線剛度法也不見得就是最好的能解決所有問題的算法，但是它是在程序開發(fā)難度不高和精度方面較高的條件下相對來說最好的 本文采用的本構關系是同性硬化彈塑性模型[14]，采用Mises屈服準則，下面將根據(jù)J2理論[15]，分別推導常剛度法和切線剛度法計算該問題的的算法公式 4.3. 增量理論常剛度法公式推導 由應力應變關系得： (4-8) (4-9) 其中是彈性矩陣，它的表達式為 是塑性應變增量，它的表達式為 其中為等效塑性應變增量，它的表達式為 為切線模量 對于3維空間問題，流動方向： 等效應力 應力偏量 4.4. 增量理論切線剛度法公式推導 本文采用的是一種切線剛度法，其應力應變關系為 (4-10) (4-11) 3-14兩端左右同乘 ，得到： 對于強化材料，Mises準則 代入可得 其中，為等效塑性應變增量，它的表達式為 由于： 由流動法則可知 應用上式得 所以得到 (4-12) 這就是切線剛度法的矩陣表達式 其中為彈塑性矩陣，在ABAQUS里面稱為雅可比矩陣，它的表達式為 其中是彈性矩陣，它的表達式為 為塑性矩陣，它的表達式為 對于3維空間問題 流動方向 所以 最后推導可得，彈塑性矩陣的表達式 其中 為切線模量，對本構關系求導得到 總結推導過程：上面的推導過程看似復雜，其實核心的問題只有一個，即兩者對于塑性階段應力更新的算法不同。 常剛度法采用的是： 而切線剛度法采用的是： 把握好了兩者的本質(zhì)上的區(qū)別，對于兩者的算法設計和程序開發(fā)問題便迎刃而解 5. UMAT程序設計和編碼 本章將嚴格按照前一章推導的公式展開程序設計和編碼，為了便于編程，本文將本構關系做了抽象化處理，即將其描述成一個含參數(shù)的表達式，改變參數(shù)即可應用于不同的模型，這樣做的好處是能保證程序的復用性，這也是本文反復強調(diào)的使用UMAT的原則。 5.1. 本構關系描述 本文采用各向同性硬化彈塑性材料，材料參數(shù)如下： (5-1) 彈性部分：，彈性模量E=200000Mpa，泊松比Mu=0.3 圖5-1 彈性部分本構關系 塑性部分：，為了研究方便，取A=700，B=0.5，C=400 圖5-2 塑性部分本構關系 將2個曲線統(tǒng)一到同一個坐標系（為方便顯示，x軸標注時擴大了1000倍） 圖5-3 本構關系 由此可以求出兩條線的交點即初始屈服點的應力Yield0=400Mpa（注意兩條曲線相差了一個屈服應變，因為兩者其實不是一個坐標系） 綜上定義的材料常數(shù)見表5-1： 表5-1 彈性模量E 200000 泊松比Mu 0.3 屈服應力Yield0 400 A 700 B 0.5 C 400 注意：上面A,B,C的取值只是為了便于理解和分析本文的材料模型，為了保證程序的 通用性，本文的參數(shù)在程序中的A,B,C一律用變量表示。 5.2. 常剛度法程序設計 算法設計 1：定義程序需要用到的常數(shù)和變量 2：讀取ABAQUS定義的材料常數(shù)和狀態(tài)變量(這里只定義了一個狀態(tài)變量)，材料常數(shù)為，彈性模量E，泊松比Mu，屈服應力Yield0，參數(shù)A,B,C,并且計算出剪切模量G,狀態(tài)變量為等效塑性應變EQPLAS 3：讀取應力分量，計算平均應力，應力偏量以及Mises等效應力 平均應力： 應力偏量： Mises等效應力： 4：根據(jù)3計算的Mises等效應力和2讀取的屈服應力Yield0比較，如果Mises等效應力小于屈服應力，表明此時材料未屈服，那么轉(zhuǎn)到5，否則轉(zhuǎn)到6 5：雅可比矩陣,初始化為0，計算彈性矩陣，按照彈性理論更新應力 6：雅可比矩陣,初始化為0 1） ：計算切線模量H ,注意到當?shù)刃苄詰獣r對應于本構關系的屈服點，此時的H不能通過上式計算，可以取此時的H為彈性模量 2） ：計算等效塑性應變增量并更新 EQPLAS=DEQPLAS+EQPLAS更新狀態(tài)變量 3）計算流動方向 4） 計算塑性應變增量 5）更新應力 算法流程圖 圖5-4 常剛度法算法流程圖 5.3. 常剛度法程序編碼 根據(jù)算法流程，用FORTRAN77 固定格式編制了常剛度法的計算程序如下： SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,DDSDDT, 1 DRPLDE,DRPLDT,STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED, 2 CMNAME,NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT, 3 PNEWDT,CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC) INCLUDE ABA_PARAM.INC CHARACTER*8 CMNAME DIMENSI