王学智 汤雨婷 车军伟 令狐佳珺 侯兆阳
1) (长安大学理学院应用物理系,西安 710064)
2) (湖南大学机械与运载工程学院,汽车车身先进设计制造国家重点实验室,长沙 410082)
3) (西安交通大学物理学院,教育部物质非平衡合成与调控重点实验室,西安 710049)
具有非晶状热导率的固体材料在热能转换和热管理应用中备受青睐.因此,揭示晶体材料的非晶状热传导机理对于开发和设计低热导率材料至关重要.本文运用原子模拟方法揭示了萤石结构二元简单晶体Yb3TaO7 的非晶状低热导率的物理机理.研究发现,萤石Yb3TaO7 的低热导率主要是由O-Yb 和O-Ta 之间的原子间结合力相差较大引起的.这种相差较大的原子键可以极大地软化声子模式,从而抑制声子输运.振动模式分解显示,萤石Yb3TaO7 中的大多数声子模式位于Ioffe-Regel 极限以下,表现出强烈的扩散特征.萤石Yb3TaO7 中绝大部分(＞ 90%)的热流是通过扩散模式而不是传播模式传输.因此,萤石Yb3TaO7 中的热传导表现出独特的类非晶特性.同时发现,萤石Yb3TaO7 中的光学声子模式在热传导中发挥着重要的作用.本文对于原子间结合力与低热导率之间关系的认识,以及开发和设计低热导率材料提供了新思路.
热障涂层(TBC)广泛应用于航空发动机和燃气轮机的高温热端部件表面,以防止部件被高温腐蚀并提高热效率和服役寿命[1].随着航空发动机和燃气轮机向更高温度和高推重比发展,其内部温度也越来越高.因此,开发和设计具有超低热导率的TBC 材料是提高航空发动机和燃气轮机高温服役性能的有效途径.由于TBC 材料为绝缘体,热能在TBC 材料中主要以晶格振动(声子)为热载体进行传导[2].目前,基于声子的缺陷工程和高熵工程是降低材料固有热导率的两种主要方式.通过缺陷工程在基体材料中引入缺陷可显著阻碍声子输运,进而降低材料的热导率.例如,Y2O3稳定的ZrO2(YSZ) 中的大量氧空位有效地降低了ZrO2的热导率,甚至当Y2O3的摩尔质量掺杂比为12%时出现了非晶热传导现象[3].利用高熵工程在基体中引入外来原子增加材料的结构复杂性亦可使声子发生“瑞利散射”,进而获得低热导率.Zhu 等[4]报道了高熵TBC 材料(La0.2Nd0.2Y0.2Er0.2Yb0.2)2Zr2O7和(La0.2Nd0.2Sm0.2Gd0.2Yb0.2)2Zr2O7具有类似非晶的超低热导率.虽然缺陷和高熵工程可有效地减少热导率,但它们也带来一些科学和工程上的挑战.例如,当Y2O3的摩尔质量掺杂比超过12%时,YSZ 的热导率将达到极限[3].同时,过多的氧空位又会降低YSZ 的力学稳定性和断裂韧性[5].同样,高熵工程法会提高成本以及技术复杂性.例如,利用物理气相沉积(PVD)或等离子体喷雾(PS)制备高熵涂层的过程中,由于各元素的蒸发度不同,所获得试样的成分通常会偏离所需的成分,甚至有可能无法制备出所需的试样[6].因此,在简单体系的材料中实现低热导率就显得非常重要.
相对于晶体材料,完全无序或非晶材料具有非常低的热导率[7].稀土钽酸盐Yb3TaO7通常具有两种晶体结构,即无序的萤石结构(fluorite,F)和有序的冰镁石结构(weberite,W)[8].近年来,萤石结构的Yb3TaO7成为一种重要的TBC 材料,并且由于其低热导率和优异的相位稳定性而引起了研究人员极大的兴趣[9].有趣的是,尽管这种晶体材料成分简单,却拥有反常的非晶状热传导特性.然而,目前对这种低热导率的微观机理仍不清晰.特别是,对其表现出的非晶状热传导的机理还缺乏深入理解.揭示这些机制有利于低热导率材料的开发和设计.另一方面,实验测得的热导率通常因试样密度和热辐射的存在而变化[10,11].因此,仅依靠实验手段研究萤石Yb3TaO7在更广泛温度范围内的固有热输运性质是具有挑战性的.为此,采用原子模拟来研究Yb3TaO7的固有热导率及其机理.
在本文中,首先利用分子动力学(MD)计算了Yb3TaO7在300—1500 K 温度范围内的热导率.接着,利用声子谱能量密度和晶格动力学揭示了Yb3TaO7中的声子模式.最后,利用双通道声子模型定量计算了Yb3TaO7中各声子模式对热导率的贡献.研究发现,波动的原子间相互作用所驱动的扩散声子模式导致简单成分的萤石Yb3TaO7表现出类似非晶的低热导率.该研究表明,可以通过增大化学键之间的差异性来软化声子模式,进而可在简单材料中实现低热导率.
图1 所展示的是萤石型和冰镁石型Yb3TaO7的晶体结构.对于萤石型Yb3TaO7,其晶体结构的空间群为Fm3m,Yb 和Ta 原子随机占据4a(0,0,0)位点,氧原子和1/8 的氧空位随机占据8c(1/4,1/4,1/4)位点.对于冰镁石型Yb3TaO7,其晶体结构的空间群为C2221,Yb 原子占据4b(0,1/2,0) 和8g(x,y,1/2)位点,Ta 原子占据4a(0,0,0)位点,O 原子占据4c(1/2,y,1/4),16h(x,y,z)和8g(x,y,1/4)位点.本文中的Yb3TaO7超胞模型包含8×8×8 个原胞,共有6144 个原子.
图1 Yb3TaO7 (a) 萤石结构;(b)冰镁石结构Fig.1.Yb3TaO7: (a) Fluorite-type;(b) weberite-type.
Yb3TaO7的原子间相互作用以基于刚体的波恩模型描述,其表达式为
其中第一项是长程库仑势能,后面两项为短程排斥势能.ε0为真空介电常数,rij表示第i个原子和第j个原子之间的距离,qi和qj表示相应的电荷,Aij,ρij和Cij是可调参数,如表1 所列.同时,采用精度为1.0×10—6的Ewald 近似算法来计算长程库仑势.所有的计算都是通过LAMMPS 软件来实现,其中在三维空间的每个方向都是周期性边界,并利用Nose-Hoover 热浴法来控制体系的温度和压力.
表1 Yb3TaO7 的力场参数[12,13]Table 1.Force field parameters for Yb3TaO7[12,13].
Yb3TaO7的热导率(κ)由基于波动耗散理论的平衡分子动力学(equilibrium molecular dynamics,EMD)计算而来,具体表达式如Green-Kubo (G-K)公式所示[14]:
其中kB表示玻尔兹曼常数;T表示温度;V表示模拟体系的体积,〈 J(0)·J(t)〉 表示热流自相关函数(heat current autocorrelation function,HCACF).EMD-GK 方法可消除模拟体系的尺寸效应,即得到材料的固有热导率.热电流J计算表达式为
其中υi,ri,mi,vi和Ui分别表示第i个原子的速度、位置、质量、位力和势能.本文对体系进行了10 次独立计算以计算出平均热导率和误差.
本文采用声子谱能量密度(spectral energy density,SED)法计算Yb3TaO7的声子色散关系.该方法直接将原子的振动轨迹投影到振动模式上,以确定振动能量在波矢-频率空间的分布[15-18]:
晶格振动模式(声子)的分析基于经典晶格动力学(lattice dynamics,LD)展开[19]:
其中uiλα表示质量为mi的原子在α方向的振动位移;eλ,iα表示声子本征矢量;ωλ表示模式λ的声子本征值(频率);t表示时间.为了获得本征值ωλ,需求解波动方程:
其中D表示动力学矩阵.D可根据下式计算而来:
式中,Ep表示总势能.
3.1 热导率
为了对比研究,本文以冰镁石结构的Yb3TaO7为参照物.在计算热导率前,首先需要对模型和力场的可靠性进行验证.为此,将晶格常数、格林艾森常数和弹性模量的计算值与文献报道的实验值[8]进行了对比,如表2 所列.从表2 可以看出,晶格常数的计算值与实验值非常吻合,故本文的力场参数可以正确描述Yb3TaO7的晶体结构.格林艾森常数可以定量地描述固体材料中声子间的非简谐相互作用强度[20],故可用来验证力场描述声子的可靠性.此外,格林艾森常数的计算值与实验值也相当吻合.最后,对比了弹性模量,因为弹性模量也与声子热输运有关.从表2 还可以看出,弹性模量的计算值与实验值吻合较好.总之,通过对比Yb3TaO7的结构参数和热力学性质的计算值与实验值,证实了我们所采用力场的可靠性.
表2 计算的La2Zr2O7 和Yb3TaO7 的晶格常数、格林艾森常数和弹性模量Table 2.Calculated lattice constants,Grüneisen constants,and elastic modulus for La2Zr2O7 and Yb3TaO7.
基于验证的原子模型和力场,利用EMD-GK法计算了热导率.为了获得Yb3TaO7的固有热导率,需要对模拟时间进行验证.图2(a)和图2(b)显示的是萤石Yb3TaO7在1500 K 下的HCACF 和热导率随关联时间的变化关系.从图2(a)和图2(b)中可以看出,HCACF 随着关联时间的增加而迅速衰减,同时热导率趋于收敛,这与波动-耗散理论是一致的.因此,本文中使用的20 ps 的关联时间足以获得收敛的热导率.接着,研究了超胞的大小对计算热导率的影响,如图2(c)所示.从图2(c)可以看出,热导率随着超胞的大小的增加依然保持收敛,这是因为本文中选取的超胞尺寸远大于声子平均自由程(phonon mean free path,PMFP)(0.41—0.33 nm).该结果同时与EMD-GK 方法可以消除尺寸效应的原则是一致的.总之,上述结果说明本研究中选择的模拟时间和超胞尺寸足以获得固有的热导率.
图2 1500 K 温度下的热导率计算结果 (a) 归一化热流关联函数随关联时间的变化关系;(b) 热导率随关联时间的变化关系;(c) 不同超胞的热导率;(d) 热导率随温度的变化关系Fig.2.Calculated thermal conductivity at 1500 K: (a) Normalized HCACF versus correlation time;(2) thermal conductivity versus correlation time;(c) calculated thermal conductivity with different supercell;(d) temperature dependence of thermal conductivity.
图2(d)显示的是萤石型和冰镁石型Yb3TaO7在300—1500 K 的温度范围内热导率的计算结果.同时,将热导率的计算结果与文献报道的实验结果进行了对比.从图2(d)可以看出,萤石型Yb3TaO7的热导率的计算结果与实验结果非常吻合,但目前还没有冰镁石型Yb3TaO7的热导率的实验报道.计算结果表明,萤石型Yb3TaO7的热导率随温度的增大而增大,表现出反常的非晶热传导特征.相反,冰镁石型Yb3TaO7的热导率随温度的增大而减小,表现出经典的晶体热传导特征.热导率对温度的依赖关系可定量地描述为
其中A,B和n是材料参数.根据(8)式,萤石结构和冰镁石结构的n值分别为—0.08 和1.032.更重要的是,在低温下,萤石结构的热导率远低于冰镁石结构的热导率.上述结果表明,萤石型二元氧化物Yb3TaO7具有类非晶的低热导率.因此,揭示这一现象背后的机理对于开发和设计低热导率的材料具有重要的意义.
3.2 声子散射
尽管EMD-GK 计算给出了和实验值相一致的热导率,但它并没给出关于声子层面的机理描述.声子气模型(PGM)是理解微纳尺度热输运的经典模型[21].假设晶体萤石Yb3TaO7中所有声子模式都是正常模式,即声子有明确的波矢和周期性的本征矢量(速度场).基于PGM,萤石Yb3TaO7的微观热传导可描述为
其中V是晶胞体积;C(ωλ)是模式比热容;υg(ωλ)是声子群速度;τ(ωλ)是声子寿命(弛豫时间).C(ωλ)由玻色-爱因斯坦统计学计算:
其中 ℏ 是约化普朗克常数.基于声子色散关系,群速度υg可根据公式υg=∂ω/∂ q计算而来.利用洛伦兹函数拟合SED 可计算出声子寿命τ(ωλ)[15]:
其中I表示SED 的峰值;ωc表示峰中心对应的频率;γ表示SED 峰的半高宽(HWHM).据此,声子寿命可表示为τλ=1/2γ.以300 K 为例,图3 显示的是萤石Yb3TaO7的模式比热容、声子群速度和声子寿命的计算结果.基于PGM,Yb3TaO7的热导率为0.36 W·m—1·K—1,该值远小于EMD-GK 值(1.25 W·m—1·K—1)以及实验值(1.10 W·m—1·K—1).如此大的差异表明,基于PGM 的单通道模型无法准确描述萤石Yb3TaO7的声子热传导性质.事实上,这并不是PGM 第一次在描述声子热传导时失效.例如,PGM 无法准确处理超晶格[22]和声子晶体中的相干声子[23],也无法处理La2Zr2O7[24]和α-DX[25]中的扩散声子.上述结果意味着萤石Yb3TaO7中包含不明确的声子模式,而这些模式无法用PGM描述.
图3 在300 K 温度下,萤石Yb3TaO7 的热输运性质 (a) 模式比热容;(b) 声子寿命;(c)声子群速度Fig.3.Thermal transport properties of F-Yb3TaO7 at 300 K: (a) Mode capacity;(b) phonon lifetime;(c) phonon group velocity.
基于此,利用SED 考察了声子色散关系,如图4所示.从图4 可以看出,萤石Yb3TaO7的色散关系只有在1.5 THz 以下的声学支部分是清晰的,而高频光学模式有很大的声子线宽且色散曲线极其模糊.相比于萤石Yb3TaO7,冰镁石Yb3TaO7的声子色散关系就显得相对清晰明确.接着,将SED投影到频率空间,每个峰对应一支色散关系,其峰宽和声子散射率有关.从图4 可以看到,萤石Yb3TaO7的SED 峰在低频时较尖锐,而在高频时则明显变宽,说明萤石Yb3TaO7中存在很强的声子散射.然而,冰镁石Yb3TaO7的SED 峰相比于萤石Yb3TaO7更尖锐,该结果与其相对清晰的色散关系相一致.同时可以看到,萤石Yb3TaO7的SED 峰值远小于冰镁石Yb3TaO7的SED 峰值,说明萤石Yb3TaO7中的声子模式具有非常低的振动能量和传热率,进而表现为低热导率.综上,SED 分析表明,萤石Yb3TaO7由于强声子散射而包含非正常的声子模式.
图4 Yb3TaO7 的SED 计算结果 (a) 萤石;(b) 冰镁石Fig.4.Calculated SED plots of Yb3TaO7: (a) F-type;(b) W-type.
萤石型和冰镁石型的Yb3TaO7具有完全相同的元素成份和氧空位浓度,故离子半径差、离子质量差和氧空位无法解释萤石型Yb3TaO7的强声子散射.声子是晶格的集体振动,故强声子散射可能和原子间相互作用有关.目前,还没有可以定量地描述声子散射率和原子间相互作用的理论.因此,分别计算了萤石型和冰镁石型Yb3TaO7的原子间相互作用力,以给出定性解释.如图5(a)所示,萤石型Yb3TaO7中O-Ta 和O-Yb 之间的结合力差值远大于冰镁石结构中的差值,即萤石型Yb3TaO7的原子间相互作用是高度不均匀的.为了解释这种不均匀相互作用力对声子的散射,以一维双原子链为例.一维双原子链的声子色散关系为
图5 (a) 原子间相互作用力;(b) 一维双原子链的声子色散关系Fig.5.(a) Calculated interatomic bonding force;(b) phonon dispersion relationship of for a one-dimensional diatomic chain.
其中M和m为相邻原子的质量;G和g为相应的力常数(详细推导过程请见补充材料).如图5(b)所示,随着力常数G和g之间差异的增大,声子模式变得软化,声子群速度降低,表明声子散射增强.因此,萤石型Yb3TaO7的强声子散射是由其高度不均匀的原子间相互作用所导致.
3.3 声子模式
上述内容揭示了萤石Yb3TaO7中存在非正常的声子模式及其原因.为了直观地理解这些声子模式,将声子本征矢量归一化到单位极化球上[3]:
其中eiα,λ表示第i个原子在模式λ中α方向的本征矢量;e* 为共轭本征矢量.如果一个给定的声子模式λ是正常的声子模式,即它具有类似于平面波的周期性波矢和本征矢量,那么它的极化矢量将聚焦于极化球表面.这种正常的模式即传播模式,也叫作传播子(propagons).如果极化矢量随机地散布在极化球表面,那么这种模式是不明确的,对应于局域化模式和扩散模式,也叫作局域子(locons)和扩散子(diffusons).如图6 所示,当频率范围小于1.0 THz 时,萤石和冰镁石型Yb3TaO7的声子极化均聚焦于球面上点,这对应于具有明确色散关系的声学模式.当频率范围扩大到2.5 THz时,萤石结构的声子极化开始随机分布在球面上,声子本征矢量的极化表现出明显的局域化特点.然而,冰镁石结构的声子极化仍聚焦于球面,且可以分辨.萤石结构中随机分布的模式对应于非正常模式.最后,我们将频率范围扩大到25 THz.从图6可以看出,萤石结构中绝大多数的声子极化随机分布于球面而无法区分,而冰镁石结构的大量声子聚焦于极化球的赤道区域.这一显著差异是因为萤石结构的原子间相互作用是高度不均匀的,而冰镁石结构的原子间相互作用是相对均匀的.总之,声子本征矢量的极化分析直观地证实了Yb3TaO7同时具有局域化和非局域化的声子模式.
图6 Yb3TaO7 在不同频率区间内的声子极化 (a) 萤石;(b) 冰镁石Fig.6.Phonon polarization of Yb3TaO7 in different frequency domains: (a) F-type;(b) W-type.
接下来,对整个声子谱进行了模式分解,以精确确定局域子、传播子和扩散子的分布区间.局域子可直接通过声子参与率(participation ratio,PR)来确定[26]:
其中N是体系的总原子数.PR 定量地衡量了体系中有多少原子对给定的声子模式λ有贡献.一般地,如果PR ＜ 0.1,则该模式是局域子.如图7(a)和图7(b)所示,萤石Yb3TaO7在高频有大量的局域子,而冰镁石Yb3TaO7在低频和高频均存在少量的局域子.这些局域子是局域化的振动模式,对热传导的贡献可以忽略不计.同时可以看到,萤石结构的声子PR 小于冰镁石结构的声子PR,说明非均匀的原子间相互作用抑制了声子模式的激发.总的声子PR 可以进一步投影到各原子上,以确定不同的原子对声子模式的贡献:
如图7(c)所示,萤石Yb3TaO7低于7 THz 的低频声子主要来源于Yb 原子的振动,高于7 THz的高频声子主要来源于O 原子的振动,而Yb 的振动对声子的贡献相当小.此外,Ta 原子的PR 在7 THz 处达到最大值,而Yb 原子的PR 在4 THz处达到最大值.这一差异可归因于O-Ta 的原子间结合力远大于O-Yb 的原子间结合力.如图7(d)所示,冰镁石Yb3TaO7低于7 THz 的低频声子也主要来源于Yb 原子的振动,介于7—8 THz 的声子主要来源于Ta 原子的振动,高于8 THz 的声子则主要来源于O 原子的振动.与萤石Yb3TaO7不同的是,冰镁石Yb3TaO7中Ta 原子对声子的贡献明显较大.
图7 (a) 萤石Yb3TaO7 和(b) 冰镁石Yb3TaO7 的总参与率;(c) 萤石Yb3TaO7 和(d) 冰镁石Yb3TaO7 中各元素的参与率Fig.7.Total phonon participation ratio of (a) F-Yb3TaO7 and (b) W-Yb3TaO7;atomic phonon participation ratio of (c) F-Yb3TaO7 and (d) W-Yb3TaO7.
传播子和扩散子可利用Ioffe-Regel 准则区分[27].根据该准则,若给定声子的平均自由程(MFP)大于体系的原子间最小距离,则该声子是传播子,否则是扩散子.如图8 所示,萤石Yb3TaO7中 90%的声子模式落在Ioffe-Regel 区间,表现出扩散子特征.然而,冰镁石Yb3TaO7存在大量的传播子.此外,通过对比可发现,萤石Yb3TaO7中声子的MFP 整体上远小于冰镁石Yb3TaO7,这说明高度不均匀的原子间相互作用导致萤石Yb3TaO7的声子表现出强烈的扩散特征.
图8 Yb3TaO7 的声子平均自由程 (a) 萤石,(b)冰镁石;Yb3TaO7 中传播子和扩散子的热导率 (c) 萤石,(d) 冰镁石Fig.8.Phonon mean free paths of (a) F-Yb3TaO7 and (b) F-Yb3TaO7;thermal conductivity of propagons and diffusons in (c) FYb3TaO7 and (d) F-Yb3TaO7.
上述模式分解表明,热能在Yb3TaO7中是由传播子和扩散子传导的.为了定量地描述传播子和扩散子对热导率的贡献,我们对热导率进行了谱分解.如前所述,传播子具有周期性的波矢和本征矢量,其对热导率的贡献(κP)可根据(9)式计算.扩散子无周期性的波矢和本征矢量,且以随机扩散的方式进行热传导,故扩散子对热导率的贡献(κD)可根据Allen-Feldman (A-F)扩散理论计算[28]:
式中,D(ωλ)是声子模式λ的扩散率:
(15)式中的S λλ′是非对角化的热流矩阵:
其中Q为布里渊区内的波矢量之和;λ和λ′代表不同的声子模式.Θ Qλλ′可根据晶格动力学矩阵计算:
式中,e是声子本征矢量;α和β表示笛卡尔方向;Rs表示周期性超胞内每个单胞(标记为s)相对原胞(标记为0)的位移矢量;表示单胞内原子k相对原子k′的位移矢量,是根据(7)式计算的动力学矩阵.图8(a)和图8(b)展示的是萤石型和冰镁石型Yb3TaO7的模式热导率.从图8(a)和图8(b)中可看出,萤石Yb3TaO7的热导率主要由扩散子贡献,故萤石Yb3TaO7的热导率整体表现出非晶热传导特性.虽然当温度高于500 K时,冰镁石Yb3TaO7的热导率主要由扩散子贡献,但热导率随温度的变化关系仍由传播子主导,故冰镁石Yb3TaO7的热导率表现出晶体热传导特性.
综上,PR 和Ioffe-Regel 准则准确给出了萤石Yb3TaO7的局域子、传播子和扩散子的分布.因此,在Yb3TaO7中存在双通道声子热传导现象.热传导的模式分解表明,萤石Yb3TaO7的热导率主要是由高频的扩散子贡献的.为了预测或计算材料的热导率,目前已经建立了一些简化的物理模型,即通过一些简单的物理参数即可得到材料的热导率,且这些模型大都建立在高频光学声子对热导率可忽略的假设上[29-32].然而,本文中萤石Yb3TaO7的结果证实,高频的光学声子在热传导中起着不可忽视甚至是决定性的作用.这一结果为调控热传导提供了新视角.
本文从原子间结合和晶格振动模式的角度出发,研究了萤石Yb3TaO7晶体的非晶状热传导机理.萤石Yb3TaO7的低热导率主要是由O-Ta 和O-Yb 之间较大的差异驱动的软化声子导致的.由于存在大量的扩散振动模式,基于PGM 的单通道模型无法合理描述萤石Yb3TaO7的热传导.振动模式分解发现,萤石Yb3TaO7中的绝大部分热量是由扩散模式传输的,导致其热传导具有类似非晶的特性.同时发现,萤石Yb3TaO7中的光学声子模式在热传导中具有重要的作用,这可为调控热传导提供新的物理维度.该研究启发我们,增强化学键的不均匀性可以激发更多的扩散振动模式,从而有利于在简单材料中实现非常低的热导率.
猜你喜欢热传导声子萤石半无限板类声子晶体带隙仿真的PWE/NS-FEM方法汽车实用技术(2022年5期)2022-04-02一类三维逆时热传导问题的数值求解数学物理学报(2022年1期)2022-03-16冬天摸金属为什么比摸木头感觉凉？红蜻蜓·中年级(2021年12期)2021-12-19纳米表面声子 首次实现三维成像今日农业(2021年7期)2021-07-28声子晶体覆盖层吸声机理研究舰船科学技术(2021年12期)2021-03-29构造对河南南召县观音庙萤石矿的控制作用矿产勘查(2020年7期)2020-12-25具有非线性边界条件的瞬态热传导方程的二择一结果数学物理学报(2020年5期)2020-11-26基于声子晶体理论的导线防舞方法及数值验证三峡大学学报(自然科学版)(2017年1期)2017-03-20红山文化萤石质鱼形饰的鉴别和探讨文物保护与考古科学(2016年4期)2016-05-17某萤石矿开拓系统优化铜业工程(2015年4期)2015-12-29