1843阅读
5回复

[RFEDA原创]FDTD迷你代码-2D CPML

UID ：90389

注册：2012-03-16
登录：2014-09-11
发帖：17
等级：仿真新人

0楼发表于: 2014-06-11 09:00:31

clear;
clc;
e0=1e-9/(36*pi);%epsilon0
m0=1e-7*(4*pi);%mu0
c=3e8;%light speed
f=1e9;
lambda=c/f;%wavelength
omega=2*pi*f;%omega
nd=20;%space grid ratio to lambda
d=lambda/nd;%grid size
s=0.5;%courantsteadycoefficient
dt=s*d/c;
nx=120;%grid x number
ny=120;%grid y number
na=10;%absorber layers thickness
ez=zeros(nx+1,ny+1);
hx=zeros(nx-1,ny);
hy=zeros(nx,ny-1);
psixHy=zeros(nx-1,ny-1);
psiyHx=zeros(nx-1,ny-1);
psixEz=zeros(nx,ny-1);
psiyEz=zeros(nx-1,ny);
kappaxHy=ones(nx-1,ny-1);
kappayHx=ones(nx-1,ny-1);
kappaxEz=ones(nx,ny-1);
kappayEz=ones(nx-1,ny);
expAlphaxHy=ones(nx-1,ny-1);
expAlphayHx=ones(nx-1,ny-1);
expAlphaxEz=ones(nx,ny-1);
expAlphayEz=ones(nx-1,ny);
cxHy=zeros(nx-1,ny-1);
cyHx=zeros(nx-1,ny-1);
cxEz=zeros(nx,ny-1);
cyEz=zeros(nx-1,ny);
m=2;
sigmaMax=(m+1)/(150*pi*d);%lost tangent max
kappaMax=1;
aMax=0.05;
sigma=sigmaMax*((na-1:-1:0)/(na-1)).^m;
kappa=1+(kappaMax-1)*((na-1:-1:0)/(na-1)).^m;
a=aMax*((0:na-1)/(na-1));
cb=dt/(e0*d);
cq=dt/(m0*d);
kappaH=kappa;
alphaH=sigma/e0./kappa+a/e0;
expAlphaH=exp(-alphaH*dt);
cH=sigma/e0./kappa.^2./alphaH.*(expAlphaH-1);
kappaE=(kappaH(1:end-1)+kappaH(2:end))/2;
alphaE=(alphaH(1:end-1)+alphaH(2:end))/2;
expAlphaE=(expAlphaH(1:end-1)+expAlphaH(2:end))/2;
cE=(cH(1:end-1)+cH(2:end))/2;
kappaxHy(1:na-1,:)=repmat(kappaE',1,ny-1);
kappaxHy(end-na+2:end,:)=repmat(fliplr(kappaE)',1,ny-1);
kappayHx(:,1:na-1)=repmat(kappaE,nx-1,1);
kappayHx(:,end-na+2:end)=repmat(fliplr(kappaE),nx-1,1);
kappaxEz(1:na,:)=repmat(kappaH',1,ny-1);
kappaxEz(end-na+1:end,:)=repmat(fliplr(kappaH)',1,ny-1);
kappayEz(:,1:na)=repmat(kappaH,nx-1,1);
kappayEz(:,end-na+1:end)=repmat(fliplr(kappaH),nx-1,1);
expAlphaxHy(1:na-1,:)=repmat(expAlphaE',1,ny-1);
expAlphaxHy(end-na+2:end,:)=repmat(fliplr(expAlphaE)',1,ny-1);
expAlphayHx(:,1:na-1)=repmat(expAlphaE,nx-1,1);
expAlphayHx(:,end-na+2:end)=repmat(fliplr(expAlphaE),nx-1,1);
expAlphaxEz(1:na,:)=repmat(expAlphaH',1,ny-1);
expAlphaxEz(end-na+1:end,:)=repmat(fliplr(expAlphaH)',1,ny-1);
expAlphayEz(:,1:na)=repmat(expAlphaH,nx-1,1);
expAlphayEz(:,end-na+1:end)=repmat(fliplr(expAlphaH),nx-1,1);
cxHy(1:na-1,:)=repmat(cE',1,ny-1);
cxHy(end-na+2:end,:)=repmat(fliplr(cE)',1,ny-1);
cyHx(:,1:na-1)=repmat(cE,nx-1,1);
cyHx(:,end-na+2:end)=repmat(fliplr(cE),nx-1,1);
cxEz(1:na,:)=repmat(cH',1,ny-1);
cxEz(end-na+1:end,:)=repmat(fliplr(cH)',1,ny-1);
cyEz(:,1:na)=repmat(cH,nx-1,1);
cyEz(:,end-na+1:end)=repmat(fliplr(cH),nx-1,1);
h=imagesc(ez,0.05*[-1,1]);%handle for image
axis equal tight;
for n=1:500
    ez(nx/2+1,ny/2+1)=ez(nx/2+1,ny/2+1)+sin(omega*n*dt);
    psixHy=expAlphaxHy.*psixHy+cxHy.*diff(hy,1,1);
    psiyHx=expAlphayHx.*psiyHx+cyHx.*diff(hx,1,2);
    ez(2:end-1,2:end-1)=ez(2:end-1,2:end-1)+cb*(diff(hy,1,1)./kappaxHy+psixHy-diff(hx,1,2)./kappayHx-psiyHx);
    psiyEz=expAlphayEz.*psiyEz+cyEz.*diff(ez(2:end ..

未注册仅能浏览部分内容，查看全部内容及附件请先登录或注册

描述:请输入描述

图片:untitled.jpg

共1条评分，rf币+5

chenjiabao1989

rf币 +5

十分期待您分享下一贴！

2014-12-12

离线八云

所有的懦弱皆源自于没有爱或爱得不够彻底，这两者一样。

UID ：122

注册：2006-10-30
登录：2024-06-13
发帖：1971
等级：七级仿真大师

1楼发表于: 2014-12-12 00:52:07

同学，你的PML基础知识来源于什么文献啊？能共享一下吗？

所谓的大学，是为了人类的进步，与世俗权力和常识斗争的崇高之地。
无数科学家牺牲自己与宗教、权力斗争，才有了今天。

离线zly0427

UID ：118577

注册：2015-10-28
登录：2017-12-07
发帖：151
等级：仿真二级

2楼发表于: 2016-03-06 17:41:16

大神厉害啊

离线zly0427

UID ：118577

注册：2015-10-28
登录：2017-12-07
发帖：151
等级：仿真二级

3楼发表于: 2016-03-30 11:33:43

回 zly0427 的帖子

zly0427:大神厉害啊 (2016-03-06 17:41) V=dOeuYd

你下载到这个文件了么我的为什么看不到
912028709qq.com 可以发给我么

离线zly0427

UID ：118577

注册：2015-10-28
登录：2017-12-07
发帖：151
等级：仿真二级

4楼发表于: 2016-03-30 11:34:04

回八云的帖子

八云:同学，你的PML基础知识来源于什么文献啊？能共享一下吗？ (2014-12-12 00:52) .,xyE--;d

你下载到这个文件了么我的为什么看不到
912028709qq.com 可以发给我么

离线zly0427

UID ：118577

注册：2015-10-28
登录：2017-12-07
发帖：151
等级：仿真二级

5楼发表于: 2016-03-30 11:34:38

你下载到这个文件了么我的为什么看不到
912028709qq.com 可以发给我么

发帖回复