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

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

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

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

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

QAQ
QAQ
QAQ

0、问题定义

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

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

其中介电常数是非均匀的，其表达式为：
$${ϵ}^{\mathrm{r}}=4+(2^{\mathrm{I}}-\mathrm{j}0.1)(1-x/L{)}^2$$
磁导率为固定的复数：
$${\mu }_{\mathrm{r}}=2-\mathrm{j}0.1$$
介质板的厚度为5个波长：
$$L=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({ϵ}_{\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 基础参数定义

1. 定义波长，随意设置对最终结果没有影响
2. 自由空间波数计算公式
3. 基板的厚度为5倍波长-这是给定的条件
4. 有限元剖分-100个单元
5. 单元节点坐标x1到xm
6. 求解对应角度下的反射数据
7. 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 材料信息定义

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

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

%材料参数
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节点信息，其实际上是个结构体，主要包含下面信息：

其中：

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

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

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

2.2.4 节点数组赋值

代码如下：

fems=[];%清空fems节点数组
% fem节点数组赋值，并和fems合并为全局节点数组
for j=1:length(x)-1
    str=['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}\varphi }{\mathrm{d}x}\right)+\beta \varphi =f0<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({ϵ}_{\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};
    else
        error('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}\varphi }{\mathrm{d}x}+\gamma \varphi \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 ] { ϕ } = { b } [K]\{\phi\}=\{b\} [K]{ϕ}={b}

$$K=\left[\begin{array}{cccc}{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{array}\right]$$
{ 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）：
$$\begin{array}{rl}& 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{array}$$
由此写出的代码如下（其中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.28
    a(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.29
    b(1)=data(1).f*data(1).l/2; %非齐次项第一项-f为已知的源或者激励
    b(eleNum)=data(num).f*data(num).l/2; %非齐次项最后一项3.30
    
    for j=2:eleNum-1
        a(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.38
        b(j)=data(j-1).f*data(j-1).l/2+data(j).f*data(j).l/2;%非齐次项剩余项3.41
        c(j)=-data(j).alpha/data(j).l+data(j).beta*data(j).l/6;%非对角元的公式3.29
    end
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:num
    K(j,j)=a(j);
end
for j=1:num-1
    K(j,j+1)=c(j);
    K(j+1,j)=c(j);
end

对于狄利克雷边界条件，实际的代码会根据下式进行修正（金-式3.62）：
$$[\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}]\left\{\begin{array}{c}{\varphi }_1\\ {\varphi }_2\\ {\varphi }_3\\ {\varphi }_4\end{array}\right\}=\left\{\begin{array}{c}p\\ b_2-K_{21}p\\ b_3-K_{31}p\\ b_4-K_{41}p\end{array}\right\}$$

%% 狄利克雷边界条件
if boundaryCondition1.type=='d'
    K(1,1)=1;
    b(1)=boundaryCondition1.p;
    for j=2:num
        K(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:num
        a(j)=K(j,j);
    end
    for j=1:num-1
        c(j)=K(j,j+1);
    end
    for j=2:num
        a(j)=a(j)-c(j-1)^2/a(j-1);
        b(j)=b(j)-b(j-1)*c(j-1)/a(j-1);
    end
    
    results=zeros(num,1);
    results(num)=b(num)/a(num);
    for j=num-1:-1:1
        results(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 }+R{E}_0{\mathrm{e}}^{-\mathrm{j}{k}_{0}x\cos \theta }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)-1
        str=['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

    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))));
    
    [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)))