2004阅读
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; 0zkT8'v
        hy=0; f%{ ag
        obj_parameters[0]=1; WG!;,~f>o
        obj_parameters[1]=1; -n$rKEC4
        obj_parameters[2]=0; y*TNJJ|
        obj_parameters[3]=0; M,dzf
        CA=ca; 9h{:!
        CB=cb; hghto \G5Y
        CP=cp; [n@!=T
        CQ=cq; =<27qj
    } Yi! >8
*p<5(-J3
     "YJ;-$rb
Hi 0df3t
    // 开始主循环 @R Yb-d
$9u
    Nsteps=1; B7QtB3bn
@):NNbtA
    while(Nsteps>0) 3[g++B."pC
    { 3Tte8]0
        printf("Nsteps-->"); GzE3B';g
        scanf("%d",&Nsteps); 113x9+w[
"u>sS
        for(n=0;n<Nsteps;n++) r:\5/0(
        { kl<B*:RqH
            T+=1; b"3T(#2<*
                         =;?PVAdu%#
            // 计算 Ex 5@r Zm4U
            for(i=e_low+1;i<e_high;i++) U=<.P;+f9
            { E)Epr&9S
                                                                 // 采用葛德彪书中的一维FDTD 公式 #K~j9DuR
                ex=CA*ex-CB*(hy-hy[i-1])/ddz; QcVtv7+*v
            } 7Mbt*[n
     XIW:Nk!S
            if(T<200){ Xcpm?aTo
R;m0eG`
            ex[pos_driv_source]+=gauss_pulse(T,40,12); j`&i4K:
                 E;Ftop
            } _=M'KCL*)
[LEh
             Ej3hdi)
            //Mur 边界吸收 ;{C{V{
             ,(Hmk(,
            ex_low_m2=ex_low_m1; |@9I5Eg)iE
            ex_low_m1=ex[e_low]; ?F), 4Q
             (6y3"cbe
            ex_low_s2=ex_low_s1; zNTu j p
            ex_low_s1=ex[e_low+1]; h+d;`7Z>
\36;csu
            ex_high_m2=ex_high_m1; [";5s&)q
            ex_high_m1=ex[e_high]; :^L]Da3
             x!Y@31!Dy
            ex_high_s2=ex_high_s1; >Ezwl5b
            ex_high_s1=ex[e_high-1]; %Yn)t3d
.*ovIU8
           r4>I?lD
                             ex[e_low]=ex_low_s2+var_border*(ex[e_low+1]-ex_low_m2); sWCm[HpG
                             ex[e_high]=ex_high_s2+var_border*(ex[e_high-1]-ex_high_m2); Q]'!FmXf
~O7(0RsCN
                                          //                行波时延法 %H~gN9Vn#@
        //    ex[e_high]=0.5*(ex_high_m2+ex_high_s2); \@t5S
        //    ex[e_low]=0.5*(ex_low_m2+ex_low_s2); <;Z3 5{
13az[
            // 计算 Hy ]<?)(xz
1^>g>bn_"
            for(i=e_low;i<e_high;i++) `ha:Gf
            { RL)3k8pk
                hy=CP*hy-CQ*(ex[i+1]-ex)/ddz; 'i-O
            } G&H"8REm
.v&h>@'m
        } Oydmq,sVe(
0e](N`
        // 存储数据 1iW9?=a"
U:jf9L2
        fp=fopen("Ex","w"); vj$6
        for(i=0;i<NUM_of_ZAXIS;i++) N9|.D.#MF
        { w[G_w:$a
        fprintf(fp,"%6.2f \n",ex); qP'g}Pc
        } m =F@CA~C
        fclose(fp); Ax*~[$$~%
YYh_lAS>
        fp=fopen("Hy","w"); "kg$s5o
        for(i=0;i<NUM_of_ZAXIS;i++) F}DD;K
        { Xkm2C)
        fprintf(fp,"%6.5f \n",hy); wdV?&W+
        } _7.y4zQJ
        fclose(fp); O;sQPG,v
}]?G"f t K
        printf("T=%6.2f",T); Y"%o\DS*
*?"{T;4u~O
e;[8GE.
    } 3) 0~:
a$O]'}]`
} Z!eq/
r,KK%B
float gauss_pulse(float T,float t0,float spread) W-QPO
{ _?c.3+;s
    float pulse; .)zISa*Xy
    pulse=exp(-0.5*(pow((T-t0)/spread,2))); wO%:WL$5
    return pulse; /CE d14.
=lD]sk
} ,4,c-
[0#hgGO]P
// 运行这个程序，总是很大的反射波，有没有哪位朋友帮我看看程序有没问题啊？

静思而后动

离线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