ABAQUS非線性有限元分析與實例umat

1、第19章ABAQUS用戶材料子程序（UMAT）雖然ABAQUS為用戶提供了大量的單元庫和求解模型，使用戶能夠利用這些處理絕大多數(shù)的問題；但是現(xiàn)實世界畢竟十分復雜*ABAQUS不可能把所有可能的問題都包含進去。所以ABAQUS提供了大量的用戶子程序（user覘子程序允許用戶在找不到合適模型的情況下自行定義符合自己問題的模型*這些用程序涵蓋了建模從載荷到單元的幾乎各個部分。ABAQUS為用戶提供的這個接口，允許用戶通過自定義的子程序定制ABAQ以實現(xiàn)特定的功能。用戶子程序具有以下的功能和特點；1）如果ABAQUS的一些固有選項模型功能有限，用戶子程序可以提高ABA，中這些選項的功能。2）通常用戶子

2、程序是用Fortran語言的代碼寫成。3）它可以以幾種不同的方式包含在模型中。4）由子它們沒有存儲在restart文件中，如果需要的話，可以在重新開始運行I改它。5）在某些情況下它可以利用ABAQUS允許的已有程序要在模型中包含用戶子程序，可以利用ABAQUS執(zhí)行程序，在ABAQUS執(zhí);序中應用user選項指明包含這些子程序的Fortran源程序或者目標程序的名字提示：ABAQUS的輸入文件除了可以通過ABAQUS/CAE的作業(yè)摸塊提交運古還可以在ABAQUSCommand口中輸入ABAQUS執(zhí)疔程序直接運疔：ABAQUSjob二輸入文件名user-用戶予程序的Fortran文件名ABAQUS

3、/Standard和ABAQUS/Explicit都支持用戶子程序功能，但是它們所的用戶子程序種類不盡相同，讀者在使用時請注意查詢手冊。在接下來的兩章中，我們將討論兩種常用的用戶子程序用戶材料子程序和單元子程序。本章將通過在ABAQUS/Standard中創(chuàng)建Johnson-Cook的材料模型，對編寫Sdard的用戶材料子程序UMAT進行簡單介紹。ABAQUS/Explic讓中的用戶材料子1VUMAT的思想與之相似，但是由于隱式和顯式兩種方法本身的差異”它們之間一些不同，請讀者在自己具體使用前首先仔細查閱ABAQUS手冊中的相關內(nèi)容。19-1引言用戶材料子程序是ABAOUS握供給陽戶白宗義材料

4、屬件的Fnrtrn稈序榕口ao=A+B(19ao=A+B(19Eee=E(Eep)(過ABAQUS/Sundard的用戶材料子程序UMAT編程實現(xiàn)。在UMAT編程中使用相關塑性理論以及完全隱式的應力更新算法。19.2模型的數(shù)學描述19*2.1Johnson-Cook強化模型簡介Johnson-Cook模型用來模擬在沖擊載荷作用下的變形。Johnson-Cook強化模型(JC)表示為三項的乘積，分別反映了應變硬化、應變化和溫度軟化。這里使用JC模型的修正形式：(1a=(A+Ben)1+Cln1+(1一模型中包含A.n.C、椚五個參數(shù)，需要通過實驗來確定。使參考應變率&這樣公式中的A即為材料的靜態(tài)

5、屈服應力。公式中的:T為無量綱化的溫度：T.r-rrTm-TT其中丁為室溫,Tm為材料的熔點。Johnson-Cook模型在溫度從室溫到材料熔點的范圍內(nèi)都是有效的。高應變率的變形經(jīng)常伴有溫升現(xiàn)象，這是因為材料變形過程中塑性功轉(zhuǎn)化為?對于大多數(shù)金屬，90%-100%的塑性變形將耗散為熱量。所以JC模型中溫度的喪以用如下的公式計算：丁=2血任(1其中，丁為溫度的增量；盤為塑性耗散比，表示塑性功轉(zhuǎn)化為熱量的比例；為的比熱，p為材料的密度口公式(19-2)考慮的是一個絕熱過程，即認為溫度的？|全起因于塑性耗散。19.2.2率相關塑性的基本公式JohnsomCook本構(gòu)模型考慮率相關塑性，塑性變形是關聯(lián)

6、的，即塑性流動沿看面的法線方向，并采用Mises屈服面。將應變的增量分解為彈性部分和塑性部分：de=dee+dep(1將上式兩端同時對時間的增量df微分得到率形式：右二于+尹(1在率相關塑性中，材料的塑性反應取決于加載率，以率的形式給岀材料的弾應；上式即為流動法則，其中誓為塑性流動勢能，入為塑性率參數(shù)*“為運動硬化時的力。對于各項同性硬化，不存在背應力，因此有0=0,此時有(r誓hI7I二crsign(d-)(I1與率無關塑性不同的是,而是直接通過經(jīng)驗定律給出，sign(7)率相關塑性中等效塑性應變率不能通過一致性條件老成為過應力模型：A(F式中0是過應力，可為粘性。在過應力模型中，等效塑性應

7、變率取決于超過了多少應力。上而一維的率相關塑性公式可以很方便地推廣到三維情況。對于小應變的情形力度量之間無須區(qū)分，這里采用Cauchy應力6塑性率參數(shù)由應力和內(nèi)變量的經(jīng):數(shù)給出。對照一維情況，三維情況下分解應變率為彈性和塑，性部分：(7=了).WF)(19s=ee+(19應力率和彈性應變率之間的關系為F=0,則=/;如果/0,則(f)二0擊為Mises等:力,云為等效應變，神為率敏感系數(shù)。對于Johnson-Cook模型，可以得到等效塑性應變率的表達式為(19ao=A+B(19ao=A+B(19Eee=E(Eep)(ao=A+B(19ao=A+B(19Eee=E(Eep)(其中為靜態(tài)屈眼應力:

8、Eee=E(Eep)(TOC o 1-5 h zsn+l二兔+壯(1屠+1=叮十為+1G+1(1盛+Q=騎+九+1.九十1(1%+1二C：(畸+1礁十1)(1fni-/(rn+in+i)=0(1在時刻肚給出一組（馮，?,qn）和應變增量爲公式（19-17）是一組關（為+“瑞務I）的非線性代數(shù)方程。將公式（19-17b）代入（19-17d）得到為十1=C：（程n+lr畔+1）=C：務+醇一碼j二C：（en一咄）+C：Ac口公十1=(r+C：Ac)-C：典+itrailj-6+1一二典+t=AA(0)=O,(f(o)=Cj略+護)應力在第k次迭代時為ff(A=rt-co)一AA(i)C護(勵且在整

9、個算法的塑性修正狀態(tài)過程中始終保持不變，因此塑性應變的更新是久的參函數(shù)。在k次迭代時將檢查屈服條件：(19-：(19；嚴）二存的一aY（&）二（笄）-3#AA）卻3）若收斂，則迭代完畢，增量步結(jié)束。否則將計算塑性參數(shù)的增量：加_育一卄；（）令（巨）并對塑性應變和內(nèi)變量進一步更新:(19-2：(19-2；(19-2：應+3二Fa)-2闕叫號nFu+n=W+加AAU+1)=3+虻=ptk)+厶嚴)（19-2；（19-2：然后將更新的變量返回屈服條件進行檢査，整個過程將重復直至收斂為止。這就是埠步中應力更新的過程。19.3ABAQUS用戶材料子程序用戶材料子程序（user-definedmateri

10、almechanicalbehavior簡稱UMAT）通迥ABAQUS主求解程序的接口實現(xiàn)與ABAQUS的數(shù)據(jù)交流。在輸入文件中，使用關鋌*USERMATERIAL表示定義用戶材料屬性口19.3.1子程序概況與接口UMAT子程序具有強大的功能使用UMAT子程序有以下幾點需注意:1）可以定義材料的本構(gòu)關系，使用ABAQUS材料庫中沒有包含的材料進行計擴充程序功能。2）幾乎可以用于力學行為分析的任何分析過程，幾乎可以把用戶材料屬性賦ABAQUS中的任何單元。3）必須在UMAT中提供材料本構(gòu)模型的雅可比（Jacobian）矩陣，即應力增量應變增量的變化率。4）可以和用戶子程序“USDFIJT聯(lián)合使用

11、，通過USDFUT重新定義單元一物質(zhì)點上傳遞到UMAT中場變量的數(shù)值。由于主程序與UMAT之間存在數(shù)據(jù)傳遞*甚至共用一些變量，因此必須遵守有TTMAT園爭宜構(gòu)才TTMAT中譽田的恥畐肚方枠#乩豐劉審暫滴直施倉嗎.3NDI,NSHR,WTENS,NSTATV,PROPS.NPROPS,COORDS,IWT.PNEDT,LENT,DFGRDO,HFGRDl,NOEL,NET,LAYER,KSPT,KSTEP,KING)INCLUDE1ABA_PARAM.INCCHARACTER*80CMNAMEDIMENSIONSTRESS(NTENS),STATEV(NSTATV),IXJSDDE(WTENS,

12、tTTEHS),MSDDT(1WNS),DRPLDE(NTENS)fSTRAN(MENS),DSTRAN(WIENS),TIME(2),PREDEF(1)TDPRED，PROPS(NPROPS),COORDS(3),DROT(3,3),LFGRDO(3.3),HFGRDl(3f3)usercoding七0defineDDSDDE,STRESS,STATEV,SSE,SPD、SCDand,ifnecessary,RPL,DDSDDT,DRHBE,imOT.PNEBJCTRTURHUMAT中的應力矩陣、應變矩陣以及矩陣DDSDDE、DDSDDT.DRPLDE是直接分量存儲在前，剪切分量存儲在后。直

13、接分量有個DNI個，剪切分量有個。各分量之間的順序根據(jù)單元自由度的不同有一些差異，所以編寫UMAT時到所使用單元的類別。下面對UMAT中用到的一些變量進行說明。DDSDDE(NTENS,NTENS):是一個NTEN5維的方陣，稱作雅可比矩P尬是應為的增量，虛是應變的增量，DDSDDE(/,J)表示增量2時第J個應變分量的改變引起的第個應力分量的變化通常雅可比是一個對稱除非在“詵USERMATERIAL”語句中加入了UNSYMM參數(shù)乜STRESS(NTENS)：應力張量矩陣，對應NDI個直接分量和NSHR個剪切在增量步的開始，應力張量矩陣中的數(shù)值通過UMAT和主程序之間的接口申UMAT中；在增量

14、步的結(jié)束，UMAT將對應力張量矩陣更新。對于包含剛體轉(zhuǎn)Z限應變問題，一個增量步調(diào)用UMAT之前就已經(jīng)對應力張量的進行了剛體轉(zhuǎn)動,在UMAT中只需處理應力張量的共旋部分。UMAT中應力張量的度量為柯西(應力。STATEV(NSTATEV):用于存儲狀態(tài)變量的矩陣，在增量步開始時將數(shù)勺到UMAT中。也可在于程序USDFLD或UEXPAN中先更新數(shù)據(jù)，然后在增量戈時將更新后的數(shù)據(jù)傳遞到UMAT中。在增量步結(jié)束時必須更新狀態(tài)變量矩陣匚據(jù)。和應力張量矩陣不同的是：對于有限應變問題，除了材料本構(gòu)行為引起的數(shù)寸以外，狀態(tài)變量矩陣中的任何矢量或者張量都必須通過旋轉(zhuǎn)來考慮材料的剛體運2NSTATEV：狀態(tài)變量矩

15、陣的維數(shù)，尊于關鍵字rDEPVAR中定義的數(shù)STANTS”常數(shù)設定的值PROPS(NPROPS)：材料常數(shù)矩陣，矩陣中元素的數(shù)值對應于關鍵字u*US：MATERIAL下面的數(shù)據(jù)行。SSE.SPD.SCD：分別是義每一增量步的彈性應變能、塑性耗散和蠕變耗散。們對計算結(jié)果沒有影響，僅僅作為能量輸出。其他變量：STRAN(NTENS)：應變矩陣。DSTRAN(NTENS)：應變增量矩陣。DTIME：增量步的時間增量。NDI：直接應力分量的個數(shù)。NSHR；剪切應力分量的個數(shù)。NTENS：總應力分量的個數(shù)，NTENS=NDI+NSHRo使用UMAT時需要注意單元的沙漏控制剛度和橫向剪切別度。通常減縮積分

16、単的沙漏控制剛度和板、殼、梁單元的橫向剪切剛度是通過材料屬性中的彈性性質(zhì)定的。這些剛度基于材料初始剪切模量的值，通常在材料定義中通過ELASTIC，項定義。但是使用UMAT時，ABAQUS對程序輸入文件進行預處理的時候得不到變模量的數(shù)值，所以這時候用戶必須使用JHOURGLASSSTIFFNESS選項來定5!有沙漏模式的單元的沙漏控制剛度，使用*TRANSVERSESHEARSTIFFNESS項來定義板、殼、梁單元的模向剪切剛度。19.3.2編程基于上面所述的率相關材料公式和應力更新算法參照ABAQUS用戶材料子稍的接口規(guī)范*進行UMAT的編程。有限元模擬結(jié)果將在下一節(jié)給岀，在最后一節(jié)中給出了

17、相應的程序源代碼。由于UMAT在單元的積分點上調(diào)用，增量步開始時，主程序路徑將通過UM的接口進入UMAX單元當前積分點必要變量的初始值將隨之傳遞給UMAT的相応量。在UMAT結(jié)束時，變量的更新值將通過接口返回主程序。整個UMAT的流磚圖19-2所示。一共有8個材料常數(shù)需要給定，并申請了一個13維的狀態(tài)變量矩陣，它們表示物理含義如表19-1所示。19-1UMAT材料常數(shù)PROPS12345678物理性質(zhì)楊氏模量泊松比塑性耗散比ABnCmSTATEV1671213變量意義彈性應變塑性應變等效鱉性應變512512$11ra19-2UMATJt稈囹19.4SHPB實驗的有限元模擬F面將建立SHPB實驗

18、的有限元模型.井把前面所建立的UMAT接入ABAQUS/Standard進行有PH元棋擬，進行有眼元模擬的目的不繪單地為了禪現(xiàn)SHPB實驗的過程.同時也足為了對選擇的ttffim和建立的IJMAT進行評價，19.4.1分離式Hopkhison壓桿(SHPB)奕監(jiān)分離式Hopkinson壓桿(splitHopkirvonpressurebar.面稱SHPB)實驗是從經(jīng)典blopkinson實臉系砒之上發(fā)展而來的一種實驗技術.用來測供材料的動念應力應變行為。該實驗技術的理論慕礎是堆應力波理論.通過測傲兩根壓桿上的應變來推尋試件上的應力-應變關系。分離式Hopkizn壓桿實驗的示jg圖見圖193。I

19、U19-35/：HcvkinwimiF*W19.4.2有限元連模有限元模熨主要杲倉照前血介紹的SHPB實驗裝程.if過載荷、邊界條件爭的定義在有限元中棋擬SHPB實臉的殲境.盡螢在較少的機時耗費F達到更髙的榊度.L援型的尙億與有限無岡站為f不使棋熨過干龐大對箱型進行了一些簡化。首先.改變?nèi)肓U和出力桿的尺寸.長度由K(來的MMUrwn減小為1(MWmm,直位堆加到25nun,試件的長度和點徑分別變化為22mmfU18mmv這樣不僅優(yōu)化了網(wǎng)格的.還成倍地咸小了復中的規(guī)模其帶來的負面影響就是試件能達到的應變將降低另外.由于撞擊桿僅僅起到產(chǎn)生應力脈神的作用莊數(shù)值模型中沒必姿有慮曲擊桿.取代的方法足點

20、接在入力桿的輸入端魏加均布的應力關沖著慮到實夠裝世的對稱性.也作了一觀簡化。菸個實驗裝匱以及戟荷弄柿是關于桿的中心線輔對稱的所以可以使用軸對稱血元進行二維分析。刃外也建遼了四分之一橫嚴面的維模型作為補充一蠢軸對休模型和三維棋熨分別如圖194、圖195所示。在棋乜中對試件以及入力桿、出力杼和試件接粧的部分進行了局85阿格加密.這樣的網(wǎng)格IM分可以取得比校經(jīng)濟的結(jié)果。團194元禪51fflN5堆有阪元鍥型$11單元類型上，選擇一階常規(guī)單元，由于沒有使用減縮積分單元，所以使用時無須指定單元的沙漏控制剛度。最后的模型中，二維網(wǎng)格單元的總數(shù)為122C模型網(wǎng)格的單元總數(shù)為17160.關于單元的詳細信息如表

21、19-2所示。表1頭2模型信息-型模L”X-單-一數(shù)個元單總二籐型桿力入-一Z-I一5-n41-1-31桿力出三養(yǎng)型一桿力入r勺r-5007-一711一-桿力出-750在二維模型局部網(wǎng)格存在疏密連接的部位，一個單元邊要同時和兩個單元逍在通常的有限元網(wǎng)格中是不易實現(xiàn)的。本例在這里使用了ABAQU5的多點約(multi-pointsconstrain,簡稱MFC)來解決，線性多點約束方式如圖19-6所示口圖中點p使用線性多點約束后，其節(jié)點自由度均由旁邊的節(jié)點撫和方線性到。所以使用多點約束方式可以很好地連接模型中網(wǎng)格疏密不同的部位，劃分出細的網(wǎng)格。2材料定義入力桿和岀力桿使用線彈性材料，彈性模量和泊

22、松比分別為200GPa和0二為7,85xio3試件采用用戶在UMAT中的自定義材料*材料參數(shù)如表表193試件的材料定義性質(zhì)密度/(kg/m5)楊氏模量/MPa泊松比Johnson-Cook模型參數(shù)A/MPaB/MPanCm數(shù)值27*068.0X1030.3366x562108.8530.2380.0290J3邊界條件為了保證SHPB實驗的要求，在二維模型和三維模型中均施加了必要的邊界條在對稱軸或?qū)ΨQ面上施加了對稱性邊界條件，同時保證壓桿和試件可以沿軸線方向E無約束地運動。壓桿和試件之間的接觸為硬接觸，光滑無摩擦*為了確定輸入應力脈沖的時間，進行了簡單的計算.彈性材料中縱波波速的計爭式為$11其

23、中E為材料彈性模量，卩為材料密度。由此可以計算輸入應力波在壓桿中的傳*度為Cd=5048m/so要求在入力桿應力波的輸入端不能出現(xiàn)入射波和反射波的重疊，也就是說在輸丿力脈沖的時間內(nèi)，應力波的傳播距離不應超過兩倍的桿長，即TsAtJ4=TEtef*/Cr.上茲丄urznixL-tAstiz上八宀ZnI+rfIAtJ4=TEtef*/Cr.上茲丄ur$11圖1916四種應變率下的應變-時間曲線70伽頂ss應應應應二二J-J:-80604020-04mm0.00000.0001O.OOO20,00030.000400005圖19J7四種應變率下的應力時間曲線生，引起單元溫度的變化，但是相鄰單元之間沒

24、有熱量的傳遞。在應變率250s下件最終的溫度場如圖19-21所示。取塑性耗散比為0.95,即95%的塑性耗散將轉(zhuǎn)化為熱量。這一比例對于大多數(shù)屬是合適的.從圖19.21中看出，早期應力波沒有到達試件時.試件上的溫度沒有:$11T應變率100應變寧250圖19-18四種應變率下的應力-應變曲線Edw、R圏圖19-19四種設計應變率下的真實應變率-時間曲線將考慮溫度效應時單元的應力應變與前面沒有考慮溫度效應的情況作一下比件表面中間單元的溫度曲線如圖1922所示。曲線中發(fā)生塑性變形后，單元的溫上升。右祖詹隔毘響卞單三萌廟*審亦兩由站緯伽園IQg霸條肚閆酚苗木窘$11z1a*鍥sJu魚篆：SWM一z-6

25、-s#但qsY2xofss*0仝君?-*辺&Hs-E:*tiqe(設I8-$11sVIS*整連驢裂ST嘗1VM355SE寫猱J國$&班七SH徉另Bgm亠89U?_EOWW鬼壬譽整房藝糊廻注藝ww9sfe潑生E4姿題如犖aJYlte5-$1128$11$11270耳000000.0001o6oo200003o.6oo4o.doos時間蟲圖19-22單元溫度的變化$11$110.04-003-制圏0020.01aoo0,0000o.oooi$11（町應變時閭曲錢$11$11時同滄800800524urr1-r-t1*r-OQOO0850.010001$n.(i200.02500300035(XH0

26、OXHS800800524800800524圖19-24專球催度效應與不(Sfi度效應的比較19.4.4三維動態(tài)分析本例IS進行的二維動態(tài)分析也相當成功試件往康擬過程中各處應力都比較均勻議大應力和廉小應力的建別不超過IMPao圖19-25為試件的Mises應力分布。800800524800800524800800524變形過穆中單元的應力、應變以及應變來歷史桶與二維的悄況非常吻合如圖1926所示。圖1927所示為應力應變曲線.這是對編寫的UMAT用于三維實體單元的-個笑證。UMATM然見從-tftJohn-son-Cook施熨中建立起來的但是在UMAT中率相關也性公式將它擴展列了三維情況.所以

27、能用于三維實體單元也在意料之中。雖怡取得的結(jié)果是一致的.但是在計算時間上卻冇很大差別在個人計算機上使用軸對稱-.維網(wǎng)格完成一次計篦只需要10min.然后完成-次濰網(wǎng)恪的計韓需翌3h。圖1926三維有限元模擬的應力、應變率歷史圖19-27三維有限元模擬中的應力-應變曲線相比而言使用軸對稱二維模型要經(jīng)濟得多。195UMAT的Fortran程序19.5.1UMATSUBROUTINEUMAT(STRESS,STATEV,DDSME,SSE,SEO,9CD,RPL,EDSDOT,DRPLEffi.EREWT,STRAN,DSTFiAN,524INCLUDE7ABAPARAM-INCCHARACTERB

28、OMATERLDIMENSIONSTBESS(WIE陽)，STAIEV旳/TV.DDSXe(NTENS,!nENS),DDSDET(inENS,DRElDE(iniMS人STRAMdnENShDSTRMinENShTIMEdhPFEWdhEREDd),PK5(NPROPS.COORDS.lWr(32),DFDO(3,3.DFGRD1(33)DIMEWSIONLAS，EELAS(6),F1QW524524PARAMETER(ONE=l.OIX),TTO=2,0D0tTHREE=3.ODO,SIX=+ODO,HALF=0.5d0)DATANEKTON,TOLER/40覧,96/CCUMATFORJ

29、OHNSON-COCKMODELCCPROPSYANGSMCDULUSCPFCPSPOISSONRATIOCFROPS(3)INELASTICHEATFRACTIONPARAMETOFJOHNSOJUDCKMCQEL：524524CPROPSATOC o 1-5 h zCPP0PS(5)七BCPROPS-nCFR0PS7)-CCPROPS(0)mIF(NDLNE.3)THENWRTrE(6,l)FQEMAT(/r3QXSl粉ERROR-THISUMATMAYCMLYHEUSEDFOR1ELEMEKTSWTIHTHREEDIRECTSTRESSCCMPONEJnS，)ENDIFCCELASTIC

30、PROPERTIESCEtW=FRQPS(l)ENUPE?OPS(2)EG3=THREE*EGELAM二(EBUTO-EG2)/THREECCELASTICSTIFFNESSCDO20Kl=l,NTEKSDO10K2=1,NTENSDDSM(K23KI)=0R10COHTINUE20OGMTINUECDO40Kl=l.NDIDO30K2=1?NDIDDSDDE(K2,K1)二ELAM30COKTINUEEDSCCEtKl,KI)=+ELAM40CONTINUEDO50KI=NDI+l.NTENS50CONTINUECCCALCULMESTRESSFROMELASTICSTRAINSCDO70K

31、l=lsNTENSDO60K2=1,NTEHSSTRESS(K2)=STRESS(K2)+DDSEC(K2,Kl)并DSIRAN(Kl)60CONTINUE70CONTINUECCREOOVERELASTICANDPLASTICSTRAINSCDO80Kll.NTENSLRS(Kl)-STATEV(Kl)+DSTRfiN(Kl)EPLAS(Kl)=STATEV(K1+WIENS)80(XMINUEEQPLRS=STATEVC1+2供NTENS)CCALCULATEMISESSTRESSDO90KI=NDI+l.NTENSSHISES=SMISES+SIX*STRESS(K1)*STRESS(K

32、1)90OCSTINUESHISES=SQRTfSMlSES/TWO)CCCALLUSEPHAK3SUBROUTINE,GETHARDENINGRATEANDYIELDSTRESSCCCALLUSERHMD(SYIELO,HMD,EQPLAS,H?0PS)CDEIERHINEIFACTIVELYYIELDINGCIF(SMISES.GT_(1.0+TOLER)*SYIELO)THENcCMATERIALRESPONSEISPLASTIC,DETERMINEFLOWDIRECTIONCSHDRO=(STRESS(l)+STRESS(2)+STRESS(3)/THREEONESY=ONE/SMIS

33、ESDO110KI=1,NDIFTjCW(KI)=ONESY*(STRESS(K1)-SHYDRO)110CONTOJUEDO120KI-NDI+1,NTENSFLOW(Kl)=STRESS(K1)*ONESY120CONTINUElaREADPARAMETERSCFJOHNS(潞COOKFCA=TOPS(4)B=PROPS(5)EN=PR0PS(6)C=PR0PS(7)EM=PR0PS(6)CCNEWTONITERATIONCSyiELD=SYIELODECPL=(SMISES-SYIEED)/ES3DSTRES=IOLER*SYIELO/EG3DEQMIN=HALF*nriME*EXP(1

34、.OD-4/C)DO130KEWTON=1FNEWTONTTfTtviriThlTtvirnrnrETtr524HARD1=HM3TVRL+SYIELX/EEQPLSYIELD=SYIELD*TVP1TTfTtviriThlTtvirnrnrETtr524TTfTtviriThlTtvirnrnrETtr524RHS-SMISES一EG3#FL-SYIELDTTfTtviriThlTtvirnrnrETtr524EEQPL二DEQRL+RHS/(EG3+HARD1)IF(ABS(RHS/EG3).LE.DSTRES)GOTO140130CONTINUEWRITE(6,2)NEWTON2FORM

35、AT/,30X,*WARNINGPLASTICITYALGORITHMDIDNOT1CONVERTAFTER13,ITERATICSS)140CONTINUEEFFHRD=EG3*HARD1/(EG3+HARD1)CCCALCULATESTRESSANDUPDATESTRAINSDO150KI=1#NDISTRESS(Kl)=FLOW(K1)*SYIELD+SHYDROEPLfiS(Kl)=EPLRS(K1)+IHREE餐FLO鍬K1”DEQPL/TO150CONTINUEDO160K1=NDI+1,NTENSSTRESS(Kl)=FLOW(SYIELO+SYIELD)/IWORPL=PROPS星SPD/DTIMECCJACOBIANEFFG=EOSYIELD/31ISESEETG2=TH(E