当前位置：首页 > news >正文

金属衬底介质片对平面波的反射-问题的解析求解和FEM求解

news 2025/9/23 2:05:44

金属衬底介质片对平面波的反射-问题的解析求解和FEM求解

参考有限元从零单排系列4
代码参考了上面大佬文章提供的，但是部分计算系数错了，我改了下加了许多注释，便于大家理解。

书籍参考的电磁场有限元方法(金建铭)，所用的公式都是这本书的。

下载地址：金属衬底介质片对平面波的反射-问题的解析求解和FEM求解-MATLAB代码

QAQ
QAQ
QAQ

0、问题定义

一个均匀平面波斜入射从自由空间入射到电介质中，电介质中的 $\varepsilon_{r}$ 和 $\mu_{r}$ 随着其深度 $L$ 而变化，电介质材料的衬底为理想金属PEC。

入射波沿着 $\theta$ 角入射，并发生反射，我们需要求解其反射系数。

其中介电常数是非均匀的，其表达式为：
$\epsilon^{\mathrm{r}}=4+(2^{\mathrm{I}}-\mathrm{j}0.1)(1-x/L)^{2}$
磁导率为固定的复数：
$\mu_{\mathrm{r}}=2-\mathrm{j}0.1$
介质板的厚度为5个波长：
${\lambda}$
在这里插入图片描述

1、解析求解方法

2、FEM求解方法

2.1、FEM求解目标方程

任意电磁平面波都能分解成只有 $E_{z}$ 分量的电场极化波和只有 $H_{z}$ 分量的磁场极化波，因此这个问题只需要分别考虑这两种情况即可。对于极化波，入射波的电场强度 $E_{z}$ 表达式为（金-式3.81）：
$E_{z}^{\mathrm{inc}}(x,y)=E_{0}\mathrm{e}^{\mathrm{j}k_{0}x\cos\theta-\mathrm{j}k_{0}y\sin\theta}$
其中 ${E_{0}}$ 为入射场大小。

此处实际要求解方程为（金-式3.82）：
${\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$

推导得到的边界条件为：
$\left[{\frac{\mathrm{d}E_{z}}{\mathrm{d}x}}+\mathrm{j}k_{0}\cos\theta E_{z}(x)\right]_{x=L+0}=2\mathrm{j}k_{0}\cos\theta E_{0}\mathrm{e}^{\mathrm{j}k_{0}L\cos\theta}$

2.2、FEM求解Matlab编程

2.2.1 基础参数定义

定义波长，随意设置对最终结果没有影响
自由空间波数计算公式
基板的厚度为5倍波长-这是给定的条件
有限元剖分-100个单元
单元节点坐标x1到xm
求解对应角度下的反射数据
E0为入射场大小，对反射计算结果无影响

%基本物理参数
lambda=1;%定义波长，随意设置对最终结果没有影响
k0=2*pi/lambda;%自由空间波数计算公式
L=5*lambda;%基板的厚度为5倍波长-这是给定的条件
element_num=100;%有限元剖分-100个单元
x=linspace(0,L,element_num+1);%单元节点坐标x1到xm
theta=linspace(0,pi/2,20);%求解对应角度下的反射数据
% theta=[0,pi/4];%求解45度下的反射数据
E0=1;%E0为入射场大小，对反射计算结果无影响

2.2.2 材料信息定义

依旧题目条件定义的材料信息如下所示：

介质板内材料的非均匀介电常数
介质板外自由空间介电常数
介质板内材料的磁导率
介质板外自由空间磁导率

%材料参数
epsilon_r=4+(2-0.1*1j)*(1-x/L).^2;%介质板内材料的非均匀介电常数
epsilon_r(element_num+1)=1;%介质板外自由空间介电常数
mu_r=ones(1,length(epsilon_r))*(2-0.1j);%介质板内材料的磁导率
mu_r(element_num+1)=1;%介质板外自由空间磁导率

2.2.3 FEM节点信息定义

下面需要定义fem节点信息，其实际上是个结构体，主要包含下面信息：
在这里插入图片描述
其中：

M为节点的标记，若使用100个单元，则M的范围为1-100
alpha为（金-式3.1）中的 $\alpha$ （如下所示），也就是求解的通用微分方程的已知系数
beta为（金-式3.1）中的 $\beta$ （如下所示），也就是求解的通用微分方程的已知系数
f为（金-式3.1）中的 $f$ （如下所示），也就是求解的通用微分方程的激励项
l是有限元节点之间的距离

alpha、 beta、f对应（金-式3.1）
$-{\frac{\mathrm{d}}{\mathrm{d}x}}\left(\alpha{\frac{\mathrm{d}\phi}{\mathrm{d}x}}\right)+\beta\phi=f\quad0<x<L$
代码：

%存放fem的所有节点信息(结构体信息)
fems=[];
%存放计算得到的反射系数信息
plotData=[];

2.2.4 节点数组赋值

代码如下：

fems=[];%清空fems节点数组
% fem节点数组赋值，并和fems合并为全局节点数组
for j=1: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(j+1));'];eval(str);str2=['fems=[fems;fem',num2str(j),'];'];eval(str2);
end

其中inputElementData函数主要给每个fem节点赋值：

%% 生成单个有限元单元的元胞(单元编号，alpha，beta，f，单元区间)
function fem=inputElementData(M,alpha,beta,f,x1,x2)fem.M=M;fem.alpha=alpha;fem.beta=beta;fem.f=f;fem.l=x2-x1;
end

对比下面两个方程，alpha、beta、f的表达式非常清楚了，可以看到和代码完全对应（金-式3.1和式3.82）：
$-{\frac{\mathrm{d}}{\mathrm{d}x}}\left(\alpha{\frac{\mathrm{d}\phi}{\mathrm{d}x}}\right)+\beta\phi=f\quad0<x<L$
${\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$

2.2.5 初始化边界条件

先给出代码，解释如下：

boundary1=inputBoundaryData('d','min',0);
boundary2=inputBoundaryData('n','max',1j*k0*cos(theta(i)),2j*k0*cos(theta(i))*E0*exp(1j*k0*L*cos(theta(i))));

其中：

%% 边界条件，d为狄利克雷，n为纽曼
function boundaryContent=inputBoundaryData(type,value,varargin)if type == 'd'boundaryContent.type=type;boundaryContent.value=value;boundaryContent.p=varargin{1};boundaryContent.gamma=0;boundaryContent.q=0;elseif type == 'n'boundaryContent.type=type;boundaryContent.value=value;boundaryContent.p=0;boundaryContent.gamma=varargin{1};boundaryContent.q=varargin{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的值也非常显而易见，这也和上面的编程实现一致：
$\left[\alpha{\frac{\mathrm{d}\phi}{\mathrm{d}x}}+\gamma\phi\right]_{x=L}=q$

$\left[{\frac{\mathrm{d}E_{z}}{\mathrm{d}x}}+\mathrm{j}k_{0}\cos\theta E_{z}(x)\right]_{x=L+0}=2\mathrm{j}k_{0}\cos\theta E_{0}\mathrm{e}^{\mathrm{j}k_{0}L\cos\theta}$

2.2.6 依据FEM节点生成求解矩阵

（金-式3.38）附近，实际的求解的矩阵如下所示：
$[K]\{\phi\}=\{b\}$

$K=\begin{bmatrix}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)}\end{bmatrix}$
$\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\}$

[a,b,c]=createMatrixElements(fems);

基于上面给出的矩阵，我们下一步是计算其中的元素，计算公式如下（金-式3.28~3.30）：
$\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}$
由此写出的代码如下（其中a为对角上的元素值，c为次对角元素的值，b就对应上面的非齐次项）：

%% data里面是元胞数组，每个元胞里面都有 alpha，beta，f和l
function [a,b,c]=createMatrixElements(data)sizes=size(data); num=sizes(1);%有限元单元数eleNum=num+1;%有限元单元数+1为节点数%% 生成K矩阵对角元,非对角元和非齐次项的过程a=zeros(eleNum,1);b=zeros(eleNum,1);c=zeros(num,1);a(1)=data(1).alpha/data(1).l+data(1).beta*data(1).l/3;%对角元第一个的公式3.28a(eleNum)=data(num).alpha/data(num).l+data(num).beta*data(num).l/3;%对角元最后一个的公式c(1)=-data(1).alpha/data(1).l+data(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 j=2:eleNum-1a(j)=data(j-1).alpha/data(j-1).l+data(j-1).beta*data(j-1).l/3+data(j).alpha/data(j).l+data(j).beta*data(j).l/3;%中间的公式3.38b(j)=data(j-1).f*data(j-1).l/2+data(j).f*data(j).l/2;%非齐次项剩余项3.41c(j)=-data(j).alpha/data(j).l+data(j).beta*data(j).l/6;%非对角元的公式3.29end
end

2.2.7 强加边界条件

加入边界条件会对上面的k矩阵和b矩阵进行修正，主要在下面代码进行：

[K,b]=createMatrixK([boundary1,boundary2],a,b,c);%强加边界条件

在此函数中，先把之前的abc数组组合成实际的k矩阵：

sizes=size(a);
num=sizes(1);
K=zeros(num);
for j=1:numK(j,j)=a(j);
end
for j=1:num-1K(j,j+1)=c(j);K(j+1,j)=c(j);
end

对于狄利克雷边界条件，实际的代码会根据下式进行修正（金-式3.62）：
$\left.\left[\begin{array}{cccc}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}\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}$

%% 狄利克雷边界条件
if boundaryCondition1.type=='d'K(1,1)=1;b(1)=boundaryCondition1.p;for j=2:numK(1,j)=0;b(j)=b(j)-b(1)*K(j,1);K(j,1)=0;end
end

对于纽曼边界条件，原方程组使用下式进行修正（只在边界上对应的元素进行修正）（金-式3.58~3.59）：
$K_{NN}=\frac{\alpha^{(M)}}{l^{(M)}}+\beta^{(M)}\frac{l^{(M)}}{3}+\gamma$
$b_{N}=f^{(M)}\frac{l^{(M)}}{2}+q$

%% 纽曼边界条件
if boundaryCondition2.type=='n'K(num,num)=K(num,num)+boundaryCondition2.gamma;b(num)=b(num)+boundaryCondition2.q;
end

2.2.8 方程组求解

还在等什么，线性方程组不会求解嘛，此处用高斯消元法求解：

results=solveWithGuassianMethod(K,b);%高斯消元法求解

function results=solveWithGuassianMethod(K,b)%生成a和c向量sizes=size(K);num=sizes(1);a=zeros(num,1);c=zeros(num-1,1);for j=1:numa(j)=K(j,j);endfor j=1:num-1c(j)=K(j,j+1);endfor j=2:numa(j)=a(j)-c(j-1)^2/a(j-1);b(j)=b(j)-b(j-1)*c(j-1)/a(j-1);endresults=zeros(num,1);results(num)=b(num)/a(num);for j=num-1:-1:1results(j)=(b(j)-c(j)*results(j+1))/a(j);end
end

2.2.9 求解反射系数

反射系数就是下面的R（金-式3.94）：
$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 x>L$
我们求解得到results的是在各个节点的电场强度，将位于介质与空气交界处的电场强度带入上式中计算，即可得到反射系数，公式如下：

R=(results(element_num+1)-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;
%基本物理参数
lambda=1;%定义波长，随意设置对最终结果没有影响
k0=2*pi/lambda;%自由空间波数计算公式
L=5*lambda;%基板的厚度为5倍波长-这是给定的条件
element_num=100;%有限元剖分-100个单元
x=linspace(0,L,element_num+1);%单元节点坐标x1到xm
theta=linspace(0,pi/2,20);%求解对应角度下的反射数据
% theta=[0,pi/4];%求解45度下的反射数据
E0=1;%E0为入射场大小，对反射计算结果无影响
%材料参数
epsilon_r=4+(2-0.1*1j)*(1-x/L).^2;%介质板内材料的非均匀介电常数
epsilon_r(element_num+1)=1;%介质板外自由空间介电常数
mu_r=ones(1,length(epsilon_r))*(2-0.1j);%介质板内材料的磁导率
mu_r(element_num+1)=1;%介质板外自由空间磁导率%存放fem的节点信息(结构体信息)
fems=[];
%存放计算得到的反射系数信息
plotData=[];for i=1:length(theta)%对每个theta进行一次有限元求解fems=[];%清空fems节点数组% fem节点数组赋值，并和fems合并为全局节点数组for j=1: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(j+1));'];eval(str);str2=['fems=[fems;fem',num2str(j),'];'];eval(str2);endboundary1=inputBoundaryData('d','min',0);boundary2=inputBoundaryData('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);%强加边界条件results=solveWithGuassianMethod(K,b);%高斯消元法求解R=(results(element_num+1)-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)))