2748阅读
9回复

[已解决]Mur二阶边界条件求助

UID ：18905

注册：2008-10-10
登录：2011-04-03
发帖：59
等级：仿真一级

0楼发表于: 2008-10-31 18:51:54

前面编了一阶的,能很好的运行,现在想编一个二阶的,于是换了个公式,但是当波进入边界后,马上不正常了(出现超过极限的反射),下面就是使用的二阶Mur的边界条件(按葛书编的),望大家指教.
% 3.二阶Mur边界条件
k0=8.854e-12;%真空介电常数
h1=(c0*dt-ds)/(c0*dt+ds);
h2=c0^2*k0*dt/2/(c0*dt+ds);

for i=1
    for j=2:JE-1
        Hz(i,j)=Hz(i+1,j)+h1*(Hz(i+1,j)-Hz(i,j))+h2*(Ex(i,j)-Ex(i,j-1)+Ex(i+1,j)-Ex(i+1,j-1));%左侧边界
    end
end
for i=IE
    for j=2:JE-1
        Hz(i,j)=Hz(i-1,j)+h1*(Hz(i-1,j)-Hz(i,j))+h2*(Ex(i,j)-Ex(i,j-1)+Ex(i-1,j)-Ex(i-1,j-1));%右侧边界
    end
end
for j=1
    for i=2:IE-1
        Hz(i,j)=Hz(i,j+1)+h1*(Hz(i,j+1)-Hz(i,j))+h2*(Ey(i,j)-Ey(i-1,j)+Ey(i,j+1)-Ey(i-1,j+1));%下侧边界
    end
end
for j=JE
    for i=2:IE-1
        Hz(i,j)=Hz(i,j-1)+h1*(Hz(i,j-1)-Hz(i,j))+h2*(Ey(i,j)-Ey(i-1,j)+Ey(i,j-1)-Ey(i-1,j-1));%上侧边界
& ..

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

离线gwzhao

方恨少

UID ：17098

注册：2008-08-24
登录：2019-01-09
发帖：1374
等级：荣誉管理员

1楼发表于: 2008-10-31 21:31:00

好像公式不是这样的吧。

逆流而上

离线aoso

<a href="http://b

UID ：4796

注册：2007-09-04
登录：2025-06-23
发帖：592
等级：八级仿真大师

2楼发表于: 2008-10-31 21:36:18

找本葛德彪老师那本书后面有个2维2阶Mur的例子，fortran编的，和matlab语法差不多

离线xiaoweike

UID ：18905

注册：2008-10-10
登录：2011-04-03
发帖：59
等级：仿真一级

3楼发表于: 2008-10-31 21:48:44

能不能给我一下matlab的二阶Mur边界条件,我手头上现在没书.

离线cem-uestc

UID ：9061

注册：2008-03-07
登录：2019-01-05
发帖：2575
等级：荣誉管理员

4楼发表于: 2008-10-31 21:49:30

Mur2阶的公式有4个系数，公式有问题

欢迎光临
http://www.mwtee.com/home.php?mod=space&uid=13535

离线gwzhao

方恨少

UID ：17098

注册：2008-08-24
登录：2019-01-09
发帖：1374
等级：荣誉管理员

5楼发表于: 2008-10-31 22:57:08

% 本程序实现2维TM波FDTD仿真
% 此程序用PML设置吸收边界条件
% FDTD_2D_PML
% 仅含有Ez,Hx,Hy分量

clear;
clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 1.初始化
T=200; % 迭代次数
IE=100; %
JE=100;
npml=8; % PML的网格数量

c0=3*10^8; % 波速
f=1.5*10^(9); % 频率
lambda=c0/f; % 波长

wl=10; %数值波长
dx=lambda/wl;
dy=lambda/wl;
pi=3.14159;

dt=dx/(2*c0); % 时间间隔

epsz=1/(4*pi*9*10^9); % 真空介电常数
epsilon=1; % 相对介电常数
sigma=0; % 电导率

spread=6; % 脉冲宽度
t0=20; % 脉冲高度

ic=IE/2; % 源的X位置
jc=JE/2; % 源的Y位置
%初始化
for i=1E+1;
for j=1:JE+1;
dz(i,j)=0; % z方向电荷密度
ez(i,j)=0; % z方向电场
hx(i,j)=0; % x方向磁场
hy(i,j)=0; % y方向磁场
ihx(i,j)=0;%
ihy(i,j)=0;
iz(i,j)=0; % z方向求和参量，频域卷积转化为时域求和
end;
end;

for i=2E; %
for j=2:JE;
ga(i,j)=1;
end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%PML参数的设置
for i=1E;
gi2(i)=1;
gi3(i)=1;
fi1(i)=0;
fi2(i)=1.0;
fi3(i)=1.0;
end

for j=1:JE;
gj2(j)=1;
gj3(j)=1;
fj1(j)=0;
fj2(j)=1;
fj3(j)=1;
end

for i=1:npml+1; %设置PML层中的参数
xnum=npml+1-i;
xn=0.33*(xnum/npml)^3;
gi2(i)=1.0/(1+xn);
gi2(IE-1-i)=1/(1+xn);
gi3(i)=(1-xn)/(1+xn);
gi3(IE-1-i)=(1-xn)/(1+xn);
xn=0.25*((xnum-0.5)/npml)^3;
fi1(i)=xn;
fi1(IE-2-i)=xn;
fi2(i)=1.0/(1+xn);
fi2(IE-2-i)=1/(1+xn);
fi3(i)=(1-xn)/(1+xn);
fi3(IE-2-i)=(1-xn)/(1+xn);
end

for i=1:npml+1;
xnum=npml+1-i;
xn=0.33*(xnum/npml)^3;
gj2(i)=1.0/(1+xn);
gj2(JE-1-i)=1/(1+xn);
gj3(i)=(1-xn)/(1+xn);
gj3(JE-1-i)=(1-xn)/(1+xn);
xn=0.25*((xnum-0.5)/npml)^3;
fj1(i)=xn;
fj1(JE-2-i)=xn;
fj2(i)=1.0/(1+xn);
fj2(JE-2-i)=1/(1+xn);
fj3(i)=(1-xn)/(1+xn);
fj3(JE-2-i)=(1-xn)/(1+xn);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 2.迭代求解电场和磁场
for t=1:T;
for i=2E; % 为了使每个电场周围都有磁场进行数组下标处理
for j=2:JE;
dz(i,j)=gi3(i)*gj3(j)*dz(i,j)+gi2(i)*gj2(j)*0.5*(hy(i,j)-hy(i-1,j)-hx(i,j)+hx(i,j-1));
end;
end; % 电场循环结束

pulse=sin(2*pi*f*t*dt); % 正弦波源
dz(ic,jc)=dz(ic,jc)+pulse; % 软源

for i=1E; % 为了使每个电场周围都有磁场进行数组下标处理
for j=1:JE;
ez(i,j)=ga(i,j)* dz(i,j); %反映煤质的情况都是放到这里的
% iz(i,j)=iz(i,j)+gb(i,j)*ez(i,j) ;
end;
end; % 电荷密度循环结束

for j=1:JE;
ez(1,j)=0;
ez(IE,j)=0;
end

for i=1E;
ez(i,1)=0;
ez(i,JE)=0;
end;

for i=1E; % 为了使每个磁场周围都有电场进行数组下标处理
for j=1:JE-1;
curl_e=ez(i,j)-ez(i,j+1);
ihx(i,j)=ihx(i,j)+fi1(i)*curl_e;
hx(i,j)=fj3(j)*hx(i,j)+fj2(j)*0.5*(curl_e+ihx(i,j));
end;
end; % 磁场HX循环结束

for i=1E-1; % 为了使每个磁场周围都有电场进行数组下标处理
for j=1:JE;
curl_e=ez(i+1,j)-ez(i,j);
ihy(i,j)=ihy(i,j)+fj1(j)*curl_e;
hy(i,j)=fi3(i)*hy(i,j)+fi2(i)*0.5*(curl_e+ihy(i,j));
end;
end; % 磁场HY循环结束
end;
end;

逆流而上

离线xiaoweike

UID ：18905