5930阅读
31回复

[求助]3d-UPML程序,大家都来看看顺便给我解答一下疑问(有些疑问解决的地方我给出了注释)

UID ：20925

注册：2008-11-06
登录：2021-04-22
发帖：227
等级：仿真三级

0楼发表于: 2009-01-02 16:24:53

   学习这个程序有短时间了，得到了论坛上不少朋友的指点感谢herosword师兄和gwzhao版主的帮助，还有网上一些朋友帮助，现在我把一些别人的讲解在程序后面注释，当然还有一些细节问题需要大家来讨论，帮助解决。
%***********************************************************************
%     3-D FDTD code with UPML absorbing boundary conditions
%***********************************************************************
%
%     Program author: Keely J. Willis, Graduate Student
%                     UW Computational Electromagnetics Laboratory
%                           Director: Susan C. Hagness
%                     Department of Electrical and Computer Engineering
%                     University of Wisconsin-Madison
%                     1415 Engineering Drive
%                     Madison, WI 53706-1691
%                     [url=mailto:kjwillis@wisc.edu]kjwillis@wisc.edu[/url]
%
%     Copyright 2005
%
%     This MATLAB M-file implements the finite-difference time-domain
%     solution of Maxwell's curl equations over a three-dimensional
%     Cartesian space lattice comprised of uniform cubic grid cells.
%
%     The dimensions of the computational domain are 8.2 cm
%     (x-direction), 3.4 cm (y-direction), and 3.2 cm (z-direction).
%     The grid is terminated with UPML absorbing boundary conditions.
%
%     An electric current source comprised of two collinear Jz components
%     (realizing a Hertzian dipole) excites a radially propagating wave.
%     The current source is located in the center of the grid. The
%     source waveform is a differentiated Gaussian pulse given by
%          J(t)=J0*(t-t0)*exp(-(t-t0)^2/tau^2),
%     where tau=50 ps. The FWHM spectral bandwidth of this zero-dc-
%     content pulse is approximately 7 GHz. The grid resolution
%     (dx = 2 mm) was chosen to provide at least 10 samples per
%     wavelength up through 15 GHz.
%
%     To execute this M-file, type "fdtd3D_UPML" at the MATLAB prompt.
%
%     This code has been tested in the following Matlab environments:
%     Matlab version 6.1.0.450 Release 12.1 (May 18, 2001)
%     Matlab version 6.5.1.199709 Release 13 Service Pack 1 (August 4, 2003)
%     Matlab version 7.0.0.19920 R14 (May 6, 2004)
%     Matlab version 7.0.1.24704 R14 Service Pack 1 (September 13, 2004)
%     Matlab version 7.0.4.365 R14 Service Pack 2 (January 29, 2005)
%
%     Note: if you are using Matlab version 6.x, you may wish to make
%     one or more of the following modifications to this code:
%       --uncomment line numbers 485 and 486
%       --comment out line numbers 552 and 561
%
%***********************************************************************
clear
%***********************************************************************
%     Fundamental constants
%***********************************************************************
cc=2.99792458e8;
muz=4.0*pi*1.0e-7;
epsz=1.0/(cc*cc*muz);
etaz=sqrt(muz/epsz);%真空中波阻抗匹配
%***********************************************************************
%     Material parameters
%***********************************************************************
mur=1.0;
epsr=1.0;
eta=etaz*sqrt(mur/epsr);
%***********************************************************************
%     Grid parameters
%
%     Each grid size variable name describes the number of sampled points
%     for a particular field component in the direction of that component.
%     Additionally, the variable names indicate the region of the grid
%     for which the dimension is relevant. For example, ie_tot is the
%     number of sample points of Ex along the x-axis in the total
%     computational grid, and jh_bc is the number of sample points of Hy
%     along the y-axis in the y-normal UPML regions.
%
%***********************************************************************
ie=41;          % Size of main grid
je=17;
ke=16;
ih=ie+1;
jh=je+1;
kh=ke+1;
upml=10;        % Thickness of PML boundaries
ih_bc=upml+1;
jh_bc=upml+1;
kh_bc=upml+1;
ie_tot=ie+2*upml;          % Size of total computational domain
je_tot=je+2*upml;
ke_tot=ke+2*upml;
ih_tot=ie_tot+1;
jh_tot=je_tot+1;
kh_tot=ke_tot+1;
is=round(ih_tot/2);         % Location of z-directed current source
js=round(jh_tot/2);
ks=round(ke_tot/2);
%***********************************************************************
%     Fundamental grid parameters
%***********************************************************************
delta=0.002;
dt=delta*sqrt(epsr*mur)/(2.0*cc);%?
nmax=200;
%***********************************************************************
%     Differentiated Gaussian pulse excitation
%***********************************************************************
rtau=50.0e-12;
tau=rtau/dt;
ndelay=3*tau;%这个延迟我个人认为是电磁波在场中传播达到时谐状态所要的时间
J0=-1.0*epsz;%
%***********************************************************************
%     Initialize field arrays
%***********************************************************************
ex=zeros(ie_tot,jh_tot,kh_tot); ％这个地方电场量的初始赋值，X方向比y,z方向多一个（这个地方我还不是很理解）后面的一次类推。
ey=zeros(ih_tot,je_tot,kh_tot);
ez=zeros(ih_tot,jh_tot,ke_tot);
dx=zeros(ie_tot,jh_tot,kh_tot);
dy=zeros(ih_tot,je_tot,kh_tot);
dz=zeros(ih_tot,jh_tot,ke_tot);
hx=zeros(ih_tot,je_tot,ke_tot);

hy=zeros(ie_tot,jh_tot,ke_tot);
hz=zeros(ie_tot,je_tot,kh_tot);
bx=zeros(ih_tot,je_tot,ke_tot);
by=zeros(ie_tot,jh_tot,ke_tot);
bz=zeros(ie_tot,je_tot,kh_tot);
%***********************************************************************
%     Initialize update coefficient arrays
%***********************************************************************
C1ex=zeros(size(ex)); ％迭代系数赋初始值
C2ex=zeros(size(ex));
C3ex=zeros(size(ex));
C4ex=zeros(size(ex));
C5ex=zeros(size(ex));
C6ex=zeros(size(ex));
C1ey=zeros(size(ey));
C2ey=zeros(size(ey));
C3ey=zeros(size(ey));
C4ey=zeros(size(ey));
C5ey=zeros(size(ey));
C6ey=zeros(size(ey));
C1ez=zeros(size(ez));
C2ez=zeros(size(ez));
C3ez=zeros(size(ez));
C4ez=zeros(size(ez));
C5ez=zeros(size(ez));
C6ez=zeros(size(ez));
D1hx=zeros(size(hx));
D2hx=zeros(size(hx));
D3hx=zeros(size(hx));
D4hx=zeros(size(hx));
D5hx=zeros(size(hx));
D6hx=zeros(size(hx));
D1hy=zeros(size(hy));
D2hy=zeros(size(hy));
D3hy=zeros(size(hy));
D4hy=zeros(size(hy));
D5hy=zeros(size(hy));
D6hy=zeros(size(hy));
D1hz=zeros(size(hz));
D2hz=zeros(size(hz));
D3hz=zeros(size(hz));
D4hz=zeros(size(hz));
D5hz=zeros(size(hz));
D6hz=zeros(size(hz));
%***********************************************************************
%     Update coefficients, as described in Section 7.8.2.
%
%     In order to simplify the update equations used in the time-stepping
%     loop, we implement UPML update equations throughout the entire
%     grid. In the main grid, the electric-field update coefficients of
%     Equations 7.91a-f and the correponding magnetic field update
%     coefficients extracted from Equations 7.89 and 7.90 are simplified
%     for the main grid (free space) and calculated below.
%
%***********************************************************************
C1=1.0; ％自由空间的电导率为0带入公式可以得到，下面的系数也是如此
C2=dt;
C3=1.0;
C4=1.0/2.0/epsr/epsr/epsz/epsz;
C5=2.0*epsr*epsz;
C6=2.0*epsr*epsz;
D1=1.0;
D2=dt;
D3=1.0;
D4=1.0/2.0/epsr/epsz/mur/muz;
D5=2.0*epsr*epsz;
D6=2.0*epsr*epsz;
%***********************************************************************
%     Initialize main grid update coefficients
%***********************************************************************
C1ex(:,jh_bc:jh_tot-upml,:)=C1; ％总场区域的系数定义
C2ex(:,jh_bc:jh_tot-upml,:)=C2;
C3ex(:,:,kh_bc:kh_tot-upml)=C3;
C4ex(:,:,kh_bc:kh_tot-upml)=C4;
C5ex(ih_bc:ie_tot-upml,:,:)=C5;
C6ex(ih_bc:ie_tot-upml,:,:)=C6;
C1ey(:,:,kh_bc:kh_tot-upml)=C1;
C2ey(:,:,kh_bc:kh_tot-upml)=C2;
C3ey(ih_bc:ih_tot-upml,:,:)=C3;
C4ey(ih_bc:ih_tot-upml,:,:)=C4;
C5ey(:,jh_bc:je_tot-upml,:)=C5;
C6ey(:,jh_bc:je_tot-upml,:)=C6;
C1ez(ih_bc:ih_tot-upml,:,:)=C1;
C2ez(ih_bc:ih_tot-upml,:,:)=C2;
C3ez(:,jh_bc:jh_tot-upml,:)=C3;
C4ez(:,jh_bc:jh_tot-upml,:)=C4;
C5ez(:,:,kh_bc:ke_tot-upml)=C5;
C6ez(:,:,kh_bc:ke_tot-upml)=C6;
D1hx(:,jh_bc:je_tot-upml,:)=D1;
D2hx(:,jh_bc:je_tot-upml,:)=D2;
D3hx(:,:,kh_bc:ke_tot-upml)=D3;
D4hx(:,:,kh_bc:ke_tot-upml)=D4;
D5hx(ih_bc:ih_tot-upml,:,:)=D5;
D6hx(ih_bc:ih_tot-upml,:,:)=D6;
D1hy(:,:,kh_bc:ke_tot-upml)=D1;
D2hy(:,:,kh_bc:ke_tot-upml)=D2;
D3hy(ih_bc:ie_tot-upml,:,:)=D3;
D4hy(ih_bc:ie_tot-upml,:,:)=D4;
D5hy(:,jh_bc:jh_tot-upml,:)=D5;
D6hy(:,jh_bc:jh_tot-upml,:)=D6;
D1hz(ih_bc:ie_tot-upml,:,:)=D1;
D2hz(ih_bc:ie_tot-upml,:,:)=D2;
D3hz(:,jh_bc:je_tot-upml,:)=D3;
D4hz(:,jh_bc:je_tot-upml,:)=D4;
D5hz(:,:,kh_bc:kh_tot-upml)=D5;
D6hz(:,:,kh_bc:kh_tot-upml)=D6;
%***********************************************************************
%     Fill in PML regions
%
%     PML theory describes a continuous grading of the material properties
%     over the PML region. In the FDTD grid it is necessary to discretize
%     the grading by averaging the material properties over a grid cell
%     width centered on each field component. As an example of the
%     implementation of this averaging, we take the integral of the
%     continuous sigma(x) in the PML region
%
%         sigma_i = integral(sigma(x))/delta
%
%     where the integral is over a single grid cell width in x, and is
%     bounded by x1 and x2. Applying this to the polynomial grading of
%     Equation 7.60a produces
%
%         sigma_i = (x2^(m+1)-x1^(m+1))*sigmam/(delta*(m+1)*d^m)
%
%     where sigmam is the maximum value of sigma as described by Equation
%     7.62.
%
%     The definitions of x1 and x2 depend on the position of the field
%     component within the grid cell. We have either
%
%         x1 = (i-0.5)*delta, x2 = (i+0.5)*delta
%
%     or
%
%         x1 = (i)*delta,      x2 = (i+1)*delta
%
%     where i varies over the PML region.
%
%***********************************************************************
rmax=exp(-16); %desired reflection error, designated as R(0) in Equation 7.62
orderbc=4;      %order of the polynomial grading, designated as m in Equation 7.60a,b
%   x-varying material properties
delbc=upml*delta; ％PML的厚度
sigmam=-log(rmax)*(orderbc+1.0)/(2.0*eta*delbc); ％参见Taflove书的7。62式
sigfactor=sigmam/(delta*(delbc^orderbc)*(orderbc+1.0));
kmax=1; ％这个地方不是很理解
kfactor=(kmax-1.0)/delta/(orderbc+1.0)/delbc^orderbc; ％这个式子应该是Taflove书的7.60(b)式。
for i=1:upml

    % Coefficients for field components in the center of the grid cell
   ％x方向我们只考虑跟x有关的系数，把迭代公式的系数中含有x的系数挑出：c5ex,c6ex d5hx,d6hx,c3ey,c4ey,c1ez,c2ez,d1hz,d2hz,d3hy,d4hy。根据Taflove书中Yee网格的形式，从X方向看去ex,hy,hz是在一个平面的上，而ey,ez,hx与ex，hy,hz相差半个网格。所以这样可以分网格中心和网格边界2部分来定义。同样Y和Z方向是一样考虑的
    x1=(upml-i+1)*delta;
    x2=(upml-i)*delta;
    sigma=sigfactor*(x1^(orderbc+1)-x2^(orderbc+1));
    ki=1+kfactor*(x1^(orderbc+1)-x2^(orderbc+1));
    facm=(2*epsr*epsz*ki-sigma*dt);
    facp=(2*epsr*epsz*ki+sigma*dt);
    C5ex(i,:,:)=facp;
    C5ex(ie_tot-i+1,:,:)=facp;
    C6ex(i,:,:)=facm;
    C6ex(ie_tot-i+1,:,:)=facm;
    D1hz(i,:,:)=facm/facp;
    D1hz(ie_tot-i+1,:,:)=facm/facp;
    D2hz(i,:,:)=2.0*epsr*epsz*dt/facp;
    D2hz(ie_tot-i+1,:,:)=2.0*epsr*epsz*dt/facp;
    D3hy(i,:,:)=facm/facp;
    D3hy(ie_tot-i+1,:,:)=facm/facp;
    D4hy(i,:,:)=1.0/facp/mur/muz;
    D4hy(ie_tot-i+1,:,:)=1.0/facp/mur/muz;
    % Coefficients for field components on the grid cell boundary
    x1=(upml-i+1.5)*delta;
    x2=(upml-i+0.5)*delta;
    sigma=sigfactor*(x1^(orderbc+1)-x2^(orderbc+1));
    ki=1.0+kfactor*(x1^(orderbc+1)-x2^(orderbc+1));
    facm=(2.0*epsr*epsz*ki-sigma*dt);
    facp=(2.0*epsr*epsz*ki+sigma*dt);
    C1ez(i,:,:)=facm/facp;
    C1ez(ih_tot-i+1,:,:)=facm/facp;
    C2ez(i,:,:)=2.0*epsr*epsz*dt/facp;
    C2ez(ih_tot-i+1,:,:)=2.0*epsr*epsz*dt/facp;
    C3ey(i,:,:)=facm/facp;
    C3ey(ih_tot-i+1,:,:)=facm/facp;
    C4ey(i,:,:)=1.0/facp/epsr/epsz;
    C4ey(ih_tot-i+1,:,:)=1.0/facp/epsr/epsz;
    D5hx(i,:,:)=facp;
    D5hx(ih_tot-i+1,:,:)=facp;
    D6hx(i,:,:)=facm;
    D6hx(ih_tot-i+1,:,:)=facm;

end
%   PEC walls
％这个还需高人来解释，因为我解释不清楚，希望有人来补充
C1ez(1,:,:)=-1.0;
C1ez(ih_tot,:,:)=-1.0;
C2ez(1,:,:)=0.0;
C2ez(ih_tot,:,:)=0.0;
C3ey(1,:,:)=-1.0;
C3ey(ih_tot,:,:)=-1.0;
C4ey(1,:,:)=0.0;
C4ey(ih_tot,:,:)=0.0;
%   y-varying material properties
delbc=upml*delta;
sigmam=-log(rmax)*epsr*epsz*cc*(orderbc+1.0)/(2.0*delbc);
sigfactor=sigmam/(delta*(delbc^orderbc)*(orderbc+1.0));
kmax=1.0;
kfactor=(kmax-1.0)/delta/(orderbc+1.0)/delbc^orderbc;
for j=1:upml

    % Coefficients for field components in the center of the grid cell
    y1=(upml-j+1)*delta;
    y2=(upml-j)*delta;
    sigma=sigfactor*(y1^(orderbc+1)-y2^(orderbc+1));
    ki=1+kfactor*(y1^(orderbc+1)-y2^(orderbc+1));
    facm=(2*epsr*epsz*ki-sigma*dt);
    facp=(2*epsr*epsz*ki+sigma*dt);

    C5ey(:,j,:)=facp;
    C5ey(:,je_tot-j+1,:)=facp;
    C6ey(:,j,:)=facm;
    C6ey(:,je_tot-j+1,:)=facm;
    D1hx(:,j,:)=facm/facp;
    D1hx(:,je_tot-j+1,:)=facm/facp;
    D2hx(:,j,:)=2*epsr*epsz*dt/facp;
    D2hx(:,je_tot-j+1,:)=2*epsr*epsz*dt/facp;
    D3hz(:,j,:)=facm/facp;
    D3hz(:,je_tot-j+1,:)=facm/facp;
    D4hz(:,j,:)=1/facp/mur/muz;
    D4hz(:,je_tot-j+1,:)=1/facp/mur/muz;

    % Coefficients for field components on the grid cell boundary
    y1=(upml-j+1.5)*delta;
    y2=(upml-j+0.5)*delta;
    sigma=sigfactor*(y1^(orderbc+1)-y2^(orderbc+1));
    ki=1+kfactor*(y1^(orderbc+1)-y2^(orderbc+1));
    facm=(2*epsr*epsz*ki-sigma*dt);
    facp=(2*epsr*epsz*ki+sigma*dt);

    C1ex(:,j,:)=facm/facp;
    C1ex(:,jh_tot-j+1,:)=facm/facp;
    C2ex(:,j,:)=2*epsr*epsz*dt/facp;
    C2ex(:,jh_tot-j+1,:)=2*epsr*epsz*dt/facp;
    C3ez(:,j,:)=facm/facp;
    C3ez(:,jh_tot-j+1,:)=facm/facp;
    C4ez(:,j,:)=1/facp/epsr/epsz;
    C4ez(:,jh_tot-j+1,:)=1/facp/epsr/epsz;
    D5hy(:,j,:)=facp;
    D5hy(:,jh_tot-j+1,:)=facp;
    D6hy(:,j,:)=facm;
    D6hy(:,jh_tot-j+1,:)=facm;
end
%   PEC walls
C1ex(:,1,:)=-1;
C1ex(:,jh_tot,:)=-1;
C2ex(:,1,:)=0;
C2ex(:,jh_tot,:)=0;
C3ez(:,1,:)=-1;
C3ez(:,jh_tot,:)=-1;
C4ez(:,1,:)=0;
C4ez(:,jh_tot,:)=0;
%   z-varying material properties
delbc=upml*delta;
sigmam=-log(rmax)*epsr*epsz*cc*(orderbc+1)/(2*delbc);
sigfactor=sigmam/(delta*(delbc^orderbc)*(orderbc+1));
kmax=1;
kfactor=(kmax-1)/delta/(orderbc+1)/delbc^orderbc;
for k=1:upml
    % Coefficients for field components in the center of the grid cell
    z1=(upml-k+1)*delta;
    z2=(upml-k)*delta;
    sigma=sigfactor*(z1^(orderbc+1)-z2^(orderbc+1));
    ki=1+kfactor*(z1^(orderbc+1)-z2^(orderbc+1));
    facm=(2*epsr*epsz*ki-sigma*dt);
    facp=(2*epsr*epsz*ki+sigma*dt);

    C5ez(:,:,k)=facp;
    C5ez(:,:,ke_tot-k+1)=facp;
    C6ez(:,:,k)=facm;
    C6ez(:,:,ke_tot-k+1)=facm;
    D1hy(:,:,k)=facm/facp;
    D1hy(:,:,ke_tot-k+1)=facm/facp;
    D2hy(:,:,k)=2*epsr*epsz*dt/facp;
    D2hy(:,:,ke_tot-k+1)=2*epsr*epsz*dt/facp;
    D3hx(:,:,k)=facm/facp;
    D3hx(:,:,ke_tot-k+1)=facm/facp;
    D4hx(:,:,k)=1/facp/mur/muz;
    D4hx(:,:,ke_tot-k+1)=1/facp/mur/muz;

    % Coefficients for field components on the grid cell boundary
    z1=(upml-k+1.5)*delta;
    z2=(upml-k+0.5)*delta;
    sigma=sigfactor*(z1^(orderbc+1)-z2^(orderbc+1));
    ki=1+kfactor*(z1^(orderbc+1)-z2^(orderbc+1));
    facm=(2*epsr*epsz*ki-sigma*dt);
    facp=(2*epsr*epsz*ki+sigma*dt);

    C1ey(:,:,k)=facm/facp;
    C1ey(:,:,kh_tot-k+1)=facm/facp;
    C2ey(:,:,k)=2*epsr*epsz*dt/facp;
    C2ey(:,:,kh_tot-k+1)=2*epsr*epsz*dt/facp;
    C3ex(:,:,k)=facm/facp;
    C3ex(:,:,kh_tot-k+1)=facm/facp;
    C4ex(:,:,k)=1/facp/epsr/epsz;
    C4ex(:,:,kh_tot-k+1)=1/facp/epsr/epsz;
    D5hz(:,:,k)=facp;
    D5hz(:,:,kh_tot-k+1)=facp;
    D6hz(:,:,k)=facm;
    D6hz(:,:,kh_tot-k+1)=facm;
end
%   PEC walls
C1ey(:,:,1)=-1;
C1ey(:,:,kh_tot)=-1;
C2ey(:,:,1)=0;
C2ey(:,:,kh_tot)=0;
C3ex(:,:,1)=-1;
C3ex(:,:,kh_tot)=-1;
C4ex(:,:,1)=0;
C4ex(:,:,kh_tot)=0;
%figure
%set(gcf,'DoubleBuffer','on')
%***********************************************************************
%     Begin time stepping loop
%***********************************************************************
for n=1:nmax

    % Update magnetic field
    bstore=bx;
    bx(2:ie_tot,:,:)=D1hx(2:ie_tot,:,:).* bx(2:ie_tot,:,:)-...
                     D2hx(2:ie_tot,:,:).*((ez(2:ie_tot,2:jh_tot,:)-ez(2:ie_tot,1:je_tot,:))-...
                                          (ey(2:ie_tot,:,2:kh_tot)-ey(2:ie_tot,:,1:ke_tot)))./delta;
    hx(2:ie_tot,:,:)= D3hx(2:ie_tot,:,:).*hx(2:ie_tot,:,:)+...
                      D4hx(2:ie_tot,:,:).*(D5hx(2:ie_tot,:,:).*bx(2:ie_tot,:,:)-...
                                           D6hx(2:ie_tot,:,:).*bstore(2:ie_tot,:,:));
    bstore=by;
    by(:,2:je_tot,:)=D1hy(:,2:je_tot,:).* by(:,2:je_tot,:)-...
                     D2hy(:,2:je_tot,:).*((ex(:,2:je_tot,2:kh_tot)-ex(:,2:je_tot,1:ke_tot))-...
                                          (ez(2:ih_tot,2:je_tot,:)-ez(1:ie_tot,2:je_tot,:)))./delta;
    hy(:,2:je_tot,:)= D3hy(:,2:je_tot,:).*hy(:,2:je_tot,:)+...
                      D4hy(:,2:je_tot,:).*(D5hy(:,2:je_tot,:).*by(:,2:je_tot,:)-...
                                           D6hy(:,2:je_tot,:).*bstore(:,2:je_tot,:));
    bstore=bz;
    bz(:,:,2:ke_tot)=D1hz(:,:,2:ke_tot).* bz(:,:,2:ke_tot)-...
                     D2hz(:,:,2:ke_tot).*((ey(2:ih_tot,:,2:ke_tot)-ey(1:ie_tot,:,2:ke_tot))-...
                                          (ex(:,2:jh_tot,2:ke_tot)-ex(:,1:je_tot,2:ke_tot)))./delta;
    hz(:,:,2:ke_tot)= D3hz(:,:,2:ke_tot).*hz(:,:,2:ke_tot)+...
                      D4hz(:,:,2:ke_tot).*(D5hz(:,:,2:ke_tot).*bz(:,:,2:ke_tot)-...
                                           D6hz(:,:,2:ke_tot).*bstore(:,:,2:ke_tot));

    % Update electric field
    dstore=dx;
    dx(:,2:je_tot,2:ke_tot)=C1ex(:,2:je_tot,2:ke_tot).* dx(:,2:je_tot,2:ke_tot)+...
                            C2ex(:,2:je_tot,2:ke_tot).*((hz(:,2:je_tot,2:ke_tot)-hz(:,1:je_tot-1,2:ke_tot))-...
                                                        (hy(:,2:je_tot,2:ke_tot)-hy(:,2:je_tot,1:ke_tot-1)))./delta;
    ex(:,2:je_tot,2:ke_tot)=C3ex(:,2:je_tot,2:ke_tot).*ex(:,2:je_tot,2:ke_tot)+...
                            C4ex(:,2:je_tot,2:ke_tot).*(C5ex(:,2:je_tot,2:ke_tot).*dx(:,2:je_tot,2:ke_tot)-...
                                                        C6ex(:,2:je_tot,2:ke_tot).*dstore(:,2:je_tot,2:ke_tot));
    dstore=dy;
    dy(2:ie_tot,:,2:ke_tot)=C1ey(2:ie_tot,:,2:ke_tot).* dy(2:ie_tot,:,2:ke_tot)+...
                            C2ey(2:ie_tot,:,2:ke_tot).*((hx(2:ie_tot,:,2:ke_tot)-hx(2:ie_tot,:,1:ke_tot-1))-...
                                                        (hz(2:ie_tot,:,2:ke_tot)-hz(1:ie_tot-1,:,2:ke_tot)))./delta;
    ey(2:ie_tot,:,2:ke_tot)=C3ey(2:ie_tot,:,2:ke_tot).*ey(2:ie_tot,:,2:ke_tot)+...
                            C4ey(2:ie_tot,:,2:ke_tot).*(C5ey(2:ie_tot,:,2:ke_tot).*dy(2:ie_tot,:,2:ke_tot)-...
                        &nb ..

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

共2条评分

gwzhao	技术分 +1	有自己的观点	2009-01-02
gwzhao	rf币 +10	感谢您的资料	2009-01-02

离线wq_463

UID ：20925

注册：2008-11-06
登录：2021-04-22
发帖：227
等级：仿真三级

1楼发表于: 2009-01-02 16:30:05

% Coefficients for field components in the center of the grid cell
    x1=(upml-i+1)*delta;
    x2=(upml-i)*delta;
    sigma=sigfactor*(x1^(orderbc+1)-x2^(orderbc+1));
    ki=1+kfactor*(x1^(orderbc+1)-x2^(orderbc+1));
    facm=(2*epsr*epsz*ki-sigma*dt);
    facp=(2*epsr*epsz*ki+sigma*dt);
    C5ex(i,:,:)=facp;
    C5ex(ie_tot-i+1,:,:)=facp;
    C6ex(i,:,:)=facm;
    C6ex(ie_tot-i+1,:,:)=facm;
    D1hz(i,:,:)=facm/facp;
    D1hz(ie_tot-i+1,:,:)=facm/facp;
    D2hz(i,:,:)=2.0*epsr*epsz*dt/facp;
    D2hz(ie_tot-i+1,:,:)=2.0*epsr*epsz*dt/facp;
    D3hy(i,:,:)=facm/facp;
    D3hy(ie_tot-i+1,:,:)=facm/facp;
    D4hy(i,:,:)=1.0/facp/mur/muz;
    D4hy(ie_tot-i+1,:,:)=1.0/facp/mur/muz;
    % Coefficients for field components on the grid cell boundary
    x1=(upml-i+1.5)*delta;
    x2=(upml-i+0.5)*delta;
    sigma=sigfactor*(x1^(orderbc+1)-x2^(orderbc+1));
    ki=1.0+kfactor*(x1^(orderbc+1)-x2^(orderbc+1));
    facm=(2.0*epsr*epsz*ki-sigma*dt);
    facp=(2.0*epsr*epsz*ki+sigma*dt);
    C1ez(i,:,:)=facm/facp;
    C1ez(ih_tot-i+1,:,:)=facm/facp;
    C2ez(i,:,:)=2.0*epsr*epsz*dt/facp;
    C2ez(ih_tot-i+1,:,:)=2.0*epsr*epsz*dt/facp;
    C3ey(i,:,:)=facm/facp;
    C3ey(ih_tot-i+1,:,:)=facm/facp;
    C4ey(i,:,:)=1.0/facp/epsr/epsz;
    C4ey(ih_tot-i+1,:,:)=1.0/facp/epsr/epsz;
    D5hx(i,:,:)=facp;
    D5hx(ih_tot-i+1,:,:)=facp;
    D6hx(i,:,:)=facm;
    D6hx(ih_tot-i+1,:,:)=facm;
1.C5,C6为什么要跟C1,C2,C3,C4分开定义呢,同样的还有D的参数是不是C5ex,C6ex,D1hz,D2hz,D3hy,D4hy含有X方向的参数呢?
2.还有为什么X1,X2两次定义的不一样呢
Coefficients for field components in the center of the grid cell
Coefficients for field components on the grid cell boundary
怎么理解这2句话呢
3.C1ez(i,:,:)=facm/facp;
    C1ez(ih_tot-i+1,:,:)=facm/facp;
ih_tot-i+1这个怎么理解呢?

离线hughsq

UID ：2114

注册：2007-04-24
登录：2009-07-25
发帖：70
等级：仿真二级

2楼发表于: 2009-01-03 15:49:36

这段程序我研究了不少时间，你可对照Taflove的2000年版的书看看，程序基本和书上一致
这是我研究时的一段注释，希望对你有帮助。
弄懂了这段程序，就明白了UPML方法的好处就是：问题域和UPML区同时迭代处理，简洁高效。

R_0=exp(-16);   %设定的反射误差, designated as R(0) in Equation 7.62
m_bc=4;         %order of the polynomial grading, designated as m in Equation 7.60a,b

%%   X方向,即前后UPML层参数填充
d_bc=UPML*delta;%UPML层的物理厚度
sigmam=-log(R_0)*(m_bc+1.0)/(2.0*Z0r*d_bc); %第二版的式7.57
sigfactor=sigmam/(delta*(d_bc^m_bc)*(m_bc+1.0));
kmax=1;%???????????
kfactor=(kmax-1.0)/delta/(m_bc+1.0)/d_bc^m_bc;

for I=1:UPML%前后两面UPML层

    % Coefficients for field components in the center of the grid cell
    x1=(UPML-I+1)*delta;%x1=x2+delta
    x2=(UPML-I)*delta;
    sigma=sigfactor*(x1^(m_bc+1)-x2^(m_bc+1));%第二版的式7.55a,b
    ki=1+kfactor*(x1^(m_bc+1)-x2^(m_bc+1));
    facm=(2*Epsr*Eps0*ki-sigma*dt);%factor minus,分子
    facp=(2*Epsr*Eps0*ki+sigma*dt);%factor plus ,分母
%UPML区系数赋值
    C5Ex(I,:,:)=facp;
    C5Ex(IE_all-I+1,:,:)=facp;
    C6Ex(I,:,:)=facm;
    C6Ex(IE_all-I+1,:,:)=facm;
    D1Hz(I,:,:)=facm/facp;
    D1Hz(IE_all-I+1,:,:)=facm/facp;
    D2Hz(I,:,:)=2.0*Epsr*Eps0*dt/facp;
    D2Hz(IE_all-I+1,:,:)=2.0*Epsr*Eps0*dt/facp;
    D3Hy(I,:,:)=facm/facp;
    D3Hy(IE_all-I+1,:,:)=facm/facp;
    D4Hy(I,:,:)=1.0/facp/Mur/Mu0;
    D4Hy(IE_all-I+1,:,:)=1.0/facp/Mur/Mu0;

    % Coefficients for field components on the grid cell boundary
    x1=(UPML-I+1.5)*delta;
    x2=(UPML-I+0.5)*delta;
    sigma=sigfactor*(x1^(m_bc+1)-x2^(m_bc+1));
    ki=1.0+kfactor*(x1^(m_bc+1)-x2^(m_bc+1));
    facm=(2.0*Epsr*Eps0*ki-sigma*dt);
    facp=(2.0*Epsr*Eps0*ki+sigma*dt);

    C1Ez(I,:,:)=facm/facp;
    C1Ez(IH_all-I+1,:,:)=facm/facp;
    C2Ez(I,:,:)=2.0*Epsr*Eps0*dt/facp;
    C2Ez(IH_all-I+1,:,:)=2.0*Epsr*Eps0*dt/facp;
    C3Ey(I,:,:)=facm/facp;
    C3Ey(IH_all-I+1,:,:)=facm/facp;
    C4Ey(I,:,:)=1.0/facp/Epsr/Eps0;
    C4Ey(IH_all-I+1,:,:)=1.0/facp/Epsr/Eps0;
    D5Hx(I,:,:)=facp;
    D5Hx(IH_all-I+1,:,:)=facp;
    D6Hx(I,:,:)=facm;
    D6Hx(IH_all-I+1,:,:)=facm;

end

%   PEC walls
C1Ez(1,:,:)=-1.0;
C1Ez(IH_all,:,:)=-1.0;
C2Ez(1,:,:)=0.0;
C2Ez(IH_all,:,:)=0.0;
C3Ey(1,:,:)=-1.0;
C3Ey(IH_all,:,:)=-1.0;
C4Ey(1,:,:)=0.0;
C4Ey(IH_all,:,:)=0.0;

学有用的知识，热爱从事的事业，做有用的创新。

离线wq_463

UID ：20925

注册：2008-11-06
登录：2021-04-22
发帖：227
等级：仿真三级

3楼发表于: 2009-01-03 19:10:33

引用第2楼hughsq于2009-01-03 15:49发表的 : ^w/_hY!4/
这段程序我研究了不少时间，你可对照Taflove的2000年版的书看看，程序基本和书上一致 JBo/<W#|
这是我研究时的一段注释，希望对你有帮助。 SxdH%agM
弄懂了这段程序，就明白了UPML方法的好处就是：问题域和UPML区同时迭代处理，简洁高效。
R_0=exp(-16); %设定的反射误差, designated as R(0) in Equation 7.62 Po#;SG#Ee
.......

感谢楼上的解答问题是这段程序中的系数是为什么这么定义,而且是分开定义呢?
截断边界为什么要这么设置呢?
顺便问一下楼上的qq多少,能否告知,我想跟你交流一下
还有我手头这本书版本不对,不知你能把第二版给我传一下

离线wq_463

UID ：20925

注册：2008-11-06
登录：2021-04-22
发帖：227
等级：仿真三级

4楼发表于: 2009-01-04 20:34:30

to gwzhao师兄或者哪位高人给我解释一下
C5ez(:,:,k)=facp;
    C5ez(:,:,ke_tot-k+1)=facp;
    C6ez(:,:,k)=facm;
    C6ez(:,:,ke_tot-k+1)=facm;
    D1hy(:,:,k)=facm/facp;
    D1hy(:,:,ke_tot-k+1)=facm/facp;
    D2hy(:,:,k)=2*epsr*epsz*dt/facp;
    D2hy(:,:,ke_tot-k+1)=2*epsr*epsz*dt/facp;
    D3hx(:,:,k)=facm/facp;
    D3hx(:,:,ke_tot-k+1)=facm/facp;
    D4hx(:,:,k)=1/facp/mur/muz;
    D4hx(:,:,ke_tot-k+1)=1/facp/mur/muz;

    % Coefficients for field components on the grid cell boundary
    z1=(upml-k+1.5)*delta;
    z2=(upml-k+0.5)*delta;
    sigma=sigfactor*(z1^(orderbc+1)-z2^(orderbc+1));
    ki=1+kfactor*(z1^(orderbc+1)-z2^(orderbc+1));
    facm=(2*epsr*epsz*ki-sigma*dt);
    facp=(2*epsr*epsz*ki+sigma*dt);

    C1ey(:,:,k)=facm/facp;
    C1ey(:,:,kh_tot-k+1)=facm/facp;
    C2ey(:,:,k)=2*epsr*epsz*dt/facp;
    C2ey(:,:,kh_tot-k+1)=2*epsr*epsz*dt/facp;
    C3ex(:,:,k)=facm/facp;
    C3ex(:,:,kh_tot-k+1)=facm/facp;
    C4ex(:,:,k)=1/facp/epsr/epsz;
    C4ex(:,:,kh_tot-k+1)=1/facp/epsr/epsz;
    D5hz(:,:,k)=facp;
    D5hz(:,:,kh_tot-k+1)=facp;
    D6hz(:,:,k)=facm;
    D6hz(:,:,kh_tot-k+1)=facm;
为什么要分开定义参数,按这个字面意思就是位于网格中心的场量和位于网格边界的场量分别定义的参数,可是我不是很理解啊

离线herosword

UID ：3621

注册：2007-07-01
登录：2019-11-03
发帖：30
等级：仿真二级

5楼发表于: 2009-01-05 15:55:48

你看程序需要看他程序所涉及的文章,因为那里才有相关的离散迭代公式,不然看程序会是一头雾水
个人认为

离线herosword

UID ：3621

注册：2007-07-01
登录：2019-11-03
发帖：30
等级：仿真二级

6楼发表于: 2009-01-05 16:00:25

我有第二版的书,你加我qq,我传给你,22348572

离线wq_463

UID ：20925

注册：2008-11-06
登录：2021-04-22
发帖：227
等级：仿真三级

7楼发表于: 2009-01-05 19:31:38

引用第5楼herosword于2009-01-05 15:55发表的 : <j"}EEb^
你看程序需要看他程序所涉及的文章,因为那里才有相关的离散迭代公式,不然看程序会是一头雾水 *c'nPa$+|S
个人认为

首先谢谢你的解答,我不是这个迭代公式不会推或者看不明白,这几个参数的公式表达是没有问题我推出来也是这个结果,我是说他为什么要把系数分开定义,而且C1EX,C2EX,C3EX,C4EX,C5EX,C6EX这几个参数他定义的也不全啊只有C5EX,C6EX.
还有D的系数,

离线cem-uestc

UID ：9061

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

8楼发表于: 2009-01-05 20:09:40

这应该是编程问题，不同人的编程思路不一样。
分开写的目的是方便计算各场分量在不同的PML区域场量计算，各场量系数在不同PML区域（各区域的PML材料参数不同）造成的。
只要把各PML区域材料参数搞清楚，对应哪一个场量分析清楚，我想PML的思路就清楚了。

共1条评分

gwzhao

rf币 +5

积极参与

2009-01-05

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

离线wq_463

UID ：20925

注册：2008-11-06
登录：2021-04-22
发帖：227
等级：仿真三级

9楼发表于: 2009-01-06 16:13:41

引用第8楼cem-uestc于2009-01-05 20:09发表的 : %(4G[R[
这应该是编程问题，不同人的编程思路不一样。 tPO\e]
分开写的目的是方便计算各场分量在不同的PML区域场量计算，各场量系数在不同PML区域（各区域的PML材料参数不同）造成的。 ibDMhW$n
只要把各PML区域材料参数搞清楚，对应哪一个场量分析清楚，我想PML的思路就清楚了。

Coefficients for field components in the center of the grid cell是说位于网格中央的场量的系数是否可以理解成磁场量HX,HY,HZ是位于网格中央的所以这么定义?
Coefficients for field components on the grid cell boundary是说位于网格边界上的场量的系数是否可以理解成为电场量EX,EY,EZ是位于网格愣边的所以这么定义?

发帖回复