Applications of random sample consistency algorithm on laser spectroscopy
-
摘要: 为了解决波长调制激光光谱技术探测大气痕量气体浓度中信号处理算法的不足，提出了一种基于随机抽样一致性算法的气体浓度反演算法。以大气甲醛分子的仿真信号和实际测量信号为例，进行了理论分析和实验研究，并与传统的最小二乘法相比较。结果表明，该算法具有较强的抗噪声和异常点干扰能力，尤其是在低信噪比的条件下，精确度可提高1个量级，体现出较高的可靠性和优越性。Abstract: In order to solve the insufficient of signal process algorithms during the detection of atmospheric trace gas concentration by wavelength modulation laser spectroscopy technique, a new method of gas concentration inversion based on the random sample consistency (RANSAC) algorithm was proposed. By choosing the simulation signal and the actual measurement signal of formaldehyde in the atmosphere as examples, theoretical analysis and experimental study were carried out and compared with the traditional least square method. The results show that the proposed algorithm has better immunity to noises and outliers. Especially under the conditions of low signal-to-noise ratio (SNR), the measurement accuracy can be improved by one order of magnitude. The algorithm shows better reliability and superiority.
-
引言
激光打孔过程既存在熔化又存在气化蒸发，是一个复杂的多态多物理场耦合过程。国内外研究人员围绕激光参量对激光打孔质量的影响规律与激光打孔工艺参量优化已经开展了大量的实验研究[1-9]，如WANG[1]等人基于单因素法的激光打孔实验，分析了激光器电压、脉冲宽度、重复频率、聚焦条件、辅助气体等对不锈钢激光打孔的影响规律。FU[2]等人利用光纤激光打孔实验分析了激光功率、占空比、切割速率、重复频率等参量对光纤激光打孔质量的影响。QIAN[3]等人基于正交实验方法进行的激光打孔实验，得出了SUS304不锈钢激光打孔的最优实验参量组合。
但是，采用传统的实验手段难以解释和分析激光打孔的机理以及激光打孔过程中材料的相变过程。采用模拟仿真的方法有利于研究激光打孔的物理机制和瞬态过程，并为实际的激光打孔加工提供理论指导。
早期针对激光打孔过程的模拟仿真主要基于解析法[10-11]，如SWIFT-HOOK等人引入高斯分布的激光热源建立了激光热作用物理模型，分析了激光移动速率对熔池大小及形状的影响[10]，但未考虑相变潜热的影响。WU等人以无限大均匀介质中热扩散方程的解为理论工具，给出了激光打孔过程的近似数学表达式，模拟了激光打孔的孔洞形成过程[11]。但利用解析法求解时，无法处理材料热物理参量随温度变化的实际情况[12]。
目前，针对激光打孔过程的数值模拟主要基于热传导理论[13-17]，即通过求解热传导方程得到激光打孔过程的温度场。ZHANG等人建立了1维稳态热传导模型，研究了长脉冲激光打孔的效能比，计算了不同激光能量密度下的打孔深度[13]。CHU等人基于有限元法建立了2维激光打孔模型，考虑了激光束空间分布和材料相变潜热对孔的影响，计算了激光打孔过程温度场和孔型演化过程[14]。BEGIC-HAJDAREVIC等人基于有限体积法计算了激光打孔过程中温度的分布和孔的几何形状[15]。SONG等人利用ANSYS软件中的单元生死技术对激光打孔过程的温度场进行模拟仿真，得到了小孔的孔深、孔径的时间特性以及随激光能量的变化曲线[16]，但该模型忽略了打孔过程中产生的熔化现象, 只考虑材料的气化。
基于热传导理论的模型大多忽略了重力、反冲压力和粘滞力等对激光打孔的影响，且仅考虑材料的固-液相变过程，忽略了材料的液-气相变过程，因此仅能模拟激光打孔过程中的温度分布情况。另外，激光打孔过程中材料发生相变会导致材料的热物理属性发生阶跃跳变，需要对随温度变化的材料参量进行平滑处理。材料的熔化和气化蒸发过程还涉及到粘滞力和反冲压力的作用，以及相变潜热和热对流的影响，因此基于热传导理论建立准确的激光打孔模型较为困难。
本文中以航空涡轮发动机叶片常用的GH4037镍基高温合金为研究对象，基于流体传热和流体力学理论建立了控制方程组，模型中考虑了重力、反冲压力和粘滞力的作用，也考虑了材料的固/液相变和气/液相变过程，提高了激光打孔模型的准确性。通过计算得到了激光打孔相变过程中不同时刻的温度场分布、熔池流动速率以及气化蒸发速率，为进一步开展激光打孔的研究奠定了理论基础。
1. 数学模型
1.1 有限元模型
本文中有限元模型的建立、参量设置、边界条件设置、求解以及后处理等均基于COMSOL MULTIPHY-SICS多物理场仿真软件。为便于计算，有限元模型采用2维模型，分为上下两个矩形域。上矩形域为空气，长4mm，高4mm；下矩形域为GH4037镍基高温合金，长4mm，高1mm。考虑到模型精确性，同时减少计算量，在空气以及材料表面以下0.5mm左右的深度构造较细的网格，其它区域构建较粗的网格，网格类型采用三角形网格。几何模型的网格划分效果如图 1所示。
1.2 控制方程
假设激光打孔过程中主要受到重力、粘滞力和反冲压力的作用，熔化和气化材料的传热过程主要受到热传导和热对流的作用。根据质量守恒、动量守恒和能量守恒原则，给出气、液、固三相统一的控制方程如下：
\frac{{\partial \rho }}{{\partial t}} + \nabla \left( {\rho \mathit{\boldsymbol{u}}} \right) = 0 (1) \begin{array}{l} \rho \frac{{\partial \mathit{\boldsymbol{u}}}}{{\partial t}} + \rho \left( {\mathit{\boldsymbol{u}} \cdot \nabla } \right)\mathit{\boldsymbol{u = }}\nabla \cdot \left[ { - pI + } \right.\\ \left. {\mu \left( {\nabla \mathit{\boldsymbol{u}} + {{\left( {\nabla \mathit{\boldsymbol{u}}} \right)}^\prime }} \right) - \frac{2}{3}\mu \left( {\nabla \cdot \mathit{\boldsymbol{u}}} \right)I} \right] + \rho \cdot g \end{array} (2) \rho {c_p}\frac{{\partial T}}{{\partial t}} + \rho {c_p}\mathit{\boldsymbol{u}}\nabla T = \nabla \cdot \left( {k\nabla T} \right) + S(x, t) (3) 式中，ρ为密度，u为速度矢量，p为压强，μ为动态粘度，g为重力加速度，T为绝对温度，cp为比定压热容，κ为热导率，S(x, t)为热源项。
激光热源采用高斯面热源模型，热源项的数学表达式如下:
S(x, t) = \frac{{2P}}{{{\rm{ \mathsf{ π} }}{r^2}}}\left( {1 - R} \right)\exp \left( { - 2\frac{{{x^2}}}{{{r^2}}}} \right) (4) 式中，激光功率P=2000W，反射率R=0.74，光斑半径r=0.25mm。
1.3 初始条件和边界条件
(1) 初始值速率u=0m/s，压力p=1×105Pa，T为环境温度，则有T=300K。
(2) 温度边界条件设置为绝热，边界方程为：
- \kappa \cdot T = 0 (5) (3) 速度边界条件设置为无滑移边界，边界方程为:
\mathit{\boldsymbol{u}} = 0 (6) 1.4 空气与GH4037镍基高温合金物理参量
Table 1. Physical properties of airT/K cp/(J·kg-1·K-1) κ/(W·m·-1K-1) μ /(Pa·s) 300 1.005×103 0.0262 1.983×10-5 500 1.029×103 0.0403 2.671×10-5 1000 1.141×103 0.0675 4.152×10-5 1500 1.230×103 0.0946 5.400×10-5 2000 1.338×103 0.1240 6.500×10-5 2500 1.688×103 0.1750 7.670×10-5 随温度变化的空气密度ρ可由以下表达式计算：
\rho = p \cdot M/(\alpha \cdot T) (7) 式中，空气摩尔质量M=0.0289kg/mol，通用气体常数α=8.314J/(mol·K)。
Table 2. Physical properties of GH4037 Nickel-based superalloypropertiy value specific heat of solid phase cp, s 440J/(kg·K) thermal conductivity of solid phase κs 13.8W/(m·K) thermal conductivity of melting phase κm 23.9W/(m·K) dynamic viscosity μ 0.006Pa·s latent heat of melting Hm 296kJ/kg latent heat of vaporization Hv 6423kJ/kg melting temperature Tm 1585K vaporization temperature Tv 3005K 采用等效热熔法处理材料熔化和气化过程的相变潜热的影响，等效热熔法的表达式如下[20]：
\begin{array}{l} {c_p} = {c_{p, s}} + \left( {\frac{{{H_{\rm{m}}}}}{{{T_{\rm{m}}}}}} \right)H(\left( {T - {T_{\rm{m}}}), \Delta T} \right) + \\ {\delta _{\rm{m}}}{H_{\rm{m}}} + \left( {\frac{{{H_{\rm{v}}}}}{{{T_{\rm{v}}}}}} \right)H((T - {T_{\rm{v}}}), \Delta T) + {\delta _{\rm{v}}}{H_{\rm{v}}} \end{array} (8) {\delta _{\rm{m}}} = \frac{{{\rm{exp}}\left[ {{{(T - {T_{\rm{m}}})}^2}/\Delta {T^2}} \right]}}{{\Delta T\sqrt {\rm{ \mathsf{ π} }} }} (9) {\delta _{\rm{v}}} = \frac{{{\rm{exp}}\left[ { - {{(T - {T_{\rm{v}}})}^2}/\Delta {T^2}} \right]}}{{\Delta T\sqrt {\rm{ \mathsf{ π} }} }} (10) 式中，H((T-Tm), ΔT)和H((T-Tv), ΔT)均为Hea-viside平滑函数，分别用于处理材料在固-液相变和液-气相变过程中的相变潜热; δm和δv为高斯函数，高斯函数的中心分别为材料的熔点和沸点; ΔT为相转变温度范围，本文中取ΔT=50K。
在材料发生固-液相变和液-气相变过程中，材料的密度ρ也会发生两次阶跃性跳变，容易导致求解过程出现奇异。本文中也采用Heaviside平滑函数对其进行处理，改善了模型的收敛性。经平滑处理后的GH4037镍基高温合金的密度ρ如图 2所示。
2. 模拟结果与分析
2.1 温度场模拟结果
图 3为激光作用0.80ms, 1.20ms, 1.60ms和1.70ms左右的温度场模拟结果。图 3中的横纵坐标分别代表几何区域的长和高，单位为mm，色柱表示温度，单位为K。如图 3a、图 3b所示，在激光加热早期阶段，激光束主要用于加热基体，激光束下方的材料达到融化温度后，开始形成小的熔池，此时激光能量以菲涅尔吸收的方式被材料表面吸收，其余能量被反射，传热过程以热传导方式为主。随着加热时间增加，热传导现象加剧，熔池的深度随之增加。如图 3c所示，激光加热1.60ms左右，材料表面发生气化现象，在小孔内外压力差的作用下，气化材料沿垂直方向向上喷出，并与周围空气形成了热对流。如图 3d所示，激光加热1.70ms左右，材料气化蒸发现象更加明显，气化后的材料进一步上升，并向周围扩散，形成了类似“蘑菇云”的形状。
![Figure 3. Simulation results of temperature field]()
a—t=0.80ms b—t=1.20ms c—t=1.60ms d—t=1.70msFigure 3. Simulation results of temperature field2.2 速度场模拟结果
图 4a、图 4b分别为熔池表面和熔池内部不同位置(x, y)的水平流动速率。从图 4a可以看出，熔池中心区域的水平流动速率较快，两端的水平流动速率较慢，水平流动速率最高为1.7m/s。从图 4b可以看出，熔池表面的水平速率较高，越往下水平速率越低。
结合图 4a、图 4b可以看出，在激光加热早期阶段，熔池表面和内部的水平流动速率都较低，随着激光作用时间增加，熔池表面和内部的水平流动速率开始增加，但增速均较为缓慢。在激光加热1.70ms左右，熔池表面和内部的水平流动速率都迅速上升。
图 5为熔池表面的垂直流动速率。从中可以看出，在激光加热0.50ms左右，熔池开始出现垂直流动。随着时间的增加，熔池的垂直流动速率不断增加，在激光加热1.70ms左右，熔池的垂直速率也迅速上升。从图 5中还能看出，靠近熔池中心区域的垂直流动速率较快，熔池两端的垂直流动速率较慢，垂直流动速率最高为1.1m/s。
图 6为材料表面不同位置的压力随时间变化的曲线。从图 6可以看出，在激光加热初始阶段，此时的压力与大气压力相同，然后随着激光加热时间增加，压力开始缓慢的上升，在激光加热1.70ms左右，压力开始迅速上升，中心点压力最大达到1.38×105Pa，而且距离中心区域越近最大压力越大。
对比熔池水平和垂直流动速率与压力的模拟结果可以看出，在激光加热1.70ms左右压力迅速上升，距离中心点越近最大压力越大，此时速率场模拟中的熔池流速也在迅速上升，且越靠近熔池中心区域流动速率较快，表明蒸发气体的反冲压力可以加快熔池的流动。
图 7为激光作用1.60ms, 1.64ms, 1.66ms和1.70ms时的材料气化蒸发速率场模拟结果。图中的横纵坐标分别代表几何区域的长和高，单位为mm，色柱表示速率，单位为m/s。从图上可以看出, 在材料发生气化后，材料的最大气化蒸发速率由1.60ms时的35m/s迅速增加到1.70ms时的250m/s。靠近材料表面中心区域的材料蒸发速率最高，向周围扩散的气化蒸发材料在空气的阻力作用下速率逐渐递减。
![Figure 7. Simulation results of vertical field]()
a—t=1.60ms b—t=1.64ms c—t=1.66ms d—t=1.70msFigure 7. Simulation results of vertical field3. 结论
基于流体力学和流体传热理论建立了GH4037镍基高温合金激光打孔相变模型，计算了激光打孔相变过程的温度场分布、熔池流动速率以及气化蒸发速率，模拟结果较好地展示了材料熔化和气化蒸发的相变过程。
(1) 随着激光作用时间增加，材料相继发生熔化和气化现象，熔化材料初期在粘滞力和重力的阻碍作用下流动较为缓慢，然后在气化蒸发材料的反冲压力下，熔池的流动速率迅速增加，气化蒸发材料沿垂直方向从小孔喷射到空气中，并向周围空气中扩散。
(2) 在激光功率为2000W、脉宽为1.70ms时，熔池的水平流动最高速率为1.7m/s，垂直流动最高速率为1.1m/s，材料最大气化蒸发速率可以达到250m/s。
-
![]()
Figure 3. a~c—the simulated second harmonic signal of formaldehyde with noise level A=10, 20, 50, respectively e~f—the corresponding fitted results by using LSM and RANSAC algorithms, respectively
![]()
Figure 4. a—the experimentally measured second harmonic signal of formaldehyde with different concentrations b—the fitting results by using LSM and RANSAC algorithms, respectively
![]()
Figure 5. a—the experimentally measured second harmonic signal of formaldehyde with same concentration b—the corresponding fitted results by using LSM and RANSAC algorithm, respectively
Table 1 Fitting results of two harmonic signal of formaldehyde simulation with different SNR (adding Gaussian noise with different amplitude A to Y, while keeping X unchanged)
experimental results R2 slope LSM RANSAC LSM RANSAC A=0(RSNR=∞) 1.0 1.0 2.0 2.0 A=5(RSNR=8.87) 0.9833 0.9935 1.9845 2.0059 A=10(RSNR=6.00) 0.9371 0.9831 2.0183 1.9997 A=20(RSNR=3.54) 0.8045 0.9412 2.0252 2.0003 A=50(RSNR=0.64) 0.2654 0.7482 1.6848 1.9083 
下载: 导出CSV
Table 2 Linear fitting results of formaldehyde spectra under different experimental conditions
formaldehyde sample 1 sample 2 LSM RANSAC LSM RANSAC actual ratio 1.20 1.20 1.0 1.0 fitted value 1.1818 1.1818 0.5342 1.0272 correlation coefficient R2 0.8923 0.9853 0.2111 0.9743 error/% 1.517 -0.18 46.58 -2.72
下载: 导出CSV
-
[1] PHILIPPE L C, HANSON R K. Laser diode wavelength-modulation spectroscopy for simultaneous measurement of temperature, pressure, and velocity in shock-heated oxygen flows[J]. Applied Optics, 1993, 32(30):6090-6103. DOI: 10.1364/AO.32.006090
[2] LI J S, YU B L, ZHAO W X, et al. A review of signal enhancement and noise reduction techniques for tunable diode laser absorption spectroscopy[J]. Applied Spectroscopy Reviews, 2014, 49(8):666-691. DOI: 10.1080/05704928.2014.903376
[3] REID J, LABRIE D. Second-Harmonic detection with tunable diode lasers-comparison of experiment and theory[J]. Applied Physics, 1981, B26(3):203-210. DOI: 10.1007/BF00692448
[4] LI J, REIFFS A, PARCHATKA U, et al. In situ measurements of atmospheric CO and its correlation with NOx and O3 at a rural mountain site[J]. Metrology and Measurement Systems, 2015, 22(1):25-38. DOI: 10.1515/mms-2015-0001
[5] CAI Y, WU Sh Q, WU A, et al. Study on calculation method of detection limit based on wavelength modulation spectroscopy[J]. Laser Technology, 2012, 36(3):390-393(in Chinese). http://en.cnki.com.cn/Article_en/CJFDTotal-JGJS201203028.htm
[6] XU Y Z, GUO J Q, GAO X R, et al. Effect of temperature on absorption spectral lines of carbon monoxide[J]. Laser Technology, 2010, 34(6):778-780(in Chinese). http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=jgjs201006016
[7] PLACKETT R L. The discovery of the method of least squares[J]. Biometrika, 1972, 59(2):239-251. http://d.old.wanfangdata.com.cn/OAPaper/oai_arXiv.org_1310.3198
[8] JIA X Y, XU C S, BAI X. The invention and way of thinking on least squares[J]. Joumal of Northwest University, 2006, 36(3):507-511(in Chinese). http://www.ams.org/mathscinet-getitem?mr=2254113
[9] FISCHLER M A, BOLLES R C.Random sample consensus:a paradigm for model fitting with applications to image analysis and automated cartography[J]. Communications of the ACM, 1981, 24(6):381-395. DOI: 10.1145/358669.358692
[10] ZHOU C L, ZHU H H, LI X J. Research and application of robust plane fitting algorithm with RANSAC[J]. Computer Engineering and Applications, 2011, 47(7):177-179(in Chinese). http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=jsjgcyyy201107051
[11] CAO Y, FENG Y, YANG Y T, et al. Application of estimation algorithm based on RANSAC in road points cloud optimization[J]. Infrared and Laser Engineering, 2012, 41(11):3108-3112(in Chinese). http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=hwyjggc201211047
[12] WEI Y Z, LIU X L. Robust plane fitting of clouds based on RANSAC[J]. Journal of Beijing University of Technology, 2014, 40(3):400-403(in Chinese). http://en.cnki.com.cn/Article_en/CJFDTOTAL-BJGD201403015.htm
[13] ZHEN Y, LIU X J, WANG M Zh. An improved RANSAC of fundamental matrix estimation method[J]. Bulletin of Surveying and Mapping, 2014(4):39-43(in Chinese). http://d.old.wanfangdata.com.cn/Periodical/chtb201404011
[14] ZHANG H M, ZHENG Z. An improvement of the adjacent probability random sampling consistency algorithm[J]. Laser Journal, 2013, 34(5):29-30(in Chinese). http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=jgzz201305013
[15] HAST A, NYSJÖ J, MARCHETTI A. Optimal RANSAC-towards a repeatable algorithm for finding the optimal set[J]. Journal of WSCG, 2013, 21(1):21-30.
[16] LI J S, PARCHATKA U, FISCHER H. A formaldehyde trace gas sensor based on a thermoelectrically cooled CW-DFB quantum cascade laser[J]. Analytical Methods, 2014, 6(15):5483-5488. DOI: 10.1039/C3AY41964A
[17] LI J, PARCHATKA U, FISCHER H. Development of field-deployable real time QCL spectrometer for simultaneous detection of ambient N2O and CO[J]. Sensors and Actuators, 2013, B182(3):659-667. http://www.sciencedirect.com/science/article/pii/S0925400513003547
计量
- 文章访问数: 7
- HTML全文浏览量: 1
- PDF下载量: 5