光学自由曲面铣床重力/惯性力误差补偿

Gravity/Inertia Force Error Compensation for Optical Free-Form Surface Milling Machines

摘要

为了提高光学自由曲面的加工精度，针对传统机床在热和力方面的不足，构建了一种新的五轴机床拓扑结构，该拓扑结构可以有效减少热变形对机床精度的影响，并提出了重力/惯性力误差同步补偿方法.基于所建立的误差补偿模型对不同误差补偿方法的误差补偿效果进行了比较研究.结果表明，重力平衡结合力学分流立柱可以提高铣床精度，变力（variable force，VF）的误差补偿效果显著优于定力（constant force，CF），二维运动重力平衡（two-dimensional motion gravity balance，TMGB）的误差补偿效果显著优于固定重力平衡（fixed gravity balance，FGB）.最差的补偿方法为固定重力平衡结合定力，最好的补偿方法为二维运动重力平衡结合变力.固定重力平衡结合变力与二维运动重力平衡结合定力在不同条件下各有优缺点.

Abstract

In order to improve the machining accuracy of optical free-form surfaces， given the shortcomings of traditional machine tools in heat and force， a new five-axis machine tool topology is constructed， which can effectively reduce the influence of thermal deformation on the accuracy of machine tools， and a synchronous compensation method of gravity/inertial force error is proposed. The error compensation effect of different error compensation methods are compared and analyzed. The results show that the accuracy of milling machines can be greatly improved by using gravity balance combined with mechanical shunt column， and the error compensation effect of variable force （VF） is significantly better than that of constant force （CF）， the error compensation effect of two-dimensional motion gravity balance （TMGB） is significantly better than that of fixed gravity balance （FGB）. The worst compensation method is FGB combined with CF， while the best compensation method is TMGB combined with VF. The combination of FGB combined with VF and TMGB combined with CF have their own advantages and disadvantages under different conditions.

Graphical abstract

关键词

光学自由曲面 / 热误差 / 惯性力误差 / 重力误差 / 误差补偿 / 误差预防

Key words

optical free-form surface / thermal error / inertia force error / gravity error / error compensation / error prevention

光学系统中的自由曲面通常是指不能用球面或非球面系数表示的曲面，也指非传统的、不对称的或只能用参数向量表示的曲面，光学自由曲面是一种精密光学元件，需要用五轴联动精密机床进行加工.机床误差的控制方法主要有误差预防方法和误差补偿方法.机床加工装配误差、力载荷误差和热误差，极大地影响了机床的加工精度，因此，许多学者对机床的几何误差和热误差进行了研究，包括对机床的几何误差和热误差进行建模，对机床的关键误差进行识别，以及对机床误差进行补偿.为了消除重力引起的误差，Nguyen等^［1］设计了一种基于齿轮弹簧模块的三角形并联机器人重力补偿机构.Han等^［2］基于有限差分法，计算并实验验证了考虑材料缺陷的大跨度重型立式车床的重力变形曲线.Huang等^［3］提出了一种重力补偿主从控制方法.Miao等^［4］研究了温度敏感点变化对数控机床热误差补偿模型的影响.Andolfatto等^［5］估计了热扰动下五轴机床误差识别中自适应蒙特卡罗的不确定性.Qiao等^［6］基于指数乘积（POE）公式，建立公共垂线（CPL）变换的五轴机床几何误差校准模型.Huang等^［7］基于Abbe和Bryan原理对三轴机床的体积误差进行建模.Abdulshahed等^［8］利用灰色神经网络模型对龙门五轴机床的热误差进行了建模.Fan等^［9］对考虑重力变形的五轴机床进行了误差优化分配.Jiang等^［10］提出一种使用触发探针对五轴机床旋转轴上的12个误差进行单次夹紧识别的方法.Zhong等^［11］提出一种基于空载条件下双球杆仪测量五轴机床体积误差的评估标准.Xing等^［12］提出一种基于体积误差和矢量相似性测量的五轴机床精度状态监测方法.Li等^［13］提出一种用于机床热误差建模的重构变量回归方法.Xiang等^［14］对五轴机床的体积误差进行了建模和补偿.Blaser等^［15］研究了五轴机床热误差补偿的自适应学习控制方法.Li等^［16］构建一种带有误差补偿机构的五轴机床，并进行了动力学研究.综上所述，众多学者对机床的热误差和重力误差补偿进行了研究，但对惯性力误差补偿的研究还比较少.重力为固定值，重力误差为静态误差，会使机床产生静变形；惯性力为变化值，惯性力误差为动态误差，其除了会使机床产生静变形，还会使机床产生振动.

鉴于传统重力误差补偿机构补偿力的大小、方向和作用点都不可变，本文提出一种可以同时补偿重力和惯性力的误差补偿方法，并建立光学自由曲面铣床误差补偿模型（半解析模型），研究将惯性力的“动”问题转化为“静”问题的方法，并对不同误差补偿方法的误差补偿效果进行比较分析.铣床的上部重力和惯性力从Z轴外立柱传出，不影响直线轴的精度，铣床结构的对称设计可以减少热量对加工精度的影响.研究铣床的误差补偿，对于分析具有这种拓扑结构机床的极限加工精度具有重要意义.

1 光学自由曲面铣床

图1展示了光学自由曲面铣床主体结构图和结构拓扑图，其中Z轴电机安装在顶板上，Z轴由4个电机驱动，此时铣床为七轴五联动机床，2个运动轴为重力、惯性力误差平衡轴.当重力平衡未补偿时，铣床大部分重力和惯性力从内立柱传递，一小部分从外立柱传递；当重力平衡补偿时，铣床大部分重力和惯性力从外立柱传递，一小部分从内立柱传递.图2展示了光学自由曲面铣床的重力和惯性力传递路径，其中Z轴电机安装在立柱上.当重力平衡未补偿时，铣床的全部重力和惯性力都从内立柱传递；当重力平衡补偿时，铣床的全部重力和惯性力都从外立柱传递.图3a展示了固定重力平衡的工作原理，固定重力平衡只能平衡铣床悬空部分的重力，但不能消除重力矩对铣床精度的影响；图3b展示了二维运动重力平衡的工作原理，由于力臂长度为零，因此二维运动重力平衡既能平衡铣床悬空部分的重力，又能消除重力矩对铣床精度的影响.固定重力平衡补偿时由铣床的重力产生的重力矩（M₁）和惯性力产生的惯性力矩（M₂）分别为

M_{\mathrm{1}}=\sqrt[]{（G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}{）}^{\mathrm{2}}·x^{\mathrm{2}}+（G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}+G_Y{）}^{\mathrm{2}}·y^{\mathrm{2}}}；（1）

M_{\mathrm{2}}=\pm \frac{a}{g}\sqrt[]{（G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}{）}^{\mathrm{2}}·x^{\mathrm{2}}+（G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}+G_Y{）}^{\mathrm{2}}·y^{\mathrm{2}}}.（2）

其中：G_X为X轴运动机构重力；G_箱为主轴箱重力；G_轴为主轴部件重力；G_Y为Y轴运动机构重力；a为Z轴的加速度，向上取正值，向下取负值；g为重力加速度.因重力和惯性力共同作用产生的倾覆力矩为

M_{\mathrm{合}}=\left(\mathrm{1}\pm \frac{a}{g}\right)\times \sqrt[]{（G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}{）}^{\mathrm{2}}·x^{\mathrm{2}}+（G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}+G_Y{）}^{\mathrm{2}}·y^{\mathrm{2}}}.（3）

G为重力和惯性力的合力，铣床在加工热态时，立柱热源设置应保证4个立柱对溜板4个角的X轴和Y轴的热力分力分别指向远离铣床原点的方向，使溜板受到张力，以达到最佳的热设计效果.

2 误差补偿原理

重力平衡装置由二维主运动平台、二维随动平台和平衡缸组成.实现最佳误差补偿效果的重力平衡要求：当主轴轴线通过机床原点时，不仅要求溜板、X轴运动机构、Y轴运动机构、主轴箱和主轴的总重心投影为（0，0），同时溜板、Y轴运动机构、X轴运动机构、主轴箱和主轴的投影重心也均为（0，0）.重力平衡装置的平衡力为溜板重力G_{\mathrm{溜}}、X轴运动机构重力G_X、Y轴运动机构重力G_Y、主轴箱重力G_{\mathrm{箱}}和主轴部件重力G_{\mathrm{轴}}之和，二维主运动平台与铣床的X轴行程比{\gamma }_X为

{\gamma }_X=\frac{G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}}{G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}+G_{\mathrm{溜}}+G_Y} .（4）

二维主运动平台与铣床的Y轴行程比{\gamma }_Y为

{\gamma }_Y=\frac{G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}+G{}_Y{}^{}}{G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}+G_{\mathrm{溜}}+G_Y} .（5）

二维主运动平台原点与铣床原点的连线平行于Z轴.建立二维运动重力平衡的运动控制算法：当铣床在原点位置时，铣床悬空部分的重心坐标为（0，0），当铣床从原点位置移动到坐标（x，y）时，二维主运动平台也从原点位置移动到坐标

\left(\frac{G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}}{G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}+G_{\mathrm{溜}}+G_Y}x\right.，\left.\frac{G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}+G_Y}{G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}+G_{\mathrm{溜}}+G_Y}y\right)

位置.由二维运动重力平衡运动控制算法可知，二维主运动平台X轴坐标如式（6）所示：

X_{\mathrm{二维}}=\frac{G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}}{G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}+G_{\mathrm{溜}}+G_Y}{X}_{\mathrm{铣床}} .（6）

式（6）两端对时间求导，得到二维主运动平台X轴运动速度，如式（7）所示：

v_{\mathrm{二维}X}=\frac{G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}}{G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}+G_{\mathrm{溜}}+G_Y}{v}_{\mathrm{铣床}X} .（7）

式（7）两端再对时间求导，得到二维主运动平台X轴运动加速度，如式（8）所示：

a_{\mathrm{二维}X}^{}=\frac{G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}}{G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}+G_{\mathrm{溜}}+G_Y}{a}_{\mathrm{铣床}X} .（8）

式（6）~（8）中，X_{\mathrm{铣床}}，v_{\mathrm{铣床}X}，a_{\mathrm{铣床}X}分别为铣床的X轴坐标、X轴速度、X轴加速度.

由二维运动重力平衡运动控制算法可知，二维主运动平台Y轴坐标如式（9）所示：

Y_{\mathrm{二维}}=\frac{G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}+G_Y}{G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}+G_{\mathrm{溜}}+G_Y}{Y}_{\mathrm{铣床}}.（9）

式（9）两端对时间求导得到二维主运动平台Y轴运动速度，如式（10）所示：

v_{\mathrm{二维}Y}=\frac{G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}+G_Y}{G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}+G_{\mathrm{溜}}+G_Y}{v}_{\mathrm{铣床}Y} .（10）

式（10）两端再对时间求导，得到二维主运动平台Y轴运动加速度，如式（11）所示：

a_{\mathrm{二维}Y}=\frac{G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}+G_Y}{G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}+G_{\mathrm{溜}}+G_Y}{a}_{\mathrm{铣床}Y} .（11）

式（9）~（11）中，Y_{\mathrm{铣床}}，v_{\mathrm{铣床}Y}，a_{\mathrm{铣床}Y}分别为铣床的Y轴坐标、Y轴速度、Y轴加速度.

从式（6）至式（11）可以看出，当X轴和Y轴运动已知时，二维运动平衡的运动规律也是确定的.重力平衡的动态平衡力如式（12）所示：

F=（G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}+G_{\mathrm{溜}}+G_Y）·\left(\mathrm{1}\pm \frac{a}{g}\right) .（12）

动态平衡力可简化为

F=G_{\mathrm{1}}+G_{\mathrm{2}}.（13）

其中：G_{\mathrm{1}}为铣床Z轴重力，如式（14）所示，

G_{\mathrm{1}}=G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}+G_{\mathrm{溜}}+G_Y；（14）

G_{\mathrm{2}}为铣床Z轴惯性力，如式（15）所示，

G_{\mathrm{2}}=\pm \frac{a}{g}（G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}+G_{\mathrm{溜}}+G_Y）.（15）

3 误差补偿模型建立与结果

铣床的误差补偿模型建立可分为3个阶段.首先，建立铣床的简化三维几何模型；然后，进行材料特性定义、网格划分和接触单元生成（与网格划分质量无关），以生成等效的有限元模型；最后，根据铣床的工作状态，加载边界约束条件和载荷，生成等效力学模型.铣床的主体结构材料为铸铁.负载情况：切削载荷F_X =100 N，F_Y =100 N，F_Z =100 N.图4展示了光学自由曲面铣床的误差补偿模型，其中4个Z轴电机分别安装在4个立柱上.

图5a~5c分别展示了当重力平衡未补偿、固定重力平衡补偿和二维运动重力平衡补偿时，Z=0，X轴和Y轴处于不同工作位置时，在重力载荷下的铣床主轴端的空间误差.由于结构的对称性，铣床主轴端的空间误差沿X轴和Y轴对称分布.主轴端最大空间误差与最小空间误差的差值分别为3.5，3.4，2.2 μm，该差值将影响铣床非加工工况时的重复定位精度.可以看出固定重力平衡可以显著减小主轴端部最大空间误差和最小空间误差的差值，二维运动重力平衡可以进一步减小该差值.图6a~6c分别展示了当重力平衡未补偿、固定重力平衡补偿和二维运动重力平衡补偿时，在重力载荷和切削载荷下，Z=0，X轴和Y轴处于不同的工作位置时，铣床主轴端的空间误差.在原点位置时，主轴端部3个方向的刚度分别为K_X =92 N/μm，K_Y =53 N/μm和K_Z =32 N/μm.最大空间误差和最小空间误差之差分别为4.5，4.8，2.8 μm，该差值将影响铣床加工工况时的重复定位精度.

4 误差补偿方法对比研究

当铣床在原点位置时，图7为在Z轴不同加速度下使用不同补偿方法的Z轴导轨空间误差.

当铣床在最大行程时，图8展示了在Z轴不同加速度下使用不同补偿方法时Z轴导轨的直线度.可以看出，误差补偿可以显著提高Z轴导轨直线度.

当铣床Z轴加速度在1g~-1g范围内变化时，根据是否采用补偿方法，Z轴导轨直线度在0.045~3.100 μm之间变化.当铣床在原点位置时，图9展示了在不同Z轴加速度下使用不同补偿方法时Z轴丝杠的空间误差.当铣床在最大行程时，图10展示了在不同Z轴加速度下使用不同补偿方法时Z轴丝杠的直线度.可以看出，使用误差补偿时可以显著提高Z轴丝杠直线度.当铣床Z轴加速度在1g~-1g范围内变化时，根据是否采用补偿方法，Z轴丝杠直线度在0.3~4.8 μm之间变化.可以看出，二维运动重力平衡补偿法优于固定重力平衡补偿法，变力补偿法优于定力补偿法，固定重力平衡结合变力补偿与二维运动重力平衡结合定力补偿在不同条件下互有胜负，所有补偿方法均优于无补偿方法.

当铣床在原点位置时，图11展示了在Z轴不同加速度下使用不同补偿方法时Y轴导轨的空间误差.当铣床在最大行程时，图12展示了Z轴不同加速度下使用不同补偿方法时Y轴导轨直线度.可以看出，使用误差补偿时可以显著提高Y轴导轨直线度.当铣床Z轴加速度在1g~-1g范围内变化时，根据是否采用补偿方法，Y轴导轨直线度在0.5~29.0 μm之间变化.

当铣床在原点位置时，图13展示了在Z轴不同加速度下使用不同补偿方法时X轴导轨空间误差.当铣床在最大行程位置时，图14展示了在Z轴不同加速度下使用不同补偿方法时X轴导轨直线度.可以看出，使用误差补偿时可以显著提高X轴导轨直线度.当铣床Z轴加速度在1g~-1g范围内变化时，根据是否采用补偿方法，X轴导轨直线度在0.6~12.0 μm之间变化.

当采用固定重力平衡补偿法结合变力补偿法进行误差补偿时，铣床在X轴和Y轴最大行程时的重力矩为

M_{\mathrm{1}}=\sqrt[]{（G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}{）}^{\mathrm{2}}·x_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{2}}+（G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}+G_Y{）}^{\mathrm{2}}·y_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{2}}}，

惯性力矩为

M_{\mathrm{2}}=\pm \frac{a}{g}\sqrt[]{（G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}{）}^{\mathrm{2}}·x_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{2}}+（G_X+G_{\mathrm{箱}}+G_{\mathrm{轴}}+G_Y{）}^{\mathrm{2}}·y_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{2}}}，该倾覆力矩将降低铣床的稳定性.采用二维运动重力平衡补偿法结合变力补偿法进行误差补偿时，由重力和惯性力共同作用产生的倾覆力矩为零，此时铣床的稳定性较好.在重力载荷作用下，重力平衡未补偿的铣床空间误差由重力和惯性力引起，固定重力平衡补偿的铣床空间误差由重力矩和惯性力矩引起；在重力载荷及切削载荷共同作用下，重力平衡未补偿的铣床空间误差由重力、惯性力和切削力引起，固定重力平衡补偿的铣床空间误差由重力矩、惯性力矩和切削力引起.在重力载荷及切削载荷共同作用下，二维运动重力平衡补偿的铣床空间误差由切削力引起.固定重力平衡结合定力补偿法只能补偿重力误差，而无法补偿惯性力和重力矩误差；固定重力平衡结合变力补偿法能补偿重力和惯性力误差，但无法补偿重力矩和惯性力矩误差；二维运动重力平衡结合定力补偿法只能补偿重力和重力矩误差，而无法补偿惯性力误差；二维运动重力平衡结合变力补偿法对上述4种误差都可以补偿.

5 结 论

1）构建了一种完全中心对称的五轴机床拓扑结构——新型机床拓扑结构，受到几何约束，铣床各个方向的热变形相互抵消.同时，在机床的上部设置直线轴，可以减少直线轴发热产生的热量传入床身.与传统的机床拓扑结构相比，这种拓扑结构可以有效地减少热变形对机床精度的影响.

2） 提出了基于力学三要素的重力平衡结合力学分流立柱的重力和惯性力误差同步补偿方法，并进一步建立了重力平衡动态平衡力计算方程和二维运动重力平衡运动控制算法.与传统误差补偿方法相比，这种补偿方法可以提高误差补偿效果.

3） 误差补偿可以显著提高直线轴的导向精度和驱动精度，最差的补偿方法为固定重力平衡结合定力（传统误差补偿方法），最好的补偿方法为二维运动重力平衡结合变力（重力误差和惯性力误差同步补偿法），固定重力平衡结合变力与二维运动重力平衡结合定力在不同条件下各有优缺点.二维运动重力平衡相比固定重力平衡在变力补偿相比定力补偿时的误差补偿效果的总体提升比例显著增加.

参考文献

[1] Nguyen V L， Lin C Y， Kuo C H. Gravity compensation design of delta parallel robots using gear-spring modules［J］. Mechanism and Machine Theory， 2020， 154： 104-135.

[2] Han Z Y， Wang H， Shao Z X， et al. Finite difference method-based calculation of gravity deformation curve for the large-span beam of heavy-duty vertical lathe［J］. Advances in Mechanical Engineering， 2016， 8（4）： 1-8.

[3] Huang L T， Yamada H， Ni T， et al. A master-slave control method with gravity compensation for a hydraulic teleoperation construction robot ［J］. Advances in Mechanical Engineering， 2017， 9（7）： 1-11.

[4] Miao E M， Liu Y， Liu H， et al. Study on the effects of changes in temperature-sensitive points on thermal error compensation model for CNC machine tool［J］. International Journal of Machine Tools and Manufacture， 2015， 97（7）： 50-59.

[5] Andolfatto L， Mayer J R R， Lavernhe S. Adaptive Monte Carlo applied to uncertainty estimation in five axis machine tool link errors identification with thermal disturbance［J］. International Journal of Machine Tools and Manufacture， 2011， 51（7/8）： 618-627.

[6] Qiao Y， Chen Y P， Yang J X， et al. A five-axis geometric errors calibration model based on the common perpendicular line （CPL） transformation using the product of exponentials （POE） formula［J］. International Journal of Machine Tools and Manufacture， 2017， 118： 49-60.

[7] Huang Y B， Fan K C， Lou Z F， et al. A novel modeling of volumetric errors of three-axis machine tools based on Abbe and Bryan principles［J］. International Journal of Machine Tools and Manufacture， 2020， 151： 103527.

[8] Abdulshahed A M， Longstaff A P， Fletcher S， et al. Thermal error modelling of a gantry-type 5-axis machine tool using a grey neural network model［J］. Journal of Manufacturing Systems， 2016， 41（8）： 130-142.

[9] Fan J W， Tao H H， Pan R， et al. Optimal tolerance allocation for five-axis machine tools in consideration of deformation caused by gravity［J］. The International Journal of Advanced Manufacturing Technology， 2020， 111（1）： 13-24.

[10] Jiang Z X， Song B， Zhou X D， et al. Single setup identification of component errors for rotary axes on five-axis machine tools based on pre-layout of target points and shift of measuring reference［J］. International Journal of Machine Tools and Manufacture， 2015， 98（2）： 1-11.

[11] Zhong L， Bi Q Z， Wang Y H. Volumetric accuracy evaluation for five-axis machine tools by modeling spherical deviation based on double ball-bar kinematic test［J］. International Journal of Machine Tools and Manufacture， 2017， 122（4）： 106-119.

[12] Xing K L， Achiche S， Mayer J R R. Five-axis machine tools accuracy condition monitoring based on volumetric errors and vector similarity measures［J］. International Journal of Machine Tools and Manufacture， 2019， 138（1）： 80-93.

[13] Li Y， Zhao J， Ji S J. A reconstructed variable regression method for thermal error modeling of machine tools［J］. The International Journal of Advanced Manufacturing Technology， 2017， 90（9）： 3673-3684.

[14] Xiang S T， Altintas Y. Modeling and compensation of volumetric errors for five-axis machine tools［J］. International Journal of Machine Tools and Manufacture， 2016， 101（5）： 65-78.

[15] Blaser P， Pavliček F， Mori K， et al. Adaptive learning control for thermal error compensation of 5-axis machine tools［J］. Journal of Manufacturing Systems， 2017， 44（2）： 302-309.

[16] Li Q， Yu T B， Wang Z X， et al. Design of optical free-form surface milling machine based on mechanical shunt and dynamic analysis［J］. Applied Sciences， 2021， 11（24）： 11764.

基金资助

国家重点研发计划项目(2017YFA0701201)