Helmholtz 散射（DtN 边界处理）¶

examples/wave/helmholtz_dtn 中的两个自包含脚本分别处理全空间散射和水平周期层散射。同一标量形式同时适用于电磁 TE/TM 极化和声压波。

如果需要吸收层（PML）或端口边界，而不是精确的 DtN 映射，请参见开放域波动（PML 与端口）。

采用时间约定 \(e^{-i\omega t}\)，并将方程写成

\[\nabla\!\cdot\!\bigl(a(x)\nabla u\bigr)+\omega^2b(x)u=0.\]

这里 \(a,b\in L^\infty\) 是实值、一致正的标量函数。两个系数都可以随空间变化；下文解析示例采用的分片常数材料只是其中一种特殊情形。

表 1 标量系数的物理含义¶
模型	标量场	\(a(x)\)	\(b(x)\)
`TM_EZ`	\(u=E_z\)	\(\mu(x)^{-1}\)	\(\varepsilon(x)\)
`TE_HZ`	\(u=H_z\)	\(\varepsilon(x)^{-1}\)	\(\mu(x)\)
声压	\(u=p\)	\(\rho(x)^{-1}\)	\(K(x)^{-1}=1/[\rho(x)c(x)^2]\)

电磁模型沿用本项目约定：TM_EZ 表示面外电场，TE_HZ 表示面外磁场。示例采用归一化单位。使用相对材料参数时，TM_EZ 的系数为 \((a,b)=(\mu_r^{-1},\varepsilon_r)\)，TE_HZ 的系数为 \((a,b)=(\varepsilon_r^{-1},\mu_r)\)。在任一均匀区域 \(j\) 中，相应波数为

\[k_j=\omega\sqrt{b_j/a_j}.\]

运行示例

python examples/wave/helmholtz_dtn/circular_dtn.py
python examples/wave/helmholtz_dtn/periodic_dtn.py

全空间散射¶

一般散射问题¶

设 \(D\Subset\mathbb R^2\) 是边界为 \(\Gamma\) 的有界障碍物或可穿透夹杂。在某个有界区域之外，材料为常数 \(a_0,b_0>0\) 的均匀介质，外部波数为 \(k_0=\omega\sqrt{b_0/a_0}\)。

在选定具体入射波之前，令 \(u^{\rm inc}\) 是外部方程的任意整解

\[\Delta u^{\rm inc}+k_0^2u^{\rm inc}=0.\]

散射体产生出射场 \(u^{\rm sc}\)。在外部区域中，总场为

\[u=u^{\rm inc}+u^{\rm sc}.\]

系数本身并非紧支，因为波必须在外部介质中传播；具有紧支集的是材料反差，

\[\operatorname{supp}(a-a_0)\cup\operatorname{supp}(b-b_0) \Subset\mathbb R^2.\]

对于障碍物，总场未知量定义在 \(D\) 的外部。对于可穿透夹杂，总场定义在整个平面上，且系数可以在紧支反差区域内变化。后面的圆形验证算例进一步将夹杂内部取为常数 \((a_1,b_1)\)。

物理边界与界面条件¶

三个算例只在物理圆周上的条件不同。令 \(\nu\) 从 \(D\) 指向外部，它们分别满足

\[\begin{split}\begin{aligned} u&=0 &&\text{on }\Gamma &&\text{(Dirichlet obstacle)},\\ a_0\partial_\nu u&=0 &&\text{on }\Gamma &&\text{(Neumann obstacle)},\\ [u]&=0,\qquad [a\partial_\nu u]=0 &&\text{on }\Gamma &&\text{(penetrable inclusion)}. \end{aligned}\end{split}\]

透射算例中的跳跃定义为内侧迹减去外侧迹。因此，场和加权法向通量均连续。

全空间辐射条件¶

散射场满足出射条件。对于 \(e^{-i\omega t}\) 时间约定，它满足二维 Sommerfeld 条件

\[\lim_{r\to\infty}\sqrt r\, \left(\partial_ru^{\rm sc}-ik_0u^{\rm sc}\right)=0\]

且关于方向一致成立。该条件选取第一类 Hankel 函数；在上述材料一致正的假设下，它唯一确定出射解。

圆周 DtN 截断¶

选择人工圆周 \(\Gamma_R=\{r=R\}\)，使其包围物理散射体和全部材料反差。若

\[\widehat\phi_m=\frac{1}{2\pi}\int_0^{2\pi} \phi(\theta)e^{-im\theta}\,d\theta,\]

则精确的外部 Fourier--Hankel DtN 映射为

\[T_R\phi =\sum_{m\in\mathbb Z}\tau_m\widehat\phi_m e^{im\theta}, \qquad \tau_m=k_0\frac{H_m^{(1)\prime}(k_0R)}{H_m^{(1)}(k_0R)}.\]

由于 \(\partial_n-T_R\) 消去每个出射 Hankel 模态，无界域辐射条件等价于如下总场边界条件

\[\partial_nu-T_Ru=g_R, \qquad g_R=\partial_nu^{\rm inc}-T_Ru^{\rm inc} \quad\text{on }\Gamma_R,\]

其中 \(n=e_r\) 是截断区域的外法向。

全空间弱形式¶

对障碍物令 \(\Omega_R=B_R\setminus\overline D\)，对透射问题令 \(\Omega_R=B_R\)。Dirichlet 障碍物的试验空间为

\[V_D=\{v\in H^1(\Omega_R):v|_\Gamma=0\}\]

Neumann 和透射算例采用 \(V=H^1(\Omega_R)\)。在相应仿射试探空间中求 \(u\)，使得对其齐次对应空间内任意试验函数 \(v\)，都有

\[\int_{\Omega_R}a\nabla u\cdot\nabla\overline v -\omega^2\int_{\Omega_R}bu\overline v -a_0\langle T_Ru,v\rangle_{\Gamma_R} =a_0\langle g_R,v\rangle_{\Gamma_R}.\]

Dirichlet 迹通过 Condenser 强制施加。齐次 Neumann 条件是自然边界条件，因而不产生边界积分。对透射问题，共形网格和逐单元材料系数 \(a,b\) 通过弱形式自然满足 \([u]=[a\partial_\nu u]=0\)。

周期散射¶

周期材料与竖直方向紧支性¶

设 \(x_1\) 方向的周期为 \(L\)。系数满足

\[a(x_1+L,x_2)=a(x_1,x_2), \qquad b(x_1+L,x_2)=b(x_1,x_2).\]

非平凡周期系数不可能在整个平面内紧支。正确的开放域假设是：在一个周期胞内，材料反差沿竖直方向紧支，即在某个有界水平带的上方和下方，系数均变为常数。在本示例中，

\[\exists\,y_-<y_+:\qquad a(x_1,x_2)=a_0,\quad b(x_1,x_2)=b_0 \quad\text{for }x_2\notin[y_-,y_+],\ 0<x_1<L.\]

slab 算例在 \(0<x_2<h\) 内采用材料 \((a_1,b_1)\)，上下外部半空间均采用材料 \((a_0,b_0)\)。wall 算例只保留 \(x_2>0\)：同一有限层在 \(x_2=0\) 处由物理边界封闭，并在上方开放。

Bloch 条件¶

给定水平 Bloch 波数 \(\alpha\)，总场满足

\[u(L,x_2)=e^{i\alpha L}u(0,x_2), \qquad \partial_{x_1}u(L,x_2) =e^{i\alpha L}\partial_{x_1}u(0,x_2).\]

由于 \(a\) 具有周期性，水平方向材料通量也具有相同相位。试探和试验函数均满足这一拟周期性，因此 \(\overline v\) 带有共轭相位，弱形式中左右边界项相互抵消。

一般入射场、散射场与总场¶

向下传播的平面波从顶部参考线 \(\Gamma_+=\{x_2=y_t\}\) 入射：

\[u^{\rm inc} =Ae^{i\alpha x_1-i\beta_{0,0}(x_2-y_t)}.\]

对 \(m\in\mathbb Z\)，定义

\[\alpha_m=\alpha+\frac{2\pi m}{L}, \qquad \beta_{j,m}=\sqrt{k_j^2-\alpha_m^2}.\]

假设入射阶为传播模态，即 \(|\alpha|<k_0\)，因此 \(\beta_{0,0}>0\) 且入射功率非零。对于每个保留阶，传播模态的出射分支满足 \(\beta_{j,m}>0\)，倏逝模态满足 \(\operatorname{Im}\beta_{j,m}>0\)；掠射模态则有 \(\beta_{j,m}=0\)。示例参数避开了保留模态满足 \(\alpha_m^2=k_j^2\) 的 Wood 异常；在这类截止频率处，Rayleigh 展开和功率归一化需要单独处理。结构上方的散射场采用如下出射 Rayleigh 展开

\[u^{\rm sc} =\sum_{m\in\mathbb Z}r_m e^{i\alpha_mx_1+i\beta_{t,m}(x_2-y_t)}, \qquad x_2\geq y_t.\]

因此，顶部总场为

\[u=u^{\rm inc}+u^{\rm sc}.\]

对于双端口 slab，底部没有入射波，因此底部总场就是透射出射场

\[u=u^{\rm tr} =\sum_{m\in\mathbb Z}t_m e^{i\alpha_mx_1-i\beta_{b,m}(x_2-y_b)}, \qquad x_2\leq y_b.\]

这些 Rayleigh 展开是 Sommerfeld 辐射条件在周期问题中的对应形式。

物理界面与墙面边界条件¶

三个配置分别为：

slab：两个平面材料界面上的场和加权通量连续，即在 \(x_2=0,h\) 上满足 \([u]=0\) 和 \([a\partial_\nu u]=0\)，且上下端口均开放；
dirichlet：\(x_2=h\) 处的材料界面满足相同的两个透射条件，墙面 \(x_2=0\) 上满足 \(u=0\)；
neumann：\(x_2=h\) 处满足界面条件，\(x_2=0\) 上满足 \(a\partial_\nu u=0\)。

Rayleigh DtN 截断¶

令 \(\Gamma_-\) 表示 slab 的底部参考线。以拟周期基 \(e^{i\alpha_mx_1}\) 表示模态系数，定义

\[T_\pm\phi =\sum_{m\in\mathbb Z}i\beta_{\pm,m}\widehat\phi_m e^{i\alpha_mx_1}.\]

两个算子的符号都是 \(i\beta_{\pm,m}\)，因为它们分别使用截断周期胞的外法向：顶部为 \(n=+e_2\)，底部为 \(n=-e_2\)。顶部总场条件为

\[\partial_nu-T_+u=g_+, \qquad g_+=-2i\beta_{0,0}Ae^{i\alpha x_1} \quad\text{on }\Gamma_+.\]

对 slab，底部场只含出射波，因此

\[\partial_nu-T_-u=0 \quad\text{on }\Gamma_-.\]

周期弱形式¶

令 \(\Omega_L\) 为一个截断周期胞，并定义

\[V_\alpha=\left\{v\in H^1(\Omega_L): v(L,x_2)=e^{i\alpha L}v(0,x_2)\right\}.\]

Dirichlet 墙面采用 \(V_{\alpha,D}=\{v\in V_\alpha:v=0\text{ on }x_2=0\}\)；slab 和 Neumann 墙面采用 \(V_\alpha\)。对 slab 令 \(\delta_{\rm slab}=1\)，对 wall 令其为零。在相应试探空间中求 \(u\)，使得

\[\int_{\Omega_L}a\nabla u\cdot\nabla\overline v -\omega^2\int_{\Omega_L}bu\overline v -a_t\langle T_+u,v\rangle_{\Gamma_+} -\delta_{\rm slab}a_b\langle T_-u,v\rangle_{\Gamma_-} =a_t\langle g_+,v\rangle_{\Gamma_+}\]

对任意齐次试验函数 \(v\) 均成立。当前 slab 中 \(a_t=a_b=a_0\)。Neumann 墙面是自然边界，不产生底部边界项。该弱形式确定实现顺序：先装配体积和 DtN 算子，再用 BlochReducer 形成 \(T^HAT\) 和 \(T^Hb\)，最后只对 Dirichlet 墙面应用 Condenser。

用于验证的分层精确解¶

一般周期问题可能耦合多个 Rayleigh 阶。为了验证代码，平界面保持入射模态的阶数不变。令 \(\beta=\beta_{0,0}\)、\(\gamma=\beta_{1,0}\)，则入射场、散射场和顶部总场为

\[\begin{split}\begin{aligned} u^{\rm inc}&=Ae^{i\alpha x_1-i\beta(x_2-y_t)},\\ u^{\rm sc}&=re^{i\alpha x_1+i\beta(x_2-y_t)},\\ u_{\rm top}&=u^{\rm inc}+u^{\rm sc}. \end{aligned}\end{split}\]

介质层内总场和底部透射总场为

\[\begin{split}\begin{aligned} u_{\rm layer}&=e^{i\alpha x_1} \left[pe^{-i\gamma x_2}+qe^{i\gamma x_2}\right],\\ u_{\rm bottom}&=u^{\rm tr} =te^{i\alpha x_1-i\beta(x_2-y_b)}. \end{aligned}\end{split}\]

在 \(x_2=h\) 和 \(x_2=0\) 处令 \(u\) 与 \(a\partial_{x_2}u\) 连续，可得到脚本中关于 \((r,p,q,t)\) 的四阶线性系统。wall 算例移除底部端口，并改用介质层内总场

\[u_D=d\,e^{i\alpha x_1}\sin(\gamma x_2), \qquad u_N=d\,e^{i\alpha x_1}\cos(\gamma x_2).\]

正弦函数满足 Dirichlet 迹，余弦函数在 \(x_2=0\) 处法向导数为零；在 \(h\) 处匹配可得到关于 \((r,d)\) 的二阶线性系统。

周期求解使用 TensorMesh 的公共约束算子：reduce() 构造 \(T^H A T\)，reduce_rhs() 构造 \(T^H b\)，reduce_mask() 则把墙面自由度映射到约化后的自由度上，供 Condenser 使用：

列表 14 examples/wave/helmholtz_dtn/periodic_dtn.py（约束处理核心代码）¶

reducer = BlochReducer(mesh.points, [[period, 0.0]], sign=+1)
wavevector = torch.tensor([alpha, 0.0], dtype=torch.float64)
reduced_matrix = reducer.reduce(matrix, wavevector)       # T^H A T
reduced_rhs = reducer.reduce_rhs(rhs, wavevector)         # T^H b

if case == "dirichlet":
    master_mask = reducer.reduce_mask(bottom_mask)
    zero_trace = torch.zeros(int(master_mask.sum()), dtype=torch.complex128)
    condenser = Condenser(master_mask, dirichlet_value=zero_trace)
    inner_matrix, inner_rhs = condenser(reduced_matrix, reduced_rhs)
    inner_solution = inner_matrix.solve(inner_rhs, backend="scipy", method="lu")
    reduced_solution = condenser.recover(inner_solution)
else:                                                     # natural Neumann wall
    reduced_solution = reduced_matrix.solve(
        reduced_rhs, backend="scipy", method="lu"
    )

solution = reducer.recover(reduced_solution, wavevector)

Dirichlet 背衬介质层的数值场、精确场与误差场。 — 图 57 Dirichlet 背衬层。¶

Neumann 背衬介质层的数值场、精确场与误差场。 — 图 58 Neumann 背衬层。¶

小写的 \(r_m,t_m\) 表示复振幅。反射率与透射率是功率比例：

\[\mathcal R_m= \frac{\operatorname{Re}(a_t\beta_{t,m})} {\operatorname{Re}(a_t\beta_{t,0})} \left|\frac{r_m}{A}\right|^2, \qquad \mathcal T_m= \frac{\operatorname{Re}(a_b\beta_{b,m})} {\operatorname{Re}(a_t\beta_{t,0})} \left|\frac{t_m}{A}\right|^2.\]

壁面没有下方辐射端口，因此脚本只报告 \(r\) 和 \(\mathcal R\)，不会人为给出零透射量。

TensorMesh 装配细节

阅读顺序与官方 Helmholtz 示例一致：

网格几何、基函数、Jacobian 和求积均保持为 float64；
将 HelmholtzAssembler 提升为 complex128，并用 ElementAssembler 装配 \(a\nabla u\cdot\nabla v-\omega^2buv\)；
用 FacetAssembler 计算模态迹；
将非局部 DtN 块加入 TensorMesh 稀疏 COO 矩阵；
施加约束，求解一般复数系统，并与真解比较。

加权体积变分形式是唯一针对本问题定义的单元装配器。基函数计算、数值求积、Jacobian 变换、单元索引和全局稀疏散布均由 TensorMesh 完成：

列表 15 加权复 Helmholtz 装配（共享模式）¶

class HelmholtzAssembler(ElementAssembler):
    def forward(self, gradu, gradv, u, v, a, b, omega_sq):
        return a * (gradu @ gradv) - omega_sq * b * u * v

assembler = HelmholtzAssembler.from_mesh(mesh, quadrature_order=4)
assembler.type(torch.complex128)
volume = assembler(
    points=mesh.points,
    element_data={
        "a": a.to(torch.complex128),
        "b": b.to(torch.complex128),
        "omega_sq": torch.full_like(a, omega**2, dtype=torch.complex128),
    },
)

TensorMesh 0.2.0 还提供 FacetBilinearAssembler 来装配局部 Robin 或阻抗矩阵。但 DtN 映射是非局部的，因此这些示例有意使用 FacetAssembler 计算模态矩，并在积分完成后再形成低秩边界块。

面装配路径为实数，因此每个指数模态先按余弦和负正弦两个通道积分，再组合为 complex128。若 \(\Psi\) 为模态矩矩阵，则两个示例都使用

Psi.conj() @ torch.diag(symbols) @ Psi.T

最后一个因子采用普通转置，因为 \(\Psi^Tu\) 给出模态系数。这里不假设矩阵为 Hermitian 或正定。脚本调用 TensorMesh 的 SparseMatrix.solve，并显式选择 SciPy LU 后端，而不是请求 SPD 求解器。

对于带周期 Dirichlet 壁面的算例，操作顺序至关重要：

装配完整的体积矩阵和 DtN 矩阵；
使用 BlochReducer(sign=+1) 形成 \(T^HAT\) 和 \(T^Hb\)；
将完整底部迹映射到唯一的 Bloch 主自由度；
应用 Condenser 并求解；
先通过 Condenser 恢复，再通过 BlochReducer 恢复。

Neumann 壁面是自然弱边界，因此跳过凝聚。测试将 Bloch 缩减与显式矩阵 \(T\) 比较，把阶数 \(-2,0,3\) 代回圆形边界条件，检验 P1 收敛性，并验证振幅、功率、Bloch 迹和非零模态泄漏。

绘图与解析验证

两个脚本均使用 plot()，没有另行搭建绘图流程：

列表 16 数值场、真解场与逐点误差¶

numerical = solution.detach().cpu()
reference = torch.as_tensor(exact, dtype=torch.complex128)
fields = {
    "Re(u_h)": numerical.real,
    "Re(u_exact)": reference.real,
    "error": abs(numerical - reference),
}
mesh.plot(fields, save_path=filename, show_mesh=False)

在默认 CPU complex128 设置下，三个圆形算例的相对节点误差均小于 \(5.6\times10^{-2}\)，拟合得到的 P1 收敛阶位于 \(1.71\) 到 \(1.74\) 之间。对于周期算例，解析振幅的相对误差小于 \(2.0\times10^{-2}\)，能量平衡残差小于 \(2.0\times10^{-15}\)，恢复后的 Bloch 迹在机器精度范围内一致。提交的图片展示了每个算例的数值场、解析场及其逐点误差。