您的位置: 首页  >  文章  >  电磁波在柱坐标和球坐标下的本征解

电磁波在柱坐标和球坐标下的本征解

分类: 文章 • 2024-11-01 22:24:10

一、背景

课程和习题中，我们通常接触的都是平面电磁波，但现实生活中却常常碰到柱面波和球面波，如通电导线辐射场、手机信号等。而且工程上，平面波也可以按柱面波和球面波展开。

那么电磁波在柱坐标和球坐标下的本征解是什么样的形式呢？什么样的波源可以产生这样的电磁场呢？这就是本文要讨论的问题。

二、真空中的时谐场

首先列出真空中的麦克斯韦方程：
{∇⋅E⃗=0∇×E⃗=−∂B⃗∂t∇⋅B⃗=0∇×B⃗=1c2∂E⃗∂t \begin{cases} \nabla \cdot \vec{E}=0 \\ \nabla \times \vec{E}=-\frac{\partial \vec{B}}{\partial t} \\ \nabla \cdot \vec{B}=0 \\ \nabla \times \vec{B} =\frac{1}{c^2} \frac{\partial \vec{E}}{\partial t} \\ \end{cases}
分别表示空间中无电荷，法拉第电磁感应定律，无磁单极子，安培定律。第四项对时间求导可得：
1c2∂2E⃗∂t2=∇×∂B⃗∂t=∇×(−∇×E⃗)=−[∇(∇⋅E⃗)−∇2E⃗]=∇2E⃗ \frac{1}{c^2} \frac{\partial ^2\vec{E}}{\partial t^2} = \nabla \times \frac{\partial \vec{B}}{\partial t} = \nabla \times (-\nabla \times \vec{E})=-[\nabla(\nabla \cdot \vec{E})-\nabla ^2 \vec{E}]=\nabla ^2 \vec{E}
第三步用到了矢量微分的运算公式，可以用直角坐标展开来证明，过程比较繁琐。

在考虑时谐场，即E⃗=E⃗(x,y,z)e−iωt\vec{E}=\vec{E}(x,y,z)e^{-i \omega t}，带入上式可得到：
∇2E⃗+k2E⃗=0 \nabla ^2 \vec{E}+k^2\vec{E}=0
k≡ω/ck \equiv \omega /c这就得到了亥姆霍兹方程。可以容易得到直角坐标系下本征解：
E⃗=C⋅ei(ky−ωt)z^ \vec{E}=C \cdot e^{i(ky-\omega t)} \hat{z}
电磁波在柱坐标和球坐标下的本征解

三、柱坐标下的本征解

亥姆霍兹方程求解的一个难点是，Laplace算符作用的是带有方向的矢量，直角坐标系下，基矢方向不变，但柱坐标和球坐标下，基矢也会随求导改变。
电磁波在柱坐标和球坐标下的本征解

为了计算简单，假设波矢方向k⃗\vec{k}沿r^\hat{r}方向，电场方向沿zz轴，大小只与r,ϕr,\phi有关，与zz无关。即E⃗=E(r,ϕ)z^\vec{E}=E(r,\phi)\hat{z}，参考附录中的矢量微分公式可得：
∇2E⃗=(∇2E⃗)z=∇2Ezz^=∇2E(r,ϕ)z^=[1r∂∂r(r∂E∂r)+1r2∂2E∂ϕ2+∂2E∂z2]z^=(∂2E∂r2+1r∂E∂r+1r2∂2E∂ϕ2)z^ \nabla ^2 \vec{E}=(\nabla ^2 \vec{E})_z=\nabla ^2 E_z \hat{z}=\nabla ^2 E(r,\phi) \hat{z}=[\frac{1}{r} \frac{\partial}{\partial r}(r\frac{\partial E}{\partial r})+\frac{1}{r^2} \frac{\partial^2 E}{\partial \phi^2}+\frac{\partial^2 E}{\partial z^2}]\hat{z}\\ =(\frac{\partial^2 E}{\partial r^2}+\frac{1}{r}\frac{\partial E}{\partial r}+\frac{1}{r^2} \frac{\partial^2 E}{\partial \phi^2})\hat{z}
带入亥姆霍兹方程，消去方向z^\hat{z}得：
∂2E∂r2+1r∂E∂r+1r2∂2E∂ϕ2+k2E=0 \frac{\partial^2 E}{\partial r^2}+\frac{1}{r}\frac{\partial E}{\partial r}+\frac{1}{r^2} \frac{\partial^2 E}{\partial \phi^2}+k^2E=0
这个方程需要用分离变量法求解：E=R⋅ΨE=R\cdot\Psi
{d2Ψdϕ2+m2Ψ=0d2Rdr2+1rdRdr+(k2−m2r2)R=0 \begin{cases} \frac{ d^2\Psi}{d\phi^2}+m^2 \Psi=0 \\ \frac{d^2 R}{dr^2}+\frac{1}{r}\frac{d R}{d r}+(k^2-\frac{m^2}{r^2})R=0 \\ \end{cases}
解得：
E=Hm(kr)[Acos⁡(mϕ)+Bsin⁡(mϕ)] E=H_m(kr)[A\cos(m\phi)+B\sin(m\phi)]
这里R满足的方程是贝塞尔方程，不同的m可以解出不同的表达式，一般的求解方法是多项式待定系数法，解起来非常繁琐，姑且用一个整体表示它即Hm(kr)H_m(kr)称为汉克函数，它包含实部和虚部：
Hm(kr)=Jm(kr)+iNm(kr) H_m(kr)=J_m(kr)+iN_m(kr)
Jm(kr)J_m(kr)为贝塞尔函数，Nm(kr)N_m(kr)为诺伊曼函数，定性上只要知道这两个函数在方向上震荡衰减，衰减的振幅正比于1r\frac {1}{\sqrt{r}}，即E⃗\vec{E}在较远处(大于一个波长就有相当好的近似)可简化为：
E⃗=[Acos⁡(mϕ)+Bsin⁡(mϕ)]Hm(kr)e−iωtz^≈A′sin⁡(mϕ+B′)1krei(kr−ωt)z^ \vec{E}=[A\cos(m\phi)+B\sin(m\phi)]H_m(kr)e^{-i\omega t}\hat{z}\approx A'\sin(m\phi+B')\frac{1}{\sqrt{kr}}e^{i(kr-\omega t)}\hat{z}
我们可以看出，电场的大小绕着z轴(沿ϕ\phi方向)周期震荡，沿着r方向呈现1r\frac {1}{\sqrt{r}} 递减，相位关系与平面波类似ei(kr−ωt)e^{i(kr-\omega t)}。电场大小随1r\frac {1}{\sqrt{r}} 递减也反映了能流守恒，因为单位高度的同轴椭圆的侧面积S正比于r，单位时间流出单位侧面积的能量正比于E2∝1/rE^2\varpropto 1/r，因此单位高度的同轴椭圆总能量流出守恒。

四、球坐标下的本征解

电磁波在柱坐标和球坐标下的本征解
为了简化问题，假设波矢k⃗\vec{k}的方向沿r^\hat{r}，电场的大小只与r,θr,\theta有关，而与ψ\psi无关，即E⃗=E(r,θ)θ^\vec{E}=E(r,\theta)\hat{\theta}，参考附录中的矢量微分公式可得：
∇2E⃗=(∇2E⃗)θ=∇2Eθθ^+2r2[0−Eθ2sin⁡2θ+0]θ^=[∂2E∂r2+2r∂E∂r+1r2sin⁡θ∂∂θ(sin⁡θ∂E∂θ)−1r2sin⁡2θE]θ^ \nabla ^2 \vec{E}=(\nabla ^2 \vec{E})_{\theta}=\nabla ^2 E_{\theta} \hat{\theta}+\frac{2}{r^2}[0-\frac{E_{\theta}}{2\sin^2\theta}+0]\hat{\theta}\\ =[\frac{\partial^2 E}{\partial r^2}+\frac{2}{r}\frac{\partial E}{\partial r}+\frac{1}{r^2 \sin \theta} \frac{\partial}{\partial \theta}(\sin \theta \frac{ \partial E}{\partial \theta})-\frac{1}{r^2\sin^2 \theta}E]\hat{\theta}
带入亥姆霍兹方程，消去方向θ^\hat{\theta}得：
∂2E∂r2+2r∂E∂r+1r2sin⁡θ∂∂θ(sin⁡θ∂E∂θ)−1r2sin⁡2θE+k2E=0 \frac{\partial^2 E}{\partial r^2}+\frac{2}{r}\frac{\partial E}{\partial r}+\frac{1}{r^2 \sin \theta} \frac{\partial}{\partial \theta}(\sin \theta \frac{ \partial E}{\partial \theta})-\frac{1}{r^2\sin^2 \theta}E+k^2E=0
用分离变量法求解：E=R⋅ΘE=R\cdot\Theta
{1sin⁡θddθ(sin⁡θdΘdθ)−1sin⁡2θΘ+l(l+1)Θ=0d2Rdr2+2rdRdr+[k2−l(l+1)r2]R=0 \begin{cases} \frac{1}{ \sin \theta} \frac{d }{d \theta}(\sin \theta \frac{ d \Theta}{d \theta})-\frac{1}{\sin^2 \theta}\Theta+l(l+1)\Theta=0 \\ \frac{d^2 R}{dr^2}+\frac{2}{r}\frac{d R}{d r}+[k^2-\frac{l(l+1)}{r^2}]R=0 \\ \end{cases}
解得：
E=C⋅hl(kr)Pl(cos⁡θ) E=C\cdot h_l(kr)P_l(\cos \theta)
其中CC为常数，Pl(cos⁡θ)P_l(\cos \theta)是Θ\Theta所满足方程的解，它是ll阶m=1m=1的关联勒让德函数，这个函数在解氢原子波函数的时候也会用到。径向的方程RR满足球贝塞尔方程，其解hl(kr)h_l(kr)是球汉克函数：
hl(kr)=jl(kr)+inl(kr) h_l(kr)=j_l(kr)+in_l(kr)
jl(kr)j_l(kr)表示ll阶的球贝塞尔函数，nl(kr)n_l(kr)表示ll阶的球诺伊曼函数，与柱坐标系的结果类似，我们只需要知道这个表达式在较远处的行为：
E⃗=C⋅Pl(cos⁡θ)hl(kr)e−iωtθ^≈C⋅Pl(cos⁡θ)1krei(kr−ωt)θ^∝1rei(kr−ωt)θ^ \vec{E}= C\cdot P_l(\cos \theta)h_l(kr)e^{-i\omega t}\hat{\theta} \approx C \cdot P_l(\cos \theta)\frac{1}{kr}e^{i(kr-\omega t)}\hat{\theta}\varpropto \frac{1}{r}e^{i(kr-\omega t)}\hat{\theta}
因此可以定性看出，电场的方向沿θ^\hat{\theta}，幅度与ψ\psi无关，随θ\theta变化，随1/r1/r递减，电磁波整体的方向沿r^\hat{r}传播。

振幅呈1/r1/r衰减也反应了能量守恒，因为球面的面积正比于r2r^2，球面上单位面积流出的能量正比于E2∝1/r2E^2 \varpropto 1/r^2，因此半径为rr的球面流出的能量是守恒的。

Tip:

求解亥姆霍兹方程时要利用拉普拉斯算符∇2\nabla ^2，与作用于不带方向的标量(如氢原子波函数)相比，当它作用于带方向的电场时，会多出几项(查看附录的矢量微分公式)，这是因为这里的“方向”也是会随坐标改变的，求偏导时会多出因方向变化导致的添加项。

五、本征解对应的波源可能是什么？

现在我们已经清楚了电场波在直角坐标、柱坐标和球坐标下的本征解，假设空间中真的有这种波存在，那么它的波源是什么呢？总不可能凭空产生电磁场吧？

我们先列出这些波在基态的形式：
平面波：E⃗=C⋅ei(ky−ωt)z^\vec{E}=C \cdot e^{i(ky-\omega t)} \hat{z}
柱面波：E⃗=H0(kr)e−iωtz^≈C′ei(kr−ωt)krz^\vec{E}=H_0(kr)e^{-i\omega t}\hat{z}\approx C'\frac{e^{i(kr-\omega t)}}{\sqrt{kr}}\hat{z}(m=0m=0)
球面波：E⃗=P0(cos⁡θ)h0(kr)e−iωtθ^≈C′sin⁡θei(kr−ωt)krθ^\vec{E}= P_0(\cos \theta)h_0(kr)e^{-i\omega t}\hat{\theta} \approx C' \sin \theta \frac{e^{i(kr-\omega t)}}{kr}\hat{\theta}(l=0l=0)

由空间对称性可以猜想，平面波是由无限大平面的震荡电流产生的；柱面波是由无限长导线的震荡电流产生的；而球面波的形式就是偶极子。下面具体推导验证这些波源产生的波是否和本征解一致。

1、无限大平面的电流源

电磁波在柱坐标和球坐标下的本征解

如图，无限大平面在oxz平面上，电流线密度α=α0e−iωt\alpha=\alpha_0 e^{-i\omega t}，现在计算(0,y,0)处的电场强度，为此需要先求磁矢势A⃗\vec{A}，再对其求旋度得到磁场强度B⃗\vec{B}，最后再求电场强度E⃗\vec{E}，oxz平面上，距离远点l=x2+z2l=\sqrt{x^2+z^2}处的一小段电流为α0dx⋅dz\alpha_0 dx \cdot dz，因为它会随时间变化，因此该点到(0,y,0)处的相位延时为ky2+l2k\sqrt{y^2+l^2}，可以得到磁矢势的积分形式：
A⃗=μ04π∫α0e−i(ωt−ky2+l2)dxdzy2+l2z^ \vec{A}=\frac{\mu_0}{4\pi}\int\frac{\alpha_0 e^{-i(\omega t-k\sqrt{y^2+l^2})}dxdz}{\sqrt{y^2+l^2}}\hat{z}
x、z的平面积分转为极坐标积分，因为小段电流的相位和大小与角度无关，可直接积分：
A⃗=μ0α02e−iωt∫0∞eiky2+l2ldly2+l2z^ \vec{A}=\frac{\mu_0\alpha_0}{2}e^{-i\omega t}\int_0^{\infty}\frac{e^{i k\sqrt{y^2+l^2}}ldl}{\sqrt{y^2+l^2}}\hat{z}
对e指数泰勒展开可得：
A⃗=μ0α02e−iωt∫0∞∑n=0∞(iky2+l2)nn!y2+l212dl2z^ \vec{A}=\frac{\mu_0\alpha_0}{2}e^{-i\omega t}\int_0^{\infty}\sum_{n=0}^\infty \frac{({i k\sqrt{y^2+l^2}}) ^n}{n!\sqrt{y^2+l^2}}\frac{1}{2} dl^2\hat{z}
将l2l^2看作积分变量，分别对各项积分得：
A⃗=μ0α02e−iωt∑n=0∞(ik)n(n+1)![(y2+L2)n+12−yn+1]z^ \vec{A}=\frac{\mu_0\alpha_0}{2}e^{-i\omega t}\sum_{n=0}^\infty \frac{({i k}) ^n}{(n+1)!}[(y^2+L^2)^{\frac{n+1}{2}}-y^{n+1}]\hat{z}
这里的LL应该趋于无穷大，这样会导致∣A∣|A|发散，这是因为磁矢势的大小与零势点的位置有关系，对磁矢势求旋度可以消除零势点位置的影响，我们可以先求出磁场，然后再对L取无穷：
B⃗=∇×A⃗=∂A∂yx^=μ0α02e−iωt∑n=0∞(ik)n(n+1)![(n+1)2(y2+L2)n−122y−(n+1)yn]x^ \vec{B}=\nabla\times \vec{A}=\frac{\partial A}{\partial y} \hat{x}=\frac{\mu_0\alpha_0}{2}e^{-i\omega t}\sum_{n=0}^\infty \frac{({i k}) ^n}{(n+1)!}[\frac{(n+1)}{2}(y^2+L^2)^ {\frac{n-1}{2}}2y-(n+1)y^{n}]\hat{x}
将求和级数凑成指数，容易得到：
B⃗=μ0α02e−iωt∑n=0∞[(iky2+L2)nn!yy2+L2−(iky)nn!]x^=−μ0α02ei(ky−ωt)x^+μ0α02e−iωteiky2+L2yy2+L2x^ \vec{B}=\frac{\mu_0\alpha_0}{2}e^{-i\omega t}\sum_{n=0}^\infty[\frac{(ik\sqrt{y^2+L^2})^n}{n!}\frac{y}{\sqrt{y^2+L^2}}-\frac{(iky)^n}{n!}]\hat{x}=-\frac{\mu_0\alpha_0}{2}e^{i(ky-\omega t)}\hat{x}+\frac{\mu_0\alpha_0}{2}e^{-i\omega t}\frac{e^{ik\sqrt{y^2+L^2}}y}{\sqrt{y^2+L^2}} \hat{x}
容易知道，第二项在L→∞L\rightarrow \infty时为0，因此最后结果为：
B⃗=−μ0α02ei(ky−ωt)x^ \vec{B}=-\frac{\mu_0\alpha_0}{2}e^{i(ky-\omega t)}\hat{x}
又可求出电场：
E⃗=−wk^k×B⃗=−μ0α0c2ei(ky−ωt)z^ \vec{E}=-\frac{w\hat{k}}{k}\times\vec{B}=-\frac{\mu_0\alpha_0c}{2}e^{i(ky-\omega t)}\hat{z}
这和平面波的表达式一致。

2、无限长导线电流源

电磁波在柱坐标和球坐标下的本征解
假设在z轴上有变化的电流，I=I0e−iωtI=I_0e^{-i\omega t}。接下来求解，oxy平面内，离原点r处的磁矢势A⃗\vec{A}。同样要考虑不同位置处小段电流对r处的相位不同：
A⃗=μ0I04πe−iωt∫−∞∞eikr2+l2dlr2+l2z^ \vec{A}=\frac{\mu_0I_0}{4\pi}e^{-i\omega t}\int_{-\infty}^{\infty}\frac{e^{ik\sqrt{r^2+l^2}}dl}{\sqrt{r^2+l^2}}\hat{z}
这个积分很难处理，先求电场表达式：
E⃗=−∂A⃗∂t=−μ0I0ω4πe−iωt(−i)∫−∞∞eikr2+l2dlr2+l2z^ \vec{E}=-\frac{\partial\vec{A}}{\partial t}=-\frac{\mu_0I_0\omega}{4\pi}e^{-i\omega t}(-i) \int_{-\infty}^{\infty}\frac{e^{ik\sqrt{r^2+l^2}}dl}{\sqrt{r^2+l^2}}\hat{z}
对比上一节柱坐标系下的本征解形式，只能先猜想积分项就是m=0时的汉克函数，实际上，我在一本参考书的课后习题上发现确实有这种关系(奚定平.贝塞尔函数.北京：高等教育出版社，1999)
电磁波在柱坐标和球坐标下的本征解
我们需要用欧拉公式展开指数项，得到：
{∫−∞∞cos⁡(kr2+l2)r2+l2dl=−πN0(kr)i∫−∞∞sin⁡(kr2+l2)r2+l2dl=iπJ0(kr) \begin{cases} \int_{-\infty}^\infty \frac{ \cos(k \sqrt{r^2+l^2})}{\sqrt{r^2+l^2}}dl= -\pi N_0(kr)\\ i\int_{-\infty}^\infty \frac{ \sin(k \sqrt{r^2+l^2})}{\sqrt{r^2+l^2}}dl=i\pi J_0(kr) \end{cases}
因此电场可化为：
E⃗=−∂A⃗∂t=−μ0I0ω4πe−iωt(−i)[−πN0(kr)+iπJ0(kr)]z^=−μ0I0ω4e−iωt[J0(kr)+iN0(kr)]z^ \vec{E}=-\frac{\partial\vec{A}}{\partial t}=-\frac{\mu_0I_0\omega}{4\pi}e^{-i\omega t}(-i) [-\pi N_0(kr)+i \pi J_0(kr)]\hat{z}=-\frac{\mu_0I_0\omega}{4}e^{-i\omega t} [J_0(kr)+i N_0(kr)]\hat{z}
为了保险起见，我们可以用数值的方式验证上述积分是否正确：(matlab里面有现成的贝塞尔和诺伊曼函数)
电磁波在柱坐标和球坐标下的本征解

六、软件模拟

最后，为了验证这种源确实可以产生对应的电磁波，我打算用电磁波软件comsol3.5模拟。

1、平面波

初始条件及环境：

下图是无限长空腔的俯视图：
电磁波在柱坐标和球坐标下的本征解

尺寸及模块：长X宽=8mX4m，RF模式下的TE模式；
region 1 为真空；region 2为PML系数物质；
border 1：垂直于纸面的面电流源α=1⋅e−iωtA\alpha=1·e^{-i\omega t} A，频率f=2π/ω=0.5GHzf=2\pi /\omega=0.5G Hz，及波长λ=0.6m\lambda=0.6m
border 2: 连续边界条件，即磁场平行分量相等：n⃗×(H1⃗−H2⃗)=0\vec{n}\times(\vec{H_1}-\vec{H_2})=0

仿真结果：

电磁波在柱坐标和球坐标下的本征解
颜色表示电磁场的值，红色为正，蓝色为负。容易看出，无限大平面电流产生的确实是平面波，波长为0.6m。

选取y=0处的数据点磁场强度H−xH-x进一步验证：
电磁波在柱坐标和球坐标下的本征解
重新推导磁场强度，带入数值可以得到：
H⃗=(α02)ei(kx−ωt)y^=12cos⁡(kx)y^ \vec{H}=(\frac {\alpha_0}{2})e^{i(kx-\omega t)}\hat{y}= \frac {1}{2}\cos(kx)\hat{y}
这里令t=0t=0，取实部，可以发现理论计算比实际小一半，什么原因呢？

因为理论计算时考虑了波向y的正方向和负方向两边扩散，而模拟时相当于负方向的波和正方向叠加，因此会有二倍的效果。

二、柱面波

下图是无限长圆柱的俯视图：
电磁波在柱坐标和球坐标下的本征解

尺寸及模块：三个园半径分别为0.2、5、6m，RF模式下的TE模式；
region 1 为真空；region 2为PML系数物质；
border 1：垂直于纸面的面电流源I=2πrα⋅e−iωtA=1e−iωtI=2\pi r \alpha · e^{-i\omega t} A=1e^{-i\omega t}，频率f=2π/ω=0.5GHzf=2\pi /\omega=0.5G Hz，或波长λ=0.6m\lambda=0.6m
border 2: 连续边界条件。
其他边界条件：PEC完美电导体，即电场无垂直分量Ez=0E_z=0;

值得注意的是，这里的电流源是有一定粗细的导线，这和理论推导时的假设不一样，仿真结果会在数值上于理论不一致，但仍然是贝塞尔函数。

模拟结果：

电磁波在柱坐标和球坐标下的本征解
初步可以看出这是柱面波，波长接近0.6m。

选取y=0的半径上获取数据，做出电场随半径变化图E−rE-r：
电磁波在柱坐标和球坐标下的本征解
将数据带入理论表达式有：(t=0，取实部)
E⃗=−μ0I0ω4e−iωt[J0(kr)+iN0(kr)]z^=−987J0(kr)z^ \vec{E}=-\frac{\mu_0I_0\omega}{4}e^{-i\omega t} [J_0(kr)+i N_0(kr)]\hat{z}=-987J_0(kr)\hat{z}
两者为什么会不一样呢？这就是源的问题：如果把comsol模拟的源"0.2m的圈"改成点源，可以得到与理论一致的结果。

总结

  • 拉普拉斯算符对矢量的作用不同于标量，前者还要考虑不同坐标系下单位矢量随位置的变化。
  • 平面波的波源是无限大平面的震荡电流（二维电流），柱面波的波源是无线长导线的震荡电流（一维波源），球面波的波源是偶极震荡（零维电流）。
  • 虽然不会证明，还是列出零阶贝塞尔函数和诺伊曼函数的积分表达式：
    {N0(r)=−1π∫−∞∞cos⁡(r2+l2)r2+l2dlJ0(r)=1π∫−∞∞sin⁡(r2+l2)r2+l2dl \begin{cases} N_0(r)=-\frac{1}{\pi}\int_{-\infty}^\infty \frac{\cos(\sqrt{r^2+l^2})}{\sqrt{r^2+l^2}}dl\\ J_0(r)=\frac{1}{\pi}\int_{-\infty}^\infty \frac{ \sin( \sqrt{r^2+l^2})}{\sqrt{r^2+l^2}}dl \end{cases}
  • 仿真模拟需要在数值上考虑更多细节，对理论的理解会更深一层。

附录：

电磁波在柱坐标和球坐标下的本征解

相关推荐