2005阅读
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; Ep^B,;~
        hy=0; vqrBRlZ
        obj_parameters[0]=1; *)j@G:
        obj_parameters[1]=1; "_nX5J9
        obj_parameters[2]=0; 2^zg0!z
        obj_parameters[3]=0; gzl%5`DBw
        CA=ca; S[-.tvI;Q
        CB=cb; ]sX7%3P
        CP=cp; Bv;I0i:_
        CQ=cq; H"Q(2I
    } fl!mYCPv
?,WUJH?^
     X?KGb{
G[|3^O>P
    // 开始主循环 euRCBzc
21.YO]Et
    Nsteps=1; umJay/>
CMC?R,d
    while(Nsteps>0) hWe}'L-
    { {@Blj3;w}
        printf("Nsteps-->"); mcvDxjk,h
        scanf("%d",&Nsteps); pO\S#GnX
w0#%AK
        for(n=0;n<Nsteps;n++) R!sNg
        { Vil@?Y"
            T+=1; Wl,%&H2S<
                         H"2U)HJl
            // 计算 Ex 4qqF v?O[r
            for(i=e_low+1;i<e_high;i++) IetCMp
            { ~wfoK7T}
                                                                 // 采用葛德彪书中的一维FDTD 公式 y' 2<qj
                ex=CA*ex-CB*(hy-hy[i-1])/ddz; WEno+Z~=1'
            } 0v;ve
     F?!FD>L{`
            if(T<200){ n(Qj||:
UP\8w#~
            ex[pos_driv_source]+=gauss_pulse(T,40,12); K3La9O)>
                 8&i;hZm
            } ^%-NPo<
!.9l4@z#
             wG_4$kyj
            //Mur 边界吸收 aLV~|$:2
             .=?Sz*3
            ex_low_m2=ex_low_m1; ?A|zRj{
            ex_low_m1=ex[e_low]; HTxB=Q|
             %(fL?
            ex_low_s2=ex_low_s1; a a4$'8s
            ex_low_s1=ex[e_low+1]; ZQ@3P7T
:RPVT,O}
            ex_high_m2=ex_high_m1; ,Yo: &>As
            ex_high_m1=ex[e_high]; ty':`)
             X>2? `8M
            ex_high_s2=ex_high_s1; mG X\wta
            ex_high_s1=ex[e_high-1]; _\p`4-.V
8a7YHUL<3i
           A8J?A#R*{q
                             ex[e_low]=ex_low_s2+var_border*(ex[e_low+1]-ex_low_m2); $H4=QVj6
                             ex[e_high]=ex_high_s2+var_border*(ex[e_high-1]-ex_high_m2); bC6X?m=
b7Yq_%+
                                          //                行波时延法 SzRL}}I
        //    ex[e_high]=0.5*(ex_high_m2+ex_high_s2); JAN|aCzD
        //    ex[e_low]=0.5*(ex_low_m2+ex_low_s2); ,Ak ^nX
6s'[{Ov
            // 计算 Hy yj>){NcX
HP#ki!'
            for(i=e_low;i<e_high;i++) HTw#U2A;+
            { Lg8]dBXu
                hy=CP*hy-CQ*(ex[i+1]-ex)/ddz; A5+q^t}
            } S45'j(S=
ACgt" M.3F
        } dZF8R
.hxin[Y
        // 存储数据 .a {QA
X@cSP7b
        fp=fopen("Ex","w"); ZHz^S)o\[s
        for(i=0;i<NUM_of_ZAXIS;i++) }&mj.hGv
        { \{lE0j7}h
        fprintf(fp,"%6.2f \n",ex); 4yhcK&
        } ) 9xX
        fclose(fp); " l.!Ed
4AJ9`1d4
        fp=fopen("Hy","w"); ZH% we
        for(i=0;i<NUM_of_ZAXIS;i++) D$ej+s7
        { #iiwD|
        fprintf(fp,"%6.5f \n",hy); 8*vFdoE_oO
        } y!F:m=x<
        fclose(fp); %,XI]+d
-W9gH
        printf("T=%6.2f",T); ATo}FL 2
l4zw]AYk+X
_n7%df
    } ad9EG#mD#
D6Dn&/>Zp
} gUspGsfr
f0OgK<.>T
float gauss_pulse(float T,float t0,float spread) TFYw
{ $<:'!#%
    float pulse; KA?v.s
    pulse=exp(-0.5*(pow((T-t0)/spread,2))); ]/a g*F
    return pulse; YKNb59k
ZxI]I1)
} s88y{o
Zct!/u9 Q
// 运行这个程序，总是很大的反射波，有没有哪位朋友帮我看看程序有没问题啊？

静思而后动

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