摘要:利用光声光谱技术进行痕量气体的检测具有独特的优势, 光声池是系统装置中最为重要的核心部件, 它决定着整机性能的优劣. 以一圆柱形共振型光声池为研究对象, 基于声学与吸收光谱学的基本理论, 建立了光声池声场激发的数学模型; 利用数值模拟方法对光声池空腔结构进行了声学模态仿真, 获得了前8阶声学模态值以及声压可视化振型; 在考虑热黏性声学损耗的作用下, 对光声池进行了热-声耦合多物理场仿真计算;将仿真结果与解析计算和实验结果进行对比, 明确了利用数值模拟方法来计算光声池有关指标的可靠性与可行性; 针对光声池的优化问题, 提出了一种将响应面代理模型与遗传算法相结合的优化算法, 在将原光声池中的谐振腔两端形貌更改为喇叭口形的情况下, 通过优化算法获得了以光声池品质因数Q及池常数C_cell为最大值寻优的Pareto最优解集; 选取一组解进行考察, 结果表明, 代理模型预测值与数值模拟值指标最大误差仅为1.3%, 优化后的新型光声池Q较之前增长了48.9%, C_cell增长了34.4%. 研究方法可为光声光谱中光声池的优化设计提供参考借鉴.
关键词: 光学/
光声光谱/
光声池/
数值计算

English Abstract

全文HTML

--> --> -->

随着现代化工业的突飞猛进, 痕量气体检测技术的重要性不言而喻, 目前痕量气体检测技术已广泛应用于农业生态^[1]、环境监测^[2]、生物医疗^[3]、高压变电^[4]、化工环保^[5]等领域. 光声光谱(photoacoustic spectroscopy, PAS)是以吸收光谱为基础理论, 基于光热声复杂效应来实现对气体定性与定量测量的一种检测技术. 1880年, Bell在研究无线通话技术的过程中发现固体的光声效应, 紧接着又发现了气体与液体的光声光谱效应, 由此拉开了光声光谱科学研究的序幕. 然而, 限于当时技术的制约, 研究进展甚为缓慢, 直到20世纪60年代, 光声光谱技术的研究才重新活跃起来, 80年代以后, 出于光声光谱检测技术的独特性能, 光声光谱技术进入了快速发展期, 各类商业化产品也相继问世. 光声光谱技术主要有高灵敏度、高选择性、快速响应、在线监测和不消耗气样等优点, 其应用与科学研究前景十分广阔.
当前, 国内外****从理论计算、技术应用、实验设计等方面对光声光谱的检测装置做出了研究, 并取得了丰硕的成果. Kreuzer^[6]最早利用3.39 $\mu$, 15 mW He-Ne激光器建立了CH₄光声光谱探测系统, 检测灵敏度达到10 × 10^–9, 并从理论上分析了光声光谱检测精度可达到0.1 × 10^–12; Besson等^[7]搭建了多组分光声光谱检测装置, 并从理论和实验两方面分析了缓冲气体对光声信号的影响; Tavakoli等^[8]提出了一种LC电路模型对一维声谐振器声学共振频率、品质因素进行了模拟, 并计算出了传输噪声最小化的几何结构; Pernau等^[9]设计制作了4.3${\text{μ}{\rm m}} $的中红外发光二极管(LED)光源式CO₂共振型光声光谱检测装置; Baumann等 ^[10]提出了针对光声池光声信号进行了仿真计算, 仿真结果与实验结果较为吻合. 国内陈伟根研究组^[11]构建了基于分布反馈半导体激光器便携式光声光谱检测装置, 探讨分析了气体光声信号多种因素的影响特性; 高晓明研究组^[12,13]建立了7.6 ${\text{μ}{\rm m}}$中红外量子级联激光的N₂O光声光谱实验系统与单探测器同时探测H₂O, CO₂和CH₄的三通道光声光谱检测装置; 于清旭研究组^[14]研制了基于可调谐掺铒光纤激光器的共振式光声光谱乙炔气体检测系统, 并对多组分气体光声光谱的独立成分分析方法进行了研究; 董磊研究组^[15]报道了拍频石英音叉增强型光声探测技术, 可用于超灵敏无标定痕量气体检测和快速光谱扫描应用; 马欲飞研究组^[16]采用波长为4.6${\text{μ}{\rm m}} $的新颖中红外高功率分布反馈量子级联激光器作为光源, 利用石英增强光声光谱技术实现了对CO气体基频吸收带的激发与测量, 探测极限达1.95 × 10^–9.
综上, 当前文献主要是围绕光声光谱检测装置的设计和应用展开的, 有关光声光谱数值模拟与优化设计方面的研究内容偏少, 实际设计中, 针对圆柱型甚至球型光声池^[17]的解析计算结果尚可以接受, 然而针对复杂形体的光声池声场偏微分方程却难以建立与求解. 采用数值模拟的方法不仅可以获得更高精度的计算结果, 同时还可以直观地表征物理量的分布, 对于光声光谱产品的开发而言, 利用仿真模拟可有效地减少设计成本、缩短研发周期和增加产品设计的可靠性与优化性.
本文以一个圆柱形共振型光声池为研究对象, 按照1︰1比例建立了其求解域模型, 利用多物理场数值计算方法对光声光谱中光热声效应中声腔模态与声场分布进行计算分析, 并与实验结果进行了比较, 验证了多物理场数值计算方法的可靠性与准确性. 在此基础上, 对原圆柱形共振型光声池进行了形状上的探索改造及优化, 并与原光声池进行了对比, 结果表明, 优化后的光声池性能得到一定幅度的提升, 研究方法可为光声光谱光声池的性能改善与优化工作提供设计思路.

2

-->

2.1.基本理论

光声光谱检测技术的基本理论为: 当待测气体分子吸收激光能量后会由基态跃迁至激发态, 处于激发态的气体分子通过无辐射弛豫退激后会释放能量, 释放的能量使得待测气体分子及其空间介质按激光调制的频率产生周期性加热, 在体积一定的情况下, 气体周期性的温度变化会引发相应的压力变化, 由此产生声场效应, 通过高灵敏度麦克风便可以感知声学信号, 以此规律来反演气体的浓度等其他物态信息, 检测原理图如图1所示.
图 1 光声光谱气体检测原理示意图
Figure1. Schematic diagram of photoacoustic spectroscopy gas detection principle.

2 -->

2.2.系统结构

基于光声光谱的基本理论, 本文构建了一套用于气体检测的光声光谱实验装置, 并以甲烷气体为检测对象. 检测系统主要由分布式反馈激光器(distributed feedback laser, DFB)激光器, DFB半导体激光器(NEL, 1653 nm DFB)、激光控制器(THORLABS, CLD1015)、函数发生器(南京盛普, SPF05)、锁相放大器(Stanford Research System, SR830)、自动配气系统(Environics, S4000)、光声池(自制)、高灵敏度麦克风(北京声望, MP201)和采集卡(NI, USB6212)等组成, 实验装置实物图如图2所示.
图 2 光声光谱气体检测实验装置图
Figure2. Photoacoustic spectroscopy gas detection device.

2

-->

3.1.光声信号数学模型

基于声学基本方程与光谱吸收定律, 假定待测气体为理想性、各向同性、均匀且初始静止的连续性流体介质, 并假设介质中声波传输中的声压与气体密度变化为微扰动量, 根据流体质量守恒方程、动量守恒方程以及气体的物态方程, 可得到无激励源线性化气体声压波动方程:

${\nabla ^2}p - \frac{1}{{{v^2}}}\frac{{{\partial ^2}p}}{{\partial {t^2}}} = 0.$

针对线性声场在谐激励下的稳态时域响应, 将时间项从(1)式中消去, 可得

$ \left({\nabla ^2} + \frac{{{\omega _j}^2}}{{{v^2}}} \right){p_j} = 0.$

当激光光源通过光声光谱系统中光声池时, 由于激励源的缘故, 气体声压波动方程为^[18]

${\nabla ^2}p - \frac{1}{{{v^2}}}\frac{{{\partial ^2}p}}{{\partial {t^2}}} = - \left[ {\frac{{ \left(\gamma - 1 \right)}}{{{v^2}}}} \right]\frac{{\partial H}}{{\partial t}}.$

对(3)式进行傅里叶变换, 可得

$\left({\nabla ^2} + \frac{{{\omega ^2}}}{{{v^2}}} \right)p(\omega ) = - \left[ {\frac{{(\gamma - 1)}}{{{v^2}}}} \right]{\rm{i}}\omega H(\omega ).$

将$p(\omega )$按简正模态${p_j}$的叠加形式展开:

$p(\omega ) = \sum {{A_j}(\omega ){p_j}} , $

其中, $p$为$p(x, y, z, t)$, 描述的是空间与时间状态下的声压; ${p_j}$为简正模态; ${A_j}(\omega )$为第$j$阶模态的振幅; $v$为待测气体中的声速; $\omega $为圆频率; $\gamma $为热容比, $\gamma = {C_p}/{C_v}$, ${C_p}$为定压热容, ${C_v}$为定容热容; $H$为热功率密度源.
本文初始研究的光声光谱装置中光声池形状为圆柱形, 因而可利用柱坐标系$(r, \theta , z)$来表示(2)式:

$\frac{1}{r}\frac{\partial }{{\partial r}}\left(r\frac{{\partial {p_j}}}{{\partial r}} \right) + \frac{1}{{{r^2}}}\frac{{{\partial ^2}{p_j}}}{{\partial {\theta ^2}}} + \frac{{{\partial ^2}{p_j}}}{{\partial {z^2}}} + {\left(\frac{{{\omega _j}}}{v} \right)^2} = 0, $

${p_j} = {}_{\sin }^{\cos }(m\theta )[A{{\rm{J}}_m}({k_r}r)] \cdot [B\sin ({k_z}z) + C\cos ({k_z}z)], $

其中, ${{\rm{J}}_m}$为第一类贝塞尔函数; A, B, C分别由声学边界条件决定的系数; ${k_z} = ({\text{π}}/{L_{\rm{c}}})q$, ${k_r} =$$ ({\text{π}}{\alpha _{mn}})/{R_{\rm{c}}}$, [qmn]为描述${p_j}$模态纵向、径向和角向的振动模式特征数; ${\alpha _{mn}}$为含m阶贝塞尔函数方程的第n个根.

$p({{r}}, \omega ) = \sum {{A_j}(\omega ){p_j}({{r}})} , $

${A_j}(\omega ) = \frac{{ - {\rm{i}}\omega }}{{\omega _{_j}^2}}\frac{{[(\gamma - 1)/{V_{\rm{c}}}]\displaystyle\iiint {p_{_j}^ * ({{r}})H({{r}}, \omega ){\rm{d}}V}}}{{1 - {{\left( {\displaystyle\frac{\omega }{{{\omega _j}}}} \right)}^2}}}.$

实际中, 需考虑热传导与黏滞损耗:

${A_j}(\omega ) = \frac{{ - {\rm{i}}\omega }}{{\omega _{_j}^2}}\frac{{[(\gamma - 1)/{V_{\rm{c}}}]\displaystyle\iiint {p_{_j}^ * ({{r}})H({{r}}, \omega ){\rm{d}}V}}}{{1 - {{\left( {\displaystyle\frac{\omega }{{{\omega _j}}}} \right)}^2} - i\left( {\displaystyle\frac{\omega }{{{\omega _j}{Q_j}}}} \right)}}, $

$H({{r}}, \omega ) = \alpha I({{r}}, \omega ) = N\sigma cI({{r}}, \omega ) = N\sigma c{P_{{\rm{laser}}}}g({{r}}, \omega ).$

令: $U = \displaystyle\iiint {p_{_j}^ * ({{r}})g({{r}}, \omega ){\rm{d}}V}$,
若调制频率$\omega = {\omega _j}$时, 光声池中谐振腔将发生声学共振, 由(8)和(10)式可得

$p({{r}}, {\omega _j}) = - (\gamma - 1)\frac{{{Q_j}}}{{{\omega _j}}}\frac{{N\sigma c{P_{{\rm{laser}}}}U}}{{{V_{\rm{c}}}}}.$

令: ${C_{{\rm{cell}}}} = - (\gamma - 1)\displaystyle\frac{{{Q_j}}}{{{\omega _j}}}\displaystyle\frac{U}{{{V_{\rm{c}}}}}$,

$p({{r}}, {\omega _j}) = {C_{{\rm{cell}}}}N\sigma c{P_{{\rm{laser}}}}.$

对于两端开口的圆柱形光声池的简正模态频率计算公式为:

$f = \frac{{{\omega _j}}}{{2{\text{π}} }} = \frac{v}{2}\sqrt {\left[{{\left(\frac{{{a_{mn}}}}{{{R_c}}}\right)}^2} + {{\left(\frac{q}{{{L_{{\rm{eff}}}}}}\right)}^2}\right]} , $

${L_{{\rm{eff}}}} = {L_{\rm{c}}} + \frac{{16}}{{3{\text{π}}}}{R_{\rm{c}}} , $

其中, $p_{_j}^ * ({{r}})$为${p_j}({{r}})$复数共轭模式; ${Q_j}$为${p_j}({{r}})$模式下的品质因数; $N$为气体总分子数密度; $\sigma $为气体分子的吸收截面; $c$为待测气体分子的体积浓度; $I({{r}}, \omega )$为激光入射光强; ${V_{\rm{c}}}$为光声池谐振腔体积; ${R_{\rm{c}}}$为谐振腔半径; ${L_{\rm{c}}}$为谐振腔长度; ${L_{{\rm{eff}}}}$为考虑谐振腔端部效应后的修正长度; ${C_{{\rm{cell}}}}$为池常数; $g({{r}}, \omega )$为归一化光强度分布函数.
通常情况下, 圆柱形共振式光声池一般工作在一阶纵向简正模式, 即[qmn]为[100]模式, 因而可计算出:

${f_{100}} = \frac{v}{{2{L_{{\rm{eff}}}}}}, $

${Q_{100}} = \frac{{{R_c}}}{{{d_{{\rm{visc}}}} + (\gamma - 1){d_{{\rm{therm}}}}(1 + 2{R_{\rm{c}}}/{L_{\rm{c}}})}}, $

${C_{100}} = \frac{{2(\gamma - 1){L_{{\rm{eff}}}}{Q_{100}}}}{{{\text{π}}{V_{\rm{c}}}{\omega _{100}}}}, $

其中, ${Q_{100}}$, ${C_{100}}$为一阶纵向简正模式的品质因数与池常数; ${d_{{\rm{visc}}}}$为黏性边界层厚度, ${d_{{\rm{therm}}}}$为热边界层厚度.
2 -->

3.2.物理模型与基本参量

光声池是光声光谱气体检测装置中最重要的核心部分, 它决定着整个检测系统的性能的优劣. 本实验装置中所采用的光声池为圆柱形, 其内部空间是一个封闭的空腔, 空腔内填充待测气体, 光声池主要由谐振腔、缓冲腔、入射窗口、出射窗口、进、出气口等部分组成. 初始结构参量如下: 光声池外壳材质为硬铝, 谐振腔体材质为紫铜, 窗口采用透过率大于90%的石英窗片密封, 谐振腔半径为5 mm, 长度为100 mm, 缓冲腔半径为20 mm, 长度为50 mm, 忽略一些微小细节, 如较小倒角、圆角、用于装配螺钉、螺母及其对应的孔等, 在基本假设的前提下, 对光声池物理模型作出一定简化, 建立的光声池空腔声学物理模型如图3所示.
图 3 光声池空腔声学物理模型
Figure3. Acoustic physical model of photoacoustic cell cavity.

2

-->

4.1.空腔声学模态数值计算

光声池空腔结构具有声学模态频率与相应的振型, 声学模态可以表征其内部空腔的压力分布特性, 当声波以某一模态频率在空腔内传播时, 光声腔将随之发生声学共振, 求解光声池空腔模态, 即解出(5)式中p_j(r), 根据(1)式及界面法向速度及法向力连续条件, 并采用有限元离散化方式, 建立声场的数值方程^[19]如(19)式:

$[{{M}}_e^p]\left\{ {{{\ddot p}_e}} \right\} + [{{C}}_e^p]\left\{ {{{\dot p}_e}} \right\} + [{{K}}_e^p]\left\{ {{{{p}}_e}} \right\} = \left\{ 0 \right\}, $

空腔壁设定为硬声场边界, 进一步可得

$[{{M}}_e^p]\left\{ {{{\ddot p}_e}} \right\} + [{{K}}_e^p]\left\{ {{{{p}}_e}} \right\} = \left\{ 0 \right\}, $

其中, $[{{M}}_e^p]$, $[{{C}}_e^p]$, $[{{K}}_e^p]$分别为等效流体质量矩阵、阻尼矩阵和刚度矩阵; $\left\{ {{{{p}}_e}} \right\}$为单元节点声压矩阵.
进一步进行频域下求解:

$ \left( \left[{{K}}_e^p \right] - {\omega ^2} \left[{{M}}_e^p \right] \right)\left\{ {{{{p}}_e}} \right\} = \left\{ 0 \right\}.$

由(21)式可得出声腔共振圆频率${\omega _i}$和声压分布${p_j}$ (i = 1, 2, 3…, n).
利用多物理场仿真程序进行数值求解, 定义光声池空腔中流体材料为氮气, 密度为1.25 g/L, 声速为349 m/s, 采用仿真平台中物理场智能控制网格方式进行网格划分, 四面体网格单元数为28711, 三角形单元5302, 光声池声学模态仿真求解结果如图4所示.
图 4 光声池空腔声学模态仿真云图
Figure4. Acoustic mode simulation of photoacoustic cell cavity: (a)The first mode (265 Hz); (b) the second mode (1659 Hz); (c) the third mode (3125 Hz); (d) the fourth mode (3490 Hz); (e) the fifth mode (3680 Hz); (f) the sixth mode (4966 Hz) ; (g) the seventh mode (5123 Hz); (h) the eighth mode (6207 Hz).

从图4中光声池空腔各阶声学模态振型中可知, 在前8阶模态中, 第1—6阶为纵向模态振型, 第7阶为角向模态振型, 第8阶为纵向与角向的混合振型, 由于结构的对称性, 第7阶与8阶还存在模态频率值相同且为对称性的模态振型(对称模态振型本文只选择其一来说明), 但可以预见, 随着模态阶数的增加, 光声池声学模态振型将越来越复杂, 在仿真范围内没有出现径向模态振型. 由于受缓冲腔的影响, 光声池纵向声学模态与单纯圆柱形光声池的计算结果稍有不同, 仿真结果中出现低频声学共振是由于缓冲腔造成的, 但这不影响后续的分析研究; 图4仿真云图中颜色从红到蓝依次表示声压的相对强弱, 颜色接近于深蓝的地方声压为0, 其连线称为光声池声腔模态的节线, 当光声池受到外界激励源时, 红色区域处声压响应相对大, 蓝色区域处声压响应相对小, 节线位置处没有响应. 因而通过模态分析可知, 在设计光声麦克风采音位置时, 应尽量选择靠近声腔模态中声压相对较强的区域, 设计进、出气口时, 应尽量选择靠近声腔模态中节线的位置. 本文光声池的信号采样选取的是仿真结果中第2阶模态频率与振型, 其重要特征为麦克风安装于谐振腔的中部位置, 为区别有无缓冲腔带来的效应影响, 本文将光声池工作在上述状态下的声学模态称之为“工作纵向声学简正模态”.
2 -->

4.2.光声信号数值计算与分析

4.1节中光声池声学模态频率与模态振型的求解结果是基于未施加激励源的条件, 所获得的声压分布均为相对值, 计算光声池通过激光调制并由待测气体吸收后产生的绝对声压值, 还需要将激光激励源、待测气体物理参数和热黏性声学损耗等添加于仿真模拟中.
3 -->

4.2.1.计算参数设置与网格划分

利用多物理场仿真程序对光声池热黏性声学进行数值求解, 黏性边界层厚度与热边界层厚度解析式计算公式^[18]如下:

${d_{{\rm{visc}}}} = \sqrt {\frac{{2\mu }}{{\omega {\rho _0}}}} , $

${d_{{\rm{therm}}}} = \sqrt {\frac{{2k}}{{\omega {\rho _0}{C_p}}}} , $

其中, $\mu $为动力黏度, ${\rho _0}$为气体密度, $k$为气体的导热系数.
在 1659 Hz, 20 ℃, 1 atm环境中, 通过粗略分析, 光声池中的黏性边界层及热边界层厚度数量级为10^–2 mm, 因此在光声池三维声学模型的情况下, 直接利用“热黏性声学”接口程序仿真求解时, 计算机运行成本太高, 结合光声池对称式的结构特点, 本文将光声池三维模型转化为二维轴对称模型来进行热黏性模拟仿真计算, 如此在保持一定求解精度的情况下可以节省计算机仿真成本与计算时间, 仿真模型如图5所示.
仿真中参量与仿真条件如下: 待测气体为甲烷, 体积浓度为100 × 10^–6; 激光功率为4.5 mW, 激光光源分布为高斯分布, 束腰半径为0.9 mm, 中心波长为1653 nm; 选用甲烷2v3带R(3)支吸收谱线, 通过HITRAN数据库查找谱线的吸收截面, 根据 (11) 式计算出热功率密度源H; 在光声池壁处, 对速度场施加无滑移边界条件, 对温度场施加等温边界条件; 对热黏性声学模型进行恰当的网格剖分非常重要, 这也是仿真求得准确解的关键之一, 黏滞边界层和热边界层在仿真中必须被网格解析, 本文采用手动控制网格单元尺寸与边界层网格的混合控制剖分技术, 预定义普通物理场极细化网格划分, 边界层采用10层网格, 第一层厚度设为0.01 mm, 边界层拉伸因子为1, 网格统计为, 单元总数为3462, 网格顶点2641, 其中四边形单元数为1500, 三角形单元数为1962, 谐振腔中部剖分细节如图6所示.
图 5 光声池二维轴对称物理模型
Figure5. Two dimensional axisymmetric physical model of photoacoustic cell.

图 6 谐振腔中部混合网格剖分细节图
Figure6. Detail drawing of mixed meshes in the middle of resonator.

3 -->

4.2.2.扫频结果

在模拟计算中, 提取光声池谐振腔正交中心线上最大声压值, 首先以粗频率(range 50—6500 Hz; step 50 Hz)对光声池进行频率扫描, 粗频率响应曲线如图7所示. 从曲线上可知, 粗频扫描情况下在1650 Hz与4950 Hz处出现了声压峰值, 说明在此频率上出现了声学共振, 这与4.1节光声池第2阶与6阶的声学模态计算结果非常吻合, 共振频率值略为减小, 原因在于声学共振增强了热声效应, 从而降低了其共振频率. 在黏性损耗和热损耗效应的作用下, 模态中其他振型并没有被有效激发, 由于研究的光声池光声信号采样选取的是第2阶模态频率与振型, 因而在1650 Hz附近对光声池进行频率细扫描(range 1550—1750 Hz; step 5 Hz), 并对扫频结果进行Lorentz函数拟合, 细频率响应曲线如图8所示, 从曲线上可知, 细扫频结果的共振峰值对应频率为1650 Hz, 拟合结果为1648 Hz, 拟合程度非常显著(R² = 0.9968), 可以计算出半功率处的全频线宽为25.9 Hz, 利用拟合结果可计算出该光声池的品质因数Q = 63.7.
图 7 光声池粗频响应仿真曲线
Figure7. Simulation curve of photoacoustic cell's coarse frequency response.

图 8 光声池细频响应仿真曲线
Figure8. Simulation curve of photoacoustic cell's fine frequency response.

3 -->

4.2.3.声场与温度场分布结果

频扫均设置为(range 1550—1750Hz; step 5Hz), 扫描探测路径上的声压和温度场分布及大小, 声压探测路径为光声池空腔中沿轴纵向上与紫铜接触的边界线, 即谐振腔的母线, 其分布特征可以表征光声池光声信号的转化能力. 如图9仿真结果所示, 当扫频值等于声学共振频率时, 声压值较非共振情况得到显著增益, 因而在系统实际检测中, 应该将光源调制频率准确地落在光声池简正频率上, 这样可以获得较高的系统检测灵敏度; 在不同探测频率下沿光声池谐振腔母线的声压特性曲线基本相同, 均呈现二次抛物线形状, 沿谐腔轴线方向的两侧声压呈对称分布, 在谐振腔的中部为最大声压值, 此位置也正是高灵敏麦克风的探测位置.
图 9 光声池谐振腔母线声压特性曲线
Figure9. Sound pressure characteristic curve of cavity geometry of photoacoustic cell.

温度场探测路径为光声池谐振腔轴线正交中心线, 即安装的麦克风的轴线上, 仿真结果如图10所示. 从图10可知, 在不同探测频率下光声池谐振腔轴线正交中心线的温度场特性曲线趋势基本相同, 在r = 0—1.5 mm区间上, 温升趋于减弱状态, 其原因在于光声池轴线上的激光光源引起的, 所设激光光源的束腰直径为1.8 mm且光强为高斯分布, 接近激光中心的温度较强, 远离激光中心的温升稍微减弱, 这与仿真曲线较为吻合, r = 2.1—4.9 mm区间上, 温升趋于平稳状态, r = 4.9—5 mm区间上, 温升发生突变, 其接近外壁处温升梯度较大, 此区间为黏性边界层与热边界层的厚度, △r ≈ 0.06 mm, 该区域内黏性损耗与热损耗较大, 通过(22)和(23)式可计算两类边界层厚度约为0.06 mm, 这与仿真结果较为吻合, 说明了数值模拟具有较高的可信度.
图 10 光声池谐振腔轴线正交中心线温升特性曲线
Figure10. Temperature rise characteristic curve of orthogonal center line of photoacoustic cavity resonator axis.

2 -->

4.3.实验与计算结果比较

根据3.1节光声信号解析计算与4.2节光声信号数值模拟结果, 基于搭建的光声光谱气体实验装置(如图2所示)进行相关的实验验证. 实验中选用氮气为缓冲气体, 由自动配气系统配成100 × 10^–6浓度的甲烷气体, 缓慢充入光声池腔内, 并关紧阀门形成密闭的光声池空腔, 保持腔内气压为1 atm, 调节激光器频率由1400 Hz缓慢增至1700 Hz, 通过粗调与精调相结合的方法, 每个频率粗调为1次采样光声信号值, 细调为20次采样平均. 对实验中光声光谱气体检测系统进行相关性能测试, 重点测量光声池的相关性能指标, 计算对比结果如图11所示. 从图11(a)中可知, 对于光声池一阶共振频率值, 其解析结果、数值模拟与实验结果三者较为吻合, 误差较小; 由图11(b)中可知, 对于光声池品质因素, 其理论解析与数值模拟结果非常接近, 但实验结果与上述两者存在一定误差, 原因在于光声池存在加工误差、谐振腔与缓冲腔的光洁度不够以及声速随温度与湿度影响的变化等, 理论计算存在一定的理想化误差; 由图11(c)中可知, 对于光声池池常数, 理论解析与数值模拟、实验结果具有一定偏差, 这是由于光声池存在各种误差导致的共振频率与各类损耗计算误差综合累积引起的, 同时光声池中缓冲腔的内壁也存在一定的表面损耗, 理想计算下并没有考虑; 图11(d)中可知, 对于光声池黏性边界层和热边界层厚度而言, 实验测量困难, 从图中给出的解析与数值模拟的结果可知, 声波在光声池壁周围的气体中传播时, 导热效应与粘滞性效应的影响程度相当, 光声池内腔表面热黏边界层仿真与解析计算结果几乎吻合, 误差较小; 图11(b)—(d)结果还需考虑光声池结构三维到二维的简化. 上述对比结果显示, 本文所建立的光声池物理模型与边界条件可以适用, 数值模拟方法可以为后续光声池的优化工作提供设计参考.
图 11 计算结果比较
Figure11. Comparison of calculation results: (a) Resonance frequency; (b) quality factor; (c) pool constant; (d) boundary layer thickness.

光声光谱气体检测系统在很大程度上主要取决于光声池的性能, 形状对于光声池内的声场分布影响较大, 合理设计光声池的形状可以获得优秀的增益效果, 但光声池的形貌经过改变后, 解析计算将变得更为复杂, 甚至得不到合理的解析结果, 此时利用数值模拟方法对复杂形状的光声池进行计算就显得更为方便、有效. 上述研究已经验证, 数值模拟技术可以很好地替代解析计算与物理样机实验, 它可以获得较为直观的物理场分布、发现物理量的特性规律, 在产品设计初期可提供较好的设计决策. 因此, 本文利用数值模拟技术与优化算法相结合的方法对光声池进行优化设计探讨.
2

-->

5.1.优化模型设定及流程

考虑到光声池的加工制造, 对于本实验中的光声池结构改进方案, 根据参考文献[20], 探索将原光声池中的谐振腔两端拟更改为喇叭口形, 如图12所示, 蓝色为原始结构剖面图, 红色为拟更改的内腔半剖图. 具体优化策略如下: 以光声池中的喇叭形状底圆半径r_c及其高度h_c, 谐振腔半径R_c, 谐振腔的长度L_c为变量因子, 并约束缓冲腔的形貌与几何尺寸不变、光声池轴向总长度不再增大, 利用DOE(design of experiment)试验设计方法选择分析样本, 建立设计变量与优化目标之间参数映射的响应面代理模型(surrogate model method), 最后通过多目标优化算法求得优化解.
图 12 光声池参数优化设计选取
Figure12. Optimum design parameters of photoacoustic cell.

2 -->

5.2.试验设计与代理模型

DOE采用Box-Behnken试验设计采样方法, Box-Behnken是一种具有高效率的响应曲面设计方法, 比较适应于二阶响应面模型的采样, 该采样方法包括中心点和设计空间的中心. 根据光声池原加工尺寸及制作经验确定设计变量的取值范围, 设计4 因素3水平的分析实验, 总计25组试验组合, 选取的试验设计水平表如表1所列.

因素	编码水平
	–1	0	1
A: 底圆半径rc/mm	5	6	7
B: 圆台高度hc/mm	8	10	12
C: 谐振腔半径Rc/mm	3	4	5
D: 谐振腔长度Lc/mm	60	70	80

表1因素水平表
Table1.The factors and levels graph.

由于光声池设计参数与其品质因素及池常数之间函数关系较为复杂, 因而代理模型选择二阶响应面模型来进行表示, 函数关系如下^[21]:

${f_u}(x) = {a_0} + \sum\limits_{i = 1}^n {{b_i}{x_i} + } \sum\limits_{i = 1}^n {{c_{ii}}x_{_i}^2 + } \sum\limits_{i < j}^n {{d_{ij}}{x_i}} {x_j} + \delta , $

式中, ${f_u}(x)$(u = 1, 2, 3)分别为光声池品质因素Q、池常数C_cell和第一阶纵向声学固有频率值f, 三者均为响应值; ${x_i}$为设计变量; $n$为设计变量个数; ${a_0}$, ${b_i}$, ${c_{ii}}$, ${d_{ij}}$分别为多项式回归系数; $\delta $为随机误差.
在保证材料属性、网格控制方法和边界条件不变的情况下, 通过数值模拟方法对上述25组试验组合下的光声池进行仿真求解并拟合代理模型, 拟合结果如表2所列. 为提高拟合模型的精度和简化模型公式, 代理模型中已筛选去除了影响不显著的部分因素组合. 表2结果表明, 三者响应面代理模型高度显著且拟合程度高, 可以满足后续优化设计的需求.

目标响应	相关系数R2	校正系数R2adj	P值
f1(x): Q	0.9997	0.9993	< 0.0001
f2(x): Ccell/(Pa·cm) · W–1	0.9999	0.9982	< 0.0001
f3(x): f/Hz	0.9998	0.9996	< 0.0001

表2代理模型拟合结果
Table2.Fitting results of surrogate models.

光声池品质因素f₁(x): Q, 池常数f₂(x): C_cell和工作纵向声学模态值f₃(x): f的代理模型公式如下:

$\begin{split}&{f_1}\left( x \right){\rm{ }} = {\rm{ }}16.93 + 2.39{\rm{A}} + 3.26{\rm{B}} + 29.19{\rm{C}}-\\&\quad 1.26{\rm{D}} + 0.73{\rm{AC}} + 0.137{\rm{BC}}-0.0137{\rm{BD}}-\\&\quad 0.145{\rm{CD}}-0.417{{\rm{A}}^2}-0.110{{\rm{B}}^2}-\\&\quad 0.804{{\rm{C}}^2} + 9.21 \times {10^{3}}{{\rm{D}}^2}, \end{split}$

$\begin{split}&{f_2}\left( x \right){\rm{ }} = {\rm{ }}3729.06-142.47{\rm{A}} + 324.62{\rm{B}}-\\&\quad 1746.72{\rm{C}} + 63.42{\rm{D}} + 18.51{\rm{AC}}-6.94{\rm{BC}}-\\&\quad 1.95{\rm{BD}}-3.22{\rm{CD}}-5.26{{\rm{B}}^2} +\\&\quad 163.55{{\rm{C}}^2}-0.196{{\rm{D}}^2}, \end{split}$

$\begin{split}&{f_3}\left( x \right){\rm{ }} = {\rm{ }}5635.08 + 183.5{\rm{A}}-86.54{\rm{B}}-\\&\quad 332.83{\rm{C}}-62.21{\rm{D}} + 3.87{\rm{AB}} + 11.50{\rm{AC}}-\\&\quad 1.25{\rm{AD}}-4.0{\rm{BC}} + 0.7{\rm{BD}} + 2.62{\rm{CD}}-\\&\quad 10.0{{\rm{A}}^2} + 2.5{{\rm{C}}^2} + 0.217{{\rm{D}}^2}.\end{split}$

2 -->

5.3.多目标优化求解结果

优化设计是在给定的设计区域和约束条件下求解目标值最优化的过程, 一般包括设计变量、设计目标以及约束条件三要素. 针对光声池的优化, 其品质因素及池常数在一定内范围内应尽可能大, 根据上述优化试验的设定, 光声池的优化数学模型可以表示为:

$\begin{split}&{\max\{ {f_1}\left( {{X}} \right), {f_2}\left( {{X}} \right)\} ;}\\&{{{X}} = {{[A, B, C, D]}^{\rm{T}}};}\\&{{\rm{s.t}}.}\\&{{f_3}\left( {{X}} \right) \geqslant 1000\;{\rm{Hz}};}\\&{5\;{\rm{mm }} \leqslant A \leqslant 7\;{\rm{mm }};8\;{\rm{mm }} \leqslant B \leqslant 12\;{\rm{mm }};}\\&{3\;{\rm{mm }} \leqslant C \leqslant 5\;{\rm{mm }};60\;{\rm{mm }} \leqslant D \leqslant 80\;{\rm{mm }};}\end{split}$

式中, X为设计变量, 表示一个4维列向量; f₁(X), f₂(X)为设计目标函数; f₃(X)为约束函数; 以上代号意义与前文相同.
光声池的优化是一个多目标优化问题, 通过上述响应面模拟结果可知光声池品质因素与池常数不能保证同时取到最优值, 设计目标之间相互冲突并且与设计变量呈现较为复杂的非线性关系. 本文采用多目标遗传算法MOGA (multi objective genetic algorithm)来解决该优化模型. 遗传算法是模拟生物在自然界进化过程形成的自适应的全局优化搜索算法, 通过群体的选择、交叉与变异等技术, 使群体进化到包含或接近最优解的状态. 优化模型采用Matlab中多目标遗传算法工具箱进行求解, 算法计算采用二进制编码, 遗传代数为 500, 种群数量为 60, 交叉率为0.8, 变异率为 0.01, 计算得到的优化目标Pareto 最优解集, 如图13所示.
从图13中可知, 光声池品质因数及其池常数两者之间是互为矛盾的, 品质因数取大时, 池常数就会较小, 反之亦然, 这使得光声池品质因数与池常数同时都取最大是不可能的. 根据Pareto最优解集, 考虑到改进后的光声池要优于之前性能指标, 对光声池两设计目标选择进行折衷处理, 从Pareto解集方案中选取3组进行考察, 同时希望光声池的品质因数与池常数相比较更加均衡增长, 因此选择方案2为最终设计方案并对参数值进行圆整, 即为方案4, 如表3所列, 初始结构参数值见本文3.2节.
图 13 Pareto最优解前沿分布图
Figure13. Pareto optimal solution frontier distribution map.

方案	rc/mm	hc/mm	Rc/mm	Lc/mm
1	5.96	11.99	3.55	61.43
2	5.6	11.99	4.78	60.28
3	7.0	11.99	5.0	60.0
4	5.6	12.0	4.7	60.0

表3优化后设计变量值
Table3.Optimized scheme value.

依据方案4所给出的结构参数值, 同时为了验证响应面代理模型的准确性, 对光声池进行模型重构与数值模拟计算, 结果如表4所列. 表4中的变化率以数值模拟结果为参照, 从结果中可知, 代理模型预测与数值模拟结果误差非常小, 光声池三种指标的最大误差仅为1.3%, 可见针对代理模型来进行优化计算的方法可行, 结果可靠; 在一定约束条件下, 优化后的光声池的品质因素Q较初始值增长了48.9%、池常数C_cell增长了34.4%, 光声池相关指标得到了一定改善.

指标	初始值	优化后			变化率/%
		代理 模型	数值 模拟	误差 率/%
Q	63.7	95.2	94.9	0.32	+ 48.9
Ccell/( Pa·cm) ·W–1	1750	2321	2353	1.3	+ 34.4
f/Hz	1648	2000	2002	0.1	+ 21.4

表4相关指标结果对比
Table4.Comparison of index results.

1)以声学方程与光谱吸收定律为基础, 建立了圆柱形光声池激发声场的数学模型, 确定了以工作纵向声学简正频率、品质因数和池常数为描述光声池特性的三项刻画指标及其解析公式;
2)利用数值模拟技术对光声池空腔结构进行了声学模态计算, 获得了光声池各阶声学模态值与其振型, 结合热黏性声学模拟方法对光声池的光热声效应进行了多物理场数值求解, 通过与解析计算和实验的对比, 明确了利用数值模拟方法来研究光声池特性的可行性;
3)提出了一种基于响应面代理模型与遗传算法相结合的光声池多目标优化求解方法, 探索将原光声池中的谐振腔两端形貌更改为喇叭口形的情况下, 通过求解获得了以光声池品质因数及池常数为最大值寻优的Pareto最优解集, 经过比照, 代理模型预测值与数值模拟值指标最大误差仅为1.3%, 优化后的新型光声池品质因素Q较之前增长了48.9%, 池常数C_cell增长了34.4%.