电磁波时域有限差分方法（FDTD）-C语言-FDTD和激励源和MUR边界

参考：

葛德彪, 闫玉波. 电磁波时域有限差分方法[M]. 西安电子科技大学出版社, 2011.

Schneider J B. Understanding the Finite-Difference Time-Domain Method[J]. 2013.

最近要用FDTD做仿真，所以找了以上两本书，FDTD确实是有点难的。

想了想流程：

1：FDTD

2：激励源

3：吸收边界

4：散射

现在完成了FDTD和激励源，的TM部分，虽然都是抄书的，但对于一个不懂电磁波电磁场，数学极差的我来说，光是看懂代码就已经竭尽全力了。

何况还遇到了C语言宏定义的坑。

然后做了一部分的MUR吸收，效果不咋的，不知道就是这样的还是我代码有问题。

上代码，基本就是抄书的，但是简洁了一点。。

就不先做注释了，等到最后全部写完再做注释。

// Boundary.cpp: 定义控制台应用程序的入口点。
//

#include "stdafx.h"
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define _CRT_SECURE_NO_WARNINGS

#define ALLOC_2D(PNTR, NUMX, NUMY, TYPE)					\
	PNTR = (TYPE *)calloc((NUMX) * (NUMY), sizeof(TYPE));	\
	if (!PNTR) {												\
		perror("ALLOC_2D");									\
		fprintf(stderr, "Allocation failed for" #PNTR ". Terminating.\n");	\
		exit(-1);											\
	}
#define SizeX g->sizeX
#define SizeY g->sizeY
#define CDTDS g->cdtds
#define Time g->time
#define MaxTime g->maxTime

#define Ez(M, N) g->ez[M * SizeY + N]
#define CA(M, N) g->ceze[M * SizeY + N]
#define CB(M, N) g->cezh[M * SizeY + N]

#define Hx(M, N) g->hx[M * (SizeY) + N]
#define CP(M, N) g->chxh[M * (SizeY) + N]
#define CQ(M, N) g->chxe[M * (SizeY) + N]

#define Hy(M, N) g->hy[M * SizeY + N]
//#define Chyh(M, N) g->chyh[M * SizeY + N]
//#define Chye(M, N) g->chye[M * SizeY + N]

struct Grid {
	double *hx, *chxh, *chxe;
	double *hy, *chyh, *chye;
	double *ez, *ceze, *cezh;

	int sizeX, sizeY;
	int time, maxTime;
	int type;
	double cdtds;	//库朗数
};
typedef struct Grid Grid;

void gridInit(Grid *g);
double driveSource(Grid *g, double time, double locX, double locY);
void writeFile(Grid *g);
void updateE(Grid *g);
void updateH(Grid *g);
void updateBoundary(Grid *g);
void EHInit(Grid *g);

float c = 3e8;
float PI = 3.1415926535;
float MU0 = PI * 4e-7;
float EPS0 = (1 / (36 * PI)) * 1e-9;
float SIGMA0 = 0;
float SIGMAm0 = 0;
float DELTA = 5e-3;
float DELTAT = DELTA / (sqrt(2.0) * c);

int main()
{
	int m, n;
	double source;
	Grid *g;
	ALLOC_2D(g, 100, 100, Grid);

	EHInit(g);
	gridInit(g);

	for (Time = 0; Time < 150; Time++)
	{
		updateH(g);

		updateBoundary(g);

		updateE(g);
		source = driveSource(g, Time, 0, 0);

		Ez(SizeX / 2, SizeY / 2) = source;
	}

	writeFile(g);

	getchar();
	return 0;
}
void EHInit(Grid *g)
{
	int n, m;
	for (m = 0; m < SizeX; m++)
	{
		for (n = 0; n < SizeY; n++)
		{
			Hy(m, n) = 0;
			Hx(m, n) = 0;
			Ez(m, n) = 0;
		}
	}
}

void gridInit(Grid *g)
{
	int m, n;

	SizeX = 100;
	SizeY = 100;
	MaxTime = 300;
	CDTDS = 1.0 / sqrt(2.0);
	ALLOC_2D(g->hx, SizeX, SizeY, double);
	ALLOC_2D(g->chxh, SizeX, SizeY, double);
	ALLOC_2D(g->chxe, SizeX, SizeY, double);
	ALLOC_2D(g->hy, SizeX, SizeY, double);
	ALLOC_2D(g->chyh, SizeX, SizeY, double);
	ALLOC_2D(g->chye, SizeX, SizeY, double);
	ALLOC_2D(g->ez, SizeX, SizeY, double);
	ALLOC_2D(g->ceze, SizeX, SizeY, double);
	ALLOC_2D(g->cezh, SizeX, SizeY, double);

	for (m = 0; m < SizeX; m++)
	{
		for (n = 0; n < SizeY; n++)
		{
			//Ez(m, n) = 0;
			CA(m, n) = (1 - 0.5 * SIGMA0 * DELTAT / MU0) / (1 + 0.5 * SIGMA0 * DELTAT / MU0);
			CB(m, n) = (DELTAT / EPS0) / (1 + (SIGMA0 * DELTAT / (2 * EPS0)));
 		}
	}
	for (m = 0; m < SizeX; m++)
	{
		for (n = 0; n < SizeY; n++)
		{
			//Hx(m, n) = 0;
			CP(m, n) = (1 - (SIGMAm0 * DELTAT) / (2 * MU0)) / (1 + (SIGMAm0 * DELTAT) / (2 * MU0));
			CQ(m, n) = (DELTAT / MU0) / (1 + (SIGMAm0 * DELTAT / (2 * MU0)));
		}
	}

	printf("The CA is :%f\n", CA(50, 50));
	printf("The CB is :%f\n", CB(50, 50));
	printf("The CP is :%f\n", CP(50, 50));
	printf("The CQ is :%f\n", CQ(50, 50));
}

double driveSource(Grid *g, double time, double locX, double locY)
{
	double result;
	double ricker_result;
	double ppw = 20;
	double cdtds = CDTDS;
	double w0 = 1000000000.0; //天线的中心频率

	double arg;

	result = time * time * sin(2 * 3.1415 * w0 * time) * exp(-0.93 * w0 * time);

	arg = 3.14159 * ((cdtds * time - 0.0) / ppw - 1.0);
	arg = arg * arg;
	ricker_result = (1.0 - 2.0 * arg) * exp(-arg);

	return ricker_result;
}

void updateE(Grid *g)
{
	int m, n;

	for (m = 1; m < (SizeX - 1); m++)
	{
		for (n = 1; n < (SizeY - 1); n++)
		{
			Ez(m, n) = CA(m, n) * Ez(m, n) + CB(m, n) * ((Hy(m, n) - Hy((m - 1), n)) / DELTA - (Hx(m, n) - Hx(m, (n - 1))) / DELTA);
		}
	}
}

void updateH(Grid *g)
{
	int m, n;

	for (m = 0; m < (SizeX); m++)
	{
		for (n = 0; n < (SizeY - 1); n++)
		{
			Hx(m, n) = CP(m, n) * Hx(m, n) - CQ(m, n) * (Ez(m, (n + 1)) - Ez(m, n)) / DELTA;
		}
	}
	for (m = 0; m < (SizeX - 1); m++)
	{
		for (n = 0; n < (SizeY); n++)
		{
			Hy(m, n) = CP(m, n) * Hy(m, n) + CQ(m, n) * (Ez((m + 1), n) - Ez(m, n)) / DELTA;
		}
	}
}


void updateBoundary(Grid *g)
{
	int m, n;
	int firstX = 0, firstY = 0, lastX = 99, lastY = 99;

	Grid *g1;
	ALLOC_2D(g1, 100, 100, Grid);
	memcpy(g1, g, sizeof(Grid));
	gridInit(g1);
	//updateH(g1);
	updateE(g1);

	//修正Hy左边界
	m = firstX;
	for (n = 1; n < lastY; n++)
	{
		Ez(m, n) = Ez((m + 1), n) + ((c * DELTAT - DELTA) / (c * DELTAT + DELTA)) * (g1->ez[(m+1) * SizeY + n] - Ez(m, n))
				//-((c * c * DELTAT * MU0) / (2 * (c * DELTAT + DELTA))) * (g1->hx[m*SizeY+(n+1)] - g1->hx[m*SizeY + (n - 1)] + g1->hx[(m+1)*SizeY + (n + 1)] - g1->hx[(m + 1)*SizeY + (n - 1)]);
		-((c * c * DELTAT * MU0) / (2 * (c * DELTAT + DELTA))) * (Hx(m, (n + 1)) - Hx(m, (n - 1)) + Hx((m + 1), (n + 1)) - Hx((m + 1), (n - 1)));
	}

	
	//修正Hy右边界
	m = lastX;
	for (n = 1; n < lastY; n++)
	{
		Ez(m, n) = Ez((m - 1), n) + ((c * DELTAT - DELTA) / (c * DELTAT + DELTA)) * (g1->ez[(m - 1) * SizeY + n] - Ez(m, n))
			- ((c * c * DELTAT * MU0) / (2 * (c * DELTAT + DELTA))) * (Hx(m, (n + 1)) - Hx(m, (n - 1)) + Hx((m - 1), (n + 1)) - Hx((m - 1), (n - 1)));
	}
	
	/*
	//修正Hx下边界


	//修正Hx上边界
	*/
}

void writeFile(Grid *g)
{
	int m, n;
	FILE *f = fopen("C:\\Users\\Administrator\\Desktop\\data.txt", "a");
	if (f == NULL) printf("open file error");
	for (n = SizeY - 1; n >= 0; n--)
	{
		for (m = 0; m < SizeX; m++)
		{
			fprintf(f, "%f,", Ez(m, n));
		}
		fprintf(f, "\n");
	}

	fclose(f);
}

然后用matlab对生成的东西成像：

imshow(data50,[]);

时间步为50：

时间步为100：

时间步为150：

会发现150的有反射波了，这是因为没加上吸收边界。

如果加上MUR吸收边界是这样的：

可以看到反射波确实有变好一点，但是边框上出现了黑边，还得对边上的Ez进行修正才行

电磁波时域有限差分方法（FDTD）-C语言-FDTD和激励源和MUR边界

猜你喜欢