一种多相复合材料力学行为的多尺度模拟方法，属于计算材料学、模拟仿真、高通量计算等多学科交叉领域。通过纳观尺度的第一性原理计算、微观尺度的分子动力学模拟、介观尺度的热力学计算可以获得宏观尺度的有限元模拟所需的各种物性参数，包括复合材料中各个相在不同温度、不同晶粒尺寸下的弹塑性物理参数等。利用聚焦离子束实验和图像处理引入获得真实的材料组织结构。通过各个尺度计算结果之间的参数耦合与参数传递，结合材料的组织结构，可以模拟多相复合材料在复杂应力和不同温度下的应力应变关系、应力分布及其演变、塑性变形等力学行为。
1.一种多相复合材料力学行为的多尺度模拟方法，其特征在于，包含以下步骤：(1)纳观尺度的第一性原理计算
纳观尺度的第一性原理计算可以分别获得多相复合材料中单独的金属相和陶瓷相在
0K下的弹性性质以及不同体积的晶体结构对应的能量，计算步骤如下：
1.1对0K下的多相复合材料中的金属相和陶瓷相的杨氏模量、体弹模量、剪切模量、泊松比进行计算：首先分别对金属相和陶瓷相材料的单晶结构进行结构弛豫；利用弛豫后的晶体结构对材料的弹性矩阵C的矩阵元，即弹性常数，进行计算；
其中σiεicij(i,j＝1,2,3,4,5,6)分别为应力、应变和单晶材料的弹性常数；
进而利用Voigt-Reuss-Hill近似，将计算得到的弹性常数转化为杨氏模量、体弹模量、剪切模量和泊松比；
B＝(Bv+Br)/2 (2)G＝(Gv+Gr)/2 (3)E＝9B·G/(3·B+G) (4)υ＝0.5·(3B-2G)/(3B+G) (5)其中B为体弹模量、G为剪切模量、E为杨氏模量、υ为泊松比；泊松比近似认为不随温度变化；体弹模量、剪切模量、杨氏模量会随温度变化而变化，随温度变化规律的计算方法将在步骤(2)中详述；式(2)(3)中下角标v和r分别代表Voigt近似和Reuss近似计算方法；
Bv＝(c11+c22+c33)/9+2(c12+c23+c13)/9 (6)Gv＝(c11+c22+c33-c12-c13-c23)/15+(c44+c55+c66)/5 (7)Br＝1/(s11+s22+s33+2s12+2s13+2s23) (8)Gr＝15/(4s11+4s22+4s33-4s12-4s13-4s23+3s44+3s55+3s66) (9)其中sij(i,j＝1,2,3,4,5,6)为材料的柔度矩阵S的矩阵元，即柔度系数；柔度矩阵S与弹性矩阵C互为逆矩阵；
1.2对0K下的多相复合材料中的金属相和陶瓷相的晶体体积-能量关系分别进行计算：将步骤1.1中弛豫后的金属相和陶瓷相的晶体结构的晶胞体积分别进行-10％至+10％范围内的收缩或者膨胀，计算对应不同晶胞体积的能量；
介观尺度的热力学计算，基于准谐德拜模型，可以分别获得多相复合材料中的金属相和陶瓷相在不同温度下的弹性性质和热膨胀系数，计算步骤如下：
2.1对不同温度下的多相复合材料中的金属相和陶瓷相的杨氏模量分别进行计算不同温度下的杨氏模量可以通过不同温度下的体弹模量计算得到，而计算不同温度下的体弹模量需要首先计算非平衡态Gibbs自由能；
G*＝E(V)+pV+Avib(θ,T) (10)其中E(V)表示步骤1.2中计算得到的晶胞总能量，p和V分别为外部压强和晶胞体积，T为温度，Avib为振动Helmholtz自由能，θ为德拜温度；振动Helmholtz自由能Avib可表示为：其中D(θ/T)为德拜积分， n为每个晶胞的原
其中M为晶胞质量，υ为泊松比， 为约化普朗克常数，f(υ)为泊松比函数Bs为绝热体弹模量，可使用步骤1.2中获得的体积-能量关系数据获得，其计算方法如下：
对于某一给定的(p,T)，非平衡Gibbs自由能对体积求极小值，即：
由式(14)可获得体系热力学状态，即获得体系的普适状态方程，明确体系的p、T、V三者的关系；进一步的，某一温度的体弹模量可表示为：通过该步计算可将温度效应引入体弹模量，计算特定温度下的体弹模量BT；
ET＝3BT(1-υ) (16)可以获得不同温度下的杨氏模量；利用公式(16)计算得到的0K下的杨氏模量ET,0与步骤1.1计算得到的0K下的杨氏模量E存在差值ΔE＝ET,0-E；此差值作为利用公式(16)计算得到的其他任意温度的杨氏模量的零点修正项，以此计算任意温度的杨氏模量E(T)＝ET-ΔE；
2.2对不同温度下的多相复合材料中的金属相和陶瓷相的热膨胀系数分别进行计算从另一方面，恒压热容Cv和Gruneisen参数γ可分别描述为：基于式(15)、(17)、(18)，可获得热膨胀系数α(T)：
由此，可将温度效应引入热膨胀系数，计算不同温度下的热膨胀系数；
微观尺度的分子动力学模拟可以获得多相复合材料中的金属相在特定的温度、特定晶粒尺寸下的塑性性质；硬脆陶瓷相的塑性相比金属相极其微小，可忽略不计；
首先通过LAMMPS软件获得多相复合材料中金属相的单晶模型，随后通过AtomEye软件中的Vironoi方法构造相对窄范围内的不同晶粒尺寸的多晶模型，相对窄范围内的晶粒尺寸为30到40nm之间；对此相对窄范围内的晶粒尺寸多晶模型在特定温度下进行压缩模拟，绘制应力应变曲线，从应力应变曲线读取屈服强度和硬化系数；假定硬化系数与晶粒尺寸无关，为了减少误差，针对所有晶粒尺寸的多晶模型的硬化系数采用多个不同小晶粒尺寸的硬化系数的平均值；
对于任意不同晶粒尺寸模型的屈服强度σs，利用上述从应力应变曲线读取的相对窄范围内的多晶模型的屈服强度数据和Hall-Petch关系进行拟合；
σs＝σ0+kd-1/2 (20)其中σs为屈服强度，σ0为单个位错移动时晶格摩擦阻力，k为常数，d为可取的任意晶粒尺寸；
拟合过程中，σ0采用派纳力τp的计算结果，σ0≈τp；派纳力的取值为金属相的主滑移系的派纳力数值，计算公式如下：其中υ为泊松比，取步骤1.1中的计算结果；a为滑移面的面间距，b为滑移方向上的原子间距，可以利用步骤1.1中弛豫的晶体结构通过几何关系计算得到；G(T)为特定温度下的剪切模量，取步骤2.1中的计算结果；
(4)构建多相复合材料的真实微观组织的三维有限元模拟几何模型
通过聚焦离子束双束系统，利用离子束对多相复合材料加工获得边长为微米级的长方体，以20nm的间隔对实验区域进行断层扫描；通过聚焦离子束实验获得图像堆栈后，采用Avizo软件进行图像处理，图像处理过程如下：首先应用FIB Stack Wizard模块对图像进行对齐、剪切、变形修正和阴影修正，进而采用中值滤波器和保持边界平滑的高斯滤波器进行去噪处理，最后，采用阈值分割方式对图片进行二值化处理，区分图像中金属相和陶瓷相，获得可用于构建有限元几何模型的复合材料组织图片；然后构建有限元几何模型；
(5)利用有限元模拟方法对多相复合材料力学行为进行模拟
根据步骤(1)中计算得到的泊松比、步骤(2)中计算得到的特定温度下的金属相和陶瓷相的杨氏模量与热膨胀系数，步骤(3)中计算得到的特定温度下特定晶粒尺寸的金属相的屈服强度和硬化系数，构建成陶瓷相和金属相的特定温度有限元参数模型；
读取步骤(4)中构建的复合材料组织的有限元几何模型，定义边界条件和约束条件，将陶瓷相和金属相的特定温度有限元参数模型赋予到复合材料组织的有限元几何模型中进行求解。
2.按照权利要求1所述的一种多相复合材料力学行为的多尺度模拟方法，其特征在于，步骤(4)构建复合材料组织的有限元几何模型过程如下：(a)通过Matlab提取上述所述的复合材料组织图片中金属相所有像素的坐标即金属相像素坐标；(b)在有限元软件中建立复合材料组织图片的有限元几何模型即总有限元几何模型，即将金属相和陶瓷相都等同是陶瓷相，建立有限元几何模型；并获取总有限元几何模型中所有单元的坐标；(c)然后结合(b)中总有限元几何模型中所有单元的坐标，根据步骤(a)中的金属相像素坐标确定总有限元几何模型所有单元中应属于金属相的单元，并将该类单元的材料属性修改为金属相材料属性，完成复合材料组织的有限元几何模型建立。
3.按照权利要求1所述的一种多相复合材料力学行为的多尺度模拟方法，其特征在于，步骤(5)还包括以下步骤：定义边界条件和约束条件，然后对复合材料组织的有限元几何模型在不同温度下进行复杂应力的加载，进一步量化分析复合材料在不同温度、不同复杂应力下的应力应变关系、应力分布及其演变规律、塑性变形。
一种多相复合材料力学行为的多尺度模拟方法
技术领域
[0001] 本发明属于计算材料学、模拟仿真、高通量计算等多学科交叉领域，具体涉及基于材料真实组织结构的适用于描述多相复合材料力学行为的计算机模拟分析方法。
背景技术
[0002] 在“材料基因组计划(MGI)”的革命性思维的引领下，材料计算在指导新材料设计研发中的作用越来越重要，可为新材料设计开发提供准确的指导，有效改变传统材料设计开发的“试错”过程，大大提升新材料的研发效率，从而有效推动高性能新材料开发及相关先进制造业的快速发展。目前常用的材料计算模拟方法包括纳观、微观尺度的第一性原理计算和分子动力学模拟等方法，可以对材料成分和结构进行筛选设计；宏观的有限元模拟等方法能够对材料的服役行为和使用性能进行计算分析。然而，目前国内外公认的普遍问题是，不同尺度下计算模拟方法所用的计算模型、典型算法、输入输出数据和参数各不相同，虽然各个尺度下的计算可能有丰富的方法和工具，但是严重缺乏跨尺度、多尺度的计算模型和模拟方法；不同尺度的计算方法之间很难进行数据传输、无法实现同一类参量在不同尺度之间数据的融合贯通。因此，当前迫切需要开发材料成分-组织-工艺-性能之间内禀关联的多尺度计算模型，以此实现材料全链条一体化设计和高效、准确的材料性能预测和分析。
[0003] 多相复合材料有着广泛的工程应用，粉末冶金方法是制备该类材料的典型途径。对于多相复合材料而言，由于各相之间的物性参数差异较大，尤其，热膨胀系数的差异会造成粉末烧结过程中产生较大的热残余应力，而制备态热应力会对材料性能产生很大影响。
进一步，热应力与外加应力的交互作用会极大地影响材料服役过程中的力学行为和使用性能。在多相复合材料中，通常会出现应力分布不均匀的现象，服役过程中各相内会出现差异较大的应力响应。结合真实微观组织的三维有限元模拟能够准确揭示相内不均匀应力响应规律，对材料力学行为与性能的研究具有重要意义。然而三维有限元模拟需要大量材料物性参数，很多参数缺乏实验数据，因此构建多尺度计算模型和算法，结合其他尺度的计算方法，为有限元模拟提供参数数据，可解决单一有限元模拟数据匮乏的问题。而且，由于多尺度模型涉及纳观、微观、宏观多种计算方法，能够同时完成对材料成分筛选、组织结构优化和使用性能的预测分析，因此开发基于真实材料组织结构的多相复合材料力学行为的多尺度模拟计算方法具有极为重要的意义和应用价值。
发明内容
[0004] 本发明提供了一种基于材料真实微观组织，对多相复合材料的力学行为进行多尺度模拟的方法。通过纳观尺度的第一性原理计算、微观尺度的分子动力学模拟、介观尺度的热力学计算可以获得宏观尺度的有限元模拟所需的各种物性参数，包括复合材料中各个相在不同温度、不同晶粒尺寸下的弹塑性物理参数等。通过各个尺度计算结果之间的参数耦合与参数传递，可以模拟多相复合材料在复杂应力下的应力应变关系、应力分布及其演变规律、塑性变形等力学行为。
[0005] 为实现上述目的，本发明一种多相复合材料力学行为的多尺度模拟方法，其特征在于，包含以下技术方案和步骤：
[0006] (1)纳观尺度的第一性原理计算
[0007] 纳观尺度的第一性原理计算可以分别获得多相复合材料中金属相和陶瓷相在0K下的弹性性质以及不同体积的晶体结构对应的能量，计算方案如下：
[0008] 1.1对0K下的多相复合材料中的金属相和陶瓷相的杨氏模量、体弹模量、剪切模量、泊松比进行计算：
[0009] 首先分别对金属相和陶瓷相材料的单晶结构进行结构弛豫。利用弛豫后的晶体结构对材料的弹性矩阵C的矩阵元，即弹性常数，进行计算。
[0010]
[0011] 其中σiεicij(i，j＝1，2，3，4，5，6)分别为应力、应变和单晶材料的弹性常数；
[0012] 进而利用Voigt-Reuss-Hill近似，将计算得到的弹性常数转化为杨氏模量、体弹模量、剪切模量和泊松比。
[0013] B＝(Bv+Br)/2 (2)
[0014] G＝(Gv+Gr)/2 (3)
[0015] E＝9B·G/(3·B+G) (4)
[0016] υ＝0.5·(3B-2G)/(3B+G) (5)
[0017] 其中B为体弹模量、G为剪切模量、E为杨氏模量、υ为泊松比；泊松比近似认为不随温度变化；体弹模量、剪切模量、杨氏模量会随温度变化而变化，随温度变化规律的计算方法将在步骤(2)中详述；式(2)(3)中下角标v和r分别代表Voigt近似和Reuss近似计算方法；
[0018] Bv＝(c11+c22+c33)/9+2(c12+c23+c13)/9 (6)
[0019] Gv＝(c11+c22+c33-c12-c13-c23)/15+(c44+c55+c66)/5 (7)
[0020] Br＝1/(s11+s22+s33+2s12+2s13+2s23) (8)
[0021] Gr＝15/(4s11+4s22+4s33-4s12-4s13-4s23+3s44+3s55+3s66) (9)
[0022] 其中sij(i，j＝1，2，3，4，5，6)为材料的柔度矩阵S的矩阵元，即柔度系数；柔度矩阵S与弹性矩阵C互为逆矩阵；
[0023] 1.2对0K下的多相复合材料中的金属相和陶瓷相的晶体体积-能量关系分别进行计算：
[0024] 将步骤1.1中弛豫后的金属相和陶瓷相的晶体结构的晶胞体积分别进行-10％至+10％范围内的收缩或者膨胀，计算对应不同晶胞体积的能量；
[0025] (2)介观尺度的热力学计算
[0026] 介观尺度的热力学计算，基于准谐德拜模型，可以分别获得多相复合材料中的金属相和陶瓷相在不同温度下的弹性性质和热膨胀系数，计算方案如下：
[0027] 2.1对不同温度下的多相复合材料中的金属相和陶瓷相的杨氏模量分别进行计算[0028] 非平衡态Gibbs自由能可描述为：
[0029] G*＝E(V)+pV+Avi6(θ，T) (10)
[0030] 其中E(V)表示步骤1.2中计算得到的晶胞总能量，p和V分别为外部压强和晶胞体积，T为温度，Avib为振动Helmholtz自由能，θ为德拜温度；振动Helmholtz自由能Avib可表示为：
[0031]
[0032] 其中D(θ/T)为德拜积分， n为每个晶胞的原子数，kB为玻尔兹曼常数，德拜温度可表示为：
[0033]
[0034] 其中M为晶胞质量，υ为泊松比， 为约化普朗克常数，f(υ)为泊松比函数Bs为绝热体弹模量，可使用步骤1.2中获得的体积-能量关系数据获得，其计算方法如下：
[0035]
[0036] 对于某一给定的(p，T)，非平衡Gibbs自由能对体积求极小值，即：
[0037]
[0038] 由式(14)可获得体系热力学状态，即获得体系的普适状态方程，明确体系的p、T、V三者的关系。进一步的，某一温度的体弹模量可表示为：
[0039]
[0040] 通过该步计算可将温度效应引入体弹模量，计算特定温度下的体弹模量BT。
[0041] 通过公式(16)
[0042] ET＝3BT(1-υ) (16)
[0043] 可以获得不同温度下的杨氏模量。利用公式(16)计算得到的0K下的杨氏模量ET，0与步骤1.1计算得到的0K下的杨氏模量E存在差值ΔE＝ET，0-E。此差值作为利用公式(16)计算得到的其他任意温度的杨氏模量的零点修正项，以此计算任意温度的杨氏模量E(T)＝ET-ΔE；
[0044] 进而利用公式 计算特定温度下的剪切模量G(T)；
[0045] 2.2对不同温度下的多相复合材料中的金属相和陶瓷相的热膨胀系数分别进行计算
[0046] 从另一方面，恒压热容Cv和Gruneisen参数γ可分别描述为：
[0047]
[0048]
[0049] 基于式(15)、(17)、(18)，可获得热膨胀系数α(T)：
[0050]
[0051] 由此，可将温度效应引入热膨胀系数，计算不同温度下的热膨胀系数；
[0052] (3)微观尺度的分子动力学模拟
[0053] 微观尺度的分子动力学模拟可以获得多相复合材料中的金属相在特定的温度(即可以取不同的温度)、特定晶粒尺寸下的塑性性质；硬脆陶瓷相的塑性相比金属相极其微小，可忽略不计。
[0054] 首先通过LAMMPS软件获得多相复合材料中金属相的单晶模型，随后通过AtomEye软件中的Vironoi方法构造相对窄范围内的不同晶粒尺寸的多晶模型(晶粒尺寸在30到40nm之间)；对不同相对窄范围内的晶粒尺寸多晶模型在特定温度下进行压缩模拟，绘制应力应变曲线，从应力应变曲线读取屈服强度和硬化系数。假定硬化系数与晶粒尺寸无关，为了减少误差，针对所有晶粒尺寸的多晶模型的硬化系数采用多个不同小晶粒尺寸的硬化系数的平均值；
[0055] 对于任意不同晶粒尺寸模型的屈服强度σs，利用上述从应力应变曲线读取的相对窄范围内的多晶模型的屈服强度数据和Hall-Petch关系进行拟合；
[0056] σs＝σ0+kd-1/2 (20)
[0057] 其中σs为屈服强度，σ0为单个位错移动时晶格摩擦阻力，k为常数，d为可取的任意晶粒尺寸；
[0058] 拟合过程中，σ0采用派纳力τp的计算结果，σ0≈τp；派纳力的取值为金属相的主滑移系的派纳力数值，计算公式如下：
[0059]
[0060] 其中υ为泊松比，取步骤1.1中的计算结果；a为滑移面的面间距，b为滑移方向上的原子间距，可以利用步骤1.1中弛豫的晶体结构通过几何关系计算得到；G(T)为特定温度下的剪切模量，取步骤2.1中的计算结果；
[0061] (4)构建多相复合材料的真实微观组织的三维有限元模拟几何模型
[0062] 通过聚焦离子束双束系统，利用离子束对多相复合材料加工获得边长为微米级的长方体，如边长约为15μm的立方体，以20nm的间隔对实验区域进行断层扫描；通过聚焦离子束实验获得图像堆栈后，采用Avizo软件进行图像处理，图像处理过程如下：
[0063] 首先应用FIB Stack Wizard模块对图像进行对齐、剪切、变形修正和阴影修正，进而采用中值滤波器和保持边界平滑的高斯滤波器进行去噪处理，最后，采用阈值分割方式对图片进行二值化处理，区分图像中金属相和陶瓷相，获得可用于构建有限元几何模型的复合材料组织图片；然后构建有限元几何模型；
[0064] 构建复合材料组织的有限元几何模型过程如下：(a)通过Matlab提取上述所述的复合材料组织图片中金属相所有像素的坐标即金属相像素坐标；(b)在有限元软件中建立复合材料组织图片的有限元几何模型即总有限元几何模型，即将金属相和陶瓷相都等同是陶瓷相，建立有限元几何模型；并获取总有限元几何模型中所有单元的坐标；(c)然后结合(b)中总有限元几何模型中所有单元的坐标，根据步骤(a)中的金属相像素坐标确定总有限元几何模型所有单元中应属于金属相的单元，并将该类单元的材料属性修改为金属相材料属性，完成复合材料组织的有限元几何模型建立；
[0065] (5)利用有限元模拟方法对多相复合材料力学行为进行模拟
[0066] 步骤(1)-(3)分别针对的是多相复合材料中独立或单独的金属相和陶瓷相的性能进行计算。根据步骤(1)中计算得到的泊松比、步骤(2)中计算得到的特定温度下的金属相和陶瓷相的杨氏模量与热膨胀系数，步骤(3)中计算得到的特定温度下特定晶粒尺寸的金属相的屈服强度和硬化系数，构建成陶瓷相和金属相的不同温度有限元参数模型；
[0067] 读取步骤(4)中构建的复合材料组织的有限元几何模型，将陶瓷相和金属相的不同温度有限元参数模型赋予到复合材料组织的有限元几何模型中，定义边界条件和约束条件；然后对复合材料组织的有限元几何模型进行复杂应力的加载和求解，进一步量化分析复合材料在不同温度、不同复杂应力条件下的应力应变关系、应力分布及其演变规律、塑性变形等力学行为。
[0068] 本发明的特色在于：
[0069] (1)本发明建立的多尺度模型和计算方法，以统一的研究对象和变量组合将不同尺度的模型和方法紧密连接、耦合，进行多重跨尺度数据流的传输、贯通，首次实现了纳观-微观-介观-宏观的连续性计算和模拟。
[0070] (2)本发明建立的多尺度模型和计算方法，可同时实现对多元体系材料成分的设计筛选、组织结构优化和服役性能预测分析。
[0071] (3)多尺度计算模型的建立，引入多相复合材料的真实微观组织和温度、压力等实际参量，可充分阐释多相复合材料在服役过程中的力学行为规律。
[0072] (4)首次实现了多尺度连续计算和服役过程模拟，解决了传统有限元分析方法对多相复合材料因严重缺乏实验数据而导致无法准确模拟其力学行为的共性难题。
[0073] (5)本发明建立的多尺度模型和计算方法，可实现并发式高通量计算，并可以推广应用于多种多元多相复杂材料体系的计算模拟，为新型高性能材料的设计研发开拓了新的准确、高效的途径。
附图说明
[0074] 图1为多尺度模拟方法和模型的计算流程
[0075] 图2为实施例用于构建有限元几何模型的二值化组织图片示例
[0076] 图3为实施例重现的真实材料三维微观组织
[0077] 图4为实施例计算的体积-能量关系
[0078] 图5为实施例模拟结果验证曲线
[0079] 图6为实施例应力场分布模拟结果。
具体实施方式
[0080] 下面结合实施例对本发明作进一步说明，但本发明并不限于以下实施例。以下实施例应用本发明的多尺度模拟计算方法研究复合材料WC-Co在热应力与压缩外力的交互作用下，应力应变规律和应力及其分布的演变规律。
[0081] 实施例1：
[0082] 通过聚焦离子束实验获得WC-Co图像堆栈，对图片进行二值化处理，获得可用于构建有限元几何模型的组织图片，如图2所示。采用Avizo软件进行图像处理，获得真实微观组织的WC-Co三维有限元模拟几何模型，建立的WC-8Co有限元模型尺寸为1.77μm(轴向)×2.35μm(横向)×0.40μm(厚向)，如图3所示。
[0083] 通过第一性原理计算，获得WC-Co复合材料中的金属相Co和陶瓷相WC在0K下的杨氏模量、体弹模量、剪切模量、泊松比以及体积-能量关系。体积-能量关系如图4所示。
[0084] E(GPa) B(GPa) G(GPa) υ
WC 679.6 378.0 283.1 0.200
Co 276.1 205.6 108.2 0.276
[0085] 利用第一性原理计算得到的0K下的杨氏模量、体弹模量、剪切模量、泊松比以及体积-能量关系数据，通过热力学计算获得金属相Co和陶瓷相WC在300K、700K、1100K下的经过零点修正的杨氏模量和热膨胀系数。
[0086]
[0087] 利用分子动力学模拟，对晶粒尺寸分别为34nm、36nm、38nm的金属相Co的三个多晶模型，在300K、700K和1100K下进行压缩模拟计算，获得应力应变曲线，从应力应变曲线读取屈服强度和硬化系数的数值。
[0088] 并对金属相Co在300K、700K和1100K下的派纳力进行计算。
[0089]
[0090] 温度(K) 派纳力(GPa)300 0.237
700 0.200
1100 0.151
[0091] 进而利用从应力应变曲线读取的屈服强度数值和派纳力的计算结果，应用Hall-Petch关系进行拟合，获得晶粒尺寸为400nm的金属相Co在300K、700K、1100K下的屈服强度。假定金属相Co的硬化系数与晶粒尺寸无关，400nm的金属相Co的硬化系数取34nm、36nm、
38nm计算的硬化系数的平均值。
[0092] 温度(K) 屈服强度(GPa) 硬化系数(GPa)300 1287 168.7
700 1204 131.9
1100 927 139.2
[0093] 将前述计算所得到的金属相Co的泊松比，以及在300K、700K、1100K下杨氏模量、热膨胀系数、屈服强度、硬化系数，和陶瓷相WC的泊松比，以及在300K、700K、1100K下杨氏模量、热膨胀系数输入有限元参数模型，并读取利用前述通过聚焦离子束实验和图像处理获得的真实组织构建的有限元几何模型。采用自由边界条件，在三个相邻的面上设置对称边界条件，在另外三个对立的面上设置自由边界条件。设置模型从1100K冷却至300K模拟烧结冷却过程。初始温度1100K时设为无应力状态，构建随时间变化的弹塑性模型。烧结态复相组织单轴压缩模拟中，加载过程为从0MPa至-3000MPa，每一载荷步增加100MPa，卸载过程为从-3000MPa至0MPa，每一载荷步降低100MPa。
[0094] 对压缩过程中Co相内应力应变进行统计，并与文献中中子衍射法测试结果进行对比，如图5所示。本模型模拟结果与实验测试结果中应力应变曲线具有相似的趋势，验证了本模型的合理性。在此基础上，可实现对压缩过程中应力应变规律、应力演变规律、塑性变形行为进行基于真实微观组织的量化表征与分析，如图6所示为压缩前、加载后和卸载后的应力场分布。继而研究微观组织-塑性变形行为-宏观性能之间的联系，为复相材料微观组织的设计与调控提供可靠指导。
