Size quantification of non-spherical bubbles by ultrasound

Zhang Ya-Jing; Li Fan; Lei Zhao-Kang; Wang Ming-Hao; Wang Cheng-Hui; Mo Run-Yang

doi:10.7498/aps.72.20222074

Abstract

Authors and contacts

FullText HTML

1. 引　言
2. 实验装置与测量方法
2.1 实验装置
2.2 基本理论
3. 结果与讨论
3.1 高速摄影-气泡形态和运动轨迹分析
3.2 运动气泡高频超声检测
4. 结　论
附录A　校准因子Cm的确定方法

References

Abstract

Ultrasonic detection is an effective method to quantify bubbles in opaque liquid, and acoustic scattering model is the key in ultrasonic inversion technique. Classical scattering models are usually based on the spherical assumption, and ka is much less than 1. However, these conditions are not always satisfied in practical applications. In this study, a quantitative strategy of ultrasonic inversion is proposed for non-spherical bubbles and ka deviation assumption. A series of solution models for a spherical gas bubble is established without considering the ka constraint, and it is compared with the classical Medwin (

$ka\ll1$ ) and Anderson (ka ≈ 1) models. The difference in scattering cross section σ_bs betweem them is only at the higher order formants of scattering, so the fitted line can be used to solve the multi-valued problem between σ_bs and ka. For a non-spherical bubble, σ_bs is determined by the frequency domain backscattering signal, the size is characterized by the equivalent radius a^*, and the inversion is performed by fitted curve from series solution model. Ultrasonic quantitative results are examined by high-speed photography. Results show that during the bubbles rising along a zigzag path, they develop non-spherical bubbles, their scattering cross sections are measured by the frequency domain scattering signal obtained at a position of ultrasonic measurement, and the equivalent radius is inverted by the series solution fitting curve. The deviation of the result from the actual result r₀ is about 1mm (relative error less than 45%) when 9≤kr₀≤35. This method can be used for implementing the acoustic inversion of non-spherical bubbles in a certain range of measurement accuracy.

Keywords:

PACS:

43.20.+g	(General linear acoustics)
43.30.+m	(Underwater sound)
43.58.+z.

Authors and contacts

Shaanxi Key Laboratory of Ultrasonics, Shaanxi Normal University, Xi’an 710119, China

Corresponding author: Mo Run-Yang, mmrryycn@snnu.edu.cn

Funds: Project supported by the National Natural Science Foundation of China (Grant Nos. 12074238, 11974232, 11774212)

FullText HTML

1. 引　言

液体中气泡在生物应用、渔业科学、工业选矿、纸浆生产及声空化等工程和环境应用中至关重要. 如在深水盆地和大陆边缘等地质环境中进行水下甲烷气体泄漏监测^[1,2], 工业选矿中利用气泡大小分布监测浮选过程^[3,4], 复合材料^[5,6]或推进剂^[7]生产中对气泡是否混入及气泡含量进行监测以提高质量; 造纸工序中O₂泡大小分布是决定纸浆脱木素环节是否顺利进行的关键^[8].

气泡表征常用方法有光学和声学两类. 光学方法如高速摄影^[9]、粒子图像测速系统^[10](particle image velocimetry, PIV)、Mie散射^[11]等, 定量精度高但在光线不充足或血液、深海及磁流体等不透明环境中则无能为力. 超声作为液体中最佳的信息载体, 由气泡引起的声能量衰减、声速变化、目标强度和散射截面等均与气泡尺寸和数量等信息直接相关. Ren等^[12]设计了脉冲超声和双光纤探头组合测量系统, 利用声衰减测量气泡直径, 并对三相流中气泡含量进行定量研究. Leighton等^[1]通过测量气泡引起的声速变化, 对水下沉积物表面气泡大小分布进行量化. Padilla等^[13]实验建立了水下气体通量与声背散射强度之间的关系, 对南加州海岸圣巴巴拉海峡油田天然气的通量进行评估, 并对该区天然气分布予以描绘.

利用声散射信息反演气泡大小的技术, 大都基于气泡散射解析模型. 较为经典的散射模型主要有Anderson模型、Medwin模型、Ainslie-Leighton 模型和Yeh模型^[14]等, 每个模型都是在一套特定假设条件下的数学表达. 散射不仅是入射声波频率和泡内外介质特性的函数, 而且要求泡形状为球形, 尺寸需满足特定假设条件, 若a为气泡半径, k为入射波的波数, Medwin模型^[15]假设 $ka\ll1$ , Anderson模型^[16]假设ka≈1. 然而, 这些基于气泡形状为球形的假设使得这些模型对非球形气泡难以应用. 事实上, 超声造影剂气泡半径通常在0.5—5 μm间, 形状通常为球形, 而液体介质中的气泡若其半径超过1 mm则很难保持球形^[17]. 海洋、工业过程最为常见的是半径在1—10 mm的非球形气泡, 而且用于气泡监测的声波频率一般从千赫兹到兆赫兹直接导致ka > 1^[14]. 气泡形状的非球形偏离及测量中ka偏离使得上述经典解析模型在用于气泡特性测量如半径反演时出现偏差. Padilla和Weber^[14]研究了非球形变对气泡背散射截面的影响, 并对基于Anderson, Medwin, Ainslie-Leighton和Yeh等经典模型进行对比, 发现ka在[0.03, 0.5]间模型与实验结果一致性很好, 而当ka在[0.5, 4.4]范围时理论与实验产生明显差异, 原因可能是由于这几种模型受ka远小于1或近似为1的假设条件约束.

本研究首先建立不受ka约束的球形气泡级数散射模型, 通过与Medwin ( $ka\ll1$ )和Anderson (ka≈1)模型对比, 分析ka偏离引发的散射截面效应; 实验测量水中上升的非球形气泡的背散射信号确定其散射截面, 通过散射等价对非球形气泡引入当量半径, 同时利用散射级数模型拟合关系对气泡当量半径进行反演. 反演结果与高速摄影图像分析结果进行对比, 研究非球形气泡超声定量方法.

2. 实验装置与测量方法

2.1 实验装置

以管内上升气泡为对象, 采用光学和声学两种方法同时进行. 高速摄影记录气泡运动过程并分析气泡尺寸、形变及轨迹等特征, 以这些数据作为对照组, 用以检验超声散射信号反演结果. 实验装置如图1所示, 充满水的圆柱形水槽高300 mm, 底面直径200 mm, 在其底部安插一个橡胶管, 管出口处接入一个微型针管并用以向水槽中注入气泡. 通过调节针管直径和气体流量等控制气泡大小和初速度, 气泡半径远远小于容器半径. 同一种针管产生的气泡在出口处具有相同的半径, 称为气泡出口半径或气泡初始半径并用r₀表示. 为研究半径不同气泡的运动, 采用9种不同规格的针管在当气体流量一定时, 获得出口半径r₀分别为0.5, 1.0, 1.3, 2.5, 3.2, 4.2, 4.5, 4.8, 5.5 mm的气泡.

$图 1 气泡测量装置(1-高速摄影机, 2-水槽, 3-流量控制阀, 4-进气口, 5-超声换能器, 6-声卡, 7-工控机)\r\nFig. 1. Bubble measuring device. 1-high-speed camera, 2-tank, 3-flow control valves, 4-air inlet, 5-ultrasonic transducer, 6-sound card, 7-controller.$

图 1 气泡测量装置(1-高速摄影机, 2-水槽, 3-流量控制阀, 4-进气口, 5-超声换能器, 6-声卡, 7-工控机)

Fig. 1. Bubble measuring device. 1-high-speed camera, 2-tank, 3-flow control valves, 4-air inlet, 5-ultrasonic transducer, 6-sound card, 7-controller.

下载: 全尺寸图片幻灯片

多功能超声信号发射接收卡 (型号: SUT20080125S, 中国科学院)产生电脉冲信号并激励固定在水槽壁上的8个超声换能器, 采集和存储气泡背散射信号 (采样频率为100 MHz). 以换能器声束中心线所在水平位置为超声测量位置, 以气泡离开管口时刻为计时起点. 高速摄影机 (型号: 千眼狼X213, 拍摄帧率: 1000帧/秒)安放在水槽左侧方, 记录气泡上升过程.

2.2 基本理论

2.2.1 气泡的声散射级数模型

声学反演技术的理论基础是散射模型. 球形气泡有多种散射模型, 其中以Medwin模型和Anderson模型最为常用. Anderson模型给出泡半径尺寸a与入射声波长λ相近 (ka≈1)的球形气泡散射模型, 未考虑热阻尼和液体黏性的影响, 仅考虑传播距离造成的声能衰减. Medwin模型综合考虑了传播距离、热阻尼和周围液体黏度对声能量的影响, 描述了半径远小于波长 (即 $ka\ll1$ )球形气泡的声散射. 然而, 这两个模型中球形假设以及ka≤1这些限制条件, 在很多实际应用需求中并不总是能得以满足. 那么, 若ka > 1以及非球形气泡的散射截面, 因偏离上述假设所产生的效应以及这些模型的偏离程度是值得研究的问题. 为此, 首先构建球形泡的声散射级数模型, 然后与上述两个模型对比, 利用偏离程度调整级数模型, 使其能用于非球形上升气泡半径的反演.

如图2所示, 球形气泡半径为a, 泡内气体为空气, 泡外的介质为水. 在幅值为1的平面波声场中, 为求解在泡外任意M点的散射声场, 借鉴Anderson^[16]模型中入射、散射及透射声波的表达式和边界条件, 在以r, θ和 $\varphi$ 为坐标变量的球坐标系中, 入射波p_i沿z轴方向, 省略时间因子, 则p_i展开为

$图 2 球形气泡散射模型\r\nFig. 2. Scattering model of spherical bubble.$

图 2 球形气泡散射模型

Fig. 2. Scattering model of spherical bubble.

下载: 全尺寸图片幻灯片

${p_{\text{i}}} = \sum\limits_{n = 0}^\infty {{{\rm{i}}^n}(2n + 1){{\text{j}}_n}({k_0}r)} {{\text{P}}_n}(\cos \theta ),$

(1)

其中k₀ = ω/c₀为水中波数, ω为入射声波角频率, c为水的密度, f为入射声波频率且ω = 2πf, n为阶数, j_n(k₀r)为球贝塞尔函数, P_n(cosθ)为勒让德函数.

泡外散射声压p_s可表示为

${p_{\text{s}}} = \sum\limits_{n = 0}^\infty {{{\rm{i}}^n}} (2n + 1){a_n}{\text{h}}_n^{(1)}({k_0}r){{\text{P}}_n}(\cos \theta ),$

(2)

其中, ${\text{h}}_n^{(1)}({k_0}r)$ 为第一类汉克尔函数, a_n为散射系数且可利用边界条件确定.

泡内的透射波声压p₁可表示为

${p_{\text{1}}} = \sum\limits_{n = 0}^\infty {{{\rm{i}}^n}} (2n + 1){b_n}{{\text{j}}_n}({k_1}r){{\text{P}}_n}(\cos \theta ),$

(3)

其中k₁= ω/c₁为空气中的波数, c₁为空气中声速, b_n为待定系数. 在r=a处满足边界条件:

${p_{\text{i}}} + {p_{\text{s}}} = {p_1},$

(4)

$\frac{{\rm{i}}}{{{\rho _0}{c_0}}}\frac{{\partial {p_{\text{i}}}}}{{\partial r}} + \frac{{\rm{i}}}{{{\rho _0}{c_0}}}\frac{{\partial {p_{\text{s}}}}}{{\partial r}} = \frac{{\rm{i}}}{{{\rho _1}{c_1}}}\frac{{\partial {p_1}}}{{\partial r}}.$

(5)

将(1)—(3)式代入(4)式和(5)式得散射系数a_n. 考虑到a_n与气泡形态函数 ${f_\infty }$ 存在关系:

${f_\infty } = \frac{1}{{{k_0}}}\sum\limits_{n = 0}^\infty {{{( - 1)}^n}(2n + 1){{\rm{i}}^n}} {{\text{e}}^{ - {\text{i}}(n + 1)/2}}{a_n}.$

散射截面σ_bs可表示为^[18]

${\sigma _{{\text{bs}}}} = {\left| {{f_\infty }} \right|^2}.$

(6)

需要注意, 利用级数法构建球形气泡散射模型时, 忽略热阻尼和液体黏度的影响, 因为当入射声波频率远高于气泡的共振频率时, 热阻尼和液体黏度对散射的影响很小^[19]. 其次, 构建级数模型时仅进行远场近似, 且未对ka进行任何限定, 故级数模型并不受ka限制.

2.2.2 ka偏离引发的散射截面效应

关注ka偏离假设后的效应, 将上述3个模型进行对比. 图3(a)和图3(b)给出ka在10^–3—1和1—50范围内无量纲散射截面(σ_bs/(πa²))的变化规律. 其中蓝色实线表示级数模型, Anderson模型和Medwin模型分别用红色虚线和绿色点划线表示. ka ≤ 0.1时 (图3(a)), 3条曲线完全重合且具有单值对应性, 此时利用散射截面σ_bs反演泡半径a时3个模型完全一致; 0.1 < ka ≤ 1时, 三者开始出现差异, 级数模型与Anderson模型仍保持高度一致性但Medwin模型产生偏离.

$图 3 气泡散射模型对比　(a) 0.001≤ka≤1; (b) 1<ka≤50\r\nFig. 3. Comparison of scattering models for bubbles: (a) 0.001≤ka≤1; (b) 1<ka≤50.$

图 3 气泡散射模型对比　(a) 0.001≤ka≤1; (b) 1<ka≤50

Fig. 3. Comparison of scattering models for bubbles: (a) 0.001≤ka≤1; (b) 1<ka≤50.

下载: 全尺寸图片幻灯片

ka在1—50间级数模型和Anderson模型高阶共振特点见图3(b). ka在1—10区间, 两模型曲线包括稀疏尖峰位置基本重合, 说明ka偏离对Anderson模型影响较小. 随ka增大, 两者高阶峰均更加密集且位置存在错落, ka偏离产生的散射截面效应主要体现在高阶共振上. ka > 1后高阶共振峰的出现意味着 (σ_bs/(πa²))-ka单调关系遭到破坏, 利用散射截面σ_bs反演泡半径a则无可能, 除非散射波干涉影响可忽略. 所幸的是, Sage等^[20]在实验中并未观察到高阶共振现象, 得到一条近似Anderson模型的 σ_bs/(πa²)-ka单值关系, 认为可忽略高阶共振继续采用Anderson模型进行反演. 事实上, 高阶共振未被实验观察到有两种可能: 一是散射模型均是基于球形假设, 不能分辨气泡的二维形状; 另一种可能是, 测量中数据点离散恰好与共振点错过. 然而, 当ka>10甚至更高为50时, 理论上高阶峰密集程度更高, 要使半径反演成为可能须对散射信号进行处理并设法消除高阶影响, 具体方法在实验部分3.2.2节详述. 另外, 气泡偏离球形假设产生的散射截面效应体现在依据σ_bs反演出的泡尺寸与实际大小的差距上, 同样将在3.2.2节中进行评价.

3. 结果与讨论

3.1 高速摄影-气泡形态和运动轨迹分析

3.1.1 上升气泡的形态变化

图4(a)—(c)分别为r₀ = 0.5, 1.3和2.5 mm的气泡从深度h = 120 mm水下上升时不同时刻的形态变化. 在实验过程中, 为保证气泡上升条件的一致性, 调节气体流量使气泡离开管口时的初速度极小几乎为零, 认为所有气泡均从静止开始上升. 由图4可见, 气泡刚离开管口时均为球形, 上升过程中形态发生了非球形变化, 对比发现, r₀越小形变越小. 通常认为半径在1 mm内的气泡基本保持球形形状^[17], 而较大的气泡则难以继续保持球形. 上升过程中气泡发生形变, 主要归因于泡内外压力差、黏滞阻力和浮力不断变化. 另外, 出口半径不同的气泡, 在上升运动路径中各自的力学过程不同, 到达同一水平位置时形态不同, 如在图4中超声测量位置处, 3个气泡在竖直方向上的形状完全不同.

$图 4 上升气泡的形态变化 (h = 120 mm)　(a) r0 = 0.5 mm; (b) r0 = 1.3 mm; (c) r0 = 2.5 mm\r\nFig. 4. Shape change of rising bubbles (h = 120 mm): (a) r0 = 0.5 mm; (b) r0 = 1.3 mm; (c) r0 = 2.5 mm.$

图 4 上升气泡的形态变化 (h = 120 mm)　(a) r₀ = 0.5 mm; (b) r₀ = 1.3 mm; (c) r₀ = 2.5 mm

Fig. 4. Shape change of rising bubbles (h = 120 mm): (a) r₀ = 0.5 mm; (b) r₀ = 1.3 mm; (c) r₀ = 2.5 mm.

下载: 全尺寸图片幻灯片

考虑到在测量位置处气泡形态不规则, 难以对其大小进行量化评定. 为此, 将这些非球形气泡近似看作为椭球形, 若该椭球与半径为r^*的球具有相同体积, 则气泡大小可用等效半径r^*表示^[21]:

${r^*} = \frac{1}{2}{(m{n^2})^{{1 \mathord{\left/ {\vphantom {1 3}} \right. } 3}}},$

(7)

其中, m和n分别为椭球的短轴和长轴, χ = n/m为气泡形变率. 通过高速摄影所得测量位置处气泡的长轴和短轴, 按(7)式对9种出口半径r₀不同气泡, 在测量位置处的等效半径r^*进行计算, 结果见表1.

表 1 测量位置处各泡的等效半径

$r^*$

Table 1. Equivalent radius r^* of each bubble at the measurement position.

	r₀/mm
	0.5	1.0	1.3	2.5	3.2	4.2	4.5	4.8	5.5
r^*/mm	0.5	1.1	1.4	2.4	3.2	4.6	4.4	4.9	5.9
χ	1.5	2.7	3.2	2.7	1.9	1.5	3.6	1.6	2.2

下载: 导出CSV

| 显示表格

从表1可以看出, r₀不同的气泡在上升过程中虽然大小形态变化各异, 却仍有共同的规律, 即到达测量位置时气泡的等效半径r^*略大于r₀, 个别意外可能是测量误差所致, 整体上两者差别非常小. 这符合力学规律, 表明可用等效半径r^*对不规则形状气泡进行描述.

3.1.2 气泡上升轨迹

高速摄影记录了气泡上升全过程, 并利用tracker软件分析其运动轨迹. 为便于对比, 将追踪软件所得每个泡的轨迹数据重构到一张轨迹图中. 以针管出口中心为坐标原点, 建立坐标系, 其中X表示气泡中心偏离原点O的水平位移, Y表示上升高度. 图5(a)和图5(b)分别为r₀ = 4.2和4.5 mm两个气泡, 从水深h = 120 mm (紫实线)和h = 160 mm (蓝点线)位置上升的轨迹. 由图5(a)和图5(b)可见, 气泡上升轨迹呈现“之”字形曲线, 水平偏移量与水深h和r₀有关. 除此之外, 水平偏移可能还有一定随机性. 图5(c)为r₀ = 4.2 mm两个等大气泡, 从同一深度h = 160 mm处上升的轨迹追踪结果. 开始时两泡轨迹基本重合, 但随着上升高度增大出现水平偏离, 两个等大、上升条件完全相同的气泡轨迹偏离表明随机性存在的可能.

$图 5 气泡“之”字形上升轨迹　(a) r0 = 4.2 mm, h = 120, 160 mm; (b) r0 = 4.5 mm, h = 120, 160 mm; (c) r0 = 4.2 mm, h = 160 mm\r\nFig. 5. Zigzag rising trajectory of the bubbles: (a) r0 = 4.2 mm, h = 120, 160 mm; (b) r0 = 4.5 mm, h = 120, 160 mm; (c) r0 = 4.2 mm, h = 160 mm.$

图 5 气泡“之”字形上升轨迹　(a) r₀ = 4.2 mm, h = 120, 160 mm; (b) r₀ = 4.5 mm, h = 120, 160 mm; (c) r₀ = 4.2 mm, h = 160 mm

Fig. 5. Zigzag rising trajectory of the bubbles: (a) r₀ = 4.2 mm, h = 120, 160 mm; (b) r₀ = 4.5 mm, h = 120, 160 mm; (c) r₀ = 4.2 mm, h = 160 mm.

下载: 全尺寸图片幻灯片

气泡上升轨迹发生水平偏移的原因与其形状变化密不可分, 而气泡形变则与浮力、黏滞阻力以及泡内外压力差和压差所致的内外气体交换等有关. 气泡上升的“之”字形轨迹以及水平偏移量与其出口半径r₀和水深h有关. 为进一步探究当h一定时水平偏移量与r₀间关系. 令h = 120 mm, 图6(a)—(c)为高速摄影直接追踪所获r₀ = 0.5, 1.3, 2.5 mm三个气泡的运动轨迹. 可见, r₀ = 0.5 mm气泡轨迹近似直线, 而后两者出现明显水平偏移. 可以推测, r₀在1 mm内的气泡上升过程中非球形变很小、轨迹近乎直线, 大泡则一定会发生非球形变.

$图 6 自同一水深处上升气泡轨迹 (h = 120 mm)　(a) r0 = 0.5 mm; (b) r0 = 1.3 mm; (c) r0 = 2.5 mm\r\nFig. 6. Rising bubble trajectory with same depth (h = 120 mm): (a) r0 = 0.5 mm; (b) r0 = 1.3 mm; (c) r0 = 2.5 mm.$

图 6 自同一水深处上升气泡轨迹 (h = 120 mm)　(a) r₀ = 0.5 mm; (b) r₀ = 1.3 mm; (c) r₀ = 2.5 mm

Fig. 6. Rising bubble trajectory with same depth (h = 120 mm): (a) r₀ = 0.5 mm; (b) r₀ = 1.3 mm; (c) r₀ = 2.5 mm.

下载: 全尺寸图片幻灯片

3.2 运动气泡高频超声检测

气泡在测量位置处的大小、形状与其运动过程密切相关. 气泡散射一般用散射截面σ_bs或远场背散射声压P_f表示, 利用测量位置处气泡的声散射特性, 可对其散射截面进行测量, 对形状及半径进行反演.

3.2.1 基于时域回波信号幅值的气泡截面形状轮廓

在图7(a)所示圆周每间隔45^°采集散射时域信号, 整个圆周共计采集8个信号. 图7(b)为在信号采集硬件、参数设置及其他条件完全相同时, r₀= 1.3 mm气泡在不同角度测量位置所获散射信号的时域波形图, 其中以红、蓝、绿和黑色表示的波形分别对应图7(a)中的1, 3, 5, 7号换能器, 因气泡形状为非球形, 在各测量角度所得回波幅值P_fj (j = 1, 2,···, 8)和回波时间t_j (j = 1, 2,···, 8)并不相同 (图7(b)中虚线对应不同位置处脉冲回波前沿), 利用回波到达时间和幅度差异对非球形气泡的水平截面进行描绘.

$图 7 气泡散射时域信号及采集方法　(a)信号采集示意图; (b)不同角度的气泡时域散射信号\r\nFig. 7. Time domain scattering signal of bubble and data acquisition method: (a) Schematic diagram of signal acquisition; (b) time domain scattering signals of bubbles at different angles.$

图 7 气泡散射时域信号及采集方法　(a)信号采集示意图; (b)不同角度的气泡时域散射信号

Fig. 7. Time domain scattering signal of bubble and data acquisition method: (a) Schematic diagram of signal acquisition; (b) time domain scattering signals of bubbles at different angles.

下载: 全尺寸图片幻灯片

勾画泡的水平截面形状时, 借用超声无损检测中当量半径概念, 即若一个非球形气泡的散射强度与同距离处球形气泡的散射强度相等, 则认为该球形泡的半径即为该非球形气泡的当量半径, 并用a^*表示. 换能器在8个角度采集到的气泡远场回波声压幅值P_fj与其当量半径 $a_j^*$ 间关系近似为^[22]

${P_{{\text{f}}j}} = A \times \frac{{a_j^*}}{{d_j^2}} , \quad A = \frac{{{P_0}{F_{\text{s}}}}}{{2\lambda }},$

(8)

其中d_j为泡表面到换能器端面的距离且d_j=(ct_j)/2, 水中声速c = 1480 m/s, P₀为入射波声压, F_s为压电晶片面积, λ为水中声波波长, A是与测量系统相关的系数且对同一套测试系统A为常数, 此时P_fj与a_j^*成线性关系.

具体步骤是, 首先将换能器在不同位置处采集到的回波信号幅值P_fj (j =1, 2, ···, 8)和回波时间t_j (j =1, 2,···, 8), 按关系式(8)求出其对应的等效半径a_j^*(j =1, 2, ···, 8), 然后分别以a_j^*为半径、以45°圆心角做8个圆弧 (由于散射信号是间隔45°采集, 故取45°为圆心角), 最后将这8段圆弧拼接得到非球形泡的水平截面形状. 当然, 若信号采集单元更多、圆心角分割更小, 则勾画出的水平截面图则更详尽. 表2第1行为出口半径r₀分别为0.5, 1.3和2.5 mm三个气泡在测量位置处的水平截面形状, 三者几乎都为椭圆. 为便于全面了解气泡形态, 表2第2行列出了高速摄影所获泡的竖直面图像. 为进一步分析气泡在上升过程中的形变, 将气泡在水平面上的形变率用χ^*表示, 上述3个气泡的χ^*依次为1.2, 1.2, 1.3, 对应于表1中竖直面上形变率χ分别为1.5, 3.2, 2.7. 可见, 气泡在上升过程中的形变主要发生在竖直方向上, 水平方向上的形变程度很小, 以至于气泡的水平面长短轴之比依然接近1, 更接近于圆形.

表 2 测量位置处气泡的形状及形变率

Table 2. Shape and deformation rate of bubbles at the measurement position.

	r₀/mm
	0.5	1.3	2.5
水平(χ^*)
竖直(χ)

下载: 导出CSV

| 显示表格

3.2.2 气泡散射截面σ_bs测量及当量半径a^*反演

用σ_bs描述气泡散射并通过频域信息对其进行测量. 具体步骤是: 采集并提取气泡散射时域信号, 对其进行傅里叶变换并将结果记为g₁(f ), 如图8(a), 提取某一频率所对应的g₁, 再按(9)式计算气泡散射截面σ_bs^[14]:

$图 8 频域散射信号(a)和散射截面校准因子(b)\r\nFig. 8. Frequency domain scattering signal (a) and calibration factor of scattering cross section (b).$

图 8 频域散射信号(a)和散射截面校准因子(b)

Fig. 8. Frequency domain scattering signal (a) and calibration factor of scattering cross section (b).

下载: 全尺寸图片幻灯片

${\sigma _{{\text{bs}}}}(f) = \frac{{{{\left| {{g_1}(f)} \right|}^2}{d^4}}}{{{C_{\text{m}}}(f)\exp ( - 4\alpha d)}},$

(9)

其中d为超声换能器端面到气泡表面的距离; α为水中的衰减系数; C_m^[23]为图8(b)所示校准因子(具体确定方法见附录A), 是一个与气泡大小无关而仅与收发系统、发射信号及泡周围介质有关的常量. 实验采用中心频率为4.5 MHz (–6 dB带宽, 2.5—6.5 MHz)的换能器激励和接收散射信号.

图9中为实验所测r₀=0.5, 1.3, 2.5, 3.2 mm四个气泡对应于各频率点上的散射截面σ_bs. 对比4条曲线发现, σ_bs具有频率依赖性; 相同频率下, 气泡r₀越大σ_bs越大、散射越强. 利用σ_bs进行a^*反演时不同频率下所得a^*值不同.

$图 9 –6 dB带宽范围σbs测量结果\r\nFig. 9. Measurements of σbs in the –6 dB bandwidth.$

图 9 –6 dB带宽范围σ_bs测量结果

Fig. 9. Measurements of σ_bs in the –6 dB bandwidth.

下载: 全尺寸图片幻灯片

为使反演结果更接近实际, 需根据气泡尺度合理选择频率或频率范围. 在此首先在探头–6 dB带宽范围 (2.5—6.5 MHz), 按照100 kHz间隔逐点对气泡当量半径a^*进行反演. 在此以f = 4.5 MHz为例介绍具体反演过程. 图10黑色点线为按(6)式计算级数模型当f = 4.5 MHz时的σ_bs-a^*曲线, 显然因高阶共振影响σ_bs与a^*并非单值对应, 在此对该曲线进行二次方拟合且拟合方程为 $\sigma_{\rm bs} = 3 \times 10^{-7}a^{*2} - 1.1\times 10^{-7}a^* + 5.8 \times 10^{-8}~\rm m^2$ (红色实线所示), 将超声所测σ_bs = 9.7×10^–8 m²与拟合曲线对照, 反演出其当量半径a^* = 0.6 mm(绿色所示).

$图 10 非球形气泡当量反演示意图\r\nFig. 10. Inversion schematic for non-spherical bubbles.$

图 10 非球形气泡当量反演示意图

Fig. 10. Inversion schematic for non-spherical bubbles.

下载: 全尺寸图片幻灯片

图11给出了上述4种气泡在探头–6 dB带宽内a^*的反演结果, 其中用δ = |a^*－r₀|/r₀表示其与高速摄影测量结果的相对误差. 对比r₀= 1.3, 2.5和3.2 mm三种气泡发现, r₀= 1.3 mm时a^*偏离r₀最小(用绿色表示, 18%≤δ≤42%), r₀= 2.5 mm次之(用红色表示, 35% ≤ δ ≤ 61%), r₀= 3.2 mm偏差最大(用蓝色表示, 32%≤δ≤67%). 由此可见, 气泡初始半径r₀越小、a^*偏离越小, 主要是由于反演的级数模型是建立在球形假设上, 气泡越小其形状越接近球形, 反之气泡尺寸越大非球形变越严重, 理论反演时偏差越大. r₀= 0.5 mm的气泡在9 ≤ kr₀≤ 14范围内的结果依然表现出和其他3种气泡相同的规律(用黄色表示), 但当5 ≤ kr₀＜9时却出现反常, 这可能是由于小气泡散射较弱使测量结果不准确所致, 而且对高速摄影图像处理得到的r₀与气泡实际大小也存在误差. 从kr₀范围上分析发现, 9 ≤ kr₀≤ 35范围内, 非球形气泡当量半径反演结果与高速摄影定量结果间的相对误差基本可控制在45%以内, 另外经计算发现实验中的4种气泡在这一范围里a^*和r₀的绝对误差|a^*－r₀|都基本可以控制在1 mm以内, 这种定量精度在很多工业应用中基本可以满足.

$图 11 超声与高速摄影定量结果对比\r\nFig. 11. Comparison of quantitative results of ultrasound and high-speed camera.$

图 11 超声与高速摄影定量结果对比

Fig. 11. Comparison of quantitative results of ultrasound and high-speed camera.

下载: 全尺寸图片幻灯片

4. 结　论

改进了偏离ka ≤ 1假设及发生非球形形变的气泡的散射反演技术, 给出了一种基于超声量化评定气泡形状和尺度的方法. 建立气泡的散射级数模型解除经典散射模型中ka ≤ 1这一约束条件, 同时通过引入当量半径概念对不规则的非球形气泡大小进行定量. 对于ka > 1时σ_bs-a曲线出现的高阶共振峰, 采用对级数解进行二次拟合的方法进行消除并将拟合曲线用于非球形气泡当量半径的反演. 在9 ≤ ka ≤ 35范围, 该方法在对实验所用出口半径分别为0.5, 1.3, 2.5, 3.2 mm气泡的测量中, 非球形气泡当量半径反演结果与高速摄影定量结果间的绝对误差基本可控制在1 mm以内, 相对误差可控制在45%以内. 此外, 还根据时域回波信号勾画出了非球形气泡的水平切面轮廓.

附录A　校准因子C_m的确定方法

校准因子C_m^[23]是一个与探头频率、发射信号类型及泡所处的液体介质有关的量, 和所测气泡半径并无直接关系, 因此在C_m确定中可任意选择某一尺度的气泡进行. 在此任选出口半径r₀= 1.0 mm气泡确定C_m(f ). 由于气泡上升过程中的形变和水平偏移具有随机性, 故通过对8组数据求平均的方式来降低这种随机性的影响.

将(9)式变形可得 $\small{C_{{\text{m}}i}}(f) = \dfrac{{{{\left| {{g_{2 i}}(f)} \right|}^2}{d^4}}}{{\sigma _{{\text{bs}}}^{\prime}(f)\exp ( - 4\alpha d)}}$ , 其中g_2i(f )是r₀= 1.0 mm气泡1组时域信号的傅里叶变换结果, $\sigma _{{\text{bs}}}^{\prime}(f)$ 是由(6)式所得级数散射截面. 重复处理r₀ = 1.0 mm气泡的8组散射信号得到了8组C_mi(i = 1, 2,···, 8), 最终本实验环境下的校准因子为 ${C_{\text{m}}}(f) = \dfrac{1}{8}\displaystyle\sum\limits_{i = 1}^8 {{C_{{\text{m}}i}}(f)}$ , 见图8(b)所示. 在C_m(f )的基础上, 取待测气泡的散射回波, 按(9)式计算σ_bs(f ), 即可通过回波频域求得待测气泡的散射截面, 在此过程中同样需要考虑到随机性带来的影响, 因而仍然取8组数据的平均值作为该气泡最终的背散射截面测量结果.

References

[1]	Leighton T G, Dogan H, Fox P, Mantouka A, Best A I, Robb G B R, White P R 2021 J. Acoust. Soc. Am. 150 2705 Google Scholar
[2]	Judd A G, Hovland M, Dimitrov L I, García Gil S, Jukes V 2002 Geofluids 2 109 Google Scholar
[3]	Kracht W, Moraga C 2016 Miner. Eng. 98 122 Google Scholar
[4]	Liu J, Gao Q, Tang Z, Xie Y, Gui W, Ma T, Niyoyita J P 2020 IEEE Trans. Instrum. Meas. 69 9618 Google Scholar
[5]	Guédra M, Inserra C, Gilles B 2017 Ultrason. Sonochem. 38 298 Google Scholar
[6]	Buckey J C, Knaus D A, Alvarenga D L, Kenton M A, Magari P J 2005 Acta Astronaut. 56 1041 Google Scholar
[7]	Wen W, Zong G H, Bi S S 2014 Rev. Sci. Instrum. 85 065106 Google Scholar
[8]	Jarmo I, Tuomas E, Heikki M, Lasse L, Jari K, Heikki K 2014 19th Iberoamerican Congress Puerto Vallarta, Mexico Puerto Vallarta, Mexico, November 2−5, 2014 p38
[9]	Ilonen J, Juránek R, Eerola T, Lensu L, Dubská M, Zemčík P, Kälviäinen H 2018 Pattern Recognit. Lett. 101 60 Google Scholar
[10]	Ohta J, Doyama N, Wakabayashi D, Suzuki H 2013 T. Jpn. Soc. Mech. Eng. Part B 79 2397 Google Scholar
[11]	Bradley P B, Seth J P 1992 Phys. Rev. Lett. 69 03839 Google Scholar
[12]	Ren W, Jin N, Zhang J 2022 Ultrasonics 124 106740 Google Scholar
[13]	Padilla A M, Loranger S, Kinnaman F S, Valentine D L, Weber T C 2019 J. Geophys. Res. Ocean. 124 2472 Google Scholar
[14]	Padilla A M, Weber T C 2021 J. Acoust. Soc. Am. 149 2504 Google Scholar
[15]	Clay C S, Medwin H 1977 Acoustical Oceanography: Principles and Applications (New York: Wiley) pp461–466
[16]	Anderson V C 1950 J. Acoust. Soc. Am. 22 426 Google Scholar
[17]	Zheng B L, Poojitha D Y M 2000 J. Hydraul. Eng. 126 852 Google Scholar
[18]	Dong X, Su M, Cai X 2012 Particuology 10 117 Google Scholar
[19]	Spiekhout S, Voorneveld J, van Elburg B, Renaud G, Segers T, Lajoinie G P R, Versluis M, Verweij M D, de Jong N, Bosch J G 2022 J. Acoust. Soc. Am. 151 3993 Google Scholar
[20]	Sage K A, George J, Oberall H 1979 J. Acoust. Soc. Am. 65 1413 Google Scholar
[21]	Sam A, Gomez C O, Finch J A 1996 Int. J. Miner. Process. 47 177 Google Scholar
[22]	郑晖, 林树青 2008 超声检测 (北京: 中国劳动社会保障社) 第75页 Zheng H, Lin S Q 2008 Ultrasonic Inspection (Beijing: China Labor and Social Security Publishing House) p75 (in Chinese)
[23]	Weber T C, Ward L G 2015 J. Acoust. Soc. Am. 138 2169 Google Scholar

图 1 气泡测量装置(1-高速摄影机, 2-水槽, 3-流量控制阀, 4-进气口, 5-超声换能器, 6-声卡, 7-工控机)

Figure 1. Bubble measuring device. 1-high-speed camera, 2-tank, 3-flow control valves, 4-air inlet, 5-ultrasonic transducer, 6-sound card, 7-controller.

DownLoad: Full-Size Img PowerPoint

图 2 球形气泡散射模型

Figure 2. Scattering model of spherical bubble.

DownLoad: Full-Size Img PowerPoint

图 3 气泡散射模型对比　(a) 0.001≤ka≤1; (b) 1<ka≤50

Figure 3. Comparison of scattering models for bubbles: (a) 0.001≤ka≤1; (b) 1<ka≤50.

DownLoad: Full-Size Img PowerPoint

图 4 上升气泡的形态变化 (h = 120 mm)　(a) r₀ = 0.5 mm; (b) r₀ = 1.3 mm; (c) r₀ = 2.5 mm

Figure 4. Shape change of rising bubbles (h = 120 mm): (a) r₀ = 0.5 mm; (b) r₀ = 1.3 mm; (c) r₀ = 2.5 mm.

DownLoad: Full-Size Img PowerPoint

图 5 气泡“之”字形上升轨迹　(a) r₀ = 4.2 mm, h = 120, 160 mm; (b) r₀ = 4.5 mm, h = 120, 160 mm; (c) r₀ = 4.2 mm, h = 160 mm

Figure 5. Zigzag rising trajectory of the bubbles: (a) r₀ = 4.2 mm, h = 120, 160 mm; (b) r₀ = 4.5 mm, h = 120, 160 mm; (c) r₀ = 4.2 mm, h = 160 mm.

DownLoad: Full-Size Img PowerPoint

图 6 自同一水深处上升气泡轨迹 (h = 120 mm)　(a) r₀ = 0.5 mm; (b) r₀ = 1.3 mm; (c) r₀ = 2.5 mm

Figure 6. Rising bubble trajectory with same depth (h = 120 mm): (a) r₀ = 0.5 mm; (b) r₀ = 1.3 mm; (c) r₀ = 2.5 mm.

DownLoad: Full-Size Img PowerPoint

图 7 气泡散射时域信号及采集方法　(a)信号采集示意图; (b)不同角度的气泡时域散射信号

Figure 7. Time domain scattering signal of bubble and data acquisition method: (a) Schematic diagram of signal acquisition; (b) time domain scattering signals of bubbles at different angles.

DownLoad: Full-Size Img PowerPoint

图 8 频域散射信号(a)和散射截面校准因子(b)

Figure 8. Frequency domain scattering signal (a) and calibration factor of scattering cross section (b).

DownLoad: Full-Size Img PowerPoint

图 9 –6 dB带宽范围σ_bs测量结果

Figure 9. Measurements of σ_bs in the –6 dB bandwidth.

DownLoad: Full-Size Img PowerPoint

图 10 非球形气泡当量反演示意图

Figure 10. Inversion schematic for non-spherical bubbles.

DownLoad: Full-Size Img PowerPoint

图 11 超声与高速摄影定量结果对比

Figure 11. Comparison of quantitative results of ultrasound and high-speed camera.

DownLoad: Full-Size Img PowerPoint

表 1 测量位置处各泡的等效半径 $r^*$

Table 1. Equivalent radius r^* of each bubble at the measurement position.

r₀/mm
0.5 1.0 1.3 2.5 3.2 4.2 4.5 4.8 5.5

r^*/mm 0.5 1.1 1.4 2.4 3.2 4.6 4.4 4.9 5.9
χ 1.5 2.7 3.2 2.7 1.9 1.5 3.6 1.6 2.2

DownLoad: CSV

表 2 测量位置处气泡的形状及形变率

Table 2. Shape and deformation rate of bubbles at the measurement position.

r₀/mm
0.5 1.3 2.5

水平(χ^*)
竖直(χ)

DownLoad: CSV

[1]	Leighton T G, Dogan H, Fox P, Mantouka A, Best A I, Robb G B R, White P R 2021 J. Acoust. Soc. Am. 150 2705 Google Scholar
[2]	Judd A G, Hovland M, Dimitrov L I, García Gil S, Jukes V 2002 Geofluids 2 109 Google Scholar
[3]	Kracht W, Moraga C 2016 Miner. Eng. 98 122 Google Scholar
[4]	Liu J, Gao Q, Tang Z, Xie Y, Gui W, Ma T, Niyoyita J P 2020 IEEE Trans. Instrum. Meas. 69 9618 Google Scholar
[5]	Guédra M, Inserra C, Gilles B 2017 Ultrason. Sonochem. 38 298 Google Scholar
[6]	Buckey J C, Knaus D A, Alvarenga D L, Kenton M A, Magari P J 2005 Acta Astronaut. 56 1041 Google Scholar
[7]	Wen W, Zong G H, Bi S S 2014 Rev. Sci. Instrum. 85 065106 Google Scholar
[8]	Jarmo I, Tuomas E, Heikki M, Lasse L, Jari K, Heikki K 2014 19th Iberoamerican Congress Puerto Vallarta, Mexico Puerto Vallarta, Mexico, November 2−5, 2014 p38
[9]	Ilonen J, Juránek R, Eerola T, Lensu L, Dubská M, Zemčík P, Kälviäinen H 2018 Pattern Recognit. Lett. 101 60 Google Scholar
[10]	Ohta J, Doyama N, Wakabayashi D, Suzuki H 2013 T. Jpn. Soc. Mech. Eng. Part B 79 2397 Google Scholar
[11]	Bradley P B, Seth J P 1992 Phys. Rev. Lett. 69 03839 Google Scholar
[12]	Ren W, Jin N, Zhang J 2022 Ultrasonics 124 106740 Google Scholar
[13]	Padilla A M, Loranger S, Kinnaman F S, Valentine D L, Weber T C 2019 J. Geophys. Res. Ocean. 124 2472 Google Scholar
[14]	Padilla A M, Weber T C 2021 J. Acoust. Soc. Am. 149 2504 Google Scholar
[15]	Clay C S, Medwin H 1977 Acoustical Oceanography: Principles and Applications (New York: Wiley) pp461–466
[16]	Anderson V C 1950 J. Acoust. Soc. Am. 22 426 Google Scholar
[17]	Zheng B L, Poojitha D Y M 2000 J. Hydraul. Eng. 126 852 Google Scholar
[18]	Dong X, Su M, Cai X 2012 Particuology 10 117 Google Scholar
[19]	Spiekhout S, Voorneveld J, van Elburg B, Renaud G, Segers T, Lajoinie G P R, Versluis M, Verweij M D, de Jong N, Bosch J G 2022 J. Acoust. Soc. Am. 151 3993 Google Scholar
[20]	Sage K A, George J, Oberall H 1979 J. Acoust. Soc. Am. 65 1413 Google Scholar
[21]	Sam A, Gomez C O, Finch J A 1996 Int. J. Miner. Process. 47 177 Google Scholar
[22]	郑晖, 林树青 2008 超声检测 (北京: 中国劳动社会保障社) 第75页 Zheng H, Lin S Q 2008 Ultrasonic Inspection (Beijing: China Labor and Social Security Publishing House) p75 (in Chinese)
[23]	Weber T C, Ward L G 2015 J. Acoust. Soc. Am. 138 2169 Google Scholar

[1]	Li Jiong-Yuan, Meng Ju, Wang Ke-Dong. Low-energy electron elastic scattering of ${\mathbf{C}}_{4}^{-}$ anions: Resonance states and conformers. Acta Physica Sinica, 2024, 73(24): 243401. doi: 10.7498/aps.73.20241377
[2]	Liu Yu, Tian Qiang, Wang Xin-Yan, Guan Xue-Fei. Efficient grain size evaluation based on single direction measurement of ultrasonic backscattering coefficient. Acta Physica Sinica, 2024, 73(7): 074301. doi: 10.7498/aps.73.20231959
[3]	Zhang Ya-Jing, Wang Ming-Hao, Lei Zhao-Kang, Shen Wen-Jie, Ma Yan-Qiang, Mo Run-Yang. Acoustic scattering properties of multilayer membrane structured magnetic microbubbles. Acta Physica Sinica, 2022, 71(18): 184302. doi: 10.7498/aps.71.20220847
[4]	Wu Xue-You, Liang Jin-Fu. Translation and nonspherical oscillation of single bubble in ultrasound field. Acta Physica Sinica, 2021, 70(18): 184301. doi: 10.7498/aps.70.20210513
[5]	Ma Yan, Lin Shu-Yu, Xu Jie, Tang Yi-Fan. Influnece of nonspherical effects on the secondary Bjerknes force in a strong acoustic field. Acta Physica Sinica, 2017, 66(1): 014302. doi: 10.7498/aps.66.014302
[6]	Di Ma-Ke, Shen Guang-Xian, Zhao Yun-Qiang, Zeng Ruo-Sheng, Wang Rong-Kai. Theoretical calculation of the vib-rotational interaction potential and the scattering cross section for the Ar-H2 (D2, T2) collision system. Acta Physica Sinica, 2015, 64(13): 133101. doi: 10.7498/aps.64.133101
[7]	Xu Shi-Long, Hu Yi-Hua, Zhao Nan-Xiang, Wang Yang-Yang, Li Le, Guo Li-Ren. Impact of metal target’s atom lattice structure on its quantum radar cross-section. Acta Physica Sinica, 2015, 64(15): 154203. doi: 10.7498/aps.64.154203
[8]	Han Jin-Zhong, Qin Zhen, Wang Xue-Lei. Joint production of Z boson and charged top-pion pair at the ILC. Acta Physica Sinica, 2012, 61(4): 041201. doi: 10.7498/aps.61.041201
[9]	Hu Ya-Hua, Ye Dan-Dan, Qi Yue-Ying, Liu Xiao-Ju, Liu Ling. Research on collision ionization process of Be atom by proton. Acta Physica Sinica, 2012, 61(24): 243401. doi: 10.7498/aps.61.243401
[10]	Gong Ming-Yan. Rotational excitation of He-BH collision system. Acta Physica Sinica, 2011, 60(7): 073401. doi: 10.7498/aps.60.073401
[11]	Dai Bing, Luo Xiang-Dong, Wang Ya-Wei. Multiple light scattering of non-spherical particles with elliptical cross section. Acta Physica Sinica, 2009, 58(6): 3864-3869. doi: 10.7498/aps.58.3864
[12]	Huang Jin-Shu, Luo Peng-Hui, Lu Gong-Ru. Bottom quark pair production in γγ collision. Acta Physica Sinica, 2009, 58(12): 8166-8173. doi: 10.7498/aps.58.8166
[13]	Shen Guang-Xian, Wang Rong-Kai, Linghu Rong-Feng, Yang Xiang-Dong. Study of isotope effect in 3He(4He)-H2 collisions. Acta Physica Sinica, 2008, 57(6): 3452-3457. doi: 10.7498/aps.57.3452
[14]	Shen Guang-Xian, Wang Rong-Kai, Linghu Rong-Feng, Yang Xiang-Dong. Theoretical calculation of the partial cross section in He-H2(D2,T2) collisions. Acta Physica Sinica, 2008, 57(1): 155-159. doi: 10.7498/aps.57.155
[15]	Li Ying-Le, Huang Ji-Ying, Wang Ming-Jun. Investigation of electromagnetic complex scattering for spherical targets in Ku wave band. Acta Physica Sinica, 2008, 57(12): 7630-7634. doi: 10.7498/aps.57.7630
[16]	Yu Chun-Ri, Song Xiao-Shu, Cheng Xin-Lu, Yang Xiang-Dong, Shen Chuan-Sheng. The interaction potential of Ne-HF system and close coupling calculations of the scattering cross sections. Acta Physica Sinica, 2008, 57(6): 3446-3451. doi: 10.7498/aps.57.3446
[17]	Wang Rong-Kai, Shen Guang-Xian, Yu Chun-Ri, Yang Xiang-Dong. Theoretical calculation of the scattering cross section for He-HF(DF,TF) collision system. Acta Physica Sinica, 2008, 57(11): 6932-6938. doi: 10.7498/aps.57.6932
[18]	Li Min-Quan, Tao Xiao-Jun, Zhao Jin, Wu Xian-Liang. Radar cross section computation using symplectic Runge-Kutta-Nystrom method. Acta Physica Sinica, 2007, 56(4): 2115-2118. doi: 10.7498/aps.56.2115
[19]	Yu Chun-Ri, Wang Rong-Kai, Cheng Xin-Lu, Yang Xiang-Dong. Theoretical study of the effect of potential models on scattering cross sections for He-HF system. Acta Physica Sinica, 2007, 56(5): 2577-2584. doi: 10.7498/aps.56.2577
[20]	SHEN HAO, CHENG HUAN-SHENG, TANG JIA-YONG, YANG FU-JIA. THE STUDY OF NON-RUTHERFORD BACKSCATTERING OF 4He FROM CARSON. Acta Physica Sinica, 1994, 43(10): 1569-1575. doi: 10.7498/aps.43.1569

Catalog

Vol. 72, Issue 3 - 2023-02-05

Metrics

Abstract views: 4488
PDF Downloads: 93

姓名
邮箱
手机号码
标题
留言内容
验证码

Search

Article

留言板