晶体取向对镍基单晶DD6合金蠕变性能的影响

贺鹏飞，吴 宸，张诚江

（同济大学航空航天与力学学院，上海200092）

提高涡轮进口温度是提升发动机功率和效率的根本技术途径，为了满足高温下的服役需求，国内外所有的先进航空发动机涡轮叶片都采用镍基单晶高温合金制造。其中，镍基单晶合金DD6由于具有良好的抗氧化，抗热腐蚀性能，以及优异的蠕变和疲劳性能，因此被用于装备国产第四代主力战机的WS-1X发动机以及装备国产大飞机C919的CJ1000AX发动机的涡轮叶片制造［1-3］。

单晶涡轮叶片在高温、高转速（高离心力）的环境下长时服役，叶片将发生蠕变变形甚至断裂破坏，危及发动机的安全。由于单晶材料的晶胞模量的各向异性，单晶合金在蠕变过程中表现出明显的各向异性与取向敏感性，对此，国内外学者已经开展了大量的研究。Han等人［4］通过对两种温度下不同取向的镍基高温合金开展蠕变试验，发现在760oC温度下［001］取向寿命最长，而［011］取向寿命最短；在1040oC下寿命最短的仍为［011］取向，而［111］取向试样的寿命最长。同样的，Chatterjee和Yu等人［5-9］分别通过蠕变试验和微观组织观测的方法，研究了镍基单晶合金蠕变寿命和取向的相关性，同样发现不同取向的单晶合金的蠕变持久性能存在差异性。很多学者将蠕变的取向敏感性归结于滑移系的开动，Li［10］对第三代镍基单晶开展了850oC条件下3种取向的蠕变试验，认为3种取向下表现出蠕变各向异性主要与所激活的滑移系有关，［111］取向开动的是＜110＞｛111｝滑移系，而［001］取向首先开动的则是＜112＞｛111｝滑移系。此外，Latief等人［11-15］提出，镍基单晶合金的蠕变各向异性在中温（750oC）条件下比高温（900oC）条件下更显著，这主要和不同取向下应力的不同和开动的滑移系的不同有关。在有限元模拟方面，西北工业大学岳珠峰课题组［16-20］基于晶体滑移理论建立了单晶的晶体取向相关性模型，认为在蠕变第一阶段晶格错配度、热不协调性和微观结构的定向特点这3个因素影响了该阶段镍基单晶合金的取向相关性，而这种取向相关性又直接导致晶体在不同取向下寿命和力学性能有显著差异。

目前对单晶高温合金不同取向下蠕变性能差异的机理与寿命模型间的关联尚未完全建立。本文基于第二代镍基单晶高温合金DD6，采用试验与有限元分析相结合的方法，研究［001］、［011］与［111］3种典型取向试样在980oC下的蠕变性能差异，并对各取向微观组织进行表征，系统的研究了镍基单晶合金的微结构演化对蠕变性能的影响。基于试验观测所得到的结果，揭示蠕变变形与损伤机制，建立基于晶体塑性理论的各向异性蠕变粘塑性本构模型，同时编写Abaqus有限元软件的Umat子程序实现模拟工作。

1 DD6单晶合金蠕变试验方案

1.1 试验设备

本项目的蠕变试验于陕西省航空发动机结构强度与可靠性重点实验室完成，所采用的试验机为RDL100微控电子式蠕变持久试验机，采用恒定力控制，同时采用上中下三段式加热使高温炉内温度均匀分布，使试样受热均匀，每段配备Q型热电耦进行测温，达到校正试件温度的目的。试验采用接触式高温引伸计测量试样标距段的应变，获取完整蠕变曲线。

1.2 试验材料与试件

DD6单晶高温合金主要成分见表1，其铸造过程为：真空条件下将严格配比的原材料熔炼为DD6母合金，其后在定向凝固炉上采用籽晶法将母合金定向浇注为［001］、［011］和［111］3个取向的单晶试棒，再经过冷却凝固，最后对铸态的单晶进行热处理，其条件为1 290oC×1h+1 300oC×2h+1 315oC×4h/空 冷+1 120oC×4h/空冷+870oC×32h/空冷。

表1 DD6镍基单晶高温合金主要化学成分Tab.1 Main chemical compositions of DD6 nickelbased single crystal superalloy

镍基单晶合金DD6的微观组织是由宽度为0.5～1μm的立方状γ’强化相及分布在其周围的γ基体相组成。经过电化学腐蚀的微观组织形貌如图1所示，其中γ相为无序的面心立方结构的Ni固溶体，包含Re、Mo、W等大量难溶元素，起固溶强化作用；而γ’相为有序的L12型Ni3Al面心立方结构，起沉淀强化作用。γ’相以共格的方式镶嵌于基体相中，呈规则的正方形状排列，γ’相有较好的稳定性，是维持合金高温下优异力学性能的重要保障，对合金起到主要的强化作用。

为完成不同取向下DD6单晶的蠕变试验并获取应变-时间曲线，本研究依照GB/T 2039-2012设计圆棒试样，如图2所示。其中用来测量应变的标距段长度为19mm，直径为4mm。

图2 980oC/300MPa蠕变试件尺寸（单位：mm）Fig.2 Creep specimen size at 980oC/300MPa（unit：mm）

1.3 试验方法

每个取向采用两件试样开展试验。将试样温度升至980oC后保温30min，后增加载荷至目标载荷300MPa，保持总拉力恒定不变，直至试样断裂。试验结束后，绘制蠕变曲线，拍摄断口形貌，选取试样特定位置取样切块，打磨抛光后用20gCuSO4+5ml H2SO4+100ml HCL+80ml H2O配比的腐蚀液对试样表面进行腐蚀，采用ZEISS场发射电子扫描显微镜（SEM）观察试样的微观组织形貌。

2 DD6单晶合金蠕变试验结果

2.1 蠕变寿命与蠕变曲线

试验获得的不同取向下的蠕变寿命与延伸率如表2所示。由试验结果可知DD6单晶合金在980oC下3种典型取向试样的蠕变寿命与延伸率均为［111］取向＞［001］取向＞［011］取向。

表2 不同取向下的蠕变寿命与延伸率Tab.2 Creep life and elongation under different orientations

图3为DD6合金在980oC及300MPa下［001］、［011］和［111］取向试样的蠕变变形曲线，表明DD6合金的蠕变过程具有明显的取向相关性。3种取向的试样在980oC下都没有明显的减速蠕变阶段。［111］取向试样蠕变稳态阶段的时间较长，加速蠕变阶段的伸长率最长，总寿命最长，接近180h；［011］取向试样蠕变稳态阶段较短，合金迅速进入加速蠕变阶段，表现出较高的蠕变速率，并迅速产生断裂，寿命仅为40h左右；［001］取向试样在蠕变稳态阶段中，前30h变形量与［011］取向试样基本一致，但在30h之后较［011］取向试样变形速率更缓慢，后进入加速蠕变阶段，直至断裂，寿命为100h左右。综上所述，［111］取向蠕变稳态时间最长，寿命与延性最优；［011］取向的蠕变稳态时间相对最短，其寿命与延性最低；［001］取向的寿命与延性处于［111］与［011］取向之间。

图3 3种取向下的蠕变曲线图Fig.3 Creep curves under three orientations

2.2 不同取向DD6合金蠕变断口形貌

图4a～4c分别为3种取向下蠕变断口的解理面形貌，3种取向下的试样断口解理面呈现较大差异，其中［001］方向的断口有大量规则且相互平行的正方形解理面（小刻面），该断裂平面垂直于应力方向，其解理平面边缘存在撕裂棱，解理面中心即为微空洞形成处，如图4a所示箭头所指处，孔洞区域为该解理面裂纹源开动处，随着蠕变的进行，由孔洞附近的微裂纹逐渐扩展至解理面边缘。可知每一个凹陷的立方状或者球状的解理面位于断口的不同高度处，故而表明，孔洞产生处位于不同平面，进而，其解理断裂也在不同平面处发生，因此最终断裂时产生了不同高度间的断裂台阶，使得该取向的蠕变断裂呈现类解理断裂特征。

［011］取向的解理面如图4b所示。裂纹沿不同晶面扩展，形成不同平面的解理面和解理台阶，台阶高度大约与1～2个强化相的宽度接近，与其他两个取向不同的是，该取向的解理台阶较为平坦，且都平行于椭圆断面的长轴。图中上半部分的台阶较密集和细小，而下半部分的台阶较平坦稀疏，层片状区域较宽大，表明随着蠕变的进行，裂纹由断口上部向下逐渐扩展。箭头所指处为该断面的孔洞和裂纹源产生的位置。由图5b断口整体形貌可知，该取向下试样断裂后发生较大变形，断口呈椭圆形。

图5 3种取向下断裂试样示意图Fig.5 Schematic diagram of fractured sample under three orientations

［111］取向下的单晶合金蠕变断裂后断口形貌如图4c所示。该取向断口表面粗糙，分布大量清晰平整的正三角形的解理面，解理面与解理面之间交线形成的棱清晰可见，且相互平行，形成一个个正三角形的韧窝。与［001］取向类似，［111］取向解理面的微空洞被包围处于构成韧窝的3个面的中心，如箭头所示，这些在解理面与解理面夹角的中心处的孔洞随着蠕变的进行，逐渐产生微裂纹，进而再随着蠕变时间的延长，发生裂纹扩展，形成这种正三角形特征的解理面，并导致试样最终发生断裂，该取向也表现出了明显的解理断裂特征。

图4 3种取向试样断口的局部特征图Fig.4 Local characteristics of fracture of three orientations

3种取向下的试件伸长量如图5所示。单晶试样材料的伸长率与不同取向下材料的微观结构的破坏机理相关，DD6合金［001］和［111］取向的断口都有较为明显的韧窝状解理形貌，韧窝中心是孔洞和裂纹源产生的位置，从微空洞或者微裂纹开始，向外扩展成解理面，解理面和解理面之间相互连接形成韧窝状的断裂形貌，最终发生断裂。在整个扩展过程中，由于内部的微空洞和微裂纹的逐步扩展受阻，需要更多能量，导致试样积累了大量的塑性变形，因此［001］和［111］取向下的蠕变第二阶段更长，材料的应变更大，伸长率更高。而［011］取向断面平整，断口呈椭圆形，表明该取向下参与变形的滑移系较少，且变形有轴向不对称性的特征［21］。有文献通过（EBSD）电子背散射衍射技术［22-23］，测量［011］取向断裂后试样的晶格旋转方向，研究发现该取向蠕变时晶格先沿［001］－［111］边界方向旋转然后转动至［111］极点直到试样断裂。这种晶格扭转畸变使得裂纹沿椭圆长轴方向迅速扩展，很快贯穿整个截面，因此呈现一个台阶状的解理面，累计的塑性变形较少，最终导致该取向下试样的伸长率较低，稳态蠕变时间较短，寿命也相应较短。

这些不同的解理面表明，不同取向下试样断裂时孔洞产生的位置，空洞附近的裂纹的扩展方向不同，使得这些不同取向下的试样断口的解理面形态各异，最终导致蠕变寿命与延伸率的不同。

2.3 不同取向DD6合金蠕变筏化形貌

为观测蠕变后不同取向的单晶合金微结构演化，采用图6所示位置进行取样。取标距段点A和点B处取样用以观测其靠近断口和远离断口的微观组织两相结构，取点C处取样用以观测非标距段处（低应力区域）的微观两相组织结构。试样组织的观察面为纵截面，即与蠕变载荷方向平行。

图6 微观形貌观测取样位置Fig.6 Sampling position for microscopic morphology

图7为DD6合金［001］、［011］与［111］取向试样的不同区域在扫描显微镜下的两相结构的微观组织形貌，化学腐蚀后的黑色区域为强化γ’相，主要成分是Ni3Al，白色区域为基体γ相。有文献指出，γ’相的体积分数对合金的性能也有较大影响，当γ’相体积分数较高时，该相稳定性会降低且易于连接，进而导致蠕变速率增加，蠕变寿命缩短［24］。

图7a～7c为［001］取向试样不同位置的微观组织形貌。不同位置筏化的差异是由应力的不同造成的。图c为取样点C处，横截面面积较大，因此应力更小，约为192MPa，与图7a、图7b相比，该处γ’相基本呈规则立方体状，但此时已经有部分立方状的γ’相相互连接，基体通道与蠕变前原始形貌相比变宽，形貌较不规则，粗细不一。在取样点B处，如图7b所示，其横截面面积比非标距段小，应力大，该区域可以观察到大量的γ’相在高温和应力作用下沿横向筏化，形成了长条状的筏状组织。与此同时，［001］方向基体通道被切断，［010］方向的基体相呈长条状，零星分布着点状的基体相结构。图7a所示区域为断口附近的取样点A，由于该区域发生颈缩变形，横截面积急剧减小，受到的应力最大，立方状γ’相发生了较大的形态的变化，与图7b相比，该区域已经到了筏化完成后的解筏阶段，基体通道的宽度大致相同，大部分的［001］和［010］方向的γ’相相连接，随着筏化阶段的完成，可以看到长条状的基体相断裂，形成短粗状的结构。

图7d～7f为［011］取向试样不同位置的微观组织形貌。与［001］取向类似，由于不同取样点处的应力不同，各位置的筏化程度也不相同。其中，图7f为取样点C的微观组织形貌，在［011］取向下的γ’相结构的晶轴坐标与拉伸方向呈45°角，立方状的γ’相排列均匀，只有极少量的γ’相相互连接，少部分浅色的基体通道断开。在取样点B处，截面积小，使得该处应力大，筏化程度更高，部分立方状的γ’相在拉应力的作用下形成与［001］方向呈45°角的筏状结构，如图7e中箭头所指方向，其垂直于筏化方向的γ’相相互连接尺寸变大，只有少量的γ’相还保持立方状结构。整体来看，筏状的γ’相结构相互连接，夹角大约为90°，少量的点状的γ相分布在深色的γ’相中，但与［001］取向相比，点状基体相更少。图7d所示为［011］方向取样点A的纵截面微观组织形貌。与图7e相比，该区域的筏化方向更加明显，如图中箭头所指方向。箭头方向上的γ’相连接形成筏形结构，基体通道沿箭头方向连接成更长的条带状，几乎看不到垂直箭头所指方向的γ相结构，且γ相基体通道变宽，粗化现象明显，且只有极少量的点状γ相分布于γ’相中间。

图7g～7i为［111］取向试样不同位置微观组织形貌。可以看到该取向下取样点C的两相组织形貌与［001］和［011］不同，其粗化现象更加明显，这是由于［111］取向的蠕变寿命更长所导致，γ’相沿箭头所指方向排列，形筏取向与应力轴方向接近35°角，且部分横向和纵向的基体相相互连接形成新的γ’相，基体通道也不再是排列整齐的交叉结构，部分横向和纵向的基体通道断开，基体通道粗化程度较为明显。图7h所示取样点B处的微观形貌中呈立方状的γ’相消失，白色的相互连接的基体通道发生断裂且逐渐变粗，部分基体通道相互连接呈波浪状的筏形结构，γ’相形筏取向与应力轴方向接近35°角。图7g所示取样点A位置的微观组织图像，在宏观上有明显的颈缩现象，该处横截面积减小，所受应力最大，所以筏化现象也更加明显，相互连接的基体通道断裂溶解在γ’相中，大部分相互连接的基体通道都发生了断裂，形成更加短粗状的形貌，可以看到在强化相中还溶解了部分点状的基体相组织，此时γ’相相互连接，尺寸增大，筏化程度也逐渐进入到解筏阶段。

图7 ［001］、［011］和［111］取向断口纵截面微观结构形貌Fig.7 Microstructure of longitudinal section of[001][011]and[111]orientation fracture

3 基于晶体塑性理论的粘塑性本构方程

单晶的蠕变变形可归结于滑移系的开动，基于一个晶体塑性理论的模型［25-26］，将各滑移系的蠕变剪应变率用̇(α)表示出来，其中α代表不同的滑移系，A、n是与温度大小有关的参数，则

τ(α)为滑移系α的分切应力，可表示为

式中：P(α)是取向因子，可表示为

式中：m(α)和n(α)是两个不同的向量，分别表示滑移系开动的方向和该滑移系对应的滑移面单位法向量的方向。

用̇表示宏观应变率，已知在蠕变变形中可分弹性变形和非弹性变形两种，那么用̇表示弹性应变率，用表示非弹性应变率，可得

弹性部分的应力应变关系可用胡克定律来表示，非弹性部分则可得

如果将蠕变变形用不同滑移系的蠕变变形来表示，则可得

由于镍基单晶高温合金特殊的两相组织形式，对于合金蠕变的强化作用主要来源于两方面，一方面是在基体相和强化相中间形成的相界面对位错产生的阻碍作用，把这种相界面处的位错网络对蠕变产生的阻碍应力用表示［27］，Zhao等人［28-30］在对高温合金的蠕变行为的研究中，认为合金的蠕变强化作用与其位错的密度的平方根成正比，即

式中：b是伯格斯矢量的模；G为剪切模量；c2是一个常数；表示γ基体相不同取向的位错密度，其演化规律为

式中：k1，k2都是材料常数。

另一方面在蠕变过程中，强化相对于蠕变中的位错运动也有阻碍作用，位错要越过强化相需要一个较大的门槛应力，称为Orowan应力，记为τor。蠕变阶段位错需要克服强化相的阻碍作用，这种行为被称为Orowan绕过机制。Tinga等学者［31-33］认为这种Orowan应力与基体通道的宽度（即强化相间距）成反比，表示为

式中：κ表示基体通道的宽度；G是剪切模量；b表示伯格斯矢量的模；λ是材料常数。通过对蠕变中断试验后得到的处于筏化不同阶段的试样显微组织的分析中不难看出，基体通道的粗化程度与蠕变时间呈显一个抛物线关系，表示为

式中：κ0为初始基体通道宽度；c1为表征材料筏化速率的参数；可由微结构观测得到；t为蠕变时间。

高温合金材料在拉伸时发生破坏的方式是直接瞬断还是长时蠕变断裂和外部荷载大小与材料屈服应力有一定关系。张诚江等人［1］运用一个判据定义了材料发生瞬断还是蠕变断裂的失效判据，即

可分为以下3种情况：

（1）当R(α)≤0时，即分切应力非常小，小于位错进入相界面位错网络的临界值时，此时令

（2）当R(α)＞0时，可知，此时材料分切应力小于临界分切应力，进入蠕变变形的状态。基于Kachanov和Ravbot-nov［34-35］提出的连续损伤模型加上Yeh等［36］提出的损伤演化率的基础上，建立了分切应力与剪切应变率有关的蠕变损伤模型。其中为初始蠕变率；χ和φ跟n一样都是与温度相关的参数；β是取值为2.5的常数；ω表示材料的损伤；为稳态蠕变率，由蠕变第二阶段的速率得到；为初始损伤率，其中是温度和应力的函数，即

采用葛庭燧-Dorn（K-D）公式与Arrhenius表达式相结合的形式，则

式中：T为绝对温度；R为气体常数；Q为激活能。

（3）当R(α)≥0且τ(α)≥τc时，也就是，此时分切应力大于材料的临界分切应力，材料将发生瞬时拉伸破坏，位错直接越过相界面的位错网络和强化相，上述所推出的蠕变本构方程和损伤方程将不再适用。

可以将式（13）从0～1的积分，得到损伤从初始的0～1时的断裂时间，即推出材料的寿命预测模型为

式中：N为某滑移系下的滑移面开动方向的个数，八面体滑移系（Oct1，＜110＞［111］）与十二面体滑移系（Oct2，＜112＞［111］）的N为12，六面体滑移系（Cub，＜110＞［100］）的N为6。式（1）～（15）中涉及的弹性参数与蠕变参数分别见表3与表4。通过蠕变过程中位错的观测［1］，最终确定了不同取向下滑移系的开动方式，为［001］与［011］取向为八面体滑移系开动，而［111］取向为六面体滑移系开动。

表3 980℃下DD6镍基单晶高温合金弹性参数Tab.3 Elastic parameters of DD6 nickel based single crystal superalloy at 980℃

表4 980℃下DD6镍基单晶高温合金蠕变参数Tab.4 Creep parameters of DD6 nickel based single crystal superalloy at 980℃

4 模拟讨论

单晶合金在蠕变时随蠕变时间的增长，一般呈现出3个不同的阶段，即蠕变第1阶段，蠕变第2阶段和蠕变第3阶段。在蠕变第1阶段应变量持续增加，但蠕变速率不断降低，材料发生硬化，又称减速蠕变阶段；而在蠕变第2阶段，蠕变速率维持在最小值，并保持相当长的一段时间，又称稳态蠕变阶段；在蠕变进行到第3阶段时，伴随材料的颈缩，蠕变速率迅速上升，蠕变变形急剧增长，随后试样发生断裂破坏，又称蠕变加速阶段。对于DD6单晶合金材料来说，在高于1 000oC的蠕变试验过程中，蠕变的3阶段特征较为明显［37］；而低于1 000oC蠕变过程中蠕变第1阶段不明显［38］。因此该试验条件下选取的3个取向的试样均只呈现出蠕变的第2、第3阶段。

把蠕变损伤本构模型通过ABAQUS有限元软件，并将用户子程序导入实现模型的计算，得到如图8所示的3种取向980°C/300MPa条件下的模拟曲线。由3个取向下的模拟取向和试验曲线的对比可以看到，模拟曲线和试验曲线较为吻合，这说明该3种取向下的蠕变损伤本构模型理论基本符合实际情况。图9为3个取向在蠕变第3阶段的变形模拟。

图8 3种取向980°C/300MPa试验曲线与模拟曲线Fig.8 Test and simulation curves for three orientations at 980°C/300MPa

由图9的有限元仿真的应变模拟结果与图5的真实蠕变断裂位置对比可知，除了寿命的模拟结果与试验一致，每个取向下的圆棒试验件的蠕变变形模拟结果与破坏位置均能够与蠕变试验一一对应。具体来说，由于变截面处的多轴应力状态，［001］取向的颈缩截面在标距段的1/4处，［111］取向的颈缩截面在标距段的1/6处，而［011］取向的颈缩截面在标距段1/2处，断口呈椭圆形，断面光滑平整，与拉伸方向呈45°，试验结果与模拟一致。由于试样的微观组织形貌与单晶合金的晶体取向、应力和蠕变的时间均有关，其中，晶体取向决定了筏化的方向，应力大小决定了筏排化的速度，蠕变的时间决定了远离标距段位置的筏排化程度。而蠕变的取向敏感性主要是因为外部载荷分解到滑移面上的分切应力大小不同，导致沿滑移方向的变形速率以及损伤累计的速率不同。综上，说明基于晶体塑性理论的粘塑性本构方程能够模拟由于蠕变带来的大变形，能够适用于真实叶片服役过程的蠕变变形与寿命预测。

图9 3种取向980°C/300MPa有限元仿真结果Fig.9 Finite element simulation results of three orientations at 980°C/300MPa

5 结论

（1）对3种取向下的试验件进行了980°C/300MPa条件下的蠕变试验，［011］取向合金的蠕变第二阶段的蠕变速率最高，寿命最短；［001］取向合金的蠕变速率和寿命居中；［111］取向合金的稳态蠕变速率也较低，蠕变寿命最长。同时晶体取向与试样的蠕变伸长率也有一定关系，［001］和［111］取向的蠕变伸长率较高，延性较好，而［011］取向蠕变伸长率较低。

（2）通过对［001］、［011］和［111］3个取向试样980°C/300MPa下的蠕变断口形貌观测，［011］取向的颈缩现象较明显，断口呈椭圆形，整体呈台阶状的层片形，层与层之间有明显的解理台阶；［001］和[111]取向呈明显的解理断裂特征，其中，［001］方向的断口有正方体形或者圆形的解理面；［111］取向断口则分布大量的正三角形的解理面。

（3）在SEM扫描电镜观测下，可以看到3种取向下试样的非标距段部分纵截面的两相结构未发生明显的筏化现象。［001］取向远离断口处纵截面两相结构呈长条状的筏状组织，而靠近断口处已经处于解筏阶段；［011］取向远离断口处纵截面两相结构中立方状的γ’相与［001］方向呈45°角的筏状结构，而靠近断口处沿45°角的筏化趋势更为明显；［111］取向远离断口处纵截面两相结构可以看到γ’相形筏取向与应力轴方向接近35°角，靠近断口处沿该方向的筏化趋势更加明显，且基体通道断裂，进入解筏阶段。

（4）采用试验和有限元相结合的方式，基于原有的晶体塑性理论下的蠕变本构模型，建立了考虑晶体各向异性特征的蠕变本构模型和蠕变寿命预测模型，通过Abaqus子程序实现有限元的模拟工作，并与试验结果进行对比，发现模型能较好的反映试样的蠕变寿命和变形，对于工程实际问题有较大的参考价值。

作者贡献说明：

贺鹏飞：提供研究的构思，设计，及方案的提出，论文审阅；

吴宸：参与论文的试验内容，数据的收集，模型的建立，起草论文；

张诚江：数据分析，论文结果讨论，论文的审阅及定稿。