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

一、获取代码方式

获取代码方式1：
完整代码已上传我的资源：【通信】基于matlab FDTD法研究移动通信终端电磁辐射对人体的影响【含Matlab源码 761期】

获取代码方式2：
通过订阅紫极神光博客付费专栏，凭支付凭证，私信博主，可获得此代码。

备注：
订阅紫极神光博客付费专栏，可免费获得1份代码（有效期为订阅日起，三天内有效）；

二、简介

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);

  

 

  
1
  
2
  
3
  
4
  
5
  
6
  
7
  
8
  
9
  
10
  
11
  
12
  
13
  
14
  
15
  
16
  
17
  
18
  
19
  
20
  
21
  
22
  
23
  
24
  
25
  
26
  
27
  
28
  
29
  
30
  
31
  
32
  
33
  
34
  
35
  
36
  
37
  
38
  
39
  
40
  
41
  
42
  
43
  
44
  
45
  
46
  
47
  
48
  
49
  
50
  
51
  
52
  
53
  
54
  
55
  
56
  
57
  
58
  
59
  
60
  
61
  
62
  
63
  
64
  
65
  
66
  
67
  
68
  
69
  
70
  
71
  
72
  
73
  
74
  
75
  
76
  
77
  
78
  
79
  
80
  
81
  
82
  
83
  
84
  
85
  
86
  
87
  
88
  
89
  
90
  
91
  
92
  
93
  
94
  
95
  
96
  
97
  
98
  
99
  
100
  
101
  
102
  
103
  
104
  
105
  
106
  
107
  
108
  
109
  
110
  
111
  
112
  
113
  
114
  
115
  
116
  
117
  
118
  
119
  
120
  
121
  
122
  
123
  
124
  
125
  
126
  
127
  
128
  
129
  
130
  
131
  
132
  
133
  
134
  
135
  
136
  
137
  
138
  
139
  
140
  
141
  
142
  
143
  
144
  
145
  
146
  
147
  
148
  
149
  
150
  
151
  
152
  
153
  
154
  
155
  
156
  
157
  
158
  
159
  
160
  
161
  
162
  
163
  
164
  
165
  
166
  
167
  
168
  
169
  
170
  
171
  
172
  
173
  
174
  
175
  
176
  
177
  
178
  
179
  
180
  
181
  
182
  
183
  
184
  
185
  
186
  
187
  
188
  
189
  
190
  
191
  
192
  
193
  
194
  
195
  
196
  
197
  
198
  
199
  
200
  
201
  
202
  
203
  
204
  
205
  
206
  
207
  
208
  
209
  
210
  
211
  
212
  
213
  
214
  
215
  
216
  
217
  
218
  
219
  
220
  
221
  
222
  
223
  
224
  
225
  
226
  
227
  
228
  
229
  
230
  
231
  
232
  
233
  
234
  
235
  
236
  
237
  
238
  
239
  
240
  
241
  
242
  
243
  
244
  
245
  
246
  
247
  
248
  
249
  
250
  
251
  
252
  
253
  
254
  
255
  
256
  
257
  
258
  
259
  
260
  
261
  
262
  
263
  
264
  
265
  
266
  
267
  
268
  
269
  
270
  
271
  
272
  
273
 

四、运行结果

四、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 沈再阳.精通MATLAB信号处理[M].清华大学出版社，2015.
[2]高宝建,彭进业,王琳,潘建寿.信号与系统——使用MATLAB分析与实现[M].清华大学出版社，2020.
[3]王文光,魏少明,任欣.信号处理与系统分析的MATLAB实现[M].电子工业出版社，2018.

文章来源: qq912100926.blog.csdn.net，作者：海神之光，版权归原作者所有，如需转载，请联系作者。

原文链接：qq912100926.blog.csdn.net/article/details/115610809