分享自己编写的一维有限元计算实例“金属衬底介质片对平面波的反射”代码: 本帖最后由 tingo 于 2012-1-4 15:19 编辑
附件中的FEM1D.m、computKeBe.m是求解一维有限元问题的通用接口代码，该实例是在磁场极化下的求解程序，有问题大家可以相互讨论
自己先顶起
{:7_1234:}
支持原创
ernest0512
楼上啥意思？？
本帖最后由 kerbcurb 于 2010-6-1 11:28 编辑
曾经用另外一种方法弄过这个，结果还不错，实际上那个问题就是一个微分方程的边值问题，用有限元练练手还是很不错的，E和H的结果见下面的两个图
曾经用另外一种方法弄过这个，结果还不错，实际上那个问题就是一个微分方程的边值问题，用有限元练练手还是 ...
kerbcurb 发表于 2010-6-1 11:01
这个问题金教授的书上有详细的讲解，但它的解析公式好像有问题，计算的曲线与书上曲线和FEM曲线有很大的出入
{:7_1273:}
本帖最后由 kerbcurb 于 2010-6-2 10:57 编辑
这个问题金教授的书上有详细的讲解，但它的解析公式好像有问题，计算的曲线与书上曲线和FEM曲线有很大 ...
cem-uestc 发表于 2010-6-1 11:57
我发的这两个图上的小方块就是从金建铭书上用软件[GetData]扣下来的数据，图形不是很规整，我感觉我算的结果应该更精确，我使用的是高阶方法。
印象中，书上关于有限元部分说到的公式没有问题，书上介绍的解析法计算公式我没有验证过。
虽然金先生的书看了好几遍，但是从来没有用有限元解决过问题，我是用自己的方法求解的，金德书上若干个例子都求解过，结果符合的不错，最近较忙，一直没有整理。
下面的图形是69页上的那个例子，假设微带线上的电势等于100
{:7_1247:}{:7_1273:}
要时常过来客串一下
我发的这两个图上的小方块就是从金建铭书上用软件[GetData]扣下来的数据，图形不是很规整，我感觉我算的结 ...
kerbcurb 发表于 2010-6-2 09:30
金老师书上的公式计算是对应不上解析解的值，你以前采用什么算法计算？
我发的这两个图上的小方块就是从金建铭书上用软件[GetData]扣下来的数据，图形不是很规整，我感觉我算的结 ...
kerbcurb 发表于 2010-6-2 09:30
感谢分享，你这个微带线问题是用有限差分法做的？
本帖最后由 kerbcurb 于 2010-6-2 16:08 编辑
金建铭书上的公式没有问题，上传很慢啊,下面的附件是一个说明：
感谢分享，你这个微带线问题是用有限差分法做的？
tingo 发表于 2010-6-2 15:24
不是用有限差分，那个精度应该没有这么高
金建铭书上的公式没有问题，上传很慢啊,下面的附件是一个说明：
kerbcurb 发表于 2010-6-2 15:45
您传的附件我看了，有限元求解的公式没问题，但是这个反射系数的理论解析解，就是金老师书中36到37页的那个递推公式求解金属介质反射系数，这儿有错，我编程计算了，没有得到正确的结果。
本帖最后由 kerbcurb 于 2010-6-2 20:13 编辑
我没用那个计算过，按照哪种方法，我感觉应该分的很密才能得到比较精确的结果，你把厚度分成200份试一下，不行就再分的密一些。
本帖最后由 tingo 于 2010-6-4 09:13 编辑
这个是我用解析解计算的结果，当所分的薄片数比较少的时候，其解是不稳定的，当数目足够大的时候就趋向于我给出的那个结果，都是在0.36~0.37之间
为了解决这个问题，今天上午我把那些公式重新推导了一遍，没有问题，我把推导过程写下来了，你看看，或许有帮助
太感谢你了，我看看哈
交流的氛围不是不能获得的！只要用心去交流，才能收获更多！
下午编程序，一直没有得到正确结果，不知道是怎么回事。
其实，所谓解析方法，有多种方案，我给出的那个结果，实际上是解析的，但是我得利用另一个边界条件（x=L时的那个条件），但是这里的方法好像没有利用x=L时的那个条件
下午编程序，一直没有得到正确结果，不知道是怎么回事。
其实，所谓解析方法，有多种方案，我给出的那个结 ...
kerbcurb 发表于 2010-6-4 20:03
我刚看了你的推导，自己也重新推了一下，书中给出的公式没错，但是还是得不到正确的结果，我刚才试着改动了参数，但是结果仍是错的。目前还是不知道问题出在哪了。非常很感谢你亲自做了这么多，我这几天在忙着写论文，等论文写完后，我再看看这个问题。
X=L的条件也要加上，具体见附件，可能是程序有点问题，算出来的结果，只是在一些特殊的角度符合，不知道是怎么回事，推导是否有错误，我自己也检查不出来，你看看，或许有帮助。
这才是交流，这才是社区的价值所在，大力支持，我们也要努力做好这一方面的工作了
回复 24# kerbcurb 的帖子
谢谢你，我先下下来，这几天在赶论文，完了后，好好看看。
金书上的解析法是可以的，有一点是为什么我说不清楚，就是上面的那个文档的最后一个公式中exp(jk0Lcos(sita))需要改成：exp(jk0L)，结果见附件图形，其中红色三角是从书上扣下来的数据，篮圈是做上面介绍的改动以后以后，所谓的解析法的解，黑实线是用我自己的方法解出来的。以下是相关代码：由于这个网站不能正确显示C++代码，代码放在一个pdf文件中
- #include <iostream.h>
- #include <complex.h>
- #include <math.h>
- using namespace std;
- int main()
- {
- double L = 5.0, Lambda = 1.0;
- double x, dx, ds,dt, s, k0 = 2 * M_PI / Lambda ;
- const int N = 200, M = 91;
- complex<double> MIU[N+1], EPS[N+1], K[M][N + 1], R[M][N + 1], YETA, X1, X2,P, P1;
- dx = L / N;
- ds = L / ( 2 * N + 2 );
- dt = 0.5 * M_PI / ( M ‐ 1 );
- for ( int i = 0; i < N; i++ )
- {
- x = ( i + 0.5 ) * dx;
- MIU[i] = complex<double> ( 2 , ‐0.1 );
- EPS[i] = complex<double> ( 4 + 2.0 * pow ( 1 ‐ x / L, 2 ), ‐0.1 * pow ( 1 ‐ x / L, 2 ) );
- }
- MIU[N] = 1.0;
- EPS[N] = 1.0;
- for ( int i = 0; i < M; i++ )
- {
- s = dt * i;
- for ( int j = 0; j < N; j++ )
- {
- K[i][j] = k0 * sqrt ( MIU[j] * EPS[j] ‐ sin ( s ) * sin ( s ) );
- }
- K[i][N] = k0;
- }
- for ( int i = 0; i < M; i++ )
- {
- R[i][0] = ‐1;
- s = dt * i;
- for ( int j = 0; j < N ‐ 1; j++ )
- {
- x = j * dx;
- P = complex<double> ( 0.0, 2.0 ) * K[i][j] * x;
- P1 = complex<double> ( 0.0, 2.0 ) * K[i][j + 1] * x;
- X1 = exp ( P );
- X2 = exp ( P1 );
- YETA = ( K[i][j+ 1] ‐ K[i][j] ) / ( K[i][j+1] + K[i][j] );
- R[i][j + 1] = ( YETA + R[i][j] / X1 ) * X2 / ( 1.0 + YETA * R[i][j] / X1 );
- }
- P = complex<double> ( 0.0, 2.0 ) * K[i][N‐1] * L;
- P1 = complex<double> ( 0.0, 2.0 ) * k0 * L;
- X1 = exp ( P );
- X2 = exp ( P1 );
- R[i][N] = 2 * k0 * ( X1 + R[i][N‐1] ) * cos ( s ) * exp ( complex<double> ( 0, k0 * L ) ) /
- ( ( K[i][N‐1] / MIU[N‐1] + k0 * cos ( s ) ) * X1 ‐ ( K[i][N‐1] / MIU[N‐1] ‐ k0 * cos ( s ) ) *
- R[i][N‐1] ) *
- exp ( complex<double> ( 0, k0 * L ) ) ‐ X2 ;
- cout<<R[i][N]<<endl;
- }
- cin.get();
- return 0;
- }
- chuguozhen@yahoo.com.cn
不知道这个问题是否是类似于FDTD吸收边界的问题，或者吸波材料，不过我感觉，这一类问题如果能从物理的角度去分析，或许能得到更好的解答，比如磁导率的实部虚部分别对应什么物理机制，相对介电常数的实部虚部又分别对应什么物理机制，改变这些参数后，会有什么表现，这些或许更重要。
不知道计算电磁学领域目前是否存在为计算而计算，为了发论文而计算的现象，此为题外话。
不知道这个问题是否是类似于FDTD吸收边界的问题，或者吸波材料，不过我感觉，这一类问题如果能从物理的角度 ...
kerbcurb 发表于 2010-6-9 17:06
这个问题很好，两者都有，呵呵
网站是可以显示各类代码的，只要在发帖回帖编辑器上点击“代码”按钮，将代码复制过去就可以了
你最后一层的反射系数是如何得到？
本帖最后由 kerbcurb 于 2010-6-9 23:14 编辑
回复 31# cem-uestc 的帖子
最后的公式，就是计算介质区域以外的发射率，前面给出的Doc_1那个描述，最后一个公式，需要修改，不过，我自己的方法好像不需要改变那个，我的方法也是解析法
本帖最后由 kerbcurb 于 2010-6-13 11:23 编辑
下面的是Hz模式所谓的计算结果，没有采取金的书上迭代求解R的方法，是求得全部的Am,Bm后，计算R。实际上，分为N段后，每一段有A,B两个未知数，总共2N个未知数，分为N段提供N + 1个匹配点，两个端点【边界点】提供两个方程，其余N - 1个点【连接点】提供2N - 2个方程，因此总共2N个方程，因此可以解出所有的Am,Bm。蓝色圆点是金的书上的数据，红色+号是用这种方法计算的结果，迭代方法速度要快一些，因为没有关于矩阵的计算，不过本方法速度也不慢
多谢分享，谢谢
O(∩_∩)O哈哈~，沙发是我的啦
00d44 于 2010-6-16 11:06 使用 抢沙发 抢夺本帖沙发
学习下有限元
多谢分享了
很好的东西，看看！！
ddddddddddddddddddddddddd
Thanks LZ for your sharing
正好学习这个 FEM
谢谢楼主分享
程序这方面的哦资料不多啊
我是程序菜鸟
辛勤的回复你的劳动成果
太感谢了，我找了好久了。
正在学习中
我也在做，参考下
回复 00d44 的帖子
管理员大哥，拜托你件事，我技术分不够，你能不能把帖子附件发下给我，qq邮箱：68213412.
灰常感谢啊！！！~~~~
回复 tingo 的帖子
我也遇到这个问题，公式好像不对，用手大概算下数量级都对不上么，你当时是怎么解决的。
学习一下
好好学习一下，非常感谢哦！
辛苦楼主，收下了,谢谢了!
顶一下 谢谢分享
谢谢分享
一直没算出来呀，对照下我哪里错了
水泥工 发表于 2011-12-23 16:50
一直没算出来呀，对照下我哪里错了
呵呵，有问题欢迎在此讨论！
看一下源程序
很不错的源程序
Thanks for sharing. I was looking into this problem too.
好贴，可惜现在看不了，有点遗憾
自己编了，得不到正确的结果
泪奔~
yanie1 发表于 2012-2-10 14:06
自己编了，得不到正确的结果
泪奔~
有问题，欢迎在计算电磁学子版块交流！
想看看
我也在研究这个东西 很有用的资料~~
谢谢楼主的分享
有没有金建铭的其他程序呢
有c语言的代码吗
谢谢楼主分享
谢谢楼主分享
谢谢楼主分享
谢谢楼主分享
谢谢楼主分享
谢谢您的分享！
通过学习您的代码，自己在角度的输入和有限元的实现过程中修正下您的计算过程，很感谢您的程序
磁场极化.rar (2.05 KB) 
Doc7.pdf (1.1 MB) 
反射系数求解.bmp
