面向机器人喷涂的多变量涂层厚度分布模型
王国磊 1 , 伊强 1 , 缪东晶 1 , 陈恳 1 , 王力强 2     
1. 清华大学 机械工程系, 北京 100084;
2. 成都飞机工业 (集团) 有限责任公司, 成都 610091
摘要:为了解决传统机器人喷涂模型在喷涂工艺参数改变时会失效的问题,该文将喷涂工艺参数作为模型变量,研究多变量喷涂模型的建模方法。首先,提出了一种基于β分布的涂层生长速率函数,并通过对其进行积分推导出涂层厚度分布方程;其次,通过分析、拟合喷涂实验数据,分别建立喷枪流量、喷涂距离与涂层生长率最大值的关系式,以及喷枪流量、空气压力与喷幅宽度的关系式,并将其代入到涂层厚度分布方程中,建立了以5种常变喷涂工艺参数为自变量的涂层厚度分布泛化模型;最后,通过实验对模型进行验证。结果表明:该模型能够根据工艺参数的变化预测相应的涂层厚度分布,且平均预测误差小于4.3%。
关键词工业机器人    机器人喷涂    涂层生长速率    涂层厚度分布    可变喷涂参数    多变量模型    
Multivariable coating thickness distribution model for robotic spray painting
WANG Guolei1, YI Qiang1, MIAO Dongjing1, CHEN Ken1, WANG Liqiang2     
1. Department of Mechanical Engineering, Tsinghua University, Beijing 100084, China;
2. Chengdu Aircraft Industrial (Group) Co., Ltd., Chengdu 610091, China
Abstract: A multivariable robotic spray painting model was developed for a range of painting parameters to improve the restricted traditional model. A β distribution based coating growth rate function was used with a coating thickness distribution formula then deduced from the integral of the growth rate function. The maximum coating growth rate was related to the paint flow rate and painting distance with the paint flow rate related to the painting air pressure and painting width from experimental data. Then, a generalized coating thickness distribution model was developed with five painting parameters as independent variables by substituting these relations into the coating thickness distribution equations. The model was validated through experiments with the results showing that it can predict the coating thickness distribution for various painting parameters with an average forecasting error of less than 4.3%.
Key words: industrial robot     robotic spray painting     coating growth rate     coating thickness distribution     variable painting parameter     multivariable model    

喷涂模型的建立是喷涂机器人离线编程的基本和关键问题之一,有了准确的模型才能够对涂层厚度进行可靠的预测,从而优化喷涂工艺参数、规划喷涂路径轨迹和控制涂层厚度等。

现有的喷涂模型可以分为无限范围模型和有限范围模型2大类:无限范围模型有Cauchy分布模型[1]和Gauss分布模型[2]等;有限范围模型有椭圆形分布模型[3]、抛物线分布模型[4]β分布模型[5-6]、椭圆双β分布模型[7]、分析沉积模型[8]和组合模型[9]等。但是,这些模型的共性问题是缺乏普适性和泛化能力,即模型仅对特定工艺参数下的涂层厚度分布预测精度很高,一旦喷涂工艺参数改变,涂层厚度分布发生变化时,模型就会失准甚至失效。因为这些模型是基于曲线拟合思想建立的,没有把工艺参数视为变量,模型中不包含工艺参数项,无法反映工艺参数的变化对涂层厚度分布的影响。

因此,本文拟对喷涂工艺参数与涂层厚度分布的关系进行研究,建立一种以喷涂工艺参数为自变量的涂层厚度分布泛化模型。

1 模型自变量确定

空气喷涂是利用压缩空气将涂料雾化成微小颗粒并使之沉积到工件表面形成连续涂层的一种涂装工艺方法。涂层的形成过程非常复杂,涉及到流体流动、雾化、挥发、沉积等很多物理过程,因此影响涂层厚度分布的因素很多,本文将这些影响因素分为4类,如表 1所示。

表 1 影响涂层厚度分布的因素
序号 类别 因素
1 喷涂装置 喷枪、喷嘴、针阀、空气帽、供气装置、供料装置等
2 外部环境 涂料特性、稀释剂比例、涂料温度、工件特性、喷涂环境等
3 喷涂参数 针阀位置、雾化压力、喷幅压力、供料压力、喷枪流量等
4 路径参数 喷涂距离、喷涂速率、喷涂角度等

这4类参数中,喷涂装置和外部环境对于同一台喷涂机器人或者同一批喷涂作业来说一般是不变的,因此可不作为喷枪模型的自变量。

喷涂设定参数中,针阀位置设定好后一般不再改变;雾化压力、喷幅压力和喷枪流量则均属于喷涂作业中的常调节量,需作为模型自变量,而供料压力的变化直接体现在流量上,故无需重复选择。

路径参数中,喷涂距离和喷涂速率都经常根据不同喷涂作业加以调整,因此需要作为喷涂模型的自变量。而喷涂角度则可不必,原因是大多数喷涂作业中喷枪轴线都是严格垂直于工件表面的,而且“与喷涂速度相比,微小的角度改变 ( < 20%) 对喷涂质量的影响是细微的”[10]

因此,最终确定喷枪模型中自变量包括喷涂距离d、喷涂速度v、喷枪流量q、喷幅空气压力pf和雾化空气压力pw

2 实验设计

喷涂实验基本原理见图 1。其中雾化压力、喷幅压力和喷枪流量通过供漆系统调整,喷涂速率和喷涂距离则由喷涂机器人控制。

图 1 喷涂实验

为了进行建模,设计了2种实验:1) 喷幅实验,即喷枪轴线垂直于平面试件且相对于工件保持不动,然后在很短时间内 (0.3~0.5 s) 完成开关枪动作,此实验可获得一个椭圆状的涂层,如图 2a所示,便于分析不同参数对喷幅的影响;2) 直行实验,即喷枪在试件平板上沿直线行程进行喷涂,此实验可以得到一条中间厚两边薄的直线状涂层,如图 2b所示,便于分析不同参数对涂层厚度分布的影响。

图 2 喷涂实验照片

实验所用喷涂机器人为自行研制的THPT-1型多自由度手臂型喷涂机器人[11-12],实验设备与材料有铝制平板试件、FT-150型测厚仪、量杯、电子称重计和锌黄环氧脂漆料等。

3 多变量模型的建立 3.1 基本模型

当喷枪垂直于工件表面并保持静止进行喷涂时,漆雾颗粒附着在工件表面上形成一个椭圆形区域,设其长短轴参数分别是ab。喷涂作业中通常让喷枪沿长轴垂向移动以提高作业效率,因此通常所说的喷幅宽度为2a

椭圆形区域内的涂层的厚度分布呈类椭球状,即喷嘴正下方的涂层最厚,其余各处沿XY轴方向逐渐变薄,并且在X向或Y向断面上涂层的厚度曲线形状是类似的。因此,假设在X向和Y向断面上,涂层生长速率曲线都服从β分布模型,并且在相互平行的断面上β值是相同的。

那么,在y=0的X向断面上,涂层生长率分布满足高度为tmax,开口宽度为2aβ分布为

$t\left( {x, 0} \right) = {t_{\max }}{\left[{1-\frac{{{x^2}}}{{{a^2}}}} \right]^{{\beta _1} -1}}, -a \le x \le a.$ (1)

其中: tmax为椭圆区域内的生长率最大值,即喷枪喷嘴正下方处的生长率,β1β分布系数。

x=kY向断面上,涂层生长率也服从β分布,易知其开口宽度为$2b\sqrt {1-{{\left( {k/a} \right)}^2}} $,再设其高度为tmaxx=k,则该断面上的涂层生长率为

$\begin{array}{c} t\left( {k, y} \right) = t_{\max }^{x = k}{\left[{1-\frac{{{y^2}}}{{{b^2}\left( {1-{k^2}/{a^2}} \right)}}} \right]^{{\beta _2} -1}}, \\ {\rm{s}}{\rm{.t}}{\rm{.}}\\ -a \le k \le a, \\ -b\sqrt {1 - {{\left( {k/a} \right)}^2}} \le y \le b\sqrt {1 - {{\left( {k/a} \right)}^2}} . \end{array}$ (2)

其中β2β分布系数。

事实上,式 (2) 中tmaxx=k就是y=0的X向断面上的x=k点的涂层生长率,可根据式 (1) 求得,因此椭圆区域内任意一点的生长率方程如下所示。

$\begin{array}{c} t\left( {x, y} \right) = \\ {t_{\max }}{\left[{1-\frac{{{x^2}}}{{{a^2}}}} \right]^{{\beta _1} - 1}}{\left[{1-\frac{{{y^2}}}{{{b^2}\left( {1-{x^2}/{a^2}} \right)}}} \right]^{{\beta _2} -1}}, \\ {\rm{s}}{\rm{.t}}{\rm{.}}\\ -a \le x \le a, \\ -b\sqrt {1 - {{\left( {x/a} \right)}^2}} \le y \le b\sqrt {1 - {{\left( {x/a} \right)}^2}} . \end{array}$ (3)

式 (3) 是喷枪垂直于工件表面喷涂时涂层的瞬时生长速率函数,图形化以后如图 3所示。

图 3 涂层生长速率分布

由于机器人喷涂作业时通常使喷枪沿着工件表面匀速喷涂,故需要基于涂层生长速率函数进一步推导喷枪沿工件表面匀速喷涂时的涂层厚度分布。很明显,喷涂区域内一点O(x0, y0) 的涂层厚度是点O被喷炬扫掠过程中涂层不断累积的结果,那么可以认为点O的涂层厚度就是涂层生长速率在点O被喷炬所覆盖的时间段内进行积分的结果,如下所示:

$T\left( {{x_0}, {y_0}} \right) = \int_0^{\frac{{2b\sqrt {1-{{\left( {{x_0}/a} \right)}^2}} }}{v}} {t\left( {{x_0}, {y_0}} \right){\rm{d}}\tau } .$ (4)

其中τ是喷涂时间。变换积分限,得到

$T\left( {{x_0},{y_0}} \right) = 2\int_0^{b\sqrt {1 - {{\left( {{x_0}/a} \right)}^2}} } {\frac{{t\left( {{x_0},y} \right)}}{v}{\rm{d}}y} .$ (5)

为了简化计算,令β1=β2=2后将式 (3) 代入式 (5),得到喷枪沿工件表面匀速喷涂时的涂层剖面厚度分布

$T\left( x \right) = \frac{{4b}}{{3\upsilon }}{t_{\max }}{\left( {1-\frac{{{x^2}}}{{{a^2}}}} \right)^{3/2}}.$ (6)

式 (6) 即是匀速喷涂时涂层厚度分布的基本模型,显然,对于不同的工艺参数,涂层生长率最大值tmax以及喷幅长短轴ab的值是不同的,因此如何将工艺参数作为自变量引入到喷涂模型中的问题就转化为如何将tmaxab表达为喷涂工艺参数。

3.2 雾锥角的引入

由于从喷枪喷出的漆雾呈锥形喷射到工件表面,因此喷幅的长轴和短轴长度ab随喷涂距离d的增大而增大。为了表征这种关系,本文提出了X向雾锥角αxY向雾锥角αy的概念,那么不同喷涂距离下的椭圆长、短轴长度为

$\left\{ \begin{array}{l} a = d \cdot \tan {\alpha _x}, \\ b = d \cdot \tan {\alpha _y}. \end{array} \right.$ (7)
3.3 空气压力的引入

进一步分析,雾锥角又与空气帽结构、针阀开度、喷幅压力和雾化压力有关,在喷涂过程中通常通过喷幅压力和雾化压力进行调节,因此雾锥角是空气压力的因变量。

通过大量喷涂实验发现如下规律:1)Xαxpw增大或pf减小时增大,反之减小;2)Yαypwpf同时增大时增加,反之减小;3)q增大时,两向雾锥角均有所增大。因此本文提出了雾锥角、空气压力、喷枪流量三者之间的关系为

$\left\{ \begin{array}{l} \tan {\alpha _x} = K_\alpha ^x{\left( {\frac{{{p_{\rm{f}}}}}{{{p_{\rm{w}}}}}} \right)^{{a_1}}}{q^{{a_2}}}, \\ \tan {\alpha _y} = K_\alpha ^y{\left( {{p_{\rm{f}}} + {p_{\rm{w}}}} \right)^{{b_1}}}{q^{{b_2}}}, \end{array} \right.{p_{\rm{w}}} > 0, {p_{\rm{f}}} > 0.$ (8)

其中KαxKαy为系数项。

经实验数据拟合与参数优化,可知当a1=1,a2=1/3,b1=1/2,b2=1/3,Kαx=0.133,Kαy=0.069时,式 (8) 可准确预测不同流量、喷幅压力、雾化压力下的X向和Y向雾锥角,如图 45所示。

图 4 X向雾锥角预测曲线与实测数据

图 5 Y向雾锥角预测曲线与实测数据

3.4 喷涂流量的引入

q增加时,单位时间内沉积到工件上的漆雾粒子相应增加,涂层生长速率随之增大。此外,由式 (8) 可知此时喷幅宽度也随之增大,因此假设涂层断面仍然服从相同的β分布。经拟合实验数据,建立qtmax的关系式,如下所示:

$\frac{{{t_{\max }}\left( {{q_1}} \right)}}{{{t_{\max }}\left( {{q_2}} \right)}} = {\left( {\frac{{{q_1}}}{{{q_2}}}} \right)^{1/3}}.$ (9)
3.5 喷涂距离的引入

在雾锥角不变的情况下,喷幅的长短轴长度abd成正比,即

$\frac{{{a_1}}}{{{a_2}}} = \frac{{{b_1}}}{{{b_2}}} = \frac{{{d_1}}}{{{d_2}}}.$ (10)

此外,由于不同喷涂距离下漆雾飞散和稀释剂挥发的程度不同,喷涂距离d还对上漆率 (干燥涂层质量与喷枪喷出漆料质量的比值) 有很大影响,进行上漆率实验可绘制出上漆率ε相对d的变化趋势线,如图 6所示。

图 6 上漆率与喷涂距离关系

那么,可认为d在一定范围内和上漆率ε存在如下关系:

$\varepsilon \left( d \right) = \frac{{{k_\varepsilon }}}{d}.$ (11)

其中kε为上漆率系数。

若在其他喷涂工艺参数相同而仅喷涂距离不同的情况下进行平板直行喷涂,由于喷枪在相同时间内所喷出的涂料流量是相同的,故可以建立等式关系,即

$\frac{{S\left( {{d_1}} \right)vt}}{{\varepsilon \left( {{d_1}} \right)}} = \frac{{S\left( {{d_2}} \right)vt}}{{\varepsilon \left( {{d_2}} \right)}}.$ (12)

其中S(d1) 和S(d2) 分别是喷涂距离d1d2下进行平板直行喷涂所得到涂层的横剖面面积,可计算如下:

$S = 2\int_0^a {T\left( x \right){\rm{d}}x = \frac{{\pi ab}}{{2v}}{t_{\max }}.} $ (13)

联立式 (10)—(13),可推出

$\frac{{{t_{\max }}\left( {{d_1}} \right)}}{{{t_{\max }}\left( {{d_2}} \right)}} = {\left( {\frac{{{d_2}}}{{{d_1}}}} \right)^3}.$ (14)
3.6 喷涂速率的引入

对漆雾粒子的飞行速度来说,喷枪的移动速率是很小的,因此喷涂速率对涂层生长速率几乎没有影响。但是喷涂速率影响工件表面被喷枪喷炬所覆盖的时间,因此喷涂速率越大,涂层生长时间越短;反之,涂层生长时间越长。换个角度来说,喷涂速率越大,喷枪在单位时间内扫掠的区域越大,工件表面的涂层厚度自然也就越薄。

喷涂速率v已存在于式 (4) 的积分上限中,因此无需再行引入。

3.7 泛化模型的建立

联立式 (9) 和 (14),可知椭圆区域内的生长率最大值tmax

${t_{\max }} = K\frac{{{q^{1/3}}}}{{{d^3}}}.$ (15)

其中K为生长率系数。

将式 (7)、(8)、(15) 代入式 (3) 即可得到雾锥覆盖区域内的涂层生长率多变量泛化模型。但该模型表达式过于复杂不便应用,实际应用中可以将式 (7)、(8)、(15) 代入式 (6) 得到喷枪沿工件表面匀速喷涂的涂层剖面厚度

$\begin{array}{c} T\left( x \right) = \\ K'\frac{{{{\left( {{p_{\rm{f}}} + {p_{\rm{w}}}} \right)}^{1/2}}{q^{2/3}}}}{{v{d^2}}}{\left( {1-{{\left( {\frac{{{p_{\rm{w}}}x}}{{K_\alpha ^x{p_{\rm{f}}}{q^{1/3}}d}}} \right)}^2}} \right)^{3/2}}, \\ {\rm{s}}{\rm{.t}}{\rm{.}}\\ -K_\alpha ^x\frac{{{p_{\rm{f}}}}}{{{p_{\rm{w}}}}}{q^{1/3}}d \le x \le K_\alpha ^x\frac{{{p_{\rm{f}}}}}{{{p_{\rm{w}}}}}{q^{1/3}}d. \end{array}$ (16)

其中K′为系数项,可由实验数据拟合获得。

4 参数优化与模型验证 4.1 参数优化

模型确定之后,系数K′的取值决定了模型的预测精确度。为了对K′值进行优化,同时避免过多的实验次数,本文对每种工艺参数在其常用取值范围内按照较小、常用、较大3个原则选取3个不同值,列于表 2之中,然后按照L18(37) 正交实验方法进行平板直行喷涂实验。

表 2 实验参数取值及测量数据点
参数 单位 取值
d cm 15, 17.5, 20
v cm/s 10, 16, 22
q L/min 0.13, 0.16, 0.19
pf MPa 0.22, 0.26, 0.30
pw MPa 0.18, 0.21, 0.24

涂层厚度分布的测量方法是:在平面试件上选取5个垂直于喷枪运动方向的测量截面,在每个截面上分别以1 cm的测量间隔测量涂层厚度,为降低随机因素影响,取5个截面上相同位置测量点的平均膜厚作为该位置的膜厚,形成涂层厚度分布。

实验数据的处理方法是:根据式 (16) 对每组实验的厚度分布进行预测,然后求取测量厚度与预测厚度的误差值E,并绘制K′与E的关系曲线。

$E = \sum\limits_{i = 1}^{18} {\sum\limits_{j = 1}^n {\left| {{T_{ij}}-{{T'}_{ij}}} \right|} } .$ (17)

其中:i表示实验编号;j表示每个喷幅的测量点数目;n的值取决于喷幅宽度,在22~36之间;Tij表示实测厚度,T′ij表示预测厚度。

最终,得到K′与E的关系曲线如图 7所示,可知最优K′值为3 377.00。

图 7 K′与E的关系

4.2 模型验证

为了验证模型的预测精度和泛化能力,取一组不同于表 2的实验参数 (q=145 mL/min, v=14.5 cm/s, d=16 cm, pw=0.2 MPa, pf=0.24 MPa) 进行平板直行喷涂实验,并将测量得到的涂层厚度分布与根据式 (16) 计算所得涂层厚度分布曲线进行比较,实验参数与结果对比见图 8

图 8 预测涂层厚度分布与实测数据对比

实测涂层厚度与模型预测的涂层分布曲线拟合得很好,其中全部21个测量点的测量值与预测值的绝对平均误差仅为0.524 μm,完全可以满足实际工程应用要求。

为了进一步测试模型的泛化能力,根据表 3所示参数,按L8(27) 正交实验方案进行模型验证实验,实验方法同前,然后统计每组实验中涂层厚度测量值与预测值的平均误差,结果列于表 4

表 3 实验参数取值及测量数据点
参数 单位 取值
d cm 16, 18
v cm/s 14.5, 16.5
q L/min 0.145, 0.175
pf MPa 0.24, 0.28
pw MPa 0.2, 0.22

表 4 模型预测结果与实测数据对比
实验序号 幅宽/cm 膜厚绝对平均误差/μm
预测 实测
1 26.8 27.0 0.524
2 28.5 28.0 0.829
3 30.2 30.0 0.709
4 32.0 32.0 0.634
5 37.5 37.0 0.613
6 29.2 29.0 0.722
7 33.3 33.0 0.655
8 25.9 26.0 0.613

所有实验中,模型预测的喷幅宽度与实测喷幅宽度的绝对误差小于0.5 cm,相对误差小于1.8%,膜厚测量值与预测值的绝对平均误差小于0.92 μm,相对平均误差小于4.3%。

通过以上实验,说明本文所提出的多变量涂层厚度分布模型具有良好的泛化能力,可以准确预测不同工艺参数下的涂层厚度分布。在实际应用中,该模型可以为工艺参数的优化和调整提供理论依据,避免不必要的喷涂模型重建,有效减少喷涂实验次数和喷涂机器人离线编程的时间。

5 结论

传统喷枪模型的主要问题是一旦喷涂参数改变,就需要重新建模。本文所建立的模型则以喷涂过程中的可控与常调参数为自变量,因此在喷涂装置、涂料、工件等外部条件不改变的情况下,无需重新实验、拟合就可以获得不同喷枪速度、喷涂距离、喷涂流量、喷幅压力、雾化压力下的涂层厚度分布,较好地解决了传统模型存在的主要问题。

需要指出的是,部分外部因素 (如涂料自身性质、稀释剂含量、涂料的黏度、外界环境温度、工件形状等) 也对涂层厚度分布有很大的影响,如能将这些因素引入到模型中会进一步提高模型的普适性。但是每增加一种可变因素都会使喷涂实验数量呈爆炸性的增长,建模也更加困难、复杂,因此还有待在未来进一步研究和探讨。

参考文献
[1] Antonio J K. Optimal trajectory planning for spray coating[C]//Proc IEEE International Conf on Robotics and Automation. Piscataway, NJ, USA:IEEE Press, 1994:2570-2577.
[2] Freund E, Rokossa D, Rossmann J. Process-oriented approach to an efficient off-line programming of industrial robots[C]//Proc 24th Conf of IEEE Industrial Electronics Society. Piscataway, NJ, USA:IEEE Press, 1998:208-213.
[3] Suh S H, Woo I K, Noh S K. Development of an automatic trajectory planning system (ATPS) for spray painting robots[C]//Proc IEEE International Conf on Robotics and Automation. Piscataway, NJ, USA:IEEE Press, 1991:1948-1955.
[4] Balkan T, Arikan M A S. Surface and process modeling and off-line programming for robotic spray painting of curved surfaces[J]. Journal of Robotic Systems, 2000, 17(9): 479–494. DOI:10.1002/(ISSN)1097-4563
[5] Balkan T, Arikan M A S. Modeling of paint flow rate flux for circular paint sprays by using experimental paint thickness distribution[J]. Mechanics Research Communications, 1999, 26(5): 609–617. DOI:10.1016/S0093-6413(99)00069-5
[6] 陈雁, 颜华, 王力强, 等. 机器人匀速喷涂涂层均匀性分析[J]. 清华大学学报 (自然科学版), 2010, 50(8): 1210–1213. CHEN Yan, YAN Hua, WANG Liqiang, et al. Coating uniformity with a uniform robotic spray gun velocity[J]. J Tsinghua Univ (Sci and Tech), 2010, 50(8): 1210–1213. (in Chinese)
[7] 张永贵, 黄玉美, 高峰, 等. 喷漆机器人空气喷枪的新模型[J]. 机械工程学报, 2006, 42(11): 226–233. ZHANG Yonggui, HUANG Yumei, GAO Feng, et al. New model for air spray gun of robotic spray painting[J]. Chinese Journal of Mechanical Engineering, 2006, 42(11): 226–233. DOI:10.3901/JME.2006.11.226 (in Chinese)
[8] Conner D C, Greenfield A, Atkar P N, et al. Paint deposition modeling for trajectory planning on automotive surfaces[J]. IEEE Transactions on Automation Science and Engineering, 2005, 2(4): 381–392. DOI:10.1109/TASE.2005.851631
[9] 冯川, 孙增圻. 机器人喷涂过程中的喷炬建模及仿真研究[J]. 机器人, 2003, 25(4): 353–358. FENG Chuan, SUN Zengqi. Models of spray gun and simulation in robotics spray painting[J]. Robot, 2003, 25(4): 353–358. (in Chinese)
[10] Pal J F, Jan T G. A real-time algorithm for determining the optimal paint gun orientation in spray paint applications[J]. IEEE Transactions on Automation Science and Engineering, 2010, 7(4): 803–816. DOI:10.1109/TASE.2009.2033567
[11] Chen Y, Shao J Y, Zhang C Q, et al. Redundant-robot-based painting system for variable cross-section S-shape pipe[C]//Proc 2009 ASME/IFToMM International Conf on Reconfigurable Mechanisms and Robots. London, UK:IEEE Press, 2009:743-749.
[12] 陈雁, 邵君奕, 张传清, 等. 复杂管道喷涂系统研制[J]. 机械设计与制造, 2009(11): 1–3. CHEN Yan, SHAO Junyi, ZHANG Chuanqing, et al. Development of spraying system for complex duct[J]. Machinery Design & Manufacturing, 2009(11): 1–3. (in Chinese)