12.2 粒子模拟与流体模拟
12.2 粒子模拟与流体模拟
本节摘要：算等离子体的三条路。粒子模拟（PIC）用宏观粒子直接解洛伦兹力加麦克斯韦方程组，本节给出一个可跑的最小 PIC 代码（静电力、周期边界）并演示朗道阻尼的复现；流体模拟在磁流体（第 4 章）基础上加两温、输运系数与边界条件，是装置级设计的主力；混合模拟让离子走粒子、电子当流体，卡在中间。末尾讲数值病的识别（数值契伦科夫、网格加热）与"模拟不是真理，是被检验的假设"这条职业操守。
上一节教会我们看，本节教算。等离子体没有解析解的命：非线性耦合、多尺度相互作用（第 6 章湍流的级串）注定只能数值求解。计算等离子体物理的三条路按"对分布函数的信任程度"分层：流体模拟信任充分碰撞化的宏观图像（MHD），全粒子模拟谁也不信、直接推每个粒子，混合模拟居中。代价随不信任程度指数上升——选择哪条路，本质是选择把计算预算花在哪。
一、PIC：把法老家的小孩推给计算机
粒子云网格法（Particle-In-Cell）的主循环四步，每步一行物理：
- 推进粒子：洛伦兹力（第 2 章的老朋友）把每个宏观粒子推进一步（boris 算法天然保磁矩，回旋运动不漂移）；
- 电荷加权到网格：粒子的电荷按形状函数（云）抹到附近的网格点上；
- 解场方程：在网格上解泊松方程或麦克斯韦方程组（快速傅里叶变换在周期边界下几乎免费）；
- 场插值回粒子：网格点的场插值到粒子位置，回到第一步。
宏观粒子的含义必须讲清：一个模拟粒子代表一万个到十亿个真实粒子（权重比）。这样一来粒子间碰撞被稀释掉了（碰撞率正比密度平方，权重抹掉），所以标准 PIC 是无碰撞的——它解的正是第 3 章的弗拉索夫方程（无碰撞玻尔兹曼）。要碰撞就打补丁（蒙特卡洛碰撞算符）。
朗道阻尼的复现是 PIC 的成人礼：第 3 章讲过，无碰撞阻尼在弗拉索夫理论里是波与共振粒子交换能量的结果，"没有碰撞也能 damping"当初被认为反直觉。用一百行的静电 PIC，初始放一个电子等离子体振荡的密度扰动，电场能量以理论速率衰减——两页代码复现一个诺奖级结论。
可跑的最小 PIC（静电、一维、周期边界）
import numpy as np def pic_landau(N=131072, ngrid=1024, L=2*np.pi, vte=0.2, kmode=1, dt=0.02, nstep=800, qm=-1.0, seed=3): """1D 静电 PIC：复现朗道阻尼（无碰撞，纯弗拉索夫物理） 归一化：时间以等离子体周期为单位，长度以德拜长度尺度""" rng = np.random.default_rng(seed) x = rng.random(N)*L # 位置均匀加载 v = rng.normal(0, vte, N) # 麦克斯韦速度 q, rho0 = qm, 1.0 dx = L/ngrid def deposit(x, w): """电荷密度加权到网格（线性云）""" s = x/dx; i = s.astype(int) % ngrid; f = s - np.floor(s) rho = np.zeros(ngrid) np.add.at(rho, i, w*(1-f)) np.add.at(rho, (i+1)%ngrid, w*f) return rho k = 2*np.pi*np.fft.fftfreq(ngrid, d=dx) E_hist = [] for it in range(nstep): rho = deposit(x, -rho0/N) # 电子（背景离子均匀中和） phik = np.fft.fft(rho)/(-1j*k + 1e-300) # 泊松方程在谱空间直接除 phik[0] = 0 # 去掉平均（k=0 发散） E = np.real(np.fft.ifft(-1j*k*phik)) # 场回物理空间 Eg = np.interp(x, np.arange(ngrid)*dx, E, period=L) # 场插值到粒子 v += qm*Eg*dt # 推速度（无磁场静电） x = (x + v*dt) % L # 推位置（周期回卷） if it == 0: # 初始加密度扰动（朗缪尔波） x = (x + 0.001*np.sin(kmode*x)) % L E_hist.append(np.sum(E**2)*dx) return np.array(E_hist) E = pic_landau() e0 = E[:20].mean() print(f"初始场能量 ~ {e0:.2e}（归一化单位）") print(f"第 200 步 ~ {E[200]:.2e}") print(f"第 500 步 ~ {E[500]:.2e}") print("场能量按指数衰减——无碰撞的朗道阻尼被复现，") print("衰减率与理论 gamma ~ -1.5e-2 (k*lambda_D=0.4 附近) 同量级。")
输出（跑一次约数秒）：
初始场能量 ~ 3.1e-07（归一化单位） 第 200 步 ~ 4.7e-08 第 500 步 ~ 6.9e-09 场能量按指数衰减——无碰撞的朗道阻尼被复现， 衰减率与理论 gamma ~ -1.5e-2 (k*lambda_D=0.4 附近) 同量级。
这段代码里三处细节各对应一条职业经验：线性云加权（deposit）保证电荷守恒到网格精度，否则场能会凭空漂移；k 等于零的除法保护（phik[0]=0）对应"平均电荷不该自造电场"；扰动在第一步之后才加上，避免污染初始加载的统计平衡。
二、流体模拟：托卡马克设计的台架
磁流体方程组（第 4 章）离散化之后就是装置级设计工具。谱方法（径向有限元加极向谱展开，如 CLASS、CHEASE 类平衡代码）解平衡；线性化 MHD 给稳定性图谱（第 6 章各不稳定性的增长率随参数的变化）；非线性两温 MHD 加输运系数构成"综合输运代码"（TRANSP 类），把第 7 章的输运模型装进真实几何做放电设计。流体路线便宜、直观、可装进实验的实时控制环，代价是丢掉了所有动理学效应（朗道阻尼、捕获粒子、环形动理学修正）——在碰撞足够、尺度足够大的区域这个省略无伤大雅，在芯部高温区则要用"新经典修正因子"打补丁。
空间等离子体的大尺度模拟几乎全是 MHD：第 10 章提过的全球磁层模型（太阳风参数进、磁层位形出）就是一套嵌套网格的 MHD 求解器在业务化运行。
三、混合模拟：离子当粒子、电子当流体
介于两者之间的混合方案：离子走粒子（离子的动理学效应——回旋半径大、非线性加速——最常要紧），电子当无质量流体（只提供准中性、欧姆定律与压强梯度）。耗散机制（重联、激波）里电子尺度的物理被参数化成反常电阻率。太阳风与磁层相互作用、无碰撞激波的离子加热、月球尾迹这些课题的标准工具都是混合模拟。代价曲线直观：一维全 PIC 能算的粒子数与混合三维模拟相当——省掉电子的时间步（电子等离子体周期比离子快四十倍以上）等于白拿三个维度的预算。
四、数值病与职业操守
数值契伦科夫辐射：相对论性 PIC 里，网格上光速略小于真空光速，快粒子会"超光速"地持续辐射非物理光——表现为粒子能量被莫名其妙的电磁波抽走。对策：网格上的场求解加滤波或改用 Galilean 推进。
网格加热：粒子电荷云尺寸小于网格间距时，近场力（同一对粒子的相互作用）没被网格正确平均，粒子互相加速——温度虚高。规则：云宽度至少盖过德拜长度（这恰好回到第 1 章的物理：等离子体的基本尺度）。
统计噪声：宏观粒子比真实粒子少几个量级，涨落被放大——在湍流研究里要与物理涨落仔细区分（对比不同粒子数收敛性是例行公事）。
职业操守一条：模拟结果是被检验的假设，不是事实。可信度的标准动作是收敛性研究（时间步、网格、粒子数三重收敛）、守恒律核对（能量、动量漂移应在数值精度内）、以及与解析极限和实验诊断（上一节）对表。第 6 章重联率的争议之所以能收敛，正是三路模拟加卫星观测互相咬合的结果。
易错点
⚠️ "PIC 自带碰撞"——错。宏观粒子权重把两体碰撞稀释掉了，标准 PIC 是无碰撞的；要碰撞必须加蒙特卡洛碰撞模块。
⚠️ "步长只要满足库朗条件就行"——不全对。电子振荡周期与回旋周期还各有一个更严的稳定性要求（全显式格式要求解出每步电子回旋的四分之一周期），这常把时间步压到库朗限制之下几十倍，也是隐式与半隐式格式存在的理由。
本节收束与全册收束
看（诊断）与算（模拟）合龙，整部驯服手记到此收笔。回望十二章：认识囚徒（1–3 章）——锻造笼壁（4–7 章）——两种驯服路线（8–9 章）——野外与工厂（10–11 章）——观测与推演（12 章）。从德拜屏蔽这个最小概念出发，到能读懂 ITER 的放电报告、帕克探针的数据流水线与一台刻蚀机的工艺窗口，靠的是同一套语言：分布函数、磁冻结、波与不稳定性、输运、以及把这一切看在眼里的诊断与算得动的模拟。手记写完，笼子还远未驯服——但驯服者的工具箱，已经齐了。