【数据分析】基于有限差分时域（FDTD）方法实现微带结构的全波分析附Matlab代码

matlab科研助手

于 2022-09-27 17:03:17 发布

阅读量1k

点赞数 1

CC 4.0 BY-SA版权

分类专栏：图像处理文章标签：数据分析 matlab 数据挖掘

本文链接：https://blog.youkuaiyun.com/matlab_dingdang/article/details/127074502

图像处理专栏收录该内容

872 篇文章

订阅专栏

✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，matlab项目合作可私信。

🍎个人主页：Matlab科研工作室

🍊个人信条：格物致知。

更多Matlab仿真内容点击👇

智能优化算法神经网络预测雷达通信无线传感器

信号处理图像处理路径规划元胞自动机无人机

⛄ 内容介绍

将直接三维有限差分时域（FDTD）方法应用于各种微带结构的全波分析。该方法被证明是对复杂的微带电路元件和微带天线进行建模的有效工具。从时域结果计算出线馈矩形贴片天线的输入阻抗以及低通滤波器和支线耦合器的频率相关散射参数。制作了这些电路，并将对它们进行的测量与 FDTD 结果进行了比较，结果表明它们非常吻合。

⛄ 部分代码

%% Microstrip low-pass filter analysis using 3D FDTD code with UPML

%% absorbing borders (ABC)

% Here we use FDTD 3D with UPML to calculate scattering coefficients S_{11}

% and S_{21} for planar microstrip low-pass filter following by the original

% paper by D. Sheen, S. Ali, M. Abouzahra, J. Kong "Application of the

% Three-Dimensional Finite-Difference Time-Domain Method to the Analysis of

% planar Microstrip Circuits", IEEE Trans. on Microwave Theory and Techniques

% (http://dx.doi.org/10.1109/22.55775).

% Also |S_{21}| dependence can be taken from the paper "Computational

% electromagnetic method for interconnects and small structures" by C.

% Balanis, A. Policarpou and S. Georgakopoulos

% (http://dx.doi.org/10.1006/spmi.2000.0865)

function FDTD_3D_Lowpass

close all; clear; clc;

%% Physical constants

epsilon0 = 8.85418782e-12; mu0 = 1.25663706e-6;

c = 1.0/sqrt(mu0*epsilon0);

%% Gaussian half-width

t_half = 15.0e-12;

%% Microstrip transmission lines parameters

lineW = 2.413e-3;

lineH = 1.0e-3;

lineEr = 2.2;

Z0 = 49.2526;

%% End time

t_end = 1.5e-9;

%% Total mesh dimensions and grid cells sizes (without PML)

nx = 80; ny = 100; nz = 16;

dx = 0.4064e-3; dy = 0.4233e-3; dz = 0.2650e-3;

%% Number of PML layers

PML = 5;

%% Matrix of material's constants

number_of_materials = 4;

% For material of number x = 1,2,3... :

% Material(x,1) - relative permittivity, Material(x,2) - relative permeability,

% Material(x,3) - specific conductivity

% Vacuum

Material(1,1) = 1.0; Material(1,2) = 1.0; Material(1,3) = 0.0;

% Metal (Copper)

Material(2,1) = 1.0; Material(2,2) = 1.0; Material(2,3) = 5.88e+7;

% Substrate material (RT/Duroid 5880)

Material(3,1) = lineEr; Material(3,2) = 1.0; Material(3,3) = 0.0;

% Matched load material is calculated from transmission line parameters

Material(4,1) = 1.0; Material(4,2) = 1.0; Material(4,3) = lineH/(Z0*lineW*dy);

% Add PML layers

nx = nx + 2*PML; ny = ny + 2*PML; nz = nz + 2*PML;

% Calculate dt

dt = (1.0/c/sqrt( 1.0/(dx^2) + 1.0/(dy^2) + 1.0/(dz^2)))*0.9999;

number_of_iterations = ceil(t_end/dt);

%% 3D array for geometry

Index = ones(nx, ny, nz);

%% Define of low-pass filter geometry

% Ground plane

Index((1+PML):(nx-PML), (1+PML):(ny-PML), PML+1) = 2;

% Rectangular patch (one cell thickness)

Index((nx/2-25):(nx/2+25), (ny/2-3):(ny/2+3), PML+5) = 2;

% Transmission line from port 1 to rectangular patch

Index((nx/2-10):(nx/2-5), (PML+1):ny/2, PML+5) = 2;

% Transmission line from rectangular patch to port 2

Index((nx/2+5):(nx/2+10), ny/2:(ny-PML), PML+5) = 2;

% Dielectric substrate between ground plane and filter patch

Index((1+PML):(nx-PML), (1+PML):(ny-PML), (PML+2):(PML+4)) = 3;

% Matched load before port 1

Index((nx/2-10):(nx/2-5), PML+1, (PML+2):(PML+4)) = 4;

% Matched load after port 2

Index((nx/2+5):(nx/2+10), ny-PML, (PML+2):(PML+4)) = 4;

%% 3D FDTD physical (fields) and additional arrays are defined as 'single'

%% to increase performance

Ex = zeros(nx, ny+1, nz+1, 'single');

Gx = zeros(nx, ny+1, nz+1, 'single');

Fx = zeros(nx, ny+1, nz+1, 'single');

Ey = zeros(nx+1, ny, nz+1, 'single');

Gy = zeros(nx+1, ny, nz+1, 'single');

Fy = zeros(nx+1, ny, nz+1, 'single');

Ez = zeros(nx+1, ny+1, nz, 'single');

Gz = zeros(nx+1, ny+1, nz, 'single');

Fz = zeros(nx+1, ny+1, nz, 'single');

Hx = zeros(nx+1, ny, nz, 'single');

Bx = zeros(nx+1, ny, nz, 'single');

Hy = zeros(nx, ny+1, nz, 'single');

By = zeros(nx, ny+1, nz, 'single');

Hz = zeros(nx, ny, nz+1, 'single');

Bz = zeros(nx, ny, nz+1, 'single');

%% FDTD PML coefficients arrays. Here they are already filled with values

%% corresponding to free space

m = 4; ka_max = 1.0; R_err = 1.0e-16;

eta = sqrt(mu0/epsilon0*Material(1,1)/Material(1,2));

k_Ex_c = ones(nx, ny, nz, 'single')*2.0*epsilon0;

k_Ex_d = ones(nx, ny, nz, 'single')*(-2.0*epsilon0);

k_Ey_a = ones(nx+1, ny, nz, 'single');

k_Ey_b = ones(nx+1, ny, nz, 'single')/(2.0*epsilon0);

k_Gz_a = ones(nx+1, ny, nz, 'single');

k_Gz_b = ones(nx+1, ny, nz, 'single');

k_Hy_a = ones(nx, ny, nz, 'single');

k_Hy_b = ones(nx, ny, nz, 'single')/(2.0*epsilon0);

k_Hx_c = ones(nx+1, ny, nz, 'single')*2.0*epsilon0/mu0;

k_Hx_d = ones(nx+1, ny, nz, 'single')*(-2.0*epsilon0/mu0);

k_Bz_a = ones(nx, ny, nz, 'single');

k_Bz_b = ones(nx, ny, nz, 'single')*dt;

k_Gx_a = ones(nx, ny+1, nz, 'single');

k_Gx_b = ones(nx, ny+1, nz, 'single');

k_Ey_c = ones(nx, ny, nz, 'single')*2.0*epsilon0;

k_Ey_d = ones(nx, ny, nz, 'single')*(-2.0*epsilon0);

k_Ez_a = ones(nx, ny+1, nz, 'single');

k_Ez_b = ones(nx, ny+1, nz, 'single')/(2.0*epsilon0);

k_Bx_a = ones(nx, ny, nz, 'single');

k_Bx_b = ones(nx, ny, nz, 'single')*dt;

k_Hy_c = ones(nx, ny+1, nz, 'single')*2.0*epsilon0/mu0;

k_Hy_d = ones(nx, ny+1, nz, 'single')*(-2.0*epsilon0/mu0);

k_Hz_a = ones(nx, ny, nz, 'single');

k_Hz_b = ones(nx, ny, nz, 'single')/(2.0*epsilon0);

k_Ex_a = ones(nx, ny, nz+1, 'single');

k_Ex_b = ones(nx, ny, nz+1, 'single')/(2.0*epsilon0);

k_Gy_a = ones(nx, ny, nz+1, 'single');

k_Gy_b = ones(nx, ny, nz+1, 'single');

k_Ez_c = ones(nx, ny, nz, 'single')*2.0*epsilon0;

k_Ez_d = ones(nx, ny, nz, 'single')*(-2.0*epsilon0);

k_Hx_a = ones(nx, ny, nz, 'single');

k_Hx_b = ones(nx, ny, nz, 'single')/(2.0*epsilon0);

k_By_a = ones(nx, ny, nz, 'single');

k_By_b = ones(nx, ny, nz, 'single')*dt;

k_Hz_c = ones(nx, ny, nz+1, 'single')*2.0*epsilon0/mu0;

k_Hz_d = ones(nx, ny, nz+1, 'single')*(-2.0*epsilon0/mu0);

%% General FDTD coefficients

I = 1:number_of_materials;

K_a(I) = (2.0*epsilon0*Material(I,1) - Material(I,3)*dt)./...

(2.0*epsilon0*Material(I,1) + Material(I,3)*dt);

K_b(I) = 2.0*dt./(2.0*epsilon0*Material(I,1) + Material(I,3)*dt);

K_c(I) = Material(I,2);

Ka = single(K_a(Index)); Kb = single(K_b(Index)); Kc = single(K_c(Index));

%% PML coefficients along x-axis

sigma_max = -(m + 1.0)*log(R_err)/(2.0*eta*PML*dx);

for I=0:(PML-1)

sigma_x = sigma_max*((PML - I)/PML)^m;

ka_x = 1.0 + (ka_max - 1.0)*((PML - I)/PML)^m;

k_Ey_a(I+1,:,:) = (2.0*epsilon0*ka_x - sigma_x*dt)/...

(2.0*epsilon0*ka_x + sigma_x*dt);

k_Ey_a(nx-I,:,:) = k_Ey_a(I+1,:,:);

k_Ey_b(I+1,:,:) = 1.0/(2.0*epsilon0*ka_x + sigma_x*dt);

k_Ey_b(nx-I,:,:) = k_Ey_b(I+1,:,:);

k_Gz_a(I+1,:,:) = (2.0*epsilon0*ka_x - sigma_x*dt)/...

(2.0*epsilon0*ka_x + sigma_x*dt);

k_Gz_a(nx-I,:,:) = k_Gz_a(I+1,:,:);

k_Gz_b(I+1,:,:) = 2.0*epsilon0/(2.0*epsilon0*ka_x + sigma_x*dt);

k_Gz_b(nx-I,:,:) = k_Gz_b(I+1,:,:);

k_Hx_c(I+1,:,:) = (2.0*epsilon0*ka_x + sigma_x*dt)/mu0;

k_Hx_c(nx-I,:,:) = k_Hx_c(I+1,:,:);

k_Hx_d(I+1,:,:) = -(2.0*epsilon0*ka_x - sigma_x*dt)/mu0;

k_Hx_d(nx-I,:,:) = k_Hx_d(I+1,:,:);

sigma_x = sigma_max*((PML - I - 0.5)/PML)^m;

ka_x = 1.0 + (ka_max - 1.0)*((PML - I - 0.5)/PML)^m;

k_Ex_c(I+1,:,:) = 2.0*epsilon0*ka_x + sigma_x*dt;

k_Ex_c(nx-I-1,:,:) = k_Ex_c(I+1,:,:);

k_Ex_d(I+1,:,:) = -(2.0*epsilon0*ka_x - sigma_x*dt);

k_Ex_d(nx-I-1,:,:) = k_Ex_d(I+1,:,:);

k_Hy_a(I+1,:,:) = (2.0*epsilon0*ka_x - sigma_x*dt)/...

(2.0*epsilon0*ka_x + sigma_x*dt);

k_Hy_a(nx-I-1,:,:) = k_Hy_a(I+1,:,:);

k_Hy_b(I+1,:,:) = 1.0/(2.0*epsilon0*ka_x + sigma_x*dt);

k_Hy_b(nx-I-1,:,:) = k_Hy_b(I+1,:,:);

k_Bz_a(I+1,:,:) = (2.0*epsilon0*ka_x - sigma_x*dt)/...

(2.0*epsilon0*ka_x + sigma_x*dt);

k_Bz_a(nx-I-1,:,:) = k_Bz_a(I+1,:,:);

k_Bz_b(I+1,:,:) = 2.0*epsilon0*dt/(2.0*epsilon0*ka_x + sigma_x*dt);

k_Bz_b(nx-I-1,:,:) = k_Bz_b(I+1,:,:);

end

%% PML coefficients along y-axis

sigma_max = -(m + 1.0)*log(R_err)/(2.0*eta*PML*dy);

for J=0:(PML-1)

sigma_y = sigma_max*((PML - J)/PML)^m;

ka_y = 1.0 + (ka_max - 1.0)*((PML - J)/PML)^m;

k_Gx_a(:,J+1,:) = (2.0*epsilon0*ka_y - sigma_y*dt)/...

(2.0*epsilon0*ka_y + sigma_y*dt);

k_Gx_a(:,ny-J,:) = k_Gx_a(:,J+1,:);

k_Gx_b(:,J+1,:) = 2.0*epsilon0/(2.0*epsilon0*ka_y + sigma_y*dt);

k_Gx_b(:,ny-J,:) = k_Gx_b(:,J+1,:);

k_Ez_a(:,J+1,:) = (2.0*epsilon0*ka_y - sigma_y*dt)/...

(2.0*epsilon0*ka_y + sigma_y*dt);

k_Ez_a(:,ny-J,:) = k_Ez_a(:,J+1,:);

k_Ez_b(:,J+1,:) = 1.0/(2.0*epsilon0*ka_y + sigma_y*dt);

k_Ez_b(:,ny-J,:) = k_Ez_b(:,J+1,:);

k_Hy_c(:,J+1,:) = (2.0*epsilon0*ka_y + sigma_y*dt)/mu0;

k_Hy_c(:,ny-J,:) = k_Hy_c(:,J+1,:);

k_Hy_d(:,J+1,:) = -(2.0*epsilon0*ka_y - sigma_y*dt)/mu0;

k_Hy_d(:,ny-J,:) = k_Hy_d(:,J+1,:);

sigma_y = sigma_max*((PML - J - 0.5)/PML)^m;

ka_y = 1.0 + (ka_max - 1.0)*((PML - J - 0.5)/PML)^m;

k_Ey_c(:,J+1,:) = 2.0*epsilon0*ka_y+sigma_y*dt;

k_Ey_c(:,ny-J-1,:) = k_Ey_c(:,J+1,:);

k_Ey_d(:,J+1,:) = -(2.0*epsilon0*ka_y-sigma_y*dt);

k_Ey_d(:,ny-J-1,:) = k_Ey_d(:,J+1,:);

k_Bx_a(:,J+1,:) = (2.0*epsilon0*ka_y-sigma_y*dt)/...

(2.0*epsilon0*ka_y+sigma_y*dt);

k_Bx_a(:,ny-J-1,:) = k_Bx_a(:,J+1,:);

k_Bx_b(:,J+1,:) = 2.0*epsilon0*dt/(2.0*epsilon0*ka_y+sigma_y*dt);

k_Bx_b(:,ny-J-1,:) = k_Bx_b(:,J+1,:);

k_Hz_a(:,J+1,:) = (2.0*epsilon0*ka_y-sigma_y*dt)/...

(2.0*epsilon0*ka_y+sigma_y*dt);

k_Hz_a(:,ny-J-1,:) = k_Hz_a(:,J+1,:);

k_Hz_b(:,J+1,:) = 1.0/(2.0*epsilon0*ka_y+sigma_y*dt);

k_Hz_b(:,ny-J-1,:) = k_Hz_b(:,J+1,:);

end

%% PML coefficients along z-axis

sigma_max = -(m + 1.0)*log(R_err)/(2.0*eta*PML*dz);

for K=0:(PML-1)

sigma_z = sigma_max*((PML - K)/PML)^m;

ka_z = 1.0 + (ka_max - 1.0)*((PML-K)/PML)^m;

k_Ex_a(:,:,K+1) = (2.0*epsilon0*ka_z - sigma_z*dt)/...

(2.0*epsilon0*ka_z + sigma_z*dt);

k_Ex_a(:,:,nz-K) = k_Ex_a(:,:,K+1);

k_Ex_b(:,:,K+1) = 1.0/(2.0*epsilon0*ka_z + sigma_z*dt);

k_Ex_b(:,:,nz-K) = k_Ex_b(:,:,K+1);

k_Gy_a(:,:,K+1) = (2.0*epsilon0*ka_z - sigma_z*dt)/...

(2.0*epsilon0*ka_z + sigma_z*dt);

k_Gy_a(:,:,nz-K) = k_Gy_a(:,:,K+1);

k_Gy_b(:,:,K+1) = 2.0*epsilon0/(2.0*epsilon0*ka_z + sigma_z*dt);

k_Gy_b(:,:,nz-K) = k_Gy_b(:,:,K+1);

k_Hz_c(:,:,K+1) = (2.0*epsilon0*ka_z + sigma_z*dt)/mu0;

k_Hz_c(:,:,nz-K) = k_Hz_c(:,:,K+1);

k_Hz_d(:,:,K+1) = -(2.0*epsilon0*ka_z - sigma_z*dt)/mu0;

k_Hz_d(:,:,nz-K) = k_Hz_d(:,:,K+1);

sigma_z = sigma_max*((PML - K - 0.5)/PML)^m;

ka_z = 1.0 + (ka_max - 1.0)*((PML - K - 0.5)/PML)^m;

k_Ez_c(:,:,K+1) = 2.0*epsilon0*ka_z + sigma_z*dt;

k_Ez_c(:,:,nz-K-1) = k_Ez_c(:,:,K+1);

k_Ez_d(:,:,K+1) = -(2.0*epsilon0*ka_z - sigma_z*dt);

k_Ez_d(:,:,nz-K-1) = k_Ez_d(:,:,K+1);

k_Hx_a(:,:,K+1) = (2.0*epsilon0*ka_z - sigma_z*dt)/...

(2.0*epsilon0*ka_z + sigma_z*dt);

k_Hx_a(:,:,nz-K-1) = k_Hx_a(:,:,K+1);

k_Hx_b(:,:,K+1) = 1.0/(2.0*epsilon0*ka_z + sigma_z*dt);

k_Hx_b(:,:,nz-K-1) = k_Hx_b(:,:,K+1);

k_By_a(:,:,K+1) = (2.0*epsilon0*ka_z - sigma_z*dt)/...

(2.0*epsilon0*ka_z + sigma_z*dt);

k_By_a(:,:,nz-K-1) = k_By_a(:,:,K+1);

k_By_b(:,:,K+1) = 2.0*epsilon0*dt/(2.0*epsilon0*ka_z + sigma_z*dt);

k_By_b(:,:,nz-K-1) = k_By_b(:,:,K+1);

end

%% Main 3D FDTD+UPML routine (operates with 'singles' to increase speed)

hhh = waitbar(0, 'Calculations in progress...');

tic;

for T=0:(number_of_iterations-1)

%% Calculate Fx -> Gx -> Ex

I = 1:nx; J = 2:ny; K = 2:nz;

Fx_r = Fx(I,J,K);

Fx(I,J,K) = Ka(I,J,K).*Fx(I,J,K) + Kb(I,J,K).*...

((Hz(I,J,K) - Hz(I,J-1,K))/dy - (Hy(I,J,K) - Hy(I,J,K-1))/dz);

Gx_r = Gx(I,J,K);

Gx(I,J,K) = k_Gx_a(I,J,K).*Gx(I,J,K) + k_Gx_b(I,J,K).*(Fx(I,J,K) - Fx_r);

Ex(I,J,K) = k_Ex_a(I,J,K).*Ex(I,J,K) + k_Ex_b(I,J,K).*...

(k_Ex_c(I,J,K).*Gx(I,J,K) + k_Ex_d(I,J,K).*Gx_r);

%% Calculate Fy -> Gy -> Ey

I = 2:nx; J = 1:ny; K = 2:nz;

Fy_r = Fy(I,J,K);

⛄ 运行结果

⛄ 参考文献

[1] Sheen D M , Ali S M , MD Abouzahra, et al. Application of the three-dimensional finite-difference time-domain method to the analysis of planar microstrip circuits[J]. IEEE Transactions on Microwave Theory & Techniques, 1990, MTT-38(7):849-857.

❤️ 关注我领取海量matlab电子书和数学建模资料

❤️部分理论引用网络文献，若有侵权联系博主删除