[发明专利]一种多相复合材料力学行为的多尺度模拟方法

权利要求书

1.一种多相复合材料力学行为的多尺度模拟方法，其特征在于，包含以下步骤：

(1)纳观尺度的第一性原理计算

纳观尺度的第一性原理计算可以分别获得多相复合材料中单独的金属相和陶瓷相在0K下的弹性性质以及不同体积的晶体结构对应的能量，计算步骤如下：

1.1对0K下的多相复合材料中的金属相和陶瓷相的杨氏模量、体弹模量、剪切模量、泊松比进行计算：

首先分别对金属相和陶瓷相材料的单晶结构进行结构弛豫；利用弛豫后的晶体结构对材料的弹性矩阵C的矩阵元，即弹性常数，进行计算；

其中σ_iε_ic_ij(i,j＝1,2,3,4,5,6)分别为应力、应变和单晶材料的弹性常数；

进而利用Voigt-Reuss-Hill近似，将计算得到的弹性常数转化为杨氏模量、体弹模量、剪切模量和泊松比；

B＝(B_v+B_r)/2 (2)

G＝(G_v+G_r)/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近似计算方法；

B_v＝(c₁₁+c₂₂+c₃₃)/9+2(c₁₂+c₂₃+c₁₃)/9 (6)

G_v＝(c₁₁+c₂₂+c₃₃-c₁₂-c₁₃-c₂₃)/15+(c₄₄+c₅₅+c₆₆)/5 (7)

B_r＝1/(s₁₁+s₂₂+s₃₃+2s₁₂+2s₁₃+2s₂₃) (8)

G_r＝15/(4s₁₁+4s₂₂+4s₃₃-4s₁₂-4s₁₃-4s₂₃+3s₄₄+3s₅₅+3s₆₆) (9)

其中s_ij(i,j＝1,2,3,4,5,6)为材料的柔度矩阵S的矩阵元，即柔度系数；柔度矩阵S与弹性矩阵C互为逆矩阵；

1.2对0K下的多相复合材料中的金属相和陶瓷相的晶体体积-能量关系分别进行计算：

将步骤1.1中弛豫后的金属相和陶瓷相的晶体结构的晶胞体积分别进行-10％至+10％范围内的收缩或者膨胀，计算对应不同晶胞体积的能量；

(2)介观尺度的热力学计算

介观尺度的热力学计算，基于准谐德拜模型，可以分别获得多相复合材料中的金属相和陶瓷相在不同温度下的弹性性质和热膨胀系数，计算步骤如下：

2.1对不同温度下的多相复合材料中的金属相和陶瓷相的杨氏模量分别进行计算

不同温度下的杨氏模量可以通过不同温度下的体弹模量计算得到，而计算不同温度下的体弹模量需要首先计算非平衡态Gibbs自由能；

非平衡态Gibbs自由能可描述为：

G^*＝E(V)+pV+A_vib(θ,T) (10)

其中E(V)表示步骤1.2中计算得到的晶胞总能量，p和V分别为外部压强和晶胞体积，T为温度，A_vib为振动Helmholtz自由能，θ为德拜温度；振动Helmholtz自由能A_vib可表示为：

其中D(θ/T)为德拜积分，n为每个晶胞的原子数，k_B为玻尔兹曼常数，德拜温度可表示为：

其中M为晶胞质量，υ为泊松比，为约化普朗克常数，f(υ)为泊松比函数B_s为绝热体弹模量，可使用步骤1.2中获得的体积-能量关系数据获得，其计算方法如下：

对于某一给定的(p,T)，非平衡Gibbs自由能对体积求极小值，即：

由式(14)可获得体系热力学状态，即获得体系的普适状态方程，明确体系的p、T、V三者的关系；进一步的，某一温度的体弹模量可表示为：

通过该步计算可将温度效应引入体弹模量，计算特定温度下的体弹模量B_T；

通过公式(16)

E_T＝3B_T(1-υ) (16)

可以获得不同温度下的杨氏模量；利用公式(16)计算得到的0K下的杨氏模量E_T,0与步骤1.1计算得到的0K下的杨氏模量E存在差值ΔE＝E_T,0-E；此差值作为利用公式(16)计算得到的其他任意温度的杨氏模量的零点修正项，以此计算任意温度的杨氏模量E(T)＝E_T-ΔE；

进而利用公式计算特定温度下的剪切模量G(T)；

2.2对不同温度下的多相复合材料中的金属相和陶瓷相的热膨胀系数分别进行计算

从另一方面，恒压热容C_v和Gruneisen参数γ可分别描述为：

基于式(15)、(17)、(18)，可获得热膨胀系数α(T)：

由此，可将温度效应引入热膨胀系数，计算不同温度下的热膨胀系数；

(3)微观尺度的分子动力学模拟

微观尺度的分子动力学模拟可以获得多相复合材料中的金属相在特定的温度、特定晶粒尺寸下的塑性性质；硬脆陶瓷相的塑性相比金属相极其微小，可忽略不计；

首先通过LAMMPS软件获得多相复合材料中金属相的单晶模型，随后通过AtomEye软件中的Vironoi方法构造相对窄范围内的不同晶粒尺寸的多晶模型，相对窄范围内的晶粒尺寸为30到40nm之间；对此相对窄范围内的晶粒尺寸多晶模型在特定温度下进行压缩模拟，绘制应力应变曲线，从应力应变曲线读取屈服强度和硬化系数；假定硬化系数与晶粒尺寸无关，为了减少误差，针对所有晶粒尺寸的多晶模型的硬化系数采用多个不同小晶粒尺寸的硬化系数的平均值；

对于任意不同晶粒尺寸模型的屈服强度σ_s，利用上述从应力应变曲线读取的相对窄范围内的多晶模型的屈服强度数据和Hall-Petch关系进行拟合；

σ_s＝σ₀+kd^-1/2 (20)

其中σ_s为屈服强度，σ₀为单个位错移动时晶格摩擦阻力，k为常数，d为可取的任意晶粒尺寸；

拟合过程中，σ₀采用派纳力τ_p的计算结果，σ₀≈τ_p；派纳力的取值为金属相的主滑移系的派纳力数值，计算公式如下：

其中υ为泊松比，取步骤1.1中的计算结果；a为滑移面的面间距，b为滑移方向上的原子间距，可以利用步骤1.1中弛豫的晶体结构通过几何关系计算得到；G(T)为特定温度下的剪切模量，取步骤2.1中的计算结果；

(4)构建多相复合材料的真实微观组织的三维有限元模拟几何模型

通过聚焦离子束双束系统，利用离子束对多相复合材料加工获得边长为微米级的长方体，以20nm的间隔对实验区域进行断层扫描；通过聚焦离子束实验获得图像堆栈后，采用Avizo软件进行图像处理，图像处理过程如下：

首先应用FIB Stack Wizard模块对图像进行对齐、剪切、变形修正和阴影修正，进而采用中值滤波器和保持边界平滑的高斯滤波器进行去噪处理，最后，采用阈值分割方式对图片进行二值化处理，区分图像中金属相和陶瓷相，获得可用于构建有限元几何模型的复合材料组织图片；然后构建有限元几何模型；

(5)利用有限元模拟方法对多相复合材料力学行为进行模拟

根据步骤(1)中计算得到的泊松比、步骤(2)中计算得到的特定温度下的金属相和陶瓷相的杨氏模量与热膨胀系数，步骤(3)中计算得到的特定温度下特定晶粒尺寸的金属相的屈服强度和硬化系数，构建成陶瓷相和金属相的特定温度有限元参数模型；

读取步骤(4)中构建的复合材料组织的有限元几何模型，定义边界条件和约束条件，将陶瓷相和金属相的特定温度有限元参数模型赋予到复合材料组织的有限元几何模型中进行求解。