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

《ABAQUS-Fortran二次開(kāi)發(fā).doc》由會(huì)員分享，可在線閱讀，更多相關(guān)《ABAQUS-Fortran二次開(kāi)發(fā).doc（73頁(yè)珍藏版）》請(qǐng)?jiān)谘b配圖網(wǎng)上搜索。

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