简介:面向光纤通信、光电子学与数值模拟领域的研发人员和研究生,这份MATLAB程序围绕单模光纤与LED光源之间的耦合效率计算展开,解决光学系统设计中如何快速评估能量传输损失的问题。程序基于FDE有限差分方程建立光传播数值模型,可调节光源发射角、光纤芯径、数值孔径、透镜焦距等关键参数,完整模拟光场进入光纤后的分布变化,并给出耦合效率结果,便于对不同设计方案进行敏感性分析与优化,在光纤通信、激光技术、光学传感等场景中均有参考价值。压缩包体积仅668KB,共包含2个文件,其中m脚本为核心计算代码,PDF文档为算法原理和操作流程说明,代码与文档结合能够让读者既掌握FDE建模思路,又可直接修改参数复现仿真结果。已有1529人浏览学习,适合需要快速搭建光纤耦合仿真模型、验证理论公式或继续扩展为更完整光学设计工具的研究者。
1. LED 光耦合进单模光纤,难的不是对准而是模场匹配
把一块发光面积几百微米的 LED 对准一根芯径 9 μm 的单模光纤,很多人第一反应是“只要对准就行”。实际测量下来,端面对准时耦合效率通常只有几个百分点,再往上推,问题都不在横向对准,而在模场匹配。LED 按朗伯体发光,角度发散接近 ±90°,而单模光纤的接收角被数值孔径限制在十几度量级,大部分功率在到达纤芯前就漏掉了。更麻烦的是,这类失配没法靠简单对准消除。用 MATLAB 仿真的好处是能分别建模,把光源尺寸、透镜倍率、NA 和光纤基模拆开算;压缩包里的 p183_exam4_8.m 正是沿着这条链路做单模光纤耦合效率计算。下面按实际调程序的顺序,把 FDE 建模、重叠积分、参数调节和批量扫描逐层拆开讲。
2. FDE 建模与模场重叠积分:先弄清楚耦合效率在算什么
2.1 为什么要用离散方程而不是解析公式
单模光纤的基模在弱导近似下可以用高斯函数近似,直接套公式也能得到粗略的重叠积分结果。但真实光纤有折射率差、截止波长和模场边界效应,解析式在芯径偏小、数值孔径偏大时误差会明显增大。FDE 的做法是把光纤横截面切进网格,把亥姆霍兹方程里的偏导数替换成差分,原问题就变成一个稀疏矩阵特征值问题,解出来的特征向量就是模式场分布。
这种方式在规则光纤截面上比有限元更轻,网格分辨率可以直接控制,不需要复杂的前处理。p183_exam4_8.m 这类教材程序选用 FDE 而不是 COMSOL 类工具,就是为了让参数透明、步骤可查。在光纤光学语境里,FDE 可以理解为基于有限差分方程的数值建模,商用软件中类似算法也被称为 Finite-Difference Eigenmode,本质是一回事。
2.2 网格生成与折射率分布
常见做法是先定义计算窗口和网格数,再按坐标填充折射率分布。下面这段代码建立一个边长为 20 μm 的方形窗口,把半径为 4.2 μm 的阶跃折射率光纤放进去。
% 网格参数:窗口大小和网格数共同决定空间分辨率 L = 20; % 计算窗口边长,单位 μm N = 200; % 每个方向上的采样点数 x = linspace(-L/2, L/2, N); dx = x(2) - x(1); [X, Y] = meshgrid(x, x); % 光纤参数,先按阶跃折射率设置 n_core = 1.4681; % 纤芯折射率 n_clad = 1.4628; % 包层折射率 a = 4.2; % 纤芯半径,单位 μm % 折射率分布矩阵,用平方量参与计算 epsr = n_clad^2 * ones(N, N); epsr(X.^2 + Y.^2 <= a^2) = n_core^2;这里L和N的比值就是空间分辨率dx。单模光纤基模的特征尺寸通常在 4~5 μm,dx至少要小于特征尺寸的五分之一,否则模场边界会被离散误差抹圆。epsr存放的是介电常数而不仅是折射率,因为波动方程最终要写成关于介电常数的矩阵形式;如果只在后面对填折射率值,边界误差会在模式求解时被放大。
网格参数与精度的对应关系见表:
| 参数 | 作用 | 推荐范围 |
|---|---|---|
L | 计算窗口边长 | 芯径的 3~5 倍,太小会截断模场 |
N | 单方向采样数 | 100~300,越大越接近连续解 |
dx | 空间步长 | 小于模场尺寸的 1/5 |
| 边界条件 | 窗口边缘处理 | 零边界简单但需配合大窗口,或用 PML |
2.3 拉普拉斯算子与基模求解
模式求解的核心是纵向传播常数对应的特征值问题。在二维横截面上,用稀疏矩阵构造拉普拉斯算子,再把介电常数矩阵加到对角线上,就能用eigs提取基模。
lambda0 = 1.55; % 工作波长,单位 μm k0 = 2 * pi / lambda0; % 一维二阶差分算子,三点中心差分格式 e = ones(N, 1); D2 = spdiags([e, -2*e, e], -1:1, N, N) / dx^2; % 二维拉普拉斯算子:x 方向与 y 方向叠加 Lap = kron(speye(N), D2) + kron(D2, speye(N)); % 波动方程离散矩阵:A = 横向拉普拉斯 + k0^2 * epsr % 特征值最大的一项对应传播常数最大的基模 A = Lap + k0^2 * spdiags(epsr(:), 0, N*N, N*N); [Evec, ~] = eigs(A, 1, 'largestabs'); E_fiber = reshape(Evec, N, N);spdiags在这里把epsr矩阵转成 N²×N² 的对角稀疏矩阵,内存占用远小于diag展开。kron组合出二维算子后,矩阵规模虽然很大,但绝大多数元素为零,eigs只用其中前几个特征向量,速度很快。零边界条件会人为截断泄漏场,所以窗口不能贴得太紧;如果后续要算更具泄漏性的结构,需要在边界加 PML 层,这一项在教材程序里通常预留了接口。
2.4 耦合效率的物理定义:重叠积分
解出光纤端面模场后,耦合效率由入射场与光纤模场的重叠积分决定。归一化后的表达式为:
% E_fiber 与 E_source 均为 N x N 复场分布 nume = abs(sum(sum(conj(E_fiber) .* E_source)))^2; deno = sum(sum(abs(E_fiber).^2)) * sum(sum(abs(E_source).^2)); eta = nume / deno;分子是两场的相关叠加,分母把两场各自能量归一,eta介于 0 到 1 之间。这里必须对E_fiber取共轭,因为光纤模式是沿 z 方向传播的复场,横向分布可能带相位;如果省略共轭,遇到波前倾斜或相位失配时效率会被明显高估。这个积分公式是整个程序后续所有参数扫描的度量标准。
3. 拆解 p183_exam4_8.m:从光源参数到效率输出
3.1 程序结构如何对应实验链路
从命名来看,p183_exam4_8.m是教材示例程序,编号里的 p183 通常指页码,exam4_8 指第 4 章第 8 个例子。这类程序的特点是参数集中在文件头部,链路按“光源 → 耦合系统 → 光纤”的顺序展开。压缩包里的“光纤耦合样稿2.pdf”一般对应推导说明或需求文档,用来反查公式和参数来源。
程序真正做三件事:确定 LED 光源在光纤端面处的等效场,求解单模光纤基模,计算两场重叠积分。光源部分最常见做法是按朗伯体建模,再用一个带束腰半径的高斯场近似端面处的入射光斑。
3.2 参数区里放什么
下面表格列出这类程序中几乎必然出现的物理量,命名习惯大同小异:
| 变量 | 物理含义 | 典型取值 |
|---|---|---|
lambda0 | 工作波长 | 1.31 或 1.55 μm |
a_core | 纤芯半径 | 4~5 μm |
n_core/n_clad | 纤芯与包层折射率 | 1.468 / 1.463 |
NA | 数值孔径 | 0.1~0.14 |
w_source | LED 经透镜后在端面的等效束腰 | 2~8 μm |
z_offset | 光源到光纤端面距离 | 0~几十 μm |
这些参数不是随便填的。lambda0决定归一化频率 V 值,NA决定光纤能接收的角度范围,w_source则直接控制重叠积分里的模场匹配程度。程序中常见问题是把芯径写成直径而不是半径,a_core一旦翻倍,V 值立刻翻倍,单模条件可能被破坏。
3.3 核心计算流程
核心流程可以归纳成下面这段骨架代码,实际文件里的函数封装方式可能不同,但逻辑顺序一致。
% Step1 参数区 lambda0 = 1.55; a_core = 4.2; NA = sqrt(n_core^2 - n_clad^2); % Step2 生成网格与折射率分布 % 同上文 2.2 节内容 % Step3 用 eigs 求解光纤基模 E_fiber = solve_fiber_mode(epsr); % Step4 构造 LED 光场 % 简化:用束腰为 w_source 的高斯场近似 E_source = exp(-(X.^2 + Y.^2) / w_source^2); % Step5 重叠积分与效率输出 eta = abs(sum(sum(conj(E_fiber) .* E_source)))^2 ... / (sum(sum(abs(E_fiber).^2)) * sum(sum(abs(E_source).^2))); fprintf('LED -> SMF 耦合效率 = %.2f%%\n', eta * 100);E_source的构造是整个程序的软肋。把 LED 当成一束高斯光,等价于把它看作空间相干光源,这并不严格等于真实 LED 的非相干辐射。严格处理需要把发光面分解成多个互不相干的子源,对每个子源做功率求和,再取平均效率。在初版程序里,先以高斯近似跑通链路,后续替换成子源阵列,是一种工程上合理的迭代路径。
3.4 与样稿 PDF 配合反查参数
拿到“光纤耦合样稿2.pdf”时,先看两处:一是坐标系定义,PDF 里的 z 轴正方向必须与 MATLAB 中光束传播方向一致;二是孔径公式,确认程序里用的是数值孔径 NA 还是半发散角。样稿里一旦出现过角度的半角全角混用,程序结果会整体偏移一个倍数关系。
检查 V 数是一个快速校验手段。单模条件要求 V < 2.4048,按下式计算:
V = 2 * pi / lambda0 * a_core * NA; assert(V < 2.4048, '当前参数不满足单模条件');程序若在单模条件之外运行,eigs提取出的“基模”可能混入高阶模,重叠积分算出的效率就会失真。样稿 PDF 如果没有明确给出截止波长,建议先跑这个断言。
4. 耦合效率上不去的典型原因与参数调整
4.1 数值孔径、芯径与光源尺寸的匹配关系
耦合效率的高低取决于光纤模场与入射光场在横向和角度两个维度的重合程度。横向维度看束腰:光纤基模束腰w_f与入射光斑w_s越接近,效率越高。假设两者都是高斯分布且完全对准,最大耦合效率可以写成:
% wf 为光纤模场半径,ws 为光源束腰半径 eta_max = 4 / (wf/ws + ws/wf)^2;当ws = wf时效率为 1;当ws是wf的 3 倍时,效率只剩约 44%。这个公式提示一个方向:对于发散角大、光斑大的 LED,直接端面耦合不可能做到高耦合,必须通过透镜系统在光纤端面处重新整形光斑。程序中如果只调光源功率而忽略ws,效率曲线不会出现实质变化。
4.2 加透镜:用倍率同时换光斑和角度
透镜耦合的核心是横向倍率M。光束经过透镜后,光斑尺寸变为原来的M倍,角度发散角变为原来的1/M倍。所以可以用一个放大倍率把宽角度光源压缩成细光束,代价是光斑变大;再配合焦距选择,让光斑尺寸落到光纤模场附近。先算目标倍率:
M_opt = w_fiber / w_source; % 光斑匹配对应的横向倍率 NA_after = sin(theta_div) / M_opt; % 经过倍率后的等效数值孔径 if NA_after > sqrt(n_core^2 - n_clad^2) warning('角度失配仍偏大,需减小光源发散角或增大倍率'); end这个计算里M_opt只解决了光斑匹配,角度约束是否满足要看NA_after。LED 全发散角可能超过 100°,直接除以 10 倍倍率后仍在 0.2 左右,对单模光纤 0.1 上下的 NA 来说仍然偏大,所以实际工程中光源尺寸本身就有硬约束。调程序时如果发现计算效率低于预期,优先检查这一行而不是去优化透镜数值。
4.3 常见异常与排错路径
模式不收敛、效率超过 1、结果随网格抖动是三类最常见问题。下表给出典型现象与排查方向:
| 现象 | 可能原因 | 处理方式 |
|---|---|---|
eigs不收敛 | 窗口太小,边界把模场截断 | 增大L,或加 PML 边界 |
| 效率大于 1 | 光源场或模式场做了单点归一化 | 检查分母是否用矩阵求和,而不是某个峰值 |
效率随N明显变化 | 网格太粗,模场形状失真 | 增加N,直到效率变化小于 1% |
| 扫描曲线有锯齿 | 光源场未对准,重叠积分相位振荡 | 检查网格原点是否在纤芯中心 |
| 偶发 NaN | epsr含有复数且对角矩阵出错 | 确认epsr展开后是 N²×N² 稀疏对角阵 |
其中效率大于 1 最容易出现在初学者代码里。重叠积分分母必须是两场各自总能量,写成sum(abs(E).^2),而不是max(abs(E(:)).^2)。后者只取了峰值功率,会丢掉能量归一化,效率结果直接翻倍甚至更高。
5. 批量扫描与解析验证:让单次计算变成设计工具
5.1 对光源束腰做参数扫描
单次耦合效率只能说明一个点上的表现。把束腰半径、轴向距离或横向偏移设成数组循环求解,就能得到效率随参数的变化曲线,这是把 MATLAB 程序变成设计工具的第一步。
ws_list = linspace(1, 10, 10); eta_list = zeros(size(ws_list)); for k = 1:numel(ws_list) E_source = exp(-(X.^2 + Y.^2) / ws_list(k)^2); eta_list(k) = abs(sum(sum(conj(E_fiber) .* E_source)))^2 ... / (sum(sum(abs(E_fiber).^2)) * sum(sum(abs(E_source).^2))); end plot(ws_list, eta_list * 100, '-o'); xlabel('光源束腰半径 ws (μm)'); ylabel('耦合效率 (%)');这段循环不需要对模式场重复求解,因为驱动光纤的E_fiber不随光源变化。把eigs放在循环外是提速的关键,否则每个参数点都要重算一次模式,10 个点还能接受,100 个点就明显卡顿。
5.2 用解析高斯式兜底验证
仿真结果需要一条参考线。用高斯-高斯重叠解析公式作为对照,验证程序里模式求解和归一化逻辑是否正确。
wf = 5.2; % 从模式场二阶矩提取的光纤模场半径 eta_analytic = 4 ./ (wf./ws_list + ws_list./wf).^2; hold on; plot(ws_list, eta_analytic * 100, '--'); legend('FDE 数值结果', '高斯解析参考');当数值曲线与解析参考线接近但不完全重合时,说明模式场不是完美高斯,这是正常现象;如果趋势偏移方向相反,则应该回到重叠积分公式确认共轭和归一化。wf应该用二阶矩提取而不是从光纤芯径直接取,更稳妥的做法是用模式场强度分布计算等效半径。
5.3 网格收敛性检查
参数扫描之前,先对网格数做一次收敛性测试:分别用 N = 100、150、200、300 计算同一组参数下的效率,观察结果是否稳定。效率随 N 增大单调逼近一个定值时,说明网格分辨率已经足够;如果 300 点的结果还在明显变化,就要检查窗口边界或模式求解精度。扫描和验证做完后,把通过校验的参数组合保存成.mat文件,再交给下一步的透镜参数优化循环,这套 MATLAB 程序就不只是算单个效率,而是能支撑完整的光路设计迭代。
本文还有配套的精品资源,点击获取