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

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