晶体塑性有限元：从基础理论到工程应用的全套学习指南与资料汇编

晶体塑性有限元学习资料晶体塑性有限元这东西说难也难满篇的张量推导能把人看晕说简单也简单核心就是抓住“晶体变形靠滑移”这一点——说白了晶体不是像橡皮泥那样随便捏的得沿着特定的滑移面和滑移方向动这就是它和常规弹塑性最大的区别。咱们拿最基础的单晶体拉伸模拟开刀直接看怎么把理论撸成代码以ABAQUS的VUMAT为例毕竟工程里用得最多。首先得定义晶体的取向也就是欧拉角这是晶体相对于试样坐标系的旋转角度。代码里直接给三个数就行C --- 定义单晶体取向欧拉角ZXZ顺序 PARAMETER (EULER10.0D0, EULER20.0D0, EULER30.0D0)别小看这三个数我第一次学的时候把欧拉角顺序搞反了ABAQUS默认ZXZ我写成了XYZ结果模拟出来的[001]取向铜拉伸屈服强度跟[111]的一模一样查了半天才发现是取向转错了方向——晶体塑性里取向就是“性格”性格错了行为当然不对。接下来是核心计算每个滑移系的Schmid因子判断哪个滑移系先激活。铜是面心立方有12个滑移系咱们先把这些滑移系的面法向量m和滑移方向s定义好C --- 面心立方滑移系定义12个 INTEGER, PARAMETER :: N_SLIP 12 REAL*8 M(3,N_SLIP), S(3,N_SLIP) C --- 第一个滑移系(111)[1-10] M(:,1) [1.0D0/SQRT(3.0D0), 1.0D0/SQRT(3.0D0), 1.0D0/SQRT(3.0D0)] S(:,1) [1.0D0/SQRT(2.0D0), -1.0D0/SQRT(2.0D0), 0.0D0] C --- 剩下11个滑移系以此类推比如(111)[-110]、(111)[10-1]...这段代码就是把每个滑移系的两个关键向量存起来。然后计算Schmid因子和分切应力——这直接决定了哪个滑移面先“启动”C --- 计算每个滑移系的分切应力 REAL*8 TAU(N_SLIP), TAU020.0D0 ! TAU0是初始临界分切应力MPa DO IS1,N_SLIP TAU(IS) 0.0D0 DO I1,3 DO J1,3 ! 应力张量STRESS(3,3)和m⊗s的点积 TAU(IS) TAU(IS) STRESS(I,J)*M(I,IS)*S(J,IS) ENDDO ENDDO ENDDO“分切应力”其实就是应力在滑移方向上的“有效分力”当它的绝对值超过TAU0这个滑移系就开始干活了。我刚开始写的时候直接把TAU0设成常数结果模拟出来的应力达到屈服后就一直平着完全是理想塑性后来才想起要加硬化——晶体滑移后会产生位错堆垛后续滑移得更费力所以临界分切应力得随滑移量变大C --- 线性硬化模型 REAL*8 GAMMA(N_SLIP), H50.0D0 ! H是硬化模量MPa DO IS1,N_SLIP TAU_CR(IS) TAU0 H*GAMMA(IS) ! 当前临界分切应力 IF (ABS(TAU(IS)) .GT. TAU_CR(IS)) THEN ! 滑移系激活计算滑移率简单用粘性正则化 D_GAMMA(IS) (ABS(TAU(IS)) - TAU_CR(IS))/1.0D-3 D_GAMMA(IS) SIGN(D_GAMMA(IS), TAU(IS)) ! 保持和分切应力同符号 ELSE D_GAMMA(IS) 0.0D0 ENDIF ENDDO这里的GAMMA是累积滑移量每次计算都要更新。我试过把H设得太大模拟出来的铜硬得像钢设太小又软得离谱后来查了铜的实验数据才调到合适值——别光靠理论参数得结合实验调。晶体塑性有限元学习资料最后就是更新应力应变了晶体的塑性应变是所有激活滑移系贡献的总和C --- 计算塑性应变率张量 REAL*8 D_EPS_P(3,3)0.0D0 DO IS1,N_SLIP IF (D_GAMMA(IS) .NE. 0.0D0) THEN DO I1,3 DO J1,3 ! 塑性应变率是滑移率乘以对称化的(m⊗s) D_EPS_P(I,J) D_EPS_P(I,J) 0.5D0*(M(I,IS)*S(J,IS)M(J,IS)*S(I,IS))*D_GAMMA(IS) ENDDO ENDDO ENDIF ENDDO C --- 用胡克定律更新弹性应力弹性部分还是常规的 CALL UPDATE_ELASTIC_STRESS(STRESS, D_EPS_P, E110.0D3, NU0.34, DDTIME)我当初忘记对称化(m⊗s)了结果模拟出来的应力有非对称分量直接报错不收敛——塑性应变张量必须是对称的这步绝对不能省把这些代码拼进VUMAT再在ABAQUS里建一个单晶体拉伸试样材料设置一定要选“Crystal Plasticity”别选常规的弹塑性我第一次就犯了这错跑出来的结果跟普通材料没区别尴尬到抠脚设置好拉伸边界条件跑一下出来的应力应变曲线应该是先直线弹性阶段然后屈服之后应力随应变上升因为硬化这就对了。要是曲线一直线性那肯定是晶体塑性没生效——要么材料设置错了要么代码里的滑移系激活判断写砸了。其实学晶体塑性有限元别一开始就死啃张量推导先跑通最简单的单晶体例子再慢慢加东西从单晶体到多晶体从拉伸到复杂加载从线性硬化到考虑滑移系间的相互作用硬化。遇到问题就调参数看结果模拟的滑移带不对就查取向应力曲线不对就查Schmid因子或硬化模型。最后说句掏心窝子的晶体塑性的乐趣在于你能看到微观变形机制怎么影响宏观性能——同样是铜[001]取向和[111]取向的屈服强度差好几倍模拟曲线一眼就能看出来这种从微观到宏观的对应感比死推导公式有意思多了。