UMAT子程序在復(fù)合材料強度分析中的應(yīng)用

UMAT子程序在復(fù)合材料強度分析中的應(yīng)用UMAT子程序在復(fù)合材料強度分析中的應(yīng)用UMAT子程序在復(fù)合材料強度分析中的應(yīng)用UMAT子程序在復(fù)合材料強度分析中的應(yīng)用編制僅供參考審核批準(zhǔn)生效日期地址：電話：傳真:郵編:UMAT子程序在復(fù)合材料強度分析中的應(yīng)用本例使用UMAT用戶子程序進行復(fù)合材料單層板的應(yīng)力分析和漸進損傷壓縮強度分析，介紹UMAT用戶子程序編寫方法及在Abaqus/CAE中的設(shè)置。本章使用最大應(yīng)變強度理論作為復(fù)合材料單層板的失效準(zhǔn)則，相應(yīng)的Fortran程序簡單易讀，便于理解UAMT子程序的工作原理。知識要點：強度分析UMAT用戶子程序最大應(yīng)變理論剛度折減講師：孔祥宏版本：Abq難度：關(guān)鍵詞：強度分析，UMAT

本章內(nèi)容簡介本章通過兩個實例介紹UMAT用戶子程序在復(fù)合材料單層板的應(yīng)力分析和強度分析中的應(yīng)用。在第一個實例中，對一個簡單的復(fù)合材料單層板進行應(yīng)力分析，UMAT子程序主要計算應(yīng)力，不進行強度分析，本例用于驗證UMAT子程序的計算精度。在第二個實例中，對復(fù)合材料單層板進行漸進損傷強度分析，UMAT子程序用于應(yīng)力計算、強度分析和剛度折減。本章所用復(fù)合材料為T700/BA9916，材料屬性如表&-1所示。表&-1T700/BA9916材料屬性參數(shù)值強度值E1/GPa114XT/MPa2688E2/GPaXC/MPa1458E3/GPaYT/MPaμ12YC/MPa236μ13ZT/MPaμ23ZC/MPa175G12/GPaSXY/MPa136G13/GPaSXZ/MPa136G23/GPaSYZ/MPa實例一：UMAT用戶子程序應(yīng)力分析在使用UMAT用戶子程序進行高級應(yīng)用之前，應(yīng)該先了解UMAT子程序，熟悉UMAT子程序的工作原理，了解UMAT中的參數(shù)、變量的含義。為了便于讀者快速了解和使用UMAT，本例通過復(fù)合材料單層板的應(yīng)力分析來介紹一個簡單的UMAT子程序。讀者可將本例中的單層板替換為層壓板，進行對比分析。&.問題描述復(fù)合材料單層板幾何尺寸為15mm×10mm×，纖維方向為45°，單層板的3D實體模型如圖&-1所示，X軸方向為0°方向，左側(cè)面施加X軸向?qū)ΨQ邊界條件，下側(cè)面施加Y軸向?qū)ΨQ邊界條件，垂直于Z軸且Z=0的平面施加Z軸向?qū)ΨQ邊界條件，右側(cè)面施加100MPa的拉力。圖&-1單層板邊界條件及加載情況本例中單位系統(tǒng)為mm、MPa。&.UMAT用戶子程序本例使用的UMAT用戶子程序的全部代碼如下，字母C及“!”之后為注釋內(nèi)容。SUBROUTINEUMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,1RPL,DDSDDT,DRPLDE,DRPLDT,2STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,3NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,4CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)CINCLUDE''CCHARACTER*80CMNAMEDIMENSIONSTRESS(NTENS),STATEV(NSTATV),1DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),2STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),3PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),4JSTEP(4)DIMENSIONEG(6),XNU(3,3),STRAND(6),C(6,6),STRESS0(6)C****************************CEG.....E1,E2,E3,G12,G13,G23CXNU.....NU12，NU21,NU13，NU31，NU23，NU32CSTRAND.....STRAINTATTHEENDOFTHEINCREMENTCC.....6X6STIFFNESSMATRIXCSTRESS0.....STRESSATTHEBEGINNINGOFTHEINCREMENTC****************************CINITIALIZEXNU&CMATRIXXNU=0C=0CGETTHEMATERIALPROPERTIES---ENGINEERINGCONSTANTSEG(1)=PROPS(1)!E1,YOUNG'SMODULUSINDIRECTION1EG(2)=PROPS(2)!E2,YOUNG'SMODULUSINDIRECTION2EG(3)=EG(2)!E3,YOUNG'SMODULUSINDIRECTION3XNU(1,2)=PROPS(3)!POISON'SRATIOPOI_12XNU(2,1)=XNU(1,2)*EG(2)/EG(1)!POISON'SRATIOPOI_21XNU(1,3)=XNU(1,2)!POISON'SRATIOPOI_13XNU(3,1)=XNU(1,3)*EG(3)/EG(1)!POISON'SRATIOPOI_31XNU(2,3)=PROPS(4)!POISON'SRATIOPOI_23XNU(3,2)=XNU(2,3)*EG(3)/EG(2)!POISON'SRATIOPOI_32EG(4)=PROPS(5)!G12,SHEARMODULUSIN12PLANEEG(5)=EG(4)!G13,SHEARMODULUSIN13PLANEEG(6)=PROPS(6)!G23,SHEARMODULUSIN23PLANEC****************************CFILLTHE6X6STIFFNESSMATRIXC(6,6)RNU=1/(1-XNU(1,2)*XNU(2,1)-XNU(1,3)*XNU(3,1)-1XNU(3,2)*XNU(2,3)-2*XNU(1,3)*XNU(2,1)*XNU(3,2))CSTIFFNESSMATRIXC(6,6)C(1,1)=EG(1)*(1-XNU(2,3)*XNU(3,2))*RNUC(2,2)=EG(2)*(1-XNU(1,3)*XNU(3,1))*RNUC(3,3)=EG(3)*(1-XNU(1,2)*XNU(2,1))*RNUC(4,4)=EG(4)C(5,5)=EG(5)C(6,6)=EG(6)C(1,2)=EG(1)*(XNU(2,1)+XNU(3,1)*XNU(2,3))*RNUC(2,1)=C(1,2)C(1,3)=EG(1)*(XNU(3,1)+XNU(2,1)*XNU(3,2))*RNUC(3,1)=C(1,3)C(2,3)=EG(2)*(XNU(3,2)+XNU(1,2)*XNU(3,1))*RNUC(3,2)=C(2,3)C****************************CCALCULATESTRAINDOI=1,6STRAND(I)=STRAN(I)+DSTRAN(I)ENDDOCCALCULATESTRESSDOI=1,6STRESS0(I)=STRESS(I)STRESS(I)=0DOJ=1,6STRESS(I)=STRESS(I)+C(I,J)*STRAND(J)ENDDOENDDOCCALCULATESSEDOI=1,6SSE=SSE+*(STRESS0(I)+STRESS(I))*DSTRAN(I)ENDDOC****************************CUPDATEDDSDDEDOI=1,6DOJ=1,6DDSDDE(I,J)=C(I,J)ENDDOENDDORETURNEND第1到14行及第81、82行為UMAT子程序固定格式，其中，第1到5行括號內(nèi)的變量為UMAT子程序中可以使用的變量，第10到14行定義各變量數(shù)組的維數(shù)和長度。部分主要變量的含義如表&-2所示。表&-2UMAT部分變量名及其含義STRESS增量步開始時的應(yīng)力(S11,S22,...)，用增量步結(jié)束時的應(yīng)力計算結(jié)果對其更新STATEV(NSTATV)狀態(tài)變量（狀態(tài)變量個數(shù)），如果在材料中定義了狀態(tài)變量，則在UMAT中需要對其更新STRAN增量步開始時的應(yīng)變(E11,E22,...)DSTRAN當(dāng)前增量步的應(yīng)變增量(ΔE11,ΔE22,...)NDI,NSHR,NTENS應(yīng)力、應(yīng)變的個數(shù)，NDI為正應(yīng)力或正應(yīng)變的個數(shù)，NSHR為剪應(yīng)力或剪應(yīng)變的個數(shù)，NTENS=NDI+NSHRPROPS,NPROPS材料參數(shù)、材料參數(shù)的個數(shù)DDSDDE雅克比矩陣，SSE,SPD,SCD特定的彈性應(yīng)變能、塑性耗散、蠕變耗散，只對能量輸出有影響，對其他計算結(jié)果無影響，在UMAT中需要對其更新CELENT單元特征長度第15到83行為用戶自己編寫的固定格式的Fortran程序，用于計算剛度矩陣、應(yīng)力、應(yīng)變能、雅克比矩陣。由于本例中沒有使用狀態(tài)變量，因此不需要更新STATEV，只需要更新STRESS、DDSDDE和SSE即可。第16行定義了5個數(shù)組，其中EG、STRAND、STRESS0為一維數(shù)組，XNU、C為二維數(shù)組。第18到22行為注釋部分，EG存放材料的3個彈性模量和3個剪切模量；STRAND存放當(dāng)前增量步結(jié)束時的應(yīng)變(E11,E22,...)；STRESS0存放增量步開始時的應(yīng)力(S11,S22,...)；XNU為3×3的二維矩陣，存放泊松比ν12、ν21、ν13、ν31、ν23、ν32；C為6×6的剛度矩陣。第25、26行初始化二維數(shù)組XNU、C，使其每個元素都為0。第28到39行，讀取材料常數(shù)，計算泊松比。泊松比的計算公式如下。(&-1)第42到56行，計算剛度矩陣。剛度矩陣的計算公式如下。(&-2)對于本例所用材料，由于E2=E3，G12=G13，ν12=ν13，所以ν21=ν31、ν23=ν32，可以將式(&-2)化簡后在UMAT中計算剛度矩陣。本例為保證的可讀性，剛度矩陣的計算沒有化簡，直接按照式(&-2)編寫。第59到61行，計算當(dāng)前增量步的應(yīng)變，計算公式如式(&-3)，式中上標(biāo)表示增量步序號。(&-3)第63到69行，將當(dāng)前分析步開始時的應(yīng)力STRESS的值賦給STRESS0，然后計算當(dāng)前增量步的應(yīng)力并賦給STRESS，當(dāng)前增量步的應(yīng)力計算如式(&-4)，式中上標(biāo)表示增量步序號，應(yīng)力σi和應(yīng)變εi為列向量。(&-4)第71到73行，計算應(yīng)變能，如式(&-5)所示，式中上標(biāo)表示增量步序號；下標(biāo)表示應(yīng)力、應(yīng)變增量的分量序號，其中下標(biāo)為1、2、3表示正應(yīng)力、正應(yīng)變增量，4、5、6表示剪應(yīng)力、剪應(yīng)變增量。(&-5)第76到80行為更新雅克比矩陣。由于沒有對剛度矩陣沒變化，應(yīng)力與應(yīng)變的關(guān)系如式(&-4)所示，所以應(yīng)力增量與應(yīng)變增量的關(guān)系如式(&-6)所示，所以雅克比矩陣就等于剛度矩陣。(&-6)第76到80行可簡寫為一行，如下所示：DDSDDE(1:6,1:6)=C(1:6,1:6)，或DDSDDE=CUMAT中的剪應(yīng)變?yōu)楣こ碳魬?yīng)變。&.復(fù)合材料單層板應(yīng)力分析創(chuàng)建部件及劃分網(wǎng)格創(chuàng)建部件在Part模塊，單擊工具區(qū)的(CreatePart)，在CreatePart對話框中，Name后面輸入Lam-C，ModelingSpace選擇3D，Type選擇Deformable，在BaseFeature區(qū)域選擇Solid、Extrusion，Approximatesize使用默認的200，單擊Continue...進入繪圖模式。單擊工具區(qū)的(CreateLines:Rectangle(4Lines))，在提示區(qū)輸入第1個點的坐標(biāo)(0,0)后按回車鍵，再輸入第2個點的坐標(biāo)(15,10)后按回車，再按Esc鍵或單擊鼠標(biāo)中鍵。單擊提示區(qū)的Done或鼠標(biāo)中鍵，在EditBaseExtrusion對話框Depth后面輸入，單擊OK完成。劃分網(wǎng)格在環(huán)境欄Module后面選擇Mesh，進入Mesh模塊。環(huán)境欄中Object選擇Part:Lam-C。單擊工具區(qū)的(SeedPart)，在GlobalSeeds對話框中Approximateglobalsize后面輸入，單擊OK。單擊工具區(qū)的(MeshPart)，單擊提示區(qū)的Yes或鼠標(biāo)中鍵，完成網(wǎng)格劃分，如圖&-2所示，板的厚度方向只劃分為1層單元。圖&-2劃分網(wǎng)格單擊工具區(qū)的(AssignElementType)，在ElementType對話框中，選擇依次選擇Standard，Linear，3DStress，在Hex標(biāo)簽頁中勾選Reducedintegration，在ElementControls區(qū)域Hourglasscontrol選擇Enhanced，即選擇C3D8R單元，單擊OK完成。單擊工具區(qū)的(AssignStackDirection)，在視圖區(qū)選擇部件平行于X-Y平面的面，單擊鼠標(biāo)中鍵或提示區(qū)的Yes完成。創(chuàng)建材料并給部件賦材料屬性創(chuàng)建材料在環(huán)境欄Module后面選擇Property，進入Property模塊。單擊工具區(qū)的(CreateMaterial)，在EditMaterial對話框中，Name后面輸入UMat-T700；單擊General→UserMaterial，在UserMaterial區(qū)域中Data區(qū)域的MechanicalConstants一欄依次輸入114000，8610，，，4160，3000，單擊OK完成。單擊工具區(qū)的(CreateMaterial)，在EditMaterial對話框中，Name后面輸入Mat-T700；單擊Mechanical→Elasticity→Elastic，在Elastic區(qū)域中Type選擇EngineeringConstants，在Data區(qū)域輸入從左到右依次輸入114000，8610，8610，，，，4160，4160，3000，單擊OK完成。材料UMat-T700用于UMAT用戶子程序，Mat-T700用于做對比分析。輸入數(shù)據(jù)時，每輸入完一行后按回車（Enter）鍵，光標(biāo)會自動移到下一行。也可以通過右鍵快捷菜單添加或刪除一行。本例UMAT子程序較簡單，不需要使用狀態(tài)變量，因此在材料UMat-T700中沒有定義Depvar。給部件賦材料屬性單擊工具區(qū)的(CreateCompositeLayup)，在打開的對話框中Name使用默認名稱，Initialplycount后面輸入1，ElementType選擇Solid，單擊Continue...；在EditCompositeLayup對話框中，LayupOrientation區(qū)域的Definition選擇Coordinatesystem，單擊Definition下一行的(CreateDatumCSYS)。在CreateDatumCSYS對話框中使用默認名稱Datumcsys-1，類型選擇Rectangular，單擊Continue...，在提示區(qū)輸入原點坐標(biāo)(0,0,0)后按回車鍵，再輸入(1,0,0)后按回車鍵，最后輸入(0,1,0)后按回車鍵，單擊CreateDatumCSYS對話框的Cancel。在EditCompositeLayup對話框中單擊(SelectCSYS...)，在視圖區(qū)選擇剛創(chuàng)建的Datumcsys-1，StackingDirection選擇Elementdirection3，Rotationaxis選擇Axis3。在Plies標(biāo)簽頁中，雙擊Region，在視圖區(qū)選擇部件后單擊鼠標(biāo)中鍵；右單擊Material，在快捷菜單中單擊EditMaterial...，在SelectMaterial對話框中選擇UMat-T700，單擊OK；右單擊ElementRelativeThickness，在快捷菜單中單擊EditThickness...，在Thickness對話框中SpecifyValue后面輸入1后單擊OK；在RotationAngle一欄輸入0；IntegrationPoints使用默認的1，單擊OK完成。在定義復(fù)合材料鋪層時，視圖區(qū)部件上會顯示鋪層方向，在EditCompositeLayup對話框中的Display標(biāo)簽頁中可以設(shè)置所需顯示的方向，在視圖區(qū)部件上白色箭頭及字母S表示StackingDirection。裝配在環(huán)境欄Module后面選擇Assembly，進入Assembly模塊。單擊工具區(qū)的(CreateInstance)，在CreateInstance對話框中選擇Parts:Lam-C，單擊OK完成。創(chuàng)建分析步、設(shè)置輸出變量創(chuàng)建分析步在環(huán)境欄Module后面選擇Step，進入Step模塊。單擊工具區(qū)的(CreateStep)，在CreateStep對話框中，在Initial分析步之后插入Static,General分析步，單擊Continue...；在EditStep對話框中使用默認設(shè)置，單擊OK完成。設(shè)置輸出變量單擊工具區(qū)的(FieldOutputManager)，在FieldOutputRequestsManager對話框中，選中F-Output-1，單擊Edit...；在EditFieldOutputRequest對話框中，設(shè)置如圖&-3所示，輸出整個模型最后一個增量步的S、E、U，單擊OK完成。圖&-3場輸出變量設(shè)置單擊工具區(qū)的(HistoryOutputManager)，在HistoryOutputRequestsManager對話框中，選中H-Output-1，單擊Edit...；在EditHistoryOutputRequest對話框中，設(shè)置輸出整個模型的內(nèi)能和應(yīng)變能，即ALLIE和ALLSE，單擊OK完成。創(chuàng)建邊界條件及施加載荷創(chuàng)建邊界條件在環(huán)境欄Module后面選擇Load，進入Load模塊。單擊工具區(qū)的（CreateBoundaryCondition），在CreateBoundaryCondition對話框中，Name后面輸入BC-X，Step選擇Initial，Category選擇Mechanical，TypesforSelectedStep選擇Symmetry...，單擊Continue...；在視圖區(qū)選擇裝配實例Lam-C-1左側(cè)端面，即垂直于X軸且X=0的側(cè)面，單擊鼠標(biāo)中鍵或提示區(qū)的Done，在EditBoundaryCondition對話框中選擇XSYMM，單擊OK完成。類似操作，選擇Lam-C-1的下側(cè)端面，即垂直于Y軸且Y=0的側(cè)面，創(chuàng)建邊界條件BC-Y，邊界類型為YSYMM；選擇Lam-C-1垂直于Z軸且Z=0的側(cè)面，創(chuàng)建邊界條件BC-Z，邊界類型為ZSYMM。邊界條件及加載情況見圖&-1。施加載荷單擊工具區(qū)的(CreateLoad)，在CreateLoad對話框中，Name使用默認的Load-1，Step選擇Step-1，Category選擇Mechanical，TypesforSelectedStep選擇Pressure，單擊Continue...；在視圖區(qū)選擇Lam-C-1的右側(cè)端面，即垂直于X軸且X=15的側(cè)面，單擊鼠標(biāo)中鍵，在EditLoad對話框中Magnitude后面輸入-100，單擊OK完成。創(chuàng)建分析作業(yè)并提交分析創(chuàng)建分析作業(yè)在環(huán)境欄Module后面選擇Job，進入Job模塊。單擊工具區(qū)的（JobManager），在JobManager對話框中單擊Create...；在CreateJob對話框中，Name后面輸入Job-Lam-Stress-Umat，Source選擇Model-1，單擊Continue...；在EditJob對話框的General標(biāo)簽頁中，單擊Usersubroutinefile后面的(Select...)，在相應(yīng)路徑下找到并選擇文件；單擊EditJob對話框中的OK完成。在EditJob對話框的Parallelization標(biāo)簽頁中可以設(shè)置多核并行計算。提交分析在JobManager對話框中，選中Job-Lam-Stress-Umat分析作業(yè)，單擊Submit提交計算。當(dāng)Job-Lam-Stress-Umat的狀態(tài)（Status）由Running變?yōu)镃ompleted時，計算完成，單擊Results進入可視化后處理模塊。保存模型單擊工具欄的File工具條中的（SaveModelDatabase），在SaveModelDatabaseAs對話框的FileName后面輸入Laminate-Umat，單擊OK完成。修改材料在環(huán)境欄Module后面選擇Property，進入Property模塊。單擊工具區(qū)的(CompositeLayupManager)，在打開的對話框中選擇CompositeLayup-1，單擊Edit...；在EditCompositeLayup對話框的Plies標(biāo)簽頁中，右單擊Material，單擊快捷菜單中的EditMaterial...；在SelectMaterial對話框中選擇Mat-T700，單擊OK；在EditCompositeLayup對話框中單擊OK完成。再次創(chuàng)建分析作業(yè)并提交分析創(chuàng)建分析作業(yè)在環(huán)境欄Module后面選擇Job，進入Job模塊。單擊工具區(qū)的（JobManager），在JobManager對話框中單擊Create...；在CreateJob對話框中，Name后面輸入Job-Lam-Stress，Source選擇Model-1，單擊Continue...；在EditJob對話框使用默認設(shè)置，單擊OK完成。在EditJob對話框的Parallelization標(biāo)簽頁中可以設(shè)置多核并行計算。提交分析在JobManager對話框中，選中Job-Lam-Stress分析作業(yè)，單擊Submit提交計算。當(dāng)Job-Lam-Stress的狀態(tài)（Status）由Running變?yōu)镃ompleted時，計算完成，單擊Results進入可視化后處理模塊。保存模型單擊工具欄的File工具條中的(SaveModelDatabase)保存模型?？梢暬筇幚盹@示云圖在視圖區(qū)顯示。長按工具區(qū)的(PlotContoursonDeformedShape)，顯示隱藏工具后單擊(PlotContoursonUndeformedShape)，在FieldOutput工具條中設(shè)置輸出S11。單擊菜單欄Result→SectionPoints...，打開SectionPoints對話框，Selectionmethod選擇Plies，在Plies區(qū)域選擇PLY-1，單擊Apply，在視圖區(qū)顯示該鋪層的S11應(yīng)力云圖。相同操作，可以顯示各鋪層的各應(yīng)力分離的云圖。和的各應(yīng)力分量的云圖如圖&-4所示。(a)S11，Job-Lam-Stress-Umat(b)S11，Job-Lam-Stress(c)S22，Job-Lam-Stress-Umat(d)S22，Job-Lam-Stress(e)S12，Job-Lam-Stress-Umat(f)S12，Job-Lam-Stress圖&-4各應(yīng)力分量對比讀者在閱讀Abaqus幫助文件AbaqusExampleProblemsGuide中Failureofbluntnotchedfibermetallaminates一例的文件時注意SSE的計算公式。輸出應(yīng)變能在視圖區(qū)顯示。單擊工具區(qū)的(XYDataManager)，在打開的對話框中單擊Create...；在CreateXYData對話框中選擇ODBhistoryoutput，單擊Continue...；在HistoryOutput對話框中選擇ALLSE，單擊SaveAs...；在SaveXYDataAs對話框中Name使用默認的XYData-1，單擊OK；在XYDataManager對話框中選擇XYData-1，單擊Edit...；在EditXYData對話框中可以看到應(yīng)變能數(shù)據(jù)，分析結(jié)束時整個模型的應(yīng)變能為。同樣操作，觀察中整個模型的應(yīng)變能數(shù)據(jù)，同樣為。通過對比，驗證了中第72行SSE計算的準(zhǔn)確性。單擊菜單欄Result→HistoryOutput...，可以直接打開HistoryOutput對話框。&.Inp文件解釋本例中節(jié)選如下。*Heading**Jobname:Job-Lam-Stress-UmatModelname:Model-1**Generatedby:Abaqus/CAE*Preprint,echo=NO,model=NO,history=NO,contact=NO****PARTS**部件Lam-C-1的節(jié)點、單元數(shù)據(jù)*Part,name=Lam-C*Node1,15.,10.,……節(jié)點編號及坐標(biāo)*Element,type=C3D8R1,64,65,23,22,43,44,2,1……單元編號及節(jié)點編號*Orientation,name=Ori-11.,0.,0.,0.,1.,0.3,0.**Section:CompositeLayup-1-1**定義復(fù)合材料鋪層*SolidSection,elset=CompositeLayup-1-1,composite,orientation=Ori-1,controls=EC-1,stackdirection=3,layup=CompositeLayup-11,1,UMat-T700,0.,Ply-1*EndPart****ASSEMBLY*Assembly,name=Assembly**使用部件Lam-C創(chuàng)建裝配實例Lam-C-1*Instance,name=Lam-C-1,part=Lam-C*EndInstance*EndAssembly**ELEMENTCONTROLS**單元控制，使用沙漏控制*SectionControls,name=EC-1,hourglass=ENHANCED1.,1.,1.****MATERIALS**使用工程常數(shù)創(chuàng)建材料Mat-T700*Material,name=Mat-T700*Elastic,type=ENGINEERINGCONSTANTS114000.,8610.,8610.,,,,4160.,4160.3000.,**使用用戶材料創(chuàng)建材料UMat-T700*Material,name=UMat-T700*UserMaterial,constants=6114000.,8610.,,,4160.,3000.****BOUNDARYCONDITIONS**創(chuàng)建對稱邊界條件BC-X，BC-Y和BC-Z略**Name:BC-XType:Symmetry/Antisymmetry/Encastre*Boundary_PickedSet4,XSYMM**----------------------------------------------------------------**STEP:Step-1**創(chuàng)建分析步Step-1*Step,name=Step-1,nlgeom=NO*Static1.,1.,1e-05,1.****LOADS**在Step-1創(chuàng)建Pressure類型的載荷**Name:Load-1Type:Pressure*Dsload_PickedSurf6,P,-100.****OUTPUTREQUESTS**設(shè)置輸出變量*Restart,write,frequency=0****FIELDOUTPUT:F-Output-1**場輸出變量U、E、S*Output,field,frequency=99999*NodeOutputU,*ElementOutput,directions=YESE,S****HISTORYOUTPUT:H-Output-1**設(shè)置歷史輸出變量*Output,history*EnergyOutputALLIE,ALLSE*EndStep&.應(yīng)用UMAT子程序應(yīng)力分析小結(jié)本例所介紹的應(yīng)力分析是UMAT用戶子程序最簡單的應(yīng)用，讀者在了解UMAT子程序的基礎(chǔ)上，可以結(jié)合相應(yīng)的強度理論、剛度折減方法對UMAT子程序中的剛度矩陣、雅克比矩陣進行計算和更新，從而達到剛度折減和漸進損傷強度分析的目的。實例二：UMAT用戶子程序漸進損傷強度分析本例在實例一使用的UMAT子程序的基礎(chǔ)上，通過增加失效判定，對UMAT子程序中的剛度矩陣進行折減，達到漸進損傷強度分析的目的。為了使本例的UMAT子程序簡單易讀，本例對復(fù)合材料單層板進行纖維方向的拉伸強度分析，僅考慮纖維方向拉伸破壞。本例使用本章實例一的模型，稍做修改，用于本例的強度分析。&.問題描述復(fù)合材料單層板的3D實體模型的幾何尺寸及邊界條件同本章的實例一的模型。將右側(cè)面100MPa的拉力替換為X軸正方向的位移載荷。本例中單位系統(tǒng)為mm、MPa。&.UMAT用戶子程序本例使用的UMAT用戶子程序的全部代碼如下，字母C及“!”之后為注釋內(nèi)容。SUBROUTINEUMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,1RPL,DDSDDT,DRPLDE,DRPLDT,2STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,3NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,4CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)CINCLUDE''CCHARACTER*80CMNAMEDIMENSIONSTRESS(NTENS),STATEV(NSTATV),1DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),2STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),3PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),4JSTEP(4)DIMENSIONSTRESS0(6),STRAND(6),C(6,6),CD(6,6),DCDE(6,6)PARAMETER(ALPHA=,LAMBDA=,DMAX=,DRND=3)C****************************CSTRESS0.....STRESSATTHEBEGINNINGOFTHEINCREMENTCSTRAND.....STRAINATTHEENDOFTHEINCREMENTCC.....6X6STIFFNESSMATRIXCCD.....6X6DAMAGEDSTIFFNESSMATRIXCDCDE.....DCD/DECSTATEV(1).....DAMAGEVARIABLEDC****************************CGETTHEMATERIALPROPERTIESE1=PROPS(1)!E1,YOUNG'SMODULUSINDIRECTION1E2=PROPS(2)!E2=E3,YOUNG'SMODULUSINDIRECTION2&3XNU12=PROPS(3)!POISON'SRATIOPOI_12,XNU13=XNU12XNU21=XNU12*E2/E1!POISON'SRATIOPOI_21,XNU31=XNU21XNU23=PROPS(4)!POISON'SRATIOPOI_23,XNU32=XNU23G12=PROPS(5)!G12=G13,SHEARMODULUSIN12&13PLANEG23=PROPS(6)!G23,SHEARMODULUSIN23PLANESTH=PROPS(7)!FAILURESTRESSIN1DIRECTIONINTENSIONC****************************CSTIFFNESSMATRIXC(6,6)RNU=1/(1-2*XNU12*XNU21-XNU23**2-2*XNU12*XNU21*XNU23)C=0C(1,1)=E1*(1-XNU23**2)*RNUC(2,2)=E2*(1-XNU12*XNU21)*RNUC(3,3)=C(2,2)C(4,4)=G12C(5,5)=G12C(6,6)=G23C(1,2)=E1*(XNU21+XNU21*XNU23)*RNUC(2,1)=C(1,2)C(1,3)=C(1,2)C(3,1)=C(1,2)C(2,3)=E2*(XNU23+XNU12*XNU21)*RNUC(3,2)=C(2,3)CCALCULATETHESTRAINATTHEENDOFTHEINCREMENTDOI=1,6STRAND(I)=STRAN(I)+DSTRAN(I)ENDDOC****************************CCALCULATETHEFAILURECOEFFICIENTSTRANF=STH/E1IF(STRAND(1)>0)THENF=STRAND(1)/STRANFELSEF=0ENDIFCCALCULATED,DAMAGEVARIABLED=STATEV(1)DDDE=0IF(F>1)THENDV=1-EXP(ALPHA*(1-F))IF(DV>D)THEND=D*LAMBDA/(LAMBDA+1)+DV/(LAMBDA+1)DDDE=ALPHA*(1-DV)/STRANF/(1+LAMBDA)D=ANINT(D*10**DRND)/10**DRNDDDDE=ANINT(DDDE*10**DRND)/10**DRNDENDIFIF(D>DMAX)THEND=DMAXENDIFENDIFSTATEV(1)=D!UPDATEDCDAMAGEDSTIFFNESSMATRIXCD(6,6)CD=CCD(1,1)=(1-D)*C(1,1)CD(1,2)=(1-D)*C(1,2)CD(1,3)=CD(1,2)CD(2,1)=CD(1,2)CD(3,1)=CD(1,2)CD(4,4)=(1-D)*C(4,4)CD(5,5)=(1-D)*C(5,5)C****************************CCALCULATESTRESSDOI=1,6STRESS0(I)=STRESS(I)STRESS(I)=DOJ=1,6STRESS(I)=STRESS(I)+CD(I,J)*STRAND(J)ENDDOENDDOCCALCULATESSE,ELASTICSTRAINENERGYDOI=1,NTENSSSE=SSE+*(STRESS0(I)+STRESS(I))*DSTRAN(I)ENDDOC****************************CUPDATETHEJACOBIANDCDE=0DCDE(1,1)=-DDDE*CD(1,1)DCDE(1,2)=-DDDE*CD(1,2)DCDE(1,3)=DCDE(1,2)DCDE(2,1)=DCDE(1,2)DCDE(3,1)=DCDE(1,2)DCDE(4,4)=-DDDE*CD(4,4)DCDE(5,5)=-DDDE*CD(5,5)DDSDDE=CDDOI=1,6DOJ=1,6ATEMP=DCDE(I,J)*STRAND(J)ENDDODDSDDE(I,1)=DDSDDE(I,1)+ATEMPENDDORETURNEND第1到15行及第118、119行為UMAT子程序固定格式，第16到117行為本程序的主體部分。第16行定義了5個數(shù)組，其含義見第18到25行之間的注釋部分。由于本例僅考慮纖維方向拉伸破壞，且使用最大應(yīng)變強度理論，所以DCDE為折減后的剛度矩陣Cd對纖維方向正應(yīng)變ε1的偏導(dǎo)矩陣，即矩陣Cd中所以元素分別對ε1求偏導(dǎo)。本例只使用一個狀態(tài)變量STATEV(1)記錄單元纖維方向拉伸破壞的破壞變量d。第17行定義了4個常數(shù)參數(shù)，其中ALPHA(即α)、LAMBDA(即λ)用于計算破壞變量d，DMAX(即dmax)為破壞變量d的最大值，DRND(即dr)為破壞變量d保留小數(shù)的位數(shù)。第27到34行，讀取材料參數(shù)。其中，第34行STH為材料纖維方向的拉伸強度。第37到50行，計算材料無損傷時的剛度矩陣C。第52到54行，計算當(dāng)前增量步結(jié)束時的應(yīng)變ε=(ε1,ε2,ε3,γ12,γ13,γ23)，其中剪應(yīng)變?yōu)楣こ碳魬?yīng)變。第57到62行，使用最大應(yīng)變強度理論計算纖維方向拉伸失效系數(shù)f。纖維方向的拉伸失效應(yīng)變ε0及失效系數(shù)f的計算公式如下：(&-7)式(&-7)中為材料纖維方向拉伸強度，E1為材料纖維方向彈性模量，ε1為當(dāng)前分析步結(jié)束時纖維方向的拉伸正應(yīng)變。第64到78行，計算破壞變量D（即d）和破壞系數(shù)對應(yīng)變的偏導(dǎo)數(shù)（導(dǎo)數(shù)）DDDE（即）。破壞變量計算公式如下：(&-8)式(&-8)中dv為當(dāng)前增量步結(jié)束時的破壞變量，di為前一增量步結(jié)束時（或當(dāng)前增量步開始時）的破壞變量，di+1為當(dāng)前增量步結(jié)束時的破壞變量。由于僅考慮纖維方向拉伸破壞，式(&-7)中失效系數(shù)f僅為ε1的函數(shù)，破壞變量d僅為f的函數(shù)，因此DDDE可以寫為：(&-9)式(&-9)中的d對應(yīng)式(&-8)中的di+1。如果d對為ε的函數(shù)，則為一個向量，即。第71、72行，按照參數(shù)DRND的值對D和DDDE保留指定的小數(shù)位數(shù)。由于使用狀態(tài)變量STATEV(1)保存破壞變量，因此要用式(&-8)中的di+1對STATEV(1)進行更新。第80到87行，計算破壞后的剛度矩陣Cd，由于近考慮纖維方向拉伸破壞，因此剛度折減如式(&-10)所示。(&-10)第90到100行，計算應(yīng)力STRESS和應(yīng)變能SSE。第103到110行，計算DCDE，即，對于本例，實際上只計算。第111到117行，計算雅克比矩陣。雅克比矩陣的基本計算公式如下：(&-11)式(&-11)中J為NTENS階方陣，其中Jij可表示為：(&-12)在式(&-11)中，(&-13)式(&-13)中，(i=1,2,...)為一個列向量，該向量各元素分別加到雅克比矩陣J的第i列的相應(yīng)元素上。&.復(fù)合材料單層板強度分析修改材料及復(fù)合材料鋪層修改材料打開本章實例一創(chuàng)建的。在環(huán)境欄Module后面選擇Property，進入Property模塊。單擊工具區(qū)的(MaterialManager)，在MaterialManager對話框中選擇UMat-T700，單擊Edit...；在EditMaterial對話框中單擊General→Depvar，在Numberofsolution-dependentstatevariable后面輸入1；單擊MaterialBehaviors區(qū)域的UserMaterial，勾選Useunsymmetricmaterialstiffnessmatirx，在Data區(qū)域MechanicalConstants一欄的最后添加2688，即材料纖維方向的拉伸強度2688MPa；單擊OK完成。修改復(fù)合材料鋪層單擊工具區(qū)的(CompositeLayupManager)，在打開的對話框中選擇CompositeLayup-1，單擊Edit...；在Plies標(biāo)簽頁中，將Ply-1的Material改為UMat-T700，RotationAngle改為0，單擊OK完成。創(chuàng)建參考點和約束創(chuàng)建參考點在環(huán)境欄Module后面選擇Interaction，進入Interaction模塊。單擊工具區(qū)的(CreateReferencePoint)，在提示區(qū)輸入(15,5,，按回車鍵完成，單擊單擊工具區(qū)的退出創(chuàng)建參考點。創(chuàng)建集合單擊菜單欄Tools→Set→Manager...，在SetManager對話框中單擊Create...；在CreateSet對話框中Name后面輸入Set-RP，Type選擇Geometry，單擊Continue...；在視圖區(qū)選擇RP-1，單擊鼠標(biāo)中鍵完成。再單擊SetManager對話框中單擊Create...；在CreateSet對話框中Name后面輸入Set-Nds，Type選擇Node，單擊Continue...；在視圖區(qū)選擇部件右側(cè)面（垂直于X軸且X=15的側(cè)面）上的所有節(jié)點，單擊鼠標(biāo)中鍵完成。選擇一個面上的所有節(jié)點時，可以在提示區(qū)將選擇方法設(shè)為byangle。創(chuàng)建約束單擊工具區(qū)的(CreateConstraint)，在打開的對話框中Name使用默認名稱，Type選擇Equation，單擊Continue...；在EditConstraint對話框設(shè)置如圖&-5所示，單擊OK完成。圖&-5設(shè)置約束修改分析步、輸出變量修改分析步在環(huán)境欄Module后面選擇Step，進入Step模塊。單擊工具區(qū)的(StepManager)，在StepManager對話框中選擇Step-1，單擊Edit...；在EditStep對話框的Basic標(biāo)簽頁中，Nlgeom選擇On；在Incrementation標(biāo)簽頁設(shè)置如圖&-6所示，單擊OK完成。圖&-6增量步設(shè)置修改輸出變量單擊工具區(qū)的(FieldOutputManager)，在FieldOutputRequestsManager對話框中，選中F-Output-1，單擊Edit...；在EditFieldOutputRequest對話框中，設(shè)置如圖&-7所示，單擊OK完成。圖&-7場輸出變量設(shè)置單擊工具區(qū)的(HistoryOutputManager)，在HistoryOutputRequestsManager對話框中，選中H-Output-1，單擊Edit...；在EditHistoryOutputRequest對話框中設(shè)置如圖&-8，單擊OK完成。圖&-8歷史輸出變量設(shè)置修改載荷刪除舊載荷在環(huán)境欄Module后面選擇Load，進入Load模塊。單擊工具區(qū)的（LoadManager），在LoadManager對話框中選擇Load-1，單擊Delete...，單擊提示框中的Yes，關(guān)閉LoadManager對話框。可以再模型樹Loads中選擇Load-1，按Delete鍵刪除。施加位移載荷單擊工具區(qū)的(CreateBoundaryCondition)，在打開的對話框中，Name后面輸入U1，Step選擇Step-1，Category選