FDTD参数选择估计程序
针对二阶精度的时域有限差分程序.
现可直接调用的源信号是:一个周期的正弦信号,高期脉冲,ricker子波.
其它信号可手动修改源信号接口,或源生成函数.
---------------
请函数.
%************************************************************
% 1. determine maximum possible spatial field discretization.
% (in order to avoid numerical dispersion).(5 grid points per
% minimum wavelength are needed to avoid dispersion).
% 2. find the maximum possible time step using this dx and dz.
% (in order to avoid numerical instability).
% Coded by yiling. Email: yiling@email.jlu.edu.cn
% Date: 2008
%*************************************************************************+
clear;
clc;
%--------------------------------------------------------------------------
dx=0.02; % (m)
dy=0.02; % (m)
epsilonmax=25; % Epsion. maximum relative dielectric permittivity.
mumax=1; % Mu. maximum relative magnetic permeability.
sourcetype='ricker'; % can be 'cont_sine', 'gaussian', 'ricker'.
freq=100e6; % (Hz)
amp=1; % amplitude.
thres=0.02; % threshold to determine maximum frequency in source pulse.(proposed = 0.02).
%--------------------------------------------------------------------------
Timewindows=528; % (ns)
%--------------------------------------------------------------------------
%*************************************************************************+
%--------------------------------------------------------------------------
vlight=0.3;
epsilonmin=1; % Epsion. minimum relative dielectric permittivity.
mumin=1; % Mu. minimum relative magnetic permeability.
%--------------------------------------------------------------------------
dt=1/(vlight*sqrt(1/dx^2+1/dy^2));
% minwavelength=vlight/sqrt(epsilinmax);
%--------------------------------------------------------------------------
t=0:dt:Timewindows;
dt=dt*1e-9;
t=t*1e-9;
Timewindows=Timewindows*1e-9;
source=gprmaxso(sourcetype,amp,freq,dt,Timewindows);
[dxmax,wlmin,fmax] = finddx(epsilonmax,mumax,source,t,thres);
%--------------------------------------------------------------------------
disp('----------------------------------------------------------------- ');
disp(['Maximum frequency contained in source pulse = ',num2str(fmax/1e6),' MHz']);
disp(['Minimum wavelength in simulation grid = ',num2str(wlmin),' m']);
disp(['Maximum possible electric/magnetic field discretization (dx,dy) = ',num2str(dxmax),' m']);
disp(' ');
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
dtmax = finddt(epsilonmin,mumin,dxmax,dxmax);
disp(['Maximum possible time step with this discretization = ',num2str(dtmax/1e-9),' ns']);
disp('----------------------------------------------------------------- ');
%**************************************************
子函数1
function dtmax = finddt(epmin,mumin,dx,dz);
% finddt.m
%
% This function finds the maximum time step that can be used in the 2-D
% FDTD modeling codes TM_model2d.m and TE_model2d.m, such that they remain
% numerically stable. Second-order-accurate time and fourth-order-accurate
% spatial derivatives are assumed (i.e., O(2,4)).
%
% Syntax: dtmax = finddt(epmin,mumin,dx,dz)
%
% where dtmax = maximum time step for FDTD to be stable
% epmin = minimum relative dielectric permittivity in grid
% mumin = minimum relative magnetic permeability in grid
% dx = spatial discretization in x-direction (m)
% dz = spatial discretization in z-direction (m)
%
% by James Irving
% July 2005
% convert relative permittivity and permeability to true values
mu0 = 1.2566370614e-6;
ep0 = 8.8541878176e-12;
epmin = epmin*ep0;
mumin = mumin*mu0;
% determine maximum allowable time step for numerical stability
dtmax = 6/7*sqrt(epmin*mumin/(1/dx^2 + 1/dz^2));
子函数2
function [dxmax,wlmin,fmax] = finddx(epmax,mumax,srcpulse,t,thres);
% finddx.m
%
% This function finds the maximum spatial discretization that can be used in the
% 2-D FDTD modeling codes TM_model2d.m and TE_model2d.m, such that numerical
% dispersion is avoided. Second-order accurate time and fourth-order-accurate
% spatial derivatives are assumed (i.e., O(2,4)). Consequently, 5 field points
% per minimum wavelength are required.
%
% Note: The dx value obtained with this program is needed to compute the maximum
% time step (dt) that can be used to avoid numerical instability. However, the
% time vector and source pulse are required in this code to determine the highest
% frequency component in the source pulse. For this program, make sure to use a fine
加入微信
获取电子行业最新资讯
搜索微信公众号:EEPW
或用微信扫描左侧二维码