【通信】基于matlab FDTD法研究移动通信终端电磁辐射对人体的影响【含Matlab源码 761期】

221 阅读7分钟

一、简介

1 手机电磁辐射的产生及对人体的影响
手机的使用频率为800 - 1800 MHz,当手机天线辐射的微波从空中传播并与活体表面碰撞时,传播条件发生变化。微波反射、折射和吸收生物体的表面,其中一些进入生物体。此时,如果在体内传播的方向合适,就会发生全反射,形成驻波,使机体吸收最大能量。由于手机的微波辐射频率高,波长短,当它作用于活体时,它不是一个连续的场能,而是一个脉冲波能,因此对活体的影响变大。为了使手机随时在线,手机需要与基站接触电磁波,这进一步增加了人体吸收电磁波的时间。
从本质上讲,手机可以被认为是一种低功耗的无线收发器,其中发射的高频无线电波会对人体健康产生不利影响。由于这些移动终端在使用过程中往往接近人脑,因此这是一个值得考虑的问题。长期使用手机会对眼睛、耳朵等重要器官造成一定程度的损害,没有大脑的保护。生活中就有这样的现象。当我们坐在电脑前,当手机响起的时候,会有一条水平线电子液晶屏幕上闪烁,和演讲者也会发出“滴”干扰声音,表明吸收的电磁波是液晶显示器和议长,但人体并不直接的感觉。这种效应称为非热效应。红外线能使人感到发烧,称为热效应。
热效应是指人体在电磁辐射下吸收辐射能,在体内转化为热能,加热局部组织,产生生物效应。在电场的作用下,人体的正电荷和负电荷向相反的方向运动而极化,在极化和定向交替的过程中,碰撞和摩擦产生热量,电解质也存在于体内。内部离子在电场的作用下倾向于运动,当电磁场频率高时,这些离子在平衡位置振动,将电能转化为热能。电磁场也会在体内产生局部涡流,从而产生热量。对于经常使用手机的人来说,非热效应会导致头痛、记忆力减退和潜在的生物损伤。
手机辐射的测量标准主要包括功率密度标准和比吸收率(SAR)标准。本设计将采用比吸收率标准进行研究。

2 时域有限差分法(FDTD)简介
有限差分时域法(FDTD)于1966年由KSYee在AP上发表,后被称为Yee网格空间离散法。这种方法最初是在20世纪60年代引入的,并开始引起关注。1970年代。20世纪80年代以后,它逐渐被接受,并开始在电磁计算领域发展。在这个时代,这种方法已经进入了蓬勃发展的时期。现阶段,关于FDTD方法在解决电磁散射问题、电磁兼容性问题、天线问题、生物电磁问题等方面的应用文章层出不穷。
时域有限差分法的核心思想是将含时间变量的麦克斯韦旋度方程转化为微分形式,模拟电子脉冲的时域响应和理想导体作用。青蛙青蛙算法——电场和磁场中的电场和磁场通过更新时域来模拟电磁场的变化,达到数值计算的目的。这样分析问题时,要考虑研究对象的几何参数、材料参数、计算精度、计算复杂度和计算稳定性。它的优点是可以直接模拟场的分布,具有较高的精度。这是目前数值模拟中常用的方法之一。

3 本设计主要研究内容
本设计首先介绍了移动通信终端对人体的影响、研究方法和测量标准。然后分别详细介绍了时域有限差分(FDTD)理论和SAR标准、仿真和仿真,然后介绍了电磁波对介质的穿透深度。仿真模拟并提出了一种保护手机辐射的方法。最后对设计进行了总结,并提出了相应的防护措施。

二、源代码

clear

%***********************************************************************
%     Fundamental constants
%***********************************************************************

cc=2.99792458e8;            %speed of light in free space光速
muz=4.0*pi*1.0e-7;          %permeability of free space自由空间渗透性
epsz=1.0/(cc*cc*muz);       %permittivity of free space自由空间介电常数

freq=1.0e+9;                %frequency of source excitation 源频率
lambda=cc/freq;             %wavelength of source excitation 源波长
omega=2.0*pi*freq;          %角频率

%***********************************************************************
%     Grid parameters网格参数
%***********************************************************************

ie=200;                     %number of grid cells in x-direction   x方向的网格数

ib=ie+1;

dx=lambda/20.0;             %space increment of 1-D lattice
dt=dx/cc;                   %time step
omegadt=omega*dt;

nmax=round(12.0e-9/dt);     %total number of time steps

%***********************************************************************
%     Material parameters
%***********************************************************************

eps=1.0;
sig=5.0e-3;

%***********************************************************************
%     Updating coefficients for space region with nonpermeable media用非渗透性介质更新空间区域的系数
%***********************************************************************

scfact=dt/muz/dx;

ca=(1.0-(dt*sig)/(2.0*epsz*eps))/(1.0+(dt*sig)/(2.0*epsz*eps));
cb=scfact*(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps));

%***********************************************************************
%     Field arrays
%***********************************************************************

ez(1:ib)=0.0;
hy(1:ie)=0.0;

%***********************************************************************
%     Movie initialization
%***********************************************************************

x=linspace(dx,ie*dx,ie);

subplot(2,1,1),plot(x,ez(1:ie)/scfact,'r'),axis([0 3 -1 1]);
ylabel('EZ');

subplot(2,1,2),plot(x,hy,'b'),axis([0 3 -3.0e-3 3.0e-3]);
xlabel('x (meters)');ylabel('HY');

rect=get(gcf,'Position');
rect(1:2)=[0 0];

M=moviein(nmax/2,gcf,rect);

%***********************************************************************
%     BEGIN TIME-STEPPING LOOP
%***********************************************************************

for n=1:nmax

%***********************************************************************
%     Update electric fields 电场
%***********************************************************************

ez(1)=scfact*sin(omegadt*n);

rbc=ez(ie);
ez(2:ie)=ca*ez(2:ie)+cb*(hy(2:ie)-hy(1:ie-1));
ez(ib)=rbc;

%***********************************************************************
%     Update magnetic fields
%***********************************************************************

hy(1:ie)=hy(1:ie)+ez(2:ib)-ez(1:ie);

%***********************************************************************
%     3-D FDTD code with PML absorbing boundary conditions
%***********************************************************************
function fdtd3D

clear
clc
%***********************************************************************
%     Fundamental constants
%***********************************************************************
cc=2.99792458e8;            %speed of light in free space
muz=4.0*pi*1.0e-7;          %permeability of free space
epsz=1.0/(cc*cc*muz);       %permittivity of free space
etaz=sqrt(muz/epsz);

freq=9e8;
lambda=cc/freq;
omega=2.0*pi*freq;

%***********************************************************************
%     Grid parameters
%***********************************************************************

ie=80;           %number of grid cells in x-direction
if mod(ie,2)==0
    ie=ie+1;
end       %Force it to be a odd integer!
je=80;            %number of grid cells in y-direction
if mod(je,2)==0
    je=je+1;
end       %Force it to be a odd integer!
ke=80;             %number of grid cells in z-direction
if mod(ke,2)==0
    ke=ke+1;
end       %Force it to be a odd integer!

ib=ie+1;
jb=je+1;
kb=ke+1;

is=45;          %location of  hard source in x axis
js=25;          %location of  hard source in y axis
ks=50;          %location of  hard source in z axis


ds=lambda/20;        %space increment of square lattice
dt=ds/(2.0*cc);   %time step

nmax=800;        

iebc=8;           %thickness of left and right PML region
jebc=8;           %thickness of front and back PML region
kebc=8;           %thickness of bottom and top PML region
ibbc=iebc+1;
jbbc=jebc+1;
kbbc=kebc+1;
iefbc=ie+2*iebc;
jefbc=je+2*jebc;
kefbc=ke+2*kebc;
ibfbc=iefbc+1;
jbfbc=jefbc+1;
kbfbc=kefbc+1;

rmax=1.0e-7;
orderbc=2;

%***********************************************************************
%     Material parameters
%***********************************************************************

media=2;

eps=[1.0 1.0];
sig=[0.0 1.0e+7];
mur=[1.0 1.0];
sim=[0.0 0.0];

ds = ds/sqrt(eps(1)*mur(1));

%***********************************************************************
%     Wave excitation
%***********************************************************************

source=zeros(1,nmax);

for n=1:nmax
    source(n)=15.5*sin(omega*n*dt);
end

%***********************************************************************
%     Field arrays
%***********************************************************************

ex=zeros(ie,jb,kb);           %fields in main grid
ey=zeros(ib,je,kb);
ez=zeros(ib,jb,ke);

hx=zeros(ib,je,ke);
hy=zeros(ie,jb,ke);
hz=zeros(ie,je,kb);


exybcf=zeros(iefbc,jebc,kbfbc);   %fields in front PML region
exzbcf=zeros(iefbc,jebc,kbfbc);
eyzbcf=zeros(ibfbc,jebc,kbfbc);
eyxbcf=zeros(ibfbc,jebc,kbfbc);
ezxbcf=zeros(ibfbc,jebc,kefbc);
ezybcf=zeros(ibfbc,jebc,kefbc);

hxybcf=zeros(ibfbc,jebc,kefbc);
hxzbcf=zeros(ibfbc,jebc,kefbc);
hyzbcf=zeros(iefbc,jebc,kefbc);
hyxbcf=zeros(iefbc,jebc,kefbc);
hzxbcf=zeros(iefbc,jebc,kbfbc);
hzybcf=zeros(iefbc,jebc,kbfbc);


exybcb=zeros(iefbc,jebc,kbfbc);   %fields in back PML region
exzbcb=zeros(iefbc,jebc,kbfbc);
eyzbcb=zeros(ibfbc,jebc,kbfbc);
eyxbcb=zeros(ibfbc,jebc,kbfbc);
ezxbcb=zeros(ibfbc,jebc,kefbc);
ezybcb=zeros(ibfbc,jebc,kefbc);

hxybcb=zeros(ibfbc,jebc,kefbc);
hxzbcb=zeros(ibfbc,jebc,kefbc);
hyzbcb=zeros(iefbc,jebc,kefbc);
hyxbcb=zeros(iefbc,jebc,kefbc);
hzxbcb=zeros(iefbc,jebc,kbfbc);
hzybcb=zeros(iefbc,jebc,kbfbc);


exybcl=zeros(iebc,jb,kbfbc);      %fields in left PML region
exzbcl=zeros(iebc,jb,kbfbc);
eyzbcl=zeros(iebc,je,kbfbc);
eyxbcl=zeros(iebc,je,kbfbc);
ezxbcl=zeros(iebc,jb,kefbc);
ezybcl=zeros(iebc,jb,kefbc);

hxybcl=zeros(iebc,je,kefbc);
hxzbcl=zeros(iebc,je,kefbc);
hyzbcl=zeros(iebc,jb,kefbc);
hyxbcl=zeros(iebc,jb,kefbc);
hzxbcl=zeros(iebc,je,kbfbc);
hzybcl=zeros(iebc,je,kbfbc);


exybcr=zeros(iebc,jb,kbfbc);      %fields in right PML region
exzbcr=zeros(iebc,jb,kbfbc);
eyzbcr=zeros(iebc,je,kbfbc);
eyxbcr=zeros(iebc,je,kbfbc);
ezxbcr=zeros(iebc,jb,kefbc);
ezybcr=zeros(iebc,jb,kefbc);

hxybcr=zeros(iebc,je,kefbc);
hxzbcr=zeros(iebc,je,kefbc);
hyzbcr=zeros(iebc,jb,kefbc);
hyxbcr=zeros(iebc,jb,kefbc);
hzxbcr=zeros(iebc,je,kbfbc);
hzybcr=zeros(iebc,je,kbfbc);


exybcd=zeros(ie,jb,kebc);           %fields in bottom PML region
exzbcd=zeros(ie,jb,kebc);
eyzbcd=zeros(ib,je,kebc);
eyxbcd=zeros(ib,je,kebc);
ezxbcd=zeros(ib,jb,kebc);
ezybcd=zeros(ib,jb,kebc);

hxybcd=zeros(ib,je,kebc);
hxzbcd=zeros(ib,je,kebc);
hyzbcd=zeros(ie,jb,kebc);
hyxbcd=zeros(ie,jb,kebc);
hzybcd=zeros(ie,je,kebc);
hzxbcd=zeros(ie,je,kebc);


exybct=zeros(ie,jb,kebc);          %fields in top PML region
exzbct=zeros(ie,jb,kebc);
eyzbct=zeros(ib,je,kebc);
eyxbct=zeros(ib,je,kebc);
ezxbct=zeros(ib,jb,kebc);
ezybct=zeros(ib,jb,kebc);

三、运行结果

在这里插入图片描述

四、备注

版本:2014a