管状中空芯光纤拉制参数计算与Matlab标定方法

简介：面向光纤通信与特种光纤制备方向的Matlab计算资源，用于求解管状中空芯光纤拉制过程中的关键参数，帮助研究者与相关专业学生建立起从光纤结构到拉制条件的计算路径。代码采用参数化编程，用户可方便调整几何与工艺参数，并附有可直接运行的案例数据，适配Matlab 2014a/2019a/2024a环境。资源包共20个文件，包括11个m源码文件、2个mat数据文件，以及png效果图、mlapp交互界面、pdf说明文档、exe可执行程序和md说明文件，整体约2.17MB。源码按文档、运行脚本及工具目录组织，注释明细，便于在课程设计、期末大作业或毕业设计中二次开发与复用。目前已有40人学习，适合电子信息工程、数学等专业的算法仿真实践。

1. 管状中空芯光纤的拉制参数为什么要算，而不是试

中空芯光纤的拉制最让人头疼的一点是：外径可以用激光测径仪实时监控，壁厚却要等一根光纤拉完、切样、送到显微镜下才能知道结果。而壁厚和外径根本不是同一个缩放规律，预制棒进给慢、收丝快，玻璃在两千度软化区被拉长的同时，表面张力在不断收薄管壁，管内气压又在试图撑开它，三者叠加之后，最终壁厚往往和按速度比推算的值差出一倍以上。靠试拉来逼近目标壁厚，烧掉的不只是预制棒，还有整个工艺调试周期。更可靠的路线是把拉制温度、进给速度、收丝速度和管内气压这些参数先放进 Matlab 里做正演计算，再用试拉样品反演修正模型，最后让参数计算替工程师决定下一炉怎么开。这篇按这个顺序写透：先立物理模型，再给能直接改参数的 Matlab 代码，最后讲怎么用实测数据把模型校准到现场可用。

2. 拉制参数的物理模型：粘度、表面张力与气压如何决定壁厚

2.1 反谐振条件先把壁厚目标画死

管状中空芯光纤的导光机制和实心光纤完全不同，它靠的是反谐振反射：光在管壁这个高折射率薄层里发生干涉，特定波长的光被反射回空气芯中，其他波长的光则泄漏出去。这个选择性反射直接由管壁厚度决定，对工作波长 λ 和管壁折射率 n，一阶反谐振壁厚满足：

t ≈ λ / (4 × √(n² - 1))

拿近红外 1.55 μm 和石英玻璃 n≈1.444 来算，一阶壁厚差不多在 1.3 μm 量级。工程上为了压低损耗或匹配特定波段，经常取更高阶反谐振条件，把壁厚压到几百纳米。这带来一个直接后果：壁厚偏掉 100 nm，谐振谷就会移动上百纳米，光谱特性完全变样。拉制参数计算的第一性目标，就是把壁厚控制在这个严格的窗口内，而不是“差不多就行”。

这个公式也解释了为什么中空芯光纤比普通光纤难拉得多。普通单模光纤的模场直径宽容度很大，芯包结构稍微走样，损耗从 0.2 dB/km 涨到 0.3 dB/km 而已；而管状空芯光纤的壁厚一旦偏离反谐振条件，导光机制直接失效，光纤就从“能导光”变成“完全不能导光”，中间没有缓冲。所以拉制参数的设定不是优化问题，而是约束满足问题，必须先算准壁厚，再谈其他。

2.2 软化区的三股力：粘度、表面张力、管内气压

预制棒进入高温炉后，在一个通常只有几十毫米长的软化区里完成从粗管到细管的全部变形，这个过程同时存在三种相互对抗的物理作用：

第一是轴向拉伸。预制棒以毫米每分钟的量级进给，光纤以米每分钟的量级收丝，直径缩小的同时，管壁被轴向拉长。这个拉长过程不是均匀的，离软化区入口越远，玻璃温度越高、粘度越低，变形也越剧烈。第二是表面张力收缩。高温下石英玻璃的表面张力约在 0.3 N/m 左右，它倾向于把管腔缩小、把壁厚增厚，是壁厚控制的负面因素。壁厚越薄，这个收缩压力 γ/t 越大，所以细壁管的壁厚失控风险是指数级上升的。第三是管内气压膨胀。在预制棒末端接入气压源，用几百帕到几千帕的正压撑住管腔，对抗表面张力收缩，是控制芯径和壁厚的主要调节手段。

三个力的量级关系可以用一个简单的比值判断：当 γ/t 远大于管内气压 P 时，管腔会持续塌缩，最终形成实心棒；当 P 远大于 γ/t 时，管壁会被吹胀变薄，甚至破裂。拉制参数计算的核心，就是在给定的温度场和速度场下，找到让表面张力和气压基本平衡的那个工作点。注意这里有个经常被忽略的细节：管内气压在整个软化区并不是常数，预制棒一端封闭、另一端接气压源时，气体在高温区的膨胀和流动会让有效压力远低于冷端表压，现场标定时不能直接把气压表的读数当作用在软化区上的真实压力。

2.2.1 温度决定粘度，粘度决定一切敏感度

石英玻璃的粘度随温度变化极陡，工程上常用 VFT 形式描述。从 1800℃ 到 2100℃，粘度可能下降超过两个数量级，而表面张力和气压基本不随温度变化。这意味着同一台拉丝塔，炉温设定差 20℃，壁厚可能就是另一个量级的结果。

function eta = silica_viscosity(T) % 石英玻璃粘度-温度关系，VFT形式 % T: 等效温度，单位 degC % eta: 粘度，单位 Pa.s T_K = T + 273.15; % 高温段工程近似参数，覆盖 10^3 ~ 10^8 Pa.s 区间 A = -6.37; B = 2.46e4; T0 = 0; eta = 10^(A + B ./ (T_K - T0)); end

这段代码返回的粘度直接喂给后面的微分方程。A、B、T0 这三个参数在不同文献里的取值略有出入，因为石英玻璃的粘度还受羟基含量、杂质和气氛影响。实际使用中，如果手头有高温粘度计数据，应该用供应商实测值替换掉这组默认参数。拉制工程师嘴里说的“等效温度”，指的就是把炉内复杂的温度分布压缩成一个能给出正确粘度的等效值，这个值通常和热电偶读数差 50℃ 以上，需要反向标定，后面第 4 章会专门处理。

2.3 质量守恒约束下的壁厚-半径耦合方程

要算拉制参数，必须先有一个能在电脑上反复试错的状态演化模型。中空芯光纤的软化区虽然几何复杂，但工程上常用一维准稳态近似：把软化区沿轴向分成很多薄片，每个薄片内玻璃处于局部平衡。第一个方程是质量守恒，石英玻璃在高温下可视为不可压缩流体，任意截面的玻璃面积 A(z) 乘以该处轴向速度 v(z) 必须等于体积流量 Q：

A(z) × v(z) = Q

这个式子很简单但很关键，它把直径和速度死死绑在一起。第二个方程是径向演化方程，描述表面张力和气压如何驱动管径变化。对薄壁管，径向收缩速度满足简化形式：

v × (dR/dz) = -((γ/t - P) / (4η)) × R

符号约定：γ 是表面张力，t 是壁厚，P 是管内表压，η 是粘度，R 是管外半径。表面张力项 γ/t 是单位面积上的收缩应力，内压 P 与之对抗，差值通过粘度阻尼转化为变形速度。这个方程也叫 Lamb 收缩模型在拉制坐标下的改写，严格只适用于粘度均匀的薄壁段，实际炉内温度不均匀时，要人为加一个经验修正系数 C_exp，乘在右边，系数由试拉样品反演得到。

2.3.1 壁厚变化的第二个方程来源

有了 dR/dz 之后，壁厚 t 的变化可以直接从质量守恒推导。由 A = 2πRt，取对数微分得到：

dt/dz = -t × (v'/v) - (t/R) × (dR/dz)

其中 v' = dv/dz 是轴向速度梯度。速度梯度由拉丝速度差除以软化区长度近似得到，软化区长度可以看熔化炉的加热区长，也可以在仿真里当调节参数。这个方程的好处是自动满足体积守恒，不会出现半径和壁厚各自演化、最后对不上的情况。代码实现时注意量纲统一，R 和 t 用毫米、v 用毫米每秒、η 用帕秒、γ 用牛每米，计算前统一转成标准单位。

下面这张表给出了各个参数在模型中的典型量级，方便判断什么参数影响大、什么参数影响小：

2.4 哪些参数要调，哪些参数要测

上表暴露了一个常见误区：很多人上来就想把温度、压力、速度同时设成一个最优组合，但温度的有效值（等效温度）其实是测不准的。炉内温度分布受加热丝状态、气流和保护气氛影响，热电偶读数只能作为参考。真正能调的是压力、收丝速度、进给速度这三项，温度和修正系数必须用试拉样品反演。我一般建议把拟合参数控制在两个以内，优先拟合等效温度 T_eff 和经验修正系数 C_exp，气压和速度直接采用给定值，这样反问题才稳定。

3. 用 Matlab 代码计算拉制温度、速度与气压的落地实现

3.1 主程序骨架与求解流程

把上一章的 ODE 系统写成可运行的 Matlab 代码，求解流程分五步：定义材料与过程参数、构造初始几何、在软化区长度上积分求解、提取出口外径和壁厚、判断是否满足反谐振目标。下面这段主程序体现整体结构，参数都在文件头部集中修改，方便做批量扫描。

% draw_parameter_calc.m % 管状中空芯光纤拉制参数正演计算 % 求解软化区内管外半径 R(z) 与壁厚 t(z) 的演变 clear; clc; % ---- 材料与过程参数 ---- T_eff = 2050; % 等效温度 [degC] P_gauge = 600; % 管内表压 [Pa] gamma = 0.3; % 石英表面张力 [N/m] v_p = 10; % 预制棒进给速度 [mm/min] v_draw = 15; % 收丝速度 [m/min] C_exp = 1.2; % 经验修正系数，由试拉数据标定 % ---- 预制棒几何 ---- R0 = 12.5; % 预制棒外半径 [mm] t0 = 1.5; % 预制棒壁厚 [mm] A0 = pi * (2*R0*t0 - t0^2); % 玻璃截面积 [mm^2] % ---- 速度换算与体积流量 ---- v_in = v_p / 60; % 进给速度 [mm/s] v_out = v_draw * 1000 / 60; % 收丝速度 [mm/s] Q = A0 * v_in; % 体积流量 [mm^3/s] % ---- 求解软化区演化 ---- L_soft = 60; % 软化区等效长度 [mm] zspan = [0 L_soft]; y0 = [R0; t0]; fparams = struct('Q', Q, 'v_in', v_in, 'v_out', v_out, ... 'T_eff', T_eff, 'P_gauge', P_gauge, ... 'gamma', gamma, 'C_exp', C_exp, 'L', L_soft); [z, y] = ode45(@(z, y) draw_ode(z, y, fparams), zspan, y0); % ---- 提取结果 ---- R_out = y(end, 1); % 最终外半径 [mm] t_out = y(end, 2); % 最终壁厚 [mm] fprintf('最终外径: %.3f mm\n', 2*R_out); fprintf('最终壁厚: %.3f um\n', t_out*1000);

代码逻辑上，ode45 负责从软化区入口积分到出口，每一小步都调用 draw_ode 计算半径和壁厚的变化率。fprintf 输出的两个值是整个计算的最终结果，直接和目标壁厚比对。注意 Q 在这里由进给速度和预制棒截面决定，收丝速度不直接出现在体积流量里，而是通过 v_out 影响软化区的速度梯度，进而改变轴向拉伸率。想考察不同收丝速度的效果，改 v_draw 即可，不需要动 Q。

3.2 ODE 方程实现：注意量纲和发散保护

子函数 draw_ode 把上一章的物理方程翻译成 Matlab 语法。这里最容易出错的是单位换算，我在代码里用 mm 保存几何量，但物理方程里的表面张力、气压和粘度都是国际单位，所以进入计算前先转成米和帕。

function dydz = draw_ode(z, y, fparams) % 拉制软化区内管外半径与壁厚的一阶ODE % y(1): 管外半径 R [mm] % y(2): 壁厚 t [mm] R = y(1); t = y(2); t = max(t, 1e-4); % 防止壁厚趋零导致除零，下限 0.1 um % 体积流量守恒得到的局部轴向速度 [mm/s] Q = fparams.Q; A = pi * (2*R*t - t^2); v = Q / A; % 速度场简化：入口到出口线性加速 v_in = fparams.v_in; v_out = fparams.v_out; L = fparams.L; dvdz = (v_out - v_in) / L; % 转换成标准单位 R_m = R * 1e-3; t_m = t * 1e-3; v_m = v * 1e-3; dvdz_m = dvdz * 1e-3; % 每毫米变化换算成每米变化 % 粘度与经验修正 eta = silica_viscosity(fparams.T_eff); C = fparams.C_exp; P = fparams.P_gauge; g = fparams.gamma; % 径向演化方程：表面张力收缩与内压膨胀的平衡 dRdz_m = -C * (g / t_m - P) / (4 * eta) * R_m / v_m; % 壁厚方程：由质量守恒对数微分得到 dtdz_m = -t_m * (dvdz_m / v_m) - (t_m / R_m) * dRdz_m; % 转回 mm/mm dydz = [dRdz_m * 1e3; dtdz_m * 1e3]; end

参数说明：dvdz 表示单位长度上的速度变化率，它刻画了软化区拉伸的剧烈程度，速度梯度过大时壁厚会急速减薄，往往对应着拉丝断纤。C 是经验修正系数，默认为 1.0 时是理想薄壁模型，实际炉况偏“软化偏早”时 C 小于 1，偏“硬化偏早”时 C 大于 1。P 取表压，单位帕斯卡，注意几百帕的气压在模型里已经能明显改变结果，这和常识里“几百帕不算压力”的直觉完全不同。

3.3 一组能直接跑通的最小输入与预期输出

下面这张表给出三组可以直接喂给上面脚本的参数组合，对应的输出大致在中空芯光纤的典型范围，可以在调试代码时用来验证改对了没有：

注意第三组工况容易触发数值警告，因为管腔内压不足时表面张力占主导，dRdz 变成较大的负值，壁厚方程里的耦合项会让结果迅速发散。这不一定是代码 bug，而是物理上这个参数组合本来就不该做中空芯光纤，更接近拉实心棒的工况。遇到 ode45 报错，先检查是不是 P 太小或者 t 掉到了下限。

3.4 输出判读：用稳定性而不是单个数值

单次仿真跑完只得到一个壁厚数字，这个数字能不能信，要看演化曲线是否平稳。建议调用脚本后接着画两张图：R-z 曲线和 t-z 曲线。健康的曲线是入口段平缓、中间段快速收细、出口段趋于平稳。如果出口附近 t 还在以每毫米几个微米的速度变化，说明软化区长度给的太短，实际拉制中光纤还没冷却定型就被拉走了，壁厚波动会很大。我做工程师时判断一个参数组合可不可用，看的不是出口最终值，而是出口处导数的绝对值是否收敛到接近零。这个准则在实拉中对应着光纤是否在拉丝塔的冷却段之前完成定型。

% 稳定性判据：出口处壁厚变化率的绝对值 dydz_end = draw_ode(z(end), y(end, :), fparams); if abs(dydz_end(2)) < 0.05 disp('出口壁厚趋于稳定'); else disp('出口壁厚仍在快速变化，需调整软化区长度或速度差'); end

这段小脚本直接复用 draw_ode 求出口导数，不需要重新积分。判断阈值 0.05 mm/mm 是经验值，对应壁厚在最后一个毫米距离内变化小于 50 纳米，足够满足反谐振光纤的容差要求。不同拉丝塔冷却距离不同，这个阈值要按实物壁厚波动数据重新标定。

4. 从试拉数据反演拉制参数：标定 Matlab 模型

4.1 正演够了，但模型和炉子之间还差一个标定

正演模型默认“等效温度”就是炉温、软化区内气压等于冷端表压，这两个假设在实验里几乎不成立。炉内热电偶测到的是加热丝附近的气氛温度，玻璃中心的实际温度受辐射吸收和气流冷却影响，通常偏低几十度；气压从冷端到热端要经过很长的进气管路，高温段气体膨胀后实际压力也会偏离表压。硬用正演模型去算目标参数，算出来的壁厚和实测经常对不上，这时就该进入反演标定环节。

反演思路很简单：拉三到五根不同参数组合的样品，用扫描电镜或白光干涉仪实测每根光纤的壁厚和外径，然后调整模型里测不准的少数几个参数，让仿真输出和全部实测值同时吻合。这个操作在有 matlab 优化工具箱的机器上，用 lsqnonlin 就能完成，不需要额外的数据拟合工具。

4.2 构造最小二乘目标函数与 lsqnonlin 拟合脚本

下面代码把可测的输出归一化后做为残差项，防止壁厚和外径数量级差异过大导致拟合偏向某一项。拟合对象只选两个参数：等效温度 T_eff 和修正系数 C_exp，气压和速度都用给定值。

% calibrate_draw_model.m % 用多根试拉样品的实测几何数据反演模型参数 % 实测数据: [外径 um, 壁厚 um] meas_D = [256, 3.1; 212, 0.9; 230, 1.7]; meas_t = [3.1, 0.9, 1.7]; % 对应的拉制参数 exp_params = [1980, 300, 8, 10; 2050, 800, 12, 20; 2000, 400, 10, 14]; % 初始猜测: 等效温度, 经验系数 x0 = [2020, 1.0]; % 目标函数: 归一化残差平方和 fun = @(x) residual_func(x, exp_params, meas_D, meas_t); % lsqnonlin 拟合 options = optimoptions('lsqnonlin', 'Display', 'iter', ... 'Algorithm', 'trust-region-reflective'); lb = [1900, 0.3]; ub = [2150, 3.0]; [x_fit, resnorm] = lsqnonlin(fun, x0, lb, ub, options); fprintf('拟合等效温度: %.1f degC\n', x_fit(1)); fprintf('拟合经验系数: %.3f\n', x_fit(2));

参数说明：exp_params 每行对应一个试拉条件，列顺序是等效温度初猜、管内表压、进给速度、收丝速度。实测外径和壁厚要放在同一组变量里传给残差函数，残差函数内部调用正演脚本的模型部分，把仿真结果和实测结果做归一化差值。lb 和 ub 给出了拟合参数的物理可行范围，等效温度低于 1900℃ 时石英粘度太大拉不动，高于 2150℃ 时粘度太低无法成形。经验系数超过 3 说明模型结构有问题，不是调系数能救回来的。

function res = residual_func(x, exp_params, meas_D, meas_t) % 残差函数：每个样品的仿真值与实测值之差，除以实测值做归一化 res = []; for k = 1:size(exp_params, 1) % 用初猜参数跑一次正演 T_eff = x(1); C_exp = x(2); P_gauge = exp_params(k, 2); v_p = exp_params(k, 3); v_draw = exp_params(k, 4); [D_sim, t_sim] = run_forward(T_eff, P_gauge, v_p, v_draw, C_exp); res = [res; (D_sim - meas_D(k)) / meas_D(k); (t_sim - meas_t(k)) / meas_t(k)]; end end

这里 run_forward 就是把第 3 章主程序里从参数到输出的部分封装成一个只返回外径和壁厚的函数，反演脚本里不需要绘图和判据。注意每次拟合循环都会调用多次正演，如果软化区网格很密，单次求解几百毫秒，5 个样品 20 轮迭代也就几十秒，完全在可接受范围内。

4.2.1 分步拟合避免参数退化

同时拟合 T_eff 和 C_exp 时，如果初始值给得太离谱，lsqnonlin 可能收敛到一个壁厚正确但外径错误的局部解。我一般先固定 C_exp=1.0 只拟合 T_eff，让仿真外径贴近实测外径，再把外径残差放到一边，固定 T_eff 拟合 C_exp 匹配壁厚。这样两步走虽然多跑几次正演，但稳定性好得多。两个参数同时拟合只有在样品数量不少于四个、且覆盖压力大范围变化时才建议使用。

4.3 标定后怎么找目标壁厚对应的拉制参数

模型标定完成后，正演就靠谱了。下一步是用它反向搜索目标壁厚对应的温度和压力组合。这一层我通常用 fminsearch 包一个目标函数，把“仿真壁厚偏离目标壁厚的差值”作为优化目标，同时约束外径不能超出指定范围。

% 反求目标壁厚对应的气压 % 设定目标 t_target = 1.3; % 目标壁厚 [um] D_target = 220; % 目标外径 [um] P_search = fminsearch(@(P) cost_target(P, t_target, D_target, T_fit, C_fit), 600); function cost = cost_target(P, t_target, D_target, T_fit, C_fit) % 固定温度和经验系数，只扫气压 [D_sim, t_sim] = run_forward(T_fit, P, v_p, v_draw, C_fit); cost = ((t_sim - t_target)/t_target)^2 + ... ((D_sim - D_target)/D_target)^2; end

气压的搜索初始值给 600 Pa 是经验值，因为前一轮试拉在这个气压附近。搜索结束后要把目标函数里的残差数值打出来看一眼，如果 cost 还在 1e-3 量级以上，说明温度和气压的参数组合下根本达不到目标几何，需要回炉重新设计预制棒尺寸，而不是继续压搜索算法。

4.4 标定结果的常见失败与判定表

反演拟合不是每次都能成功，我列几个高频失败场景和它们的表现：

烧硅的经验判断标准是：拟合后的残差达到实测重复性以内就算合格，一般归一化残差小于 0.03 即可认为模型可用于参数指导。如果残差一直压不下去，不要继续加迭代次数，那是模型结构问题。

5. 用实测回代与参数敏感性矩阵打磨拉制参数

5.1 两种实测壁厚来源，回代方式不同

模型标定和验证依赖的壁厚数据，工程上主要来自两个渠道：扫描电镜横截面照片和白光干涉仪。两者口径不同要注意：扫描电镜直接测几何壁厚，回代仿真模型最直观；白光干涉仪测的是光纤轴向的干涉光谱，通过反谐振谷波长反推光学厚度，它反映的是 n×t，受折射率色散影响，回代时要把光谱反演出的光学厚度除以工作波长的折射率再和模型对比。用错了口径会让标定系数产生系统性偏差，在现场跟踪时特别容易踩。

5.2 参数敏感性矩阵：找出当前工况的“命门”

标定完成后，我会用有限差分法快速算一遍参数敏感性，判断当前工作点附近哪个参数对壁厚的影响最大。做法是围绕标定好的参数组合，逐个小幅扰动 T_eff、P、v_p，记录壁厚变化量，整理成归一化敏感性矩阵。

% sensitivity_sweep.m % 围绕当前拉制参数做参数敏感性分析 base = [2050, 600, 10, 15]; % T_eff, P, v_p, v_draw delta = [10, 50, 1, 1]; % 扰动步长 fields = {'T_eff', 'P_gauge', 'v_p', 'v_draw'}; Dt = zeros(1, 4); for k = 1:4 p_hi = base; p_hi(k) = p_hi(k) + delta(k); p_lo = base; p_lo(k) = p_lo(k) - delta(k); [~, t_hi] = run_forward(p_hi(1), p_hi(2), p_hi(3), p_hi(4), C_fit); [~, t_lo] = run_forward(p_lo(1), p_lo(2), p_lo(3), p_lo(4), C_fit); Dt(k) = (t_hi - t_lo) / (2 * delta(k)); end % 输出: 每变化一个单位参数，壁厚的变化量 um disp(Dt);

这个脚本会输出四个数字，单位是 μm/单位参数。通常在薄壁工况下，敏感性矩阵里温度那一项会远大于其他项，这意味着当前炉温只要抖 10℃，壁厚就超差，再想靠气压和速度去补是补不回来的。这时最优策略不是继续调参数，而是先解决炉温稳定性，比如检查加热丝老化、气流波动。

5.3 一个效果明显的实战技巧：三点标定法

最后分享一个在现场验证过的技巧。每次换预制棒批次或大修炉子后，不要直接拉目标参数，先拉三根校准样：第一根用中等压力和中等速度，第二根只把气压提高一倍，第三根只把收丝速度提高一倍，每根记录实际外径和壁厚。这三根样品的实测数据正好构成 2.4 节和 4.3 节说的分步拟合输入：第一根定 T_eff，第二根和第一根的壁厚差定 P 的有效系数，第三根和第一根的外径差定速度梯度，三根数据互不打架，一次标定就能让模型误差收敛到可用范围。这个做法比一次性拉七八根正交试验样省一半预制棒，标定出来的参数也只针对当前炉况，换炉子必须重做一遍。

本文还有配套的精品资源，点击获取