自然超空泡航行体回转运动数值研究_王瑞.pdf

返回 相似 举报
自然超空泡航行体回转运动数值研究_王瑞.pdf_第1页
第1页 / 共6页
自然超空泡航行体回转运动数值研究_王瑞.pdf_第2页
第2页 / 共6页
自然超空泡航行体回转运动数值研究_王瑞.pdf_第3页
第3页 / 共6页
自然超空泡航行体回转运动数值研究_王瑞.pdf_第4页
第4页 / 共6页
自然超空泡航行体回转运动数值研究_王瑞.pdf_第5页
第5页 / 共6页
点击查看更多>>
资源描述:
􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈 􀪈 振 动 与 冲 击 第 39 卷第 21 期JOURNAL OF VIBRATION AND SHOCKVol. 39 No.21 2020 基金项目 陕西省自然科学基金资助项目(2019JQ⁃944) 收稿日期 2019 -04 -08 修改稿收到日期 2019 -08 -08 第一作者 王瑞 男,博士,高级工程师,1984 年生 通信作者 刘传龙 男,博士,工程师,1988 年生 自然超空泡航行体回转运动数值研究 王 瑞1, 刘传龙2, 赵三飞2, 祁晓斌1 (1. 西北机电工程研究所, 陕西 咸阳 712099; 2. 中国舰船研究设计中心, 武汉 430064) 摘 要针对超空泡航行体领域亟待明确的航行体机动过程中空泡形态和流体动力特性开展数值模拟研究。 基 于 Logvinovich 独立扩张原理建立了离心力作用下的回转运动理论模型和用于求解超空泡航行体回转运动的三维数值计 算模型,通过两种计算模型下获得的计算结果进行对比,验证了数值模拟方法的有效性。 利用数值模拟方法研究了离心 力“拉直”作用对空泡内航行体流体动力特性的影响规律,分析了不同回转角速度下航行体表面压力分布特性。 结果表 明随着回转角速度增加,空泡轴线的外移量增大,离心力拉直空泡作用越强烈,将会对回转半径内外两侧空泡的发展的 产生较大影响;升力系数和力矩系数呈现“正弦”形态分布,当回转角速度小于 1 rad/ s 时,产生的俯仰力矩对航行体起到 “抬头”作用,反之,则产生的俯仰力矩对航行体起到“低头”作用。 关键词 回转运动; 离心力; 泡内压力; 数值模拟 中图分类号 TJ630. 1 文献标志码 ADOI10. 13465/ j. cnki. jvs. 2020. 21. 009 Numerical simulation for rotating motion of a natural super⁃cavitating vehicle WANG Rui1, LIU Chuanlong2, ZHAO Sanfei2, QI Xiaobin1 (1. Northwest Institute of Mechanical and Electrical Engineering, Xianyang 712099, China; 2. China Ship Development and Design Center, Wuhan 430064, China) Abstract Here, numerical simulation was performed for cavitation shape and hydrodynamic characteristics of a natural super⁃cavitating vehicle in its maneuver process. Based on Logvinovich’s principle of independent expansion, the theoretical model of its rotating motion under action of centrifugal force and the 3⁃D numerical calculation model for solving rotating motion of the super⁃cavitating vehicle were established. By comparing calculation results of these two models, the effectiveness of the numerical simulation method was verified. The numerical simulation method was used to study influence laws of centrifugal force’s “straightening” action on hydrodynamic characteristics of the vehicle in cavitation, and analyze pressure distribution characteristics on the surface of the vehicle at different rotating angular velocities. Results showed that with increase in rotating angular velocity, cavitation axis’ s external displacement increases, the centrifugal force’s straightening cavitation action becomes stronger to have greater influence on development of cavitation on both sides of the vehicle’s rotation radius; lift coefficient and moment coefficient present a sinusoidal distribution; when rotating angular velocity is less than 1 rad / s, the generated pitching moment has a “head up” effect on the vehicle; when rotating angular velocity is larger than 1 rad / s, the generated pitching moment has a " head down" effect on the vehicle. Key words rotating motion; centrifugal force; cavitation internal pressure; numerical simulation 利用超空泡技术[1⁃5]突破水下航行速度屏障已经 成为国内外学者研究的热点问题。 受实验条件和相关 技术的制约,国内外有关超空泡流动的研究大多是围 绕着水下航行体直线运动展开的不同空化器头型超 空化流状态下的实验研究[6],航行体锥段线型对自然 空泡形态及减阻特性的影响研究[7],航行体直航状态 下的超空泡流发展过程以及空泡形态变化特性研 究[8]。 对于超空泡航行体回转运动,多体现在理想状 态下的空泡形态研究[9⁃10]以及纵平面内机动运动特性 分析[11]和控制[12]的研究。 而欧、美等国在超空泡机动 航行领域[13⁃14]进行了大量研究并取得了突破,其中最 具代表性的就是德国“BARRACUDA”试验型超空泡鱼 雷[15],该鱼雷通过控制头部空化器偏转,实现高速转弯 运动,其转弯半径小于 60 m。 在超空泡状态下,航行体 回转过程中空泡形态受到离心力影响而变化以及相应 流体动力特性的分析是超空泡航行体控制与结构设计 研究的基础。 本文应用商业软件 FLUENT,采用多相流模型和 SST(Shear Stress Transport)湍流模型,建立了用于求解 超空泡航行体回转运动的三维数值计算模型,对水下 超空泡航行体定深度回转运动进行了数值模拟。 首先 基于空泡截面独立扩张原理以及伯努利方程分析了离 心力对于空泡形态的影响及通过 CFD 数值模型计算获 得同等条件下的空泡形态,相互支撑验证了计算结果 的有效性;然后基于 CFD 数值模型分析了回转运动中 超空泡航行体流体动力特性,并探讨泡内压力分布对 航行体流体动力的影响。 1 超空泡航行体回转运动理论模型 1. 1 Logvinovich 空泡截面独立扩张原理 Logvinovich 空泡截面独立扩张原理[16]定义如下 对于高速运动物体形成的超空泡,空泡的每一个固定 的横截面都相对于物体中心运动轨迹按相同的规律扩 张,该扩张规律依赖于物体通过所论截面所在平面时 刻的条件物体的尺度、速度、阻力及无限远处与空泡 内的压力差,而与物体此时刻之前或之后的运动几乎 无关。 1. 2 离心力对空泡外形影响 如图 1 所示,当超空泡航行体作回转运动时,以 O0x0y0z0为地面坐标系,坐标原点 O0位于回转运动中 心,O0x0z0平面为水平面;orl 为空泡体坐标系,坐标原 点 o 位于空泡起始截面中心,ol 坐标轴与空泡对称轴 线重合,指向空泡尾部,or 坐标轴沿空泡径向,指向按 直角坐标系右手法则确定。 根据图 1,空泡内表面(长 虚线)与外表面(短虚线)的点因回转半径及线速度不 同,内表面点速度小,外表面点速度大。 根据伯努利方 程,空泡内表面压力大,外表面压力小,内外压力差形 成沿回转半径方向合力(离心力),从而对空泡外形[17] 产生影响;空泡内气体具有一定质量,在回转过程中会 产生惯性力,方向沿回转半径方向,对空泡外形亦产生 一定影响。 此处设超空泡航行体以角速度 ω、回转半径 R 作 回转运动,用动量定理建立离心力对空泡外形影响的 计算公式。 在空泡体坐标系 orl 中考虑半径 rk(l)单位 长度空泡沿 R 方向力的平衡,其动量变化量 πρr2 kωRrH 应等于空泡离心力 Fr,即 r H(l) = Fr(l) πρωRr2 k(l) (1) drH dl = drH dt dt dl = r H 1 ωR = Fr(l) πρω2R2r2 k(l) (2) 沿 l 轴积分 rH(l) = 1 πρω2R2∫ l 0 Fr(ζ) r2 k(ζ) dζ(3) 式中rH(l)为空泡中弧线的外移量;rk(l)为 l 处空泡截 面半径;ρ 为水密度。 图 1 离心力作用下空泡变形分析 Fig. 1 Analysis of cavity deformation under the action of centrifugal force 为积分式(3)得到 rH(l),需先确定离心力 Fr表达 式。 Fr由两部分组成,一部分是由空泡内表面和外表 面压力差形成的离心力 Frp,另一部分由空泡质量产生 的惯性力 Frm。 为计算 Frp,在 rk(l)截面建立局部坐标 系 o1ξη(见图 1),o1ξ 轴垂直于航行体回转运动平面, o1η 轴沿航行体回转半径方向,指向外侧。 把空泡表面 分为内侧和外侧,o1ξ 轴之上为外侧。 空泡外、内侧表 面点在 o1ξη 坐标系中的坐标为 ξkw= rkcos θ, ηkw= rksin θ ξkn= rkcos θ, ηkn= - rksin θ 0 ≤ θ ≤ 180 (4) 空泡外侧和内侧表面点至航行体回转运动中心的 距离 LRw、LRn分别为 LRw=(R + ηkw)2+ ξ2 kw= (R + rksin θ)2+ (rkcos θ)2 LRn=(R + ηkn)2+ ξ2 kw= (R - rksin θ)2+ (rkcos θ)2 (5) 空泡内、外侧表面点的速度 vkn、vkw分别为 vkw= ωLRw vkn= ωLRn { (6) 根据伯努利方程,空泡外侧表面压力 pkw和内侧表 面压力 pkn之间关系为 pkw(l) + 1 2 ρv2 kw(l) = pkn(l) + 1 2 ρv2 kn(l) (7) 内外表面压力差 Δpk为 Δpk= pkn- pkw= 1 2 ρ(v2 kw - v2 kn) = 2ρω2sin θRrk(l)(8) 66振 动 与 冲 击 2020 年第 39 卷 沿 ξ、η 方向的压力分量分别为 Δpkξ= Δpkcos θ = ω2ρsin 2θRrk(l) Δpkη= Δpksin θ = 2ρω2sin2θRrk(l) { (9) 离心力可以通过单位长度空泡表面压力积分获 得。 沿 ξ 方向的压力左侧表面和右侧表面对应点大小 相等方向相反,沿表面积分合力为零。 沿 η 方向的压 力积分即为空泡内表面和外表面压力差形成的离 心力。 Frp(l) =∫ πrk 0 Δpkηds = 2ρω2Rr2 k∫ π 0 sin2θdθ = ρω2RVk(l)(10) 式中,Vk为半径为 rk(l)的单位长度空泡体积。 空泡质量形成的单位长度空泡离心力 Frm为单位 长度空泡质量的惯性力 Frm= ρqVk(l)ω2R(11) 式中,ρq为空泡内气体的密度。 总的离心力为式(10) 与式(11)之和 Fr= ρω2RVk(l) + ρqVk(l)ω2R = (ρ + ρq)ω2RVk(l)(12) 把式(12)代入式(3),得到回转运动离心力引发的 超空泡中弧线外移量计算公式如下 rH(l) = 1 πρR(ρ + ρq)∫ l 0 Vk(ζ) r2 k(ζ) dζ(13) 2 超空泡航行体回转运动研究方法 2. 1 计算模型及基本设置 研究中超空泡航行体模型采用了圆盘为空化器, 具体形状如图 2 所示。 航行体由的圆盘空化器、圆锥 过渡段和圆柱弹身段组成,航行体全长 L =3. 2 m,空化 器直径 Dn=75 mm,圆柱段直径 D = 213 mm。 整个计 算域均采用六面体结构化网格,如图 3 所示。 图 2 计算模型 Fig. 2 Computational model 图 3 航行体周围网格划分 Fig. 3 Grid division around vehicle 计算域及边界条件的设定如图 4 所示计算域的 进口和出口分别定义为速度入口和压力出口,外边界 采用自由滑移壁面,消除壁面黏性影响,下表面为对 称面。 图 4 计算域及边界条件 Fig. 4 Computational domain and boundary conditions 2. 2 数值计算方法验证 图 5 给出了空化器直径 Dn= 75 mm,空化数 σ = 0. 019 79,航速 v =100 m/ s,角速度分别为 ω =0. 5 rad/ s、ω =1 rad/ s、ω =1. 5 rad/ s 和 ω =2 rad/ s 条件时考虑 离心力对超空泡外形影响与数值模拟计算结果的对 比。 图 6 给出了 ω = 1 rad/ s 时,空泡中弧线外移量沿 中弧线分布。 图 7 给出了对应四种不同角速度下空泡 中弧线最大偏移量的变化规律。 结合图 5、图 6、图 7 可知,超空泡航行体回转运动 时离心力使空泡中弧线沿径向向外侧偏移,偏移量由 头部至尾部逐渐增大,在尾端达到最大值,具有使“月 牙形”空泡拉直作用,离心力所致中弧线外移量随回转 角速度增大而增大;从四种回转运动条件下空泡整体 形态对比可知,CFD 数值模拟结果与离心力影响下的 理论结果基本一致,验证了所建超空泡航行体回转运 动数学模型以及数值模拟方法的有效性。 3 超空泡航行体回转运动特性研究 3. 1 回转运动中航行体流体动力特性 在研究回转运动中超空泡航行体流体动力特性 时,分析了空化数 σ =0. 019 79,航速 v =100 m/ s,角速 度分别为 ω =0. 2 rad/ s、ω =0. 5 rad/ s、ω =1 rad/ s、ω = 1. 5 rad/ s 和 ω =2 rad/ s 时的超空泡航行体回转运动特 性,图 8 给出了不同回转角速度条件下航行体表面的 空泡形态。 对比图 8 中各图可知,随着回转角速度的增加,包 裹航行体的空泡整体弯曲程度在不断增加,但由于离 心力对空泡的拉直作用,使得空泡下表面与航行体下 表面的距离减小,空泡上表面与航行体上表面的距离 增加,航行体下表面越来越贴近空泡壁面,但未出现沾 湿情况,与图 5 中圆盘空化器空泡形态变化趋势一致。 图9和图10列出了超空泡航行体的升力系数和 76第 21 期王瑞等 自然超空泡航行体回转运动数值研究 (a) R =200 m, ω =0.5 rad/ s (b) R =100 m, ω =1 rad/ s (c) R =66. 7 m,ω =1.5 rad/ s (d) R =50 m, ω =2 rad/ s 图 5 两种方法空泡形态结果对比 Fig. 5 Comparison of two kinds of methods for cavity shape 图 6 ω =1 rad/ s 空泡中弧线外移量轴向分布曲线 Fig. 6 Axial distribution curve of Arc outward shift in cavitation 图 7 空泡中弧线最大外移量随回转角速度的变化曲线 Fig. 7 The curve of the maximum outward shift of arc in cavitation vs. the angular velocity of rotation 俯仰力矩系数随回转角速度的变化曲线,其中俯仰力 矩是取航行体升力对其浮心的力矩。 由图 9 和图 10 可知,超空泡航行体的升力及俯仰 力矩主要是由弹身产生,空化器和尾端面产生的流体 动力系数很小,这是因为回转运动过程中,空化器迎流 (a) ω =0.2 rad/ s (b) ω =0.5 rad/ s (c) ω =1 rad/ s (d) ω =1.5 rad/ s (e) ω =2 rad/ s 图 8 不同回转角速度下的空泡形态 Fig. 8 Cavity shape under different rotation angles 图 9 升力系数随回转角速度变化曲线 Fig. 9 The lift coefficient versus the rotation angles curve 面与速度方向垂直,不产攻角,而弹尾端面始终处于空 化区内,其升力主要是靠气动产生,产生的非对称升力 系数很小;自然超空泡航行体回转状态下的升力和力 矩系数随回转角速度的变化呈“正弦”变化,在 ω =0. 5 86振 动 与 冲 击 2020 年第 39 卷 图 10 力矩系数随回转角速度变化曲线 Fig. 10 The moment coefficient versus the rotation angles curve rad/ s 和 ω =1. 5 rad/ s 时,升力系数和力矩系数分别达 到正向和负向最大值,当回转角速度小于 1 rad/ s 即回 转半径大于 100 m 时,产生的俯仰力矩对航行体起到 “抬头”作用,保持航行体转“大弯”直行的趋势,而回 转半径小于 100 m 时,产生的俯仰力矩对航行体起到 “低头”作用,使得航行体呈迅速“转弯”趋势。 3. 2 泡内压力对航行体流体动力特性的影响 为研究泡内压力对航行体流体动力特性的影响。 图 11 给出了超空泡航行体回转运动状态下泡内压力 沿航行体轴向的压力测点。 其中,在航行体表面布置 65 个压力监测面,每个截面布置 5 个压力测点。 航行 体浮心位置处为坐标原点,截面 1 距航行体头部 15 mm,截面65 与航行体尾部重合,各截面间距50 mm,共 计设置了 325 个测点。 图 11 压力测点布置 Fig. 11 Pressure measuring points 图 12 列出了不同回转角速度下,不同弹身位置处 压力的分布规律。 由图 12 可知,航行体头部和尾部处 的压力为 2 367 Pa,处于完全空化状态,与图 8 中的空 泡形态相对应;弹身段(对应坐标位置区间[ - 1. 5, 1])的不同轴向位置和不同周向位置空化状态不同,所 受压力也不同,且受力主要集中在浮心前,由于航行体 上下表面压力分布的不同形成了作用在航行体上的升 力和俯仰力矩;当 ω = 0. 2 rad/ s、ω = 1 rad/ s 和 ω = 2 rad/ s 时,航行体受到的泡内压力在升力方向基本抵 消,浮心前后所受俯仰力矩也基本抵消,ω = 0. 5 rad/ s 时,弹身下部泡内压力明显大于上部泡内压力,且浮心 前泡内压力大于浮心后泡内压力,使得弹身受到正向 的升力和抬头力矩;ω =1. 5 rad/ s 时,弹身上部泡内压 力明显大于下部泡内压力,且浮心前泡内压力大于浮 心后泡内压力,使得弹身受到负向的升力和低头力矩。 (a) ω =0.2 rad/ s (b) ω =0.5 rad/ s (c) ω =1 rad/ s (d) ω =1.5 rad/ s (e) ω =2 rad/ s 图 12 空泡内压力分布曲线 Fig. 12 Curves of cavity internal pressure 4 结 论 本文采用数值模拟方法,分析航行体机动过程中 空泡形态和泡内航行体流体动力特性。 研究中,建立 了离心力作用下的回转运动理论模型和用于求解超空 泡航行体回转运动的三维数值模型,计算并获得了不 96第 21 期王瑞等 自然超空泡航行体回转运动数值研究 同回转角速度下的空泡形态和泡内压力。 所得主要结 论如下 (1) 航行体回转运动过程中,空泡外形受到离心 力的拉直作用而呈“月牙状”形态;随着回转角速度增 大,空泡轴线的外移量增大,发展趋势由空化器头部至 尾部逐渐增大,在尾部达到最大值。 (2) 由于空泡轴线的向外偏移导使得航行体回转 半径内侧表面最先出现沾湿,产生的升力和俯仰力矩 主要由弹身提供,且升力系数和力矩系数随回转角速 度的变化呈现“正弦”形态分布规律。 (3) 回转运动使得航行体两侧压力呈非对称分 布,最大压力峰值均出现在弹身质心前,形成了作用在 航行体上的升力和俯仰力矩使得航行体保持原运动轨 迹的趋势。 参 考 文 献 [ 1] ASHLEY S. Warp drive underwate[J]. Scientific American, 2001, 284(5)70⁃79. [ 2] SAVCHENKO Y N.Experimentalinvestigationof supercavitating motion of bodies[R]. Brussels VKI Special Course on Supercavitating Flows, 2001. [ 3] 曹伟,魏英杰,王聪,等. 超空泡技术现状、问题与应用 [J]. 力学进展, 2006, 36(4)571⁃579. CAO Wei, WEI Yingjie, WANG Cong, et al. Current status, problems and applications of supercavitation technology[J]. Advances in Mechanics, 2006, 36(4)571⁃579. [ 4] 杨莉,张庆明. 超空泡技术的应用现状和发展趋势[J]. 战术导弹技术, 2006(5)6⁃10. YANG Li,ZHANGQingming.Currentapplicationand perspectives on supercavitation technology research [ J ]. Tactical Missile Technology, 2006(5)6⁃10. [ 5] 肖昌美,李恒,彭佩. 国外水面舰艇鱼雷防御系统发展现 状及趋势[J]. 鱼雷技术, 2014, 22(2)150⁃156. XIAO Changmei, LI Heng, PENG Pei.Development of torpedo defence systems of foreign surface ships[J]. Torpedo Technology, 2014, 22(2)150⁃156. [ 6] 邓飞,张宇文,陈伟政, 等. 头形对细长体超空泡生成与外 形影响的实验研究[J]. 西北工业大学学报, 2004, 22 (3)269⁃273. DENG Fei, ZHANGYuwen,CHENWeizheng,etal. Experimental investigation on the incipiency and the shape of supercavity for slender bodies with different headforms[J]. Journal of Northwestern Polytechnical University, 2004, 22 (3)269⁃273. [ 7] 付英杰. 航行体锥段线型对自然空泡特性影响的数值模 拟[J]. 舰船科学技术, 2010, 32(12)126⁃130. FU Yingjie. Numerical research on influence of link section contour on natural cavity characteristics[J]. Ship Science and Technology, 2010, 32(12)126⁃130. [ 8] 易文俊,李铁鹏,熊天红,等. 水下高速航行体自然超空泡 形态特性仿真研究[J]. 南京理工大学学报(自然科学 版), 2009, 33(3)330⁃334. YI Wenjun, LI Tiepeng, XIONG Tianhong, et al. Simulation on natural supercavitation characteristics of underwater high⁃ speed vehicle[J]. Journal of Nanjing University of Science and Technology(Natural Science), 2009, 33(3)330⁃334. [ 9] ZHANG Guang, YU Kaiping, ZHOU Jingjun.Numerical research on ventilated supercavity shape and flow structure in the turning motion[J]. Journal of Ship Mechanics, 2011, 15 (12)1335⁃1343. [10] 王复峰,王国玉,张敏弟, 等. 带空化器回转体空化流场 研究[J]. 哈尔滨工程大学学报, 2015, 36(7)959⁃964. WANG Fufeng, WANG Guoyu, ZHANG Mindi, et al. Study on cavitatingflowsaroundanaxisymmetricbodywith cavitators[ J].Journal of Harbin Engineering University, 2015, 36(7)959⁃964. [11] 胡峰,罗凯,李淼, 等. 机动超空泡水下航行器纵平面运 动特性 分 析 [ J].计 算 测 量 与 控 制, 2011, 19 (7) 1643⁃1646. HU Feng, LUO Kai, LI Miao, et al. Motion characteristics on vertical plane of maneuvering supercavitating underwater vehicle[J]. Computer Measurement and Control, 2011, 19 (7)1643⁃1646. [12] LI Daijin, LUO Kai, ZHANG Yuwen, et al.Studies on fixed⁃depth control of supercavitating vehicles [ J].Acta Automatica Sinica, 2010, 36(3)421⁃426. [13] LIN G, BALACHANDRAN B, ABED E H. Dynamics and control of supercavitating vehicles[J]. Journal of Dynamic Systems, Measurement, and Control, 2008, 130 ( 2 ) 0210031⁃02100311. [14] KIRSCHNER I N, KRING D C, STOKES A W, et al. Control strategies for supercavitating vehicles[J]. Journal of Vibration and Control 2002, 8(2)219⁃242. [15] 从敏,刘乐华. 德国 BARRACUDA 超空泡高速水下导弹的 制导与控制[J]. 飞航导弹,2007,37(5)38⁃43. CONG Min, LIU Lehua. Guidance and control of high⁃speed supercavitation underwater missile of Germany[J]. Winged Missiles Journal, 2007,37(5)38⁃43. [16] 杨洪澜,史文谱. 超空泡横截面独立膨胀原理解析及其应 用[J]. 计算力学学报, 2014, 31(2)277⁃280. YANG Honglan, SHI Wenpu. Analysis of supercavity cross section independent expansion principle and its applications [J]. Chinese Journal of Computational Mechanics, 2014, 31 (2)277⁃280. [17] 李雨田,张宇文. 空泡外形对超空泡航行器机动回转运动 影响[J]. 振动与冲击,2014,33(16)165⁃170. LI Yutian, ZHANG Yuwen. The influence of cavity shape on maneuvering rotational movement of supercavitating vehicle [J].Journal of Vibration and Shock, 2014, 33 (16) 165⁃170. 07振 动 与 冲 击 2020 年第 39 卷
展开阅读全文

资源标签

最新标签

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

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

矿业文库合伙人QQ群 30735420