通过响应面方法研究竹材加工工艺参数对声级和表面质量的影响 陈海阳(Haiyang Chen)、 迪特里希·巴克(Dietrich Buck)、 丁建文(Jianwen Ding)、 郭晓蕾(Xiaolei Guo) 朱兆龙(Zhaolong Zhu)

《Forests》:Bamboo Milling Process Parameters’ Influence on Sound Level and Surface Performance via Response Surface Methodology Haiyang Chen, Dietrich Buck, Jianwen Ding, Xiaolei Guo and Zhaolong Zhu

【字体: 大 中 小 】 时间:2026年05月02日 来源:Forests 2.5

编辑推荐:

  摘要:本研究探讨了关键铣削参数如何影响层压竹材加工过程中的切削噪声和表面质量。通过使用多因素最优响应面方法,量化了纤维方向(0–135°)、主轴转速(7000–10,000 r/min)、进给速度(0.5–2.0 m/min)和铣削深度(0.5–2.0 mm)的影响,共进行了25

  摘要:本研究探讨了关键铣削参数如何影响层压竹材加工过程中的切削噪声和表面质量。通过使用多因素最优响应面方法,量化了纤维方向(0–135°)、主轴转速(7000–10,000 r/min)、进给速度(0.5–2.0 m/min)和铣削深度(0.5–2.0 mm)的影响,共进行了25次实验。切削噪声以峰值声压级(SPL)为单位进行测量,范围从86.8到95.2 dB,随着纤维角度、进给速度和铣削深度的增加而显著增加,但对主轴转速的响应呈非线性。表面粗糙度(Sa)的范围为2.6到11.7 μm,主要受铣削深度的影响,其次是纤维方向和进给速度,纤维方向与主轴转速之间存在显著交互作用。二次回归模型表现出很强的预测性能(SPL的R2 = 0.97;Sa的R2 = 0.85)。基于响应面结果,发现在中等主轴转速、低进给速度和浅铣削深度下可以实现低噪声、高质量加工。这些发现为理解竹材加工中的噪声-粗糙度耦合提供了机制基础,并为计算机数控加工、刀具选择和工业降噪策略提供了实际指导。

1. 引言
竹子是一种生长迅速且可再生的生物质材料,在建筑、家具、室内装饰和复合材料领域具有广泛的应用潜力,这得益于其良好的机械性能、独特的纹理和环境可持续性[1,2,3]。其强度、韧性和可再生性使其成为木材、金属和塑料等传统材料的理想绿色替代品。因此,竹子在这些领域的应用符合社会对资源保护和生态保护日益增长的需求[4,5]。铣削是工业竹材加工中的核心材料去除工艺,其性能直接决定了三个关键方面:加工噪声排放、加工表面质量和加工稳定性。优化铣削工艺以实现低噪声运行和高表面质量是实现竹材高价值工业化的核心前提[6]。
随着竹材利用的工业化程度不断提高,切削过程中的可加工性已成为确保产品质量和生产效率的因素[7,8]。加工引起的噪声污染较为严重,对工作环境和周边社区(包括邻近的工业活动、办公室和住宅区)产生不利影响,这限制了竹材在人口密集地区的工业化应用,并影响了公众对该行业的可持续性认知[9,10]。此外,加工质量往往缺乏稳定性,经常出现表面损伤、毛刺形成和撕裂等缺陷,降低了产品精度和附加值[11]。
现有的关于竹材加工的研究主要集中在切削力、表面粗糙度、切削过程中的功耗和刀具磨损等方面。这些研究考察了刀具几何形状、切削策略和工艺参数对加工性能的影响[12,13,14]。这些研究初步揭示了竹材独特的切削机制,这些机制源于其各向异性的结构、层状纤维微观结构和相对较高的硬度。这些材料特性是竹材加工面临的核心挑战,导致切削方向依赖性、容易的纤维撕裂、分层和毛刺缺陷,以及加剧的切削振动和噪声,从而阻碍了稳定、高质量加工的实现[15]。因此,在向低消耗、高质量制造的转型过程中,当代竹材加工实践仍面临诸多挑战[16]。
当前参数选择的主要局限在于,工业实践大多沿用传统木材加工的参数策略,未能充分考虑竹材的独特材料特性,导致参数严重不匹配。在实际铣削过程中,纤维方向、主轴转速、进给速度和铣削深度并非独立作用,而是对切削噪声和表面质量产生耦合效应。参数不匹配限制了降低峰值声压级(SPL)和改善加工质量的协调优化。除了加工参数外,刀具属性也极大地影响了竹材和木复合材料的铣削噪声和表面粗糙度。螺旋形刀具可以提高稳定性并减少振动噪声;刀片数量较多的刀具可以改善表面光洁度,但可能会增加空气动力噪声;薄刀片可以减少纤维撕裂;分段刀具可以抑制共振噪声[17]。
本研究的目的是解决竹材加工中与表面质量相关的切削噪声问题。这些问题源于参数不匹配,导致加工性能不稳定。因此,本研究关注关键工艺参数对切削噪声和加工质量的影响模式及其潜在机制。参数包括竹纤维方向、主轴转速、进给速度和铣削深度。目标是确定控制竹材铣削过程中SPL与表面粗糙度关系的多目标优化关系。

2. 材料与方法
2.1. 工件和切削工具
本研究使用的试样尺寸为120 mm × 100 mm × 12 mm,由福建品竹竹木公司生产的毛竹(Phyllostachys pubescens)层压而成。由于这些板材具有高强度、高刚度和优异的尺寸稳定性,有效弥补了天然竹材易开裂和变形的缺点,因此被广泛应用于建筑、家具制造和地板制造领域,这些试样能够很好地代表这些应用中常用的材料[15]。试样的主要物理参数列于表1中。为了更深入地研究纤维方向对切削行为的影响,准备了四种纤维方向的试样(0°、45°、90°和135°)。试样的示意图见图1。
2.2. 实验设备
在本研究中,安装了一个声压测量系统(北京盛旺盛电科技有限公司,北京),采样频率为20 kHz,靠近计算机数控(CNC)机床(东莞南兴机械有限公司,中国)(图3)。该系统用于记录加工过程中的声压变化,包括声压传感器(MPA231)、数据采集卡(USB-4431)和数据处理单元。图3显示了铣削噪声分析的实验装置,包括切削配置、声信号采集链和使用超深度场显微镜的表面粗糙度表征。声压传感器的布置遵循GB/T 3767-2016标准[20],传感器位于与工件相同的水平面上,距离声源0.5 m处。采集的声信号根据HJ 706-2014标准进行了进一步校正,以减少背景环境噪声的影响。
声信号在2 s的固定时间窗口内进行处理,对应于稳态切削阶段,排除了刀具进入和退出时的瞬态效应。所有测量均采用A加权处理,符合标准声学评估方法,报告的值以A加权声压级表示。噪声水平主要通过选定时间窗口内的等效连续声压级(Leq)来表征。声信号的采样频率设置为44.1 kHz。
在实验前,对声压测量系统进行了验证,以确保在相同操作条件下的稳定灵敏度和一致响应。在相同切削条件下重复测量显示了良好的重复性,证实了记录的SPL数据的可靠性。
本研究采用峰值声压级(SPL)作为噪声评估指标。SPL能够捕捉到切削过程中最强烈的瞬态声学事件,包括刀具进入时的冲击声、与脆性断裂相关的突然高频排放以及切削力突变引起的峰值压力峰值[10]。SPL的计算公式如下:
(1)
其中p(t)是测量得到的声压,p0 = 20 μPa是参考声压。
表面粗糙度Sa使用超深度场显微镜(VHX-7000,KEYENCE有限公司,上海)进行测量。Sa表示表面的整体微观波动,定义为采样区域内所有测量点相对于参考表面的绝对高度偏差的算术平均值[11]。因此,Sa提供了加工表面整体粗糙度的全面定量表征。
表面粗糙度测量在每个试样的中心区域采用五点采样方法进行。每个样本选择2 mm × 2 mm的五个测量区域,最终Sa值通过这些测量值的平均值获得。观测方向垂直于加工表面,以确保表面采集的一致性。三维表面地形是通过沿垂直方向获取的30张二维图像堆叠,并使用显微镜的深度合成功能自动对焦重建的。
在粗糙度评估之前,对测量到的表面数据进行了平整处理,以消除形变和倾斜。所有Sa值均基于处理后的三维表面数据计算得出。
2.3. 响应面方法
本研究应用了响应面方法(RSM),实验计划采用最优设计[22]进行开发。根据Merhar等人(2026)[23]的研究,RSM广泛用于建模和分析加工系统中工艺参数的影响,为建立输入因素与响应变量之间的经验关系提供了有效框架。与传统方案(如中心复合设计(CCD)和Box–Behnken设计)相比,最优设计在处理不规则因素范围、不可行的因素组合和具有特定约束的实验场景时具有更大的灵活性和效率[24]。在最优设计中,实验点通过数值优化算法选择,主要目的是在预定义的候选点集中最大化可用于模型估计的统计信息。标准评估标准包括D-最优性和I-最优性。
D-最优设计最大化信息矩阵(XTX)的行列式,从而最小化回归系数估计的协方差矩阵(XTX)?1。目标函数如下:
(2)
I-最优设计最小化设计区域内的平均预测方差:
(3)
其中D表示实验域。
本研究采用基于响应面方法(RSM)的实验设计(DOE)来高效研究多个工艺参数及其交互作用的影响。基于预定义的因素范围、候选点集和实际操作约束,生成了最优设计。软件根据指定的最优性标准自动选择了实验运行,从而确保了设计的稳健性和模型估计的准确性。详细的设计结果见表3。同时,根据RSM设计的黄金标准指南[25],每个参数组合进行了五次测试。表3展示了包括因素组合及其对应的铣削噪声和表面粗糙度测量结果。收集的实验数据用于拟合标准二次响应面模型,数学表达式如下:
(4)
其中y是响应变量,xi表示第i个(标准化的)因素,β0、βi、βii和βij分别是截距、线性、二次和交互作用系数,ε是实验误差。
通过分析拟合的响应面和等高线图,可以进一步分析因素之间的交互作用并确定最优参数区域。在这项研究中,定制的最优RSM设计在统计效率、模型准确性和实验成本之间取得了适当的平衡,同时满足了实际工艺约束。3. 结果与讨论 3.1. 切割噪声和加工质量的回归模型 切割噪声峰值声压级(SPL)和加工表面粗糙度(Sa)的评估指标分别由方程(5)和(6)定义。随后基于这些响应变量开发了回归模型。根据拟合统计(表4)和相关性分析(图4),两个模型都展示了令人满意的统计性能,适用于预测分析。如图4所示,数据点的颜色从红色变为蓝色,颜色越接近红色表示数值越高,而颜色越接近蓝色表示数值越低。此外,数据点越接近拟合线,模型的预测性能越好。表4. 切割噪声SPL和机器表面Sa的拟合统计。图4. (a)切割噪声和(b)表面粗糙度的实际值与预测值之间的相关性。SPL模型的拟合优度很高,决定系数(R2)为0.97,调整后的R2为0.93,表明该模型解释了实验数据中的大部分变异。变异系数(CV = 0.77%)和足够的精度值15.42进一步表明模型在实验范围内的行为稳定。同样,Sa模型的R2为0.85,反映了可接受的解释能力。足够的精度值12.08确保了合理的预测性能。如图4所示,预测值与实验结果非常吻合,数据点紧密跟随回归趋势。总体而言,两个模型都满足了在研究参数范围内进行可靠预测所需的统计标准。 (5) (6) 方程(5)和(6)中的回归系数简洁地表明了各因素的影响。在方程(5)中,主轴速度(n)的负系数(-0.97141)表明随着n的增加,SPL呈下降趋势,而相对较大的常数项(143.91740)表明SPL主要受基线切割噪声的影响,工艺参数的影响相对有限。相比之下,与纤维相关的变量(U)的绝对系数较大(-2.67081566),表明其对SPL的影响更为显著。对于方程(6),切割深度(ap)的相对较大系数(-4.20079)表明ap在决定表面粗糙度方面起主导作用。总体而言,系数的符号和大小反映了每个因素的方向和相对重要性,而高阶项则解释了非线性和交互效应。3.2. 切割噪声和加工质量的方差分析 SPL模型的方差分析(ANOVA)结果表明,该模型在95%的置信水平上总体显著(p < 0.0001),并有效解释了响应变量的变化。所有四个主要参数——纤维角度(θ)、主轴速度(n)、进给速度(U)和铣削深度(h)的主效应都非常显著(p < 0.01),证实了它们对噪声的显著单独影响。此外,四个交互项θ * U、n * U、n * h和U * h也显示出显著性,表明这些因素之间存在不可忽视的协同效应。相比之下,交互项θ * n和θ * h以及二次项θ2和h2不显著,表明它们的影响相对较小。拟合的不显著进一步支持了所提出的二次模型的适当性和可靠性。表5. SPL的方差分析。根据Sa的ANOVA结果(表6),建立的模型在统计上显著(p = 0.0003),解释了表面粗糙度Sa总变异的约85%。在主要因素中,铣削深度(h)的影响最为显著,其次是纤维角度(θ)、主轴速度(n)和进给速度(U),所有这些因素在p < 0.05的水平上均显著。此外,纤维角度和主轴速度(θ * U)之间的交互作用也显著,表明它们的组合影响了Sa。没有其他两个因素的交互作用显示出统计显著性。表6. Sa的方差分析。3.3. 铣削参数对切割噪声和加工质量的影响 图5中呈现的实验结果表明,切割噪声SPL与铣削参数之间存在明显的相关性,这些参数与材料结构密切相关。如图5所示,颜色从黄色变为蓝色,蓝色表示较高值,黄色表示较低值。此外,顶视图中的线条代表等高线。图5. (a)主轴速度和纤维角度,(b)进给速度和纤维角度,(c)铣削深度和进给角度,(d)进给速度和主轴速度,(e)铣削深度和主轴速度,(f)铣削深度和进给速度对切割噪声SPL的影响。随着纤维角度的增加,切割噪声水平相应上升。这主要是由于竹子的独特细胞结构表现出明显的各向异性,导致相对材料强度的变化。在0–90°的纤维角度范围内,即沿纤维方向切割时,纤维在铣削过程中被工具压缩并完全切断[26]。随着纤维角度的增加,相对材料强度增加,导致材料断裂时释放的声压强度增大。相反,在90–180°的纤维角度范围内,即逆纤维方向切割时,材料在工具挤压下被剥离[27]。随着纤维角度的增加,纤维细胞分层程度和严重性加剧,导致声压强度增加[11]。随着主轴速度和进给速度的增加,铣削速度相应提高,当工具接触材料时产生的瞬时冲击力增大。这增强了材料失效时产生的声压强度[28];因此,随着主轴速度和进给速度的增加,切割噪声水平也增加[29]。随着铣削深度的增加,每齿的材料去除量增加,导致加工过程中的切割阻力增大。因此,材料断裂时产生的声压强度也增加,从而在铣削深度增加时切割噪声水平提高[30]。总之,在铣削层压竹子时,选择纤维角度较小的材料、降低主轴速度和进给速度以及最小化铣削深度可以有效减少材料去除过程中释放的声压强度,从而降低切割噪声水平[31]。如图6所示,表面粗糙度Sa受到铣削参数的影响,其中铣削深度的影响最为显著。如图6所示,颜色从黄色变为蓝色,蓝色表示较高值,黄色表示较低值。此外,顶视图中的线条代表等高线。图6. (a)主轴速度和纤维角度,(b)进给速度和纤维角度,(c)铣削深度和进给角度,(d)进给速度和主轴速度,(e)铣削深度和主轴速度,(f)铣削深度和进给速度对加工质量(Sa)的影响。随着铣削深度的增加,每齿的材料去除量相应增加,导致切割阻力增大和材料失效时释放的能量增加。这导致加工表面上残留的纤维细胞碎片增多,从而增加了表面粗糙度[32]。随着纤维角度的增加,表面粗糙度增加。在0–90°的纤维角度范围内(沿纤维方向切割时),相对材料强度随着纤维角度的增加而增加,导致材料失效时释放的能量增加。因此,加工表面上残留的纤维细胞碎片增多,从而增加了表面粗糙度[26]。相反,在90–180°的纤维角度范围内(逆纤维方向切割时),材料在工具挤压下被剥离[27]。随着纤维角度的增加,纤维分层程度和严重性加剧,使得加工表面上残留的纤维细胞碎片面积增大,从而增加了表面粗糙度[11]。随着主轴速度的增加,材料去除机制从低速下的塑性变形逐渐转变为高速下的脆性断裂。与长时间的塑性变形过程相比,脆性断裂在竹子细胞中产生更光滑的剪切表面[33];因此,随着主轴速度的增加,表面粗糙度降低[34]。进给速度与表面粗糙度呈正相关。随着进给速度的增加,每齿的材料去除量也增加,使得加工表面上周期性的工具痕迹和沟槽面积增大[35]。同时,材料失效时的切割阻力和能量释放增加,导致表面上残留的纤维细胞碎片增多。因此,随着进给速度的增加,表面粗糙度增加[20]。总之,为了实现竹子铣削的低表面粗糙度,必须严格控制铣削深度和纤维方向,同时协同优化主轴速度和进给速度,特别要注意纤维方向和主轴速度的联合效应。此外,通过响应面和等高线图(图7)进一步可视化了纤维角度和主轴速度之间的交互效应,这些图提供了对这些两个因素联合影响的更深入理解。如图7所示,黄色表示较高的Z轴值,蓝色表示较低的Z轴值;Z轴在图7b中表示切割噪声,在图7d中表示表面粗糙度。如图7a,b所示,切割噪声SPL对主轴速度有轻微的非线性依赖性,同时对纤维角度有略微增加的趋势。响应面与回归模型中确定的主轴速度的显著二次效应一致。图7. 响应面和等高线图显示了纤维角度和主轴速度对(a,b)切割噪声SPL和(c,d)表面粗糙度Sa(U = 1.2 m/min, ap = 0.5 mm)的交互作用。对于表面粗糙度,图7c,d显示Sa随纤维角度呈总体增加趋势,而主轴速度的影响取决于纤维方向。在小纤维角度下,较高的主轴速度倾向于通过促进更平滑的剪切来降低Sa。这种行为与ANOVA结果中报告的纤维角度和主轴速度之间的显著交互作用一致。3.4. 切割噪声和加工质量的优化与验证 优化竹子加工中的切割噪声和表面质量对于推进其作为可持续和高性能材料的工业应用具有重要价值。过高的噪声不仅会恶化工作环境并影响职业和周围条件,还反映了能量散发的低效率和不稳定的切割动态[36]。表面质量差,通常表现为纤维撕裂、毛刺形成和粗糙度增加,会损害竹制品的结构完整性、美观质量和功能性能,从而限制了其在家具、地板和工程复合材料等精密应用中的使用[20]。因此,通过系统化的工艺参数优化实现低噪声操作和优质表面处理之间的适当平衡至关重要[37]。这种优化不仅有助于符合环境和职业标准,还提高了产品的竞争力,并促进了资源高效的制造[38]。根据Natesh等人(2025)[39]的研究,本研究采用了模型验证方法。根据本研究提出的铣削噪声SPL和表面粗糙度Sa的预测模型,确定了不同纤维角度0°、45°和90°对应的铣削参数,包括主轴速度、进给速度和铣削深度,如图8和表7所示。目标是在不同纤维角度条件下分析SPL和Sa的性能,以验证模型的可靠性。如图8所示,红色水平线代表自变量,蓝色斜线代表因变量,数据点的位置表示它们在整体参数范围内的相对位置。图8. 不同纤维角度下通过RSM优化的铣削条件:0°(a)、45°(b)和90°(c)。表7. 优化和验证结果。在这些纤维角度和铣削参数组合下,得到了预测的SPL(声压级)和Sa(表面粗糙度)值:对于0°的纤维角度,预测的SPL和Sa分别为86.83和2.91;对于45°的纤维角度,预测的SPL和Sa分别为86.83和3.79;对于90°的纤维角度,预测的SPL和Sa分别为87.43和4.01。为了验证模型在不同纤维角度条件下的可行性,进行了验证测试,并测量了实际的SPL和Sa值:对于0°的纤维角度,实际的SPL和Sa分别为82.38和3.41;对于45°的纤维角度,实际的SPL和Sa分别为83.46和3.98;对于90°的纤维角度,实际的SPL和Sa分别为92.11和4.52。分别计算了SPL和Sa的误差率:对于0°的纤维角度,SPL和Sa的误差率分别为4.45%和-14.66%;对于45°的纤维角度,误差率分别为3.37%和14.77%;对于90°的纤维角度,误差率分别为-4.68%和-11.28%。所有误差率都在合理范围内。因此,可以得出结论,建立的SPL和Sa预测模型是可靠的。

纤维角度对铣削性能及其相应参数的影响如下:(1) 当纤维角度为0°时,相应的铣削参数为主轴转速9636转/分钟、进给速度2.0米/分钟和铣削深度0.5毫米;在这种条件下,实际的SPL为82.38,实际的Sa为3.41,显示出相对较低的铣削噪声但表面粗糙度表现一般。(2) 当纤维角度为45°时,相应的铣削参数为主轴转速9297转/分钟、进给速度1.4米/分钟和铣削深度0.5毫米;在这种条件下,实际的SPL为83.46,实际的Sa为3.31,在铣削噪声和表面粗糙度方面达到了最佳综合性能。(3) 当纤维角度为90°时,相应的铣削参数为主轴转速9155转/分钟、进给速度1.0米/分钟和铣削深度0.5毫米;在这种条件下,实际的SPL为92.11,实际的Sa为4.52,导致铣削噪声和表面质量最差。

本研究建立的切削噪声和表面质量的双预测模型能够准确预测不同参数组合下的加工性能。通过多目标优化得到的不同纤维方向下的最佳参数匹配规律,可以用于指导CNC加工优化和切削工具设计,适用于地板、家具和层压竹板材等工业应用。本研究的验证基于在优化参数条件下的确认实验。预测值与测量值之间的良好一致性表明,所开发的模型在调查的参数范围内是可靠的。然而,应当注意的是,没有使用独立的验证数据集或交叉验证程序。因此,模型的预测能力主要限于定义的实验领域。

4. 结论

本研究解决了竹材切割过程中的噪声污染和表面质量不稳定问题。系统地研究了关键工艺参数对切削噪声和表面质量的影响,包括0–135°的纤维方向、7000–10,000转/分钟的主轴转速、0.5–2.0米/分钟的进给速度以及0.5–2.0毫米的铣削深度。主要结论如下:
(1) 建立了一个高精度的双响应二次回归预测模型,用于预测切削噪声峰值声压级(SPL)和加工表面粗糙度(Sa)。该模型表现出优异的预测性能,SPL的确定系数R2为0.97,Sa的确定系数R2为0.85,适用于竹材加工中的加工性能预测和参数优化。
(2) 噪声水平随着纤维角度、进给速度和铣削深度的增加而增加,并与主轴转速呈非线性关系。表面粗糙度随着铣削深度和纤维角度的增加而恶化,较高的进给速度也会加剧表面不均匀性。
(3) 纤维方向、主轴转速、进给速度和铣削深度都对切削噪声和表面质量有显著影响。其中,铣削深度是影响表面粗糙度最显著的因素,而纤维方向与主轴转速之间的相互作用在决定表面质量方面起着重要作用。
(4) 通过多目标优化获得了平衡低噪声和高质量加工的综合最佳参数组合:对于45°的纤维方向,最佳参数为主轴转速9297转/分钟、进给速度1.4米/分钟和铣削深度0.5毫米,测得的SPL为83.46分贝,Sa为3.98微米。对于0°的纤维方向,最佳参数为主轴转速9636转/分钟、进给速度2.0米/分钟和铣削深度0.5毫米,实现了最低的SPL为82.38分贝和Sa为3.41微米。对于90°垂直于纤维方向的切割,优化后的参数为主轴转速9155转/分钟、进给速度1.0米/分钟和铣削深度0.5毫米,在三种典型纤维方向中获得了最高的SPL为92.11分贝和最高的Sa为4.52微米。
(5) 未来的研究将集中在多刀头切削器、螺旋切削器、刀具磨损和涂层附着力方面。特别是,未来的研究将把当前的工作扩展到更具代表性的工业刀具配置,如多刃和螺旋铣刀,以提高在实际加工条件下的适用性。将进行包括切削力、温度和切屑变形在内的全面表征,并考虑刀具材料、结构和先进的切削方法。此外,还将探索超出二次响应表面方法(RSM)的更高级建模方法,如高阶非线性或数据驱动方法,以更好地捕捉更广泛参数范围内的复杂切削机制。
相关新闻
生物通微信公众号
微信
新浪微博

热点排行

    今日动态 | 人才市场 | 新技术专栏 | 中国科学人 | 云展台 | BioHot | 云讲堂直播 | 会展中心 | 特价专栏 | 技术快讯 | 免费试用

    版权所有 生物通

    Copyright© eBiotrade.com, All Rights Reserved

    联系信箱:

    粤ICP备09063491号