2031阅读
0回复

[求助]FDFD 频域有闲差分2D散射问题

alizwell

UID ：2851

注册：2007-05-22
登录：2017-09-27
发帖：110
等级：仿真二级

0楼发表于: 2016-10-15 20:41:34

大家好，有懂FDFD方法的朋友吗？最近我看了网上的一份教程，自己动手编写了一个求解2纬自由空间Hz模式传播的代码。虽然加上了PML和周期边界条件，可是有反射，驻波很明显。能否请达人指教？代码有三部分组成，
FDFD_2D_freespace.m是主文件，需要调用两个函数（无法用附件故把原代码写在下面）
yeeder.m 生成基于yee网个的差分算子，
calcpml2d.m生成PML边界条件

请指正
下图是计算结果，可以看到明显的驻波，SF为散射场区域，按理说应该没有场才对。

%% This block tests the 2D FDFD to visualise if standing wave exists
% positive convention must be applied because of the DEX DEY DHX DHY uses
% positive convention to account periodic boudnary conditions.
% d/dz=0; H mode with Hz Ex Ey and E mode with Ez Hx Hy
clear all
close all
clc
Degree=pi/180;
%=========== Source Parameters ===========
c0=299792458;
Theta=15*Degree;
f0=8e9; % 8GHz design frequency
f1=20e9;    % 10GHz source frequency
Lam0=c0/f0; % design wavelength
Lam1=c0/f1; % source frequency
k0=2*pi/Lam1;   % wave vector determined by the source
kx=k0*cos(Theta);
ky=k0*sin(Theta);
Mode='H';   % H mode and E mode leads to different routines
% ==============GRID PARAMETERS==============
NRES=20;    % mask for reflected wave
BUFX=1*Lam0*[1 1];  % spacer region above and below grating excluding PML
NPML=[30 30 0 0];   % PML xlo xhi, but no PML ylo yhi
%========== dx dy Nx Ny ===============
dx=min(Lam1/20, Lam0/20);
dy=min(Lam1/20, Lam0/20);
Nx=round(sum(BUFX)/dx)+NPML(1)+NPML(2);
Ny=31;   % to make odd number
% ==========================================================
% Build the Device on 2X Grids based on Yee grids
%                   2X Grids
% the dielectric assumed isotropic
% =========================================================
Nx2=2*Nx;
Ny2=2*Ny;
UR2=ones(Nx2,Ny2);  % permeability is 1 EVERYWHERE
ER2=ones(Nx2,Ny2);  % permitivity is 1 EVERYWHERE

%================Parse to 1X Grid=================
URxx=UR2(1:2:Nx2,2:2:Ny2);
URyy=UR2(2:2:Nx2,1:2:Ny2);
URzz=UR2(2:2:Nx2,2:2:Ny2);
ERxx=ER2(2:2:Nx2,1:2:Ny2);
ERyy=ER2(1:2:Nx2,2:2:Ny2);
ERzz=ER2(1:2:Nx2,1:2:Ny2);
%=============PML APPLIED TO X  =================
[sx,sy]=calcpml2d([Nx, Ny], NPML);   % 20 cell on each boudary
ERxx=ERxx.*sy./sx;
ERyy=ERyy.*sx./sy;
ERzz=ERzz.*sx.*sy;
URxx=URxx.*sy./sx;
URyy=URyy.*sx./sy;
URzz=URzz.*sx.*sy;
%==============Form Sparse Matrix      ====================
ERxx=spdiags(reshape(ERxx,Nx*Ny,1),0,Nx*Ny,Nx*Ny);
ERyy=spdiags(reshape(ERyy,Nx*Ny,1),0,Nx*Ny,Nx*Ny);
ERzz=spdiags(reshape(ERzz,Nx*Ny,1),0,Nx*Ny,Nx*Ny);
URxx=spdiags(reshape(URxx,Nx*Ny,1),0,Nx*Ny,Nx*Ny);
URyy=spdiags(reshape(URyy,Nx*Ny,1),0,Nx*Ny,Nx*Ny);
URzz=spdiags(reshape(URzz,Nx*Ny,1),0,Nx*Ny,Nx*Ny);
%===== Form the Normalised Curl Operators ============
[DEX,DEY,DHX,DHY]=yeeder([Nx Ny], [dx*k0 dy*k0], [0 -2], [kx/k0, ky/k0]);    % normalised

%==============Build Scattering Problem ===========
%   form sorce field as a column vector b
%   Mask Matrix for TF/SF regions
%   Solve the linear function
%==================================================
%=== Source field plane wave=====
x=dx.*(0:Nx-1);
y=dy.*(0:Ny-1);
[X,Y]=meshgrid(x,y);
X=X.';
Y=Y.';  % meshgrid transpose X to horizontal axis so must be corrected
fsrc=exp(1j.*kx.*X+1j.*ky.*Y);  % positive convention
b=reshape(fsrc,Nx*Ny,1); % Keep exciation as a collomn vector

% ========== Mask Q TF/SF Interface ============
Q=ones(Nx, Ny);
Q(NPML(1)+NRES:Nx,:)=0;
Q=diag(sparse(Q(:)));

% ========== solve linear function =================
if Mode=='H'
  &nb ..

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

描述:SF 为散射场区域，TF为总场区域

图片:abs_Hz.jpg

描述:PML
附件:

calcpml2d.txt (3 K) 下载次数:1

描述:YEE 差分算子
附件:

yeeder.txt (3 K) 下载次数:1

专业是电磁兼容测试和建模

发帖回复