ECHO3D尾场仿真完全教程：从物理原理到代码实现

Posted on In Views: Disqus:
前置知识：基础电磁学、加速器物理入门、Python 基础

引言：什么是尾场，为什么要仿真它

在现代粒子加速器中，一束高能带电粒子穿过真空管道、加速腔、波纹管等结构时，会在结构中激发出电磁场。这个被激发的场会反过来作用在束团自身（或后续束团）上，影响粒子的能量和轨迹——这就是尾场（Wakefield）。

尾场效应是现代加速器设计的核心问题之一：

  • 纵向尾场 $W_\parallel(s)$：使束团头部粒子减速、尾部粒子加速（或反之），导致能散增大
  • 横向尾场 $W_\perp(s)$：使偏离轴心的束团受到横向踢力，可能导致束流不稳定性（Beam Break-Up, BBU）

因此，在设计任何加速器组件时，精确计算尾场效应是必不可少的。ECHO3D 正是为此而生——它是一个基于时域有限差分法（FDTD）的三维尾场仿真程序。

物理基础：尾场、阻抗与损失因子

尾场势（Wake Potential）

考虑一个以光速 $c$ 运动的点电荷 $q$ 穿过一个任意结构。在它后面距离 $s$ 处（$s>0$ 表示在激发粒子之后），一个试验电荷感受到的纵向尾场势定义为：

$$W_\parallel(s) = -\frac{1}{q}\int_{-\infty}^{\infty} E_z\left(z, t=\frac{z+s}{c}\right) dz$$

类似地，横向尾场势定义为：

$$W_\perp(s) = \frac{1}{q}\int_{-\infty}^{\infty} \left[\mathbf{E}_\perp + c(\hat{z} \times \mathbf{B}\perp)\right]{t=\frac{z+s}{c}} dz$$

损失因子（Loss Factor）与踢因子（Kick Factor）

对于有限长度的束团，其电荷分布为 $\lambda(s)$，则：

  • 损失因子（纵向）：

$$\kappa_\parallel = -\int_{-\infty}^{\infty} \lambda(s) \cdot W_\parallel(s) , ds \quad [\text{V/pC}]$$

物理意义：束团中单位电荷的平均能量损失。

  • 踢因子（横向）：

$$\kappa_\perp = \int_{-\infty}^{\infty} \lambda(s) \cdot W_\perp(s) , ds \quad [\text{V/pC/m}]$$

物理意义：单位电荷偏移量受到的横向踢力。

阻抗（Impedance）

尾场势的傅里叶变换就是阻抗：

$$Z_\parallel(\omega) = \int_{-\infty}^{\infty} W_\parallel(s) , e^{-i\omega s/c} , \frac{ds}{c}$$

$$Z_\perp(\omega) = i \int_{-\infty}^{\infty} W_\perp(s) , e^{-i\omega s/c} , \frac{ds}{c}$$

阻抗的实部对应能量损失（电阻性），虚部对应相移（电抗性）。阻抗峰的位置对应结构的本征模式频率。

Panofsky-Wenzel 定理

一个极其重要的定理将纵向和横向尾场联系起来：

$$\frac{\partial W_\parallel}{\partial \mathbf{r}_\perp} = \frac{\partial \mathbf{W}_\perp}{\partial s}$$

这意味着：如果你知道纵向尾场在横向上的梯度，你就可以推导出横向尾场。ECHO3D 正是利用这个定理，先计算”间接尾场”（indirect wake），再通过数值微分得到横向分量。

ECHO3D 简介与时域有限差分法

什么是 ECHO3D

ECHO3D（Electromagnetic Computation Heavy Operations 3D）是一个用 Fortran 编写的三维电磁场时域求解器，专门针对加速器尾场问题优化。它的核心特点：

特性 说明
算法 时域有限差分法（FDTD），隐式格式
几何 STL 三角面片导入，支持共形网格
边界条件 PEC、PMC、Open、PML
束团模型 高斯纵向分布，点状横向
输出 二进制尾场数据 + 文本场监视器

FDTD 基本原理

FDTD 方法在时域中直接离散麦克斯韦方程组：

$$\nabla \times \mathbf{E} = -\mu \frac{\partial \mathbf{H}}{\partial t}$$

$$\nabla \times \mathbf{H} = \varepsilon \frac{\partial \mathbf{E}}{\partial t} + \mathbf{J}$$

在 Yee 网格上，电场和磁场在空间和时间上交错排列，通过蛙跳（leapfrog）格式递推求解。

ECHO3D 的”间接尾场”方法

ECHO3D 不直接计算尾场势 $W(s)$，而是计算一个中间量——间接尾场（indirect wake）。其核心思想是：

  1. 在束团路径上放置一系列”监视器”点
  2. 记录每个监视器点处的纵向电场 $E_z(t)$
  3. 通过时域积分得到间接尾场
  4. 后处理时通过 Panofsky-Wenzel 定理的数值形式计算横向尾场

ECHO3D 工作流：从几何到结果

四步流水线

ECHO3D 的仿真流程由四个独立的可执行程序组成：

flowchart LR
    A["1️⃣ Mesher.exe
网格生成"] --> B["2️⃣ InitField.exe
初始场设置"]
    B --> C["3️⃣ ECHO3D.exe
时域推进求解"]
    C --> D["4️⃣ IndirectIntegration.exe
尾场积分"]
    
    A1["input.txt
geometry.txt
*.stl"] --> A
    D --> E["wake3D.bin
wake3Dindirect.bin
Monitor_N*.txt"]

各步骤详解

步骤 1：Mesher.exe — 网格生成

读取 geometry.txt 中定义的 STL 几何和材料，生成 Yee 网格。输出放在 Mesh/ 目录下：

  • Mesh/Body_00/Eps.mat — 介电常数分布
  • Mesh/Body_00/Mue.mat — 磁导率分布
  • Mesh/Body_00/Steps.mat — 网格步长信息

步骤 2：InitField.exe — 初始场设置

在束团初始位置设置激励源（高斯束团产生的自由空间场或波导模式场）。输出放在 Fields/ 目录下。

步骤 3：ECHO3D.exe — 时域推进

这是计算的核心。在时域中推进麦克斯韦方程组，束团以光速运动，在每个时间步更新电磁场。输出：

  • OutFields/ — 各时间步的场快照（如果 DumpField=1）
  • Fields/ — 更新后的场文件
  • Monitor_N*.txt — 场监视器输出（如果配置了 FieldMonitor）

步骤 4：IndirectIntegration.exe — 尾场积分

对监视器记录的场进行时域积分，得到间接尾场。输出：

  • Results/wake3D.bin — 直接尾场
  • Results/wake3Dindirect.bin — 间接尾场（后处理的主要输入）

命令行运行方式

在 Windows 下，典型运行命令为：

1
cmd /c "cd /d ECHO3D_Minimal && Mesher.exe input.txt && InitField.exe input.txt && ECHO3D.exe input.txt && IndirectIntegration.exe input.txt"

四个程序按顺序执行，&& 确保前一步成功后才执行下一步。

输入文件详解

主输入文件 input.txt

以最小纵向尾场算例（ECHO3D_Minimal）为例：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
%%%%%%%%%%%%%% geometry %%%%%%%%%%%%%%%%%
GeometryFile = 'geometry.txt'       % 几何定义文件路径
Units = 'mm'                        % 单位：m / cm / mm
BoundaryConditionsX = [0  0]        % X方向边界：0=开放
BoundaryConditionsY = [0  1]        % Y方向边界：0=磁壁, 1=电壁
BoundaryConditionsZ = [0  1]        % Z方向边界：0=磁壁, 1=电壁

%%%%%%%%%%%%%% beam and field %%%%%%%%%%%%%%%%%
BunchSigma = 1                      % 束团RMS长度 [Units]
BunchPosition = [0  0]              % 束团横向位置 [网格步数]
InFieldDir = '-'                    % 初始场类型：'free'=自由空间, '-'=波导

%%%%%%%%%%%%%% mesh %%%%%%%%%%%%%%%%%
TimeSteps = -1                      % 时间步数，-1=自动确定
MeshLength = 200                    % 纵向网格数 nx
dY = [0 10]                         % Y方向范围 [Units]
dZ = [0 10]                         % Z方向范围 [Units]
Steps = [0.2 0.2 0.2]              % 步长 dx dy dz [Units]
Tolerance = 0.01                    % 隐式求解器容差
PMLDepth = 0                        % PML层数（0=无PML）
PMLParameters = [1 1]              % PML参数

%%%%%%%%%%%%%% solver %%%%%%%%%%%%%%%%%
SolverType = 'impl'                 % 求解器：'expl'=显式, 'impl'=隐式
Conformal = 1                       % 共形网格：0=USC, 1=Simple
Iterations = 0                      % 最大迭代次数（0=自动）
InitialIterations = 0               % 初始迭代次数
Damping = 0                         % 阻尼系数 [0, 0.5]
ThreadsNumber = 4                   % OpenMP 线程数

%%%%%%%%%%%%%% monitors %%%%%%%%%%%%%%%%%
FieldMonitor = {'Ez' 'x' 20 100 0 10 0 10 0 100 1}   % 场监视器（x-时间扫描）
FieldMonitor = {'Ez' 's' 20 100 0 10 0 10 0 100 1}   % 场监视器（s-空间扫描）
DumpMesh = 1                        % 是否导出网格
DumpField = 1                       % 是否导出快照

关键参数说明

参数 含义 典型值
BunchSigma 束团 RMS 长度 1 mm（短束团）~ 20 mm（长束团）
BunchPosition 束团在横截面中的偏移 [ny_offset, nz_offset] [0, 0] 轴上，[1, 0] 偏轴
MeshLength 纵向网格数，决定尾场计算窗口长度 200~520
Steps 三个方向的网格步长 越细越精确，但计算量立方增长
SolverType 'impl' 隐式无条件稳定，推荐；'expl' 显式需满足 CFL 条件

场监视器 FieldMonitor 语法

1
FieldMonitor = {'字段名' '扫描类型' k1 k2 k3 k4 k5 k6 k7 k8 k9}
  • 字段名：Ex, Ey, Ez, Hx, Hy, Hz
  • 扫描类型：'x' 沿 x 方向时间扫描，'s' 固定空间位置扫描
  • k1-k9：9 个整数参数，定义监视器的空间范围和采样点数

尾场监视器 WakeMonitor（用于横向尾场）

1
WakeMonitor = [ymin ymax zmin zmax]   % 以束团位置为中心的网格线偏移范围

例如 WakeMonitor = [-1 1 -1 1] 表示在束团周围 ±1 个网格线范围内放置监视器。

几何文件 geometry.txt

1
2
3
4
5
6
7
8
9
10
11
12
13
Background = 0 0 0                  % 背景材料 [eps, mu, sigma]

%%%%%%%%%%% Materials %%%%%%%%%%%%%%%%%%
MaterialsNumber = 1                  % 材料种类数
stepout.stl 1 1 0                   % STL文件 eps mu sigma

%%%%%%%%%%% Meshing parts %%%%%%%%%%%%%%
MeshParts = 1                        % 网格分区数
body -0.3 0.3                       % 分区名 x_min x_max

%%%%%%%%%%% Geometry list %%%%%%%%%%%%%%
GeometryParts = 1                    % 几何部件数
body 1                               % 部件名 材料编号
  • Background：背景材料的 $[\varepsilon_r, \mu_r, \sigma]$，0 0 0 表示真空
  • STL 行：文件名 eps mu sigma，1 1 0 表示 PEC（理想导体）
  • Meshing parts：定义纵向分段，body -0.3 0.3 表示从 $x=-0.3$ 到 $x=0.3$（单位由 input.txt 中的 Units 决定）

二进制输出格式与解析

wake3D.bin / wake3Dindirect.bin 格式

ECHO3D 的二进制尾场文件采用 Fortran 风格存储：

1
2
3
4
5
6
7
字节偏移      类型         内容
─────────────────────────────────────
0-3          int32        nx（纵向网格数）
4-7          int32        ny（Y方向网格数）
8-11         int32        nz（Z方向网格数）
12-43        float64×4    y[0], y[-1], z[0], z[-1]（坐标边界）
44-...       float64×N    W 数据，N = nx × ny × nz

数据排列顺序（Fortran 列优先）：最外层循环是 $x$（纵向），然后是 $y$，最内层是 $z$。

索引公式：
$$\text{index} = ix \times ny \times nz + iy \times nz + iz$$

用 Python 解析二进制文件

最简洁的解析方式（来自 parse_wake.py）：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
import struct
from pathlib import Path

def parse_wake_file(bin_path: Path, out_path: Path, s0_cm=-0.48, ds_cm=0.02):
    data = bin_path.read_bytes()
    
    # 读取头部：nx, ny, nz
    n  = struct.unpack_from('<i', data, 0)[0]   # 纵向网格数
    ny = struct.unpack_from('<i', data, 4)[0]   # Y方向网格数
    nz = struct.unpack_from('<i', data, 8)[0]   # Z方向网格数
    
    values_per_sample = ny * nz
    header_size = 44  # 3×int32 + 4×float64 = 12 + 32 = 44 bytes
    
    offset = header_size
    for i in range(n):
        # 读取一个纵向切片的所有 (ny×nz) 个值
        vals = struct.unpack_from(f'<{values_per_sample}d', data, offset)
        avg = sum(vals) / values_per_sample   # 取平均
        s_cm = s0_cm + i * ds_cm
        print(f"{i} {s_cm:.6f} {avg:.17g}")
        offset += values_per_sample * 8

场监视器文本格式

Monitor_N01.txt 等文件的格式：

1
2
3
4
5
% Field=Ez time=x k_ct=200 h_ct=0.02 ct0=0
% k_y=50 h_y=0.02 y0=0
% k_z=50 h_z=0.02 z0=0
% k_s=100 h_s=0.02 s0=0
[数据矩阵：ks×ky 或 kt×ky 的浮点数]

核心函数库 pyLib4ECHO.py 逐函数讲解

pyLib4ECHO.py 是所有后处理脚本的基石，提供了 10 个核心函数。

ReadInput(InputFile) — 读取输入文件

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
def ReadInput(InputFile):
    with open(InputFile, 'r') as fid:
        lines = fid.readlines()
    
    variables = {}
    for line in lines:
        line = line.strip()
        if not line or line.startswith('%'):  # 跳过空行和注释
            continue
        line = line.split('%')[0].strip()     # 去除行内注释
        
        match = re.match(r"(\w+)\s*=\s*(.*)", line)
        if match:
            var_name, var_value = match.groups()
            # 处理 MATLAB 风格数组 [1 2 3] → [1, 2, 3]
            if re.match(r"\[.*\]", var_value):
                var_value = re.sub(r'\s+', ',', var_value)
            try:
                variables[var_name] = ast.literal_eval(var_value)
            except:
                variables[var_name] = var_value
    
    # 提取关键变量并计算派生量
    nx = variables.get('MeshLength')
    hx, hy, hz = variables.get('Steps', [None, None, None])
    dY = variables.get('dY', [0, 0])
    dZ = variables.get('dZ', [0, 0])
    BunchPosition = variables.get('BunchPosition', [0, 0])
    
    ny = round((dY[1] - dY[0]) / hy) + 1 if hy else None
    nz = round((dZ[1] - dZ[0]) / hz) + 1 if hz else None
    iy0 = BunchPosition[0] + 1   # 转为 1-based 索引
    iz0 = BunchPosition[1] + 1
    
    return nx, ny, nz, hx, hy, hz, sigma, iy0, iz0, BCy, BCz, GFile, Ymin, Zmin

设计要点：

  • 使用 ast.literal_eval 安全解析 Python 字面量
  • 自动将 MATLAB 风格数组 [1 2 3] 转为 Python 风格 [1, 2, 3]
  • 从 dY、dZ 和步长自动计算 ny、nz

ReadDataB(fname) — 读取二进制尾场数据

1
2
3
4
5
6
7
def ReadDataB(fname):
    with open(fname, 'rb') as f:
        nx, ny, nz = np.fromfile(f, dtype=np.int32, count=3)
        y = np.fromfile(f, dtype=np.float64, count=ny)
        z = np.fromfile(f, dtype=np.float64, count=nz)
        W = np.fromfile(f, dtype=np.float64, count=nx * ny * nz)
    return W, nz, ny, nx, z, y

注意：返回顺序是 W, nz, ny, nx, z, y，与输入顺序不同——这是为了与 MATLAB 版本的接口兼容。

Vect2Matr3D(V, nx, ny, nz) — 向量转 3D 矩阵

1
2
3
4
5
6
7
8
9
10
def Vect2Matr3D(V, nx, ny, nz):
    A = np.zeros((nz, ny, nx))    # 注意：形状是 (nz, ny, nx)
    nynz = ny * nz
    
    for iz in range(nz):
        for iy in range(ny):
            for ix in range(nx):
                ind = ix * nynz + iy * nz + iz   # Fortran 列优先索引
                A[iz, iy, ix] = V[ind]
    return A

索引公式图解：

1
2
3
4
对于 3D 数组 A[iz, iy, ix]（对应 z, y, x 方向）：
- ix 变化最慢（最外层）
- iz 变化最快（最内层）
- 展平索引 = ix × (ny × nz) + iy × nz + iz

Matr2Vect3D(A, nx, ny, nz) — 3D 矩阵转向量

Vect2Matr3D 的逆操作，将处理后的 3D 矩阵还原为 1D 向量以便写回二进制文件。

WriteDataB(fname, W, nz, ny, nx, z, y) — 写入二进制文件

1
2
3
4
5
6
7
def WriteDataB(fname, W, nz, ny, nx, z, y):
    with open(fname, 'wb') as fid:
        fid.write(struct.pack('3i', nx, ny, nz))        # 头部：三个 int32
        fid.write(struct.pack(f'{ny}d', *y))             # y 坐标数组
        fid.write(struct.pack(f'{nz}d', *z))             # z 坐标数组
        W_flat = W.flatten()
        fid.write(struct.pack(f'{len(W_flat)}d', *W_flat))  # 数据

LongLoss3d(h, w, sigma) — 计算纵向损失因子

这是尾场分析中最核心的物理计算：

1
2
3
4
5
6
7
8
9
def LongLoss3d(h, w, sigma):
    n = len(w)
    x = np.array([-5 * sigma + h * i for i in range(1, n+1)])
    bi2 = gauss(x, sigma)    # 高斯权重 λ(s)
    
    loss = -np.sum(bi2 * w * h)                          # 损失因子
    spread = np.sqrt(np.sum(bi2 * (w + loss) ** 2 * h))  # 能散
    
    return loss, spread

物理公式：

$$\kappa_\parallel = -\sum_i \lambda(s_i) \cdot W(s_i) \cdot \Delta s$$

$$\sigma_\kappa = \sqrt{\sum_i \lambda(s_i) \cdot (W(s_i) + \kappa_\parallel)^2 \cdot \Delta s}$$

其中 $\lambda(s)$ 是归一化的高斯分布：

$$\lambda(s) = \frac{1}{\sigma\sqrt{2\pi}} \exp\left(-\frac{s^2}{2\sigma^2}\right)$$

积分区间取 $[-5\sigma, +5\sigma]$（如果 $n$ 足够大），覆盖了高斯分布 99.9999% 的面积。

IntegrTr(h, x) — 梯形累积积分

1
2
3
4
5
6
7
8
def IntegrTr(h, x):
    n = len(x)
    y = np.zeros(n)
    y[0] = 0
    for j in range(1, n):
        y[j] = y[j-1] + 0.5 * (x[j] + x[j-1])
    y = y * h
    return y

这是从间接尾场计算横向尾场的关键步骤。根据 Panofsky-Wenzel 定理：

$$W_\perp(s) \propto \frac{\partial}{\partial r_\perp} \int_{-\infty}^{s} W_\parallel(s’) , ds’$$

IntegrTr 实现的就是这个纵向积分 $\int W_\parallel(s’) ds’$。

ReadFieldMonitor_xtime / ReadFieldMonitor_stime — 读取场监视器

两个函数分别处理两种扫描模式：

  • **_xtime**：固定空间位置，沿时间扫描（time='x'）
  • **_stime**：固定时间，沿空间扫描（time='s'）

解析监视器文件头部的元数据（k_ct, h_ct, k_y, h_y, k_z, h_z, k_s, h_s 等），然后读取数值数据矩阵。

后处理脚本详解

PP_LongitudinalWake.py — 纵向尾场分析

功能：从间接尾场提取纵向尾场势 $W_\parallel(s)$ 和损失因子分布。

核心流程：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# 1. 读取输入文件和尾场数据
nx, ny, nz, hx, hy, hz, sigma, iy0, iz0, ... = ReadInput(input_file)
W0, nz1, ny1, length1, z1, y1 = ReadDataB(mydir + in_wake_file)
FW = Vect2Matr3D(W0, length1, ny1, nz1)   # 还原为 3D

# 2. 提取束团位置的纵向切片
br = FW[zpos, ypos, :]                     # 形状 (nx,)

# 3. 计算损失因子
LossL = LongLoss3d(h, br, sigma)

# 4. 计算所有横向位置的损失因子分布
LossM = np.zeros((ny, nz))
for i in range(ny):
    for j in range(nz):
        br = FW[j, i, :]
        LossM[i, j] = LongLoss3d(h, br, sigma)[0]

输出：

  • wakeL.txt：两列数据 [s, W_parallel(s)]
  • 3D 表面图：损失因子在横截面上的分布
  • 2D 曲线图：束团中心位置的纵向尾场势

8.2 PP_CreateTransverseWake.py — 生成横向尾场

功能：从间接尾场通过数值微分生成 $W_y$ 和 $W_z$。

核心算法：

1
2
3
4
5
6
7
8
9
# 步骤1：对纵向做累积积分
WI = np.zeros((nz, ny, nx))
for j in range(ny):
    for k in range(nz):
        WI[k, j, :] = IntegrTr(hx, W[k, j, :])

# 步骤2：对横向做中心差分
Wy[:, :ny-1, :] = (WI[:, :ny-1, :] - WI[:, 1:ny, :]) / hy   # ∂/∂y
Wz[:nz-1, :, :] = (WI[:nz-1, :, :] - WI[1:nz, :, :]) / hz   # ∂/∂z

物理对应：

$$W_y \approx -\frac{1}{h_y}\frac{\partial}{\partial y}\int W_\parallel ds$$

$$W_z \approx -\frac{1}{h_z}\frac{\partial}{\partial z}\int W_\parallel ds$$

输出：

  • wake3DindirectY.bin：Y 方向横向尾场
  • wake3DindirectZ.bin：Z 方向横向尾场

PP_TransDipoleWakeY.py — 偶极尾场

功能：从两个不同束团偏移位置的横向尾场计算偶极尾场。

原理：偶极尾场是横向尾场对束团偏移的导数：

$$W_y^{\text{dipole}} = \frac{W_y(y_1) - W_y(y_0)}{y_1 - y_0}$$

1
2
3
4
5
6
7
8
# 读取两个偏移位置的横向尾场
wake1 = W[iz0-1, iy0-1, :]    # 位置 y0 的尾场
wake2 = W[iz0-1, iy1-1, :]    # 位置 y1 的尾场

# 有限差分
dy = hy * (i11 - i00)
wake2[:, 1] = (wake2 - wake1) / dy
wake2[:, 1] *= 1e3   # 转换为 1/m 单位

输出：wakeDy.dat — 偶极尾场势 $W_y^{\text{dipole}}(s)$ [V/pC/m]

PP_TransMonololeWake.py — 单极横向尾场

功能：提取束团中心位置的横向尾场（不做差分）。

1
2
# 取相邻两个网格点的平均值（因为束团可能在网格点之间）
wake[:, 1] = 0.5 * (W[iz0-1, iy0-1, :] + W[iz0-1, iy0-1 - 1, :])

输出：wakeMy.dat, wakeMz.dat — 单极横向尾场 [V/pC]

PP_TransQuadWake.py — 四极尾场

功能：计算四极尾场，即横向尾场对横向坐标的二阶导数。

$$W_y^{\text{quad}} = \frac{W_y(y_0+\Delta y) - W_y(y_0-\Delta y)}{2\Delta y}$$

1
2
3
4
5
if iy0 > 1:
    wake[:, 1] = (W[iz0-1, iy0-1, :] - W[iz0-1, iy0-1+di-2, :]) / (hy * di)
else:
    wake[:, 1] = W[iz0-1, iy0-1+di-1, :] / (hy * (di - 0.5))
wake[:, 1] *= 1e3   # 转换为 1/m 单位

输出：wakeQy.dat, wakeQz.dat — 四极尾场势 [V/pC/m]

PP_FieldMonitor.py — 场监视器动画

功能：读取场监视器文件，生成场分布的实时 3D 动画，并提取特定点的场随时间变化。

关键技术：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
from scipy.interpolate import RegularGridInterpolator

# 对每个时间步
for i in range(kt):
    FF = np.reshape(F[i, :ky * kx], (ky, kx))
    
    # 2D 插值提取任意点的场值
    interp_func = RegularGridInterpolator(
        (Y, MeshPos[i] + X), FF, 
        method='linear', fill_value=0, bounds_error=False
    )
    PM_Field[i] = interp_func([PM_y, PM_x])
    
    # 实时更新 3D 表面图
    ax.plot_surface(X_mesh, Y_mesh, FF, cmap='viridis')
    plt.pause(0.01)

Monitor_InitField.py — 初始场可视化

功能：读取初始电场，与自由空间解析解对比验证。

自由空间中，线电荷的电场为 $E_y \propto 1/y$。脚本将数值解与 $20/y$ 的解析曲线对比，验证初始场设置的正确性。

8.8 Monitor_Mesh.py — 网格材料可视化

功能：读取 Eps.mat 和 Mue.mat，用键盘交互逐层浏览 3D 材料分布。

1
2
3
4
5
6
def on_key(event):
    global index
    index = (index + 1) % nx   # 按任意键切换到下一层
    update_plots()

fig.canvas.mpl_connect('key_press_event', on_key)

对比验证脚本

  • **Ref_CompareLongitudinalWake.py**：对比 ECHO3D（3D）与 ECHO2D（2D）的纵向尾场
  • **Ref_CompareDipoleWake.py**：对比偶极尾场与参考解
  • **Ref_CompareLongitudinalWake3D.py**：对比不同网格精度的结果（网格收敛性验证）
  • **Ref_CompareMonitors.py**：对比场监视器的 2D/3D 结果

完整工作流示例

示例 1：最小纵向尾场算例（轴上束团）

目录：ECHO3D_Minimal/

步骤：

1
2
3
4
5
6
7
8
9
10
11
# 1. 运行 ECHO3D 四步流水线
cmd /c "cd /d ECHO3D_Minimal && Mesher.exe input.txt && InitField.exe input.txt && ECHO3D.exe input.txt && IndirectIntegration.exe input.txt"

# 2. 解析二进制尾场
python ECHO3D_Minimal/parse_wake.py

# 3. 绘制尾场曲线
python ECHO3D_Minimal/plot_wake.py

# 4. 分析纵向阻抗和损失因子
python ECHO3D_Minimal/analyze_longitudinal.py

关键配置：

  • BunchPosition = [0, 0]：束团在轴上
  • BunchSigma = 1：1 mm RMS 长度

预期结果：

  • Loss factor ≈ -14.6 V/pC
  • $W_\parallel$ 在 $s=0$ 处达到最小值（减速最大）

示例 2：偶极尾场算例（偏轴束团）

目录：ECHO3D_Dipole/

步骤：

1
2
3
4
5
6
7
8
9
10
11
# 1. 运行 ECHO3D
cmd /c "cd /d ECHO3D_Dipole && Mesher.exe input.txt && InitField.exe input.txt && ECHO3D.exe input.txt && IndirectIntegration.exe input.txt"

# 2. 解析尾场
python ECHO3D_Dipole/parse_wake.py

# 3. 分析纵向
python ECHO3D_Dipole/analyze_longitudinal.py

# 4. 分析横向（含 kick factor 估算）
python ECHO3D_Dipole/analyze_transverse.py

关键配置：

  • BunchPosition = [1, 0]：束团在 Y 方向偏移 1 个网格步长
  • WakeMonitor = [-1 1 -1 1]：在束团周围布置尾场监视器

预期结果：

  • Loss factor ≈ -29.3 V/pC（偏轴束团损失更大）
  • Kick factor ≈ 0.47 V/pC/m

示例 3：使用官方后处理工具链

目录：ECHO3D/（官方 Python 脚本）

1
2
3
4
5
6
7
8
9
10
11
# 1. 生成横向尾场
python ECHO3D/PP_CreateTransverseWake.py

# 2. 提取偶极尾场
python ECHO3D/PP_TransDipoleWakeY.py

# 3. 提取四极尾场
python ECHO3D/PP_TransQuadWake.py

# 4. 与参考解对比
python ECHO3D/Ref_CompareDipoleWake.py

总结与进阶方向

本文涵盖的内容

层次 内容
物理 尾场势、损失因子、踢因子、阻抗、Panofsky-Wenzel 定理
建模 ECHO3D 四步流水线、输入文件配置、几何定义
代码 pyLib4ECHO.py 全部 10 个函数、8 个后处理脚本、自定义分析工具

关键公式速查

公式 表达式
损失因子 $\kappa_\parallel = -\int \lambda(s) W_\parallel(s) ds$
踢因子 $\kappa_\perp = \int \lambda(s) W_\perp(s) ds$
纵向阻抗 $Z_\parallel(\omega) = \int W_\parallel(s) e^{-i\omega s/c} ds/c$
Panofsky-Wenzel $\partial W_\parallel / \partial r_\perp = \partial W_\perp / \partial s$
间接尾场→横向尾场 $W_y \approx -\frac{1}{h_y} \frac{\partial}{\partial y} \int W_\parallel ds$

进阶方向

  1. 阻抗计算：对尾场做 FFT 得到阻抗谱 $Z(\omega)$，分析阻抗峰对应的本征模式
  2. 网格收敛性：用不同 Steps 运行同一算例，验证结果收敛
  3. PML 边界：对于开放结构（如自由空间辐射），使用 PML 吸收边界
  4. 介电材料：在 geometry.txt 中定义 $\varepsilon_r \neq 1$ 的材料
  5. 束团分布：修改高斯分布为其他分布（如抛物线分布）
  6. 多束团仿真：研究束团间的 long-range wakefield 效应

文件结构总览

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
ECHO3D/
├── pyLib4ECHO.py              ← 核心函数库（所有脚本的依赖）
├── PP_LongitudinalWake.py     ← 纵向尾场分析
├── PP_CreateTransverseWake.py ← 生成横向尾场
├── PP_TransDipoleWakeY.py     ← 偶极尾场提取
├── PP_TransMonololeWake.py    ← 单极横向尾场
├── PP_TransQuadWake.py        ← 四极尾场提取
├── PP_FieldMonitor.py         ← 场监视器动画
├── Monitor_InitField.py       ← 初始场可视化
├── Monitor_Mesh.py            ← 网格材料可视化
├── CompareFields.py           ← 场对比
├── Ref_CompareLongitudinalWake.py  ← 纵向尾场对比验证
├── Ref_CompareDipoleWake.py        ← 偶极尾场对比验证
├── Ref_CompareLongitudinalWake3D.py ← 3D 网格收敛性验证
└── Ref_CompareMonitors.py          ← 监视器对比验证

致谢：本文基于 ECHO3D 官方发布的 Python 后处理工具链编写。ECHO3D 由 DESY 等机构开发，是加速器物理领域广泛使用的开源尾场仿真工具。