二维轴对称模型中的电磁散射

具有连续旋转对称性的物体的电磁散射模拟可以在二维轴对称模型中进行，而不需要建立三维模型，从而大大减少了计算时间和资源。在这篇博文中，我们讨论了如何利用柱坐标中的平面波展开，在二维轴对称模型中模拟任意的平面波入射。

电磁模拟中的散射场公式

当传播的电磁波撞击到一个物体时，与物体的相互作用通常会改变原始波的传播特性。这种散射事件在射频、微波和光学工程中往往会引起极大的兴趣。为了准确地模拟电磁散射，COMSOL Multiphysics® 软件的附加模块——RF 模块和波动光学模块提供了一个 散射场 公式，其中总场被写为用户定义的背景场和将被模拟求解的散射场之和，背景场的常见选项包括平面波和高斯光束。还可以灵活地定义一个任意的背景场。

柱坐标中的平面波展开

在不同的背景场中，平面波是最常用的。然而，与二维和三维模型不同，我们不能简单地将二维轴对称模型中的背景场设置为沿任意方向传播的平面波，因为它破坏了旋转对称性。但是，通过接下来介绍的建模技巧，我们可以在二维轴对称模型中实现平面波背景场。

现在，让我们考虑一个由下式给出的一般平面波背景场

在不失一般性的情况下，假设 \bm{k} 在 xz 平面上， \bm{k} 与 z 轴的角度由 \theta 给出。为了简单起见，我们只考虑 p 偏振入射，但 s 偏振也可以用类似的方式实现。对于 p 偏振，\bm{E_0}=E_0 (cos \theta \bm{\hat{x}} + sin \theta \bm{\hat{z}}) ， \bm{k} = k (sin\theta \bm{\hat{x}} – cos\theta \bm{\hat{z}})，其中 k=\omega/c。因此，平面波可以被写成

其中，\phi=atan\frac{y} {x} 是方位角，并且 r=\sqrt{x^2+y^2}。这里, 为了简单起见，我们忽略了因子 e^{i\omega t} ，但是还隐含在其中。实现扩展需要两个关键成分：平面波和柱坐标中的部分波展开

以及从笛卡尔坐标到柱坐标的变换

其中，J_m 是一阶贝塞尔函数, m。结合上述方程和一些代数可得

现在，我们终于可以把平面波背景场写成一个无限的方位角模数之和，其中柱坐标模数为 m，这意味着它可以进行二维轴对称模拟。此外，正 m和负 m 只相差一个相因子，因此原则上，求和只需要从 m=0 到 m=N 进行。例如，背景场的 \bm{r} 分量可以写为

其中，\chi 是一个阶跃函数，因此对于 m=0，\chi(m)=1，对于 m>0， \chi(m)=2 。虽然求被展开到无穷大，但如果散射的大小与波长相当，那么结果将仅在几个项之后收敛。

研究长椭球体的散射

来看一个具体的例子，我们使用 COMSOL Multiphysics 来研究一个银质长椭球体在 5 THz 到 50 THz 的红外频率范围内的散射。首先，我们创建一个含半轴 a=2\ \mu m 和 b=5\ \mu m 的长椭球体。银在相关频率范围内的光学特性可以从材料库中找到。我们可以在球体的表面施加 阻抗边界条件，因为在这个频率范围内，银的导电性很高，我们仍然要考虑小的损失。在外部边界上添加完美匹配层(PMLs)，用于吸收出射的辐射。

接下来，我们定义之前在固定方位角模数 m 下得出的背景场的每个分量，并在散射公式中将它们作为背景电场使用。在 变量 子节点中，我们还定义了散射截面，它是坡印廷矢量的表面积分。(我们将在后面进一步讨论这个问题。)

背景场的每一个分量都被定义为一个变量。

使用了一个 散射场公式。阶跃函数考虑了正负 m。

最后，我们设置了频率扫描，以 0.25 THz 的增量在 5 THz ~ 50 THz 频率范围内扫描，并使用辅助扫描从 0 到 N 扫描 m 。使用频率和辅助扫描，而不是更常见的参数扫描可以提高仿真速度。

频率扫描和辅助扫描的设置

即使对单元大小进行了非常精细的设置，宽带模拟也只需大约 5 分钟就可以完成了。在后处理中，我们将使用一些技巧来可视化总散射场和总背景场，以确保它与平面波非常相似。为此，我们利用了二维镜像 数据集(详见相关模型)。我们可以在同一幅图中绘制 \phi 和 \pi-\phi 场分布。总散射场是每个展开项的散射场之和。例如，在 30 THz 及 \phi=0 处的总散射场的模可以计算为

sqrt(abs(sum(withsol('sol1', ewfd.relEz*cos(m*0), setind(m,index),setval(freq,30[THz])), index, 1, N+1))^2 + abs(sum(withsol('sol1', ewfd.relEr*cos(m*0), setind(m,index),setval(freq,30[THz])), index, 1, N+1))^2 + abs(sum(withsol('sol1', j*ewfd.relEphi*sin(m*0), setind(m,index),setval(freq,30[THz])), index, 1, N+1))^2)

这里，使用 withsol 算子挑选出频率为 30 THz ，方位角模数 index 的散射场，算子用于对每个解的贡献进行求和。

通过对每个展开项的贡献求和，在 30 THz 处的总散射场模的二维图。

另外，我们可以利用二维旋转 数据集来绘制三维散射场。请注意，在默认情况下，二维旋转 数据集只是实现了一个体旋转，其中的方位角依赖性被忽略了。因此，我们首先在二维旋转 数据集的设置中启用高级 标签下的定义变量，然后可以手动添加正确的 \phi 依赖。这样就可以为方位角启用一个名为 rev1phi 的变量。最终，三维总散射场模的正确表达就是：

sqrt(abs(sum(withsol('sol1', ewfd.relEz*cos(m*rev1phi), setind(m,index),setval(freq,30[THz])), index, 1, N+1))^2 + abs(sum(withsol('sol1', ewfd.relEr*cos(m*rev1phi), setind(m,index),setval(freq,30[THz])), index, 1, N+1))^2 + abs(sum(withsol('sol1', j*ewfd.relEphi*sin(m*rev1phi), setind(m,index),setval(freq,30[THz])), index, 1, N+1))^2)

二维旋转数据集的同一场的三维表面图。

背景场可以用与散射场相同的方式来绘制。我们可以看到，只有 N=4 展开时，散射周围的背景场才与平面波无异。

通过对每个展开项的贡献求和，得到总背景场的 z 分量。

最后，我们想将散射截面绘制为频率的函数。请注意，在变量 中，我们定义了散射截面的表达式。总散射截面可以计算为

-withsol('sol1', sigma_sc, setval(m,0),setval(freq, freq))-sum(withsol('sol1', 0.5*sigma_sc, setval(m,val),setval(freq, freq)),val,1,N)

减号是由于球状体表面的法线指向内部。 对于 m>0 的项，乘性系数 0.5 是为了补偿能量的加倍，因为每个项都是正负之和。通常，总的散射截面不是简单的每个展开项的计算之和。这是因为散射截面和坡印廷矢量直接与能量相关，而不是与场相关。能量是与场的平方成正比的。因此，当对许多项的总和进行平方时，除了平方的总和之外，还出现了交叉项。幸运的是，在我们的例子中，由于方位角模数的正交性，交叉项都消失了，也就是说，当积分到超过 \phi，交叉项都是零。因此，总的散射截面只是每个展开项的总和。我们可以看到，长椭圆球体的散射截面在质量上与球体的散射截面相似，这是著名的米氏计算的结果。

银质长椭圆球体在红外频率下的散射截面。

结论

在这篇博文中，我们介绍了一种使用二维轴对称模型模拟平面波激励下的旋转体散射特性的方法。与完整的三维模拟相比，使用这种方法进行模拟可以获得巨大的回报。计算的内存和时间成本至少要小一个数量级。因此，可以使用非常精细的网格实现非常高的模拟精度。此外，由于 RAM 的要求很小，而需要扫描的参数很多，因此使用批量扫描功能在同一时间同时运行多个模拟具有明显的优势，尽管我们在这里并没有在示例模型中演示。我们演示了散射场和散射截面的计算，但也可以得出与散射问题相关的其他数量，如远场辐射模式。最后，需要注意的是，在计算各种物理量时，必须注意追踪与模数相关的相因子。

动手尝试

点击下面的按钮，进入 COMSOL 案例库。尝试自己动手模拟平面波散射展开 模型：