A reproducible numerical laboratory

Source: supplied generic textbook, Chapter 13. The complete bilingual chapter follows, with its derivations and qualifications. Textbook source (Markdown) · Original Markdown.

来源：所给通用教材第 13。下方完整收录双语章节及推导和限定。教材原文（Markdown） · 原始 Markdown。

This lab is placed after Chapters 1–7, as the original author recommends. Its code and original data are preserved; a separate rerun report is linked below.

按原作者建议，数值实验放在第 1—7 章之后。原程序和数据完整保留；下方链接独立重跑报告。

13. A reproducible numerical laboratory 13.1 What the supplied program actually computes 13.2 Scaling the equations before integration 13.3 Running the laboratory 13.4 Verification results supplied with this edition 13.5 What not to do when the solver becomes difficult 13.6 Moving from ODEs to field simulation

13. A reproducible numerical laboratory

可复现数值实验

13.1 What the supplied program actually computes

配套程序实际计算什么

The file code/bubble_lab.py implements an incompressible Rayleigh–Plesset model, a constant-ambient Keller–Miksis model, and the empty-cavity Rayleigh benchmark. It integrates a two-state system $(R,U)$ with $U=\dot R$ . The gas follows a fixed polytropic law; vapor pressure and liquid properties are constant. The demonstration starts from an expanded bubble at rest, not from laser nucleation.

文件code/bubble_lab.py实现不可压缩Rayleigh–Plesset模型、恒定环境压力下的Keller–Miksis模型，以及空腔Rayleigh基准。它积分由 $(R,U)$ 构成的两状态系统，其中 $U=\dot R$ 。气体服从固定多方关系，蒸气压和液体性质恒定。演示从膨胀后的静止气泡开始，而不是从激光成核开始。

The default reference radius is 10 micrometers. The comparison figures use $R(0)=1.8R_0$ and $U(0)=0$ , with identical initial states for both models. This deliberately modest case stays well below the chosen Mach guard. Stronger collapse is used only to verify that the guard actually stops integration.

默认参考半径为10微米。对比图采用 $R(0)=1.8R_0$ 和 $U(0)=0$ ，两个模型使用完全相同的初态。这个有意选择的温和案例，始终明显低于所设Mach保护阈值。更强塌缩案例仅用于检验保护机制确实能够终止积分。

13.2 Scaling the equations before integration

积分前先进行尺度化

The public interface uses SI units, but the solver integrates dimensionless variables $x=R/R_0$ , $v=U/U_*$ , and $\tau=tU_*/R_0$ , where $U_*=\sqrt{(p_0-p_v)/\rho}$ . Then

程序对外接口采用SI单位，但求解器积分无量纲变量 $x=R/R_0$ 、 $v=U/U_*$ 和 $\tau=tU_*/R_0$ ，其中 $U_*=\sqrt{(p_0-p_v)/\rho}$ 。于是

\frac{dx}{d\tau}=v,\qquad \frac{dv}{d\tau}=\frac{R_0}{U_*^2}\,a(R_0x,U_*v), \tag{13.1}

13.2 Scaling the equations before integration / 积分前先进行尺度化; x: position / dimensionless radius; τ: duration / dimensionless time / dummy time; v: velocity / squared wall speed as defined; R_0: reference equilibrium radius; U: velocity / velocity scale; a: droplet radius / film span / acceleration

where $a$ is equation (3.4) or (7.3) solved for acceleration. Scaling prevents the same absolute tolerance from treating a micrometer-scale radius and a meter-per-second velocity in incompatible ways. It also makes parameter comparisons easier to interpret.

其中 $a$ 由式(3.4)或(7.3)整理为加速度形式。尺度化可以避免同一个绝对容差对微米级半径与米每秒级速度产生不相容的控制，同时也便于解释参数对比。

The implementation uses adaptive DOP853 integration with dense output. Its maximum internal step is additionally bounded relative to the linear period. Dense plotting output is not the same thing as a small integration step: drawing more interpolated points does not repair an inaccurate solution.

程序使用带稠密输出的自适应DOP853积分，并额外按线性周期限制最大内部步长。稠密绘图输出不等于较小积分步长：增加插值绘图点数，无法修复本身不准确的解。

13.3 Running the laboratory

运行实验

From the package directory, create a Python environment with NumPy, SciPy, and Matplotlib, then run the command below. The tested package versions are recorded in data/verification.json; requirements.txt records the environment used to generate the supplied results.

在交付包目录中，建立安装NumPy、SciPy和Matplotlib的Python环境，然后运行以下命令。实际测试的软件版本记录在data/verification.json中，requirements.txt记录生成本次结果所使用的环境。

python code/bubble_lab.py --output data --plots

The command writes CSV trajectories, a JSON verification report, and seven separate figures. The original paper PDFs are not needed to run the code. Running the command does not validate a physical experiment; it reproduces the teaching calculations and mathematical checks.

该命令输出CSV轨迹、JSON验证报告和七幅独立图像。运行代码不需要原始论文PDF。执行命令并不意味着验证了某个物理实验，而是复现教学计算和数学检查。

For a parameter experiment, save the following script in the package root. It changes viscosity while leaving the gas-content definition and all other parameters unchanged. The printed stop reason should always be inspected before interpreting extrema.

进行参数实验时，可将以下脚本保存在交付包根目录中。它只改变黏度，保持气体含量定义及其他参数不变。在解释任何极值之前，都应检查输出中的停止原因。

from dataclasses import replace
from code.bubble_lab import Parameters, simulate

base = Parameters()
for viscosity in (0.5e-3, 1.0e-3, 2.0e-3):
    p = replace(base, mu=viscosity)
    result = simulate(p, model="KM", radius_factor=1.8, cycles=6)
    print(viscosity, result["stop"], result["max_mach"])

The package includes code/__init__.py for this import pattern. In an interactive environment where Python's standard-library module named code is already imported, instead add the package's code directory to sys.path and import bubble_lab directly, as shown in the README.

交付包包含支持这种导入方式的code/__init__.py。若在交互式环境中已经导入了Python标准库中同名的code模块，则应按README中的方法，把交付包的code目录加入sys.path，再直接导入bubble_lab。

13.4 Verification results supplied with this edition

本版附带的验证结果

The following checks were actually executed. Errors refer to the stated normalization, not to experimental accuracy. Small numerical values demonstrate consistency with the implemented mathematical problem; they do not establish that the closure laws represent a real laser bubble.

以下检查均已实际执行。误差采用所列归一化方式，并非实验精度。较小数值说明计算与所实现数学问题相一致，但不能证明闭合关系准确代表真实激光气泡。

Check / 检查	Observed result / 实测计算结果
RP equilibrium radius error / RP平衡半径相对误差	$0$
KM equilibrium radius error / KM平衡半径相对误差	$0$
Small-amplitude period relative error / 小振幅周期相对误差	$3.97\times10^{-7}$
Inviscid RP energy drift / 无黏RP能量相对漂移	$7.29\times10^{-11}$
Viscous energy-balance residual / 黏性能量平衡相对残差	$6.33\times10^{-9}$
KM versus RP as $c\to\infty$ , radius difference divided by $R_0$ / 无限声速极限半径差除以 $R_0$	$1.26\times10^{-10}$
KM tolerance refinement, radius difference divided by $R_0$ / KM容差加密半径差除以 $R_0$	$3.59\times10^{-9}$
Rayleigh finite-cutoff time relative error / Rayleigh有限截止半径时间相对误差	$4.16\times10^{-12}$
Strong-collapse terminal Mach number / 强塌缩停止时Mach数	$0.1000$

The equilibrium checks test signs and initialization. The small-amplitude test compares a numerically measured period with equation (6.2), using roots of the velocity interpolation rather than coarse sampled maxima. The energy tests probe the full nonlinear trajectory. The infinite-sound-speed check tests the algebra of the compressible extension. The guard test verifies that invalid-regime reporting is not silently bypassed.

平衡检查检验符号和初始化。小振幅测试将数值得到的周期与式(6.2)比较，并通过速度插值的零点而不是粗采样极大值确定周期。能量测试检查完整非线性轨迹。无限声速检查检验可压缩扩展的代数实现。保护测试则确认程序不会默默越过无效区域而继续报告结果。

13.5 What not to do when the solver becomes difficult

求解困难时不要做什么

Do not replace negative trial radii by a tiny positive number and continue. Do not cap pressure at an arbitrary value merely to avoid overflow. Do not interpret a dense time trace as convergence evidence. Instead inspect the physical regime, event location, internal step size, scaling, and closure. A singular model cannot be repaired by a smaller time step alone.

不要将负的试探半径替换为极小正数后继续计算；不要只为了避免溢出而任意截断压力；也不要把稠密时间曲线当作收敛证据。应检查物理范围、停止事件、内部步长、尺度化和闭合关系。仅减小时间步长，不能修复一个奇异模型。

For a Rayleigh collapse benchmark, compare numerical and analytic times at the same finite cutoff radius. Comparing a stopped numerical trajectory against the analytic time to zero mixes two different endpoint problems. The code performs the correct finite-cutoff comparison.

对于Rayleigh塌缩基准，应在相同有限截止半径下比较数值时间和解析时间。若把提前停止的数值轨迹与解析到零半径的时间比较，就混合了两个不同终点问题。配套代码执行的是正确的有限截止比较。

13.6 Moving from ODEs to field simulation

从常微分方程走向场仿真

A practical hierarchy is: radial ODE for an isolated bubble; axisymmetric interface-resolving flow for an aligned bubble and film; fully coupled axisymmetric fluid–structure interaction for compliant loading; and three-dimensional compressible flow when symmetry breaking or shock propagation controls the target quantity. The fully compressible rigid-wall source [S418] is a reminder that jet impact and local pressure require a different numerical representation from a radius history.

一个实用层次是：用径向常微分方程研究孤立气泡；用轴对称界面解析流动研究同轴气泡与薄膜；用完整轴对称流固耦合研究柔性加载；当对称性破缺或激波传播控制目标量时，再使用三维可压缩流动。完全可压缩刚壁研究[S418]提醒我们，射流冲击和局部压力需要不同于半径历程的数值表达。

An arbitrary Lagrangian–Eulerian moving mesh can follow a smooth interface but struggles when it pinches off or changes topology. Volume-of-fluid or level-set methods can accommodate topology changes but introduce interface-resolution and conservation issues. Axisymmetry suppresses non-axisymmetric modes by construction. Every method trades one set of difficulties for another; no method name alone guarantees accuracy.

任意Lagrange–Euler移动网格可以跟踪平滑界面，但在颈缩断裂或拓扑改变时会面临困难。流体体积分数法或水平集法可以处理拓扑变化，但会引入界面分辨率和守恒问题。轴对称模型则从构造上排除了非轴对称模态。每种方法都在不同困难之间取舍，任何方法名称本身都不保证精度。

For field calculations, check mesh refinement, time refinement, domain-size sensitivity, outgoing-wave treatment, interface conservation, and structural energy transfer. Numerical verification asks whether equations were solved correctly; experimental validation asks whether those equations describe the physical system adequately. Report them separately.

对于场计算，应检查网格加密、时间加密、域尺寸敏感性、出射波处理、界面守恒和结构能量传递。数值验证关注方程是否被正确求解；实验验证关注方程是否充分描述物理系统。二者应分开报告。

Verified rerun of the original program

原程序核查重跑

The original program was rerun on 6 October 2026 without editing its equations or parameter defaults. All checks passed. Results below describe numerical consistency of the supplied model, not experimental accuracy.

2026 年 10 月 6 日重跑原程序，没有修改方程或默认参数。全部检查通过。下列结果描述模型数值一致性，并非实验精度。

Original solver / 原求解器 · New report / 新报告 · Author’s original report / 原作者报告

Check / 检查	Rerun / 重跑
Small-amplitude period / 小振幅周期	3.97007324571e-07
Inviscid energy drift / 无黏能量漂移	7.29248082186e-11
Viscous balance residual / 黏性能量残差	6.32565819262e-09
Infinite-c limit / 无限声速极限	1.2621045024e-10
Tolerance refinement / 容差加密	3.58732691786e-09
Rayleigh coefficient / Rayleigh 系数	0.91468135649
Finite-cutoff time / 有限截止时间	4.16014547544e-12

Strong collapse stopped at Mw = 0.10 as intended. The modest RP/KM ringdowns remained below Mw = 0.018. Read the exported trajectories alongside the book’s seven original plots.

强塌缩按预期在 Mw = 0.10 停止。温和 RP/KM 振荡的 Mw 均小于 0.018。可结合教材七幅原图阅读导出轨迹。

rp_ringdown.csv

km_ringdown.csv

rayleigh_dimensionless.csv

Integrated calculation and model check

综合计算与模型检查

Use the declared parameter set in each part. The values are teaching scenarios unless specifically identified as inherited project inputs. Each reference panel follows its question.

使用每小题明确给出的参数。除明确标为沿用项目输入外，数值均为教学情景。各参考解答紧随题目。

Run or inspect the saved lab. For R0 = 10 µm compute its velocity/time scales, noting that the code uses c = 1481 m/s. Explain why dense plotted points do not substitute for solver convergence.
运行或查看保存的数值实验。R0 = 10 µm 时计算速度和时间尺度，注意程序 c = 1481 m/s。解释稠密绘图点为何不能替代求解器收敛。
Reference answer / 参考答案
$U_*=\sqrt{(p_0-p_v)/\rho},\quad t_*=R_0/U_*,\quad x=R/R_0,\quad v=\dot R/U_*$
The diagram locates the physical system; its key lists the symbols in the formula. Read the definition in the surrounding derivation when a letter has different roles.
图示标出物理体系，图例列出公式中的符号。若某字母有多种作用，应以邻近推导中的定义为准。
U = 9.959 m/s*, t = 1.004 µs*. The rerun uses the original code constants without silently changing c to 1482. Dense output samples an interpolant; tolerances, internal step limits and independent benchmarks test the underlying integration.
U = 9.959 m/s*，t = 1.004 µs*。重跑保留原程序常数，没有悄悄将 c 改为 1482。稠密输出采样的是插值函数；容差、内部步长限制及独立基准检验的是实际积分。
Read the new verification report. Quote the Rayleigh coefficient, inviscid energy drift and viscous balance residual, and state what all_checks_passed proves.
阅读新验证报告，列出 Rayleigh 系数、无黏能量漂移和黏性能量平衡残差，并说明 all_checks_passed 证明什么。
Reference answer / 参考答案
$\epsilon_E=\max_t|E(t)-E(0)+\int_0^t16\pi\mu R\dot R^2dt|/|E(0)|$
The diagram locates the physical system; its key lists the symbols in the formula. Read the definition in the surrounding derivation when a letter has different roles.
图示标出物理体系，图例列出公式中的符号。若某字母有多种作用，应以邻近推导中的定义为准。
The rerun gives coefficient 0.9146813564901609, inviscid relative drift 7.29×10⁻¹¹, viscous residual 6.33×10⁻⁹. All deterministic implementation checks passed. This is numerical verification of the specified spherical closures, not experimental validation of PFC activation, jetting or film release.
重跑得到系数 0.9146813564901609、无黏相对漂移 7.29×10⁻¹¹、黏性残差 6.33×10⁻⁹。所有确定性实现检查通过。这验证了指定球形闭合的数值实现，并非 PFC 激活、射流或膜释放的实验验证。
The strong-collapse run stops at Mach 0.10. A radius plot looks smooth to that point. Should you remove the event and extrapolate a peak? Compare a cutoff benchmark with the analytic time to zero.
强塌缩计算在 Mach 0.10 停止，之前半径曲线看起来平滑。应移除事件并外推峰值吗？比较截止基准与解析到零时间。
Reference answer / 参考答案
$M_w=|\dot R|/c,\qquad t_{\mathrm{cut}}=t_*\sqrt{3/2}\int_{x_{\mathrm{cut}}}^1\frac{x^{3/2}}{\sqrt{1-x^3}}dx$
The diagram locates the physical system; its key lists the symbols in the formula. Read the definition in the surrounding derivation when a letter has different roles.
图示标出物理体系，图例列出公式中的符号。若某字母有多种作用，应以邻近推导中的定义为准。
Retain the regime guard and report its stop reason. Compare numerical and analytic times at the same finite radius. Smoothness is not evidence for the validity of post-guard physics; do not clip trial radii or pressures to manufacture a late-time maximum.
保留适用范围保护条件并报告停止原因。应在同一个有限半径比较数值与解析时间。平滑性不能证明保护条件之后物理有效；不要截断试探半径或压力来制造后期最大值。

Brief conceptual synthesis

简短概念综合

After using the reference panels, give a short explanation of the physical mechanism, the strongest assumption, and one observation that could falsify this chapter’s model. Preparation of later chapters does not require waiting for that reply.

核对参考解答后，简述物理机制、最强假设，以及一个可证伪本章模型的观测。后续章节编写无需等待该回答。