#MATLAB编写雷电波激励下空间任一点场的时域频域波形，其中雷电波基电流采用Heilder模型

#MATLAB编写雷电波激励下空间任一点场的时域频域波形其中雷电波基电流采用Heilder模型供研究分析使用 #程序包含详细注释本人在2020a版本均可运行最近在搞雷电电磁脉冲的仿真发现很多论文里的Heidler模型代码写得云里雾里。今天咱们直接上干货手把手教你用MATLAB撸出雷电波激励下的空间场分布。别被公式吓到其实就三层窗户纸——电流模型、时域计算和频域转换。先看Heidler模型的实现这个电流波形是出了名的难调参数。核心公式长这样function I heidler_current(t, I0, eta, tau1, tau2, n) % 雷电基电流Heidler模型 % t - 时间序列秒 % I0 - 峰值电流A % eta - 修正系数 % tau1 - 波头时间常数秒 % tau2 - 波尾时间常数秒 % n - 陡度系数 term1 (t/tau1).^n ./ (1 (t/tau1).^n); term2 exp(-t/tau2); I (I0/eta) * term1 .* term2; end这里有个坑tau1和tau2的单位容易搞错记得换成微秒量级。比如典型参数I030kAtau10.5e-6tau25e-6时生成的波形在0.1微秒处就会冲到峰值。#MATLAB编写雷电波激励下空间任一点场的时域频域波形其中雷电波基电流采用Heilder模型供研究分析使用 #程序包含详细注释本人在2020a版本均可运行接下来要算空间场这里用简化版偶极子模型别急后面会解释为啥能这么用function [E,H] calc_fields(r, theta, t, current, dt) % 计算空间任意点电磁场 % r - 观测点距离米 % theta - 仰角弧度 % t - 时间序列 % current - 电流时域数据 % dt - 时间步长 c 3e8; % 光速 mu0 4*pi*1e-7; eps0 8.85e-12; % 时延处理 delay floor(r/c / dt); if delay length(t) error(时延超过时间序列长度); end current_delayed [zeros(1,delay), current(1:end-delay)]; % 电场计算辐射场项 E_rad (mu0/(4*pi*r)) * sin(theta) * gradient(current_delayed, dt); % 磁场计算 H_phi (1/(4*pi*r)) * sin(theta) * (current_delayed/(c*r) gradient(current_delayed, dt)/c^2); % 这里省略了感应场和静电场的计算 E E_rad; H H_phi; end注意第14行的梯度计算gradient函数这是MATLAB自带的差分工具。但有个玄机——当dt太小时数值微分会飘建议时间步长控制在0.1微秒以上。频域分析这块很多教程都让直接fft完事其实要做窗函数处理function [f, spectrum] analyze_spectrum(signal, dt) % 频域分析带汉宁窗 N length(signal); window hann(N); % 汉宁窗抑制频谱泄露 Y fft(signal .* window); f (0:N-1)/(N*dt); spectrum 20*log10(abs(Y)/max(abs(Y))); % 归一化dB值 end这里有个工程经验雷电波的频带一般在10kHz-10MHz所以绘图时建议xlim([1e4 1e7])用log坐标更直观。最后来个主程序串联流程%% 主程序 dt 1e-8; % 10ns步长 t 0:dt:1e-5; % 10微秒时长 % 生成Heidler电流典型首次回击参数 I heidler_current(t, 3e4, 0.7, 0.5e-6, 5e-6, 5); % 计算100米外地面场thetapi/2 r 100; theta pi/2; [E, H] calc_fields(r, theta, t, I, dt); % 频域转换 [f, E_spectrum] analyze_spectrum(E, dt); %% 画图部分 figure; subplot(2,1,1) plot(t*1e6, I/1e3) xlabel(时间μs), ylabel(电流kA) title(雷电基电流时域波形) subplot(2,1,2) semilogx(f, E_spectrum) xlim([1e4 1e7]), grid on xlabel(频率Hz), ylabel(归一化场强dB) title(电场频域特性)跑这个程序时可能会遇到两个问题一是内存不足把时间序列调短点二是高频部分出现毛刺在calc_fields里加个低通滤波就行。最后出来的图应该能看到典型的双指数波衰减频域在3MHz左右有个明显拐点——这和tau1参数直接相关调参时盯着这个点就能快速验证模型准确性。