机床—主轴耦合系统动力学建模与模型修正.pdf

返回 相似 举报
机床—主轴耦合系统动力学建模与模型修正.pdf_第1页
第1页 / 共7页
机床—主轴耦合系统动力学建模与模型修正.pdf_第2页
第2页 / 共7页
机床—主轴耦合系统动力学建模与模型修正.pdf_第3页
第3页 / 共7页
机床—主轴耦合系统动力学建模与模型修正.pdf_第4页
第4页 / 共7页
机床—主轴耦合系统动力学建模与模型修正.pdf_第5页
第5页 / 共7页
点击查看更多>>
资源描述:
第 4 8 卷第 3期 2 O 1 2 年 2 月 机械工程学报 J OURNAL OF MECHANI CAL ENGI NEE G VO1 . 48 N O. 3 F e b . 201 2 DoI l O . 39 01 , JM E. 2 01 2 . 0 3 . 0 8 8 机床一主轴耦合系统动力学建模与模型修正水 曹宏瑞何正嘉 西安交通大学机械制造系统工程 国家重点实验室西安7 1 0 0 4 9 摘要复杂精密零件的加工对机床装备的性能不断提出新的挑战。为了保证零件加工精度,有必要在设计阶段对机床整机的 动态特性进行精确预测和评估。以机床主轴系统为研究对象,介绍一种机床一主轴系统耦合动力学模型,并指出该模型中存 在结合面动态特性参数的辨识问题。为了修正该耦合模型,提出~种基于频率响应函数的有限元模型修正技术。采用包含平 动自由度 De g r e e o f fr e e d o m, DO F ,转动 DO F以及平动与转动耦合 DO F的刚度矩阵来描述结合面的动态特性,通过测量结 合面处的动态响应数据,辨识出高速主轴与机床之间结合面的刚度参数,从而达到模型修正的目的。试验结果表明,修正后 的机床一主轴耦合模型反映出了机床结构对主轴动态特性的影响,能较好地预测主轴安装到机床上以后的动态特性,从而为 高性能机床的数字化设计与制造提供理论依据。 关键词机床主轴耦合模型修正结合面参数辨识 中图分类号T G 5 0 2 Dy n a mi c M o d e l i n g a n d M o d e l Upd a t i ng o f Co up l e d S y s t e m s b e t we e n ‘ M a c h i ne T o o l a n d I t s S pi n dl e C AO H o n g r u i HE Z h e n g j i a S t a t e K e y L a b o r a t o r y f o r Ma n u f a c t u r i n g S y s t e ms E n g i n e e r i n g , Xi ’ a n J i a o t o n g U n i v e r s i t y , X i ’ a n 7 1 0 0 4 9 Ab s t r a c t T h e ma c h i n i n g o f c o mp l e x a n d p r e c i s e p a r t s c o n t i n u a l l y b ri n g s n e w c h a l l e n g e s t o t h e p e rfo r ma n c e o f ma c h i n e t o o l s . T o e n s u r e t h e ma c h i n i n g a c c ura c y o f p a r t s , i t ’ S n e c e s s a r y t o p r e d i c t t h e d y n a mi c b e h a v i o r o f the wh o l e ma c h i n e t o o l i n the d e s i g n s t a g e . Ai me d a t ma c h i n e t o o l s p i n d l e s y s t e ms ,a c o u p l i n g mo d e l b e t we e n the s p i n d l e a n d ma c h i n e t o o l i s i n t r o d u c e d and t h e e x i s t e d p r o b l e m o f t h e j o i n t d y n a mi c s i d e n t i fi c a t i o n i n thi s mo d e l i s a l s o p o i n t e d o u t . I n o r d e r t o u p d a t e t h i s c o u p l i n g mo d e l , a fi n i t e mo d e l u p d a t i n g t e c h n i q u e b a s e d o n fr e q u e n c y r e s p o n s e f u n c t i o n i s p r o p o s e d . T h e j o i n t d y n a mi c s i s mo d e l e d a s a s t i f f n e s s ma t r i x w h i c h c o n t a i n s t r ans l a t i o n a l d e g r e e o f fr e e d o m D OF , r o t a t i o n a l DO F and c o u p l e d D OF . D y n a mi c r e s p o n s e d a t a a r e me a s u r e d t O i d e n t i f y t h e j o i n t d y n am i c s b e t we e n the s p i n d l e a n d ma c h i n e t o o l a n d t h e r e f o r e the a n a l y t i c a l mo d e l i s u p d a t e d . T h e e x p e r i me n t a l r e s u l t s s h o w t h a t the u p d a t e d c o u p l i n g mo d e l b e t we e n t h e s p i n d l e and ma c h i n e t o o l C a l l r e a s o n a b l y r e f l e c t the e ffe c t s o f t h e ma c h i n e t o o l s t r u c t u r e o n t h e s p i n d l e a n d p r e d i c t d yn am i c p r o p e r t i e s o f the s p i n d l e s y s t e m a f t e r mo u n t i n g , wh i c h c a ll p r o v i d e t h e o r e t i c a l p r o o f s f o r t h e d i g i tal i z e d d e s i gn a n d ma n u f a c t u r e o f h i g h p e r f o r ma n c e ma c h i n e t o o l s . Ke y wo r d s Ma c h i n e t o o l s p i n d l e Co u p l i n g Mo d e l u p d a t i n g J o i n t s P a r a me t e r i d e n t i fic a t i o n 0 前言 航空、航天、能源和国防等领域对复杂精密零 件 日 益增长的需求,促使数控机床的精度、效率、 可靠性等性能指标向更高的水平发展 。对机床装备 进行精确的数字化建模是研究其动态行为的有效手 段。 作为高性能数控机床核心部件之一的高速主轴, 国家自然科学基金 5 1 1 0 5 2 9 4 、国家重点基础研究发展计划 9 7 3计划, 2 0 0 9 C B 7 2 4 4 0 5 1 和 中央高校 基本科研 业 务费专项 资金 资助项 目。 2 0 1 1 0 4 1 3收到初稿,2 0 1 1 1 2 1 5 收到修改稿 一 直是机床制造领域内的研究重点。目前大部分针 对主轴的研究,都将主轴从机床系统 中分离出来, 单独对其进行有限元建模,研究 自由或完全固定边 界条件下转子一轴承系统的动态特性,而没有考虑 机床结构对主轴动态特性的影响 J 。 然而在实际应 用中,当主轴被安装到机床上后,随着边界条件的 改变,主轴系统的动态性能将是主轴本身、主轴与 机床之 间结合面 如螺栓联接、 导轨等 及机床这三 部分结构动态特性的综合反映。在主轴与机床组成 的相互耦合的新系统中,机床会对主轴的动态特性 产生不可忽略的影响,仅研究自由边界条件下主轴 学兔兔 w w w .x u e t u t u .c o m 学兔兔 w w w .x u e t u t u .c o m 机械工程学报 第 4 8卷第 3期 ,q 平动与转动 自由度之间的耦合 刚 度和阻尼 卜下标,结合面 相应地 , 主轴一机床耦合模型 的刚度矩阵 对称 阵 和阻尼矩阵 对称阵 分别为 f 一 ] / 1 2 l l 【2 c H十c J / I 式中 ,露 主轴、机床 的刚度矩阵 c s ,c H 主轴、机床的阻尼矩阵 结合面的等 效刚度 和 阻尼矩 阵 c J 属于未 知 参数 ,代表主轴与机床之间的耦合特性 。C AO等 J 没有系统地对这些结合面参数进行辨识,而是采用 反复试验法来获得合理的值,缺乏实用性。因此, 需要研究相关的技术对这些结合面参数进行系统的 辨识 ,从而迅速准确地建立主轴一机床耦合模型。 2 有限元模型修正理论 2 . 1 概述 对于复杂结构 ,各子结构之间的耦合特性往往 是未知的,比如焊接点及螺钉联接等。这些不确定 因素导致建立的有限元模型精度低,应用效果差。 有限元模型修正的目的是将解析模型同试验模型相 结合,导出一个能够对机械结构的动态特性进行准 确而可靠预测 的有 限元模型 。 近 2 0年来,关于有限元模型修正技术的研究 十分活跃,英 国的 MOT T E R S HE AD 等 J 、美国的 HE ME Z等 J 以及法国的 A V R I L等 著名学者分别 从不同的角度相继对该技术 的发展进行了详细的综 述和讨论。在众多的有 限元模型修正方法 中,以灵 敏度为基础的迭代修正法与有限元模型的离散本质 一 致,只需要修正系统矩阵中的某几个单元或节点 的参数,就可以使有限元模型迅速收敛于与试验数 据良好相关的修正模型,因此备受青睐。灵敏度迭 代法的原理可用式 3 描述 S A u£ 3 式中 灵敏度矩阵 △ l | 修正参数的改变量 s 残余矢量, 即解析数据与试验数据之差 式 3 通常为超定方程, 可利用最小二乘法求解 修正参数的改变量△ Ⅳ。假设修正参数的初始值为 ,每一步迭代运算中,都会得到更新的修正参数 l I I I o A u,将其代入有限元模型中计算新的解析 解 ,进而得 到新的残余矢量 £;重复迭代过程直到 残余矢量为 0或者小于某一设定的阈值。根据残余 矢量的数据组成,又可以把该方法分为基于模态参 数和基于频率响应 函数的迭代法 。基于频率响应函 数 的迭代法最早 由德 国的N A T K E E 8 1 于 1 9 7 7 年提 出, 随后该方法成为研究热点,许多学者从不同角度对 其进行发展和改进。1 9 9 3年,比利 时的 L A MME NS 等 利用动态压缩技术解决有 限元模型同试验 自由 度数不一致 问题 ,合理地利用 了试验频率响应 函数 信息,然而修正后的模型只能近似描述结构的动态 特性 ,而且没有考虑结构的阻尼。2 0 0 9年,印度学 者 A RO R A等【 J 研究了阻尼系统的有限元模型修正 方法 , 利用扩展的频率响应函数法处理复模态 问题 , 提出先修正系统质量和刚度矩 阵,再修正阻尼矩阵 两步走的办法 ,然而对于实际结构,由于很难建立 起解析的阻尼矩阵, 使这些方法的应用受到 了限制。 虽然有限元模型修正技术 已取得了长足的进 步,然而大多数的研究仍停留在对简单结构的理论 仿真和试验阶段,距 离实际应用还有差距。测量 自 由度与有限元模型中 自由度的不一致、试验模态数 据不完整、复杂的阻尼特性和测量噪声等仍然是有 限元模型修正技术成功应用的主要障碍 。本文针对 大 型复 杂 阻尼 结构 的有 限元 模 型修 正 问题 , 以 L AMME NS等 J 的修正理论为基础,提 出一种基于 频率响应 函数的修正方法,以增强有限元模型修正 技术的实用性。 2 . 2 基于频率响应函数的模型修正理论 有限元模型修正技术的理论基础是机械结构 的解析模型可以通过适当的参数改变 ,使其解析模 态参数、频率响应 函数等更能符合试验测量结果, 从而更好地描述实际结构的动态特性。假设解析模 型 用上标 a表示 和实际结构 用上标x表示 的系统 运动方程如式 4 所示 { 芝 4, lz F 国 一 式中, 为频率, 和 ,分别为输出位移矢量和输 入力矢量。解析和实际结构的系统动刚度矩阵z可 以分别表示为 { Z a 0“k a 一- 02mm a ..1_jCOC ak CO OC 5 l z 一 m j 一 将解析模型与实际结构的残余力矢量定义为 s Z 一 J 6 若激励力为参考自由度 V 的单位力矢量,由式 学兔兔 w w w .x u e t u t u .c o m 2 0 1 2年 2月 曹宏瑞等机床一主轴耦合系统动力学建模与模型修正 9 1 6 推得关于 自由度 V、频率 残余力矢量为 Z 月 7 式中,矢量 L对于参考 自由度 等于 1 ,对于其他 自由度则为 0 。定义参考 自由度 1 , 的残余频率响应 函数矢量,即解析频率响应函数与试验频率响应函 数之差为 c H 缈 一日 8 可 以证 明, 8 Ⅳ 的一阶泰勒展开式有 限。然而更稳 定的线性化可 以根据 式 7 进行 ,以 日。 左 乘式 7 ,得 H s , v 9 式 9 关于修正参数 l f 一阶泰勒展开式为 ‰ ‰ 14o 善 鲁 10 ‰ ‰ ∑ 日 69 △ “ 10 式中 。 修正参数初始值 △ “ 第 i 个修正参数增量 式 1 0 按一种合理稳定的方式将所有 的试验信 息结合在一起 ,是基于频率响应函数的有 限元模型 修正迭代法运算的基础。令 0,则有 nay 嘻 鲁 。 11 l|。 ∑ 日 。 △ 0 11 式 1 1 可以简化为标准的迭代方程 S 0 A us 1 2 式 中,灵敏度矩阵 、修正参数变化 △ 和残余 矢量 如下 △ l f △ 1 △ “ 2 ⋯ △ 灵敏度矩阵式 1 3 0 0 残余力矢量对修正参数 的 导数为 鲁 一 [嚣 一 寄 堡Ouiau 1 4 f a f a f 、 、 由于阻尼作用机理的复杂性 ,难 以建立解析的 系统阻尼矩阵。因此在实际应用中,经常令阻尼矩 阵为零,式 1 4 变为 鲁 一 盖 m IH xc 15 a a“f 一 举例来说 ,如果修正参数为第 P 自由度的局部 刚度 和第 P 自由度 与第 q自由度之 间的耦合 刚度 A k p q ,则式 1 2 可具体为 日 0 一 H a, ⋯ , x 一 日 ; H 钾a 一 H a ,x 、0 一 ; H Ⅳ 日a 一 日 一 一 ⋯ 1 6 在式 1 6 中,灵敏度矩阵 由解析和试验测量的 频率响应函数构成,每一行均包含有相同的试验信 息 g x , 。从信息冗余 的角度来说,没 有必要利用所有行来构成灵敏度矩阵,即不需要计 算结构上所有响应点的残余频率响应函数。如果选 取激励 自由度 V,响应 自由度 W的残余频率响应函 数矢量 日 一 日 作 为迭代收敛的 目标函数,那么灵敏度矩阵中的第 W行将被抽取出 来 , 在 频 率 响 应 函 数 上 任 选 个 频 率 点 q, , ⋯, 可组成新的超定方程 H a q 0 1 一 q q 一 q H a 一 一 一 % 一 /-/ 5 % Ak舶 一 q 日 1 7 式 1 7 为基于频率响应函数的有限元模型修正 算法的核心。由式 1 7 N- j 见,不管实际结构多复杂, 有限元模型的维数有多高,若修正参数为两自由度 间的耦合刚度 , 则只需要测量相应 自由 度的频 率响应函数 和 即可构成灵敏度矩 学兔兔 w w w .x u e t u t u .c o m 机械工程学报 第 4 8卷第 3 期 阵。因此该算法有效地避免了有限元模型的自由度 数与试验测量 自由度数不一致时带来的影响,能利 用最少的试验数据进行模型修正,具有很强的实用 性和应用价值 。 式 1 5 中阻尼矩阵被忽略,并不意味着提出的 模型修正算法没有考虑实际结构的阻尼 。结构的一 般粘性阻尼难 以用解析模型去描述 ,但可以采用试 验模态分析法来辨识各阶模态的阻尼参数。将辨识 出的模态阻尼参数输入到有限元模型中,即可得到 阻尼系统的解析频率响应函数 。由于解析和试验频 率响应函数均为复数,将产生复数灵敏度矩阵和残 余矢量 ,而修正参数为实数,为了便于计算,式 3 可进一步写为 f R e S 、 f R e e 1 0 S lJ Im s J h n 8 j 0 8 h n 8 3 高速主轴一机床耦合模型修正 3 . 1 模型修正 下面将应用基于频率响应函数 的有 限元模型 修正法,对主轴一机床耦合模型中螺钉联接的结合 面参数 质量、刚度和 阻尼参数 进行辨识,从而达 到模型修正之 目的。对于螺钉联接,将其质量参数 忽略 ,而阻尼参数将利用试验模态法确定,因此只 需要辨识结合面的刚度参数。在有限元模型中,主 轴与机床在结合面处的刚度参数利用等效弹簧来表 示, 弹簧的两端分别 同节点 1 7 主轴箱 和节点 2 8 主 轴座 连接,如图 4 所示。 图4 主轴一机床耦合有限元模型 节 点 垂 轴 承 I 固 定 连 接 詈 轴 向 滑 动 一 轴 承 隔 圈 将主轴端部 有限元模型第 1个节点 上平动 自 由度的残余频率响应函数 激励 自由度v l ,响应 自 由度 w 1 .。 ∞ 一 日 作为修正算法的目标函数, 其中Ha 为解析频率 响应函数, 。 为试验测量频率响应函数。系统 刚度矩阵中需要修正或辨识的刚度参数为 △l f 。 2 式 中, . 3 . 3 3 , 4 . 3 3 , 4 . 3 4 和 5 . 5 5 , 6 . 5 5 , △ 6 - 5 6 分别为节点 1 7和 2 8上需要修正的平动刚 度,平动与转动耦合刚度及转动刚度。在频率响应 函数频段范围内随机选取频率点 , , ⋯, , 根据式 1 7 建立超定方程,并利用最小二乘迭代法 进行求解从而得到结合面的刚度参数值 。在建立灵 敏度矩阵时,需要测量结合面处的试验频率响应 函 数 , , ,。 或【 5 , , , 。 由 于 转 动 自由度 动态 响应 测量 复 杂且 成本 高 昂 ,采 用 R e c e p t a n c e c o u p l i n g技术 u 来估计转动 自由度的频 率响应函数 . 或 . ,限于篇幅 , 此处不 再赘述 。 下面给出主轴一 机床有 限元耦合模型修正 的 结 果 并 进 行 分 析 。 图 5 为 残 余 矢 量 范 数 II , lJ 的 收敛情况,经过 3 0次迭代后,残余误差的变化趋 于稳定。 _、 ● 量 i j I{ 】} 图5 残余频率响应函数的收敛情况 需要辨识的结合面参数即平动刚度 、 平动与 转动耦合刚度 和转动刚度 的变化趋势如图 6 a 6 c 所示。 结合面参数的初始值均为零, 经过 3 0 次运算后,得到修正值。 为了便于 比较,将结合面刚度参数的初始值和 修正值列于下表 。 学兔兔 w w w .x u e t u t u .c o m 学兔兔 w w w .x u e t u t u .c o m 机械工程学报 第 4 8卷第 3期 4 结论 1 提出了~种基于频率响应函数的有限元模 型修正技术。该修正算法只需要测量结合面处的动 态响应数据,大大减少了所需的试验信息,有效地 避免了有限元模型的自由度数与试验测量 自由度数 不一致带来的影响。 2 利用试验模态分析法确定结构的模态阻尼 参数,避免了模型修正时对一般粘性阻尼结构的复 杂建模,试验证明这是一种行之有效的方法,具有 很强的实用性。 3 采用包含平动 自由度,转动 自由度以及平 动与转动耦合 自由度 的刚度矩阵来描述结合面的动 态特性,利用有限元模型修正技术成功提取出了主 轴与机床结合面的刚度参数,达到了修正主轴一机 床耦合模型的目的。 参考文献 [ 1 ]L I H o n g q i ,S H I N Y C .I n t e g r a t e d d y n a mi c t h e r mo me c h a n i c a l mo d e l in g o f h i g h s p e e d s p ind l e s , p a r t 1 - mo d e l d e v e l o p me n t [ J ] . J o u r n a l o f Ma n u f a c t u ri n g S c i e n c e a n d E n g i n e e ri n g , T r a n s a c t i o n s o f t h e AS ME ,2 0 0 4 ,2 6 1 1 4 8 . 1 5 8 . [ 2 】G AG NO L V,B OU GA R R O U B C,R A Y P,e t a 1 . S t a b i l i t y - b a s e d s p i n d l e d e s i g n o p t i mi z a t i o n [ J ] . J o u r n a l o f Manu f a c t u r i n g S c i e n c e a n d E n g i n e e ri n g ,T r a n s a c t i o n s o f t h e AS ME, 2 0 0 7,1 2 9 2 1 4 0 7 - 4 1 5 . 【 3 】C AO Y u z h o n g , I N T A S Y Mo d e l i n g o f s p i n d l e b e a r i n g an d ma c h i n e t o o l s y s t e ms f o r v i r t u a l s i mu l a t i o n o f mi l l ing o p e r a t i o n s [ J ] . I n t e rna t i o n a l J o u r n a l o f Ma c h i n e T o o l s Manu f a c t u r e , 2 0 0 7, 4 7 9 1 3 4 2 . 1 3 5 0 . [ 4 】DH U P I AJ S ,P 0 w L KAB,UL S oYAG,e t a 1 . E ff e c t o f a n o n l i n e a r j O i n t o n the d y n a mi c p e r f o r ma n c e o f a ma c h i n e t o o l [ J ] . J o u r n a l o f Manu f a c t u r i n g S c i e n c e and E n g i n e e r i n g , T r ans a c t i o n s o f t h e AS ME ,2 0 0 7 ,1 2 9 5 9 43 9 50 . [ 5 】MO T T E R S H E A D J E ,F R I S WE L L M I . Mo d e l u p d a t i n g in s t r u c t u r a l d yna mi c s a s u r v e y [ J ] . J o u r n a l o f S o u n d and Ⅵb r a t i o n ,1 9 9 3 ,l 6 7 2 3 4 7 3 7 5 . [ 6 ]6 H E ME Z F M,DOE B L G S W R e v i e w and a s s e s s me n t o f mo d e l u p d a t i n g for n o n - l i n e a r ,t r a n s i e n t d ynam i c s [ J ] . Me c h a n i c a l S y s t e ms a n d S i g n a l P r o c e s s i n g , 2 0 0 1 , 1 5 1 4 5 . 7 4 . [ 7 ]A V R I L S , B O N NE T M, B R E T E L L EA S , e t a 1 . O v e r v i e w o f i d e n t i fi c a t i o n me t h o d s o f me c h a n i c a l p ara me t e r s b a s e d o n f u l l fi e l d me a s u r e me n t s [ J ] . E x p e ri me n t a l Me c h a n i c s , 2 0 0 8 ,4 8 4 3 8 1 . 4 0 2 . [ 8 】N A T KE H G D i e k o r r e k t u r d e s r e c h e n mo d e l l s e i n e s e l a s t o me c h ani s c h e n s y s t e ms mi t t e l s g e me s s e n e r e r z wu n - - g e n e r s c h wi n g u n g e n [ J ] . A r c h i v e o f Ap p l i e d Me c h ani c s , 1 9 7 7 ,4 6 3 1 6 9 . 1 8 4 . NA TKE H G n l e c o r r e c t i o n o f th e c o mp u t a t i o n a l mo d e l o f a s y s t e m u s i n g the me a s ure d e l a s t o f o r c e d o s c i l l a t 一 ’ i o n s [ J ] . Ar c h i v e o f A p p l i e d Me c h a n i c s ,1 9 7 7 ,4 6 3 1 6 9 . 1 8 4 . [ 9 ]I , A MME NS S ,H E Y L E N W, S AS P Mo d e l u p d a t i n g u s i n g e x p e rime n t a l f r e q u e n c y r e s p o n s e f u n c t i o n s Ca s e s t u d i e s [ C ] .P r o c e e d i n g s o f S tr u c t u r a l D ynam i c s Mo d e l l i n g , T e s t ,An a l y s i s an d Co rre l a t i o n, 1 9 9 3 1 9 5 - 2 0 4 . [ 1 0 】AR O R AV, S 1 NG H S P ,K U ND R AT K. Da mp e d mo d e l u p d a ti n g u s i n g c o mp l e x u p dat i n g p aram e t e r s [ J ] . J o u rna l o f S o u n d a n d V i b r a t i o n ,2 0 0 9 ,3 2 0 1 - 2 4 3 8 - 4 5 1 . [ 1 1 ]P A R K S S ,C HA E J . J o i n t i d e n t i fi c a t i o n o f mo d u l a r t o o l s u s i n g a n o v e l r e c e p t anc e c o u p l i n g me t h o d [ J ] . I n t e r n a t - i o n a l J o u r n a l o f Ad v a n c e d Ma n u f a c t u r i n g T e c h n o l o g y , 2 0 0 8 ,3 5 1 1 - 1 2 1 1 2 5 1 - 1 2 6 2 . 作者简介曹宏瑞 通信作者 ,男,1 9 8 2年出生,博士,讲师。主要研 究方 向为机械设备动力学建模与故障诊断。 E ma i l c h r mail .x j t u . e d u . c n 何正嘉,男,1 9 4 2年出生,教授,博士研究生导师 。主要研究方向机械 动态分析与设计、机械设备状态监测与故障诊断。 E ma i l ma i l .x j n 1 . e d u .c n 学兔兔 w w w .x u e t u t u .c o m
展开阅读全文

资源标签

最新标签

长按识别或保存二维码,关注学链未来公众号

copyright@ 2019-2020“矿业文库”网

矿业文库合伙人QQ群 30735420