【计算电磁学】简单的C语言一维FDTD仿真程序
通过时域有限差分法(FDTD)可以利用上一时刻的磁场和同一时刻相邻空间的电场递推下一时刻的磁场以及利用上一时刻的电场和同一时刻相邻空间的磁场递推下一时刻的电场,由此可以计算出整个空间的场分布。
以下为仿真电磁波在自由空间中传播的程序,计算区域中含200个电场点和磁场点。
高斯激励源的简单一维FDTD仿真程序:
/* 硬源激励的基本FDTD仿真 */
#include <stdio.h>
#include <math.h>
#define SIZE 200
int main ()
{
double ez[SIZE] = {0.}, hy[SIZE] = {0.}, imp0 = 377.0; //{0.}保证数组初始化为0
int qTime, maxTime = 250, mm; //qTime为时间步,总步数为maxTime
FILE *fp;
fp = fopen("1.txt", "w+");
/* 开始时间循环 */
for (qTime = 0; qTime < maxTime; qTime++)
{
/* 递推磁场 */
for (mm = 0; mm < SIZE - 1; mm++)
hy[mm] = hy[mm] + (ez[mm + 1] - ez [mm]) / imp0;
/* 递推电场 */
for (mm = 1; mm < SIZE; mm++)
ez[mm] = ez[mm] + (hy[mm] - hy[mm - 1]) * imp0;
fprintf(fp, "%.6g\n", ez[50]); //记录50时刻的电场值并打印到txt文件
ez[0] = exp(-(qTime - 30.) * (qTime - 30.) / 100.); //高斯脉冲源
printf ("%g\n", ez[50]);
}
fclose(fp);
return 0;
}
作图的MATLAB程序如下:
fid3=textread('1.txt', '%n', 'whitespace', '');
[m,n]=size(fid3);
Ox=1:m;
figure;
plot(Ox,fid3,'-');
grid on;
横坐标为时间步,纵坐标为Ez[50]。由于我们设定稳定常数为1,所以场每个时间步移动一个空间步长,记录到的信号从源点延时50个时间步。
参考
[1] John B. Schneider, Understanding the Finite-Difference Time-Domain Method [M]