文章编号：1004-0609(2012)08-2333-07

粉末共注射成形充模流动过程前沿位置及场分布的数值模拟

封  娟¹，何  浩²，李益民²，王光耀²

(1. 中南大学 数学科学与计算技术学院，长沙 410075；

2. 中南大学 粉末冶金研究院，长沙 410083)

摘  要：对粉末共注射充模流动过程进行数值分析和实验验证。采用实验和数值拟合的方法确定芯/壳层界面的厚度，并运用改进的控制体积法对芯、壳层喂料前沿进行追踪。采用有限元和有限差分法对控制方程组进行数值求解，用Matlab进行程序开发，获得芯、壳层充模过程中的熔体前沿分布以及温度场和压力场的分布情况。将模拟结果与实验结果进行对比分析，发现在充填初期，模拟的喂料前沿位置与实验较为吻合，但随着充填的进行，两者偏差增大，其原因可能是在模拟过程中没有考虑注射坯的收缩。

关键词：粉末共注射成形；数值模拟；控制体积法；充模流动

中图分类号：TF124；O242.1          　　     文献标志码：A

Numerical simulation of melt front and filed profile in powder co-injection molding filling process

FENG Juan¹, HE Hao², LI Yi-min², WANG Guang-yao²

(1. School of Mathematical Science and Computing Technology, Central South University, Changsha 410075, China;

2. State Key Laboratory of Powder Metallurgy, Central South University, Changsha 410083, China)

Abstract: The powder co-injection molding filling process was analyzed and verified by the theory of numerical analysis and experiments. The thickness of core/shell interface was determined by experiments and numerical fitting. The melt front was traced by a modified control volume. The pressure equation and energy equation were solved by hybrid finite element/finite difference approaches. The simulation program was developed by means of Matlab, as a result the tracking of movement of skin front and core/skin interface and profile of temperature and pressure fields were obtained. Finally, simulation results were compared with the experimental ones. The simulated results agree well with the movement of feed’s forefront. With the filling of the feed, the deviation between simulation results and the test results increases. The reason may be that the contraction of injected parts is ignored in the process of numerical simulation.

Key words: powder co-injection molding; numerical simulation; control volume method; filling flow

粉末共注射成形是由传统粉末冶金技术与现代塑料共注射技术相结合而产生的金属零部件近净成形新技术^[1]。粉末共注射成形通过不同方式将多种组分在模具内复合注射成形，该技术整合了不同组分的优越性能，具有高性能、低成本、短流程、净成形的工艺优势。在零部件的开发设计及生产方面具有重要的应用价值。因此成为国内外学术界研究的热点。目前，国内外研究机构如英国的Cranfield大学^[2]、德国第三研究院、日本Taisei-Kogyo公司、德国的博世集团以及中国的中南大学粉末冶金国家工程研究中心^[3]等都对其产品开发和过程工艺理论进行了研究。充模过程是粉末共注射成形的关键环节，其模拟研究为优化工艺参数、改进注射模具、预测产品质量以及指导生产具有重要意义。但由于加入了大量粉末，粉末共注射流变行为比塑料共注射要复杂得多，同时喂料的微小变化也会改变其流动行为和产品最终性能，均加大了粉末共注射成形研究和控制的难度。因此，对该过程的研究显得非常必要。

目前，国内在这方面的研究起步较晚，关于共注射的研究大都是塑料共注射成形的充模模拟研究，对粉末共注射成形充模流动的研究基本是空缺的。对粉末共注射充模流动的模拟可以在注塑共注射充模流动模拟的研究上开展。关于注塑充模流动，LI等^[4-5]建立了塑料共注射成形充填过程的物理模型，并给出了数值解法。周国发等^[6]和匡唐清等^[7]将塑料共注射成 形充填过程视为多组分分层流动，采用控制体积法追踪熔体芯、壳层前沿。塑料共注射成形模型推导出的控制方程涉及压力和温度在三维方向上的梯度或散度，场的求解需要解决芯层熔体在充填方向以及厚度方向上界面变化的问题，但目前在充填方向上研究较多，对厚度方向研究却较少。由于共注射成形涉及两种喂料且存在芯/壳层界面，这个界面位置对熔体充填的影响在压力控制方程中有明显体现，仍采用上述忽略厚度方向的方法显然是不合理的。

因此，直接采用传统的塑料共注射成形模拟方法对粉末共注射成形数值模拟问题研究存在欠缺，关键在于芯/壳层界面厚度位置无法确定。本文作者通过实验与分析提出一种确定芯/壳层界面厚度位置的方法，基于改进的控制体积法，对共注射芯、壳层界面进行同时追踪，从而实现粉末共注射成形过程的动态模拟，获得芯、壳层熔体在型腔中压力场、温度场的分布。最后将模拟结果和实验结果进行对比分析。

1  数学模型

由于共注射成形过程以及共注射熔体流变性质的复杂性，本研究采用塑料共注射成型过程的研究方法，将粉末共注射喂料充填过程视为Hele-Shaw流动^[8-9]。在传统单组分注射成形充填过程分析所作假设的基础上，对三维薄壁型腔充填过程进行考虑时，针对共注射具有芯、壳层界面的特点作如下假设：1) 将实际复杂的型腔简化为薄板状型腔，型腔在x和y方向上的厚度尺寸远大于方向上的厚度尺寸，如图1所示。由于z方向上的尺寸较小，可以忽略喂料熔体在充填过在厚度方向上的速率；2) 喂料熔体中粉末与黏结剂混合均匀，喂料在充填过程中不发生两相分离，并且将喂料视为不可压缩的连续介质，即忽略密度的变化，且物性参数不变；3) 喂料在充填过程中，充填方向以对流传热为主，忽略充填方向的热传导，在厚度方向上以热传导为主，忽略厚度方向上的对流传热；4) 忽略惯性力、弹性力和重力喂料在充填过程的影响， 假定压力沿厚度方向不变，壳层前沿熔体的压力为零；5) 假设喂料的流动是层流，且流动是充分发展的，无滑移边界，忽略流动前沿处的“喷泉”效应；6) 忽略喂料充模过程中的凝固现象；7) 忽略物性参数随时间的变化；8) 熔体在芯壳层处有明显的界面存在，在界面处无明显的混合，即混合可以忽略。结合流体力学的基本方程^[10-11]，得到连续性方程：

                     (1)

其中，

                       (2)

压力控制方程：

                             (3)

式中：表示压力梯度。

图1  共注射型腔示意图

Fig. 1  Geometry of plat for co-injection molding (mm)

能量方程：

         (4)

式中：p、T与S代表熔体的压强、温度及流导；δ为芯层厚度的一半；u_i和v_i分别表示在界面处、芯层或者壳层熔体在x和y方向上的速度分量；b为型腔厚度的一半；η_c和η_s分别为芯、壳层的黏度；T为喂料熔体的温度；、c_i、λ_i、η_i和分别为熔体i的密度、比热容、热传导系数、黏度和剪切速率。

根据流动模型的简化及注射成形充填过程中的实际情况，设边界条件如下：1) 壳层流动前沿处，p=0；2)芯/壳层界面处，p_c=p_s，u_c=u_s，v_c=v_s，；3) 沿型腔侧壁，满足不渗透条件，；4) 在型腔入口处，，。其中，Q、Q_i和T₀分别表示流率、浇口处的流率和注射温度。

压力场、温度场的求解交替进行，两者通过喂料熔体黏度模型耦合。考虑熔体为非牛顿流体，其黏度模型采用Cross-WLF黏度模型^[12]：

                       (5)

其中，

                     (6)

                              (7)

                               (8)

式中：为熔体的剪切速率；T为熔体的温度；p为熔体的压强；为零剪切黏度；为熔体由牛顿流体过渡到幂率流体时的临界剪切应力；n为流动指数；D₁为非结晶体在玻璃化温度下的零剪切率黏性系数；D₂为非结晶物在零压力下的玻璃化转换温度；D₃为非结晶物玻璃化转换温度系数；A₁和A₂为决定热膨胀系 数的参数。

2  共注射成形前沿界面和芯/壳层熔体界面位置的确定方法

2.1  芯/壳层界面厚度方向的分析

前面推出的控制方程组，表面上是一个三维方程组，但由于所作出的简化在推导过程中忽略了厚度方向上的速度，从而忽略了在厚度方向上的传导，因此该方程组实际是一个介于二维和三维的方程组。虽然这种2.5维的公式推导在一定程度上简化了控制方程，但实际上并没有忽略三维空间中厚度方向上的位置函数δ的确定。对于传统单组份注射成形不存在这个厚度位置函数δ，因而可以对三维型腔采用三角单元平剖，在厚度方向分成等高的w层，各层之间的流动通过能量交换联系起来，然后采用传统的控制体积法就可以实现前沿界面的追踪。与其相比，共注射成形由

于芯层和壳层组分的不同，在芯/壳层处有一个明显的界面也需要进行追踪，如图2所示。由于芯/壳层界面厚度方向的位置对于充填情况的影响不可忽略，如果采用上述传统研究方法处理共注射成形，就无法确定位置函数δ在任意时刻的值，因此是不合理的。

图2  共注射横剖、纵剖示意图

Fig. 2  Sketch map of cross and vertical section of co-injection molding

为了分析芯/壳层界面厚度随时间的变化情况，采用316L(40%)不锈钢喂料为壳层、316L(60%)不 锈钢喂料为芯层，选择不同的芯层注射时间进行实验，观测芯层熔体厚度的变化情况。图3所示为实验观测到的穿伸位置分别为26.1、37.8、47.0、60.0和76.1 mm处芯层熔体厚度随着注射时间的变化情况。由图3可知，壳层的预填充量是影响芯层熔体厚度方向分布情况的首要因素，并且在某一固定的穿伸位置   处，芯层熔体在厚度方向上的充填时间仅需0.2~0.4 s，与芯层熔体在注射水平方向上的充填时间相比，充填时间较短。通过图3还可以发现，在固定的穿伸位置处，随着充填过程的进行，在开始注射的0.2~0.4 s内芯层熔体在厚度方向上的变化与注射时间大致呈线性关系，在厚度达到最大值后，不再随着注射时间延长发生变化。

图3  不同穿伸位置处共注射充填过程中芯层厚度的变化情况

Fig. 3  Change of thickness of core melt during filling process

因此，基于以上实验研究和分析，对型腔进行平剖在x和y方向上追踪熔体的充填情况时，通过实验先确定芯层熔体在厚度方向上的最大值，再在数值模拟分析时，假定芯、壳层充填量为定值的情况下，厚度方向上芯壳层界面位置是匀速变化的，从而芯壳层界面厚度方向的位置(δ)可以通过式(9)确定

                                      (9)

将δ代入到控制方程中，按照数值求解的方法即可获得任意时刻的芯、壳层熔体前沿的分布情况以及场的分布情况。由此便可确定任意时刻(t)厚度方向上芯/壳层界面的位置。

在壳层充填量一定的情况下，通过实验数据拟合可得芯/壳层界面在厚度方向上的最大值b_max，理论注射时间t_theory=0.3 s，则在数值模拟中假定芯层界面位置δ随时间的变化率为k=b_max/t_theory。

2.2  改进的控制体积法追踪前沿界面 

在共注射成形充填过程中，除了判断控制体积是否被充填满之外，还需要判断控制体积被何种熔体充填，为此对传统的控制体积法^[13-14]做出一些改进。设f为填充因子，其定义为

                                  (10)

式中：为节点i的控制体内已填充喂料熔体的体积；为控制体的体积。根据填充因子f值的大小判断节点的属性，即f=0，空节点(控制体内无流体)；0＜f＜1，前沿节点(控制体未被壳层熔体充满)；f=1，内部节点或入口节点(控制体完全被壳层熔体充满)；0＜f＜2，内部节点(控制体被芯、壳层体充满)；f=2，内部节点或入口节点(控制体完全被芯层熔体充满)。

首先计算所有前沿节点控制体内的体积流率，由于所有前沿节点和空节点压力为零，因此其体积流率(q_i)的计算^[15]如下：

                              (11)

式中：表示对前沿节点i周围的所有节点压力

求和。计算各个前沿节点控制体完全充满所需的时 间，取其最小值作为下一时间步长：

需要指出的是，在注射时先注入一定量的壳层熔体，再将芯壳层两种熔体一起注入。当注射壳层熔体时，最小步长的计算采用式(12a)，当注射芯层熔体时，由于型腔中存在两个前沿，即壳层前沿、芯/壳层界面。因此，最小步长的取值需要取壳层前沿、芯/壳层界面充填的最小值，计算需采用式(12a)和式(12b)。

当一个前沿节点充满而变为内部节点时，与之相邻的所有空节点成为新的前沿节点，同时更新所有其他前沿节点的填充因子：

                          (13)

即可进入下一步压力场的计算，重复这个步骤便可动态地跟踪流动前沿的发展。

这样，在数学模型、边界条件以及熔体前沿界面追踪的问题解决后，获得任意时刻熔体在型腔中的流动情况并进行分析。

3  数值求解

在充模流动过程中，对场的求解以及熔体前沿的追踪作如下处理：对时间步长的选取采用的是使每个时间刚好有一个前沿节点被充满。此时，与其相邻的所有空节点都变成前沿节点。假定流动开始时第一个控制体积已被完全充满，此时可认为熔体处于等温状态，这样便获得初始时刻的前沿位置和温度场。由此可计算出此时的剪切黏度η、流通率s和平均体积流量Q_a，然后迭代求解压力场。在获得压力场数据后计算下一时间间隔和熔体的前沿位置，以便计算下一时刻的温度场。获得温度场数据后，程序返回计算新时刻的压力场，如此循环直至型腔被熔体充满。

首先对流动区域进行剖分，沿型腔中面将其划分成一系列的三角形单元的结构性网格，在厚度方向上由中性面向模壁分成高度相同的w层。采用有限元法对压力场进行求解，设三角形单元内的压力按线性分布，单元向其节点控制体积所提供的流率q_j (j=1, 2, 3)，则

                               (14)

采用Galerkin法^[16]得到单元的局部刚度矩阵

                               (15)

式中：k_ij为线性插值函数矩阵。对每个单元的局部刚度矩阵累加，得到总体刚度矩阵，通过迭代求解确定该时刻的压力场。

求解温度场时需要考虑型腔被充填的具体情况，对于分别被芯、壳层熔体充填的控制体，能量守恒方程(4)中的物性参数、c_i、λ_i和η_i，分别取为芯、壳层熔体对应的物性系数。而对于控制体被芯、壳层熔体充填的区域，即1＜f≤2处，该处控制体遵守的能量守恒方程的物性参数值取芯、壳层熔体物性参数的加权平均值，即

                    (16)

                     (17)

                    (18)

                    (19)

能量方程的物性参数确定后，采用有限元与有限差分混合法对温度场进行求解。在能量方程中热对流项和黏性项按体积作加权平均，其中热对流项采用上风法进行差分，以保证数值计算的稳定性。能量方程的具体差分形式为

        (20)

式中：?z为z坐标变化量；λ为热传导系数；F₁为热对流项差分矩阵；F₂为黏性项差分矩阵。

4  数值模拟的实现与验证

4.1  芯/壳层界面在厚度方向上最大值的拟合

通过2.1的分析可知，为了获得芯/壳层界面δ的值，需要先确定芯/壳层界面在厚度方向上最大值，本文作者选用316L(40%)为壳层、316L(60%)为芯层进行实验，其中壳层充填量为60%。实验和模拟选用的型腔为薄板型腔，具体尺寸如图1所示，实验的工艺参数见表1所示。实验用注射成型机是海天塑机集团有限公司生产的HTSJ160B共注射成型机。研究假设喂料在充填过程中，充填方向以对流传热为主，对流传热系数通过黏度测定获得。由实验获得不同芯层穿伸距离处芯层的厚度，采用Origin软件对实验数据进行6次多项式函数拟合，其形式如式(21)所示。多项式参数值见表2，实验拟合曲线图如图4所示。拟合结果与实验结果的相关系数R=0.992 95。

    (21)

式中：x表示芯层穿伸距离，z表示芯层厚度的最小 值，即b_max。通过式(21)可获得不同芯层穿伸距离处芯层的厚度最大值b_max，然后通过本文作者所提出的芯层熔体在厚度方向上的变化与注射时间呈线性关系的假设，可得任意时刻在任意芯层穿伸距离处芯壳层界面的δ值。

4.2  界面的追踪实现以及温度场、压力场的实现

将4.1节求得的芯壳层界面δ值代入控制方程中，按照数值求解的方法进行编程获得任意时刻芯壳层熔体前沿的分布情况以及场的分布情况。所开发的Matlab程序在惠普公司生产的HPJ400工作站上运  行^[17]。图5所示为共注射成形充填过程，芯、壳层喂料在不同注射时间下的前沿运动分布图，其中，A、B、C、D、E、F为实验曲线，a、b、c、d、e、f为模拟曲线。由图5可以看出，模拟的结果大致反映出芯、壳层喂料充填过程中前沿的运动情况，模拟结果与实验结果吻合较好。对共注射充填过程的模拟，随着充填的进行，模拟结果与实验结果偏差增大。原因是模拟过程中没有考虑注射坯的收缩，造成模拟获得的注射坯与实验获得的注射坯在尺寸上存在差异(见曲线D、E、F和d、e、f)。

表1  模拟和实验工艺参数

Table 1  Injection parameters of simulation and experiment

表2  多项式拟合参数值

Table 2  Fitted values of polynomial

图4  实验拟合曲线图

Fig. 4  Fitting curve of experiment data

图5  芯、壳层熔体在不同芯层注射时间下的前沿分布图

Fig. 5  Front distribution of skin and core melt at various time (A, B, C, separately, is front position of core melt at 0.8, 1.5, 2.0 s through experiment; a, b, c, separately, is front position of core melt at 0.902, 1.653, 2.194 s through simulation; D, E, F, separately, is front position of skin melt at 0.8, 1.5, 2.0 s through experiment; d, e, f, separately, is front position of skin melt at 0.902, 1.653, 2.194 s through simulation)

图5中芯、壳层熔体在不同芯层注射时间下的前沿分布实验曲线A、B、C、D、E、F是通过如下方法获得的：设置芯层注射时间分别为0.8、1.5和2.0 s，在HTSJ160B夹层塑料注射机上进行实验，得到3个不同芯层注射时间的注射坯。由于研究假定熔体在型腔中的流动为层流，对获得的注射坯进行平剖，并沿注射方向将平剖面垂直平分N等份，如图6所示，任取其中相邻的两层，记为第i层和第i+1层，其中第i层更靠近对称轴，然后在坐标纸上相应标记出该层的芯、壳层熔体前沿位置，当所有层都标记完毕时，用曲线连接即可得到此芯层注射时间下的芯、壳层熔体前沿位置，其中Y轴方向只考虑了一半。

通过数值模拟，还可获得芯、壳层熔体在充填过程中压力场、温度场的分布情况。图7和8所示分别为芯、壳层喂料在注射时间为0.902 s时的压力场和温度场分布图。由图7和8可以看出，在充填过程中，壳层熔体在浇口处的压力和温度比较大，随着穿伸位置的增大，壳层熔体的压力和温度随之减小，这与共注射矩形模腔的生产实践相符合。

图6  共注射样品平剖面垂直平分图

Fig. 6  Midnormal drawing of athwartship plane of co-injected part

图7  芯、壳层熔体在时间为0.902 s时的压力场分布图

Fig. 7  Pressure field of core and skin feed stocks at 0.902 s

图8  芯、壳层熔体在时间为0.902 s时的温度场分布图

Fig. 8  Temperature field of core and skin feedstocks at 0.902 s

通过对芯、壳层注射过程芯/壳层界面、壳层前沿追踪的模拟结果与实验结果的对比，可以看出，模拟获得的芯/壳层界面、壳层前沿的结果与实验结果符合程度相对较好。需要指出的是，利用数值模拟技术除了能追踪前沿界面以外，还获得熔体在充填过程中的压力场、温度场等物理场的分布情况，但是，这些场的分布在实验中较难验证，需要通过间接的方法进行验证。例如，可以尝试利用注射坯的裂纹来验证压力场的分布情况，利用熔接线来验证温度场的分布等。这需要在以后的工作中进行进一步的验证。

5  结论

1) 在建立聚合物熔体充填三维薄板型腔的数学模型基础上，对传统的控制体积法进行改进并追踪芯、壳层熔体前沿的运动情况，阐述了粉末共注射成形与普通注射成形充填过程数值模拟的不同，基于芯层在厚度方向上的充填时间较短，和壳层预填充量是影响芯层厚度变化的首要因素的实验结论，提出了在固定的穿伸距离处芯层厚度在未达到最大厚度前，芯层厚度呈线性变化的假设。

2) 在充填过程中利用Malta软件进行程序开发，实现了共注射过程中芯、壳层前沿界面的追踪，以及压力场、温度场的分布的可视化，模拟结果与实验结果吻合很好。

REFERENCES

[1] ALCOCK J R, LOGAN P M, STEPHENSON D J. Surface engineering by co-injection moulding[J]. Surface and Coatings Technology, 1998, 105: 65-71.

[2] 朱 计, 类彦威, 张 杰. 共注射成型技术及其发展[J]. 工程塑料应用, 2007, 35(4): 31-35.
ZHU Ji, LEI Yan-wei, ZHANG Jie. Co-injection moulding and its development[J]. Engineering Plastics Application, 2007, 35(4): 31-35.

[3] ROSOF B H. The metal injection molding process comes of age[J]. Journal of Mater, 1989, 41(8): 13-16.

[4] LI C T, ISAYEV A I. Interface development and encapsulation in simultaneous co-injection molding of disk (Ⅰ): Two- dimensional modeling and formulation[J]. Journal of Applied Polymer Science, 2003, 88: 2300-2309.

[5] LI C T, LEE D J, ISAYEV A I. Interface development and encapsulation in simultaneous co injection molding of disk (Ⅱ): Two-dimensional simulation and experiment[J]. Journal of Applied Polymer Science, 2003, 88: 2310-2318.

[6] 周国发, 孙 懋, 柳和生, 李绅元, 罗忠民. 顺序共注成型的实验研究[J]. 中国塑料, 2003, 17(4): 53-57.
ZHOU Guo-fa, SUN Mao, LIU He-sheng, LI Shen-yuan, LUO Zhong-min. Experimental study of sequential co-injection molding[J]. China Plastics, 2003, 17(4): 53-57.

[7] 匡唐清, 柳和生, 周国发, 赖佳美. 共注射成型充填过程的运动界面追踪[J]. 塑料工业, 2005, 33(3): 32-35.

KUANG Tang-qing, LIU He-sheng, ZHOU Guo-fa, LAI Jia-mei. Tracing of moving interface during polymer mold filling in co-injection molding process[J]. China Plastics Industry, 2005, 33(3): 32-35.

[8] QU Xuan-hui. Numerical simulation of feedstock melt filling in a cylindrical cavity with solidification in powder injection molding[J]. Transactions of Nonferrous Metals Society of China, 1998, 8(4): 544-545.

[9] 曲选辉, 李益民, 黄伯云. 金属粉末注射成形技术[J]. 粉末冶金材料科学与工程, 1996, 1(2): 33-37.

QU Xuan-hui, LI Yi-min, HUANG Bo-yun. The technology of metal powder injection molding[J]. Materials Science and Engineering of Powder Metallurgy, 1996, 1(2): 33-37.

[10] LI C T, LEE D J, ISAYEV A I. Interface and encapsulation in simultaneous co-injection molding of disk: Two-dimensional simulation and experiment[C]//Annual Conference 2002, Society of Plastics Engineers, 2002, 48: 465-469.

[11] HIEBER C A, SHEN S F. A finite-element/finite-difference simulation of the injection-molding filling process[J]. Journal of Non-Newtonian Fluid Mechanica, 1987, 7: 1-32.

[12] KOSZKUL J, NABIALEK J. Viscosity models in simulation of the filling stage of the injection molding process[J]. Journal of Materials Processing Technology, 2004, 157/158: 183-187.

[13] HIRT C W, NICHOLS B D. Volume of fluid (VOF) method for the dynamics of free boundaries[J]. Journal of Computational Physics, 1981, 39: 201-225.

[14] RIDER W J, KOTHE D B. Reconstructing volume tracking[J]. Journal of Computational Physics, 1998, 141: 112-152.

[15] 刘 全, 柳和生, 水鸿寿, 张晓轶. 数值模拟界面流方法进展[J]. 力学进展, 2002, 32(2): 259-274.

LIU Quan, LIU He-sheng, SHUI Hong-shou, ZHANG Xiao-yi. Advances in numerical simulation of interfacial/free-surface flows[J]. Advances in Mechanics, 2002, 32(2): 259-274.

[16] 李庆扬. 数值分析[M]. 武汉: 华中理工大学出版社, 1982.

LI Qing-yang. Numerical analysis[M]. Wuhan: Huazhong University of Science and Technology Press, 1982.

[17] 孙 祥, 徐留美, 吴 青. MATLAB7.0基础教程[M]. 北京: 清华大学出版社, 2005.

SUN Xiang, XU Liu-mei, WU Qing. MATLAB7.0 based tutorial[M]. Beijing: Tsinghua University Press, 2005.

(编辑 龙怀中)

基金项目：国家高技术研究发展计划资助项目(2007AA03Z114)；国家自然科学基金资助项目(51075405)；中南大学博士后基金资助项目(7604130008)

收稿日期：2011-08-01；修订日期：2011-10-26

通信作者：何  浩，助理研究员，博士；电话：13973114617；E-mail: he_hao555@yahoo.com.cn