分享自己编写的一维有限元计算实例“金属衬底介质片对平面波的反射”代码

00d44

交流的氛围不是不能获得的！只要用心去交流，才能收获更多！

kerbcurb

下午编程序，一直没有得到正确结果，不知道是怎么回事。
其实，所谓解析方法，有多种方案，我给出的那个结果，实际上是解析的，但是我得利用另一个边界条件（x=L时的那个条件），但是这里的方法好像没有利用x=L时的那个条件

tingo

下午编程序，一直没有得到正确结果，不知道是怎么回事。
其实，所谓解析方法，有多种方案，我给出的那个结 ...
kerbcurb 发表于 2010-6-4 20:03

                               
登录/注册后可看大图

我刚看了你的推导，自己也重新推了一下，书中给出的公式没错，但是还是得不到正确的结果，我刚才试着改动了参数，但是结果仍是错的。目前还是不知道问题出在哪了。非常很感谢你亲自做了这么多，我这几天在忙着写论文，等论文写完后，我再看看这个问题。

kerbcurb

X=L的条件也要加上，具体见附件，可能是程序有点问题，算出来的结果，只是在一些特殊的角度符合，不知道是怎么回事，推导是否有错误，我自己也检查不出来，你看看，或许有帮助。

00d44

这才是交流，这才是社区的价值所在，大力支持，我们也要努力做好这一方面的工作了

tingo

回复 24# kerbcurb 的帖子


    谢谢你，我先下下来，这几天在赶论文，完了后，好好看看。

kerbcurb

金书上的解析法是可以的，有一点是为什么我说不清楚，就是上面的那个文档的最后一个公式中exp(jk0Lcos(sita))需要改成：exp(jk0L)，结果见附件图形，其中红色三角是从书上扣下来的数据，篮圈是做上面介绍的改动以后以后，所谓的解析法的解，黑实线是用我自己的方法解出来的。以下是相关代码：由于这个网站不能正确显示C++代码，代码放在一个pdf文件中

1. #include <iostream.h>
   
2. #include <complex.h>
   
3. #include <math.h>
   
4. using namespace std;
   
5. int main()
   
6. {
   
7. double L = 5.0, Lambda = 1.0;
   
8. double x, dx, ds,dt, s, k0 = 2 * M_PI / Lambda ;
   
9. const int N = 200, M = 91;
   
10. complex<double> MIU[N+1], EPS[N+1], K[M][N + 1], R[M][N + 1], YETA, X1, X2,P, P1;
    
11. dx = L / N;
    
12. ds = L / ( 2 * N + 2 );
    
13. dt = 0.5 * M_PI / ( M ‐ 1 );
    
14. for ( int i = 0; i < N; i++ )
    
15. {
    
16. x = ( i + 0.5 ) * dx;
    
17. MIU[i] = complex<double> ( 2 , ‐0.1 );
    
18. EPS[i] = complex<double> ( 4 + 2.0 * pow ( 1 ‐ x / L, 2 ), ‐0.1 * pow ( 1 ‐ x / L, 2 ) );
    
19. }
    
20. MIU[N] = 1.0;
    
21. EPS[N] = 1.0;
    
22. for ( int i = 0; i < M; i++ )
    
23. {
    
24. s = dt * i;
    
25. for ( int j = 0; j < N; j++ )
    
26. {
    
27. K[i][j] = k0 * sqrt ( MIU[j] * EPS[j] ‐ sin ( s ) * sin ( s ) );
    
28. }
    
29. K[i][N] = k0;
    
30. }
    
31. for ( int i = 0; i < M; i++ )
    
32. {
    
33. R[i][0] = ‐1;
    
34. s = dt * i;
    
35. for ( int j = 0; j < N ‐ 1; j++ )
    
36. {
    
37. x = j * dx;
    
38. P = complex<double> ( 0.0, 2.0 ) * K[i][j] * x;
    
39. P1 = complex<double> ( 0.0, 2.0 ) * K[i][j + 1] * x;
    
40. X1 = exp ( P );
    
41. X2 = exp ( P1 );
    
42. YETA = ( K[i][j+ 1] ‐ K[i][j] ) / ( K[i][j+1] + K[i][j] );
    
43. R[i][j + 1] = ( YETA + R[i][j] / X1 ) * X2 / ( 1.0 + YETA * R[i][j] / X1 );
    
44. }
    
45. P = complex<double> ( 0.0, 2.0 ) * K[i][N‐1] * L;
    
46. P1 = complex<double> ( 0.0, 2.0 ) * k0 * L;
    
47. X1 = exp ( P );
    
48. X2 = exp ( P1 );
    
49. R[i][N] = 2 * k0 * ( X1 + R[i][N‐1] ) * cos ( s ) * exp ( complex<double> ( 0, k0 * L ) ) /
    
50. ( ( K[i][N‐1] / MIU[N‐1] + k0 * cos ( s ) ) * X1 ‐ ( K[i][N‐1] / MIU[N‐1] ‐ k0 * cos ( s ) ) *
    
51. R[i][N‐1] ) *
    
52. exp ( complex<double> ( 0, k0 * L ) ) ‐ X2 ;
    
53. cout<<R[i][N]<<endl;
    
54. }
    
55. cin.get();
    
56. return 0;
    
57. }
    
58. chuguozhen@yahoo.com.cn

复制代码

kerbcurb

不知道这个问题是否是类似于FDTD吸收边界的问题，或者吸波材料，不过我感觉，这一类问题如果能从物理的角度去分析，或许能得到更好的解答，比如磁导率的实部虚部分别对应什么物理机制，相对介电常数的实部虚部又分别对应什么物理机制，改变这些参数后，会有什么表现，这些或许更重要。

不知道计算电磁学领域目前是否存在为计算而计算，为了发论文而计算的现象，此为题外话。

cem-uestc

不知道这个问题是否是类似于FDTD吸收边界的问题，或者吸波材料，不过我感觉，这一类问题如果能从物理的角度 ...
kerbcurb 发表于 2010-6-9 17:06

                               
登录/注册后可看大图

这个问题很好，两者都有，呵呵

00d44

网站是可以显示各类代码的，只要在发帖回帖编辑器上点击“代码”按钮，将代码复制过去就可以了