2006阅读
7回复

[讨论]关于一维FDTD Mur边界的问题

wave wave wave

UID ：55184

注册：2010-03-19
登录：2012-07-04
发帖：148
等级：仿真二级

0楼发表于: 2010-03-19 22:56:35

最近在看葛德彪的《电磁波时域有限差分方法》一书，在看到Mur边界时，自己做的一维Mur边界的吸收效果不好。
用Gauss脉冲激励（最大值为1），在边界边会有0.2的反射波不能被吸收。
程序调了很久，就是不能消除反射波，不知道有没有人做过这个简单的一维FDTD程序，有没有出现这样的情况？

还有在书中提到的一维行波延时法，当C*dt=dz/2时的吸收边界公式为
Ez（K , n+1）=1/2*{ Ez(K , n ) + Ez(K-1 , n ) )
其中 K表示空间节点， n 表示时间节点。用这个公式设置边界也不能完全吸收。
我认为这个公式是简单的线性插值，如果入射到边界处的波形在一个空间节点距离上线性度不好的话，那么就不能很好的吸 ..

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

静思而后动

离线shuangzh

wave wave wave

UID ：55184

注册：2010-03-19
登录：2012-07-04
发帖：148
等级：仿真二级

1楼发表于: 2010-03-20 00:40:38

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define NUM_of_ZAXIS 400

float gauss_pulse(float T,float t0,float spread);

int file_save(float* data,char* filename);

void main()
{
    float ex[NUM_of_ZAXIS],hy[NUM_of_ZAXIS];
    float obj_parameters[NUM_of_ZAXIS][4];  /*模型参数设置 */
    float CA[NUM_of_ZAXIS],CB[NUM_of_ZAXIS],CP[NUM_of_ZAXIS],CQ[NUM_of_ZAXIS];
    float ca,cb,cp,cq;        /*真空时的参数*/
    float var_ca,var_cb,var_cp,var_cq;        /* */
    float var_border,ex_low_m1,ex_low_m2,ex_high_m1,ex_high_m2,    ex_low_s1,ex_low_s2,ex_high_s1,ex_high_s2;    /*边界吸收参数*/
    float Epsilon,Mu,Pi,C;          /*介电系数 Epsilon 0，磁导系数 Mu 0*/
    float rel_epsz,rel_mu;   /*相对介电系数，磁导系数*/
    float e_sigma,h_sigma;   /*电导率，磁导率*/

    float dt,ddz;
    float source,T;
    int pos_driv_source;
    int i,n,Nsteps;
    int e_low,e_high;

    FILE* fp;


    /*初始化各个变量 */
    Pi=3.14159;
    Epsilon=8.85e-12;
    Mu=4*Pi*1e-7;

    C=pow( (float)(1/(Epsilon*Mu)),(float)(0.5));

    ddz=0.01;
    dt=ddz/(2*C);
    T=0;

                    // 采用葛德彪书中一维FDTD公式中的参数
    e_sigma=0;
    h_sigma=0;
    var_ca=0.5*e_sigma*dt/Epsilon;
    var_cb=dt/Epsilon;
    var_cp=0.5*h_sigma*dt/Mu;
    var_cq=dt/Mu;

    ca=(1-var_ca)/(1+var_ca);
    cb=var_cb/(1+var_ca);
    cp=(1-var_cp)/(1+var_cp);
    cq=var_cq/(1+var_cp);

    var_border=(C*dt-ddz)/(C*dt+ddz);

                   // 边界处用来存数据的临时变量
    ex_low_m1=0;
    ex_low_m2=0;
    ex_high_m1=0;
    ex_high_m2=0;
    ex_low_s1=0;
    ex_low_s2=0;
    ex_high_s1=0;
    ex_high_s2=0;

                    // 上，下边界位置
    e_low=20;
    e_high=380;

    pos_driv_source=200;

    for(i=0;i<NUM_of_ZAXIS;i++)
    {
        ex=0; o\=n4;S
        hy=0; zP)~a
        obj_parameters[0]=1; S Xr%kndS
        obj_parameters[1]=1; ,r^"#C0J}
        obj_parameters[2]=0; z5 m>H;P
        obj_parameters[3]=0; p]T"|!d
        CA=ca; ^@6q
        CB=cb; z`3( ,V
        CP=cp; 9A$m$
        CQ=cq; k%81f'H
    } !8RwO%c(
'%;\YD9
     Acm<-de
rf K8q'@
    // 开始主循环 nsuX*C7
N{Qxq>6 G
    Nsteps=1; nE W31 8
pdVQ*=c?M
    while(Nsteps>0) H)(jh
    { Gc,_v3\
        printf("Nsteps-->"); [2c{k
        scanf("%d",&Nsteps); NL"G2[e
UQ?%|y*Kc
        for(n=0;n<Nsteps;n++) dp5cDF}l
        { j<yiNHC
            T+=1; q6d~V]4:
                         T6BFX0$
            // 计算 Ex n6Z|Q@F
            for(i=e_low+1;i<e_high;i++) +cu^%CXT
            { zTm]AG|0
                                                                 // 采用葛德彪书中的一维FDTD 公式 %&<LNEiUN
                ex=CA*ex-CB*(hy-hy[i-1])/ddz; g_.^O$}
            } 8?FueAM'
     5"KlRuv%
            if(T<200){ v3[@1FQ"
iw?I
            ex[pos_driv_source]+=gauss_pulse(T,40,12); 5TKJWO.
                 wIvo"|%
            } ",qU,0
f>$``.O
             N:5[,O<m_
            //Mur 边界吸收 -7qIToO.
             <yUstz,Xu^
            ex_low_m2=ex_low_m1; LV{Q,DrP
            ex_low_m1=ex[e_low]; @* ust>7
             ts~{w;c
            ex_low_s2=ex_low_s1; clO,}Ph>
            ex_low_s1=ex[e_low+1]; ! NV#U
io2)1cE&f
            ex_high_m2=ex_high_m1; 4Ft1@
            ex_high_m1=ex[e_high]; \_6OCVil
             E2wz(,@
            ex_high_s2=ex_high_s1; y(jg#7)
            ex_high_s1=ex[e_high-1]; ~p1EF;4#
u,.3
           p<Z3tD;Z
                             ex[e_low]=ex_low_s2+var_border*(ex[e_low+1]-ex_low_m2); E>jh"|f:{
                             ex[e_high]=ex_high_s2+var_border*(ex[e_high-1]-ex_high_m2); je,}_:7
si4-3eC
                                          //                行波时延法 JH,/jR
        //    ex[e_high]=0.5*(ex_high_m2+ex_high_s2); 7E$&2U^Js
        //    ex[e_low]=0.5*(ex_low_m2+ex_low_s2); wucV_p.E
:a[Ihqfg
            // 计算 Hy g2cVZ!GIj
dtStTT
            for(i=e_low;i<e_high;i++) c7uG9
            { OR6ML-|
                hy=CP*hy-CQ*(ex[i+1]-ex)/ddz; ,~PYt*X4
            } ;F:fM!l=
g\fhp{gWB
        } Sb2v_o
l~:v (R5
        // 存储数据 b,H[I!. %
)}v3q6?_
        fp=fopen("Ex","w"); '_s}o<
        for(i=0;i<NUM_of_ZAXIS;i++) IE~%=/|
        { gwkb!#A
        fprintf(fp,"%6.2f \n",ex); kK>Xrj6
        } ,V] ]:eR
        fclose(fp); y8Xv~4qQW
>B -q@D
        fp=fopen("Hy","w"); 'vV$]/wBF
        for(i=0;i<NUM_of_ZAXIS;i++) o?Nu:&yE
        { 'Ye v}QM
        fprintf(fp,"%6.5f \n",hy); P@}Pk
        } LHCsk{3
        fclose(fp); 5MTgK=c
)+y G+
        printf("T=%6.2f",T); 7v}x?I
Wl"0m1G
bUy,5gk-
    } E|EgB33S
4'pS*v
} zoDZZ%{
oP?YA-#nc
float gauss_pulse(float T,float t0,float spread) R0Ue0pF7
{ M(q'%XL^
    float pulse; L6P1L)
    pulse=exp(-0.5*(pow((T-t0)/spread,2))); ^H'a4G3
    return pulse; 1nhtM
5~ 'Ie<Y_
} m`?MV\^
A1Y7;-D
// 运行这个程序，总是很大的反射波，有没有哪位朋友帮我看看程序有没问题啊？

静思而后动

离线gwzhao

方恨少

UID ：17098

注册：2008-08-24
登录：2019-01-09
发帖：1374
等级：荣誉管理员

2楼发表于: 2010-03-20 19:15:28

对1D 真空，可以完全吸收的啊。

逆流而上

离线shuangzh

wave wave wave

UID ：55184

注册：2010-03-19
登录：2012-07-04
发帖：148
等级：仿真二级

3楼发表于: 2010-03-21 10:25:47

回 2楼(gwzhao) 的帖子

能不能把你做的程序发给我看看。我做的怎么总是吸收不掉。
真空中当C*dt=dz/2 时用行波延迟法 E(k, n)=E(K-1, n-2)   可以在边界上完全吸收
用 Ez（K , n+1）=1/2*{  Ez(K , n ) +  Ez(K-1 , n )  ) 这个公式，我的程序也不能完全吸收。
用书上的Mur 公式 E(k, n+1)=E(k-1, n)  -  (c*dt-dz) / (c*dt + dz) * [ E(k-1, n+1) - E(k, n)  ] 设置的边界条件也不能完全吸收。

静思而后动

离线andymiao

UID ：30660