网站软件设计,wordpress 添加meta,个人做的网站,360浏览器网址金属衬底介质片对平面波的反射-问题的解析求解和FEM求解
参考有限元从零单排系列4 代码参考了上面大佬文章提供的#xff0c;但是部分计算系数错了#xff0c;我改了下加了许多注释#xff0c;便于大家理解。
书籍参考的电磁场有限元方法(金建铭)#xff0c;所用的公式都…金属衬底介质片对平面波的反射-问题的解析求解和FEM求解
参考有限元从零单排系列4 代码参考了上面大佬文章提供的但是部分计算系数错了我改了下加了许多注释便于大家理解。
书籍参考的电磁场有限元方法(金建铭)所用的公式都是这本书的。
下载地址金属衬底介质片对平面波的反射-问题的解析求解和FEM求解-MATLAB代码
QAQ QAQ QAQ
0、问题定义
一个均匀平面波斜入射从自由空间入射到电介质中电介质中的 ε r \varepsilon_{r} εr和 μ r \mu_{r} μr随着其深度 L L L而变化电介质材料的衬底为理想金属PEC。
入射波沿着 θ \theta θ角入射并发生反射我们需要求解其反射系数。
其中介电常数是非均匀的其表达式为 ϵ r 4 ( 2 I − j 0.1 ) ( 1 − x / L ) 2 \epsilon^{\mathrm{r}}4(2^{\mathrm{I}}-\mathrm{j}0.1)(1-x/L)^{2} ϵr4(2I−j0.1)(1−x/L)2 磁导率为固定的复数 μ r 2 − j 0.1 \mu_{\mathrm{r}}2-\mathrm{j}0.1 μr2−j0.1 介质板的厚度为5个波长 L 5 λ L5 {\lambda} L5λ
1、解析求解方法
2、FEM求解方法
2.1、FEM求解目标方程
任意电磁平面波都能分解成只有 E z E_{z} Ez分量的电场极化波和只有 H z H_{z} Hz分量的磁场极化波因此这个问题只需要分别考虑这两种情况即可。对于极化波入射波的电场强度 E z E_{z} Ez表达式为金-式3.81 E z i n c ( x , y ) E 0 e j k 0 x cos θ − j k 0 y sin θ E_{z}^{\mathrm{inc}}(x,y)E_{0}\mathrm{e}^{\mathrm{j}k_{0}x\cos\theta-\mathrm{j}k_{0}y\sin\theta} Ezinc(x,y)E0ejk0xcosθ−jk0ysinθ 其中 E 0 {E_{0}} E0为入射场大小。
此处实际要求解方程为金-式3.82 d d x ( 1 μ r d E z d x ) k 0 2 ( ϵ r − 1 μ r sin 2 θ ) E z 0 {\frac{\mathrm{d}}{\mathrm{d}x}}\left({\frac{1}{\mu_{\mathrm{r}}}}{\frac{\mathrm{d}E_{z}}{\mathrm{d}x}}\right)k_{0}^{2}\left(\epsilon_{\mathrm{r}}-{\frac{1}{\mu_{\mathrm{r}}}}\sin^{2}\theta\right)E_{z}0 dxd(μr1dxdEz)k02(ϵr−μr1sin2θ)Ez0
推导得到的边界条件为 [ d E z d x j k 0 cos θ E z ( x ) ] x L 0 2 j k 0 cos θ E 0 e j k 0 L cos θ \left[{\frac{\mathrm{d}E_{z}}{\mathrm{d}x}}\mathrm{j}k_{0}\cos\theta E_{z}(x)\right]_{xL0}2\mathrm{j}k_{0}\cos\theta E_{0}\mathrm{e}^{\mathrm{j}k_{0}L\cos\theta} [dxdEzjk0cosθEz(x)]xL02jk0cosθE0ejk0Lcosθ
2.2、FEM求解Matlab编程
2.2.1 基础参数定义
定义波长随意设置对最终结果没有影响自由空间波数计算公式基板的厚度为5倍波长-这是给定的条件有限元剖分-100个单元单元节点坐标x1到xm求解对应角度下的反射数据E0为入射场大小对反射计算结果无影响
%基本物理参数
lambda1;%定义波长随意设置对最终结果没有影响
k02*pi/lambda;%自由空间波数计算公式
L5*lambda;%基板的厚度为5倍波长-这是给定的条件
element_num100;%有限元剖分-100个单元
xlinspace(0,L,element_num1);%单元节点坐标x1到xm
thetalinspace(0,pi/2,20);%求解对应角度下的反射数据
% theta[0,pi/4];%求解45度下的反射数据
E01;%E0为入射场大小对反射计算结果无影响2.2.2 材料信息定义
依旧题目条件定义的材料信息如下所示
介质板内材料的非均匀介电常数介质板外自由空间介电常数介质板内材料的磁导率介质板外自由空间磁导率
%材料参数
epsilon_r4(2-0.1*1j)*(1-x/L).^2;%介质板内材料的非均匀介电常数
epsilon_r(element_num1)1;%介质板外自由空间介电常数
mu_rones(1,length(epsilon_r))*(2-0.1j);%介质板内材料的磁导率
mu_r(element_num1)1;%介质板外自由空间磁导率2.2.3 FEM节点信息定义
下面需要定义fem节点信息其实际上是个结构体主要包含下面信息 其中
M为节点的标记若使用100个单元则M的范围为1-100alpha为金-式3.1中的 α \alpha α如下所示也就是求解的通用微分方程的已知系数beta为金-式3.1中的 β \beta β如下所示也就是求解的通用微分方程的已知系数f为金-式3.1中的 f f f如下所示也就是求解的通用微分方程的激励项l是有限元节点之间的距离
alpha、 beta、f对应金-式3.1 − d d x ( α d ϕ d x ) β ϕ f 0 x L -{\frac{\mathrm{d}}{\mathrm{d}x}}\left(\alpha{\frac{\mathrm{d}\phi}{\mathrm{d}x}}\right)\beta\phif\quad0xL −dxd(αdxdϕ)βϕf0xL 代码
%存放fem的所有节点信息(结构体信息)
fems[];
%存放计算得到的反射系数信息
plotData[];2.2.4 节点数组赋值
代码如下
fems[];%清空fems节点数组
% fem节点数组赋值并和fems合并为全局节点数组
for j1:length(x)-1str[fem,num2str(j),inputElementData(j,1/mu_r(j),-k0^2*(epsilon_r(j)-1/mu_r(j)*sin(theta(i))^2),0,x(j),x(j1));];eval(str);str2[fems[fems;fem,num2str(j),];];eval(str2);
end其中inputElementData函数主要给每个fem节点赋值
%% 生成单个有限元单元的元胞(单元编号alphabetaf单元区间)
function feminputElementData(M,alpha,beta,f,x1,x2)fem.MM;fem.alphaalpha;fem.betabeta;fem.ff;fem.lx2-x1;
end对比下面两个方程alpha、beta、f的表达式非常清楚了可以看到和代码完全对应金-式3.1和式3.82 − d d x ( α d ϕ d x ) β ϕ f 0 x L -{\frac{\mathrm{d}}{\mathrm{d}x}}\left(\alpha{\frac{\mathrm{d}\phi}{\mathrm{d}x}}\right)\beta\phif\quad0xL −dxd(αdxdϕ)βϕf0xL d d x ( 1 μ r d E z d x ) k 0 2 ( ϵ r − 1 μ r sin 2 θ ) E z 0 {\frac{\mathrm{d}}{\mathrm{d}x}}\left({\frac{1}{\mu_{\mathrm{r}}}}{\frac{\mathrm{d}E_{z}}{\mathrm{d}x}}\right)k_{0}^{2}\left(\epsilon_{\mathrm{r}}-{\frac{1}{\mu_{\mathrm{r}}}}\sin^{2}\theta\right)E_{z}0 dxd(μr1dxdEz)k02(ϵr−μr1sin2θ)Ez0
2.2.5 初始化边界条件
先给出代码解释如下
boundary1inputBoundaryData(d,min,0);
boundary2inputBoundaryData(n,max,1j*k0*cos(theta(i)),2j*k0*cos(theta(i))*E0*exp(1j*k0*L*cos(theta(i))));其中
%% 边界条件d为狄利克雷n为纽曼
function boundaryContentinputBoundaryData(type,value,varargin)if type dboundaryContent.typetype;boundaryContent.valuevalue;boundaryContent.pvarargin{1};boundaryContent.gamma0;boundaryContent.q0;elseif type nboundaryContent.typetype;boundaryContent.valuevalue;boundaryContent.p0;boundaryContent.gammavarargin{1};boundaryContent.qvarargin{2};elseerror(unknown boundary condition!);end
end显而易见我们求解的是电场最左侧即x轴最小min处是PEC边界对应的电场强度为0因此左侧边界设置为狄利克雷边界金-式3.2
inputBoundaryData(d,min,0)最右侧即x轴最大max处是纽曼边界金-式3.3和式3.94,对比下面两个式子此问题的边界条件和标准式3.3的格式gamma和q的值也非常显而易见这也和上面的编程实现一致 [ α d ϕ d x γ ϕ ] x L q \left[\alpha{\frac{\mathrm{d}\phi}{\mathrm{d}x}}\gamma\phi\right]_{xL}q [αdxdϕγϕ]xLq [ d E z d x j k 0 cos θ E z ( x ) ] x L 0 2 j k 0 cos θ E 0 e j k 0 L cos θ \left[{\frac{\mathrm{d}E_{z}}{\mathrm{d}x}}\mathrm{j}k_{0}\cos\theta E_{z}(x)\right]_{xL0}2\mathrm{j}k_{0}\cos\theta E_{0}\mathrm{e}^{\mathrm{j}k_{0}L\cos\theta} [dxdEzjk0cosθEz(x)]xL02jk0cosθE0ejk0Lcosθ
2.2.6 依据FEM节点生成求解矩阵
金-式3.38附近实际的求解的矩阵如下所示 [ K ] { ϕ } { b } [K]\{\phi\}\{b\} [K]{ϕ}{b} K [ K 11 ( 1 ) K 12 ( 1 ) 0 0 K 21 ( 1 ) K 22 ( 1 ) K 11 ( 2 ) K 12 ( 2 ) 0 0 K 21 ( 2 ) K 22 ( 2 ) K 11 ( 3 ) K 12 ( 3 ) 0 0 K 21 ( 3 ) K 22 ( 3 ) ] K\begin{bmatrix}K_{11}^{(1)}K_{12}^{(1)}00\\K_{21}^{(1)}K_{22}^{(1)}K_{11}^{(2)}K_{12}^{(2)}0\\0K_{21}^{(2)}K_{22}^{(2)}K_{11}^{(3)}K_{12}^{(3)}\\00K_{21}^{(3)}K_{22}^{(3)}\end{bmatrix} K K11(1)K21(1)00K12(1)K22(1)K11(2)K21(2)00K12(2)K22(2)K11(3)K21(3)00K12(3)K22(3) { b } { b 1 ( 1 ) b 2 ( 1 ) b 1 ( 2 ) b 2 ( 2 ) b 1 ( 3 ) b 2 ( 3 ) } \left.\{b\}\left\{\begin{array}{c}{b_{1}^{(1)}}\\{b_{2}^{(1)}b_{1}^{(2)}}\\{b_{2}^{(2)}b_{1}^{(3)}}\\{b_{2}^{(3)}}\\\end{array}\right.\right\} {b}⎩ ⎨ ⎧b1(1)b2(1)b1(2)b2(2)b1(3)b2(3)⎭ ⎬ ⎫
[a,b,c]createMatrixElements(fems);基于上面给出的矩阵我们下一步是计算其中的元素计算公式如下金-式3.28~3.30 K 11 e K 22 e α e l e β e l e 3 K 12 e K 21 e − α e l e β e l e 6 b 1 e b 2 e f e l e 2 \begin{aligned}K_{11}^{e}K_{22}^{e}\frac{\alpha^{e}}{l^{e}}\beta^{e}\frac{l^{e}}{3}\\K_{12}^{e}K_{21}^{e}-\frac{\alpha^{e}}{l^{e}}\beta^{e}\frac{l^{e}}{6}\\b_{1}^{e}b_{2}^{e}f^{e}\frac{l^{e}}{2}\end{aligned} K11eK22eleαeβe3leK12eK21e−leαeβe6leb1eb2efe2le 由此写出的代码如下其中a为对角上的元素值c为次对角元素的值b就对应上面的非齐次项
%% data里面是元胞数组每个元胞里面都有 alphabetaf和l
function [a,b,c]createMatrixElements(data)sizessize(data); numsizes(1);%有限元单元数eleNumnum1;%有限元单元数1为节点数%% 生成K矩阵对角元,非对角元和非齐次项的过程azeros(eleNum,1);bzeros(eleNum,1);czeros(num,1);a(1)data(1).alpha/data(1).ldata(1).beta*data(1).l/3;%对角元第一个的公式3.28a(eleNum)data(num).alpha/data(num).ldata(num).beta*data(num).l/3;%对角元最后一个的公式c(1)-data(1).alpha/data(1).ldata(1).beta*data(1).l/6;%次对角线第一个的公式3.29b(1)data(1).f*data(1).l/2; %非齐次项第一项-f为已知的源或者激励b(eleNum)data(num).f*data(num).l/2; %非齐次项最后一项3.30for j2:eleNum-1a(j)data(j-1).alpha/data(j-1).ldata(j-1).beta*data(j-1).l/3data(j).alpha/data(j).ldata(j).beta*data(j).l/3;%中间的公式3.38b(j)data(j-1).f*data(j-1).l/2data(j).f*data(j).l/2;%非齐次项剩余项3.41c(j)-data(j).alpha/data(j).ldata(j).beta*data(j).l/6;%非对角元的公式3.29end
end2.2.7 强加边界条件
加入边界条件会对上面的k矩阵和b矩阵进行修正主要在下面代码进行
[K,b]createMatrixK([boundary1,boundary2],a,b,c);%强加边界条件在此函数中先把之前的abc数组组合成实际的k矩阵
sizessize(a);
numsizes(1);
Kzeros(num);
for j1:numK(j,j)a(j);
end
for j1:num-1K(j,j1)c(j);K(j1,j)c(j);
end对于狄利克雷边界条件实际的代码会根据下式进行修正金-式3.62 [ 1 0 0 0 0 K 22 K 23 K 24 0 K 32 K 33 K 34 0 K 42 K 43 K 44 ] { ϕ 1 ϕ 2 ϕ 3 ϕ 4 } { p b 2 − K 21 p b 3 − K 31 p b 4 − K 41 p } \left.\left[\begin{array}{cccc}1000\\0K_{22}K_{23}K_{24}\\0K_{32}K_{33}K_{34}\\0K_{42}K_{43}K_{44}\end{array}\right.\right]\begin{Bmatrix}\phi_1\\\phi_2\\\phi_3\\\phi_4\end{Bmatrix}\begin{Bmatrix}p\\b_2-K_{21}p\\b_3-K_{31}p\\b_4-K_{41}p\end{Bmatrix} 10000K22K32K420K23K33K430K24K34K44 ⎩ ⎨ ⎧ϕ1ϕ2ϕ3ϕ4⎭ ⎬ ⎫⎩ ⎨ ⎧pb2−K21pb3−K31pb4−K41p⎭ ⎬ ⎫
%% 狄利克雷边界条件
if boundaryCondition1.typedK(1,1)1;b(1)boundaryCondition1.p;for j2:numK(1,j)0;b(j)b(j)-b(1)*K(j,1);K(j,1)0;end
end对于纽曼边界条件原方程组使用下式进行修正只在边界上对应的元素进行修正金-式3.58~3.59 K N N α ( M ) l ( M ) β ( M ) l ( M ) 3 γ K_{NN}\frac{\alpha^{(M)}}{l^{(M)}}\beta^{(M)}\frac{l^{(M)}}{3}\gamma KNNl(M)α(M)β(M)3l(M)γ b N f ( M ) l ( M ) 2 q b_{N}f^{(M)}\frac{l^{(M)}}{2}q bNf(M)2l(M)q
%% 纽曼边界条件
if boundaryCondition2.typenK(num,num)K(num,num)boundaryCondition2.gamma;b(num)b(num)boundaryCondition2.q;
end2.2.8 方程组求解
还在等什么线性方程组不会求解嘛此处用高斯消元法求解
resultssolveWithGuassianMethod(K,b);%高斯消元法求解function resultssolveWithGuassianMethod(K,b)%生成a和c向量sizessize(K);numsizes(1);azeros(num,1);czeros(num-1,1);for j1:numa(j)K(j,j);endfor j1:num-1c(j)K(j,j1);endfor j2:numa(j)a(j)-c(j-1)^2/a(j-1);b(j)b(j)-b(j-1)*c(j-1)/a(j-1);endresultszeros(num,1);results(num)b(num)/a(num);for jnum-1:-1:1results(j)(b(j)-c(j)*results(j1))/a(j);end
end2.2.9 求解反射系数
反射系数就是下面的R金-式3.94 E z ( x ) E 0 e j k 0 x cos θ R E 0 e − j k 0 x cos θ x L E_{z}(x)E_{0}\mathrm{e}^{\mathrm{j}k_{0}x\cos\theta}RE_{0}\mathrm{e}^{-\mathrm{j}k_{0}x\cos\theta}\quad xL Ez(x)E0ejk0xcosθRE0e−jk0xcosθxL 我们求解得到results的是在各个节点的电场强度将位于介质与空气交界处的电场强度带入上式中计算即可得到反射系数公式如下
R(results(element_num1)-E0*exp(1j*k0*L*cos(theta(i))))/E0*exp(-1j*k0*L*cos(theta(i)));2.3 求解结果
杠杠的
plot(plotData(:,1)*180/pi,abs(plotData(:,2)))2.4 主函数代码
其余从最上面链接下载吧
%非均匀金属衬底反射问题有限元解法
clear all;
clc;
%基本物理参数
lambda1;%定义波长随意设置对最终结果没有影响
k02*pi/lambda;%自由空间波数计算公式
L5*lambda;%基板的厚度为5倍波长-这是给定的条件
element_num100;%有限元剖分-100个单元
xlinspace(0,L,element_num1);%单元节点坐标x1到xm
thetalinspace(0,pi/2,20);%求解对应角度下的反射数据
% theta[0,pi/4];%求解45度下的反射数据
E01;%E0为入射场大小对反射计算结果无影响
%材料参数
epsilon_r4(2-0.1*1j)*(1-x/L).^2;%介质板内材料的非均匀介电常数
epsilon_r(element_num1)1;%介质板外自由空间介电常数
mu_rones(1,length(epsilon_r))*(2-0.1j);%介质板内材料的磁导率
mu_r(element_num1)1;%介质板外自由空间磁导率%存放fem的节点信息(结构体信息)
fems[];
%存放计算得到的反射系数信息
plotData[];for i1:length(theta)%对每个theta进行一次有限元求解fems[];%清空fems节点数组% fem节点数组赋值并和fems合并为全局节点数组for j1:length(x)-1str[fem,num2str(j),inputElementData(j,1/mu_r(j),-k0^2*(epsilon_r(j)-1/mu_r(j)*sin(theta(i))^2),0,x(j),x(j1));];eval(str);str2[fems[fems;fem,num2str(j),];];eval(str2);endboundary1inputBoundaryData(d,min,0);boundary2inputBoundaryData(n,max,1j*k0*cos(theta(i)),2j*k0*cos(theta(i))*E0*exp(1j*k0*L*cos(theta(i))));[a,b,c]createMatrixElements(fems);%依据FEM节点生成求解矩阵[K,b]createMatrixK([boundary1,boundary2],a,b,c);%强加边界条件resultssolveWithGuassianMethod(K,b);%高斯消元法求解R(results(element_num1)-E0*exp(1j*k0*L*cos(theta(i))))/E0*exp(-1j*k0*L*cos(theta(i)));plotData[plotData;[theta(i),R]];
end
plot(plotData(:,1)*180/pi,abs(plotData(:,2)))
