+用戶材料子程序?qū)嵗?Johnson-Cook+金屬本構(gòu)模型

1、ABAQUS軟件2003年度用戶年會(huì)論文集中國(guó)CAE聯(lián)盟 工 www.CAE 資料收集 M版權(quán)歸屆原作者ABAQUS/Standard用戶材料子程序?qū)嵗齁ohnson-Cook金屬本構(gòu)模型盧劍鋒莊茁"張帆清華大學(xué)工程力學(xué)系北京100084摘要：用戶材料子程序是ABAQUS提供給用戶定義Fl己的材料屬性的Fortran程序接口，使用 戶能使用ABAQUS材料庫(kù)中沒(méi)有定義的材料模型。ABAQUS中H有的Johnson-Cook模熨只能應(yīng)用丁顯式ABAQUS/ExpUcit程序中，而我們 希與能在隱式ABAQUS/Standard程序屮更椿確的實(shí)現(xiàn)本構(gòu)積分，而且應(yīng)用Jolmson-Cook

2、模型 的修正形式。這就需更通過(guò)ABAQUS/Standard的用戶材料子程序UMAT編程實(shí)現(xiàn)。在UMAT 編程中使用了率相關(guān)蜩性理論以及完全隱式的W力更新算法。1 Jolmson-Cook強(qiáng)化模型簡(jiǎn)介Jolmson-Cook (JC)模型用來(lái)模擬高應(yīng)變率下的金屬材料。JC強(qiáng)化模型農(nóng)示為三項(xiàng)的乘積, 分別反映了應(yīng)變硬化，應(yīng)變率硬化和溫度軟化。這里使用JC模熨的修1E形式：a =（A+Bsn） 1 + Cln（l +彳（1 一廣）并使參考應(yīng)變率金=1,這樣公式中的力即為材料的莎態(tài)屈服應(yīng)力。公式中包含A,E,n,C,m五個(gè)參 數(shù)，需要通過(guò)實(shí)驗(yàn)來(lái)確定。2 ABAQUS用戶材料子程序用戶材料子程序(Us

3、eT-defined Material Mechaiucal Behavior,簡(jiǎn)稱 UM AT)通過(guò)與ABAQUS ±求解程序的接口實(shí)現(xiàn)與ABAQUS的數(shù)據(jù)交流。在輸入文件中，使用關(guān)鍵字'“USERMATERIAL"表示定義用戶材料屬性。子程序概況與接口UMAT子程序具有強(qiáng)大的功能，使用UMAT子程序：可以定義材料的木構(gòu)關(guān)系，使用ABAQUS材料咋中沒(méi)有包含的材料進(jìn)行計(jì)如擴(kuò)充程序功能。-1-BackABAQUS軟件2003年度用戶年會(huì)論文集(2) 兒乎可以用丁力學(xué)行為分析的任何分析過(guò)程，幾乎可以把用戶材料屬性賦予ABAQUS中 的任何單元：(3) 必須在UMAT中

4、提供材料本構(gòu)模型的雅可比(Jacobian)矩陣，即應(yīng)力增址對(duì)應(yīng)變?cè)隽?的變化率。(4) 可以和用戶子程序'PSDFLD聯(lián)A使用，通過(guò)rjSDFLD'亜新定義單元每一物質(zhì)點(diǎn)上傳 遞到UMAT中場(chǎng)變最的數(shù)值。由于主程序與UMAT之間存在數(shù)據(jù)傳遞，其至共用一些變帚:，因此必須遵守有關(guān)UMAT的 書(shū)寫(xiě)格式，UMAT中常用的變応在文件開(kāi)頭予以泄義，通常格式為：SUBROUTINE UMAT(STRESS,ST ATEV,DDSDDESEZSPD,SCD,1 RPL,DDSDDT,DRPLDE,DRPLDTZ2 STRAN,DSTRAN,TIlvIE,DTIME,TElvIPzDTElv

5、lP,PIxEDEF,DPRED/CMN AME,3NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDSQROT,PNEWDT,4 CELENT,DFGRDOZDFGRD1ZNOELZNPT,L AYER,KSPTZKSTEP,IONC)CINCLUDE 'ABA-PARAMINC*CCHARACTER*80 CMNAMEDIMENSION STRESS(NTENS)STATEV(NSTATV),1 DDSDDE(NTENS,NTENS ),DDSDDT(NTENS ),DRPLD 玖 NTENS),2 STRAN(NTENS)ZDSTCAN(NTENS)

6、ZTIME(2),PREDEF( 1 ),DPRED( 1),3PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRDO(3,3),DFGRD1(3,3)user coding to define DDSDDE, STRESS, STATEV, SSE, SPD, SCDand, if necessary, RPLZ DDSDDT, DRPLDE, DRPLDT, PNEWDTRETURNENDUMAT中的應(yīng)力矩陣、應(yīng)變矩陣以及矩陣DDSQDE, DDSDDT , DRPLDE等,都是H接分 最存儲(chǔ)在前，剪切分帚儲(chǔ)在后。W接分最有ND/個(gè)，剪切分最冇臨飯個(gè)。各分最之間的

7、順序 根據(jù)單元H由度的不同有一些差異，所以編寫(xiě)UMAT時(shí)?？紤]到所使用單元的類別。下面對(duì) UMAT中用到的一些變雖進(jìn)行說(shuō)明：DDSDDE (NTENS, NTENS)是一個(gè)NTENS維的方陣，稱作雅可比矩陣，沁IE, <!是應(yīng)力的増鼠，山是應(yīng)變的增暈, DDSDDE表示增帚:步結(jié)束時(shí)第丿個(gè)應(yīng)變分帚:的改變引起的第I個(gè)應(yīng)力分帚的變化。通常雅可 比是一個(gè)對(duì)稱矩陣，除非在"USER MATERIAL"語(yǔ)句中加入:T 'TJNSUNT參數(shù)。STRESS (NTENS)應(yīng)力張帚矩陣，對(duì)應(yīng)NZV個(gè)氏接分帚和冋飯個(gè)剪切分量。在增帚步的開(kāi)始，應(yīng)力張最矩陣 屮的數(shù)值通過(guò)UMAT

8、和主程序Z間的接口傳遞到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)力張帚:的 度最為柯西(真實(shí))應(yīng)力。STATEV(NSTATEV)用丁存儲(chǔ)狀態(tài)變最的矩陣，在増暈步開(kāi)始時(shí)將數(shù)值傳遞到UMAT中。也可在子程序USDFLD 或UEXPAN中先更新數(shù)據(jù)，然后增最步開(kāi)始時(shí)將更新后的數(shù)據(jù)傳遞到UMAT中。在增最步的結(jié) 朿必須蟲(chóng)新?tīng)顟B(tài)變屋矩陣中的數(shù)據(jù)°和應(yīng)力張最矩陣不同的是：對(duì)丁有限應(yīng)變問(wèn)題，除了材料木構(gòu)行為引起的數(shù)據(jù)更新以外，狀

9、態(tài)變量矩陣屮的任何矢戢或者張量都必須通過(guò)旋轉(zhuǎn)來(lái)考慮材料的剛體運(yùn)動(dòng)。狀態(tài)變砒矩陣的維數(shù)，等丁關(guān)鍵字'“DEPVAR'定義的數(shù)值。狀態(tài)變量矩陣的維數(shù)通過(guò) ABAQUS輸入文件中的關(guān)鍵字'“DEPVAR，定義，關(guān)鍵字下面數(shù)據(jù)行的數(shù)值即為狀態(tài)變帚矩陣的 維數(shù)。材料常數(shù)的個(gè)數(shù)，等丁關(guān)鍵'j"uSERMATERIAL"中"CONSTANTS'常數(shù)設(shè)定的值。PROPS(NPROPS)材料常數(shù)矩陣，矩陣中元素的數(shù)值對(duì)應(yīng)丁咲鍵字'“USERMATERIAL下而的數(shù)據(jù)行。SSE , SPD, SCD分別定義每一增最步的彈性應(yīng)變能，型性耗

10、散和疑變耗散。它們對(duì)計(jì)算結(jié)果沒(méi)有影響，僅僅 3BackABAQUS軟件2003年度用丿'*年會(huì)論文集作為能最輸出。其他變最：STRAN NTENS):應(yīng)變矩陣；DSTRAM(NTENS)： hY變?cè)鲎罹仃?；DTIME：増帚步的時(shí)間增帚:：NDI： X接應(yīng)力分量的個(gè)數(shù)：NSHR：剪切應(yīng)力分最的個(gè)數(shù)；NTENS：總應(yīng)力分斌的個(gè)數(shù)，NTENS = NDI + NSHR。使用UMAT時(shí)礙要注意單元的沙漏控制剛皮和橫向剪切剛皮。通常減縮積分單元的沙漏控 制剛度和板、売、梁?jiǎn)卧臋M向剪切剛度是通過(guò)材料屈性中的彈性性質(zhì)定義的。這些剛皮基丁材 料初始剪切模量的值，通常在材料定義小通過(guò)"ELA

11、STIU選項(xiàng)定義。但是使用UMAT的時(shí)候， ABAQUS對(duì)程序輸入文件進(jìn)行預(yù)處理的時(shí)候得不到剪切模鼠的數(shù)值。所以這時(shí)候用戶必須使用 "HOURGLASS STIFFNESS"選項(xiàng)來(lái)定義具有沙漏模式的單元的沙漏控制剛度，使用 TRANSVERSE SHEAR STIFFNESS"選項(xiàng)來(lái)定義板、殼、梁?jiǎn)卧臋M向剪切剛度。編程基丁上而所述的率相關(guān)材料公式和應(yīng)力更新算法，參照ABAQUS用戶材料子程序的接口規(guī) 范，進(jìn)行UMAT的編程。有限元模擬結(jié)果將在下一節(jié)給出，在垠后一節(jié)中還給出了相應(yīng)的程序 源代碼。由于UMAT在單.元的積分點(diǎn)上調(diào)用，增最步開(kāi)始時(shí)，主程序路徑將通過(guò)UM

12、AT的接口進(jìn)入 UMAT,單元當(dāng)詢積分點(diǎn)必要變量的初始值將隨Z傳遞給UMAT的相應(yīng)變量。在UMAT結(jié)束時(shí), 變最的更新值將通過(guò)接口返回主程序。幣個(gè)UMAT的流程如圖1所示。一共有8個(gè)材料常數(shù)需要給定，并申諳一個(gè)13維的狀態(tài)變最矩陣，它們表不的物理含義如表 1所示。下一步將使用建立的UMAT結(jié)介ABAQUS/Standard進(jìn)行袍布金森沖擊桿(SHPB)實(shí)驗(yàn)的有 限元模擬，并對(duì)結(jié)果進(jìn)行比較。圖1UMAT流程圖表1UMAT材料常數(shù)PROPS12345678物理性質(zhì)楊氏模量泊松比塑性耗散比ABnCMSTATEV1-67-1213變最意義彈性應(yīng)變塑性應(yīng)變等效塑性應(yīng)變3SHPB實(shí)驗(yàn)的有限元模擬模型的簡(jiǎn)化

13、與有限元網(wǎng)格為了不使模烈過(guò)于龐大，對(duì)模型進(jìn)行了一些簡(jiǎn)化。首先，改變?nèi)肓U和出力桿的尺寸，長(zhǎng)度 由原來(lái)的3040nmi減小為1 OOOmm, 徑增加到25mm,試件的長(zhǎng)度和H徑也分別變化為22mm 和18nmio這樣不僅優(yōu)化了網(wǎng)格的質(zhì)暈，還成倍地減小了模型的規(guī)模，其帶來(lái)的負(fù)面影響就是試 件能達(dá)到的應(yīng)變將降低。另外，由于撞擊桿僅僅起到產(chǎn)生應(yīng)力脈沖的作用，在數(shù)值模型中沒(méi)必耍 考慮撞擊桿，取代的方法是IT接在入力桿的輸入端施加均布的應(yīng)力脈沖?？紤]到實(shí)驗(yàn)裝置的對(duì)稱性，也做了一些簡(jiǎn)化。整個(gè)實(shí)驗(yàn)裝置以及載荷等都是關(guān)丁桿的中心線 軸對(duì)稱的，所以可以使用軸對(duì)稱單元進(jìn)行二維分析。二維軸對(duì)稱模型如圖2所示。在模型中

14、，對(duì)試件以及入力桿，出力桿和試件接觸的部分進(jìn)行 了局部網(wǎng)格加密，這樣的網(wǎng)格劃分可以取得比較經(jīng)濟(jì)的結(jié)果。-#-BackABAQUS軟件2003年度用戶年會(huì)論文集入力桿出力桿試件圖2二維軸對(duì)稱有限元模型表2模型信息模熨尺寸inm（ xL）單元類些單元個(gè)數(shù)總節(jié)點(diǎn)數(shù)總單元數(shù)二維模型入力桿25x1000CAX453014751220試件18x22CAX4160出力桿25x1000CAX4530材料定義入力桿和出力桿使用線彈性材料，彈件模帚:和泊松比分別為200GPa和0.3,密度為7.85x103 kg/m3o試件采用用戶在UMAT中門(mén)泄義材料，材料參數(shù)如表3所示，其屮Johnson-Cook模型 中參

15、數(shù)的數(shù)值來(lái)源丁前而的數(shù)值擬介程序。衣3試件的材料建義性質(zhì)密度Kg/irf楊氏模 量 MPa泊松比Jolmson-Cook模型參數(shù)AMPaBMPanCM數(shù)2.7x10368.0x1030.3366.562108.8530.2380.0290.5邊界條件為了保證SHPB實(shí)驗(yàn)的要求，在二維模型中施加了必要的邊界條件。在對(duì)稱軸上施加了對(duì)稱 性邊界條件，同時(shí)保證床桿和試件可以沿軸線方向IH由無(wú)約束的運(yùn)動(dòng)。床桿和試件之間的接觸為 硬接觸，光滑無(wú)摩擦。為了確定輸入應(yīng)力脈沖的時(shí)間，進(jìn)行了簡(jiǎn)單的計(jì)算。彈性材料中縱波波速的計(jì)算公式為：q =其中E為材料彈性模量，Q為材料密度。由此可以計(jì)算輸入應(yīng)力波在壓桿中的傳播速

16、度為Cd = 5048 m/So耍求在入力桿應(yīng)力波的輸入端不能出現(xiàn)入射波和反射波的重廉，也就是說(shuō)在輸入應(yīng)力脈沖的時(shí)間內(nèi)，應(yīng)力波的傳播距離不應(yīng)超過(guò)兩倍的桿長(zhǎng),即:Tt < =?« 4.0x10 (s) Cd 5048根據(jù)這一佔(zhàn)計(jì)，選擇輸入應(yīng)力脈沖的持續(xù)時(shí)間7； =2.0x10'4s,上升時(shí)間fr =3.0x10-5so經(jīng)過(guò)若干次試算，對(duì)輸入應(yīng)力脈沖的波形進(jìn)行適當(dāng)?shù)恼{(diào)格，使試件中產(chǎn)生較均勻的應(yīng)變率。最后輸入應(yīng)力脈沖的波形如圖3所示：圖3輸入應(yīng)力脈沖為了饑定用戰(zhàn)步的最大時(shí)間步長(zhǎng)，品耍先簡(jiǎn)單計(jì)算一卜單元的穩(wěn)定極限?；粋€(gè)單尤的估 算，穩(wěn)定極限可以用單元特征長(zhǎng)度厶'和材料

17、波速Q(mào)定義如下：斥桿單元的特征單元長(zhǎng)度r«10nmi.由此可以計(jì)算出應(yīng)力波金斥桿傳遞的穩(wěn)定極限為S宀0x10%)將它作為ABAQUS門(mén)動(dòng)增磺控制里而的最大時(shí)間步匕 二維動(dòng)態(tài)分析我們所對(duì)照的SHPB實(shí)驗(yàn)止是屈于這一情況，所以可以將ABAQUS/Standard結(jié)合UMAT 進(jìn)行有限元模擬的結(jié)果和實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比。下面是應(yīng)變率250 s-i下的動(dòng)態(tài)模擬過(guò)程。在時(shí)間f = 1.98x107($)左右，應(yīng)力波前沿到達(dá)試件，這一時(shí)間和前面使用彈性波波速計(jì)算的傳播時(shí)間是相同的，此前試件上的Mises應(yīng)力兒乎為零，如圖4所示。S, Mises|Ave Cert ： fil7e*05 66 60S

18、514e*05 3 6 心05 211etO5 060&-05 08be*04 571e*04 057e+04 543e*CI4 029rr*04 514e*04 801C-25圖4應(yīng)力波前沿到達(dá)試件時(shí)的Mises應(yīng)力(t=1.98xl0-4 s)在時(shí)間心3.0x107(s)，試件經(jīng)過(guò)應(yīng)力波的上升時(shí)間后達(dá)到穩(wěn)定變形的狀態(tài)，一部分入射波 反射冋入力桿，一部分應(yīng)力波經(jīng)過(guò)試件進(jìn)入出力桿，試件各點(diǎn)的變形都很均勻，如圖5(a)所示。 在圖5(b)試件的放犬圖上可以看出，各點(diǎn)Mises應(yīng)力相差不超過(guò)IMPa,這個(gè)精度是相當(dāng)可靠的。(Ave. Cr it.: 75)H. 277-05*!. L70e

19、-*-051.06*et05 -*9.574e*04 _ M 一 f383H04 5 32*04M.255e*04 >3. L91e+04 2 12 8*04 *L04e+04 9. 595e 11(a)全局視圖S, Mises(Ave Cir 丄t ：75 % |M.277e-tO5M.276etO5 M.275e-K)5 M.275e-tO5 r "274bOS r 1.273et-O5 k M.273e-K)5 k M.272e*O5 一 *1.272et-O5 i-1.271e-i-O5 *1.27Oe-t-O5 tl.27OetO5 M.2 69e-t-O5(b)試件

20、的放大視圖圖5試件經(jīng)歷均勻變形時(shí)的Mises應(yīng)力(t=3.0xl0-4 s)經(jīng)過(guò)穩(wěn)定變形階段后，反射波和傳遞波分別向入力桿和出力桿擴(kuò)散，試件上Mises應(yīng)力逐漸 減小到較低的水平，試件開(kāi)始經(jīng)歷卸載，如圖6所示。圖中Mises應(yīng)力云圖的單位為KPa。1. 189*05 1.09D-fD5 9.9X3O047.938004 +6950d 3.962e*04 H.975e*01 3.98e*04 t3.000«t04 2.O12e*04 1.02 5e>D4 3.7n4etD2圖6應(yīng)力波消退后試件時(shí)的Mises應(yīng)力(t=4.2xl0-4s)實(shí)際上有限元模擬的應(yīng)力一應(yīng)變曲線和恒定應(yīng)變率

21、下實(shí)驗(yàn)的結(jié)果也能夠很好的吻介。取出試件表面中間的一點(diǎn)，將應(yīng)變率250 s-i和200sJ下ABAQUS有限元模擬的結(jié)果與實(shí)驗(yàn)的結(jié)果對(duì)比 9BackABAQUS軟件2003年度用戶年會(huì)論文集見(jiàn)圖7和圖8。160n140120-100-80-60-20-00.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035 0.040應(yīng)變740-應(yīng)變率實(shí)驗(yàn)值-e- ABAQUS 模擬450400350300250200«1501005000.045-,500圖7應(yīng)力一應(yīng)變曲線的對(duì)比及模擬過(guò)程中真實(shí)應(yīng)變率變化(250 sJ)圖8應(yīng)力一丿應(yīng)變曲線的對(duì)比及模擬過(guò)程中冀

22、實(shí)應(yīng)變率變化(200 S-1)WRITER 6,1)BackABAQUS軟件2003年度用戶年會(huì)論文集WRITER 6,1)BackABAQUS軟件2003年度用戶年會(huì)論文集4J-CUMAT 程序SUBROUTINE UMAT(STRESS,ST ATnEVQDSDDESESPDSCD1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRANZ2HME,DTIME；rEMPQTEMP,PREDEFQPRED,MATERL,NDI,NSHR,NTENS,3NSTATVfROPS,NPROPS,COORDSQROT,PNEWDT；CELENT,4DFGRDO,DFGRD1,NO

23、EL/NPT,KSLAY/KSPT,KSTEP/KINC)CINCLUDE 'ABA-PARAM-INC*CCHARACTER*80 MATERLDIMENSION STRESS(NTENS ),ST ATE V(NST ATV),1 DDSDD E(NTENS,NTENS)QDSDDT(NTENS ),DRPLD 玖 NTENS),2 STRAN(NTENS)ZDSTR AN(NTENS ),HME(2),PREDEF( 1 ),DPRED( 1),3 PROPS(NPROPS),COORDS(3),DROT(3,3),4 DFGRDO(3,3)QFGRD1(3,3)CDIMENSI

24、ON EEL AS( 6)ZEPL AS( 6),FLOW(6)PARAMETER (ONE=1.ODO,TWO=2.0D0zTHREE=3.0D0,SIX=6.0D0, HALF =0.5d0)DATA NEWTON,TOLER/40,1 .D-6/CCCUMAT FOR JOHNSON-COOK MODELCCPROPS( 1)-YANG'S MODULUSCPROPS(2) - POISSON RATIOCPROPS(3) - INELASTIC HEAT FRACTIONCP ARAMETERS OF JOHNSON-COOK MODEL:CPROPS(4) ACPROPS(5

25、)BCPROPS(6)nCPROPS(7) CCPROPS(8) mCCIF (NDI.NE.3) THEN1FORMAT(/,30X；FERROR - THIS UMAT MAY ONLY BE USED FOR 1'ELEMENTS WITH THREE DIRECT STRESS COMPONENTS')ENDIFCC ELASTIC PROPERTIESCEMOD=PROPS(1)ENU=PROPS(2)IF(ENU.GT.0.4999.AND.ENU.LT.0.5001)ENU=0.499EBULK3=EMOD/(ONE-TWCX ENU)EG2=EMOD/(ONE+

26、ENU)EG=EG2/TWOEG3=THREE*EGEL AM=(EBULK3-EG2)/THREECC ELASTIC STIFFNESSCDO 20 Kl=l,MIENSDO 10 K2=1ZNTENSDDSDDE(K2,K1 )=0.010 CONTINUE20 CONTINUECDO 40 K1=1,NDIDO 30 K2=1,NDIDDSDDE(K2,K1)=ELAM30 CONTINUEDDSDDE(K1,K1 )=EG2+EL AM40 CONTINUEDO 50 K1=NDI+1,NTENSWRITER 6,1)BackABAQUS軟件2003年度用戶年會(huì)論文集DDSDDE(K

27、1ZK1)=EG50 CONTINUECCC ALCULATE STRESS FROM EL ASTIC STRAINSCDO 70 K1=1,NTENSDO 60 K2=1,NTENSSTRESS(K2)=STRESS(K2)+DDSDDE(K2,KirDSTRAN(Kl)60 CONTINUE70 CONTINUECC RECOVER ELASTIC AND PLASTIC STRAINSCDO 80 K1=1,NTENSEELAS(K1)=STATEV(K1)+DSTRAN(K1)EPLAS(K1)=STATEV(K1+NTENS)80 CONTINUEEQPL AS=ST ATE V(

28、 1+2*NTENS)CC C ALCUL ATE MISES STRESSCIF(NPROPS.GT,5. AND.PROPS(4).GT.O.O) THENSMISES=(STRESS(1)-STRESS(2)*(STRESS(1>STRESS(2)4-1(STRESS(2STPESS(3)r(STRESS(2STRESS(3) +1(STRESS(3>STRESS( 1 )*(STRESS(3)-STRESS( 1)DO 90 K1=NDI+1,NTENSSMISES=SMISES+SIX*STRESS(K1)*STRESS(K1)90 CONTINUESMISES=SQRT

29、(SMISES/TWO)C #BackABAQUS軟件2003年度用戶年會(huì)論文集CC ALL USERHARD SUBROUTINE, GET HARDENING RATE AND YIELD STRESSCCCALL USERH ARD(S YIELOZH ARD,EQPL AS,PROPS( 4)C DETERMINE IF ACTIVELY YIELDINGCIF (SMISES.GT.(1.0+TOLER)*SYIEL0) THENCCMATERIAL RESPONSE IS PLASTIC, DETERMINE FLOW DIRECTIONCSHWRO=(STRESS( 1 )+ST

30、RESS(2)+STBESS(3)/THREEONESY=ONE/SMISESDO 11OK1=1,NDIFLOW(K1 )=ONES Y*(STRESS(K1 >SHYDRO)110 CONTINUEDO 120 K1=NDI+1,NTENSFLOW(K1 )=SIPESS(K1 )*ONESY120 CONTINUECC READ PARAMETERS OF JOHNSON-COOK MODELCA=PROPS(4)B=PROPS(5)EN=PROPS(6)C=PROPS(7)EM=PROPS(8)CCNEWTON ITERATIONSYIELD=SYIEL0C 15BackABAQ

31、US軟件2003年度用丿'*年會(huì)論文集DEQPL=(SNnSES-SYIELD)/EG3DSTRES=TOLER*S YIEL0/EG3DEQMIN=H ALFDTTMEDOX1. OD-4/C)DO 130 KEWTON=1,NEVVTONDEQPL=MAX(DEQPLzDEQNflN)C ALL USERHARD(SYIELD,HARD,EQPL AS+DEQPL,PROPS(4)tvpylocxdeqfl/dume)TVP1=TVP+ONEHARD1=HARD*TVP1+SYIELD*C/DEQPLSYIELD=S YIELD”TVP 1RHS=SMISES-EG3*DEQPL-

32、S1ELDDEQFL=DEQPL+RHS/(EG3+HARD1)IF(ABS(RHS/EG3) ,LE. DSTRES ) GOTO 140130 CONTINUEVVRITE(6,2)NEWTON2FORMAT(/,30X/*"WARNING - PLASTICITY ALGORITHM DID NOT1'CONVERGE AFTER ',13； ITERATIONS，)140CONTINUEEFFHRD=EG3*HARD1/(EG3+HARD1)CC CALCULATE STRESS AND UPDATE STRAINSCDO 150 K1=1,NDISTRESS(K1)=FLOW(K1)*STELD+SHYDROEPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL/TWOEELAS(Kl)=EELAS(Kl)-THREEFLOW(KirDEQPL/TVVO150 CONTINUEDO 160