利用有限元法（FEM）模拟电磁场与样品的相互作用

一、引言

        电磁场与物质的相互作用是理解光学现象的基础。在实际应用中，激光光束与样品的相互作用通常涉及复杂的电磁场分布，尤其在微纳尺度结构中。因此，使用数值模拟方法如有限元法（FEM）来模拟电磁场的分布和传播，对于优化实验条件和解释实验现象具有重要意义。

二、原理分析

        有限元法(FEM)是一种求解偏微分方程（PDE）的数值方法。FEM通过将连续问题离散化为有限的离散元（通常是三角形或四面体单元）来近似求解。对于电磁场模拟，FEM可以用于求解麦克斯韦方程组，从而获得电磁场在复杂几何结构中的分布。

        FEM的基本步骤包括：

1. 几何建模：构建感兴趣区域的几何模型，并对该区域进行离散化。
2. 材料特性定义：设置模型中各部分的材料参数，如电导率、介电常数和磁导率。
3. 边界条件设定：定义模型的边界条件，如电场、磁场或电流源。
4. 求解方程：利用FEM求解器来求解离散化后的麦克斯韦方程，得到电场或磁场的分布。

三、MATLAB代码实现

        我们通过设置一个初始的电场分布（激励源），然后观察这个电磁波如何在空间内传播。用MATLAB代码实现“在一个盒子里制造电磁波并观察它们如何传播”的过程。

1.创建模型：制造一个“空盒子”，也就是一个用来装载我们的电磁场模型的容器

model = createpde();

2.定义几何结构：盒子的形状和里面的障碍物

R1 = [3, 4, 0, 1, 1, 0, 0, 0, 1, 1]'; % 矩形区域

C1 = [1, 0.5, 0.5, 0.2, 0, 0, 0, 0, 0, 0]'; % 圆形区域

gd = [R1, C1];

sf = 'R1+C1';

ns = char('R1','C1');

ns = ns';

g = decsg(gd, sf, ns);

geometryFromEdges(model, g);

• R1和C1：这些代表盒子的形状。我们这里用一个矩形和一个圆形来表示盒子的内部结构，比如矩形可能是盒子的墙壁，而圆形可以看作盒子里面的一个障碍物。
• 组合几何形状：sf = R1+C1表示我们把这个矩形和圆形放在一起，作为我们盒子的最终形状。盒子里有一个矩形区域和一个圆形区域。

3.创建网格：将盒子划分成小块，称为“网格”。就像把一个面包切成小块一样，这就是FEM的第一步：网格划分。每个小单元将用于近似原始区域内的物理场（如电场）。

generateMesh(model, 'Hmax', 0.05);

4.定义材料属性和PDE系数：盒子里面的材料特性。

epsilon = 1; % 介电常数

mu = 1; % 磁导率

c = 1/(epsilon * mu); % 波速

specifyCoefficients(model, 'm', 0, 'd', 1, 'c', c, 'a', 0, 'f', 0);

• 材料属性：我们定义了盒子内部的材料特性。epsilon 和 mu 分别表示材料的电和磁特性（比如空气、玻璃等都有不同的电磁特性）。c 是波速，由这些材料特性决定。
• PDE系数：这些系数用来描述盒子里面的电磁场如何传播。这些系数相当于在告诉MATLAB，“波在这个盒子里应该怎样运动”。

5.设置边界条件：盒子的边界是什么样的

applyBoundaryCondition(model, 'Edge', 1:model.Geometry.NumEdges, 'g', 0, 'q', 0);

• 边界条件：这些边界条件允许电磁波自由地离开盒子，类似于在开放的空间中，波可以自由传播，不会被墙壁反射回来。

6.定义初始条件：在盒子中心放一个小激励源

u0 = @(location) exp(-100*((location.x-0.5).^2 + (location.y-0.5).^2));

setInitialConditions(model, u0);

• 初始条件：我们在盒子中央放了一个“小激励源”，这就像是在盒子中央敲了一下，产生了一个初始的电磁波。这个初始激励源有点像一个脉冲或爆炸，它会在时间 t = 0 时开始传播。

7.求解PDE：看电磁波是如何传播的

tlist = linspace(0, 0.1, 100); % 时间步

results = solvepde(model, tlist);

• 求解过程：MATLAB开始计算盒子里电磁波的传播过程。我们让它从时间 t = 0 开始，一直到 t = 0.1 秒，通过100个时间步来观察这个传播过程。就像你用一个慢动作摄像机记录一个短暂的爆炸过程一样。

8.可视化结果：展示电磁波的传播动画

u = results.NodalSolution;

figure;

for i = 1:length(tlist)

    pdeplot(model, 'XYData', u(:, i), 'ZData', u(:, i), 'Mesh', 'off');

    title(['Time = ', num2str(tlist(i))]);

    zlim([-0.5 0.5]); % 设置Z轴范围以更好地展示电场分布

    drawnow;

end

• 可视化：我们通过动画来展示电磁波在盒子里的传播过程。每一帧显示的是某个时间点上的电磁波形状。你会看到电磁波从中心扩散到整个盒子，反映了电磁场如何在时间和空间中演化。

完整代码：

% 创建模型

model = createpde();


% 创建矩形区域

R1 = [3, 4, 0, 1, 1, 0, 0, 0, 1, 1]'; % 矩形的几何描述


% 创建圆形区域

C1 = [1, 0.5, 0.5, 0.2, 0, 0, 0, 0, 0, 0]'; % 圆形的几何描述, 填充至与R1相同行数


% 将几何对象合并

gd = [R1, C1];


% 定义集合函数和名称矩阵

sf = 'R1+C1';

ns = char('R1','C1');

ns = ns';


% 使用 decsg 创建几何结构

g = decsg(gd, sf, ns);

geometryFromEdges(model, g);


% 创建网格

generateMesh(model, 'Hmax', 0.05);


% 定义材料属性和PDE系数

epsilon = 1; % 介电常数

mu = 1; % 磁导率

c = 1/(epsilon * mu); % 介质中的波速

specifyCoefficients(model, 'm', 0, 'd', 1, 'c', c, 'a', 0, 'f', 0);


% 修改边界条件为 Neumann 边界条件，允许波传播

applyBoundaryCondition(model, 'Edge', 1:model.Geometry.NumEdges, 'g', 0, 'q', 0);


% 定义初始条件为高斯脉冲电场分布

u0 = @(location) exp(-100*((location.x-0.5).^2 + (location.y-0.5).^2));

setInitialConditions(model, u0);


% 求解PDE

tlist = linspace(0, 0.1, 100); % 时间步

results = solvepde(model, tlist);


% 提取结果并可视化

u = results.NodalSolution;

figure;

for i = 1:length(tlist)

    pdeplot(model, 'XYData', u(:, i), 'ZData', u(:, i), 'Mesh', 'off');

    title(['Time = ', num2str(tlist(i))]);

    zlim([-0.5 0.5]); % 设置Z轴范围以更好地展示电场分布

    drawnow;

end

实验结果