简介:本资源是面向本科及硕士阶段电磁场与天线课程教学、科研入门与仿真实践的Matlab专项工具包,聚焦无穷小偶极子天线近场与远场电磁特性的数值建模与可视化分析。通过严格依据经典电磁理论推导的电场、磁场及功率密度表达式,在Matlab中实现空间三维场分布计算与动态绘图,帮助学习者直观理解近场感应区、辐射近场与远场区的物理特征差异及过渡规律。压缩包共6个文件(1.08MB),含核心仿真脚本main.m、3幅关键结果图(png)、说明文档txt及推广引导图,结构简洁、即开即用;所有代码兼容Matlab 2014a/2019a,附带运行效果图,降低初学者环境配置与调试门槛。目前已有23人学习下载,适用于电磁场理论验证、天线基础仿真训练及课程设计参考,可直接用于教学演示或个人研究复现。 最近我把一个做了很久的Matlab工程重新整理了一遍,正好是围绕“无穷小偶极子天线近场和远场模拟”的项目。说起来这算是我读研时一个很基础、但后来反复拿出来用的练习:建模对象极其简单,简单到只有一个电流元,但真要在Matlab里把近场特性、远场方向图、场强衰减趋势全部可视化出来,并且保证每条曲线都和解析解对得上,遇到的问题比预想多得多。把整个项目打包成zip后,结构也很清晰,适合直接解压运行,所以干脆写篇文章把里面的关键细节都拆开讲一遍。
这里说的“无穷小偶极子”,也叫赫兹偶极子,是电磁场和天线理论里最经典的理想模型。网上很多讲近场远场的资料只给公式,没给出怎么在Matlab里把各个区域画出来;而很多Matlab例子又只顾着出图,没有把近场项和远场项的物理意义讲清楚。这次的项目把这两块接起来了:从球坐标系下电场磁场的完整解析式出发,在Matlab里逐项计算,最后得到近场区域的场分布图、远场的方向图以及随距离变化的衰减曲线。如果你正在学电磁场与天线,或者工作中需要用一个解析可算的标准天线去校准其他仿真软件,这个工程可以直接拿来当参考。
下面我会把项目从原理、代码实现、结果解读到常见问题完整过一遍。全文不会只贴代码然后让你自己看,而是把每个关键步骤背后为什么要这么写、算出来的东西怎么理解,都尽量讲透。
1. 项目概述:这个zip包里到底装了什么
1.1 无穷小偶极子在天线理论中的地位
无穷小偶极子可以理解成一段长度远小于工作波长、上面电流又处处等幅同相的电流元。实际天线,不管是一根偶极子天线还是复杂的阵列天线,理论上都能拆成无数个这样的小电流元来叠加分析。所以在天线教科书里,这个模型永远是最前面的一章,地位类似傅里叶分析里面的正弦基函数,是整个大厦的地基。
和贴片天线、螺旋天线这类结构不同,无穷小偶极子的场有严格闭合解,不需要建模网格,不需要数值求解器。这意味着做仿真的时候,你得到的是“标准答案”,完全可以拿来验证其他数值算法算出来的结果靠不靠谱。我在实际项目里经常遇到这种情况:用HFSS或CST仿真某根天线之前,先拿无穷小偶极子跑一遍流程,确认软件设置、边界条件、网格划分有没有人为错误。它的作用就是这一类“参考基准”。
1.2 近场与远场到底怎么划分
近场和远场的划分,物理根源在于场量随距离的衰减规律不同。对偶极子而言,电场里有随1/r³衰减的分量,也有随1/r²、1/r衰减的分量;磁场里有随1/r²、1/r衰减的分量。当距离很近时,高次项远大于远场项,场的行为接近静态场,能量在源和周围空间之间来回交换,这部分叫无功近场;到了距离足够远的地方,1/r项占绝对主导,能量真正向外辐射出去,就是远场。
工程上最常用的判据是r远大于2D²/λ,其中D是天线最大尺寸。但无穷小偶极子的D趋近于零,这个判据直接失效。更自然的划分标准是看kr这个无量纲量:当kr=1时,正好对应r=λ/2π≈0.159λ。这个位置常被当作无功近场和辐射近场的粗略分界点。实际工程里,通常取几个波长以上作为远场区,此时1/r项已经占绝对主导。
1.3 包内文件结构与运行环境
整个zip包解压后是这样一个结构:
- main_dipole_nearfar.m:主脚本,一键运行,负责计算和出图
- hertzian_fields.m:核心场量计算函数,输入半径、角度、电流参数,输出Er、Etheta、Hphi
- plot_radial_decay.m:绘制径向衰减曲线的辅助函数
- plot_direction_pattern.m:绘制E面和H面方向图的辅助函数
代码用的都是Matlab基础功能,不需要额外工具箱。唯一要求是版本别太老,建议R2016a以上,因为用了polarplot这类较新函数。如果你要二次开发,只要理清hertzian_fields这个函数,其他部分都只是数据后处理和绘图。
2. 仿真原理:写代码前必须吃透的公式
2.1 球坐标系下无穷小偶极子的完整解析式
设偶极子沿z轴放置,位于坐标原点,电流为I₀,长度为dl,自由空间波数k=2π/λ,波阻抗η=120π≈377Ω。在球坐标系下,任意观察点(r, θ, φ)处的场可以写成:
E_r = (η I₀ dl cosθ) / (2π r²) · (1 + 1/(jkr)) · e^(-jkr)
E_θ = (j η k I₀ dl sinθ) / (4π r) · [1 + 1/(jkr) − 1/(kr)²] · e^(-jkr)
H_φ = (j k I₀ dl sinθ) / (4π r) · (1 + 1/(jkr)) · e^(-jkr)
因为是理想点源模型,磁场只有φ方向一个分量,电场有r和θ两个分量,三者互相垂直。注意这里用的都是复振幅形式,e^(-jkr)表示球面波向外传播的相位滞后,j是虚数单位。实际仿真里所有量都是复数,出图时要取模或取实部,这一点新手很容易忽略。
2.2 近场区哪些项占主导,物理上意味着什么
近场区对应kr远小于1,也就是r远小于λ/2π。在这个范围内,E_r括号里的1/(jkr)项非常大,E_θ括号里的1/(kr)²项也非常大,H_φ里的1/(jkr)项同样可观。于是场量中占主导的是:
E_r ≈ (η I₀ dl cosθ) / (2π r²) · (1/(jkr)) · e^(-jkr)
E_θ ≈ −(j η k I₀ dl sinθ) / (4π r) · (1/(kr)²) · e^(-jkr)
H_φ ≈ (j k I₀ dl sinθ) / (4π r) · (1/(jkr)) · e^(-jkr)
其中E_r和E_θ都与1/r³成正比,H_φ与1/r²成正比。这个特性非常接近静态偶极子的场。正因为这些高次项的存在,近场区存在大量能量存储,并没有真正跑出去,所以近场区也叫无功近场区。近场区还有一个显著特点:E_r纵向分量和E_θ横向分量可能量级相当,导致总电场并不像远场那样垂直于传播方向,而是有明显纵向成分。
仿真中如果只关心辐射特性,近场区这些能量完全可以忽略;但如果研究天线附近的耦合、感应场,或者做近场测量,就必须把这些项都保留。这也是为什么我坚持在代码里计算完整表达式,而不是直接跳到远场近似。
2.3 远场区如何简化,方向图为什么是“8字形”
远场区对应kr远大于1,此时1/(jkr)和1/(kr)²这些项相对1可以忽略。电场和磁场简化为:
E_θ ≈ (j η k I₀ dl sinθ) / (4π r) · e^(-jkr)
H_φ ≈ (j k I₀ dl sinθ) / (4π r) · e^(-jkr)
E_r ≈ 0
这组远场近似公式信息量很大。第一,E_θ和H_φ都是按1/r衰减,是真正的辐射场;第二,两者比值E_θ/H_φ = η,等于自由空间波阻抗,说明远场中电场和磁场同相位、互相垂直,能量沿径向向外传播;第三,方向性完全包含在sinθ因子里。θ=90°时场最强,θ=0°和180°也就是偶极子轴线方向上,场严格为零。如果绘制E面方向图,得到的正是教科书里经典的那个“8字形”或说“甜甜圈”的剖面图。H面方向图则是在垂直于偶极子轴的平面内,任意方位角上场强都相等,是一个圆。
搞清楚这些之后,Matlab写起来就不容易跑偏了。因为你知道近场要重点观察哪些项,远场又要观察什么,不会拿着结果一阵乱画。
3. Matlab仿真实现全流程
3.1 工程顶层的组织方式:先搞清单文件职责
我写这个程序的时候,刻意把“计算”和“绘图”分开。很多人爱把所有代码堆在一个脚本里,前面算完后面画,最后代码一长串,想改参数就得很小心地找。我这次拆成了4个文件,主脚本只负责定义参数和调用,核心计算交给hertzian_fields函数,绘图逻辑再单独放两个辅助函数里。
这样拆分的好处是:当你换了频率,或者想改偶极子长度,只需要在主脚本里改几个变量;当你想从直角坐标观察近场分布时,只需要改网格生成方式,不需要碰场量计算公式。虽然这个项目本身不大,但这种组织习惯放到大型工程里能节约大量调试时间。
3.2 核心函数hertzian_fields的实现细节
打开hertzian_fields.m,核心就几行。我尽量用数组运算代替for循环,因为Matlab对矩阵操作做了深度优化,直接对整个矩阵求场可以省去一层甚至两层循环,计算效率高不少。以下是完整实现:
function [Er, Etheta, Hphi] = hertzian_fields(r, theta, I0, dl, k, eta) % hertzian_fields 计算无穷小偶极子在球坐标系下的场量 % 输入: % r - 观察点到偶极子的距离,可以是标量或矩阵 % theta - 极角,单位弧度,与r同尺寸 % I0 - 电流幅度,单位A % dl - 偶极子长度,单位m % k - 波数,单位rad/m % eta - 介质波阻抗,真空约120*pi % 输出: % Er, Etheta, Hphi - 复振幅场强,单位V/m或A/m kr = k * r; ejkr = exp(-1j * kr); sinT = sin(theta); cosT = cos(theta); % 电场径向分量 Er = (eta * I0 * dl ./ (2 * pi * r.^2)) .* cosT .* ... (1 + 1 ./ (1j * kr)) .* ejkr; % 电场theta分量 Etheta = (1j * eta * k * I0 * dl ./ (4 * pi * r)) .* sinT .* ... (1 + 1 ./ (1j * kr) - 1 ./ (kr).^2) .* ejkr; % 磁场phi分量 Hphi = (1j * k * I0 * dl ./ (4 * pi * r)) .* sinT .* ... (1 + 1 ./ (1j * kr)) .* ejkr; end这里有个容易踩的坑:r不能取0,否则1./r和1./kr会直接变成Inf或NaN。我在主脚本里用logspace生成从0.001λ到10λ的范围,就是为了避开原点奇点。实际计算近场时,也不用非得取到那么小,一般从0.001λ开始就足够了,再小数值意义不大,毕竟理想点源模型在这个尺度已经不“物理”了。
另一个容易出错的地方是点乘和乘法的区别。r、theta如果是矩阵,就必须用.*和./,不用的话Matlab会尝试做矩阵乘法,报错或者得到完全错误的结果。初学者在这里很容易卡一下,我刚开始也吃过亏。
3.3 主脚本:径向衰减曲线与方向图绘制
主脚本里面我做了三块内容:第一,固定theta=π/2,让r从0.001λ扫到10λ,画近场到远场过渡的场强衰减曲线;第二,固定r=10λ,让theta从0扫到2π,画E面方向图;第三,固定theta=π/2,让phi从0扫到2π,画H面方向图。
clear; clc; close all; % 基础参数 c0 = 3e8; f = 1e9; % 1 GHz lambda = c0 / f; % 波长 0.3 m k = 2 * pi / lambda; eta = 120 * pi; I0 = 1; % 电流幅度 dl = lambda / 100; % 偶极子长度,远小于波长 % 第一幅图:径向衰减曲线,观察赤道面电场 r_log = logspace(-3, 1, 500) * lambda; theta_fix = pi / 2; [Er, Etheta, Hphi] = hertzian_fields(r_log, theta_fix, I0, dl, k, eta); figure; loglog(r_log / lambda, abs(Er), '-r', ... r_log / lambda, abs(Etheta), '-b', ... r_log / lambda, abs(Hphi) * eta, '--k', 'LineWidth', 1.5); xlabel('r / \lambda'); ylabel('场强幅度 (V/m)'); legend('|E_r|', '|E_\theta|', '|H_\phi| \times \eta', 'Location', 'best'); grid on; title('赤道面上场量随距离的衰减'); % 第二幅图:E面远场方向图 r0 = 10 * lambda; theta_scan = linspace(0, 2 * pi, 721); [~, Etheta_far] = hertzian_fields(r0 * ones(size(theta_scan)), ... theta_scan, I0, dl, k, eta); norm_Etheta = abs(Etheta_far) / max(abs(Etheta_far)); figure; polarplot(theta_scan, norm_Etheta, 'b'); title('E面远场方向图 (r = 10\lambda)'); % 第三幅图:H面远场方向图 phi_scan = linspace(0, 2 * pi, 721); theta_h = pi / 2 * ones(size(phi_scan)); [~, Etheta_h] = hertzian_fields(r0 * ones(size(phi_scan)), theta_h, I0, dl, k, eta); norm_Etheta_h = abs(Etheta_h) / max(abs(Etheta_h)); figure; polarplot(phi_scan, norm_Etheta_h, 'r'); title('H面远场方向图 (r = 10\lambda)');这里有几个操作上的考虑。第一,径向衰减曲线用loglog而不是plot,是因为场量动态范围跨了好几个数量级,线性坐标下小距离处的曲线会被压缩到几乎看不见。第二,Etheta_far虽然只取了E_theta,但调用函数时还是需要用~占位符接收Er,否则语法会报错。第三,方向图做归一化是必要的,因为方向图关心的是相对大小,不关心绝对场强,归一化后不同频率、不同电流幅度下画出来的图都能直接对比。
如果你希望观察近场在某个平面上的空间分布,还可以在主脚本后面追加生成二维伪彩图。方法是建立一个直角坐标网格,比如xz平面,将网格点的坐标转换成r和θ,再调用场函数。这里贴一段我当时用的扩展代码:
% 近场xz平面二维分布图 xx = linspace(-lambda, lambda, 401); zz = linspace(-lambda, lambda, 401); [X, Z] = meshgrid(xx, zz); R = sqrt(X.^2 + Z.^2); R(R < lambda * 1e-3) = lambda * 1e-3; % 屏蔽原点奇点 THETA = acos(Z ./ R); [Er2, Et2, ~] = hertzian_fields(R, THETA, I0, dl, k, eta); Etotal = sqrt(abs(Er2).^2 + abs(Et2).^2); figure; imagesc(xx/lambda, zz/lambda, 20*log10(Etotal/max(Etotal(:)))); axis xy; axis equal; colorbar; xlabel('x / \lambda'); ylabel('z / \lambda'); title('近场电场总幅度分布 (dB)');这段代码里最关键的是R(R<阈值)=阈值。由于网格原点处r=0,不处理就会得到Inf,整个色图都会被破坏,这个处理属于参数层面的小技巧,但非常实用。用dB画色图也很有讲究,近场动态范围太大,线性色标根本看不出层次。
4. 仿真结果怎么读、怎么用
4.1 远场E面和H面方向图的实际图像
运行主脚本后,第二幅图会显示一个标准“8字形”极坐标曲线,在θ=90°方向幅度为1,在θ=0°和180°左右幅度为0。这个曲线对应偶极子天线的E面方向图。如果你看三维波瓣,它像个甜甜圈,围绕z轴旋转对称。H面方向图则是一个圆,因为赤道面上各个方向场强一样大。
很多人第一次画出来,可能会觉得“8字形”太平滑,怀疑是不是代码写错了。其实这正是解析解应该有的样子。如果以后用HFSS仿真一根有限尺寸的偶极子,方向图会和这个理想“8字形”有细微差别,主要因为振子长度不是无穷小,会带来方向图主瓣稍微变窄、轴向不那么严格为零等偏差,但大趋势仍然一致。
4.2 径向衰减曲线揭示了近场到远场的过渡
第一幅图是我认为信息量最大的。在loglog坐标系下,可以看到三个区域:r很小的时候,|Er|和|Eθ|的斜率很陡,因为场按1/r³衰减,对应loglog图的斜率约-3;再远一点,斜率逐渐趋缓,变成约-2,这是1/r²项在起作用;到了r大约λ/2π之后,斜率接近-1,就是标准辐射场的1/r衰减。
|Hφ|的情况类似,从近场的1/r²衰减逐渐过渡到远场的1/r衰减。值得注意的是,|Hφ|×η这条曲线和|Eθ|在远场区会完全重合,这就是前面说的E_θ/H_φ=η这个波阻抗关系在图形上的直接体现。我在实际看仿真曲线时,经常用这个特征判断自己设置的远场半径是不是足够大。如果两条曲线到r=10λ还没重合,说明r取得不够,要继续加大观察距离。
4.3 近场与远场特性的表格化对比
为了方便项目报告或演示,我整理了一个近场远场特性对比表,可以放到文档里:
| 特性 | 近场区 (kr << 1) | 远场区 (kr >> 1) |
|---|---|---|
| 占主导的分量 | E_r、E_θ按1/r³衰减,H_φ按1/r²衰减 | E_θ、H_φ按1/r衰减 |
| 电场纵向分量 | 明显存在,不可忽略 | 可忽略,几乎为零 |
| 电场与磁场相位 | 电场和磁场有约90°相位差 | 电场和磁场同相位 |
| 波阻抗 | 与距离和角度有关 | 恒等于η |
| 能量特性 | 无功功率为主,能量往返存储 | 有功功率为主,能量向外辐射 |
| 工程应用 | 近场探测、耦合分析、近场测量 | 天线远场方向图、增益测量 |
这张表写报告时可以直接用。近场测量这类工程场景,不能拿远场假设去近似;远场测量转台的设计、暗室尺寸计算,则要按远场条件来。模拟背后这些物理区别直接影响工程决策。
5. 常见问题与排查技巧实录
5.1 对数坐标下曲线出现断裂或NaN
很多人在loglog画图时,发现某些点没显示出来,或者曲线中间断开一段。最常见的原因是场量在某些角度上严格等于0,比如θ=0时sinθ=0,导致Hφ和Eθ等于0。loglog取对数时遇到0会变成-Inf,曲线自然断裂。E面方向图在θ=0处为0是正常现象,但画径向衰减曲线时如果选择的θ恰好是0,那整条曲线都会是0。
解决办法是:画径向衰减曲线时固定θ=π/2,不要取0;画二维分布图时,用20*log10(数据)之前先对数据做一个很小的下限裁剪,比如Etotal(Etotal<1e-10)=1e-10,这样可以避免-Inf破坏色标范围。
5.2 复数场量的模值和相位问题
Matlab里算出Er、Etheta、Hphi都是复数,直接plot会画出实部虚部两条线,很多人一时没反应过来,以为程序错了。我自己的习惯是,对幅度图一律用abs()取模,需要观察相位分布时用angle()取相位,再用unwrap使相位曲线连续。
近场区相位空间变化非常剧烈,不同位置间的相位差可能快速跳动,这是高次项占主导导致的正常结果。如果某天你在相位图上看到明显的“跳变”条纹,先别着急改代码,想想是不是这个原因。真要分析相位的话,建议对角度和距离的分辨率都提高一些,让相位变化有足够采样点。
5.3 极坐标方向图出现负值或非对称
polarplot函数要求输入的是非负数,因为极径本身不能为负。如果你直接拿abs(Etheta)去画,它天然非负,没问题。但如果你对方向图做了某种运算,比如把远场Eθ减去近场Eθ,然后再取绝对值画,图形可能出现“花瓣向内凹”的奇怪形状,很多人就蒙了。方向图应该是对θ=90°对称的,如果不对称,先检查生成theta时是否从0到2π等间隔采样,有没有漏掉某个区间。
还有一个容易犯的错:E面方向图应该用theta从0到2π扫描,但有些人在画时只用了0到π,这样只画出上半部分,看起来像一个半圆,于是误以为程序错了。实际上偶极子在下方也有辐射,必须扫满2π才能看到完整的“8字形”。
5.4 计算量过大导致脚本卡死
如果把网格设得特别密,比如xy平面网格用1000×1000,还要计算复数场再画图,Matlab可能会忽然卡住。虽然单次场量计算并不复杂,但100万个网格点上的复数运算也会让普通电脑明显变慢。建议先用401×401这类中等密度,确认结果合理后再加密。另外,hertzian_fields函数里没有用循环,所以网格再大也主要是内存和计算速度的问题,不至于指数级恶化。
注意:近场二维分布图在origin处必须屏蔽,r<λ/1000的点全部替换成阈值。否则Inf和NaN会污染整个colorbar范围,导致图像一片深蓝或一片红色,看起来完全失去层次。
6. 从仿真到工程应用的个人体会
6.1 妙用无穷小偶极子做解析基准
工程上,无穷小偶极子最常见的用途是作为解析基准天线。比如我拿到一个新的电磁仿真软件,不确定它的辐射边界条件设置有没有问题,就会先建一个偶极子模型,把仿真结果和Matlab解析解对比。如果辐射方向图、远场增益、输入阻抗都对得上,说明设置基本可靠。这种基准测试的思路,远不止这一个小项目里能用,几乎可以平移到所有场仿真软件的学习和调试中。
做近场测量或天线暗室校准的时候,某些频段甚至可以直接用一个小型探针当参考源,它产生的场结构非常接近无穷小偶极子,方便用这个理论模型反推探头位置、极化性能等参数。所以别小看这个基础项目,它背后是大量工程验证流程的共同起点。
6.2 后续可以怎么扩展
这个项目做完后,按同样套路可以延伸出不少方向。第一,把无穷小偶极子换成有限长度的对称振子,只要叠加多个电流元或者用更严格的电流分布模型,场量公式会复杂一些,但Matlab架构不用大改。第二,把单个偶极子扩展成二元阵或直线阵,通过阵因子叠加两个源的场,就可以直观看到方向图变化,这也是阵列天线仿真的入门门槛。第三,把理想自由空间介质参数换成有损耗介质,观察近场场的分布变化,这和探地雷达、生物医学电磁感应等应用都有关系。
实现这类扩展,核心都是复用hertzian_fields这个函数,然后用叠加原理把多个源的场加起来。叠加时注意场是矢量,要按球坐标分量分别叠加,而不是简单把幅度加起来,这是不少人容易忽略的点。
6.3 最后分享一个小技巧
我在实际跑这个项目时,最喜欢的一个验证方法,是把近场E_θ模值曲线在相邻两个距离点的比值取出来看。理论上如果两点都处在1/r³占主导的近场区,比值应该接近(r2/r1)³;如果两点都处在1/r主导的远场区,比值应该接近(r2/r1)。用这个简单的数值检验,就能快速判断自己选的观察点到底属于哪个区域,比单看图直观得多。这个技巧也适合以后处理其他天线仿真结果,一旦学会用“距离变化率”去判断场区,很多近场远场问题都能一眼看出门道。
本文还有配套的精品资源,点击获取