12文章编号：1004-0609(2015)10-2844-10

Al6061板在不同应变路径下粗糙度演变的有限元模拟及实验

徐戊矫，龙  江，章  磊，陈  晨

(重庆大学 材料科学与工程学院，重庆 400044)

摘  要：基于数据统计分析提出一种用于表面粗糙度仿真的有限元分析模型最佳尺寸的确定方法，采用最优模型分析多种应变路径下铝合金Al6061材料的晶粒尺寸和织构类型对表面粗糙度演变的影响。结果表明：不同的应变路径下，塑性变形后的表面粗糙度随着晶粒尺寸的增加而增大；具有较小泰勒因子的织构，在变形过程中自身变形协调性较好，可以获得较好的表面质量；而在相同的变形条件下，具有相同泰勒因子、不同剪切趋势的材料在塑性变形后会呈现出不同的表面形貌，织构之间的剪切趋势差异越大，变形后的表面就越粗糙。

关键词：Al6061板；表面粗糙度；晶体塑性；有限元模型；应变路径；晶粒尺寸；织构种类

中图分类号：TG146      　    　     文献标志码：A

Finite element simulation and experiment of surface roughness evolution for Al6061 sheet under different strain routes

XU Wu-jiao, LONG Jiang, ZHANG Lei, CHEN Chen

(College of Material Science and Engineering, Chongqing University, Chongqing 400044, China)

Abstract: The optimum crystal model size established method for surface roughness simulation was presented based on data statistic analysis. The influences of grain size and textures of Al6061 alloy on surface roughness evolution were discussed under different strain routes. The results indicate that the surface roughness after plastic deformation increases with the increase of the grain size under different strain routes. The textures with small Taylor factor have better deformation coordination, which will be beneficial to reduce the surface roughness. Under the same deformation condition, the models with the same Taylor factor and different the shear tendencies present different surface morphologies. And the larger shear tendency difference is, the worse the surface obtained after plastic deformation will be.

Key words: Al6061 sheet; surface roughness; crystal plastic; finite element model; strain route; grain size; texture type

对于大部分工业金属产品，尤其是对板料制品而言，塑性变形后的表面粗糙化现象会对成品质量产生较大的影响。一般来说，表面粗糙度受许多因素的影响，如初始的表面粗糙度、晶粒尺寸、晶体结构、晶体取向以及织构的三维空间分布等^[1]。BECKER^[2]研究发现，小尺度下的局部应变集中对表面粗糙度也有很重要的影响。因此，影响局部应变的因素，例如应变硬化等也会影响最终的表面粗糙度。ROMANOVA等^[3]通过建立三维模型来分析边界条件，在载荷和微观结构对粗糙度演变的影响，发现通过施加不同的载荷和边界条件，材料织构会很强烈地影响表面的粗糙化过程。WU等^[4]研究发现，表面粗糙度和晶粒取向的空间分布有关，并且对板的表面粗糙度来说，内部晶粒似乎比表层晶粒起着更大的作用。ZHAO等^[5]将各种织构以不同的方式进行组合，分析后发现织构的种类以及空间分布对最终的表面形貌有很大的影响。除此之外，还有很多其他相关文献也讨论表面粗糙化和晶粒取向之间的关系^[6-11]。

目前，国内外学者对晶体结构、晶体取向、载荷等因素对塑性变形后表面粗糙度的影响进行了较多研究，但大多数都是基于某一种变形条件或加载条件来分析的，而在多种变形条件下对各种因素进行系统分析和比较的研究还不多见。又因为板料成形过程中材料受力状况复杂多变，因此，有必要进一步分析和探讨各种因素在不同变形条件下对表面粗糙度演变的 影响。

1  晶体塑性力学模型

1.1  晶体变形的运动学方程

晶体变形时，总的变形梯度F可分解为如下形式^[12]：

                                    (1)

式中：F^*表示由晶格畸变和刚性转动引发的变形梯度。晶体受力后会引发晶格畸变，同时由于变形协调条件的要求和晶界的约束，会导致晶格发生刚性转动，弹性性能可以假定为不受滑移的影响。是晶体沿滑移方向的均匀剪切所产生的变形梯度，的变化率与滑移系以及滑移速率的关系如下：

                      (2)

式中：矢量和分别是参照构型中的滑移方向向量和滑移面的法向方向向量。变形构型中的滑移方向向量和法向向量如下：

                          (3)

当前状态下的速度梯度为

                               (4)

式中：变形张量D和旋转张量可分解为如下形式  (*为弹性部分，p为塑性部分)：

                         (5)

1.2  本构关系及硬化模型

晶格变形张量D^*与柯西应力焦曼率之间的关系可以表示为^[13]

                       (6)

式中：I为二阶单位张量；L为具有一组完全对称特性的弹性模量张量。

对于率相关晶体材料模型而言，滑移系上的滑移速率由其相应的分解剪应力τ^(α)来确定^[14]：

       (7)

式中：sgn(x)是一个符号函数；是滑移系上的参考剪切应变速率；m是应变速率敏感系数；是滑移系上的临界分解剪切应力。

应变强化是通过增量关系用强度的演化来表示：

                            (8)

式中：h_αβ是滑移硬化模量，且在所有开动的滑移系上进行求和。和()分别称为自硬化和潜硬化模量^[14-15]。

                       (9)

式中：h₀是初始硬化模量；τ₀为初始分切应力的临界值；τ_s为饱和分切应力的临界值；γ为所有滑移系上累积剪切应变。

2  微观有限元分析模型的建立

2.1  微观几何建模及边界条件的确定

本文作者在ABAQUS平台下基于晶体塑性力学^[12-16]和广义平面应变模型在不同应变路径下分析了晶粒尺寸和织构对塑性变形后表面粗糙度演变的影响。广义平面应变模型是平面应变模型的拓展，它允许二维平面模型单元变形时拥有垂直于面的应变(ε_zz)，当长轴方向的长度不是无限长时，这种假设的应变状态更符合实际的变形情况。同时，也可以采用二维单元更高效地模拟三维变形。

多晶体几何模型的建立方式包括扩展Taylor模型、用规则多边形或多面体代替晶粒以及结合Voronoi图的结构与算法来建立等^[17]。研究者通常以正六边形或者以每个有限元单元^[18]来建立多晶体CAE分析模型。图1所示为有限元分析模型示意图，其中图1(a)所示为板材截面示意图，图1(b)所示为一个随机建立的二维Voronoi模型用于表示多晶体材料的微观结构^[19]，其中不同颜色代表不同的晶粒。图1(b)所示的二维微观几何模型就代表了图1(a) 中所示的二维截面，如ND-RD截面。本研究中所用的材料微观参数如表1所列^[20]。

表1 材料参数

Table 1  Material parameters

图1  有限元分析模型示意图

Fig. 1  Schematic diagram of finite element model

有限元分析时，位移边界条件可以使用沿x方向施加的应变ε_xx来概括，令模型z向和x向的应变比为ρ=ε_zz/ε_xx。边界条件的设定如下所示：

         (10)

式中：下标“0”表示初始的位置或尺寸；U表示相应方向上的位移，Y₀=H₀面为模型的自由面(如图1(b)所示)。其中，不同的ρ值对应于不同的应变状态(见表2)。其中垂直于平面的应变速率等于 ，垂直于平面的应变可以通过ABAQUS文档中的广义平面模型均匀的施加在平面上，这样在所有的应变路径下都可以使Y₀=H₀保持自由状态。计算完成后的表面粗糙度用轮廓算术平均偏差R_a表示。当取样长度为L，表面轮廓到基准面的距离为δ时，R_a可表示为

                                  (11)

其离散形式为

                               (12)

表2  应变比ρ与应变路径的对应关系

Table 2  Relationship between strain ratio and strain route

2.2  有限元分析模型尺寸确定

有限元模型的尺寸与计算效率和计算精度息息相关。对于二维模型，表面粗糙度R_a的计算需要一个合适的取样长度，而最小取样长度是计算效率和精度最好的平衡点。因此，最优模型的本质就是找到二维模型表面粗糙度R_a的最小取样长度。利用最优模型不仅可以在保证计算精度的前提下提高计算效率，而且能以统计学的角度来代表材料微观组织结构的基本特 征^[21]。本文作者采用统计计算的方法来确定最优模型的长度，在此以平均晶粒尺寸70 μm的模型为例进行阐述。

由于实验中Al6061板料厚度为1.5 mm，因此，模型的高度统一设定为0.75 mm并作对称处理；而模型长度(图1(b)中的L)则从最小的0.5 mm变化到8 mm。每一个长度都对应着一个模型。每一个模型都会被拉伸到应变0.2后计算其表面粗糙度R_a，每个模型尺寸重复计算40次，每次随机赋予材料织构，最终得到的粗糙度平均值和分布标准差如图2所示(黑色方块为每次计算得到的粗糙度值)。从图2中可以看到，表面粗糙度和标准差的变化曲线几乎都在2.25 mm长度下达到稳定。因此，2.25 mm就是最小取样长度，即是在保证得到准确R_a前提下的最小模型长度。那么在平均晶粒尺寸70 μm下，最优模型的长度就可以确定为2.25 mm。

图2  模型平均表面粗糙度R_a和标准差随模型长度的变化(晶粒尺寸为70 μm)

Fig. 2  Variance of model surface roughness(R_a) and st-deviation with model length (average grain size 70 μm)

图3所示为不同晶粒尺寸下粗糙度平均值随模型尺寸的变化关系图。从图3中可以看出，单向拉伸变形后(应变为0.2)的表面粗糙度R_a和最优模型长度L与晶粒尺寸d基本成线性关系的，拟合后其表达式如下：

R_a=0.15d-1.33                                (13)

L=22d+0.774                                  (14)

因此，可以通过式(13)和(14)来近似预测变形后的表面粗糙度值以及确定分析时所需的模型尺寸大小，这为后面在不同晶粒尺度下的模拟分析奠定了基础。

图3  不同晶粒尺寸下模型长度L、晶粒尺寸d和平均表面粗糙度R_a之间的关系

Fig. 3  Relationships among model length(L), grain size(d) and average surface roughness(R_a)

3  不同应变路径下的模拟分析

3.1  晶粒尺寸的影响

以单向拉伸为例进行模拟分析。基于2.2节所确定的有限元分析模型尺寸建立晶粒尺寸为45、70、100和150 μm的模型并沿RD方向拉伸，其初始粗糙度为0，拉伸应变为0.1。为了减小晶粒取向对表面粗糙度的影响，每种晶粒尺寸下进行多次分析，每次分析都重新随机的赋予晶粒取向，将所得结果求平均值后绘出如图4所示的表面粗糙度R_a随应变的变化曲线。

图4  不同晶粒尺度在单向拉伸状态下的表面粗糙度和应变的关系

Fig. 4  Relationship between surface roughness and stretching strain for different grain-scales under uniaxial tension

从图4中可以看出，在相同应变的条件下，表面粗糙度随着晶粒直径的增大而增大。晶粒尺寸为100 μm和150 μm时的粗糙度演变曲线具有较大的斜率，说明大晶粒尺寸下的表面粗糙度的增加速度明显要高于小晶粒尺寸模型的。造成这一结果的可能原因是晶粒尺寸越细小，晶粒内部和晶界附近的应变差值越小，晶粒的变形就越均匀，因此，可以通过细化晶粒来减小变形后的表面粗糙度。

通过相同的方法可以得到其他应变状态下的表面粗糙度R_a随应变的变化曲线，如图5所示。从图5中可以看出，在这几种应变状态下，其粗糙度随应变的变化趋势和单向拉伸下是大致相同的，即是在相同的应变下，表面粗糙度和表面粗化速度都随着晶粒尺寸的增加而增大。因此，细化晶粒可以有效地提高塑性变形后的表面质量。

图5  不同应变路径下晶粒尺寸对表面粗糙度的影响

Fig. 5  Effect of grain size on surface roughness under different strain routes

3.2  织构类型的影响

对一般的冷塑性制品来说，其材料都是具有一定初始织构的。铝合金板材中常见的织构有立方(Cube)织构(0°, 0°, 0°)、黄铜(Brass)织构(35°, 45°, 0°)、戈斯(Goss)织构(0°, 45°, 0°)、S织构(60°, 32°, 65°)等^[22]。首先以单向拉伸为例来分析各种织构对塑性成形后表面粗糙度的影响。取晶粒尺寸45 μm的模型赋予不同的织构类型，然后沿RD方向进行拉伸。对每种织构类型模拟多次，取平均值后得到如图6所示结果。

从图6中可以看出，S织构和随机织构变形后表现出较大的表面粗糙度，立方织构和戈斯织构变形后表面粗糙度变化很小。由此可见，在相同的变形条件下各种织构对变形的协调能力并不相同，此种能力在一定程度上可以由泰勒因子来反映^[5]。泰勒因子越小，织构的变形越容易，因此，对变形具有良好的协调能力。各种织构沿RD方向拉伸时的泰勒因子如表3所示^{[5, 23-24]}。立方织构和戈斯织构的泰勒因子较小且相同，属于较“软”的织构，具有很好的变形协调性，当材料中具有较强的此种织构时，材料整体上变形较为均匀，局部应变集中情况较少，因此，会得到较小的表面粗糙度。

图6  单向拉伸下不同织构种类表面粗糙度与应变的关系

Fig. 6  Relationship between surface roughness and stretching strain of different texture types under uniaxial tension

表3  不同织构在RD方向上对应的泰勒因子

Table 3  Taylor factor for different textures along RD

从随机织构和S织构的粗糙度演变曲线可以看出，仅通过泰勒因子大小并不能完全解释表面粗糙度演变的规律。除泰勒因子外，不同织构剪切趋势的差异也会对变形过程中的表面质量产生影响^{[5, 25]}。这种由剪切引起的材料异常流动通常在黄铜织构、铜型织构和S织构中比较常见。由于在面心立方金属中，S织构往往会表现出较明显的剪切趋势，因此，以S织构为例探讨剪切趋势对表面形貌变化的影响。不同应变路径下S织构的U2向(即y向)位移分布云图如图7所示。由图7可知，同一种织构在不同应变状态和边界条件下会呈现出不同程度的剪切趋势，因此，对表面粗糙度的影响程度也各有不同。由图7(e)可看出，在平面应变压缩下材料发生了严重的局部失稳，因此，没有统计S织构变形后的粗糙度。

图8所示为其余不同应变路径下织构种类对表面粗糙度的影响。从图8可以看出，在不同应变条件下，各种织构的变化趋势和单向拉伸状态下是大致相似的。然而值得注意的是，它们对表面造成的粗糙化程度并不是与泰勒因子的大小严格保持一致，这再次说明泰勒因子并不是影响表面粗糙化的唯一主要因素。而导致此种现象的主要原因可能和织构剪切趋势之间的差异有关^{[2, 5]}。

由于晶体对称性，S织构通常具有几种等效形式，如S₁(60°, 32°, 65°)和S₂(300°, 148°, 245°)，它们具有相同的泰勒因子，但剪切趋势正好相反^[5]。为了能更清晰地对比这两种取向的表面质量，以平面纯剪为例进行分析(表面端部有翘曲，易于观察)。S₁和S₂模型的分析结果如图9所示，在相同的应变条件下，S₁和S₂在应变云图和应力云图上存在较大的差异：S₁具有较明显的应力和应变集中；而S₂则表现得相对均匀。此外，S₁造成了表面末端的翘曲，而S₂的表面整体较平整。

此外，织构之间剪切趋势的差异也会对表面质量产生影响。将不同织构进行组合然后模拟分析，其结果如图10所示。由于沿着RD拉伸会减弱剪切趋势的影响^[5]，为了更好地分析剪切趋势与表面质量的关系，图10的分析都是在单向拉伸下沿着TD方向进行的。图10所示为对单一的S₁和S₂织构以及S₁和S₂的混合织构拉伸后的位移云图。图11所示为3种织构组合变形后的表面轮廓图。图12所示为3种模型所得粗糙度随应变的变化曲线。

通过图10~12可以看出，单独的S₁和S₂织构经过拉伸后的表面起伏较小，而二者结合之后的表面轮廓起伏严重，但其粗糙化速度随应变的增加大致是保持一致的。因此，基于以上的分析可知，不同织构剪切趋势的差异会使表面粗糙化更加严重。

图7  不同应变路径下S型织构的U2分布云图

Fig. 7  U2 displacement contours of S texture under different strain routes

图8  不同应变路径下织构种类对表面粗糙度的影响

Fig. 8  Effect of texture type on surface roughening under different strain routes

图9  在平面纯剪下织构S₁和S₂所对应的应变和应力云图

Fig. 9  Strain((a), (b)) and stress((c), (d)) contours for S₁((a), (c)) and S₂((b), (d)) under in-plane pure shear

图10  不同织构组合下的U2位移云图

Fig. 10  U2 displacement contours of different textures

图11  不同织构组合下的U2变形形貌

Fig. 11  U2 surface profile of different textures

图12  不同织构组合对应的粗糙度与应变的关系

Fig.12  Relationships between surface roughness and stretching strain of different texture combinations along TD

4  单向拉伸实验验证

为了和模拟结果进行对比，根据《GBT16865—1997变形铝、镁及其合金加工制品拉伸试验用试样》标准设计了Al6061的单向拉伸试样，板料厚度为1.5 mm，拉伸试样尺寸示意图及照片如图13所示。单向拉伸实验在微机控制电子万能材料试验机CMT5305上进行，粗糙度采用TR200型手持式粗糙度仪测量，实验分为2组，即平均晶粒尺寸约为45和150 μm的试样，分别命名为试样1和2。其金相照片和相应的取向信息见图14和15。

图13  拉伸试样尺寸示意图及照片

Fig. 13  Schematic diagram of tensile sample(a) and photo(b) configuration of samples (Unit: mm)

图14  不同试样的金相照片

Fig. 14  Metallographs of samples

图15  试样{111}, {200}, {220}极图

Fig. 15  {111}, {200}, {220} pole figures of samples

在实验过程中，每拉伸原始标距的2%进行一次表面粗糙度的测量(沿RD方向测量5次并取平均值)，然后重新加载拉伸，直到拉伸变形量为原始标距的10%为止。同时，两种晶粒尺寸对应的最优模型在被赋予相应的取向信息后进行单向拉伸模拟。最终对比结果如图16所示。

从图16可以看到，模拟与实验得到的表面粗糙度随拉伸变形量的变化趋势基本一致，考虑到二维模型的局限性以及测量误差等原因，实验与模拟结果吻合良好。

对比结果表明，基于晶体塑性有限元模型进行表面粗糙度的模拟分析是可靠的，模拟结果是可接受的。此外，通过统计计算得到最优模型长度的方法也是可行的。

图16  实验结果与模拟结果对比

Fig.16  Comparison of experiment and simulation results

5  结论

1) 基于数据统计分析提出一种用于表面粗糙度仿真的有限元分析模型最佳尺寸的确定方法，为模型尺寸的合理选择奠定了基础，避免了尺寸选择的盲目性，此种思路也适用于其他的同类问题。

2) 在不同应变路径下，晶粒尺寸对粗糙度演变的影响规律相似，即晶粒尺寸越大表面越粗糙。因此，可以通过细化晶粒来控制塑性变形过程中的表面粗化现象。

3) 不同路径下，具有较小泰勒因子的织构可以得到较好的表面质量，剪切趋势之间的差异会使表面粗糙化更加严重。

REFERENCES

[1] RAABE D, SACHTLEBER M, WEILAND H, SCHEELE G, ZHAO Z S. Grain-scale micromechanics of polycrystal surfaces during plastic straining[J]. Acta Materialia, 2007, 51: 1539-1560.

[2] BECKER R. Effects of strain localization of surface roughness during sheet forming[J]. Acta Materialia, 1998, 46: 1385-1401.

[3] ROMANOVA V A, BALOKHONOV R R, SCHMAUDER S. Numerical study of mesoscale surface roughening in aluminum polycrystals under tension[J]. Materials Science and Engineering A, 2013, 564: 255-263.

[4] WU P D, LLOYD D J. Analysis of surface roughening in AA6111 automotive sheet[J]. Acta Materialia, 2004, 52: 1785-1798.

[5] ZHAO Z, RADOVITZKY R. A study of surface roughening in fcc metals using direct numerical simulation[J]. Acta Materialia, 2004, 52: 5791-5804.

[6] GROCHE P, SCHAFER R, JUSTINGER H, LUDWIG M. On the correlation between crystallographic grain size and surface evolution in metal forming processes[J]. International Journal of Mechanical Sciences, 2010, 52: 523-530.

[7] STOUDT M R, HUBBARD J B, LEIGH S D. On the relationship between deformation-induced surface roughness and plastic strain in AA5052 — Is it really linear[J]. Metallurgical and Materials Transactions A, 2011, 42: 2668-2679.

[8] SHIN H, AN J, PARK S, LEE D. The effect of texture on ridging of ferritic stainless steel[J]. Acta Materialia, 2003, 51: 4693-706.

[9] STOUDT M R, LEVINE L E, CREUZIGER A, HUBBARDA J B. The fundamental relationships between grain orientation, deformation-induced surface roughness and strain localization in an aluminum alloy[J]. Materials Science and Engineering A, 2011, 530: 107-116.

[10] LI H J, JIANG Z Y, WEIA D B, et al. Study on surface asperity flattening during uniaxial planar compression[J]. Wear, 2011, 271: 1778-1784.

[11] KALIDINDI S R, BRONKHORST C A, ANAND L. Crystallographic texture evolution in bulk deformation processing of FCC metals[J]. Journal of the Mechanics and Physics of Solids, 1992, 40: 552-555.

[12] ASARO R J, RICE J R. Strain localization in ductile single crystals[J]. Journal of the Mechanics and Physics of Solids, 1977, 25: 309-338.

[13] Hill R, Rice J R. Constitutive analysis of elastic-Plastic crystals at arbitrary strain[J]. Journal of the Mechanics and Physics of Solids, 1972, 20: 401-413.

[14] PIEREE D, ASARO R J, Needleman A. Material rate dependence and localized deformation in crystalline solids[J]. Acta Materialia, 1983, 31(12): 1951-1976.

[15] PIEREE D, ASARO R J, NEEDLEMAN A. An analysis of nonuniform and localized deformation in ductile single crystals[J]. Acta Materialia, 1982, 30: 1087-1119.

[16] HILL R, HAVNER K. Perspectives in the mechanics of elastoplastic crystals[J]. Journal of the Mechanics and Physics of Solids, 1982, 30: 5-22.

[17] 汪 凯.多晶体材料加工的细观塑性有限元模拟[D]. 昆明: 昆明理工大学, 2013: 25-30.

WANG Kai. The simulation of materials processing with crystal plasticity finite element method[D]. Kunming: Kunming University of Science and Technology, 2013: 25-30.

[18] 皮华春, 韩静涛, 章传国. 面心纯铝轧制织构的晶体塑性有限元模拟[J]. 北京科技大学学报, 2007, 29(9): 920-924.

PI Hua-chun, HAN Jing-tao, ZHANG Chuan-guo. Simulation of the rolling texture of face-centered pure Al using crystal plasticity finite element method[J]. Journal of University of Science and Technology Beijing, 2007, 29(9): 920-924.

[19] AURENHAMMER F. Voronoi diagrams: A survey of a fundamental geometric data structure[J]. Acm Computing Surveys, 1991, 23: 345-405.

[20] JUNG K H, K D K, IM Y T, LEE Y S. Crystal plasticity finite element analysis of texture evolution during rolling of fcc polycrystalline metal[J]. Materials Transactions, 2013, 54(5): 769-775.

[21] 章继峰, 王振清, 周建生, 张博明. 基于Python-Abaqus复合材料代表性体积元的数值模型[J]. 宇航材料工艺, 2009(3): 25-29.

ZHANG Ji-feng, WANG Zhen-qing, ZHOU Jian-sheng, ZHANG Bo-ming. Numerical modeling of composite representative volume element based on Python-Abaqus[J]. Aerospace Materials & Technology, 2009(3): 25-29.

[22] 毛卫民, 张新明. 晶体材料织构定量分析[M]. 北京: 冶金工业出版社, 1993: 65-80.

MAO Wei-min, ZHANG Xin-ming. Quantitative analysis of crystal material texture[M]. Beijing: Metallurgical Industry Press, 1993: 65-80.

[23] RAABE D, SACHTLEBER M, WEILAND H, SCHEELE G, ZHAO Z. Grain-scale micromechanics of polycrystal surfaces during plastic straining[J]. Acta Materialia, 2003, 51: 1539-1560.

[24] STARINK M J, WANG S C. A model for the yield strength of overaged Al-Zn-Mg-Cu alloys[J]. Acta Materialia, 2003, 51: 5131-5150.

[25] 刘玉振. 基于ABAQUS的材料微结构有限元计算的后处理研究[D]. 兰州: 兰州理工大学, 2008: 6-10.

LIU Yu-zhen. Finite element’s post-processing of microstructure for polycrystalline materials owing to the ABAQUS material studies[D]. Lanzhou: Lanzhou University of Technology, 2008: 6-10.

(编辑  李艳红)

基金项目：国家自然科学基金青年基金资助项目(51205427)；重庆市教委科学技术研究项目(KJ1401320)

收稿日期：2015-02-09；修订日期：2015-07-10

通信作者：徐戊矫，教授，博士；电话：13368186252；E-mail: xuwujiao_cq@163.com