疲劳加载（循环应力） 或织构演化分析

疲劳加载（循环应力）

或织构演化分析一、钛合金HCP疲劳加载适配（600℃循环应力）疲劳模拟核心是循环硬化/软化模型

+

滑移系累积损伤演化，以下是针对α-Ti600℃低周疲劳（LCF）的完整适配方案：1.疲劳模型核心扩展（PROPS新增参数）PROPS序号参数名物理意义示例值（α-Ti）18h_cycle循环硬化/软化系数-50.0（软化）19N_fatigue疲劳寿命参考循环数100020D_critical临界损伤值（失效阈值）0.821k_damage损伤演化系数1.0e-52.循环硬化/软化模块（替换单调硬化）fortran!步骤7升级：循环硬化+温度依赖+饱和硬化REAL*8::gamma_rev(NSYSTEM)!反向滑移累积剪切应变（循环加载）REAL*8::cycle_num,gamma_amplitude!循环次数、滑移应变幅值!读取循环相关状态变量（SDV25-36：反向滑移应变；SDV37：循环次数；SDV38：损伤值）IF(KINC.EQ.1.AND.JSTEP.EQ.1)THENDOalpha=1,NSYSTEMSTATEV(alpha+24)=0.0D0!gamma_rev初始化ENDDOSTATEV(37)=0.0D0!循环次数初始化STATEV(38)=0.0D0!损伤值初始化ENDIFcycle_num=STATEV(37)DOalpha=1,NSYSTEMgamma_rev(alpha)=STATEV(alpha+24)ENDDO!识别循环加载方向（应力反向判断）REAL*8::stress_amplitude=MAXVAL(STRESS)-MINVAL(STRESS)IF(stress_amplitude>1e-3.AND.SIGN(1.0D0,tau(1))/=SIGN(1.0D0,STATEV(40)))THENcycle_num=cycle_num+0.5D0!半循环计数STATEV(40)=tau(1)!记录当前滑移应力方向ENDIF!循环硬化/软化修正DOalpha=1,NSYSTEMgamma_amplitude=ABS(gamma(alpha)-gamma_rev(alpha))IF(alpha.LE.6)THEN!基面滑移系!循环软化：h_cycle为负，随循环次数降低硬化能力h_base=PROPS(11)+PROPS(18)*cycle_num/PROPS(19)tau_c(alpha)=tau_cs_base-(tau_cs_base-PROPS(9))*EXP(-h_base*gamma(alpha))ELSE!柱面滑移系h_prism=PROPS(12)+PROPS(18)*cycle_num/PROPS(19)tau_c(alpha)=tau_cs_prism-(tau_cs_prism-PROPS(10))*EXP(-h_prism*gamma(alpha))ENDIF!反向滑移应变更新（应力反向时记录）IF(SIGN(1.0D0,tau(alpha))/=SIGN(1.0D0,gamma(alpha)))THENgamma_rev(alpha)=gamma(alpha)ENDIFENDDO3.滑移系累积损伤演化（基于Manson-Coffin）fortran!步骤8：损伤演化（SDV38为总损伤值）REAL*8::D,D_alpha(NSYSTEM)!总损伤、各滑移系损伤D=STATEV(38)D_alpha=0.0D0!各滑移系损伤：D_alpha=k_damage*gamma_amplitude^2*cycle_numDOalpha=1,NSYSTEMgamma_amplitude=ABS(gamma(alpha)-gamma_rev(alpha))D_alpha(alpha)=PROPS(21)*(gamma_amplitude**2)*cycle_numENDDO!总损伤（最大滑移系损伤主导）D=D+MAXVAL(D_alpha)!损伤阈值判断（超过临界值则刚度退化）IF(D>=PROPS(20))THEN!损伤失效：弹性刚度退化50%C_elastic=C_elastic*0.5D0!输出失效信息WRITE(*,*)'Grain',NOEL,'failedatcycle',cycle_num,'Damage=',DENDIF!更新损伤状态变量STATEV(38)=DSTATEV(37)=cycle_numDOalpha=1,NSYSTEMSTATEV(alpha+24)=gamma_rev(alpha)ENDDO二、织构演化分析（EBSD对比）织构演化核心是晶粒取向旋转模型

+

ODF（取向分布函数）提取，以下是Abaqus模拟与EBSD实验对比的完整流程：1.取向旋转模型（UMAT扩展）fortran!步骤11：晶粒取向旋转（基于滑移系剪切应变）REAL*8::Rot(3,3),Rot_new(3,3),phi1,Phi,phi2!旋转矩阵、欧拉角REAL*8::omega(3)!旋转角速度（由滑移应变梯度驱动）!读取当前欧拉角（SDV39-41：φ1,Φ,φ2）IF(KINC.EQ.1.AND.JSTEP.EQ.1)THEN!初始欧拉角（Neper导入值）STATEV(39)=PROPS(22)!φ1(deg)STATEV(40)=PROPS(23)!Φ(deg)STATEV(41)=PROPS(24)!φ2(deg)ENDIFphi1=STATEV(39)*PI/180.0D0Phi=STATEV(40)*PI/180.0D0phi2=STATEV(41)*PI/180.0D0!计算旋转矩阵CALLEULER2ROT(phi1,Phi,phi2,Rot)!滑移驱动的旋转角速度（简化版：基面滑移主导旋转）omega(1)=0.0D0omega(2)=0.0D0omega(3)=SUM(d_gamma(1:6))*1e-3!绕Z轴旋转，比例系数校准!更新旋转矩阵（小变形近似）Rot_new(1,1)=Rot(1,1)*COS(omega(3)*DTIME)-Rot(1,2)*SIN(omega(3)*DTIME)Rot_new(1,2)=Rot(1,1)*SIN(omega(3)*DTIME)+Rot(1,2)*COS(omega(3)*DTIME)Rot_new(1,3)=Rot(1,3)Rot_new(2,1)=Rot(2,1)*COS(omega(3)*DTIME)-Rot(2,2)*SIN(omega(3)*DTIME)Rot_new(2,2)=Rot(2,1)*SIN(omega(3)*DTIME)+Rot(2,2)*COS(omega(3)*DTIME)Rot_new(2,3)=Rot(2,3)Rot_new(3,:)=Rot(3,:)!旋转矩阵转回欧拉角CALLROT2EULER(Rot_new,phi1,Phi,phi2)!更新欧拉角状态变量STATEV(39)=phi1*180.0D0/PISTATEV(40)=Phi*180.0D0/PISTATEV(41)=phi2*180.0D0/PI2.辅助子程序：旋转矩阵转欧拉角fortran!子程序：旋转矩阵转Z-X-Z欧拉角（适配HCP晶体）SUBROUTINEROT2EULER(Rot,phi1,Phi,phi2)REAL*8,INTENT(IN)::Rot(3,3)REAL*8,INTENT(OUT)::phi1,Phi,phi2REAL*8::eps=1e-6!计算Φ（极角）Phi=ACOS(Rot(3,3))IF(ABS(Phi)<eps.OR.ABS(Phi-PI)<eps)THEN!极角为0或π，欧拉角退化phi1=0.0D0phi2=ATAN2(-Rot(1,2),Rot(1,1))ELSE!通用情况phi1=ATAN2(Rot(1,3),-Rot(2,3))phi2=ATAN2(Rot(3,1),Rot(3,2))ENDIFENDSUBROUTINEROT2EULER三、Abaqus/CAE疲劳加载操作步骤步骤1：定义疲劳分析步进入Step模块，创建Static,General分析步，勾选“Coupledtemperature-displacement”；设置分析步参数：时间总长：1000s（对应1000次循环，1s/循环）；自动步长：最小步长1e-6，最大步长1e-3；循环加载：通过*Amplitude定义正弦波应力幅值（如±200MPa）。步骤2：施加循环载荷与温度边界温度载荷：恒定600℃（873K），通过*Temperature施加；循环应力载荷：进入Load模块，选择*Cload，施加Z向应力；关联*Amplitude,NAME=Fatigue_Amp，定义正弦波：plaintext*Amplitude,Name=Fatigue_Amp0.0,0.01.0,200.02.0,-200.03.0,200.0...!循环至1000s约束条件：固定RVE模型的X/Y平动，释放Z向位移（避免刚体转动）。步骤3：输出织构/疲劳结果变量进入Step模块，编辑分析步的*Output,Field：勾选状态变量（SDV）：SDV1-41（滑移应变、损伤、欧拉角）；勾选应力（S）、应变（E）、温度（T）；设置输出频率：每10个增量步输出一次（捕捉循环特征）。步骤4：提交作业与编译bash运行abaqusjob=ti_hcp_fatigueinput=poly_ti.inpuser=umat_hcp_fatigue.f90-fortlib=/opt/intel/mkl/lib/intel64-cpus=8四、织构演化后处理（与EBSD对比）步骤1：提取欧拉角数据进入Visualization模块，选择Tools→XYData→Create→Element/Nodal；选择“Statevariables”，提取SDV39-41（φ1,Φ,φ2），导出为txt文件。步骤2：ODF分析（与EBSD对比）将Abaqus导出的欧拉角数据导入MTEX（Matlab织构分析工具）；绘制ODF图和极图，对比EBSD实验数据：matlab%MTEX代码示例cs=crystalSymmetry('6/mmm',[2.954.68],'mineral','Ti');%HCP晶系euler=load('abaqus_euler.txt');%导入Abaqus欧拉角数据odf=calcODF(Euler(deg2rad(euler)),cs);%计算ODFplotODF(odf,'sections',6);%绘制ODF图步骤3：疲劳损伤结果提取查看SDV38（损伤值）云图，定位最先失效的晶粒；提取损伤值随循环次数的曲线，拟合Manson-Coffin公式：Nf​⋅γplc​=C（Nf​：疲劳寿命，γpl​：塑性剪切应变幅值，c：材料常数）。五、完整UMAT文件（疲劳+织构演化）fortran!UMATforHCPTifatigue+textureevolution(600℃)INCLUDE'ABA_PARAM.INC'SUBROUTINEUMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,1DDSDDT,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)!变量声明REAL*8,PARAMETER::PI=3.141592653589793D0INTEGER,PARAMETER::NSYSTEM=12REAL*8::C_elastic(6,6),schmid(NSYSTEM,6),tau(NSYSTEM),dot_gamma(NSYSTEM)REAL*8::tau_c(NSYSTEM),gamma(NSYSTEM),d_gamma(NSYSTEM),L_p(6,6),D_ep(6,6)REAL*8::C_elastic_inv(6,6),T_current,tau_cs_base,tau_cs_prismREAL*8::Q_base,Q_prism,T_ref,R,h_base,h_prismREAL*8::gamma_rev(NSYSTEM),cycle_num,gamma_amplitude,D,D_alpha(NSYSTEM)REAL*8::Rot(3,3),Rot_new(3,3),phi1,Phi,phi2,omega(3)INTEGER::i,j,alpha!步骤1：读取参数T_ref=298.0D0R=8.314D0T_current=TEMP!步骤2：初始化状态变量IF(KINC.EQ.1.AND.JSTEP.EQ.1)THEN!滑移应变、临界分切应力DOalpha=1,NSYSTEMSTATEV(alpha)=0.0D0STATEV(alpha+NSYSTEM)=MERGE(PROPS(9),PROPS(10),alpha.LE.6)STATEV(alpha+24)=0.0D0!反向滑移应变ENDDOSTATEV(37)=0.0D0!循环次数STATEV(38)=0.0D0!损伤值STATEV(39:41)=[PROPS(22),PROPS(23),PROPS(24)]!初始欧拉角STATEV(40)=0.0D0!应力方向标记ENDIF!读取状态变量cycle_num=STATEV(37)D=STATEV(38)phi1=STATEV(39)*PI/180.0D0Phi=STATEV(40)*PI/180.0D0phi2=STATEV(41)*PI/180.0D0DOalpha=1,NSYSTEMgamma(alpha)=STATEV(alpha)tau_c(alpha)=STATEV(alpha+NSYSTEM)gamma_rev(alpha)=STATEV(alpha+24)ENDDO!步骤3：构建HCP弹性刚度矩阵（损伤退化）CALLBUILD_C_ELASTIC_HCP(PROPS,C_elastic)IF(D>=PROPS(20))C_elastic=C_elastic*(1.0D0-D)!刚度退化!步骤4：初始化施密特矩阵（取向旋转）CALLINIT_SCHMID_HCP(schmid,NSYSTEM)CALLEULER2ROT(phi1,Phi,phi2,Rot)!旋转施密特矩阵schmid=MATMUL(schmid,Rot)!步骤5：计算分切应力DOalpha=1,NSYSTEMtau(alpha)=DOT_PRODUCT(schmid(alpha,:),STRESS(1:6))ENDDO!步骤6：计算剪切应变率DOalpha=1,NSYSTEMdot_gamma(alpha)=MERGE(0.0D0,PROPS(7)*(ABS(tau(alpha)/tau_c(alpha))**(1.0D0/PROPS(8)))*SIGN(1.0D0,tau(alpha)),ABS(tau_c(alpha))<1e-6)d_gamma(alpha)=dot_gamma(alpha)*DTIMEENDDO!步骤7：循环硬化+温度依赖tau_cs_base=PROPS(14)*EXP(-PROPS(16)*1000.0D0/(R*T_current)+PROPS(16)*1000.0D0/(R*T_ref))tau_cs_prism=PROPS(15)*EXP(-PROPS(17)*1000.0D0/(R*T_current)+PROPS(17)*1000.0D0/(R*T_ref))DOalpha=1,NSYSTEMgamma_amplitude=ABS(gamma(alpha)-gamma_rev(alpha))IF(alpha.LE.6)THENh_base=PROPS(11)+PROPS(18)*cycle_num/PROPS(19)tau_c(alpha)=tau_cs_base-(tau_cs_base-PROPS(9))*EXP(-h_base*gamma(alpha))ELSEh_prism=PROPS(12)+PROPS(18)*cycle_num/PROPS(19)tau_c(alpha)=tau_cs_prism-(tau_cs_prism-PROPS(10))*EXP(-h_prism*gamma(alpha))ENDIF!反向滑移更新IF(SIGN(1.0D0,tau(alpha))/=SIGN(1.0D0,gamma(alpha)))THENgamma_rev(alpha)=gamma(alpha)ENDIFENDDO!步骤8：损伤演化DOalpha=1,NSYSTEMD_alpha(alpha)=PROPS(21)*(gamma_amplitude**2)*cycle_numENDDOD=D+MAXVAL(D_alpha)D=MIN(D,1.0D0)!损伤值不超过1!步骤9：塑性本构矩阵L_p=0.0D0DOalpha=1,NSYSTEMIF(ABS(tau_c(alpha))>1e-6)THENDOi=1,6DOj=1,6L_p(i,j)=L_p(i,j)+schmid(alpha,i)*schmid(alpha,j)*dot_gamma(alpha)/tau_c(alpha)ENDDOENDDOENDIFENDDO!步骤10：弹塑性刚度矩阵CALLINV_MATRIX_6X6(C_elastic,C_elastic_inv)D_ep=C_elastic_inv+L_pCALLINV_MATRIX_6X6(D_ep,DDSDDE)!步骤11：更新应力DOi=1,6STRESS(i)=STRESS(i)+MATMUL(DDSDDE(i,:),DSTRAN(1:6))ENDDO!步骤12：取向旋转更新CALLEULER2ROT(phi1,Phi,phi2,Rot)omega(3)=SUM(d_gamma(1:6))*1e-3Rot_new(1,1)=Rot(1,1)*COS(omega(3)*DTIME)-Rot(1,2)*SIN(omega(3)*DTIME)Rot_new(1,2)=Rot(1,1)*SIN(omega(3)*DTIME)+Rot(1,2)*COS(omega(3)*DTIME)Rot_new(1,3)=Rot(1,3)Rot_new(2,1)=Rot(2,1)*COS(omega(3)