气泡与空化——Rayleigh 塌缩、Minnaert 共振与枪虾
2026 年 10 月
1 导读
一个气泡只有一个自由度——半径 \(R(t)\),却能把大气压在毫米尺度上做的功集中成上千个大气压的压强脉冲和数千开尔文以上的高温(在声致发光中可达 1 × 104 K 量级),也能像一个“液体质量挂在气体弹簧上”的振子那样唱出清脆的音调。本讲用普通物理的力学与热学(牛顿定律、能量守恒、理想气体、表面张力)推出描述气泡的 Rayleigh–Plesset 方程,再沿着它展开三条主线:塌缩、振荡、以及大量气泡共同改变液体的声速。
19 世纪 90 年代高速轮船出现后,螺旋桨的推力在高转速下出现异常下降。英国驱逐舰 HMS Daring(1893 年下水)1894 年试航时达不到设计航速,Thornycroft 与 Barnaby 把原因归于桨叶背面形成的空腔,并在 1895 年的报告中讨论了这一现象 [1], [2];“空化”(cavitation)一词由 R. E. Froude 提出,并由 Barnaby 与 Thornycroft 于 1895 年首次在文献中使用 [3]。Parsons 在为他的汽轮机快艇 Turbinia 设计螺旋桨时遇到同样的问题,并于 1895 年建造了第一座研究空化的水洞 [3], [4]。人们还发现空泡塌缩会在桨叶上“啃”出麻点(空蚀)。1917 年 Rayleigh 发表短文《论液体中球形空腔塌缩时产生的压强》:他在文中写道,读到 O. Reynolds 对水壶将沸时声音的描述后开始思考这一问题,并从 Parsons 处得知后者也因螺旋桨后的空化关心同一问题,而 S. Cook 曾在 Parsons 的建议下、基于 Besant(1859)书中的一道题计算过塌缩被突然阻止时的压强;Rayleigh 给出了更简单的推导,求出塌缩时间 \(0.915R_0\sqrt{\rho/P}\)(\(P\) 为远处压强,即本讲的 \(p_\infty\))与液体中的压强分布,指出 \(R=R_0/20\) 时液体中压强可达外压的 1260 倍 [5], [6]。这是本讲第 5 节的内容。
1933 年,比利时出生、在荷兰乌得勒支工作的天文学家 Minnaert 发表《论会唱歌的气泡与流水的声音》,用“液体惯性 + 气体弹性”解释了气泡的共振频率,他的数据要求气体按绝热(\(\gamma\approx1.4\))而非等温压缩 [7], [8](第 6 节);据 Crawford 转述,他还用其他气体检验过 \(\gamma\) 的作用 [8]。半个世纪后 Leighton 与 Walton 用氦、丙烷气泡重做了这一检验,测得 \(\kappa=1.46\pm0.09\)(氦,\(\gamma=5/3\))与 \(1.10\pm0.02\)(丙烷,\(\gamma=1.13\)) [9]。1933 年 Marinesco 与 Trillat 发现超声浴中的照相底片会感光;1934 年 Frenzel 与 Schultes 证明这是超声辐照的水在发光,即多泡声致发光 [10], [11], [12]。1990 年前后,Gaitan 等人用驻波声场把单个气泡稳定地悬浮在烧瓶中心(声压 1.1 atm–1.5 atm,纯水中为 1.1 atm–1.3 atm,较高值对应甘油水溶液;频率 21 kHz–25 kHz),每个声周期塌缩一次、发一次光,开创了单泡声致发光研究 [12], [13], [14]。
枪虾(鼓虾)的“枪声”长期被认为是螯的两指相撞;2000 年 Versluis 等用高速摄影证明它来自螯射出的高速水流造成的空化泡的塌缩,2001 年又观察到塌缩时的闪光 [15], [16](第 8.3 节)。海面下的雨声在约 15 kHz 附近有一个特征峰,来自小雨滴卷入的气泡的 Minnaert 振荡 [17], [18]。医学超声造影剂是微米级的包膜微泡,其共振频率恰在 MHz 超声频段 [19]。IPhO 2025(巴黎)理论第 3 题“Champagne!”分三部分:A 气泡的成核、长大与上升;B 气泡在液面破裂时的发声(建模为 Helmholtz 共鸣器——一段“瓶颈”中的液体质量对泡内气体弹簧振荡,可与本讲泡壁整体胀缩的 Minnaert 模型对比);C 开瓶时瓶塞下 CO\(_2\) 的绝热膨胀。可作为本讲的延伸阅读 [20]。
学习目标
从质量守恒与牛顿定律出发,导出球对称流的速度场、液体动能和非定常伯努利方程,并用能量方法与伯努利方法两种方式推出 Rayleigh–Plesset 方程。
掌握 Rayleigh 塌缩的能量积分、塌缩时间(Beta 函数)、速度奇点、液体中的压强峰,以及含气泡反弹的能量估算和它的局限。
把线性化 RP 方程理解为“质量–弹簧”振子:Minnaert 频率、有效质量、等温/绝热判据、三种阻尼与 \(Q\) 值。
推导 Wood 公式,解释含泡液体声速的极小值与热巧克力效应。
会用空化数、Blake 阈值与 Rayleigh 时间对枪虾、螺旋桨等真实空化现象做自洽的量级估算;会处理“由气泡声音测多方指数”的实验数据。
已掌握:牛顿定律与动能定理、理想气体状态方程与绝热过程、表面张力的定义、简谐振动与阻尼振动、平面声波、微积分(链式法则、分离变量、简单常微分方程)、复数。
本讲补充(第 3 节与附录):流体元的运动方程(欧拉方程)、不可压缩流的连续性、速度势与非定常伯努利方程、拉普拉斯压强、多方过程与热扩散长度、法向粘性应力、球面声波与单极子辐射、Beta 函数、泡内热传导。
建议学时
约 2 学时:前置模块与 RP 方程 30 分钟,Rayleigh 塌缩 25 分钟,Minnaert 共振 25 分钟,Wood 公式与空化/枪虾 20 分钟,例题与实验题选讲 20 分钟。附录供课后阅读。
本讲的逻辑主线
质量守恒 \(\Rightarrow u=R^2\dot R/r^2\) \(\Rightarrow\) 动能 \(2\pi\rho R^3\dot R^2\) / 伯努利方程 \(\Rightarrow\) RP 方程 \(\Rightarrow\) { 塌缩(大振幅);共振(小振幅);声辐射;含泡声速 }
2 名词与符号
| 符号 | 含义 | 单位 / 典型值 |
|---|---|---|
| \(R(t)\), \(\dot R\), \(\ddot R\) | 泡半径及其一、二阶时间导数 | m, m/s, m/s\(^2\) |
| \(R_0\) | 初始(最大)半径或平衡半径(依上下文) | m |
| \(r\) | 到泡心的距离 | m |
| \(u(r,t)\) | 液体径向速度 | m/s |
| \(\phi(r,t)\) | 速度势,\(u=\partial\phi/\partial r\) | m\(^2\)/s |
| \(\rho\) | 液体密度 | 998 kg/m3 |
| \(p_\infty(t)\), \(p_0\) | 远场压强;静压 | Pa;1.013 × 105 Pa |
| \(p_B\), \(p_g\), \(p_v\) | 泡内总压强、气体分压、饱和蒸气压 | Pa;\(p_v=2.3\,\mathrm{kPa}\) |
| \(p_L\) | 紧贴泡壁的液体压强 | Pa |
| \(\sigma\) | 表面张力系数 | 0.0728 N/m |
| \(\mu\) | 液体(动力)粘度 | 1.0 × 10−3 Pa s |
| \(\kappa\), \(\gamma\) | 多方指数;比热比 | 空气 \(\gamma=1.40\) |
| \(D\) | 气体热扩散率 \(k_{\mathrm{th}}/(\rho_gc_p)\) | 2.12 × 10−5 m2/s |
| \(c\) | 液体声速 | 1482 m/s |
| \(K\) | 液体动能 | J |
| \(m_{\mathrm{eff}}\), \(k_s\) | 泡壁振荡的有效质量、劲度 | kg, N/m |
| \(\tau\), \(U_0\) | 塌缩的时间、速度单位 \(R_0\sqrt{\rho/p_\infty}\)、\(\sqrt{p_\infty/\rho}\) | s, m/s |
| \(t_c\) | Rayleigh 塌缩时间 \(0.915\tau\) | s |
| \(z\) | 体积比 \(R_0^3/R^3\) | 1 |
| \(\omega_0\), \(f_0\) | Minnaert(角)频率 | s\(^{-1}\), Hz |
| \(\delta\), \(Q\), \(\Gamma\) | 无量纲阻尼、品质因数 \(Q=1/\delta\)、振幅衰减率 \(\Gamma\) | 1, 1, s\(^{-1}\) |
| \(\mathrm{Pe}\) | 热 Péclet 数 \(\omega R^2/D\) | 1 |
| \(\phi\)(第 7 节) | 气体体积分数(空泡率) | 1 |
| \(c_m\), \(\rho_m\), \(\beta\) | 混合物声速、密度;压缩率 | m/s, kg/m\(^3\), Pa\(^{-1}\) |
| \(\sigma_{\mathrm c}\), \(U\) | 空化数 \((p_\infty-p_v)/(\frac12\rho U^2)\);来流速度 | 1, m/s |
| \(p_{\mathrm{Bl}}\) | Blake 阈值(泡核失稳的临界液体压强) | Pa |
注意:速度势 \(\phi\) 与空泡率 \(\phi\) 用同一字母,前者只出现在第 3–5 节、第 8.2 节的绕球小练习与附录,后者只出现在第 7 节及例题 3;表面张力 \(\sigma\) 与空化数 \(\sigma_{\mathrm c}\) 以下标区分;泡内压强 \(p_B\) 与 Blake 阈值 \(p_{\mathrm{Bl}}\) 不同;\(\beta\) 只表示压缩率,阻尼振子的振幅衰减率记为 \(\Gamma\)(附录 B 的 \(\Gamma(\cdot)\) 是带自变量的 Gamma 函数,不会混淆);\(k_s\)、\(k_{\mathrm{th}}\)、\(k=\omega/c\) 分别是劲度、热导率与声波波数。
3 前置模块:描述气泡所需的流体力学与热学工具
本节把后面要用到、但普通物理课程通常没有系统讲过的工具逐一补齐。每个工具都只推导到“够用”为止,并且只针对本讲的情形——一个被无界液体包围的球形气泡(图 1)。
3.1 流体元的牛顿第二定律:压强梯度力与“随体加速度”
压强梯度力
取液体中一个边长为 \(\dd x,\dd y,\dd z\) 的小长方体。左面受力 \(p(x)\,\dd y\dd z\) 指向 \(+x\),右面受力 \(p(x+\dd x)\,\dd y\dd z\) 指向 \(-x\),合力 \[\tag{3.1} \dd F_x=\bigl[p(x)-p(x+\dd x)\bigr]\dd y\,\dd z=-\pdv{p}{x}\,\dd V .\] 即单位体积液体受到的压强合力为 \(-\pdv{p}{x}\)(三维时为 \(-\nabla p\))。在静止液体中它与重力平衡,\(\dv{p}{z}=-\rho g\),这正是普通物理里的静水压强公式。对球对称的压强场 \(p(r,t)\),在任一点把 \(x\) 轴取为该点的径向,式 (3.1) 就给出单位体积的径向力 \(-\pdv{p}{r}\)。
随体加速度
流体力学中速度用“场”描述:\(u(r,t)\) 是 \(t\) 时刻位于 \(r\) 处的那份液体的速度。牛顿定律要用的是同一份液体的加速度。设这份液体的位置为 \(r(t)\),则 \(\dot r=u(r(t),t)\),由链式法则 \[\tag{3.2} a=\dv{}{t}u\bigl(r(t),t\bigr)=\pdv{u}{t}+\dv{r}{t}\pdv{u}{r}=\pdv{u}{t}+u\pdv{u}{r} .\] 第一项是“同一地点速度随时间的变化”,第二项是“液体搬到了速度不同的地方”带来的变化(称为对流项)。纯径向流动中液体沿径线运动、方向不变,所以这个一维链式法则是精确的。
忽略粘性与重力,单位体积液体的牛顿第二定律为 \[\tag{3.3} \rho\left(\pdv{u}{t}+u\pdv{u}{r}\right)=-\pdv{p}{r}.\] 这是欧拉方程(无粘流体的运动方程)在球对称情形下的形式。粘性力在本讲的流动中体内恰好为零(附录 A),只在泡壁边界上起作用,见第 3.7 节。
重力的作用分两部分:一是建立静水压强背景 \(p_{\mathrm{atm}}+\rho g h\),我们把它直接计入“远处压强” \(p_\infty\);二是使气泡上浮。对 \(R=1\,\mathrm{mm}\) 的泡,泡上下的静压差 \(\rho g\cdot 2R\approx20\,\mathrm{Pa}\),远小于 \(p_\infty\approx10^{5}\,\mathrm{Pa}\);上浮的终端速度约为 0.2 m/s–0.3 m/s,而后面关心的塌缩、振荡发生在 1 μs–100 μs 时间内,泡在此期间只移动数十微米,远小于 \(R\)。故在球对称模型中忽略重力。
3.2 不可压缩与球对称:液体的速度场
液体密度 \(\rho\) 为常量。水的压缩率很小:压强改变 \(1\,\mathrm{atm}\),密度相对变化仅 \(\Delta p/(\rho c^2)\approx10^{5}/(10^{3}\times 1480^2)\approx4.5\times10^{-5}\)。这一假设要求液体速度远小于声速 \(c\)、过程时间远长于声波穿越气泡的时间 \(R/c\);在塌缩末期 \(\abs{\dot R}\) 可能接近 \(c\),此时该假设失效(第 5.5 节)。
取图 1 中半径 \(r_1<r_2\) 的两个同心球面(都在液体中)。两球面间液体的质量 \(\rho\cdot\frac43\pi(r_2^3-r_1^3)\) 不随时间变化(\(\rho\) 为常量、球面固定),所以单位时间从内球面流入的体积必须等于从外球面流出的体积: \[4\pi r_1^2\,u(r_1,t)=4\pi r_2^2\,u(r_2,t).\] 由于 \(r_1,r_2\) 任取,体积流量 \(q_V(t)\equiv 4\pi r^2u(r,t)\) 与 \(r\) 无关。再用泡壁条件定出 \(q_V\):
泡壁处液体紧贴泡壁运动,\(u(R,t)=\dot R\),即液体既不穿过泡壁、也不与泡壁分离。这要求忽略蒸发/凝结造成的质量流(对蒸气泡,修正量为 \(\rho_v/\rho\) 量级,常温下约 \(10^{-5}\))。
于是 \(q_V=4\pi R^2\dot R\), \[\tag{3.4} u(r,t)=\frac{R^2\dot R}{r^2}\qquad(r\ge R).\]
式 (3.4) 有一个后面反复用到的特点:\(t\) 时刻整个液体的速度场完全由 \(R\) 与 \(\dot R\) 两个数决定——气泡只有一个自由度。这也意味着远处的液体“瞬时”感受到泡壁运动,这是不可压缩近似(声速 \(\to\infty\))的结果。
3.3 液体的动能与有效质量
半径 \(r\)、厚 \(\dd r\) 的球壳(图 1 橙色)质量为 \(\rho\,4\pi r^2\dd r\),速度为 \(u(r)\)。整个液体的动能 \[\tag{3.5} K=\int_R^\infty \frac12\rho\,u^2\,4\pi r^2\,\dd r =2\pi\rho R^4\dot R^2\int_R^\infty\frac{\dd r}{r^2} =2\pi\rho R^3\dot R^2 .\] 积分收敛是因为 \(u^2r^2\propto r^{-2}\):动能主要集中在泡附近几个 \(R\) 以内(\(r>2R\) 的液体只贡献一半)。
把式 (3.5) 写成 \(K=\frac12 m_{\mathrm{eff}}\dot R^2\),则 \[\tag{3.6} m_{\mathrm{eff}}=4\pi\rho R^3=3\times\Bigl(\tfrac43\pi R^3\rho\Bigr),\] 即泡壁运动“拖动”的等效质量等于气泡所排开液体质量的 3 倍。这个“3 倍”在第 6 节的 Minnaert 共振中直接决定频率。
3.4 拉普拉斯压强
液体表面单位面积具有的表面自由能 \(\sigma\)(单位 J/m2\(=\)N/m),等价地,表面上单位长度分界线两侧相互的拉力。\(20\,{}^{\circ}\mathrm{C}\) 水–空气界面 \(\sigma=0.0728\,\mathrm{N/m}\)。
设泡内压强 \(p_{\mathrm{in}}\)、紧贴泡外的液体压强 \(p_{\mathrm{out}}\),两者保持不变,让半径增大 \(\dd R\)。压强差做功 \((p_{\mathrm{in}}-p_{\mathrm{out}})\,4\pi R^2\dd R\),全部用于增加表面能 \(\sigma\,\dd(4\pi R^2)=8\pi\sigma R\,\dd R\)(准静态、无动能)。于是 \[\tag{3.7} p_{\mathrm{in}}-p_{\mathrm{out}}=\frac{2\sigma}{R}.\]
把气泡沿赤道切开,上半球受压强差的合力 \((p_{\mathrm{in}}-p_{\mathrm{out}})\pi R^2\) 向上,赤道一圈表面张力 \(\sigma\cdot 2\pi R\) 向下,平衡即得式 (3.7)。液中气泡只有一个界面;肥皂泡有内外两个界面,压强差为 \(4\sigma/R\)。数值:\(R=1\,\mathrm{mm}\) 时 \(2\sigma/R=146\,\mathrm{Pa}\)(约 0.14 % 大气压);\(R=1\,\mu\mathrm{m}\) 时 \(2\sigma/R=1.46\times10^{5}\,\mathrm{Pa}\),已超过一个大气压。可见表面张力只对微米级气泡重要。
3.5 速度势与非定常伯努利方程
若存在标量函数 \(\phi(r,t)\) 使 \(u=\pdv{\phi}{r}\)(三维中 \(\vect{u}=\nabla\phi\)),称 \(\phi\) 为速度势,相应的流动称为势流。有速度势的流动必定无旋:流体微元在运动中只平移和变形、不自转(三维中 \(\nabla\times\vect u=\nabla\times\nabla\phi=0\));本讲的径向流显然如此。对式 (3.4),积分 \(\pdv{\phi}{r}=R^2\dot R/r^2\) 并取 \(r\to\infty\) 时 \(\phi\to0\),得 \[\tag{3.8} \phi(r,t)=-\frac{R^2\dot R}{r}.\]
把 \(u=\pdv{\phi}{r}\) 代入式 (3.3) 的左边两项,分别改写为对 \(r\) 的偏导:
交换求偏导次序:\(\pdv{u}{t}=\pdv{}{t}\pdv{\phi}{r}=\pdv{}{r}\pdv{\phi}{t}\);
乘积求导的逆用:\(u\pdv{u}{r}=\pdv{}{r}\Bigl(\frac{u^2}{2}\Bigr)\);
\(\rho\) 为常量,\(\frac1\rho\pdv{p}{r}=\pdv{}{r}\Bigl(\frac p\rho\Bigr)\)。
于是式 (3.3) 成为 \(\pdv{}{r}\Bigl[\pdv{\phi}{t}+\frac{u^2}{2}+\frac{p}{\rho}\Bigr]=0\):方括号内的量与 \(r\) 无关,只能是 \(t\) 的函数 \(C(t)\)。令 \(r\to\infty\):\(\phi\to0\) 且 \(\pdv{\phi}{t}\propto1/r\to0\),\(u\to0\),\(p\to p_\infty(t)\),故 \(C(t)=p_\infty(t)/\rho\), \[\tag{3.9} \pdv{\phi}{t}+\frac{u^2}{2}+\frac{p}{\rho}=\frac{p_\infty(t)}{\rho}.\]
式 (3.9) 成立需要:(1) 液体不可压缩;(2) 流动有速度势(本讲的径向流自动满足;一般三维流动要求无旋);(3) 体内粘性力可忽略(对势流,不可压缩牛顿流体的体内粘性力恒为零,附录 A;粘性只通过边界条件进入);(4) 重力已并入 \(p_\infty\)。若流动定常(\(\pdv{\phi}{t}=0\)),它退化为普通物理中的伯努利方程 \(p+\frac12\rho u^2=\) 常量。与定常形式相比,新出现的 \(\pdv{\phi}{t}\) 项代表“加速整段液柱所需的压强差”。
3.6 泡内气体:多方过程,等温还是绝热
理想气体满足 \(pV^\kappa=\) 常量的过程称为多方过程,\(\kappa\) 为多方指数。\(\kappa=1\) 为等温,\(\kappa=\gamma=c_p/c_V\) 为可逆绝热(空气 \(\gamma=1.40\),单原子气体如氩 \(\gamma=5/3\))。对半径为 \(R\) 的球形泡、物质的量不变, \[\tag{3.10} p_g=p_{g0}\Bigl(\frac{R_0}{R}\Bigr)^{3\kappa},\qquad T=T_0\Bigl(\frac{R_0}{R}\Bigr)^{3(\kappa-1)},\] 后一式由 \(pV=n R_{\mathrm{m}} T\) 与前一式相除得到(\(R_{\mathrm m}\) 为摩尔气体常量)。泡内总压强 \(p_B=p_g+p_v\),\(p_v\) 为液体的饱和蒸气压(20 °C 水 \(p_v\approx2.3\,\mathrm{kPa}\))。
\(\kappa\) 取 1 还是 \(\gamma\),取决于压缩过程中热量能否在泡内来得及扩散。
热传导方程 \(\pdv{T}{t}=D\pdv{^2T}{x^2}\) 中 \(D=k_{\mathrm{th}}/(\rho c_p)\) 称为热扩散率(单位 m2/s;\(k_{\mathrm{th}}\) 为热导率,\(c_p\) 为定压比热)。比较两边量纲可知:时间 \(t\) 内热量扩散的典型距离为 \[\tag{3.11} \ell_{\mathrm{th}}\sim\sqrt{Dt}.\] 20 °C、1 atm 空气 \(k_{\mathrm{th}}=0.0257\,\mathrm{W/(m\,K)}\),\(\rho=1.20\,\mathrm{kg/m^3}\),\(c_p=1005\,\mathrm{J/(kg\,K)}\),\(D=2.12\times10^{-5}\,\mathrm{m^2/s}\)。
液体的比热容与密度都远大于气体,泡壁处温度几乎保持在液体温度 \(T_0\)。若在一个振荡周期(时间尺度 \(1/\omega\))内热量能扩散过整个泡(\(\ell_{\mathrm{th}}=\sqrt{D/\omega}\gg R\)),泡内温度跟得上液体,过程近似等温;若 \(\ell_{\mathrm{th}}\ll R\),只有贴壁的薄层与液体换热,泡内主体近似绝热。判据是无量纲数(热 Péclet 数) \[\tag{3.12} \mathrm{Pe}\equiv\frac{\omega R^2}{D}=\Bigl(\frac{R}{\ell_{\mathrm{th}}}\Bigr)^2:\qquad \mathrm{Pe}\ll1\ \text{近似等温},\qquad \mathrm{Pe}\gg1\ \text{近似绝热}.\] 第 6.3 节会把这一定性判据变成定量的“有效多方指数” \(\kappa_{\mathrm{eff}}(R)\)。
3.7 粘性:泡壁上的法向粘性应力
普通物理中的牛顿粘性定律描述剪切:相邻液层速度差 \(\dd u_x/\dd y\) 引起切向应力 \(\tau_{xy}=\mu\,\dd u_x/\dd y\)。气泡胀缩时液体并不剪切,而是在径向被拉伸(或压缩)、在切向被压缩(或拉伸)。附录 A 证明:剪切与“一个方向拉伸、另一方向压缩”只是同一种形变在不同坐标轴下的样子,由此得到牛顿流体中沿某方向的法向粘性应力等于 \(2\mu\) 乘以该方向的拉伸速率。对式 (3.4),径向拉伸速率为 \[\tag{3.13} \pdv{u}{r}=-\frac{2R^2\dot R}{r^3}\quad\Longrightarrow\quad \tau_{rr}\big|_{r=R}=2\mu\pdv{u}{r}\Big|_{r=R}=-\frac{4\mu\dot R}{R}.\] 附录 A 还证明:在液体内部,这些粘性应力对每个液体微元的合力恰好为零(所以欧拉方程 (3.3) 在体内仍然精确),而整个液体的粘性耗散功率恰为 \(16\pi\mu R\dot R^2=(4\mu\dot R/R)\cdot 4\pi R^2\dot R\)。两种算法都表明:粘性的全部效果相当于在泡壁处额外施加一个阻碍运动的压强 \(4\mu\dot R/R\)。
3.8 线性声学:声速、球面波与单极子辐射
后面的阻尼(第 6.4 节)、含泡液体声速(第 7 节)和枪虾的声音(第 8.3 节)都需要“可压缩”的液体。这里在小扰动下补齐需要的三个结论。
物质的(绝热)压缩率 \(\beta\equiv-\frac1V\pdv{V}{p}=\frac1\rho\pdv{\rho}{p}\)。小扰动下压强扰动 \(p'\) 与密度扰动 \(\rho'\) 成正比:\(p'=c^2\rho'\),其中 \[\tag{3.14} c^2=\pdv{p}{\rho}=\frac{1}{\rho\beta}.\] 下面的推导将说明 \(c\) 正是声速。水 \(\rho c^2=2.2\times10^{9}\,\mathrm{Pa}\),\(c=1482\,\mathrm{m/s}\)。
设静止液体(密度 \(\rho_0\))中有小扰动 \(p',\rho',u\),只保留一阶小量:
质量守恒(对图 1 的球壳):壳内质量增加率 \(=\) 流入减流出,\(4\pi r^2\dd r\,\pdv{\rho'}{t}=-\pdv{}{r}\bigl(\rho_0u\,4\pi r^2\bigr)\dd r\),即 \(\pdv{\rho'}{t}=-\frac{\rho_0}{r^2}\pdv{}{r}(r^2u)\);
欧拉方程 (3.3) 去掉二阶的对流项:\(\rho_0\pdv{u}{t}=-\pdv{p'}{r}\);引入 \(u=\pdv{\phi}{r}\),同第 3.5 节的步骤得线性化伯努利关系 \(p'=-\rho_0\pdv{\phi}{t}\);
状态关系 \(\rho'=p'/c^2=-(\rho_0/c^2)\pdv{\phi}{t}\)。
把 3 代入 1,并用 \(u=\pdv{\phi}{r}\): \[\frac{1}{c^2}\pdv{^2\phi}{t^2}=\frac{1}{r^2}\pdv{}{r}\Bigl(r^2\pdv{\phi}{r}\Bigr)=\frac1r\pdv{^2(r\phi)}{r^2},\] 最后一步可直接展开两边验证(都等于 \(\phi''+2\phi'/r\))。于是 \(\psi\equiv r\phi\) 满足一维波动方程 \(\pdv{^2\psi}{t^2}=c^2\pdv{^2\psi}{r^2}\),其向外传播的解为 \(\psi=-F(t-r/c)\),\(F\) 为任意函数: \[\tag{3.15} \phi(r,t)=-\frac{F(t-r/c)}{r},\qquad p'(r,t)=-\rho_0\pdv{\phi}{t}=\frac{\rho_0F'(t-r/c)}{r}.\]
在泡附近 \(r\ll\lambda\)(\(\lambda\) 为声波长)处,推迟时间 \(r/c\) 可忽略,式 (3.15) 退化为 \(\phi\approx-F(t)/r\),与不可压缩解 (3.8) 比较得 \(F(t)=R^2\dot R\)。因此远处的声压为 \[\tag{3.16} p'(r,t)=\frac{\rho}{r}\,\dv{}{t}\bigl(R^2\dot R\bigr)\Big|_{t-r/c}=\frac{\rho\bigl(R^2\ddot R+2R\dot R^2\bigr)}{r}\Big|_{t-r/c}=\frac{\rho\,\ddot V(t-r/c)}{4\pi r},\] 其中 \(V=\frac43\pi R^3\) 为气泡体积(\(\dot V=4\pi R^2\dot R\),\(\ddot V=4\pi(R^2\ddot R+2R\dot R^2)\))。这种只由体积变化产生、各向同性的声源叫单极子。在远场(\(r\gg\lambda\)),球面波局部近似为平面波,\(u=p'/(\rho c)\),声强 \(I=\avg{p'u}=\avg{p'^2}/(\rho c)\)。
式 (3.16) 的近场形式 \(\rho(R^2\ddot R+2R\dot R^2)/r\) 正是不可压缩伯努利方程 (3.9) 中 \(-\rho\pdv{\phi}{t}\) 一项;不可压缩理论中另一项 \(-\rho u^2/2\propto r^{-4}\) 在远处衰减得快得多,是二阶小量。两者一致。
4 Rayleigh–Plesset 方程
4.1 泡壁处的力平衡
泡壁(气液界面)本身没有质量,作用在它上面的力必须时时平衡。取一小块面积 \(\dd A\) 的界面,以向外(指向液体)为正:
泡内气体与蒸气向外推:\(+p_B\,\dd A\);
表面张力:弯曲界面的张力合力指向曲率中心,效果相当于向内的压强 \(2\sigma/R\)(式 (3.7)),即 \(-\frac{2\sigma}{R}\dd A\);
液体对界面的法向力:液体压强向内推 \(-p_L\,\dd A\);法向粘性应力 \(\tau_{rr}\) 按“拉为正”的约定向外拉 \(+\tau_{rr}\dd A\)。
令合力为零,并代入式 (3.13): \[\tag{4.1} p_L(R,t)=p_B-\frac{2\sigma}{R}-\frac{4\mu\dot R}{R},\] 其中 \(p_L(R,t)\) 是紧贴泡壁的液体压强。静止时(\(\dot R=0\))它退化为拉普拉斯公式。
4.2 推导一:能量方法
研究对象:\(R<r<L\) 的液体,\(L\to\infty\)。先不计粘性(粘性项最后由式 (4.1) 加回,附录 A 说明这样做与“功\(-\)耗散”的算法等价)。
内边界 \(r=R\):泡壁以压强 \(p_L(R)\) 推液体,力 \(p_L\cdot4\pi R^2\) 沿 \(+r\),液体速度 \(\dot R\),功率 \(p_L\,4\pi R^2\dot R\)。
外边界 \(r=L\):外部液体以 \(p_\infty\) 向内推,功率 \(-p_\infty\,4\pi L^2u(L)=-p_\infty\,4\pi R^2\dot R\)(用了体积流量守恒)。
不可压缩液体内部压强不做净功(没有体积变化),故动能定理为 \[\dv{K}{t}=\bigl(p_L(R)-p_\infty\bigr)\,4\pi R^2\dot R .\]
对式 (3.5) 求导:\(\dv{K}{t}=2\pi\rho\bigl(3R^2\dot R^3+2R^3\dot R\ddot R\bigr)=4\pi\rho R^2\dot R\Bigl(R\ddot R+\frac32\dot R^2\Bigr)\)。
两式相等,约去 \(4\pi R^2\dot R\)(对 \(\dot R\neq0\) 的任意时刻成立,由连续性推广到全部时刻):\(R\ddot R+\frac32\dot R^2=\bigl(p_L(R)-p_\infty\bigr)/\rho\)。
4.3 推导二:非定常伯努利方程
由式 (3.8),\(\pdv{\phi}{t}=-\dfrac{R^2\ddot R+2R\dot R^2}{r}\)(\(R\)、\(\dot R\) 都是 \(t\) 的函数,乘积求导)。在 \(r=R\) 处:\(\pdv{\phi}{t}=-(R\ddot R+2\dot R^2)\),\(u^2/2=\dot R^2/2\),\(p=p_L(R)\)。代入式 (3.9): \[-R\ddot R-2\dot R^2+\frac12\dot R^2+\frac{p_L(R)}{\rho}=\frac{p_\infty}{\rho} \quad\Longrightarrow\quad R\ddot R+\frac32\dot R^2=\frac{p_L(R)-p_\infty}{\rho},\] 与能量方法结果相同。
两种推导再与式 (4.1) 合并,就得到本讲的核心方程。
\[\tag{4.2} R\ddot R+\frac32\dot R^2=\frac{1}{\rho}\left[p_B(t)-p_\infty(t)-\frac{2\sigma}{R}-\frac{4\mu\dot R}{R}\right], \qquad p_B=p_v+p_{g0}\Bigl(\frac{R_0}{R}\Bigr)^{3\kappa}.\] Rayleigh 于 1917 年给出了无表面张力、无粘性的形式 [5];Plesset 1949 年用随时间变化的远场压强研究空化泡 [21];含气体与表面张力项的完整形式由 Noltingk 与 Neppiras(1950、1951)首先写出并用于超声空化 [22], [23],粘性项由 Poritsky 1952 年首先研究 [4], [24]。
逐项物理意义
\(R\ddot R\):泡壁加速时要推动的液体惯性(与 \(m_{\mathrm{eff}}\ddot R\) 对应);
\(\frac32\dot R^2\):即使泡壁匀速,参与运动的液体也在变化——泡变小时同样的 \(\dot R\) 对应更少的液体动能,多出来的动能必须“挤”到更快的速度上。正是这一非线性项造成塌缩末期的速度奇点;
\(p_B-p_\infty\):驱动力。\(p_\infty>p_B\) 时泡被压缩;\(p_\infty(t)\) 可以是声场或流场造成的压强变化;
\(-2\sigma/R\):表面张力,总是倾向于缩小气泡;
\(-4\mu\dot R/R\):粘性,总与 \(\dot R\) 反号,阻碍运动。
左边量纲 \(\mathrm{m^2/s^2}\);右边 \(\mathrm{Pa}/(\mathrm{kg/m^3})=\mathrm{m^2/s^2}\),\(\sigma/(\rho R)\) 与 \(\mu\dot R/(\rho R)\) 亦同。静态极限 \(\dot R=\ddot R=0\):\(p_B=p_\infty+2\sigma/R\),回到拉普拉斯公式。若 \(\sigma=\mu=0\)、\(p_B=p_\infty\),则 \(R\ddot R+\frac32\dot R^2=0\) 即 \(\dv{}{t}(R^3\dot R^2)=0\),液体动能守恒——无驱动时没有能量输入,与式 (3.5) 一致。
球对称:远离壁面与其他气泡,忽略重力引起的上浮与形变;形状扰动稳定——即泡面上偶然出现的小起伏不会自发长大(剧烈塌缩末期这一点可能失效,泡会破碎或形成射流;本讲不讨论)。
液体不可压缩:\(\abs{\dot R}\ll c\)。第 5.5 节给出失效的位置。
泡内压强均匀:泡内气体的声速远大于 \(\abs{\dot R}\),且过程时间远长于 \(R/c_{\mathrm{gas}}\)。
泡内气体物质的量不变(忽略溶解气体的扩散,其时间尺度为秒量级),蒸气压 \(p_v\) 取液体温度下的饱和值(忽略蒸发、凝结的动力学)。
5 Rayleigh 塌缩
1917 年,Rayleigh 在与 Parsons 关心的螺旋桨空化问题相同的背景下(见第 1 节)研究了如下问题:液体中一个半径 \(R_0\) 的空腔(泡内无物,\(p_B=0\))在远场压强 \(p_\infty\) 作用下由静止开始塌缩,需要多长时间?液体中会产生多大的压强?[5]
\(p_B=0\)(空腔,或 \(p_v\ll p_\infty\) 的蒸气泡并把 \(p_\infty-p_v\) 记作 \(p_\infty\)),\(\sigma=\mu=0\),\(p_\infty\) 为常量,\(t=0\) 时 \(R=R_0\)、\(\dot R=0\)。
5.1 塌缩速度:能量积分
从 \(R_0\) 塌缩到 \(R\),远场压强 \(p_\infty\) 推着液体使空腔体积减少 \(\frac43\pi(R_0^3-R^3)\),做功 \(\frac43\pi p_\infty(R_0^3-R^3)\);空腔内没有物质,不做功;液体不可压缩,内部无储能。所以这部分功全部变为液体动能 (3.5): \[2\pi\rho R^3\dot R^2=\frac43\pi p_\infty\bigl(R_0^3-R^3\bigr),\] \[\tag{5.1} \dot R^2=\frac{2p_\infty}{3\rho}\left(\frac{R_0^3}{R^3}-1\right),\qquad \dot R<0 .\]
式 (5.1) 两边对 \(t\) 求导:\(2\dot R\ddot R=\frac{2p_\infty}{3\rho}\cdot\bigl(-3R_0^3R^{-4}\bigr)\dot R\),得 \(\ddot R=-\frac{p_\infty}{\rho}\frac{R_0^3}{R^4}\)。于是 \(R\ddot R+\frac32\dot R^2=-\frac{p_\infty}{\rho}\frac{R_0^3}{R^3}+\frac{p_\infty}{\rho}\Bigl(\frac{R_0^3}{R^3}-1\Bigr)=-\frac{p_\infty}{\rho}\),正是 \(p_B=0\) 时 RP 方程的右边。能量方法与运动方程一致(这本是同一定律的两种写法)。
自然的时间单位与速度单位是 \[\tag{5.2} \tau\equiv R_0\sqrt{\frac{\rho}{p_\infty}},\qquad U_0\equiv\sqrt{\frac{p_\infty}{\rho}} .\] 水中 \(p_\infty=1\,\mathrm{atm}\) 时 \(U_0=10.1\,\mathrm{m/s}\),\(R_0=1\,\mathrm{mm}\) 时 \(\tau=99\,\mu\mathrm{s}\)。
5.2 塌缩时间
由式 (5.1),\(\dd t=\dd R/\dot R=-\dd R/\abs{\dot R}\),从 \(R_0\) 积到 \(0\): \[t_c=\int_0^{R_0}\frac{\dd R}{\sqrt{\frac{2p_\infty}{3\rho}\bigl(R_0^3/R^3-1\bigr)}} =R_0\sqrt{\frac{3\rho}{2p_\infty}}\int_0^1\frac{x^{3/2}\,\dd x}{\sqrt{1-x^3}},\qquad x\equiv\frac{R}{R_0}.\] (第二步:被积函数分子分母同乘 \(x^{3/2}\),\(\sqrt{x^{-3}-1}=x^{-3/2}\sqrt{1-x^3}\)。)令 \(s=x^3\),积分化为欧拉 Beta 函数(附录 B): \[\int_0^1\frac{x^{3/2}\dd x}{\sqrt{1-x^3}}=\frac13\int_0^1 s^{-1/6}(1-s)^{-1/2}\dd s=\frac13\,\mathrm{B}\Bigl(\frac56,\frac12\Bigr)=0.74683 .\]
\[\tag{5.3} \boxed{\,t_c=\sqrt{\frac{\pi}{6}}\,\frac{\Gamma(5/6)}{\Gamma(4/3)}\,R_0\sqrt{\frac{\rho}{p_\infty}}=0.91468\,R_0\sqrt{\frac{\rho}{p_\infty}}\,}\]
把 RP 方程无量纲化(\(R\) 以 \(R_0\)、\(t\) 以 \(\tau\) 为单位):\(\ddot R=-(1+\frac32\dot R^2)/R\)。用 8 阶 Runge–Kutta 积分到 \(R=10^{-4}R_0\),剩余时间用下面的渐近律补上,得 \(t_c/\tau=0.914681\),与式 (5.3) 相差小于 \(10^{-6}\)(verify_core.py)。数值解的 \(\dot R\) 沿整条轨道与能量积分 (5.1) 的相对偏差小于 \(10^{-7}\)(图 2 右)。量纲检验:\(R_0\sqrt{\rho/p}\) 的单位是 \(\mathrm{m}\cdot\sqrt{\mathrm{kg\,m^{-3}}/\mathrm{Pa}}=\mathrm{s}\)。物理上,\(t_c\) 正比于 \(R_0\):泡越大,要加速的液体越多(\(\propto R_0^3\)),而驱动力做功也 \(\propto R_0^3\),速度标度与 \(R_0\) 无关,时间只随距离 \(R_0\) 线性增长。
5.3 速度奇点与末端渐近律
当 \(R\ll R_0\) 时式 (5.1) 中的 \(-1\) 可忽略(相对误差 \(R^3/R_0^3\)): \[\tag{5.4} \abs{\dot R}\approx\sqrt{\frac{2p_\infty}{3\rho}}\Bigl(\frac{R_0}{R}\Bigr)^{3/2}\to\infty\qquad(R\to0).\] 分离变量 \(R^{3/2}\dd R=-\sqrt{2p_\infty/3\rho}\,R_0^{3/2}\dd t\),从 \(t\) 积到 \(t_c\)(此时 \(R=0\)): \[\tag{5.5} R(t)\approx\Bigl(\frac52\Bigr)^{2/5}\Bigl(\frac{2p_\infty}{3\rho}\Bigr)^{1/5}R_0^{3/5}\,(t_c-t)^{2/5}.\] 数值解在 \(2\times10^{-4}<R/R_0<3\times10^{-2}\) 区间拟合出的指数为 \(0.400\)。物理图像:能量守恒要求 \(\frac12m_{\mathrm{eff}}\dot R^2\approx\frac43\pi p_\infty R_0^3\) 为常量,而 \(m_{\mathrm{eff}}\propto R^3\) 越来越小,速度只能越来越大——有限的能量被集中到越来越少的液体上。这就是空化“破坏力”的根源。
5.4 液体中的压强峰
塌缩中液体被加速后又要在泡壁附近“刹车”,因此液体内部的压强可以远高于 \(p_\infty\)。
由式 (3.9),\(p(r)=p_\infty-\rho\pdv{\phi}{t}-\frac12\rho u^2\),其中 \(\pdv{\phi}{t}=-(R^2\ddot R+2R\dot R^2)/r\),\(u=R^2\dot R/r^2\): \[p(r)=p_\infty+\frac{\rho\bigl(R^2\ddot R+2R\dot R^2\bigr)}{r}-\frac{\rho R^4\dot R^2}{2r^4}.\] 代入 \(\ddot R=-\frac{p_\infty}{\rho}\frac{R_0^3}{R^4}\) 与式 (5.1),记 \(z\equiv R_0^3/R^3\)(\(z\) 从 1 增大到 \(\infty\)):\(\rho(R^2\ddot R+2R\dot R^2)=p_\infty R\bigl[-z+\frac43(z-1)\bigr]=\frac{p_\infty R}{3}(z-4)\),\(\frac12\rho R^4\dot R^2=\frac{p_\infty R^4}{3}(z-1)\)。于是 \[\tag{5.6} \frac{p(r)}{p_\infty}=1+\frac{R}{3r}(z-4)-\frac{R^4}{3r^4}(z-1).\]
\(r=R\):\(1+\frac{z-4}{3}-\frac{z-1}{3}=0\),与空腔壁上 \(p_L=0\) 一致;\(r\to\infty\):\(p\to p_\infty\)。
当 \(z\le4\)(即 \(R\ge4^{-1/3}R_0=0.630R_0\))时,式 (5.6) 随 \(r\) 单调上升到 \(p_\infty\),液体中没有超过 \(p_\infty\) 的压强。当 \(z>4\) 时出现极大值:令 \(\dv{p}{r}=0\), \[-\frac{R}{3r^2}(z-4)+\frac{4R^4}{3r^5}(z-1)=0\quad\Longrightarrow\quad r_m^3=\frac{4(z-1)}{z-4}R^3 .\] 代回式 (5.6):记 \(y=R/r_m\),则 \(y^3(z-1)=(z-4)/4\),\(\frac{p_{\max}}{p_\infty}=1+\frac y3\bigl[(z-4)-y^3(z-1)\bigr]=1+\frac{y(z-4)}{4}\),即 \[\tag{5.7} \frac{p_{\max}}{p_\infty}=1+\frac{(z-4)^{4/3}}{4^{4/3}(z-1)^{1/3}}\ \xrightarrow{\ z\gg1\ }\ \frac{z}{4^{4/3}}=\frac{1}{6.35}\frac{R_0^3}{R^3},\qquad r_m\to4^{1/3}R=1.587R .\] 例如 \(R=R_0/20\) 时 \(p_{\max}\approx1260\,p_\infty\),位于距泡壁约 \(0.59R\) 处的液体中(图 3)。这就是 Rayleigh 1917 年的结论 [5]:塌缩空腔周围的液体中会出现正比于 \((R_0/R)^3\) 的高压,他据此认为这可以解释螺旋桨的金属表面被“打出”凹坑。注意这个峰只在 \(R\ll R_0\) 的末段出现、且位于液体中 \(r\approx1.59R\) 处。后来的实验与计算表明,对紧贴壁面的空泡,更重要的是不对称塌缩形成的指向壁面的高速微射流,以及反弹时向外传播的冲击波(振幅很大、波前陡峭的压缩波) [4], [25], [26]。
5.5 何时必须考虑液体可压缩性
式 (5.4) 中 \(\abs{\dot R}\) 无限增大,必然在某处达到声速 \(c\),不可压缩假设失效。令 \(\frac{2p_\infty}{3\rho}\frac{R_0^3}{R^3}=c^2\): \[\tag{5.8} \frac{R_{\mathrm{M}}}{R_0}=\Bigl(\frac{2p_\infty}{3\rho c^2}\Bigr)^{1/3}\approx0.031\qquad(p_\infty=10^{5}\,\mathrm{Pa},\ c=1480\,\mathrm{m/s}).\] 所以对空腔,RP 方程在 \(R\gtrsim R_0/30\) 时可信。更接近塌缩点时需要计入液体可压缩性的一阶修正。Keller 与 Miksis 等把单极子辐射场(第 3.8 节)与 RP 方程衔接,得到“弱可压缩 RP 方程” [4], [27],其中最重要的新项代表泡向外辐射声能,第 6.4 节将在小振幅情形下推导它。实际上,只要泡内有一点气体,气体的“缓冲”就会先让塌缩停下来——这是下一小节的内容。
5.6 含气空泡:反弹、最小半径与最高温度
\(t=0\) 时泡处于最大半径 \(R_0\) 且静止;泡内有不凝结气体,压强 \(p_{g0}\ll p_\infty\);压缩为绝热(\(\kappa=\gamma\)),忽略 \(p_v\)、\(\sigma\)、\(\mu\) 与可压缩性。绝热的判据仿照式 (3.12):在运动的时间尺度 \(t_*\) 内热扩散距离 \(\sqrt{Dt_*}\) 应远小于 \(R\),即 \(\mathrm{Pe}\sim R^2/(Dt_*)\gg1\)。对毫米泡,\(t_*\sim t_c\approx{90\,}\mathrm{\mu s}\),\(\mathrm{Pe}\sim500\),全程近似绝热;对声致发光中 \(R\approx{50\,}\mathrm{\mu m}\) 的泡,取塌缩的 Rayleigh 时间 \(t_*\approx3\)~\(5\,\mu\mathrm{s}\)(\(p_\infty\) 为 1~2.3 bar),按常压空气的 \(D\) 有 \(\mathrm{Pe}\sim30\);但此时泡内气体比常压稀薄约 \(780\) 倍,气体热扩散率 \(D=k_{\mathrm{th}}/(\rho_gc_p)\propto1/\rho_g\) 相应增大,实际 \(\mathrm{Pe}\approx0.03\ll1\),前段接近等温,只有最后极短的剧烈阶段才近似绝热(见本节末的注)。
从 \(R_0\) 到 \(R\),外压做功 \(\frac43\pi p_\infty(R_0^3-R^3)\) 一部分变为液体动能,另一部分用于压缩气体。绝热压缩时外界对气体做的功等于内能增量: \[W_g=-\int_{V_0}^{V}p_g\,\dd V'=-p_{g0}V_0^\gamma\int_{V_0}^{V}V'^{-\gamma}\dd V'=\frac{p_{g0}V_0}{\gamma-1}\Bigl[\Bigl(\frac{V_0}{V}\Bigr)^{\gamma-1}-1\Bigr].\] 能量守恒:\(2\pi\rho R^3\dot R^2=\frac43\pi p_\infty(R_0^3-R^3)-W_g\)。在最小半径 \(R_{\min}\) 处 \(\dot R=0\),记 \(x=R_{\min}/R_0\),两边除以 \(\frac43\pi R_0^3\): \[\tag{5.9} p_\infty\bigl(1-x^3\bigr)=\frac{p_{g0}}{\gamma-1}\Bigl(x^{-3(\gamma-1)}-1\Bigr).\]
当 \(p_{g0}\ll(\gamma-1)p_\infty\) 时 \(x\ll1\),式 (5.9) 左边 \(\approx p_\infty\)、右边 \(\approx\frac{p_{g0}}{\gamma-1}x^{-3(\gamma-1)}\): \[\tag{5.10} \frac{R_{\min}}{R_0}\approx\Bigl[\frac{p_{g0}}{(\gamma-1)p_\infty}\Bigr]^{\frac{1}{3(\gamma-1)}},\qquad \frac{T_{\max}}{T_0}=x^{-3(\gamma-1)}=1+\frac{(\gamma-1)p_\infty(1-x^3)}{p_{g0}}\approx\frac{(\gamma-1)p_\infty}{p_{g0}} .\] 被略去的是相对大小为 \(p_{g0}/[(\gamma-1)p_\infty]\) 的项(\(x^{-3(\gamma-1)}\) 旁的 1)和 \(x^3\);\(R_{\min}\) 近似式的相对误差约为 \(p_{g0}/[3(\gamma-1)^2p_\infty]\)。物理含义:外压对体积 \(\frac43\pi R_0^3\) 所做的功 \(p_\infty V_0\) 几乎全部集中到极少量气体的内能 \(p_{g0}V_0/(\gamma-1)\) 上,温度升高倍数就是两者之比。
对 \(\gamma=1.4\),\(p_{g0}/p_\infty=10^{-1},10^{-2},10^{-3}\),数值积分得 \(R_{\min}/R_0=0.2648,\ 0.04529,\ 0.006772\),与精确能量式 (5.9) 的根相对差 \(<10^{-5}\);近似式 (5.10) 给出 \(0.3150,\ 0.04623,\ 0.006786\),在 \(p_{g0}/p_\infty\le10^{-2}\) 时误差 \(\le2\%\)。\(T_{\max}/T_0\) 的数值结果为 \(4.9,\ 41.0,\ 401\);式 (5.10) 最右的近似给出 \(4,\ 40,\ 400\),保留“1”而只略去 \(x^3\) 的 \(1+(\gamma-1)p_\infty/p_{g0}\) 给出 \(5,\ 41,\ 401\)(图 4 右所画即此式)。\(p_{g0}/p_\infty=10^{-3}\) 时到达最小半径的时间为 \(0.916\tau\),与空腔塌缩时间 \(0.915\tau\) 几乎相同——气体只在最后极短时间内起作用。之后泡无损耗地反弹回 \(R_0\),周而复始(图 4)。实际中声辐射、热传导与粘性使反弹逐次减弱。
交互演示加载中……(需启用 JavaScript)
单泡声致发光(SBSL)的典型参数 [12]:平衡半径 \(R_a\approx{5\,}\mathrm{\mu m}\) 的(氩)气泡在频率 20 kHz–40 kHz、振幅 \(p_a\approx1.2\)~\(1.4\,\mathrm{bar}\) 的声场中,先缓慢膨胀到 \(R_{\max}\approx{50\,}\mathrm{\mu m}\),再剧烈塌缩到约 0.5 μm,并发出宽约 100 ps–200 ps 的闪光。按上面的模型:平衡时泡内气压为 \(p_0+2\sigma/R_a=1.013\times10^{5}\,\mathrm{Pa}+2.9\times10^{4}\,\mathrm{Pa}\)(表面张力项对 5 μm 的泡不可忽略),缓慢(等温)膨胀到最大时 \[p_{g0}\approx\Bigl(p_0+\frac{2\sigma}{R_a}\Bigr)\Bigl(\frac{R_a}{R_{\max}}\Bigr)^3\approx1.3\times10^{-3}p_0 .\] 取塌缩时 \(p_\infty\approx p_0\)、\(\gamma=5/3\),式 (5.10) 给出 \(R_{\min}\approx{2.2\,}\mathrm{\mu m}\)、\(T_{\max}\approx293\,\mathrm{K}\times\frac{2/3}{1.3\times10^{-3}}\approx1.5\times10^{5}\,\mathrm{K}\);若取更合理的 \(p_\infty\approx p_0+p_a\approx{2.3\,}\mathrm{bar}\)(塌缩发生在声压的正半周),则 \(R_{\min}\approx{1.5\,}\mathrm{\mu m}\)、\(T_{\max}\approx3.5\times10^{5}\,\mathrm{K}\)(verify_core.py)。而考虑传热、水蒸气等因素的精细模型给出的峰值温度在 6000 K 到略高于 2 × 104 K 之间 [12],在浓硫酸中对氩发射谱的测量给出 4000 K–15000 K [28]。
这里有一个必须正视的矛盾:模型预言 \(R_{\min}\approx1.5\)~\(2\,\mu\mathrm{m}\),实测约 0.5 μm——真实气泡被压得更小,却更冷。在绝热模型中这不可能(压得越小越热),需要分两层解释。
“更小”来自前段等温。 按上面的 Pe 判据,塌缩前段 \(\mathrm{Pe}\ll1\),气体近似等温压缩,泡内气压只按 \(R^{-3}\) 而不是 \(R^{-5}\) 增长,“气垫”软得多。用两段模型估算(verify_sbsl2.py):沿等温塌缩的轨道计算局部 \(\mathrm{Pe}=R\abs{\dot R}/D(p_g)\)(\(D\propto1/p_g\)),它在 \(R_1\approx17\)~\(19\,\mu\mathrm{m}\) 处升到 1,此后改为绝热。结果 \(p_\infty=p_0\) 时 \(R_{\min}\approx{0.76\,}\mathrm{\mu m}\),\(p_\infty={2.3\,}\mathrm{bar}\) 时 \(R_{\min}\approx{0.57\,}\mathrm{\mu m}\),与实测相当。
“更冷”却不能由此解释。 无损耗时,最高温度由外压做的总功(\(\approx p_\infty V_0\))与气体的量(\(\propto p_{g0}V_0\))之比决定。等温段流进液体的热量 \(p_{g0}V_0\ln(V_0/V_1)\) 只占总功的 \(0.2\%\)~\(0.4\%\),所以两段模型的 \(T_{\max}\) 仍为 \(1.5\times10^{5}\,\mathrm{K}\)(1 bar)与 \(3.4\times10^{5}\,\mathrm{K}\)(2.3 bar),与纯绝热估算几乎相同。降温必须来自最后几纳秒的能量出口,按重要性大致是:(1) 末期 \(\abs{\dot R}\) 接近液体声速,大部分能量以声波(冲击波)辐射掉;(2) 泡内水蒸气的离解与凝结吸收能量;(3) 末段仍不可忽略的热传导。至于气体分子本身的体积(“硬核”,按氩的 van der Waals 常量估算半径约 0.6 μm,与实测 \(R_{\min}\approx{0.5\,}\mathrm{\mu m}\) 相当,两者都只是量级值),它决定塌缩停在多大半径,并不能在能量给定时降低温度。所以本节的绝热估算只抓住了机制——把 bar 量级外压在大体积上做的功集中到极少量气体中——而定量的温度需要计入液体可压缩性、水蒸气与传热。
6 Minnaert 共振:气泡是一个“液体质量 + 气体弹簧”振子
1933 年 Minnaert 发表《论会唱歌的气泡与流水的声音》,指出溪水的潺潺声主要来自刚被卷入的气泡的体积振荡,并给出了振荡频率公式 [7]。本节从 RP 方程把它推出来。
6.1 线性化 RP 方程
远场压强恒为 \(p_0\),平衡半径 \(R_0\) 处泡内气压由拉普拉斯条件定出 \(p_{g0}=p_0+2\sigma/R_0\)(暂忽略 \(p_v\),或把它并入 \(p_0\))。令 \(R=R_0(1+\xi)\),\(\abs{\xi}\ll1\),只保留 \(\xi\) 的一阶项。气体按多方指数 \(\kappa\) 变化(\(\kappa\) 的取值见第 6.3 节)。
左边:\(R\ddot R=R_0^2(1+\xi)\ddot\xi=R_0^2\ddot\xi+\order{\xi^2}\);\(\frac32\dot R^2=\frac32R_0^2\dot\xi^2=\order{\xi^2}\),略去。
气压:由二项式展开 \((1+\xi)^{-3\kappa}=1-3\kappa\xi+\order{\xi^2}\),\(p_g=p_{g0}(1-3\kappa\xi)\)。
表面张力:\(\frac{2\sigma}{R}=\frac{2\sigma}{R_0}(1-\xi)\)。
粘性:\(\frac{4\mu\dot R}{R}=4\mu\dot\xi+\order{\xi^2}\)。
零阶项:\(p_{g0}-p_0-2\sigma/R_0=0\)(平衡条件),抵消。
收集一阶项,除以 \(\rho R_0^2\): \[\tag{6.1} \ddot\xi+\frac{4\mu}{\rho R_0^2}\dot\xi+\omega_0^2\,\xi=0,\qquad \omega_0^2=\frac{3\kappa p_{g0}-2\sigma/R_0}{\rho R_0^2}.\]
把 \(p_{g0}=p_0+2\sigma/R_0\) 代入:
\[\tag{6.2} \omega_0^2=\frac{1}{\rho R_0^2}\left[3\kappa p_0+(3\kappa-1)\frac{2\sigma}{R_0}\right],\qquad f_0=\frac{\omega_0}{2\pi}\approx\frac{1}{2\pi R_0}\sqrt{\frac{3\kappa p_0}{\rho}}\quad\Bigl(\frac{2\sigma}{R_0}\ll p_0\Bigr).\]
空气泡、绝热 \(\kappa=1.4\):取 \(p_0=1.000\times10^{5}\,\mathrm{Pa}\)、\(\rho=1000\,\mathrm{kg/m^3}\) 得 \(f_0R_0=3.26\,\mathrm{m/s}\)(文献常引的数值 [29]);取标准大气压 \(p_0=1.013\times10^{5}\,\mathrm{Pa}\)、20 °C 水 \(\rho=998\,\mathrm{kg/m^3}\) 得 \(3.29\,\mathrm{m/s}\)。若等温(\(\kappa=1\))则为 2.78 m/s。所以半径 1 mm 的气泡“唱” 3.3 kHz,1 cm 的泡唱 330 Hz,人耳可闻范围(20 Hz 至 20 kHz)对应 0.16 mm 到 16 cm 的气泡,常见流水声主要来自 0.2 mm 到 1 cm 的气泡。表面张力修正的相对大小为 \(\frac{3\kappa-1}{3\kappa}\frac{2\sigma}{R_0p_0}\),\(R_0=1\,\mathrm{mm}\) 时仅 0.1 %。
6.2 质量–弹簧图像
以泡壁位移 \(s\equiv R-R_0=R_0\xi\) 为坐标:
质量:液体动能 \(K=\frac12(4\pi\rho R_0^3)\dot s^2\),有效质量 \(m_{\mathrm{eff}}=4\pi\rho R_0^3\),即排开液体质量的 3 倍(式 (3.6))。振荡的惯性完全来自液体,泡内气体的质量(与 \(m_{\mathrm{eff}}\) 之比为 \(\rho_g/(3\rho)\approx4\times10^{-4}\))可忽略。
弹簧:泡壁受到的净向外力 \(F=\bigl(p_g-p_0-\frac{2\sigma}{R}\bigr)4\pi R_0^2\approx-\bigl(3\kappa p_{g0}-\frac{2\sigma}{R_0}\bigr)\frac{s}{R_0}\cdot4\pi R_0^2\),劲度系数 \(k_s=4\pi R_0\bigl(3\kappa p_{g0}-2\sigma/R_0\bigr)\)。弹性完全来自气体的可压缩性(表面张力略微“软化”这个弹簧:\(R\) 变大时 \(2\sigma/R\) 变小,帮助泡继续变大)。
于是 \(\omega_0^2=k_s/m_{\mathrm{eff}}\),与式 (6.1) 相同。一个 1 mm 的气泡,\(m_{\mathrm{eff}}=1.25\times10^{-5}\,\mathrm{kg}\),\(k_s\approx5.3\times10^{3}\,\mathrm{N/m}\)——一根“很硬的弹簧”挂着“一滴水”,所以频率高达 kHz。
刚性球在液体中整体平动时,被带动液体的附加质量是排开液体质量的 \(1/2\);而这里泡壁径向胀缩,速度场 \(\propto1/r^2\) 衰减得比平动的偶极场(\(\propto1/r^3\))慢,被带动的液体更多,等效质量是排开液体质量的 3 倍。二者对应不同的运动方式,不能混用。
交互演示加载中……(需启用 JavaScript)
6.3 等温还是绝热:有效多方指数
式 (6.2) 中的 \(\kappa\) 不是任意常数,而由泡内热传导决定。按式 (3.12),对 Minnaert 振荡 \[\mathrm{Pe}=\frac{\omega_0R_0^2}{D}\approx\frac{2\pi\times(3.3\,\mathrm{m/s})\times R_0}{D}=\frac{R_0}{1.0\,\mu\mathrm{m}} ,\] 即 \(R_0\gg1\,\mu\mathrm{m}\) 的泡接近绝热,\(R_0\ll1\,\mu\mathrm{m}\) 的泡接近等温。毫米泡的热扩散长度 \(\ell_{\mathrm{th}}=\sqrt{D/\omega_0}\approx32\,\mu\mathrm{m}\),只占半径的 3%。
附录 C 求解了泡内的线性化热传导方程(泡壁温度固定为液体温度),结果可以写成一个复多方指数: \[\tag{6.3} \kappa_c=\frac{\gamma}{1+(\gamma-1)G(y)},\qquad G(y)=\frac{3\,(y\coth y-1)}{y^2},\qquad y=R_0\sqrt{\frac{\ii\omega}{D}} .\] 这里沿用交流电路中的复数写法:振荡量写成“复振幅 \(\times\ee^{\ii\omega t}\)”,物理量取实部;\(\kappa_c\) 定义为压强相对扰动 \(p'/p_{g0}\) 与压缩量 \(-\delta V/V=-3\xi\) 的复振幅之比。其实部 \(\kappa_{\mathrm{eff}}=\Re\kappa_c\) 代替式 (6.2) 中的 \(\kappa\)。虚部 \(\Im\kappa_c>0\) 的含义:若压缩量 \(\propto\cos\omega t\),则压强 \(\propto\abs{\kappa_c}\cos(\omega t+\theta)\),\(\theta=\arg\kappa_c>0\),即泡内压强的变化超前于压缩(等价地,体积响应滞后于压强)。于是压缩阶段的压强高于膨胀阶段,\(p\)–\(V\) 图上的循环围出面积,每周期外界对气体做净功 \(-\oint p\,\dd V>0\),这部分振动能经热传导变成液体中的热,即热阻尼。两个极限:\(\abs{y}\to0\) 时 \(G\to1\),\(\kappa_c\to1\)(等温);\(\abs{y}\to\infty\) 时 \(G\approx3/y\to0\),\(\kappa_c\to\gamma\)(绝热)。图 5 右给出了 \(\kappa_{\mathrm{eff}}(R_0)\):1 mm 空气泡 \(\kappa_{\mathrm{eff}}=1.36\),0.1 mm 为 \(1.28\),10 μm 为 \(1.08\)。
6.4 阻尼:粘性、热传导与声辐射
把振子方程写成 \(\ddot\xi+2\Gamma\dot\xi+\omega_0^2\xi=0\),振幅按 \(\ee^{-\Gamma t}\) 衰减、能量按 \(\ee^{-2\Gamma t}\) 衰减。定义无量纲阻尼 \(\delta\equiv2\Gamma/\omega_0\),品质因数 \(Q\equiv1/\delta=\omega_0E/P_{\mathrm{loss}}\)(\(E\) 为振动能,\(P_{\mathrm{loss}}\) 为平均损耗功率)。几种互相独立的小损耗机制的 \(\delta\) 直接相加。振幅衰减到 \(1/\ee\) 的时间为 \(\tau_d=1/\Gamma=Q/(\pi f_0)\)。
粘性
由式 (6.1),\(2\Gamma_{\mathrm{vis}}=4\mu/(\rho R_0^2)\), \[\tag{6.4} \delta_{\mathrm{vis}}=\frac{4\mu}{\rho\omega_0R_0^2}.\] \(R_0=1\,\mathrm{mm}\) 时 \(\delta_{\mathrm{vis}}=2\times10^{-4}\),可忽略;它 \(\propto1/R_0\),只对微米泡重要。
声辐射
振荡的泡是一个单极子声源,会把能量辐射到远处。
设 \(s=s_a\cos\omega t\)。由式 (3.16),线性阶 \(p'(r,t)=\rho R_0^2\ddot s(t-r/c)/r\),振幅 \(\rho R_0^2\omega^2s_a/r\),均方值 \(\avg{p'^2}=\rho^2R_0^4\omega^4s_a^2/(2r^2)\)。穿过半径 \(r\) 球面的平均辐射功率 \[P_{\mathrm{rad}}=4\pi r^2\frac{\avg{p'^2}}{\rho c}=\frac{2\pi\rho R_0^4\omega^4s_a^2}{c} .\] 振动能 \(E=\frac12m_{\mathrm{eff}}\omega^2s_a^2=2\pi\rho R_0^3\omega^2s_a^2\)。于是 \[\tag{6.5} Q_{\mathrm{rad}}=\frac{\omega E}{P_{\mathrm{rad}}}=\frac{c}{\omega R_0}=\frac{1}{kR_0},\qquad \delta_{\mathrm{rad}}=\frac{\omega_0R_0}{c},\] \(k=\omega/c\) 为声波波数。
(1) 量纲:\(\omega R_0/c\) 无量纲。(2) 从泡壁受力看同一效应:把向外传播的解 (3.15) 在泡附近(\(r\ll\lambda\))对推迟时间展开到一阶,\(\phi=-F(t-r/c)/r\approx-F(t)/r+F'(t)/c\)。第二项在泡附近处处相同,不产生速度,却在泡壁处的非定常伯努利式中贡献 \(\pdv{\phi}{t}\supset F''/c\)。重做第 4.3 节的推导,RP 方程右边多出 \(\frac1c\dv{^2}{t^2}(R^2\dot R)\) 一项。线性化:\(F=R^2\dot R\approx R_0^3\dot\xi\),且对近似简谐的运动 \(\dddot\xi\approx-\omega_0^2\dot\xi\),此项成为 \(-\frac{R_0^3\omega_0^2}{c}\dot\xi\),除以 \(R_0^2\) 得附加阻尼 \(2\Gamma=\omega_0^2R_0/c\),即 \(\delta=\omega_0R_0/c\),与能流算法相同(verify_acoustic.py 还用数值积分核对了振幅衰减率)。这一修正项正是 Keller–Miksis 方程等“弱可压缩 RP 方程”中的领头项 [27]。(3) 由于 \(\omega_0R_0\approx20.5\,\mathrm{m/s}\) 几乎与 \(R_0\) 无关,\(\delta_{\mathrm{rad}}\approx20.5/1482=0.0138\),\(Q_{\mathrm{rad}}\approx72\)。\(kR_0\approx0.014\ll1\) 也说明气泡远小于它发出的声波波长——这正是“单极子”近似成立、以及在泡附近可当作不可压缩处理的理由。
热传导
由式 (6.3),\(p_g=p_{g0}(1-3\kappa_c\xi)\),其中 \(-3p_{g0}(\kappa_{\mathrm r}+\ii\kappa_{\mathrm i})\xi\) 的虚部项 \(-3\ii p_{g0}\kappa_{\mathrm i}\xi=-(3p_{g0}\kappa_{\mathrm i}/\omega)\dot\xi\)(对 \(\ee^{\ii\omega t}\) 有 \(\dot\xi=\ii\omega\xi\))与 \(\dot\xi\) 同相、方向相反,是阻尼力(\(\kappa_{\mathrm r}\)、\(\kappa_{\mathrm i}\) 为 \(\kappa_c\) 的实部与虚部)。它给出 \(2\Gamma_{\mathrm{th}}=3p_{g0}\kappa_{\mathrm i}/(\rho R_0^2\omega_0)\), \[\tag{6.6} \delta_{\mathrm{th}}=\frac{3p_{g0}\Im\kappa_c}{\rho R_0^2\omega_0^2}\approx\frac{\Im\kappa_c}{\Re\kappa_c}\qquad\Bigl(\frac{2\sigma}{R_0}\ll p_0\Bigr).\] 大泡(\(\abs{y}\gg1\))时附录 C 给出 \(\delta_{\mathrm{th}}\approx3(\gamma-1)\sqrt{D/(2\omega_0R_0^2)}\propto R_0^{-1/2}\)。
1 mm 空气泡:\(\delta_{\mathrm{th}}\approx0.026\)、\(\delta_{\mathrm{rad}}\approx0.014\)、\(\delta_{\mathrm{vis}}\approx0.0002\),\(Q\approx25\),振幅衰减时间 \(\tau_d=Q/(\pi f_0)\approx2.5\,\mathrm{ms}\),约 8 个周期——所以气泡的声音是短促的“叮”或“咕”。热阻尼的数值与 Devin 1959 年以来的经典结果一致 [30], [31]。
6.5 应用:流水声、雨声、超声造影
雨声
直径 0.8 mm–1.1 mm 的小雨滴撞击水面时会卷入气泡并激起衰减振荡,使水下雨噪声谱在约 15 kHz 出现一个显著的峰 [17], [18]。按式 (6.2) 反推,15 kHz 对应半径约 0.22 mm 的泡:\(R_0=0.22\,\mathrm{mm}\) 时 \(f_0\approx14.5\,\mathrm{kHz}\)(含热效应,\(\kappa_{\mathrm{eff}}=1.32\);这个半径是由频率反推的,不是直接测量值)。泡离水面很近,水面的影响见例题 2。这一特征峰为用水听器“听”雨提供了物理基础。
超声造影微泡
医学超声造影剂是直径约 1 μm–10 μm 的包膜微泡,它们在 MHz 超声下共振、强烈散射声波,使血流“变亮” [19]。式 (6.2) 给出半径 1.5 μm 的自由空气泡 \(f_0\approx2.4\,\mathrm{MHz}\),恰在医学超声频段。实际造影剂的包膜(脂质或蛋白)增加了弹簧的劲度与阻尼,需对式 (6.2) 加以修正。
7 含泡液体中的声速与“热巧克力效应”
用勺子敲一杯刚冲好的速溶热巧克力(或刚倒出的啤酒、刚打开的汽水),音调起初很低,随后几秒到几十秒内逐渐升高好几个八度。Crawford 1982 年把它称为“热巧克力效应”并给出了解释 [8]:他用软化的小槌敲击容器底部,测到的正是 \(1/4\) 波长的液柱纵向模式;用含溶解空气的热自来水灌满高量筒,音调在最初几秒内下降近三个八度(频率约降为 \(1/7\)),随后随气泡上浮逐渐回升;在咖啡杯里则只降约一个八度。原因是冲调或析出时形成的大量微小气泡使液体的声速大大降低,气泡上浮逸出后声速恢复。本节先推出含泡液体的声速公式。
7.1 Wood 公式
声波长远大于气泡半径与气泡间距,一个波长内有大量气泡,可把含泡液体看作均匀的“等效流体”;
声频率远低于气泡的 Minnaert 频率(\(f\ll f_0\)):气泡来得及随压强“准静态”地胀缩,泡壁惯性不起作用;
气泡随液体一起运动(无相对滑移),两相压强相等(略去 \(2\sigma/R\));
气体体积分数 \(\phi\)(空泡率)定义为单位体积混合物中气体所占体积。
密度:质量相加,\(\rho_m=\phi\rho_g+(1-\phi)\rho_l\)。
压缩率:压强增加 \(\dd p\) 时,体积 \(V\) 的混合物中气体部分 \(\phi V\) 缩小 \(\beta_g\phi V\dd p\),液体部分 \((1-\phi)V\) 缩小 \(\beta_l(1-\phi)V\dd p\),体积变化相加: \[\beta_m=-\frac{1}{V}\dv{V}{p}=\phi\beta_g+(1-\phi)\beta_l .\] 气体的压缩率由多方定律 \(pV^\kappa=\) 常量求得:\(\beta_g=-\frac1V\dv{V}{p}=\frac{1}{\kappa p}\);液体 \(\beta_l=1/(\rho_lc_l^2)\)。
声速:在一个波长内取一个包含大量气泡的小体积元。两相同速、同压,所以 (i) 它的惯性由总质量决定,线性欧拉方程为 \(\rho_m\pdv{u}{t}=-\pdv{p'}{z}\);(ii) 它的质量守恒仍是“密度变化率 \(=-\rho_m\times\) 体积流出率”,而压强变化 \(p'\) 引起的密度变化为 \(\rho'=\rho_m\beta_mp'\)(体积按 \(\beta_m\) 缩小、质量不变)。这与第 3.8 节的三步完全相同,只是 \(\rho_0\to\rho_m\)、\(1/c^2\to\rho_m\beta_m\),故 \(c_m^2=1/(\rho_m\beta_m)\)(式 (3.14))。
\[\tag{7.1} \frac{1}{\rho_mc_m^2}=\frac{\phi}{\rho_gc_g^2}+\frac{1-\phi}{\rho_lc_l^2},\qquad\rho_m=\phi\rho_g+(1-\phi)\rho_l,\] 其中气体项 \(\rho_gc_g^2=\kappa p\)。此式通常归于 A. B. Wood 1930 年的声学教科书 [32](亦有文献归于 Minnaert 1933 [4])。
“密度来自液体,弹性来自气体”
对 \(\phi\) 不太接近 0 或 1 的情形,\(\rho_g\ll\rho_l\) 使 \(\rho_m\approx(1-\phi)\rho_l\);而 \(\kappa p\approx10^{5}\,\mathrm{Pa}\) 比 \(\rho_lc_l^2\approx2.2\times10^{9}\,\mathrm{Pa}\) 小四个数量级,使压缩率由气体主导,\(\beta_m\approx\phi/(\kappa p)\)。于是 \[\tag{7.2} c_m\approx\sqrt{\frac{\kappa p}{\rho_l\,\phi(1-\phi)}},\qquad c_{\min}=c_m\Big|_{\phi=1/2}=2\sqrt{\frac{\kappa p}{\rho_l}}\approx\begin{cases}20\,\mathrm{m/s}, & \kappa=1\\ 24\,\mathrm{m/s}, & \kappa=1.4\end{cases}\quad(p=1\,\mathrm{atm}).\] 含泡水的声速可以比纯水(1482 m/s)和纯空气(343 m/s)都低得多——混合物同时拥有液体的“重”与气体的“软”。文献给出的数值为约 20 m/s [4] 与约 23 m/s [33],与上式一致。这与“质量–弹簧”图像一致:声速 \(\sim\sqrt{\text{弹性}/\text{惯性}}\)。
微量气体的巨大作用
气体项与液体项相等时 \(\phi^*=\kappa p/(\rho_lc_l^2)\approx4.6\times10^{-5}\):万分之零点五的气泡就使声速降到纯水的 \(1/\sqrt2\)。\(\phi=10^{-3}\) 时 \(c_m\approx310\,\mathrm{m/s}\),\(\phi=10^{-2}\) 时 \(c_m\approx100\,\mathrm{m/s}\)(\(\kappa=1\),图 7)。
热巧克力中的气泡半径为 10 μm–100 μm 量级,声频 \(\sim1\,\mathrm{kHz}\),\(\mathrm{Pe}=\omega R^2/D\lesssim3\),更接近等温,常取 \(\kappa\approx1\)。频率远低于这些泡的 Minnaert 频率(\(R=50\,\mu\mathrm{m}\) 时 \(f_0\approx60\,\mathrm{kHz}\)),假设 2 成立。若声频接近 \(f_0\),泡壁惯性使混合物出现强烈色散与衰减,Wood 公式不再适用 [33], [34]。
7.2 热巧克力效应的估算
敲击杯底(或杯壁)激发的主要是杯中液柱的纵向驻波:杯底(刚性)处液体速度为零,是压强波腹;液面(与空气接触,压强几乎不变)是压强波节。最低模式为 \(1/4\) 波长: \[\tag{7.3} f_1=\frac{c_m}{4h},\] \(h\) 为液深。对 \(h=8\,\mathrm{cm}\),纯水 \(f_1\approx4.6\,\mathrm{kHz}\);\(\phi=10^{-3}\) 时降到约 1.0 kHz(低 2.2 个八度);\(\phi=10^{-2}\) 时约 320 Hz(低 3.9 个八度)。气泡上浮离开时音调升高。注意这并不是“均匀的 \(\phi\) 随时间减小”:若泡大小相同、以同一速度上浮,泡群的下边界向上移动,杯底出现一层无泡的清液,而上层的空泡率不变——液柱变成两层介质。例题 3 将说明,此时基频主要由上层含泡液柱的长度决定,约为 \(c_m/[4(h-s)]\)(\(s\) 为清液层厚度),音调随含泡层缩短而升高;泡径不一时,各种泡依次离开,音调连续回升。实际杯壁并非绝对刚性,纯液体时的频率会比 \(c_l/4h\) 低,但气泡较多时 \(c_m\ll c_l\),杯壁的影响变小,式 (7.3) 的估算反而更可靠。
交互演示加载中……(需启用 JavaScript)
8 空化:从螺旋桨到枪虾
8.1 空化的起始:蒸气压、空化核与 Blake 阈值
液体中某处压强降到饱和蒸气压 \(p_v\) 以下时,液体中出现并迅速长大的充满蒸气(及少量溶解气体)的空腔。它与沸腾都是液体变成蒸气,区别在于:沸腾靠升温使 \(p_v\) 升到环境压强以上,空化靠流动或声场使局部压强降到 \(p_v\) 以下。
纯净无杂质的水理论上能承受很大的拉应力(负压)而不破裂:Frenkel 按分子间作用力的简单估计给出 \(10^9\)~\(10^{10}\,\mathrm{Pa}\) 量级的抗拉强度,而 19 世纪以来用玻璃毛细管做的实验只达到约 50~200 bar(\(5\times10^{6}\)~\(2\times10^{7}\,\mathrm{Pa}\)) [4]——差距来自液体中与器壁上的“薄弱点”。实际水中总有微米级的气泡核或固体表面缝隙中残存的气体,空化从这些“核”开始。一个核能否长大,可以用前面的工具判断。
半径 \(R_n\) 的泡核原本在压强 \(p_0\) 下平衡,泡内气压 \(p_{g0}=p_0-p_v+2\sigma/R_n\)。若周围液体压强缓慢降到 \(p_L\),泡准静态(等温)长大,平衡条件为 \[p_L=f(R)\equiv p_v+p_{g0}\frac{R_n^3}{R^3}-\frac{2\sigma}{R}.\] \(f(R)\) 在 \(R\to\infty\) 时趋于 \(p_v\),在小 \(R\) 时因气压项而很大,中间有一个极小值。令 \(f'(R)=-3p_{g0}R_n^3/R^4+2\sigma/R^2=0\): \[R_*=R_n\sqrt{\frac{3p_{g0}R_n}{2\sigma}},\qquad p_{\mathrm{Bl}}\equiv f(R_*)=p_v-\frac{4\sigma}{3R_*}.\] 若 \(p_L<p_{\mathrm{Bl}}\),就不存在平衡半径,泡将无限长大——这就是空化。\(p_{\mathrm{Bl}}\) 称为 Blake 阈值 [4]。对 \(2\sigma/R_n\gg p_0\) 的小核,\(p_{g0}\approx2\sigma/R_n\),\(R_*\approx\sqrt3R_n\),\(p_{\mathrm{Bl}}\approx p_v-\dfrac{4\sigma}{3\sqrt3R_n}=p_v-0.77\dfrac{\sigma}{R_n}\)。
数值(fig_shrimp.py 精确求极小):\(R_n=1\,\mu\mathrm{m}\) 时 \(p_{\mathrm{Bl}}=-41\,\mathrm{kPa}\)(需要拉应力);\(R_n=10\,\mu\mathrm{m}\) 时 \(p_{\mathrm{Bl}}=0.3\,\mathrm{kPa}\),略低于 \(p_v\) 即可。所以液体越“干净”(核越小),越难空化;含有较大气核的普通水只要压强降到 \(p_v\) 附近就会空化。
8.2 空化数
对速度为 \(U\)、远处压强为 \(p_\infty\) 的绕流,定常伯努利方程(式 (3.9) 去掉 \(\pdv{\phi}{t}\))给出流线上 \(p+\frac12\rho u^2=p_\infty+\frac12\rho U^2\)。流速增大处压强降低。
压强系数 \(C_p\equiv(p-p_\infty)/(\frac12\rho U^2)\) 只取决于物体形状(对给定的流动形态)。空化数 \[\tag{8.1} \sigma_{\mathrm c}\equiv\frac{p_\infty-p_v}{\frac12\rho U^2}.\] 表面最低压强 \(p_{\min}=p_\infty+C_{p,\min}\frac12\rho U^2\) 降到 \(p_v\) 的条件是 \(\sigma_{\mathrm c}\le-C_{p,\min}\):流速越高、环境压强越低,\(\sigma_{\mathrm c}\) 越小,越容易空化。
理想流体以 \(U\) 绕半径 \(a\) 的球流动,速度势 \(\phi=U\cos\theta\,(r+a^3/(2r^2))\)(第一项 \(Ur\cos\theta=Uz\) 是均匀流;第二项 \(\frac{Ua^3}{2}\frac{\cos\theta}{r^2}=\frac{Ua^3}{2}\frac{z}{r^3}=-\frac{Ua^3}{2}\pdv{}{z}\Bigl(\frac1r\Bigr)\) 是点源势 \(1/r\) 沿 \(z\) 方向的导数,称为偶极子势,相当于一对相距无限近的正、负点源。点源势在原点外满足不可压缩条件(第 3.2 节的流量守恒),对 \(z\) 求导不破坏这一点,所以叠加后的流动不可压缩。再验证球面上法向速度 \(\pdv{\phi}{r}\big|_a=U\cos\theta\,(1-a^3/a^3)=0\),即流体不穿过球面)。球面上切向速度 \(\frac1r\pdv{\phi}{\theta}\big|_a=-\frac32U\sin\theta\),赤道处最大 \(\frac32U\),\(C_{p,\min}=1-\frac94=-1.25\)。故在理想化条件下 \(\sigma_{\mathrm c}<1.25\) 时球面赤道附近开始空化。螺旋桨叶尖速度可达数十 m/s:叶尖线速度 30 m/s、水深 3 m 处 \(\sigma_{\mathrm c}=(1.3\times10^{5}-2.3\times10^{3})/(\frac12\times1000\times30^2)\approx0.28\),远低于 1,叶片低压面极易空化。
8.3 枪虾:用空化泡捕猎
鼓虾科(Alpheidae)的枪虾有一只特别巨大的螯。螯的活动指(dactyl)上的“塞子”快速插入固定指上的“凹槽”,把水从槽中挤出,形成高速射流:Versluis 等由高速图像中空化泡的推进速度估计射流约 25 m/s [15], [35],真实尺寸螯的数值模拟给出约 30 m/s,螯的闭合历时约 0.3 ms [35]。长期以来人们以为响声来自两指相撞;Versluis 等 2000 年用高速摄影与水听器同步测量证明:响声发生在空化泡塌缩的时刻,而不是螯闭合的时刻,并用 Rayleigh–Plesset 型方程定量重现了泡半径随时间的变化和声信号 [15]。距螯 4 cm 处的声压约 80 kPa——这一数值我们只见到二手转述 [36],未能确认它指正峰值还是峰峰值(原文未能获取)。2001 年他们又发现塌缩时有极短的闪光(“虾致发光”),若光来自热辐射,需要至少 5000 K 的温度 [16],每次约 \(10^4\) 个光子 [12]。按 5 倍放大的仿生螯装置得到的空化泡最大半径约 5.6 mm、空化数 \(0.14\)~\(0.30\)(与真实枪虾可比) [37]。数值模拟表明,空化主要发生在射流前端卷起的涡环(像烟圈那样首尾闭合成环的涡旋;流体绕环的截面中心高速旋转,这个截面中心附近称为涡核)核心的低压区中,核心压强亏损平均约 \(2\times10^{5}\,\mathrm{Pa}\)、最大可达 \(10^{6}\,\mathrm{Pa}\) [35]。
低压从哪里来?——不能对射流直接用伯努利方程
一个诱人的错误是:取周围静止海水(\(p_0\)、\(u=0\))与射流(速度 \(u\))列伯努利方程,得“射流区压强 \(p_0-\frac12\rho u^2\)”。但定常伯努利方程只沿同一条流线成立:射流的水来自被塞子挤压的凹槽,那里的总压高于 \(p_0\);射流离开喷口后是自由射流,其内部静压与周围相同,约为 \(p_0\)。两者不在同一条流线上,不能这样连接。真正的低压来自射流前端卷起的涡环:流体绕涡核高速旋转,维持圆周运动需要指向核心的压强梯度(\(\dd p/\dd r=\rho v_\theta^2/r\)),于是核心压强低于外部。取最简单的涡模型(Rankine 涡:半径 \(a\) 以内像刚体一样转动,\(v_\theta=\Omega r\);以外 \(v_\theta=\Omega a^2/r\)),对 \(\dd p/\dd r=\rho v_\theta^2/r\) 从核心积分到无穷远,外区贡献压强亏损 \(\frac12\rho v_{\max}^2\)、核内再贡献 \(\frac12\rho v_{\max}^2\),核心压强为 \(p_0-\rho v_{\max}^2\)(\(v_{\max}=\Omega a\);逐步推导见例题 4(1b))。取 \(v_{\max}=\lambda u\),空化(\(p_{\mathrm{core}}<p_v\))要求 \[\tag{8.2} \lambda^2>\frac{p_0-p_v}{\rho u^2}=\frac{\sigma_{\mathrm c}}{2},\qquad \sigma_{\mathrm c}\equiv\frac{p_0-p_v}{\frac12\rho u^2}.\] 取海水密度 \(\rho={1025\,}\mathrm{kg/m^3}\)、\(u={25\,}\mathrm{m/s}\),\(\sigma_{\mathrm c}\approx0.31\),空化起始只需 \(\lambda>0.39\)(若用淡水 \(\rho=998\),则 0.32 与 0.40)。但这只是核心压强低于 \(p_v\) 的静态阈值:泡要在百微秒内长到毫米级,低压还须足够深——在图 8 的模型中,\(\lambda=0.40\) 时泡只长到 0.08 mm,\(\lambda=0.5\) 时 0.39 mm,要到 \(\lambda\approx0.7\) 才长到 1.1 mm;\(\lambda=1/\sqrt2\) 时核心压强亏损恰为 \(\frac12\rho u^2\approx3.2\times10^{5}\,\mathrm{Pa}\),与模拟给出的量级一致。所以动压 \(\frac12\rho u^2\) 仍是低压的正确量级,但它出现在涡核中,而不是“射流区”。
用本讲的工具可以把整个过程拆成几步(图 8):
空化:涡核压强 \(p_0-\rho\lambda^2u^2\) 远低于 \(p_v\)(\(\lambda=1/\sqrt2\) 时约 \(-2.2\times10^{5}\,\mathrm{Pa}\),液体处于拉伸状态),由 Blake 阈值,半径大于约 0.25 μm 的核都会失稳长大。
生长:涡核低压只维持一段与螯闭合时间同量级的时间(下面的模型取 \(T_j={133\,}\mathrm{\mu s}\),是模型假设),泡在这期间被“拉”大;低压消失后泡还会凭惯性继续长大一段时间,到 \(R_{\max}\) 后在 \(p_0-p_v\) 作用下开始 Rayleigh 塌缩。
塌缩与发声:塌缩历时 \(\approx0.915R_{\max}\sqrt{\rho/(p_0-p_v)}\);泡体积的剧烈加速度按式 (3.16) 辐射一个尖锐的压强脉冲;泡内残存的“缓冲气体”(来不及凝结的蒸气与扩散进入的空气)被绝热压缩到数千 K,发出闪光。
交互演示加载中……(需启用 JavaScript)
9 决赛风格例题
以下各题的数值答案均已用脚本 verify_examples.py 重算。除特别说明,本节取水的物性为 \(\rho=998\,\mathrm{kg/m^3}\)、\(c=1482\,\mathrm{m/s}\)、\(\sigma=0.0728\,\mathrm{N/m}\)、\(\mu=1.00\times10^{-3}\,\mathrm{Pa\,s}\),空气 \(\gamma=1.40\)。
9.1 例题 1:螺旋桨空泡的塌缩、压强峰与气体缓冲
某船螺旋桨叶背的低压区中产生一个蒸气空泡,随后被水流带到压强为 \(p_\infty=2.0\times10^{5}\,\mathrm{Pa}\) 的区域。设空泡此时达到最大半径 \(R_0=2.0\,\mathrm{mm}\) 并静止,泡内蒸气压 \(p_v=2.3\,\mathrm{kPa}\) 保持不变(蒸气在塌缩中随时凝结)。忽略表面张力、粘性与重力,泡保持球形。
先设泡内只有蒸气。求塌缩时间 \(t_c\);求 \(R=R_0/10\) 时液体的动能,并说明它与哪一部分功相等。
求 \(R=R_0/10\) 与 \(R=R_0/20\) 时的泡壁速率。在多大半径处泡壁速率达到水中声速?此结论对 RP 方程的适用性意味着什么?
求 \(R=R_0/20\) 时液体中的最大压强及其位置,与常见金属的屈服强度( × 108 Pa~ × 109 Pa 量级)比较。
实际空泡中溶有空气。设 \(R=R_0\) 时泡内还有分压 \(p_{g0}=5.0\,\mathrm{kPa}\) 的空气,其压缩为绝热。(a) 写出泡壁速率与半径的关系,求最大泡壁速率及其所在半径;(b) 求最小半径 \(R_{\min}\) 与此时泡内空气的温度(初温 293 K);(c) 检验 (a) 中速率与声速之比,评价模型自洽性。
空泡往往紧挨着叶片表面。回到只有蒸气的空泡,设泡心到平直刚性壁面的距离为 \(d\)(\(R_0<d\))。(a) 说明:在“半无限液体 + 刚性壁”中,泡的流场可以用“真实泡 + 关于壁面对称的一个相同的镜像泡”在无界液体中的流场来表示,壁面上法向速度自动为零(镜像法)。(b) 在 \(R\ll2d\) 的近似下,镜像泡的速度势在真实泡附近近似为常量。写出它对真实泡壁处非定常伯努利方程的贡献,导出修正后的运动方程,并验证它等价于液体动能 \(K'=2\pi\rho R^3\dot R^2\bigl(1+\frac{R}{2d}\bigr)\) 的能量守恒。(c) 保留 \(R_0/d\) 的一阶,求塌缩时间相对于无界液体的修正,答案用 Beta 函数表示并给出数值系数。(d) 对 \(d={6.0\,}\mathrm{mm}\) 给出 \(t_c\)。
(1) 蒸气压恒定,相当于把 \(p_\infty\) 换成 \(\Delta p=p_\infty-p_v=1.977\times10^{5}\,\mathrm{Pa}\)。时间单位 \(\tau=R_0\sqrt{\rho/\Delta p}=142.1\,\mu\mathrm{s}\),由式 (5.3) \[t_c=0.9147\,\tau=130\,\mu\mathrm{s}.\] 能量:外压(扣除泡内蒸气压的反推)做功 \(\frac43\pi\Delta p(R_0^3-R^3)\) 全部变为液体动能。\(R=R_0/10\) 时 \[K=\tfrac43\pi\Delta p R_0^3\bigl(1-10^{-3}\bigr)=6.62\,\mathrm{mJ},\] 即总可用能 \(\frac43\pi\Delta pR_0^3\) 的 \(99.9\%\) 已转化为动能,但这些动能集中在泡附近越来越少的液体中。评分要点:识别 \(\Delta p=p_\infty-p_v\);能量守恒;代入数值。
(2) 由式 (5.1)(\(p_\infty\to\Delta p\)),\(\abs{\dot R}=\sqrt{\frac{2\Delta p}{3\rho}(R_0^3/R^3-1)}\):\(R=R_0/10\) 时 \(363\,\mathrm{m/s}\),\(R=R_0/20\) 时 \(1028\,\mathrm{m/s}\)。由式 (5.8),\(\abs{\dot R}=c\) 于 \(R_{\mathrm M}=R_0(2\Delta p/3\rho c^2)^{1/3}=0.039R_0=78\,\mu\mathrm{m}\)。所以在 \(R\lesssim R_0/25\) 后不可压缩的 RP 方程不再可靠,塌缩末段需计入液体压缩性(能量将以冲击波形式辐射)。
(3) \(z=R_0^3/R^3=8000\),由式 (5.7) \[p_{\max}=\Delta p\Bigl[1+\frac{(z-4)^{4/3}}{4^{4/3}(z-1)^{1/3}}\Bigr]=249\,\mathrm{MPa},\qquad r_m=\Bigl[\frac{4(z-1)}{z-4}\Bigr]^{1/3}R=1.588R=159\,\mu\mathrm{m},\] 在距泡壁约 59 μm 的液体中。它已达到金属屈服强度的量级。若塌缩发生在叶片表面附近,这样的压强脉冲(以及近壁不对称塌缩形成的微射流 [25], [26],见 (5))反复作用,就造成疲劳剥蚀。注意:(2) 表明此时已接近可压缩区,这里的数值应理解为量级估计。
(4a) 绝热压缩时外界对空气做功 \(W_g=\frac{p_{g0}V_0}{\gamma-1}\bigl(z^{\gamma-1}-1\bigr)\)(\(z=V_0/V=R_0^3/R^3\)),能量守恒除以 \(2\pi\rho R^3\): \[\dot R^2=\frac{2}{3\rho}\Bigl[\Delta p\,(z-1)-\frac{p_{g0}}{\gamma-1}\bigl(z^\gamma-z\bigr)\Bigr].\] 对 \(z\) 求导令其为零:\(\Delta p-\frac{p_{g0}}{\gamma-1}(\gamma z^{\gamma-1}-1)=0\), \[z_*^{\gamma-1}=\frac{1}{\gamma}\Bigl[\frac{(\gamma-1)\Delta p}{p_{g0}}+1\Bigr]=12.01,\qquad z_*=500,\qquad R_*=z_*^{-1/3}R_0=0.126R_0=252\,\mu\mathrm{m},\] 代回得 \(\abs{\dot R}_{\max}=141\,\mathrm{m/s}\)。数值积分 RP 方程给出相同结果(相对差 \(<10^{-3}\))。
(4b) 由式 (5.9)(\(p_\infty\to\Delta p\)),\(p_{g0}/\Delta p=0.0253\),数值求根得 \(x=R_{\min}/R_0=0.0952\),\(R_{\min}=190\,\mu\mathrm{m}\)(近似式 (5.10) 给出 \(0.100\),偏大 \(5\%\),因为 \(p_{g0}/[(\gamma-1)\Delta p]=0.063\) 已不算很小)。温度 \(T_{\max}=T_0x^{-3(\gamma-1)}=4.9\times10^{3}\,\mathrm{K}\),此时泡内气压 \(p_{g0}x^{-3\gamma}\approx97\,\mathrm{MPa}\)。
(4c) \(\abs{\dot R}_{\max}/c=0.095\),可压缩修正为 \(10\%\) 量级,模型基本自洽。对比 (2):仅仅 5 kPa 的空气就把最大速率从“超声速”压到 141 m/s,最小半径停在约 \(R_0/10\)——不凝结气体的缓冲作用极强。
(5a) 无界液体中两个相同的泡关于平面 \(z=0\) 对称地胀缩,由对称性,平面上任一点两泡产生的速度的法向分量大小相等、方向相反,合成为零——这正是刚性壁的边界条件。所以在 \(z>0\) 的半空间里,“泡 + 镜像泡”的流场就是“泡 + 壁”的流场(势流满足的方程是线性的,两个解可以叠加)。
(5b) 镜像泡的速度势 \(\phi_{\mathrm{im}}=-R^2\dot R/r'\),\(r'\) 为到镜像泡心的距离;在真实泡壁上 \(r'=2d+\order{R}\),故 \(\phi_{\mathrm{im}}\approx-R^2\dot R/(2d)\),近似为空间常量——它不贡献泡附近的速度,但贡献 \(\pdv{\phi_{\mathrm{im}}}{t}=-\frac{1}{2d}\dv{}{t}(R^2\dot R)=-\frac{R^2\ddot R+2R\dot R^2}{2d}\)。在第 4.3 节的推导中加入这一项: \[\Bigl(R+\frac{R^2}{2d}\Bigr)\ddot R+\Bigl(\frac32+\frac{R}{d}\Bigr)\dot R^2=\frac{p_L-p_\infty}{\rho}.\] 能量验证:\(\dv{}{t}\Bigl[2\pi\rho\Bigl(R^3+\frac{R^4}{2d}\Bigr)\dot R^2\Bigr]=4\pi\rho R^2\dot R\Bigl[\Bigl(R+\frac{R^2}{2d}\Bigr)\ddot R+\Bigl(\frac32+\frac Rd\Bigr)\dot R^2\Bigr]\),正好是上式左边乘 \(4\pi\rho R^2\dot R\),所以上式等价于 \(\dv{K'}{t}=(p_L-p_\infty)4\pi R^2\dot R\)。物理意义:壁面挡住了一半空间,泡附近的液体只能从一侧补充,被带动的液体更多,有效质量增加 \(R/(2d)\) 倍。
(5c) 能量守恒 \(K'=\frac43\pi\Delta p(R_0^3-R^3)\),令 \(\varepsilon=R_0/(2d)\)、\(x=R/R_0\):\(\abs{\dot R}=\sqrt{\frac{2\Delta p}{3\rho}\frac{x^{-3}-1}{1+\varepsilon x}}\), \[t_c'=\tau\sqrt{\frac32}\int_0^1\frac{x^{3/2}\sqrt{1+\varepsilon x}}{\sqrt{1-x^3}}\dd x
\approx\tau\sqrt{\frac32}\Bigl[\frac13\mathrm B\Bigl(\frac56,\frac12\Bigr)+\frac\varepsilon2\cdot\frac13\mathrm B\Bigl(\frac76,\frac12\Bigr)\Bigr],\] 第二项的积分 \(\int_0^1x^{5/2}(1-x^3)^{-1/2}\dd x\) 用与附录 B 相同的换元 \(s=x^3\) 化为 \(\frac13\mathrm B(\frac76,\frac12)\)。于是 \[\frac{t_c'}{t_c}=1+\frac{\mathrm B(7/6,1/2)}{2\,\mathrm B(5/6,1/2)}\,\varepsilon=1+0.203\,\frac{R_0}{d}.\] 数值核对:直接积分上式与数值积分修正后的运动方程一致;\(d/R_0=3\) 时精确比值 1.0655,一阶式 1.068(verify_examples.py)。
(5d) \(R_0/d=1/3\):\(t_c'=130\,\mu\mathrm{s}\times1.065\approx139\,\mu\mathrm{s}\)。近壁塌缩变慢;更重要的是镜像泡的存在破坏了球对称——泡朝向壁面一侧与背离壁面一侧受到的“补充”不同,这正是不对称塌缩与指向壁面的微射流的起源 [25], [26],已超出球对称模型。
9.2 例题 2:雨声——近水面气泡的振荡与辐射
直径约 1 mm 的雨滴落入平静水面时,常在水面下卷入一个半径 \(R_0=0.22\,\mathrm{mm}\) 的气泡,泡心离水面的深度约 \(d=1.5\,\mathrm{mm}\)(本题设定值)。气泡脱离瞬间的形变使它以径向相对振幅 \(\varepsilon=s_a/R_0=1.0\%\) 开始自由振荡。水面压强 \(p_0=1.013\times10^{5}\,\mathrm{Pa}\),取 \(\kappa_{\mathrm{eff}}=1.32\)(热传导理论值)、热阻尼 \(\delta_{\mathrm{th}}=0.050\)。
先把泡当作在无界水中。计算振荡频率(含表面张力),并与绝热、忽略表面张力的结果比较;求有效质量、等效劲度与初始振动能 \(E\)。
仍在无界水中:求距泡 \(r=1\,\mathrm{m}\) 处的声压振幅与声压级 \(L=20\lg\bigl(p_{\mathrm{rms}}/1\,\mu\mathrm{Pa}\bigr)\),以及 \(\delta_{\mathrm{rad}}\)、\(\delta_{\mathrm{vis}}\)、\(Q\) 与振幅衰减时间 \(\tau_d\)。
现在计入水面。水面上方是空气,声压几乎为零(“压强释放”边界)。(a) 由线性化伯努利关系 \(p'=-\rho\pdv{\phi}{t}\) 说明:水面处应有 \(\phi=0\),可以用关于水面对称、振动反相的镜像泡(速度势反号)来满足它。(b) 在 \(R_0\ll2d\) 下,求镜像对泡附近流场的影响,说明有效质量变为 \(m_{\mathrm{eff}}(1-R_0/2d)\),并求频率的相对变化。
(a) 在远场(\(r\gg d\)),泡与镜像到观测点的路程差为 \(2d\cos\theta\)(\(\theta\) 从竖直向下量起)。在 \(kd\ll1\)(\(k=\omega/c\))下求声压的方向分布;(b) 对下半空间积分,证明辐射功率是同一个泡在无界水中辐射功率的 \(\frac23(kd)^2\) 倍;(c) 求新的 \(\delta_{\mathrm{rad}}\)、\(Q\),以及正下方 1 m 处的声压级。
与雨滴的动能(末速约 4 m/s [38])相比,求雨滴动能转化为声能的效率,并与无界水中的结果比较。
(1) \(p_{g0}=p_0+2\sigma/R_0=p_0+662\,\mathrm{Pa}\),由式 (6.1),\(\omega_0^2=(3\kappa p_{g0}-2\sigma/R_0)/(\rho R_0^2)\),\(f_0=14.54\,\mathrm{kHz}\);绝热且忽略表面张力为 \(14.94\,\mathrm{kHz}\)。表面张力使 \(\omega_0^2\) 相对增加 \(\frac{3\kappa-1}{3\kappa}\frac{2\sigma}{R_0p_0}=0.49\%\)(频率 \(+0.24\%\)),热效应使频率降低约 \(3\%\),影响更大。\(m_{\mathrm{eff}}=4\pi\rho R_0^3=1.34\times10^{-7}\,\mathrm{kg}\),\(k_s=m_{\mathrm{eff}}\omega_0^2=1.11\times10^{3}\,\mathrm{N/m}\),\(s_a=2.2\,\mu\mathrm{m}\),\(E=\frac12k_ss_a^2=2.7\times10^{-9}\,\mathrm{J}\)。
(2) 由式 (3.16),\(p_a=\rho R_0^2\omega_0^2s_a/r=0.89\,\mathrm{Pa}\),\(L=20\lg\frac{0.89/\sqrt2}{10^{-6}}\approx116\,\mathrm{dB}\)(\(r=1\,\mathrm m\) 远大于波长 \(\approx10\,\mathrm{cm}\),是远场)。\(\delta_{\mathrm{rad}}=\omega_0R_0/c=0.0136\),\(\delta_{\mathrm{vis}}=4\mu/(\rho\omega_0R_0^2)=0.0009\),\(\delta=0.0645\),\(Q=15.5\),\(\tau_d=Q/(\pi f_0)=0.34\,\mathrm{ms}\),约 5 个周期。
(3a) 水面上 \(p'=0\) 对任意时刻成立,即 \(\pdv{\phi}{t}=0\);扰动从零开始,故 \(\phi=0\)。镜像泡的速度势与真实泡大小相等、符号相反,在对称面上两者处处抵消,\(\phi=0\) 自动满足(与例题 1 刚性壁的同号镜像相反)。
(3b) 镜像的势 \(+R^2\dot R/r'\) 在真实泡附近近似为常量 \(+R^2\dot R/(2d)\),代入伯努利方程,与例题 1(5b) 同理但符号相反:线性阶 \(R_0^2\bigl(1-\frac{R_0}{2d}\bigr)\ddot\xi=\dots\),有效质量变为 \(m_{\mathrm{eff}}(1-R_0/2d)\)——水面一侧液体“轻”,被带动的液体变少。弹簧(气体)不变,故 \[\frac{\omega'}{\omega_0}=\Bigl(1-\frac{R_0}{2d}\Bigr)^{-1/2}=1.039,\qquad f'=15.10\,\mathrm{kHz}.\] 振动能 \(E=\frac12k_ss_a^2\) 不变。
(4a) 远场中两个声源的声压幅值相等、相位相反,到达观测点的相位差为 \(k\cdot2d\cos\theta\): \[p'\propto\frac{1}{r}\Bigl[\ee^{\ii(\omega t-kr)}-\ee^{\ii(\omega t-kr-2kd\cos\theta)}\Bigr],\qquad \abs{p'}=p_{\mathrm m}\,\abs{1-\ee^{-2\ii kd\cos\theta}}\approx p_{\mathrm m}\cdot2kd\cos\theta ,\] 其中 \(p_{\mathrm m}\) 为单独一个泡的声压幅值(以泡与镜像的中点为 \(r\) 的起点,近似到 \(kd\) 的一阶)。这是偶极子辐射:正下方最强,水平方向为零。
(4b) 平均声强 \(\propto\abs{p'}^2\)。单个泡在无界水中的功率 \(\propto p_{\mathrm m}^2\cdot4\pi\)(对全部立体角);现在只有下半空间有水: \[\frac{P'}{P_{\mathrm m}}=\frac{\int_{\text{下半空间}}4k^2d^2\cos^2\theta\,\dd\Omega}{4\pi}=\frac{4k^2d^2\cdot2\pi/3}{4\pi}=\frac23(kd)^2 ,\] 用了 \(\int_0^{\pi/2}\cos^2\theta\sin\theta\,\dd\theta\cdot2\pi=2\pi/3\)。数值积分精确的 \(\abs{1-\ee^{-2\ii kd\cos\theta}}^2\) 与此相差 \(0.2\%\)(verify_examples.py)。
(4c) \(kd=\omega'd/c=0.096\),功率比 \(6.2\times10^{-3}\)。\(\delta_{\mathrm{rad}}'=\frac{\omega'R_0}{c}\cdot\frac23(kd)^2=8.7\times10^{-5}\),\(\delta_{\mathrm{vis}}'=0.0009\),\(Q'=1/(0.050+0.0009+0.0001)=19.6\)——辐射几乎不再是损耗机制。正下方 1 m:\(p_{\mathrm m}=\rho R_0^2\omega'^2s_a/r=0.96\,\mathrm{Pa}\),\(\abs{p'}=2kd\,p_{\mathrm m}=0.18\,\mathrm{Pa}\),\(L\approx102\,\mathrm{dB}\),比无界水中低约 14 dB。
(5) 雨滴质量 \(\frac43\pi(0.5\,\mathrm{mm})^3\rho=0.52\,\mathrm{mg}\),动能 \(4.2\times10^{-6}\,\mathrm{J}\)。无界水中辐射份额 \(\delta_{\mathrm{rad}}/\delta=0.21\),声能 \(5.7\times10^{-10}\,\mathrm J\),效率 \(1.4\times10^{-4}\);计入水面后份额 \(\delta_{\mathrm{rad}}'/\delta'=1.7\times10^{-3}\),声能 \(4.6\times10^{-12}\,\mathrm J\),效率约 \(10^{-6}\)。近水面的“镜像抵消”使声源效率下降两个数量级;但每个泡都在同一频率附近“唱”,大量雨滴叠加仍形成显著的谱峰。
9.3 例题 3:热巧克力效应——两层液柱
杯中液深 \(h=8.0\,\mathrm{cm}\),热饮温度约 70 °C,按同温度的水取物性:\(\rho=977.8\,\mathrm{kg/m^3}\),\(c_l=1555\,\mathrm{m/s}\),粘度 \(\mu=0.404\,\mathrm{mPa\,s}\),饱和蒸气压 \(p_v=31.2\,\mathrm{kPa}\);液面为大气压 \(p=1.013\times10^{5}\,\mathrm{Pa}\)。刚搅拌后敲杯,音调比静置很久后(无气泡时)低整整三个八度。设气泡半径均为 \(R_b=30\,\mu\mathrm{m}\),泡内为空气与饱和水蒸气,温度始终等于液温。
泡内总压强 \(p\) 中有 \(p_v\) 是蒸气分压。设压强变化时蒸气随时凝结或蒸发、使蒸气分压保持 \(p_v\),空气等温压缩。证明气体部分的压缩率为 \(\beta_g=1/(p-p_v)\),并说明它比纯空气泡“更软”的原因。
设杯壁绝对刚性。由 Wood 公式推导小 \(\phi\) 下 \(f_1(\phi)/f_1(0)\) 的表达式,并求刚搅拌后的空泡率。
实测静置后(无气泡)基频为 \(f_\infty=1.60\,\mathrm{kHz}\),远低于 (2) 中的刚性杯值。把杯壁的弹性等效为液体压缩率的增加,用 \(f_\infty\) 定义液柱的等效声速 \(c_{\mathrm{eff}}\),重新估算 \(\phi\)。哪个估算更可信?
用 Stokes 公式(半径 \(R_b\) 的刚性球以速度 \(v\) 匀速运动受阻力 \(6\pi\mu R_bv\),适用于雷诺数 \(\mathrm{Re}=2\rho vR_b/\mu\ll1\);表面活性物质使小气泡表现得像刚性球)求上浮速度。\(t\) 时刻杯底已有厚 \(s=vt\) 的清液层,其上为厚 \(h-s\)、空泡率不变的含泡层。(a) 以杯底为 \(z=0\),写出两层中压强驻波的形式(杯底 \(\pdv{p'}{z}=0\),液面 \(p'=0\)),利用界面处压强连续、法向速度连续(由线性欧拉方程,\(\frac1\rho\pdv{p'}{z}\) 连续),导出本征频率满足的方程;(b) 当清液层的 \(1/4\) 波长频率远高于含泡层时,证明 \(f\approx c_m/[4(h-s)]\);(c) 求 \(t=0\)、4、8、12、14、16 s 的基频,说明音调如何回升。
检验 Wood 公式成立的两个条件。
(1) 体积 \(V\) 的泡中空气的物质的量不变、等温:\((p-p_v)V=\) 常量,故 \(\dd V/V=-\dd p/(p-p_v)\),\(\beta_g=1/(p-p_v)=1/(70.1\,\mathrm{kPa})\)。压强升高时一部分蒸气凝结成水、体积随之消失,没有“抵抗”,只有空气提供弹性,所以比同体积的纯空气泡(\(\beta=1/p\))更易压缩。(这要求蒸发、凝结足够快;实际中这是一个近似。)
(2) 刚性杯:\(f_1(0)=c_l/(4h)=4.86\,\mathrm{kHz}\)。小 \(\phi\) 时 \(\rho_m\approx\rho\)(相对误差 \(\phi\)),Wood 公式给出 \[\frac{1}{c_m^2}\approx\rho\Bigl(\frac{1}{\rho c_l^2}+\phi\beta_g\Bigr) \ \Longrightarrow\ \frac{f_1(\phi)}{f_1(0)}=\frac{c_m}{c_l}\approx\Bigl(1+\frac{\phi}{\phi^*}\Bigr)^{-1/2},\qquad\phi^*=\frac{p-p_v}{\rho c_l^2}=2.97\times10^{-5}.\] 低三个八度即 \(c_m/c_l=1/8\),\(\phi=63\phi^*=1.9\times10^{-3}\)。
(3) 压强升高时杯子略微胀大,效果等同于液体更易压缩。把它并入液体压缩率,\(c_{\mathrm{eff}}=4hf_\infty=512\,\mathrm{m/s}\),同样的推导给出 \(\phi\approx63(p-p_v)/(\rho c_{\mathrm{eff}}^2)=0.017\)(完整 Wood 公式保留 \(1-\phi\) 因子的精确解 \(0.0175\)),比 (2) 大约 9 倍。(3) 更可信:“三个八度”是用同一只杯子的无泡频率作参照测出的比值,只有用这个实测参照频率,比值才只反映气泡的贡献。(假设杯壁的等效压缩率不随 \(\phi\) 改变。)
(4) Stokes 匀速:\(\frac43\pi R_b^3\rho g=6\pi\mu R_bv\),\(v=\frac{2\rho gR_b^2}{9\mu}=4.75\,\mathrm{mm/s}\),\(\mathrm{Re}=0.69\),基本满足;清空整杯需 \(h/v=17\,\mathrm s\)。(若误用室温粘度 \(1.0\,\mathrm{mPa\,s}\),会得到慢 2.5 倍的结果。)
(a) 记清液层(\(0<z<s\))声速 \(c_2=c_{\mathrm{eff}}\)、密度 \(\rho_2=\rho\),含泡层(\(s<z<h\))声速 \(c_1=c_m=c_{\mathrm{eff}}/8=64\,\mathrm{m/s}\)、密度 \(\rho_1\approx\rho\)。满足两端条件的驻波为 \[p_2=A\cos k_2z,\qquad p_1=B\sin k_1(h-z),\qquad k_i=\omega/c_i .\] 线性欧拉方程 \(\rho\pdv{u}{t}=-\pdv{p'}{z}\) 表明法向速度连续等价于 \(\frac1\rho\pdv{p'}{z}\) 连续。在 \(z=s\): \[A\cos k_2s=B\sin k_1h_1,\qquad \frac{Ak_2}{\rho_2}\sin k_2s=\frac{Bk_1}{\rho_1}\cos k_1h_1\qquad(h_1\equiv h-s),\] 两式相除得本征方程 \[\tan(k_1h_1)\,\tan(k_2s)=\frac{\rho_2c_2}{\rho_1c_1}\equiv\frac{Z_2}{Z_1}\approx8.1 .\] \(Z=\rho c\) 称为介质的声阻抗。
(b) 若 \(k_2s\ll1\)(清液层远短于其中的四分之一波长),\(\tan k_2s\approx k_2s\),于是 \(\tan k_1h_1=\frac{Z_2}{Z_1k_2s}\gg1\),\(k_1h_1\approx\frac\pi2\),即 \(f\approx\frac{c_m}{4(h-s)}\):清液层“又重又硬”,对含泡层来说几乎就是新的刚性杯底。
(c) 数值求本征方程的最低根(verify_examples.py):
| \(t\)/s | 0 | 4 | 8 | 12 | 14 | 16 |
|---|---|---|---|---|---|---|
| \(s\)/cm | 0 | 1.90 | 3.80 | 5.70 | 6.65 | 7.60 |
| \(f\)(两层精确)/Hz | 200 | 261 | 376 | 668 | 1073 | 1589 |
| \(c_m/[4(h-s)]\) /Hz | 200 | 262 | 381 | 695 | 1183 | 3978 |
音调先缓慢、后加速地升高,在含泡层快要消失时趋于无泡值 1.60 kHz。近似式在 \(s\lesssim h/2\) 时误差 \(<2\%\),在含泡层很薄时失效(此时条件 (b) 不再成立)。与“\(\phi\) 均匀减小”的图像不同,这里上层的空泡率始终不变,音调回升完全来自含泡层变短。
(5) (a) 长波长:基频波长 \(\ge4h=32\,\mathrm{cm}\);泡数密度 \(n=\phi/(\frac43\pi R_b^3)\approx1.6\times10^{11}\,\mathrm{m^{-3}}\)(用 \(\phi=0.0175\)),间距 \(n^{-1/3}\approx0.19\,\mathrm{mm}\),远小于波长。(b) 低频:\(R_b=30\,\mu\mathrm{m}\) 泡的 Minnaert 频率约 80 kHz(按本题 \(\kappa=1\)、气体刚度 \(p-p_v\) 估算),远高于 kHz 量级的液柱频率,气泡准静态胀缩。两条件都满足。
9.4 例题 4:枪虾的一击
海水 \(\rho=1025\,\mathrm{kg/m^3}\),浅水处 \(p_0=1.013\times10^{5}\,\mathrm{Pa}\),\(p_v=2.3\,\mathrm{kPa}\),\(\sigma=0.0728\,\mathrm{N/m}\)。枪虾螯射出的水流速度 \(u=25\,\mathrm{m/s}\)。
有人这样估计射流中的压强:对“远处静止的海水”和“射流中的水”列伯努利方程,得 \(p_0-\frac12\rho u^2\)。(a) 指出这一做法错在何处,并说明自由射流内部的静压约为多少。(b) 射流前端卷起涡环。把涡环核心附近的流动近似为平面上的 Rankine 涡:半径 \(a\) 以内像刚体一样以角速度 \(\Omega\) 转动(\(v_\theta=\Omega r\)),以外为 \(v_\theta=\Omega a^2/r\),远处压强 \(p_0\)。由流体元的径向力平衡 \(\dd p/\dd r=\rho v_\theta^2/r\) 求压强分布,证明核心压强为 \(p_0-\rho v_{\max}^2\),\(v_{\max}=\Omega a\)。(c) 设 \(v_{\max}=\lambda u\),求涡核发生空化(核心压强低于 \(p_v\))所需的 \(\lambda\),用空化数 \(\sigma_{\mathrm c}=(p_0-p_v)/(\frac12\rho u^2)\) 表示,并对 \(\lambda=1/\sqrt2\) 用 Blake 阈值估计能够失稳长大的泡核半径。
题设:空化泡最大等效半径 \(R_{\max}=1.1\,\mathrm{mm}\)(按文献中 5 倍放大仿生装置的 5.6 mm 几何相似缩小所得的估计值)。求泡在最大半径时相对于远场的“势能” \(E=\frac43\pi R_{\max}^3(p_0-p_v)\) 与 Rayleigh 塌缩时间。
水听器在距泡 \(r=4.0\,\mathrm{cm}\) 处测得声压约 80 kPa,本题设它为正峰值 \(p'_{\mathrm{pk}}\)(文献未说明是峰值还是峰峰值,见 (4))。设塌缩从 \(R_{\max}\) 静止开始,泡内有分压 \(p_{gM}\) 的“缓冲气体”并绝热压缩(\(\gamma=1.4\)),峰值出现在最小半径时刻。(a) 证明 \(p'_{\mathrm{pk}}=\frac{\rho R_{\min}^2\ddot R}{r}\big|_{R_{\min}}=\frac{R_{\min}}{r}\bigl(p_g(R_{\min})-\Delta p\bigr)\),\(\Delta p=p_0-p_v\);(b) 在 \(p_{gM}\ll(\gamma-1)\Delta p\) 的近似下,导出 \(p'_{\mathrm{pk}}\approx(\gamma-1)\Delta p\,R_{\max}/(x^2r)\),\(x=R_{\min}/R_{\max}\),由此估计 \(x\)、\(R_{\min}\) 与最高温度;(c) 不作近似,给出精确结果(可数值求解)。
与“闪光要求至少 5000 K”比较,讨论估算中的主要不确定性。
(1a) 定常伯努利方程只沿同一条流线成立。周围静止海水与射流中的水不在同一条流线上:射流的水来自被塞子挤压的凹槽,其总压高于 \(p_0\)。射流离开喷口后是自由射流,流线近似平行,横向没有加速度,横向压强梯度为零,所以射流内部静压与周围相同,约为 \(p_0\)——并不会降到 \(p_0-\frac12\rho u^2\)。
(1b) 外区:\(\dv{p}{r}=\rho\frac{\Omega^2a^4}{r^3}\),从 \(r\) 积到 \(\infty\):\(p(r)=p_0-\frac{\rho\Omega^2a^4}{2r^2}\),在 \(r=a\) 处 \(p(a)=p_0-\frac12\rho v_{\max}^2\)。核内:\(\dv{p}{r}=\rho\Omega^2r\),\(p(r)=p(a)-\frac12\rho\Omega^2(a^2-r^2)\),核心 \(p_{\mathrm{core}}=p(a)-\frac12\rho v_{\max}^2=p_0-\rho v_{\max}^2\)。压强亏损一半来自外区、一半来自核内。(外区是势流,可以验证伯努利方程沿径向成立:\(p+\frac12\rho v_\theta^2=p_0\);核内不是势流——流体刚性转动,有“旋”——伯努利方程在不同流线间不成立,必须用力平衡。)
(1c) \(p_0-\rho\lambda^2u^2<p_v\iff\lambda^2>\frac{p_0-p_v}{\rho u^2}=\frac{\sigma_{\mathrm c}}{2}\)。\(\frac12\rho u^2=3.20\times10^{5}\,\mathrm{Pa}\),\(\sigma_{\mathrm c}=0.31\),\(\lambda>0.39\) 即可开始空化——涡核的最大旋转速度只需达到射流速度的四成(但要长成毫米泡需 \(\lambda\approx0.7\),见第 8.3 节)。\(\lambda=1/\sqrt2\) 时 \(p_{\mathrm{core}}=p_0-\frac12\rho u^2=-2.2\times10^{5}\,\mathrm{Pa}\)(液体处于拉伸状态),与数值模拟的核心压强亏损(平均约 \(2\times10^{5}\,\mathrm{Pa}\))同量级。由 Blake 阈值的小核近似 \(R_n\ge\frac{4\sigma}{3\sqrt3(p_v-p_{\mathrm{core}})}=0.25\,\mu\mathrm{m}\)(精确公式给出 \(0.235\,\mu\mathrm{m}\),差约 8%)。海水中这样大小的核随处可见,空化必然发生。
(2) \(E=\frac43\pi(1.1\,\mathrm{mm})^3\times9.90\times10^{4}\,\mathrm{Pa}=0.55\,\mathrm{mJ}\)。这是塌缩中外压能做的总功,它最终变成液体动能、气体内能、声能与热。Rayleigh 时间 \(t_c=0.915R_{\max}\sqrt{\rho/\Delta p}=102\,\mu\mathrm{s}\)。
(3a) 由式 (3.16),\(p'=\rho(R^2\ddot R+2R\dot R^2)/r\)。在 \(R_{\min}\) 处 \(\dot R=0\),RP 方程给出 \(\rho R\ddot R=p_g(R_{\min})-\Delta p\),于是 \(p'_{\mathrm{pk}}=\rho R_{\min}^2\ddot R/r=R_{\min}\bigl(p_g(R_{\min})-\Delta p\bigr)/r\)。
(3b) 能量守恒式 (5.9) 写成 \(p_{gM}x^{-3(\gamma-1)}=p_{gM}+(\gamma-1)\Delta p(1-x^3)\approx(\gamma-1)\Delta p\),故 \(p_g(R_{\min})=p_{gM}x^{-3\gamma}\approx(\gamma-1)\Delta p\,x^{-3}\gg\Delta p\),代入 (3a): \[p'_{\mathrm{pk}}\approx\frac{xR_{\max}}{r}(\gamma-1)\Delta p\,x^{-3}=\frac{(\gamma-1)\Delta p\,R_{\max}}{x^2r} \ \Longrightarrow\ x=\sqrt{\frac{(\gamma-1)\Delta p\,R_{\max}}{p'_{\mathrm{pk}}r}}=0.117 .\] \(R_{\min}\approx128\,\mu\mathrm{m}\),\(T_{\max}=T_0x^{-3(\gamma-1)}\approx3.9\times10^{3}\,\mathrm{K}\),所需 \(p_{gM}\approx(\gamma-1)\Delta p\,x^{3(\gamma-1)}\approx3.0\,\mathrm{kPa}\)。
(3c) 精确地,\(p_g(R_{\min})=x^{-3}\bigl[p_{gM}+(\gamma-1)\Delta p(1-x^3)\bigr]\),\(p_{gM}=\frac{(\gamma-1)\Delta p(1-x^3)}{x^{-3(\gamma-1)}-1}\),代入 (3a) 得关于 \(x\) 的单变量方程,数值解 \(x=0.121\),\(R_{\min}=133\,\mu\mathrm{m}\),\(p_{gM}=3.4\,\mathrm{kPa}\),\(T_{\max}=3.7\times10^{3}\,\mathrm{K}\)。直接数值积分 RP 方程,取 \(p_{gM}=3.4\,\mathrm{kPa}\),得 4 cm 处峰值 80.0 kPa,验证了 (3a) 中“峰值出现在最小半径”的假设(verify_shrimp.py);而 (3b) 的近似式会使峰值高估约 20%。
(4) 估算温度 \(\sim4\times10^{3}\,\mathrm{K}\) 与“至少 5000 K”同一量级,说明“绝热压缩少量缓冲气体”的图像大体自洽。不确定性:(i) \(R_{\max}\) 是缩放估计,\(x\propto R_{\max}^{1/2}\)、\(T\propto x^{-1.2}\);(ii) 水听器的有限带宽会削平极窄(\(\sim\mu\mathrm{s}\))的脉冲,真实峰值更高,对应更小的 \(x\) 与更高的温度,故这里的温度是偏低的估计;(iii) 泡并非球形(在涡环中形成、塌缩不对称),且末段速度接近声速;(iv) 缓冲气体的来源(溶解空气、来不及凝结的蒸气)与含量本身不确定;(v) 若 80 kPa 实为峰峰值:本模型中声脉冲的负瓣只有约 \(-2.6\,\mathrm{kPa}\),正峰约为峰峰值的 97%,即 \(78\,\mathrm{kPa}\),重解得 \(x=0.123\)、\(T_{\max}=3.6\times10^{3}\,\mathrm{K}\),结论对这一口径不敏感。
10 实验与数据处理:从气泡的“叮”声测多方指数
在 20 °C 的大水池中(池壁的影响可忽略),把内径不同的 8 根针头置于水面下 \(h=1.00\,\mathrm{m}\) 处,缓慢通气,每次只释放一个气泡。水听器记录气泡脱离针头时发出的衰减正弦声信号,对其做频谱分析得峰值频率 \(f\),对包络做指数拟合得振幅衰减时间 \(\tau_d\);同时用高速相机拍照,由气泡体积换算等效球半径 \(R\)。每根针头测 5 个气泡取平均,结果见表 2。表中数据为按本讲理论模型加随机误差生成的模拟数据(生成代码与随机种子见 exp_data.py)。
已知:水 \(\rho=998\,\mathrm{kg/m^3}\),\(c=1482\,\mathrm{m/s}\),\(\mu=1.00\times10^{-3}\,\mathrm{Pa\,s}\);大气压 \(p_{\mathrm{atm}}=1.013\times10^{5}\,\mathrm{Pa}\),\(g=9.81\,\mathrm{m/s^2}\);空气 \(\gamma=1.40\),泡内空气的热扩散率 \(D=2.12\times10^{-5}\,\mathrm{m^2/s}\)。照片测半径的随机误差约 0.5%(每个平均值),频率的随机误差约 0.3%;此外相机放大率的标定有约 1% 的系统不确定度(对所有半径同比例偏大或偏小)。
写出气泡处的静压 \(p\)。由 Minnaert 理论说明,若 \(\kappa\) 为常数且表面张力可忽略,\(f\) 与 \(1/R\) 成正比;估算表面张力修正,说明可以忽略。
作 \(f\)–\(1/R\) 图。(a) 先用不加权的过原点最小二乘法求斜率 \(A\),分别由残差、由已知测量误差传递估计 \(A\) 的统计不确定度,并解释两者为何不同;(b) 再按各点误差加权拟合,求多方指数 \(\kappa\) 及其统计不确定度,计入 1% 标定系统误差给出总不确定度。用带截距的直线拟合,检验截距是否与零相容。
逐点计算 \(\kappa_i=\rho(2\pi f_iR_i)^2/(3p)\),检查它们的离散是否与测量误差相符。将 \(\kappa\) 与 \(\gamma=1.40\) 和等温值 1 比较:仅凭频率数据,能否断言气泡振荡“不是绝热的”?
由 \(\tau_d\) 求各气泡的品质因数 \(Q_i=\pi f_i\tau_{d,i}\)(由第 6.4 节的定义 \(\tau_d=Q/(\pi f_0)\))与 \(\delta_i=1/Q_i\)。扣除按式 (6.5)、(6.4) 计算的声辐射阻尼与粘性阻尼(说明在什么条件下可以对这些泡使用无界液体中的公式),得到热阻尼 \(\delta_{\mathrm{th},i}\)。理论预言大泡时 \(\delta_{\mathrm{th}}\approx3(\gamma-1)\sqrt{D/(2\omega R^2)}\)。说明为何应有 \(\delta_{\mathrm{th}}\propto R^{-1/2}\),用过原点拟合检验之,由斜率反推 \(D\),并说明这一结果受半径标定误差的影响有多大。
讨论主要的系统误差来源与改进方案。
| 编号 | \(R\) / mm | \(f\) / Hz | \(\tau_d\) / ms |
|---|---|---|---|
| 1 | 0.622 | 5458 | 1.35 |
| 2 | 0.850 | 3991 | 1.94 |
| 3 | 1.088 | 3095 | 2.58 |
| 4 | 1.382 | 2475 | 3.46 |
| 5 | 1.706 | 1998 | 4.67 |
| 6 | 2.056 | 1667 | 6.11 |
| 7 | 2.467 | 1380 | 6.93 |
| 8 | 3.022 | 1131 | 9.22 |
对数据 \((x_i,y_i)\),\(i=1..n\):
过原点 \(y=Ax\):使 \(\sum(y_i-Ax_i)^2\) 最小,\(A=\frac{\sum x_iy_i}{\sum x_i^2}\);残差 \(e_i=y_i-Ax_i\),\(s_A=\sqrt{\frac{\sum e_i^2}{(n-1)\sum x_i^2}}\)。
带截距 \(y=a+bx\):记 \(\Delta=n\sum x_i^2-(\sum x_i)^2\),\(b=\frac{n\sum x_iy_i-\sum x_i\sum y_i}{\Delta}\),\(a=\frac{\sum x_i^2\sum y_i-\sum x_i\sum x_iy_i}{\Delta}\);\(s^2=\frac{\sum(y_i-a-bx_i)^2}{n-2}\),\(s_b=\sqrt{ns^2/\Delta}\),\(s_a=\sqrt{s^2\sum x_i^2/\Delta}\)。
已知每点误差时(过原点):若 \(y_i\) 的标准误差为 \(\sigma_i\),不加权拟合的斜率误差为 \(s_A=\sqrt{\sum x_i^2\sigma_i^2}\big/\sum x_i^2\);更好的做法是按 \(w_i=1/\sigma_i^2\) 加权,使 \(\sum w_i(y_i-Ax_i)^2\) 最小:\(A=\frac{\sum w_ix_iy_i}{\sum w_ix_i^2}\),\(s_A=\bigl(\sum w_ix_i^2\bigr)^{-1/2}\)。\(x_i\) 的相对误差 \(e_x\) 可折算到 \(y\) 上:\(\sigma_i=Ax_i\sqrt{e_x^2+e_y^2}\)。拟合优度 \(\chi^2=\sum(e_i/\sigma_i)^2\) 应与自由度 \(n-1\) 相近。
不确定度传递:若 \(\kappa\propto A^2\),则 \(s_\kappa/\kappa=2s_A/A\);独立的统计与系统误差按平方和合成。
(分母中的 \(n-1\)、\(n-2\) 是拟合参数个数为 1、2 时的自由度。)
(1) \(p=p_{\mathrm{atm}}+\rho gh=1.111\times10^{5}\,\mathrm{Pa}\)。由式 (6.2),\(f=\frac{1}{2\pi R}\sqrt{\frac{3\kappa p}{\rho}}\equiv\frac AR\),\(A=\frac{1}{2\pi}\sqrt{3\kappa p/\rho}\)。表面张力使 \(\omega^2\) 相对改变 \(\frac{3\kappa-1}{3\kappa}\frac{2\sigma}{Rp}\),对最小的泡(\(R=0.62\,\mathrm{mm}\))为 \(0.16\%\),频率改变 \(0.08\%\),小于频率测量误差,可以忽略。
(2) \(x_i=1/R_i\)。(a) 不加权过原点拟合得 \(A=3.3954\,\mathrm{m/s}\);由残差 \(s_A=0.0056\,\mathrm{m/s}\),由已知误差 \(\sigma_i=Ax_i\sqrt{0.005^2+0.003^2}\) 传递得 \(s_A=0.0100\,\mathrm{m/s}\),两者相差 1.8 倍,而 \(\chi^2=8.2\) 与自由度 7 相符,说明题给误差合理。差别主要来自误差结构:本题各点的相对误差相同,\(\sigma_i\propto x_i\),小泡(大 \(x\))的绝对误差大;残差公式却默认各点误差相等。仅这一因素就使两种估计之比为 \(\sqrt{n\sum x_i^4/(\sum x_i^2)^2}=1.43\),余下的 \(1.78/1.43\approx1.24\) 才是 8 个点的样本涨落。同时,不加权拟合让误差最大的小泡主导了斜率。 (b) 正确做法是按 \(w_i=1/\sigma_i^2\) 加权:\(A=3.4042\,\mathrm{m/s}\),\(s_A=0.0070\,\mathrm{m/s}\),\(\chi^2=6.6\), \[\kappa=\frac{\rho(2\pi A)^2}{3p}=1.370,\qquad s_\kappa^{\mathrm{stat}}=2\kappa\frac{s_A}{A}=0.006,\qquad s_\kappa^{\mathrm{sys}}=2\kappa\times1\%=0.027,\] 合成 \(s_\kappa=\sqrt{0.006^2+0.027^2}=0.028\),即 \(\kappa=1.37\pm0.03\)。(对相对误差恒定的数据,加权拟合等价于直接平均 \(\kappa_i\),见 (3)。)带截距(不加权)拟合:\(b=(3.382\pm0.012)\,\mathrm{m/s}\),\(a=(13\pm10)\,\mathrm{Hz}\),截距在约 \(1.3\sigma\) 内与零相容,支持正比关系(图 9 左)。
(3) \(\kappa_i=1.363,\ 1.360,\ 1.341,\ 1.383,\ 1.374,\ 1.389,\ 1.370,\ 1.381\),平均 \(1.370\pm0.006\)(与 (2b) 的加权结果一致),标准差 \(0.015\)(1.1%);按半径 0.5%、频率 0.3% 的误差,\(\kappa_i\propto(fR)^2\) 的预期离散为 \(2\sqrt{0.5^2+0.3^2}\,\%=1.2\%\),即 0.016,两者相符。 与 1.40 比较:只计统计误差,\(1.40-1.370=0.030\) 约为 \(5.3\sigma\),似乎显著;但若相机标定有 1% 的系统误差,总误差 0.028,偏离只有 \(1.1\sigma\)。仅凭频率数据不能断言非绝热——\(\kappa\propto R^2\) 对半径标定极其敏感,必须先做标定(见 (5))。\(\kappa\) 远大于 1 则是确定的:振荡至少“接近”绝热。物理上,泡壁处厚约 \(\ell_{\mathrm{th}}=\sqrt{D/\omega}\approx\) 25 μm–55 μm 的气体层与水换热、近似等温,这一层占泡体积的比例约 \(3\ell_{\mathrm{th}}/R\sim\) 5%~12%,应把有效多方指数拉低到 1.35~1.38(附录 C),与结果相容。下一问给出不依赖半径标定的检验。
(4) \(Q_i=23.1,\ 24.3,\ 25.1,\ 26.9,\ 29.3,\ 32.0,\ 30.0,\ 32.8\)。声辐射 \(\delta_{\mathrm{rad},i}=2\pi f_iR_i/c\approx0.0144\)(几乎与 \(R\) 无关,因为 \(fR\) 近似为常数),粘性 \(\delta_{\mathrm{vis},i}=4\mu/(2\pi\rho f_iR_i^2)\le3\times10^{-4}\)。这两个公式是无界液体中的结果。粘性项只涉及泡附近的流动,只要泡离针头、壁面远大于 \(R\) 即可。辐射项还涉及远场:水面的反相镜像(例题 2)使辐射功率乘以因子 \(\frac{1}{4\pi}\int_{\text{下半空间}}\bigl(2-2\cos(2kd\cos\theta)\bigr)\dd\Omega=1-\frac{\sin2kd}{2kd}\)(例题 2(4b) 的积分对任意 \(kd\) 做出,\(\abs{1-\ee^{-\ii\alpha}}^2=2-2\cos\alpha\)),这里镜像距离取 \(d=h\)(针头深度)。本题 \(h={1\,}\mathrm{m}\)、波长 0.27 m–1.3 m,\(kd=5\)~23,该因子在 0.93~1.07 之间,\(\delta_{\mathrm{rad}}\) 有约 \(\pm0.001\) 的系统偏差;若水深只有 10 cm,低频泡的辐射将被削弱近一个数量级,无界公式完全失效。扣除后 \(\delta_{\mathrm{th},i}=0.0285,\ 0.0265,\ 0.0254,\ 0.0225,\ 0.0196,\ 0.0166,\ 0.0188,\ 0.0160\)。 \(R^{-1/2}\) 依赖的来源:热量只在厚 \(\ell_{\mathrm{th}}\) 的边界层内交换,耗散正比于边界层体积占比 \(\sim\ell_{\mathrm{th}}/R=\sqrt{D/(\omega R^2)}\);又 \(\omega R\approx2\pi A\) 为常数,所以 \(\delta_{\mathrm{th}}\propto\sqrt{D/(2\pi A\,R)}\propto R^{-1/2}\)。具体地 \(\delta_{\mathrm{th}}\approx3(\gamma-1)\sqrt{D/(4\pi A)}\,R^{-1/2}\)。过原点拟合 \(\delta_{\mathrm{th}}=bR^{-1/2}\): \[b=(7.92\pm0.24)\times10^{-4}\,\mathrm{m^{1/2}},\qquad D=4\pi A\Bigl[\frac{b}{3(\gamma-1)}\Bigr]^2=(1.86\pm0.11)\times10^{-5}\,\mathrm{m^2/s},\] 比已知值 2.12 × 10−5 m2/s 低约 12%(图 9 右)。偏差主要来自“领头阶”近似本身(其次是上面的辐射镜像因子):本数据 \(\abs{y}=25\)~55,此时它把 \(\delta_{\mathrm{th}}\) 高估约 4%~10%(附录 C 的完整式),反推的 \(D\propto b^2\) 被相应低估约 10%~20%。 关键在于:\(\delta_{\mathrm{th},i}\) 由 \(f\) 与 \(\tau_d\) 直接得到,半径只通过很小的 \(\delta_{\mathrm{rad}}\) 与横坐标 \(R^{-1/2}\) 进入,两者的效应部分抵消:1% 的标定误差只使 \(b\) 改变约 0.1%,\(D\propto Ab^2\) 改变不到 1%,远小于其 6% 的统计误差。\(b\) 偏离零 33 个标准差,且 \(R^{-1/2}\) 依赖与热传导理论相符——这是比 \(\kappa<1.40\) 稳健得多的“非绝热”证据:(3)(4) 合起来才构成完整的结论。
(5) 系统误差与改进:
半径标定:\(\kappa\propto R^2\),相机放大率 1% 的系统误差即造成 \(\kappa\) 2%(约 0.03)的系统偏差,是统计误差的 3 倍多,是 \(\kappa\) 的主要误差源。改进:拍摄已知直径的标准球标定;或收集已知个数的气泡、测总气体体积求平均半径。
近壁与针头效应:刚脱离针头的泡离针尖很近,固体边界的同号“镜像”使有效质量增大为 \(m_{\mathrm{eff}}(1+R/2d)\)、频率降低(使 \(\kappa\) 偏小,例题 1(5));离自由液面太近则相反(例题 2)。例如 \(d=5R\) 时频率改变约 5%,\(\kappa\) 改变约 10%。改进:增大水深与针头间距,或只分析泡离开针尖若干直径后的信号。
非球形与大振幅:脱离时泡有形变,形状振荡会混入其他频率;初始振幅较大时有非线性频移。改进:只拟合信号尾部的小振幅段。
水箱混响与有限尺寸:水箱壁反射使 \(\tau_d\) 偏大(\(Q\) 偏大)。改进:使用大水箱或吸声衬里,取直达声时间窗。
温度与溶解气体:水温影响 \(\rho\)、\(c\)、\(p_v\),泡内水蒸气也改变有效的 \(\gamma\) 与 \(D\)。
11 小结:推导链
运动学:不可压缩 + 球对称 + 泡壁条件 \(\Rightarrow\) \(u=R^2\dot R/r^2\)(式 (3.4))\(\Rightarrow\) 液体动能 \(K=2\pi\rho R^3\dot R^2\),有效质量为排开液体的 3 倍(式 (3.5)、(3.6))。
动力学:泡壁力平衡(拉普拉斯压强 + 法向粘性应力,式 (4.1))与 [动能定理 或 非定常伯努利方程 (3.9)] \(\Rightarrow\) RP 方程 (4.2)。
大振幅塌缩:\(p_B=0\) 时能量积分 \(\Rightarrow\) 式 (5.1) \(\Rightarrow\) 塌缩时间 (5.3)(Beta 函数)、奇点 \(\abs{\dot R}\propto R^{-3/2}\)、压强峰 (5.7);可压缩性在 \(R\sim R_0/30\) 时变重要(式 (5.8))。泡内有气体时能量守恒 (5.9) \(\Rightarrow\) 最小半径与最高温度 (5.10),它们解释了空蚀、声致发光、虾致发光的能量集中机制。
小振幅振荡:RP 方程线性化 \(\Rightarrow\) Minnaert 频率 (6.2);有效多方指数由热扩散决定(式 (3.12)、(6.3));阻尼 = 粘性 (6.4) + 声辐射 (6.5)(来自单极子辐射 (3.16))+ 热 (6.6)。
集体效应:低频、长波长下,“液体的密度 + 气体的压缩率” \(\Rightarrow\) Wood 公式 (7.1),最小声速约 20 m/s,热巧克力效应 (7.3)。
空化的发生:沿流线的伯努利方程 \(\Rightarrow\) 空化数 (8.1)(只能用于同一条流线);涡核的径向力平衡给出压强亏损 \(\rho v_{\max}^2\)(式 (8.2));泡核的 Blake 阈值决定何时失稳长大;Rayleigh 时间、单极子声压与绝热压缩共同给出枪虾的一击的自洽图像。
边界的影响(例题):镜像法——刚性壁用同号镜像(有效质量增大、塌缩变慢 \(1+0.203R_0/d\)),自由液面用反号镜像(有效质量减小、辐射退化为偶极子,功率乘以 \(1-\frac{\sin2kd}{2kd}\));分层液柱用界面上压强与法向速度连续求本征频率。
参考文献
附录 A 球对称流中的粘性应力
A.1 从剪切到拉伸:法向粘性应力
简单剪切可以分解为“纯拉伸 + 刚性转动”
考虑普通物理中的简单剪切流 \(u_x=\dot\gamma\,y\),\(u_y=0\),牛顿粘性定律给出切应力 \(\tau_{xy}=\tau_{yx}=\mu\dot\gamma\),法向粘性应力 \(\tau_{xx}=\tau_{yy}=0\)。把坐标轴转 \(45^\circ\):\(x'=(x+y)/\sqrt2\),\(y'=(-x+y)/\sqrt2\)。速度分量变为 \[u_{x'}=\frac{u_x+u_y}{\sqrt2}=\frac{\dot\gamma}{2}(x'+y'),\qquad u_{y'}=\frac{-u_x+u_y}{\sqrt2}=-\frac{\dot\gamma}{2}(x'+y') ,\] (用了 \(y=(x'+y')/\sqrt2\))。于是 \(\pdv{u_{x'}}{x'}=+\dot\gamma/2\)(沿 \(x'\) 拉伸),\(\pdv{u_{y'}}{y'}=-\dot\gamma/2\)(沿 \(y'\) 压缩);剩下的 \(\pdv{u_{x'}}{y'}=\dot\gamma/2\)、\(\pdv{u_{y'}}{x'}=-\dot\gamma/2\) 构成角速度为 \(-\dot\gamma/2\) 的刚性转动,刚性转动不使液体变形,不产生粘性应力。
同一应力在 45° 轴上的法向分量
取一个三棱柱微元:两直角面法向分别为 \(x\)、\(y\),斜面法向 \(\uvec{n}=(1,1)/\sqrt2\)(即 \(x'\) 方向)。微元很小时惯性与体积力可忽略(它们 \(\propto\) 体积,而面力 \(\propto\) 面积),面力必须平衡,由此得斜面上单位面积受力 \(t_i=\sum_j\tau_{ij}n_j\)。代入 \(\tau_{xy}=\tau_{yx}=\mu\dot\gamma\)、\(\tau_{xx}=\tau_{yy}=0\):\(\vect t=\mu\dot\gamma(1,1)/\sqrt2\),其法向分量 \[\tau_{x'x'}=\vect t\cdot\uvec n=\mu\dot\gamma=2\mu\cdot\frac{\dot\gamma}{2},\qquad\text{同理}\quad\tau_{y'y'}=-\mu\dot\gamma=2\mu\cdot\Bigl(-\frac{\dot\gamma}{2}\Bigr).\]
在纯拉伸/压缩的主轴上,法向粘性应力等于 \(2\mu\) 乘以该方向的拉伸速率。理由是流体各向同性:它对一个“沿某轴以速率 \(e\) 拉伸、沿垂直轴以速率 \(e\) 压缩”的形变的响应,不依赖于这种形变是怎样得到的。(一般形式为 \(\tau_{ij}=\mu(\partial u_i/\partial x_j+\partial u_j/\partial x_i)\),对不可压缩流体成立。)
A.2 球对称流的拉伸速率与应力
由对称性,\(r,\theta,\varphi\) 三个方向就是形变主轴。
径向:相距 \(\dd r\) 的两层液体速度差为 \(\pdv{u}{r}\dd r\),拉伸速率 \(e_{rr}=\pdv{u}{r}=-2R^2\dot R/r^3\)。
切向:半径 \(r\) 处一个随液体运动的小圆环,周长 \(2\pi r\) 的相对增长率 \(\dot r/r=u/r\),故 \(e_{\theta\theta}=e_{\varphi\varphi}=u/r=R^2\dot R/r^3\)。
检验:\(e_{rr}+e_{\theta\theta}+e_{\varphi\varphi}=0\),即体积不变(不可压缩)。
这是三个方向同时形变,而第 A.1 节只讨论了“一轴拉伸、一轴压缩”的平面情形。记 \(e\equiv R^2\dot R/r^3\),主轴形变率为 \((e_{rr},e_{\theta\theta},e_{\varphi\varphi})=(-2e,e,e)\),可以拆成两个平面纯形变之和:\((-e,e,0)+(-e,0,e)\),每一个都是第 A.1 节的情形。粘性应力与形变率成正比(牛顿流体的线性关系),所以两部分的应力可以相加: \[\tau_{rr}=2\mu(-e)+2\mu(-e)=2\mu e_{rr},\qquad\tau_{\theta\theta}=2\mu e+0=2\mu e_{\theta\theta},\qquad\tau_{\varphi\varphi}=2\mu e_{\varphi\varphi}.\] 于是 \(\tau_{rr}=-4\mu R^2\dot R/r^3\),\(\tau_{\theta\theta}=\tau_{\varphi\varphi}=2\mu R^2\dot R/r^3\)。
A.3 体内粘性合力为零
取一个以径线为轴、张角为 \(\alpha\times\alpha\)(\(\alpha\to0\))、位于 \(r\sim r+\dd r\) 的小“锥台”微元,体积 \(r^2\alpha^2\dd r\)。
两个径向端面(面积 \(r^2\alpha^2\))上的径向力之差:\(\dv{}{r}\bigl(\tau_{rr}r^2\bigr)\alpha^2\dd r\)。
四个侧面(每个面积 \(r\alpha\,\dd r\))上的法向应力 \(\tau_{\theta\theta}\)(或 \(\tau_{\varphi\varphi}\))沿侧面的外法向,而侧面外法向与微元中心径向的夹角为 \(90^\circ+\alpha/2\),故每个侧面贡献径向力 \(-\tau_{\theta\theta}\,r\alpha\,\dd r\sin(\alpha/2)\approx-\frac12\tau_{\theta\theta}r\alpha^2\dd r\)。四个面合计 \(-(\tau_{\theta\theta}+\tau_{\varphi\varphi})r\alpha^2\dd r\)。这和气球蒙皮的“环向张力产生向心合力”是同一个几何效应。
除以体积,单位体积的径向粘性力 \[f_r=\frac{1}{r^2}\dv{}{r}\bigl(r^2\tau_{rr}\bigr)-\frac{\tau_{\theta\theta}+\tau_{\varphi\varphi}}{r} =\frac{1}{r^2}\dv{}{r}\Bigl(-\frac{4\mu R^2\dot R}{r}\Bigr)-\frac{4\mu R^2\dot R}{r^4}=\frac{4\mu R^2\dot R}{r^4}-\frac{4\mu R^2\dot R}{r^4}=0 .\]
若把各向同性的压强写成 \(\tau_{rr}=\tau_{\theta\theta}=\tau_{\varphi\varphi}=-p\),上式给出 \(f_r=\frac{1}{r^2}\dv{}{r}(-r^2p)+\frac{2p}{r}=-\dv{p}{r}\),与式 (3.1) 一致,说明侧面项的处理正确。
所以粘性不改变液体内部的运动方程(欧拉方程 (3.3) 与伯努利方程 (3.9) 仍精确成立),只通过泡壁的力平衡 (4.1) 进入 RP 方程。
A.4 为什么“幂次计数”会误导
一个看似有力的反驳
速度 \(u\propto r^{-2}\),求一次导数后径向应力 \(\tau_{rr}\propto r^{-3}\),再乘以面积 \(\propto r^{2}\),球面上的径向粘性力 \(\propto r^{-1}\)——随 \(r\) 变化,内外两层球面上的力似乎不能抵消,体内粘性力似乎不为零。这个直觉对了一半,下面把账算全。
第 1 步:直觉算对的那一半
取第 A.3 节的锥台微元,张开的立体角为 \(\Omega=\alpha^2\)。半径 \(r\) 处端面(面积 \(r^2\Omega\))所受的径向粘性力为 \[F_{\mathrm{end}}(r)=\tau_{rr}\,r^2\Omega=-\frac{4\mu R^2\dot R}{r^3}\,r^2\Omega=-\frac{4\mu R^2\dot R\,\Omega}{r}\ \propto r^{-1},\] 确实随 \(r\) 变化,与直觉一致。外端面(\(r+\dd r\))受外侧液体沿 \(+r\) 方向的力 \(+F_{\mathrm{end}}(r+\dd r)\),内端面受内侧液体沿 \(-r\) 方向的力,合起来是 \[\tag{A.1} \dd F_{\mathrm{end}}=F_{\mathrm{end}}(r+\dd r)-F_{\mathrm{end}}(r)=\dv{F_{\mathrm{end}}}{r}\dd r=+\frac{4\mu R^2\dot R}{r^2}\,\Omega\,\dd r\neq0 .\] 若只有端面,微元确实受到净的径向粘性力。
第 2 步:直觉漏掉的那一半——环向应力
微元还有四个侧面。泡膨胀(\(\dot R>0\))时,液体向外流,每个液体环的周长 \(2\pi r\) 随之增大,沿切向被拉伸(\(e_{\theta\theta}=e_{\varphi\varphi}=u/r>0\)),所以侧面上的环向粘性应力是拉应力:\(\tau_{\theta\theta}=\tau_{\varphi\varphi}=2\mu R^2\dot R/r^3>0\)。如果四个侧面彼此平行,这些拉力会两两抵消;但锥台的侧面沿径向张开,相对的两个侧面夹角为 \(\alpha\),每个侧面上的拉力都向内偏转 \(\alpha/2\),合起来指向球心——这正是充气气球的蒙皮张力把气球“往里收”的几何效应。由第 A.3 节第 2 条,侧面的径向合力为 \[\tag{A.2} \dd F_{\mathrm{hoop}}=-(\tau_{\theta\theta}+\tau_{\varphi\varphi})\,r\,\Omega\,\dd r=-\frac{4\mu R^2\dot R}{r^2}\,\Omega\,\dd r .\]
第 3 步:幂次相同,系数相反
比较式 (A.1) 与 (A.2):端面差来自 \(\dv{}{r}(r^2\cdot r^{-3})\sim r^{-2}\),侧面项来自 \(r\cdot r^{-3}\sim r^{-2}\),两项是同一幂次。所以幂次计数只能告诉我们“两项同阶”,判断不了是否抵消;抵消与否完全取决于系数,而系数恰好是 \(+4\) 与 \(-4\)。直觉的毛病在于只数了端面这一项。
第 4 步:抵消不是巧合——只有 \(n=2\) 才抵消
为看清系数从哪里来,把速度形式地换成一般幂律 \(u=A/r^n\)(\(A\) 为常量),沿用同一个应力公式 \(\tau=2\mu e\)1: \[e_{rr}=\dv{u}{r}=-\frac{nA}{r^{n+1}},\qquad e_{\theta\theta}=e_{\varphi\varphi}=\frac{u}{r}=\frac{A}{r^{n+1}} .\] 代入第 A.3 节的 \(f_r\) 公式: \[f_r=\frac{1}{r^2}\dv{}{r}\Bigl(r^2\cdot\frac{-2\mu nA}{r^{n+1}}\Bigr)-\frac{2}{r}\cdot\frac{2\mu A}{r^{n+1}} =\frac{2\mu A\,n(n-1)}{r^{n+2}}-\frac{4\mu A}{r^{n+2}} =\frac{2\mu A\,(n-2)(n+1)}{r^{n+2}} .\] 它只在 \(n=2\)(以及无物理意义的 \(n=-1\))时为零。而 \(n=2\) 恰好是不可压缩条件的要求:流过半径 \(r\) 球面的体积流量 \(4\pi r^2u\) 必须与 \(r\) 无关(式 (3.4)),即 \(e_{rr}+e_{\theta\theta}+e_{\varphi\varphi}=0\)。换句话说,“拉伸速率之和为零”使得径向的压缩应力与环向的拉伸应力在体内恰好平衡。
第 5 步:与几何无关的一般证明
上面的抵消其实对任何不可压缩势流都成立,不限于球对称。用直角坐标,粘性应力为 \(\tau_{ij}=\mu\bigl(\pdv{u_i}{x_j}+\pdv{u_j}{x_i}\bigr)\)(附录 A.1 的要点框),作用在单位体积液体上的粘性力的第 \(i\) 个分量是各面应力差之和(与第 3 节由压强差求合力的做法相同,只是现在每个面上还有切向分量): \[f_i=\sum_{j}\pdv{\tau_{ij}}{x_j}=\mu\sum_j\frac{\partial^2u_i}{\partial x_j^2}+\mu\,\pdv{}{x_i}\Bigl(\sum_j\pdv{u_j}{x_j}\Bigr) =\mu\,\nabla^2u_i+\mu\,\pdv{}{x_i}(\nabla\cdot\vect u),\] 其中 \(\nabla^2\equiv\partial_x^2+\partial_y^2+\partial_z^2\),第二步交换了求导次序。对不可压缩流第二项为零。对势流 \(u_i=\pdv{\phi}{x_i}\),再交换一次求导次序: \[\nabla^2u_i=\pdv{}{x_i}\bigl(\nabla^2\phi\bigr)=\pdv{}{x_i}(\nabla\cdot\vect u)=0 .\] 所以 \(f_i\equiv0\)。本讲的 \(\phi=-R^2\dot R/r\)(式 (3.8))满足 \(\nabla^2\phi=0\),是这一结论的特例。
第 6 步:合力为零,粘性却仍在起作用
体内粘性合力为零,不等于粘性没有效果:
泡壁处应力不平衡。 体内每个微元两侧都有液体,端面差被环向拉力抵消;而泡壁外侧是液体、内侧是气体,壁面上的 \(\tau_{rr}(R)=-4\mu\dot R/R\) 没有东西与之抵消,于是出现在力平衡 (4.1) 中,成为 RP 方程的 \(4\mu\dot R/(\rho R)\) 项。
耗散不为零。 力对一个微元做功有两部分:合力乘微元的平移速度,以及应力乘微元的形变速率。合力为零只消去了前一部分;后一部分 \(\sum_i\tau_{ii}e_{ii}=2\mu\sum_ie_{ii}^2>0\) 仍把机械能变成热(下一节)。能量从泡壁输入:壁面粘性应力对液体的功率为 \(-\tau_{rr}(R)\cdot4\pi R^2\dot R=16\pi\mu R\dot R^2\),它在体内一路传递而不改变任何微元的动能,最终全部变为热,恰好等于下一节的耗散积分 \(\mathcal D\)。
一般幂律的 \(f_r=2\mu A(n-2)(n+1)/r^{n+2}\)、\(n=2\) 时端面项与环向项为 \(\pm4\mu A/r^4\)、直角坐标中 \(\mu(\nabla^2u_i+\partial_i\nabla\!\cdot\!\vect u)=0\)(对 \(\phi=-A/r\)),以及壁面粘性功率 \(=\mathcal D=16\pi\mu R\dot R^2\),均由 scripts/verify_visc.py 核对通过。
A.5 能量核对:粘性耗散
一根截面积为 \(S\)、受拉力 \(F\)、长度以 \(\dot L\) 增长的细杆,外力功率 \(F\dot L=(F/S)(\dot L/L)(SL)\),即“应力 \(\times\) 拉伸速率 \(\times\) 体积”。同理,粘性应力对单位体积液体的形变功率为 \(\sum\tau_{ii}e_{ii}=2\mu(e_{rr}^2+e_{\theta\theta}^2+e_{\varphi\varphi}^2)\),这部分机械能变为热。对全部液体积分: \[\mathcal D=\int_R^\infty2\mu\,\frac{(R^2\dot R)^2}{r^6}\,(4+1+1)\,4\pi r^2\dd r=48\pi\mu R^4\dot R^2\int_R^\infty\frac{\dd r}{r^4}=16\pi\mu R\dot R^2 .\] 泡壁对液体的推力(单位面积)为 \(p_L-\tau_{rr}=p_B-2\sigma/R\)(式 (4.1)),功率 \((p_B-2\sigma/R)4\pi R^2\dot R\);远场功率 \(-p_\infty4\pi R^2\dot R\)。机械能定理 \(\dv{K}{t}=\) 输入功率 \(-\mathcal D\): \[4\pi\rho R^2\dot R\Bigl(R\ddot R+\frac32\dot R^2\Bigr)=\Bigl(p_B-\frac{2\sigma}{R}-p_\infty\Bigr)4\pi R^2\dot R-16\pi\mu R\dot R^2 ,\] 两边除以 \(4\pi\rho R^2\dot R\) 即得 RP 方程 (4.2),且 \(16\pi\mu R\dot R^2/(4\pi R^2\dot R)=4\mu\dot R/R\)。两种算法一致。\(\mathcal D\) 的积分与体内合力为零均已用 sympy 核对(verify_core.py)。
附录 B Beta 函数与 Rayleigh 塌缩时间
\(\Gamma(a)=\int_0^\infty t^{a-1}\ee^{-t}\dd t\)(\(a>0\)),分部积分可得 \(\Gamma(a+1)=a\Gamma(a)\),\(\Gamma(n+1)=n!\),且 \(\Gamma(1/2)=\sqrt\pi\)。\(\mathrm{B}(a,b)=\int_0^1s^{a-1}(1-s)^{b-1}\dd s\)(\(a,b>0\))。
写出乘积为二重积分 \(\Gamma(a)\Gamma(b)=\int_0^\infty\!\!\int_0^\infty t^{a-1}w^{b-1}\ee^{-(t+w)}\dd t\,\dd w\)。换元 \(t=vs\),\(w=v(1-s)\)(\(v\in(0,\infty)\),\(s\in(0,1)\))。雅可比行列式 \[\abs{\frac{\partial(t,w)}{\partial(v,s)}}=\abs{\begin{matrix}s & v\\ 1-s & -v\end{matrix}}=\abs{-vs-v(1-s)}=v ,\] 于是 \(\Gamma(a)\Gamma(b)=\int_0^\infty v^{a+b-1}\ee^{-v}\dd v\int_0^1s^{a-1}(1-s)^{b-1}\dd s=\Gamma(a+b)\,\mathrm B(a,b)\)。
正文中 \(I=\int_0^1x^{3/2}(1-x^3)^{-1/2}\dd x\)。令 \(s=x^3\),\(x=s^{1/3}\),\(\dd x=\frac13s^{-2/3}\dd s\),\(x^{3/2}=s^{1/2}\): \[I=\frac13\int_0^1s^{1/2-2/3}(1-s)^{-1/2}\dd s=\frac13\mathrm B\Bigl(\frac56,\frac12\Bigr)=\frac{\Gamma(5/6)\sqrt\pi}{3\,\Gamma(4/3)} .\] 数值 \(\Gamma(5/6)=1.128787\),\(\Gamma(4/3)=0.892980\),\(I=0.746834\)。\(t_c=\sqrt{3/2}\,I\,\tau=0.914681\,\tau\);合并常数 \(\sqrt{3/2}\cdot\sqrt\pi/3=\sqrt{\pi/6}\) 即得式 (5.3)。
附录 C 泡内热传导与复多方指数
(1) 泡内压强空间均匀(泡远小于气体中的声波长),只随时间变化;(2) 泡壁温度固定为液体温度 \(T_0\)(液体热容量大、导热好);(3) 小振幅,平衡时泡内气压为 \(p_{g0}\)(含拉普拉斯压强,与正文第 6.3 节相同),所有扰动 \(\propto\ee^{\ii\omega t}\),只保留一阶;(4) 气体为理想气体,热导率 \(k_{\mathrm{th}}\) 为常量。
热力学第一定律:对一小份气体(物质的量 \(n\)、体积 \(V\)),从普通物理的形式 \(\delta Q=nC_V\,\dd T+p\,\dd V\) 出发。由 \(pV=nR_{\mathrm m}T\) 取微分,\(p\,\dd V+V\dd p=nR_{\mathrm m}\dd T\),用它消去 \(p\,\dd V\),并用 \(C_p=C_V+R_{\mathrm m}\): \[\delta Q=n(C_V+R_{\mathrm m})\dd T-V\dd p=nC_p\,\dd T-V\dd p .\] 除以体积 \(V\),\(nC_p/V=\rho_gc_p\)(\(c_p\) 为单位质量的定压比热),得单位体积吸热 \(\rho_gc_p\,\dd T-\dd p\)。这种写法的好处是:泡内压强处处相同、只随时间变化,\(\dd p\) 项与位置无关。热量由傅里叶定律(热流密度 \(=-k_{\mathrm{th}}\nabla T\))传导,对球壳做与第 3.8 节相同的“流入减流出”,单位体积的净流入功率为 \(k_{\mathrm{th}}\frac{1}{r^2}\pdv{}{r}\bigl(r^2\pdv{T}{r}\bigr)\)。线性阶下气体的对流项是二阶小量,故 \[\rho_gc_p\pdv{T'}{t}-\dv{p'}{t}=k_{\mathrm{th}}\frac{1}{r^2}\pdv{}{r}\Bigl(r^2\pdv{T'}{r}\Bigr).\]
用 \(\rho_gc_pT_0=\frac{\gamma}{\gamma-1}p_{g0}\)(理想气体 \(\rho_gc_p=nC_p/V=\frac{\gamma}{\gamma-1}\frac{nR_{\mathrm m}}{V}=\frac{\gamma}{\gamma-1}\frac{p_{g0}}{T_0}\)),记 \(\Theta=T'/T_0\),\(P=p'/p_{g0}\),\(D=k_{\mathrm{th}}/(\rho_gc_p)\),代入 \(\ee^{\ii\omega t}\): \[\ii\omega\Theta-\ii\omega\frac{\gamma-1}{\gamma}P=D\,\frac1r\dv{^2(r\Theta)}{r^2}.\]
特解为常数 \(\Theta_p=\frac{\gamma-1}{\gamma}P\)(这正是绝热时的温度扰动)。齐次方程 \(\frac1r(r\Theta_h)''=q^2\Theta_h\),\(q^2=\ii\omega/D\),在球心正则的解为 \(\Theta_h\propto\sinh(qr)/r\)。由 \(\Theta(R)=0\): \[\Theta(r)=\frac{\gamma-1}{\gamma}P\left[1-\frac{R\sinh(qr)}{r\sinh(qR)}\right].\]
理想气体 \(p=\rho_gR_{\mathrm m}T/M\)(\(M\) 为摩尔质量)线性化:\(\frac{p'}{p_{g0}}=\frac{\rho'}{\rho_g}+\frac{T'}{T_0}\),即 \(\rho'/\rho_g=P-\Theta\)。泡内气体总质量不变:\(\int\rho'\dd V+\rho_g\delta V=0\),即 \(\frac{\delta V}{V}=-\bigl(P-\avg{\Theta}\bigr)\),\(\avg{\cdot}\) 表示对泡体积的平均。利用 \(\int_0^R r\sinh(qr)\,\dd r=\frac{R\cosh(qR)}{q}-\frac{\sinh(qR)}{q^2}\): \[\avg{\frac{R\sinh(qr)}{r\sinh(qR)}}=\frac{3}{R^3}\cdot\frac{R}{\sinh(qR)}\int_0^R r\sinh(qr)\dd r=\frac{3}{(qR)^2}\bigl(qR\coth(qR)-1\bigr)\equiv G(y),\quad y=qR .\] 所以 \(\avg\Theta=\frac{\gamma-1}{\gamma}P(1-G)\),\(\frac{\delta V}{V}=-\frac{P}{\gamma}\bigl[1+(\gamma-1)G\bigr]\)。按多方定律的线性形式 \(P=-\kappa_c\,\delta V/V\) 定义复多方指数,即得式 (6.3)。
小 \(\abs y\):\(y\coth y=1+y^2/3-y^4/45+\cdots\),\(G=1-y^2/15+\cdots\),\(\kappa_c\approx1+\ii\frac{\gamma-1}{15\gamma}\frac{\omega R^2}{D}\)——实部为 1(等温),热阻尼 \(\delta_{\mathrm{th}}\approx\frac{\gamma-1}{15\gamma}\mathrm{Pe}\) 随 \(R\) 减小而趋零。大 \(\abs y\):\(\coth y\to1\),\(G\approx3/y\),\(1/y=\sqrt{D/(\omega R^2)}\,\ee^{-\ii\pi/4}\), \[\kappa_c\approx\gamma-3\gamma(\gamma-1)\sqrt{\frac{D}{2\omega R^2}}\,(1-\ii),\qquad
\delta_{\mathrm{th}}\approx\frac{\Im\kappa_c}{\Re\kappa_c}\approx3(\gamma-1)\sqrt{\frac{D}{2\omega R^2}} .\] 所以热阻尼在 \(\mathrm{Pe}\sim10\)–\(100\)(空气泡 \(R_0\sim10\,\mu\mathrm{m}\))附近最大(图 6)。verify_acoustic.py 用有限差分直接求解上面的边值问题,与解析式 (6.3) 在 \(R=1\,\mu\mathrm{m}\) 到 1 mm 范围内相对差 \(<2\times10^{-3}\)。式 (6.3) 与 Prosperetti 等人的结果 [31], [39] 相同(他们写作 \(\Phi=3\kappa_c\))。
附录 D 名词解释(中英对照)
| 中文 | 英文 | 一句话定义 |
|---|---|---|
| 空化 | cavitation | 液体中压强降到蒸气压以下时出现并长大的蒸气/气体空腔 |
| 空蚀 | cavitation erosion | 空泡反复塌缩对固体表面造成的剥蚀 |
| 空化数 | cavitation number | \((p_\infty-p_v)/(\frac12\rho U^2)\),越小越易空化 |
| Rayleigh–Plesset 方程 | Rayleigh–Plesset equation | 球形气泡半径的运动方程,式 (4.2) |
| 速度势 | velocity potential | 满足 \(\vect u=\nabla\phi\) 的标量函数 |
| 随体加速度 | material acceleration | 跟随同一流体元的加速度 \(\pdv{u}{t}+u\pdv{u}{r}\) |
| 多方指数 | polytropic exponent | \(pV^\kappa=\) 常量中的 \(\kappa\) |
| 热 Péclet 数 | thermal Péclet number | \(\omega R^2/D\),比较振荡与热扩散快慢 |
| Minnaert 频率 | Minnaert frequency | 气泡体积振荡的固有频率,式 (6.2) |
| 单极子 | monopole | 只靠体积变化发声的各向同性点声源 |
| 品质因数 | quality factor \(Q\) | \(\omega E/P_{\mathrm{loss}}\),振荡衰减快慢 |
| 空泡率 | void fraction | 混合物中气体所占体积分数 \(\phi\) |
| Wood 公式 | Wood’s equation | 含泡液体低频声速公式,式 (7.1) |
| 声致发光 | sonoluminescence | 声场驱动气泡塌缩时发出的闪光 |
| 枪虾(鼓虾) | snapping shrimp | 鼓虾科(Alpheidae)甲壳动物,用大螯闭合产生空化泡 |
| Blake 阈值 | Blake threshold | 泡核失稳长大所需的临界液体压强 \(p_{\mathrm{Bl}}\) |
| 镜像法 | method of images | 用边界外对称放置的虚拟源代替边界条件;刚性壁同号、自由液面反号 |
| 偶极子 | dipole | 一对相距很近、强度相反的点源,辐射呈 \(\cos^2\theta\) 分布 |
| 涡环 | vortex ring | 首尾闭合成环的涡旋(如烟圈),环截面中心为涡核 |
| Rankine 涡 | Rankine vortex | 核内刚性转动、核外 \(v_\theta\propto1/r\) 的简化涡模型 |
| 声阻抗 | acoustic impedance | \(Z=\rho c\),两介质界面处声波反射由其比值决定 |
| 雷诺数 | Reynolds number | \(\mathrm{Re}=2\rho vR/\mu\),惯性力与粘性力之比;\(\ll1\) 时 Stokes 公式适用 |
| Keller–Miksis 方程 | Keller–Miksis equation | 计入液体一阶可压缩性(声辐射)的 RP 方程 |
| 冲击波 | shock wave | 振幅很大、波前陡峭的压缩波 |
| 微射流 | microjet | 近壁不对称塌缩时穿过泡、指向壁面的高速液流 |
当 \(n\neq2\) 时流体会被压缩或膨胀(见下),此时完整的牛顿粘性定律还含正比于 \(\nabla\cdot\vect u\) 的项,所以这里只是一个形式计算,目的仅是看清 \(n=2\) 的特殊性。↩︎