简介一套MATLAB仿真资源面向光纤通信与光子学方向的初学者及仿真爱好者聚焦LP11模式在光纤中的电场分布与光斑形态分析可用于理解多模光纤中模式传播的基本规律。压缩包共5个文件以1个可运行的MATLAB脚本FIberiaLP11.m为核心配合4张仿真结果图像整体仅73KB便于快速下载与离线研读。脚本从光纤折射率参数出发通过数值计算求解电场与磁场分布并用三维曲面等方式呈现LP11模式的双模特征与不对称光斑图像文件则直观展示了电场强度及模式形态适合与文献对照学习。目前已有1163人学习下载。借助这套资源读者既能掌握光纤模式仿真的整体流程也可修改参数观察不同条件下光斑变化为后续光纤设计或信号优化提供可复现的参考样例。1. 从一帧双瓣光斑认识 LP11 模式仿真初用 MATLAB 仿光纤模式的人多半在 LP01 上顺利在 LP11 上卡壳费了半天劲画出来的光斑却是圆斑或者只有边缘亮圈。真正 LP11 的端面强度是两瓣中间有一条暗线两个瓣的相位相差 π电场的角向分布不再像 LP01 那样均匀而是带 cosφ 的调制。这个标题正是围绕这套仿真流程展开的从标量亥姆霍兹方程出发把 LP11 模式的本征传播常数算出来再还原光纤截面的电场分布和光斑。它面向光通信、光纤传感和光学教学的工程师也适合刚接触 COMSOL 又想快速核对结果的仿真用户。下文给模型、给代码也给三个不需要实验设备就能自查的验证手法。2. 用标量亥姆霍兹方程把 LP11 模式写进代码2.1 弱导近似下 LP11 为什么是两瓣分布光纤模式的基础方程是矢量波动方程但当纤芯与包层折射率差很小Δn/n 1%时纵向电场远小于横向电场模式可以用标量 ψ 近似描述。此时横向场满足标量亥姆霍兹方程 ∇²ψ k²n²ψ β²ψ在圆柱坐标系下分离变量令 ψ(r,φ,z) F(r)Φ(φ)exp(-jβz)。角向方程的解是 Φ cos(lφ) 或 sin(lφ)整数 l 对应方位角阶数l 0 是轴对称l 1 出现两个对称瓣这就是 LP11 的来历。两个正交解分别对应 cosφ 和 sinφ 取向实际光纤中两个简并模同时被激励时光斑的朝向取决于两者的相位关系。需要提醒的是LP11 并不是一个孤立的矢量模而是 HE21、TE01、TM01 在弱导极限下简并叠加的结果。若要做高精度近场偏振测量矢量效应会在偏振成像里显现出来但只计算光斑强度和传播常数标量模型误差在千分之一量级足够用于绝大多数工程预研。这也是后面所有 MATLAB 代码都建立在标量近似上的理由。2.2 径向方程和 U、W、V 三个无量纲数把 ψ 代入标量亥姆霍兹方程后径向方程写成d²F/dr² (1/r)dF/dr (k²n² - β² - l²/r²)F 0在纤芯 r a 和包层 r a 两段分别定义三个无量纲数U² a²(k²n1² - β²) W² a²(β² - k²n2²) V² U² W²其中 V 是归一化频率V ak√(n1² - n2²)。U 代表纤芯内横向振荡的快慢W 代表包层中衰减的速率。V 由光纤参数直接决定给定波长、纤芯半径和数值孔径就能算。V 的大小决定一根光纤支持多少模式V 2.4048 时只有 LP01 能导波2.4048 V 3.8317 时 LP11 与 LP01 共存这正是少模光纤常用的工作窗口。纤芯内 F 的解为贝塞尔函数 J_l(Ur/a)包层内要求衰减取修正贝塞尔函数 K_l(Wr/a)。在 r a 处令 F 与 dF/dr 连续得到特征方程UJ_{l1}(U)/J_l(U) WK_{l1}(W)/K_l(W)l 固定时特征方程的根从小到大排列第 m 个根对应 LP_lm 模l 1 的第 1 个根就是 LP11。这是超越方程没有闭式解后面的数值求解就是这个标题里 MATLAB 仿真的核心。2.3 截止条件表查表决定当前光纤能算哪些模式模式方位角阶数 l径向阶数 m截止 Vc说明LP01010基模无截止LP11112.4048J0 的第一个零点LP21213.8317J1 的第一个零点LP02023.8317与 LP21 同截止LP12125.5201J0 的第二个零点这张表的用法很直接先算出当前工作波长下的 V如果 V 小于模式的截止 Vc特征方程在 (0, V) 内没有该模式的根强行求解会得到虚的 W物理上没有意义。对 LP11只有当 V 2.4048 时代码才应进入求解分支。这是排错时第一件要查的事很多“仿真发散”的现场根源就是传入了不满足截止条件的参数。选型上我用标量模块而不是全矢量求解标量模型对 V 10、折射率差小于 1% 的常规石英光纤足够准实现只依赖 besselj 和 besselk 两个内置函数不需要剖网格缺点是分不开 TE01、TM01 和 HE21 的微小拍长差异。后续如果要仿真偏振串音需要升级到全矢量公式但坐标网格和可视化部分的代码可以原样保留。3. 在 MATLAB 中求 LP11 特征方程网格扫描加 fzero3.1 为什么不直接对整个区间调用 fzero特征方程里含两个贝塞尔函数的比值J_l(U) 有零点U 在这些点上残差会穿过无穷大直接在整个 (0, V) 区间上调用 fzero 很容易报错或收敛到奇点而不是物理根。我一般先用粗网格扫描残差的符号变化把根隔离到小区间里再在每个小区间上调用 fzero 精化。这样既拿到根的近似位置又避免奇点干扰。在动手前还应确认目标模式的 U 取值范围。对 LP11U 的理论区间是 (0, 3.8317)其中 3.8317 是 J1 的第一个零点对应远离截止的极限而 V 必须大于 2.4048 才有解。如果算法解出的 U 超过 3.8317说明抓到了高阶根或者错误根应该直接判为异常。3.2 完整函数find_lp.m下面这段代码保存为 find_lp.m输入 l、m 和 V返回第 m 个根的 U、W。function [U, W] find_lp(l, m, V) % FIND_LP 计算阶跃光纤 LP_lm 模的 U、W % 特征方程: U*J_{l1}(U)/J_l(U) W*K_{l1}(W)/K_l(W) % l 方位角阶数, m 径向根序号, V 归一化频率 n 5000; u linspace(1e-8, V*sqrt(1-1e-8), n); res nan(size(u)); for i 1:n w sqrt(V^2 - u(i)^2); jl besselj(l, u(i)); kl besselk(l, w); if abs(jl) 1e-10 || abs(kl) 1e-10 continue; % 跳过奇点避免假符号变化 end res(i) u(i)*besselj(l1, u(i))/jl ... - w*besselk(l1, w)/kl; end % 把连续无 NaN 的区间切出来逐段检测符号变化 valid find(isfinite(res)); breaks [1; find(diff(valid) 1) 1; numel(valid) 1]; roots_ []; for s 1:numel(breaks) - 1 idx valid(breaks(s):breaks(s1)-1); if numel(idx) 2 continue; end sg sign(res(idx)); for k find(diff(sg) ~ 0) roots_(end1) fzero((uu) ffe(l, V, uu), ... [u(idx(k)), u(idx(k1))]); %#okAGROW end end roots_ sort(roots_); if m numel(roots_) error(V %.4f 时不存在 LP%d%d, V, l, m); end U roots_(m); W sqrt(V^2 - U^2); end function e ffe(l, V, u) % 特征方程残差供 fzero 调用 w sqrt(V^2 - u^2); e u * besselj(l1, u) / besselj(l, u) ... - w * besselk(l1, w) / besselk(l, w); end逻辑说明先在整个 (0, V) 上均匀采 5000 个点逐点计算残差贝塞尔函数接近零时残差突变为无穷大这里直接记为 NaN后面按连续段切分避免把奇点误判成根。diff(sign(res))找出相邻点符号变化的位置每个变化区间对应一个根。fzero 的区间端点是网格上相邻的两个点两端残差异号所以能可靠收敛。最后对所有根排序取第 m 个保证传 (1,1) 得到的是 LP11 而不是 LP12。3.3 几个需要按场景调整的参数网格点数 n 默认 5000V 10 时足够V 到 20 以上时建议调到 12000 到 20000否则靠近截止的根容易被漏掉。u 的上限取V*sqrt(1-1e-8)而不是 V是因为 w sqrt(V² - u²) 在 u 接近 V 时趋于 0K_l(w) 数值急剧增大残差会溢出。fzero 默认容差对大多数可视化足够若要高精度传播常数可以加optimset(TolX,1e-12,TolFun,1e-14)。在调用层还需要自己写好截止判断V 2.4048 时直接提示当前光纤不支持 LP11而不是等 find_lp 报错。举例来说C 波段 1550nm、纤芯半径 5μm、n1 1.46、n2 1.45 的光纤V 大约在 3.4 附近LP11 刚好落在可导波区间是很有代表性的测试参数。3.4 怎么判断解出来的 U 是否合理对 LP11U 必然落在 (0, 3.8317) 区间内并且随 V 增大单调变大刚过截止时 U 接近 0模式严重泄漏到包层远离截止时 U 逼近 3.8317场被压进纤芯。如果打印出的 U 不在这个区间优先怀疑网格扫描点数太少、漏掉了根或者 V 计算有误。另一个常见问题是把 λ 的单位写错1550nm 直接写成 1550V 会大三个数量级特征方程的根分布完全乱掉。4. 把 LP11 电场和光斑画出来4.1 从 U、W 构造完整电场分布得到 U、W 后电场分布就是两个区域的拼接。下面这段脚本可以直接运行输出 LP11 的强度图和相位图。% lp11_demo.m —— LP11 电场强度与相位绘制 lambda 1550e-9; % 波长 1550nm a 5e-6; % 纤芯半径 5um n1 1.46; n2 1.45; % 纤芯/包层折射率 k0 2*pi/lambda; V a*k0*sqrt(n1^2 - n2^2); % 归一化频率 [U, W] find_lp(1, 1, V); N 401; x linspace(-3*a, 3*a, N); [X, Y] meshgrid(x, x); R hypot(X, Y); TH atan2(Y, X); F zeros(N); core R a; F(core) besselj(1, U*R(core)/a); F(~core) besselj(1, U) .* besselk(1, W*R(~core)/a) ./ besselk(1, W); F F / max(abs(F(:))); E F .* cos(TH); % l1 的角向项决定两瓣 I abs(E).^2; figure(Color,w); subplot(1,2,1); pcolor(x/a, x/a, I); shading interp; axis image; colorbar; title(LP11 光斑强度); subplot(1,2,2); pcolor(x/a, x/a, angle(E)); shading interp; axis image; colorbar; title(LP11 相位);参数说明F(core)和F(~core)分别对应纤芯内与包层内的径向函数系数besselj(1,U)/besselk(1,W)保证 r a 处两侧连续。cos(TH)是 l 1 的角向调制缺了它画出来就是圆形强度分布那不是 LP11。归一化只维持数值稳定不影响模式形状。窗口取 ±3a 是为了能看到包层里的衰减尾巴只取 ±1.5a 会截断模式。4.2 为什么用 pcolor 而不是 imagescimagesc 按像素中心采样对中心附近的快速变化不够平滑pcolor 配合shading interp做的是网格间线性插值更接近连续场。N 401 时内存占用约 1.3MB完全可接受。如果用 N 101中心暗线和两瓣边界会出现明显锯齿这是初学时最常见的“仿真不收敛”假象之一。实际项目里我会先用 N 256 快速预览定稿时再提到 512。4.3 LP11 仿真的关键参数速查表参数对 LP11 仿真的影响建议取值波长 λ改变 VV 增大时模式更束缚按实际光源C 波段 1550nm 常用纤芯半径 a与 V 成正比少模光纤 4~8μmn1, n2决定 NA 与 V保持 Δn/n 1%网格窗口影响包层尾巴是否被截断3a ~ 6a网格数 N决定中心暗线清晰度256 ~ 512扫描不同 V 看模式演化是判断代码是否正确的有效手段for V [2.5, 3.46, 5, 7] [U, W] find_lp(1, 1, V); fprintf(V %.2f, U %.4f, W %.4f\n, V, U, W); end运行后 U 应随 V 增大而增大W 也随之增大对应光斑从接近截止时的弥散逐渐收缩进纤芯。如果 U 出现回退或者报找不到根回去查截止条件表和 V 的输入单位。4.4 绘制阶段最容易踩的三个错第一把 K_l 当成振荡函数画出包层里的“条纹”。修正贝塞尔函数在实轴上单调衰减不会振荡出现条纹说明 W 带虚部或公式抄错。第二忘了乘角向项 cos(φ)画出来是圆斑。第三用abs(E).^2画强度时取的是复场模平方没有问题但若误写成abs(F).^2且 F 不含角向项同样会丢失两瓣结构。排错时先单独画出cos(TH)的图确认角向项对再叠加径向场。5. 三个快速验证手法和一张可发布的光斑图5.1 径向剖面检查边界连续性沿任意一条直径取径向剖面F 在 r a 处必须连续导数也应连续。导数跳变明显说明 U、W 不是同一特征方程的根常见原因是网格扫描漏根或 fzero 收敛到了奇点。验证代码很短rline linspace(0, 3*a, 2000); Fl zeros(size(rline)); coreL rline a; Fl(coreL) besselj(1, U*rline(coreL)/a); Fl(~coreL) besselj(1,U) .* besselk(1,W*rline(~coreL)/a) ./ besselk(1,W); plot(rline/a, Fl, LineWidth, 1.5);5.2 检查中心暗线两侧的 π 相位跳变LP11 的两个亮瓣相位差应为 π。取强度图两个峰值位置的复场相位做差后回绕到 [-π, π]绝对值应接近 π。这一步能同时确认角向阶数 l 1 没有写错也排除了“画出来像 LP11、实际是 LP01 加了噪声”的情况。5.3 用功率占比检查网格是否够密对强度图做数值积分计算纤芯内功率占总功率的比例。把 N 从 256 提高到 512 后这个比例的变化应小于 1%。变化过大说明网格不够密或者包层窗口取得太小截断了模式尾巴。对接近截止的 V模式泄漏严重窗口要放到 6a 以上。5.4 输出高分辨率图的小技巧出图不要用saveas用exportgraphics(gcf, lp11_spot.png, Resolution, 300)导出的是矢量级清晰度标题里的光斑图直接可以放进论文或报告。色标用parula或turbo不要用jetjet 的彩色条会在暗线附近制造虚假的对比度影响对两瓣结构的判断。本文还有配套的精品资源点击获取
