简介这份资源是一套基于 MATLAB 的气动光学效应仿真程序面向大气光学、激光传输与成像探测领域的科研人员和学生。气动光学关注大气湍流、温度梯度等非均匀条件对光传播的干扰本资源重点模拟高斯光束、涡旋光束及环形光束在湍流大气中的传输特性。涡旋光束具有螺旋相位结构环形光束中心为暗区它们在大气中的稳定性与畸变行为是当前研究热点。压缩包内共 3 个 m 脚本文件大小仅 2KB三个脚本分别承担光束传播计算、多截面光场系列分析和动态 GIF 可视化结构精简便于直接运行和二次修改。已有 299 人学习下载。运行这批脚本可直观观察不同光束在大气中的强度分布演化、相位畸变和光束扩展现象帮助理解大气对光场质量的影响机制为激光通信、遥感探测等系统设计提供初步参考适合作为课程设计或科研预研的入门模板。1. 气动光学程序在MATLAB里到底解决什么问题高超声速飞行器的光学窗口、红外导引头和激光定向器会遇到同一个问题高速流场在光窗表面形成边界层和剪切层这些区域的密度不均匀会让光束发生偏折波前产生畸变。这个现象的量化评估工程上离不开把CFD密度场或风洞实验数据转换成光程差OPD、波前参数和Strehl比这一步。气动光学程序就是搭起这条数据链路的工具。用MATLAB写气动光学效应仿真门槛不高数据处理、频谱分析和可视化都顺手修改参数再跑一遍的周期也短。这里按理论模型、程序实现、参数验证和动态扩展的顺序讲清楚整个体系的搭建路径适合正在做气动光学评估、光机热控或者想快速评估湍流对成像系统影响的工程人员。2. 气动光学效应的物理模型从密度起伏到波前畸变写程序之前先把物理模型的链条理清湍流造成密度起伏密度起伏通过Gladstone-Dale关系变换成折射率变化折射率沿光束路径积分得到光程差OPD再由波长换算成相位畸变。链条上的每一步都有对应的MATLAB表达式。哪个环节出错调试时最先暴露的就是OPD数值偏大偏小或者Strehl比不符合物理直觉。2.1 Gladstone-Dale公式是程序的数据红线从空气动力学角度看流场输出的是压力、速度和温度从光学角度看只有折射率直接作用于光传播。气体折射率与密度的关系由Gladstone-Dale公式给出n 1 K_GD · ρ对空气在0.42.0 μm波段Gladstone-Dale常数K_GD约为2.23×10⁻⁴ m³/kg随波长和气体组分略有变化。这个换算出现在程序的数据输入和输出两端把CFD密度场转成折射率场要乘K_GD把光程差反推到密度起伏时也要除K_GD。单位不统一是这条线上最常见的问题密度用g/cm³而K_GD用m³/kg光程差直接差三个数量级。% 密度场(kg/m^3)到折射率场, 逐元素运算 n ones(size(rho)) K_GD .* rho;这里的rho可以是二维或三维数组ones(size(rho))保证折射率场的维度不变。如果rho是风洞数据导出的归一化密度必须先乘回参考密度再代入公式这一步漏掉会直接污染后面所有统计量。2.2 光程差与相位畸变的换算光沿z方向穿过厚度δ的扰动区域时光程OPL的定义是沿路径的折射率积分。实际操作中扣掉孔径内的平均活塞项得到反映波前空间畸变的OPDOPL(x,y) ∫₀^δ n(x,y,z) dzOPD(x,y) OPL(x,y) − ⟨OPL⟩⟨OPL⟩是孔径上的空间平均。由OPD换相位畸变直接乘波数即可φ(x,y) (2π/λ) · OPD(x,y)同样的OPD在0.8 μm近红外下比在10.6 μm长波红外下的相位误差大一个量级这就是气动光学对短波长系统更敏感的原因。Strehl比是评价效应的核心指标小扰动近似下等于exp(−σ_φ²)更严谨的写法是复指数平均的模方% 从OPD场求相位畸变与Strehl比 phase_map 2*pi / lambda * opd; strehl abs( mean( exp(1i*phase_map), all) )^2;mean的all选项在MATLAB 2018b之后可用它会沿所有维度取平均。程序里用严格形式算Strehl不需要事先判断相位误差大小这行代码比近似式更稳。2.3 气动光学用Kolmogorov谱还是von Kármán谱湍流密度起伏的统计特性是合成随机流场模型的基础。对充分发展的湍流折射率起伏的三维功率谱常用Kolmogorov形式Φ_n(κ) 0.033·C_n²·κ^(−11/3)C_n²是折射率结构常数单位m^(−2/3)。气动光学问题里C_n²并不是常数而是随边界层动量厚度、当地密度和温度梯度变化的量通常由CFD换算或半经验公式给出。实际合成随机场时如果完全按这个纯幂律谱低频能量会发散所以程序中更稳的是von Kármán谱Φ_n(κ) 0.033·C_n²·(κ² κ₀²)^(−11/6)·exp(−κ²/κ_m²)其中κ₀2π/L₀L₀是外尺度与边界层厚度同一量级κ_m5.92/ℓ₀ℓ₀是内尺度对应Kolmogorov微尺度。外尺度控制OPD场的低频倾斜和活塞量内尺度压掉高频尾翼两者取错会让频谱形状整体变形。谱参数数学表达对OPD程序的影响外尺度L₀κ₀ 2π/L₀决定低频倾斜能量内尺度ℓ₀κ_m 5.92/ℓ₀决定高频截止位置结构常数C_n²谱幅度决定OPD总体幅值湍流层厚δ积分路径长度与OPD方差近线性2.4 从CFD密度场反推C_n²的程序化处理实际工程里C_n²并不是一个直接给定的常数。拿到CFD密度场之后需要先沿等z面做Gladstone-Dale换算得到折射率场再提取折射率结构函数D_n(r)用D_n(r) C_n²·r^(2/3)做log-log域直线拟合才能得到这个平面上的当地C_n²。这一步放到MATLAB里就是对结构函数曲线做一阶线性回归斜率应接近2/3由截距反解出C_n²。这也是气动光学程序区别于大气光学程序的地方大气湍流通常把C_n²当成路径上的慢变量气动光学则必须逐层处理边界层内的强梯度变化。3. 用MATLAB写气动光学仿真程序的基本流程程序的最小可运行版本按四段组织参数定义、OPD场生成、指标计算、可视化。下面给出一套可以直接落地的框架模块拆分保持独立后面替换CFD数据源或改成时序仿真时不至于重写主程序。3.1 主程序框架与参数定义把参数集中放在主程序顶部改动时不用去找散落各处的常数。% main_aero_optic.m % 气动光学效应最小可运行程序 clear; clc; close all; % ---- 光学参数 ---- lambda 1.064e-6; % 工作波长(m), 取Nd:YAG激光 K_GD 2.23e-4; % Gladstone-Dale常数(m^3/kg), 密度场输入时使用 % ---- 湍流层参数 ---- Cn2 5e-13; % 折射率结构常数(m^-2/3) delta 0.01; % 沿光路的湍流层厚度(m) % ---- 计算域与网格 ---- L 0.05; % 计算域边长(m), 对应50mm口径 N 256; % 每边网格数 % ---- 生成OPD场 ---- opd generate_opd_phase_screen(N, L, Cn2, delta); % ---- 波前指标 ---- phase_map 2*pi/lambda * opd; strehl abs(mean(exp(1i*phase_map), all))^2; fprintf(OPD RMS %.3f um\n, std(opd(:))*1e6); fprintf(Strehl %.4f\n, strehl);选型逻辑说明C_n²取5×10⁻¹³是临近空间高超声速边界层中常见的量级做配平和趋势研究可以先用它起步。δ取0.01 m对应薄边界层假设实际使用时应与CFD边界层厚度对齐。主程序里保留K_GD变量是给第3.3节接入密度场数据时预留的纯统计合成路径用不到它。3.2 谱反演法生成OPD屏OPD屏的生成采用频域滤波法也叫功率谱反演法。先按功率谱形状构造频域滤波函数再给每个频点配上复高斯随机数通过逆傅里叶变换回到空间域。function opd generate_opd_phase_screen(N, L, Cn2, delta) % 生成符合Kolmogorov谱的二维OPD屏 % 输入: N网格数, L计算域边长(m), Cn2结构常数(m^-2/3), delta层厚(m) % 输出: opd光程差分布(m) fx (-N/2 : N/2-1) / L; % 空间频率(cycles/m) [fx, fy] meshgrid(fx); f sqrt(fx.^2 fy.^2); f(N/21, N/21) 1e-9; % 对零频做保护 % 二维OPD功率谱密度, 薄层近似, 连续谱单位m^4 Phi_opd 0.033 * Cn2 * delta * (2*pi*f).^(-11/3); % 频域采样间隔与离散方差 df 1/L; random_amp (randn(N) 1i*randn(N)) .* sqrt(Phi_opd * df^2); % 逆FFT, 乘N^2还原傅里叶级数系数 opd real( ifft2( ifftshift(random_amp) ) ) * N^2; % 去除平均活塞 opd opd - mean(opd(:)); end这段代码有三个容易看漏的点。第一f(N/21,N/21)1e-9是保护零频用的不处理会因除零产生NaN或者出现一个很大的直流伪影。第二ifftshift把零频挪回矩阵左上角与fftshift配对才能让频率坐标和数组索引正确对应。第三sqrt(Phi_opd * df^2)把连续功率谱变成离散频率槽内的幅度ifft2自带的1/N²归一化需要用乘N²抵消。最后减均值去掉活塞项它对成像质量没有影响只改变整体光程常数。3.3 从CFD三维密度场积分得到OPD功率谱反演法的局限是没有真实空间结构。工程上拿到CFD密度场后更常用的路径是直接沿光束方向积分。这里建议单独封装一个函数后续无论换算例还是改网格都只需要动这一处。function [OPL, OPD] compute_opd_from_density(rho3d, z_grid, K_GD) % 从三维密度场计算光程和光程差 % rho3d: Nx×Ny×Nz密度场(kg/m^3) % z_grid: Nz×1, 光束方向坐标(m) % K_GD: Gladstone-Dale常数(m^3/kg) n3d 1 K_GD .* rho3d; % 折射率场 dz z_grid(2) - z_grid(1); % 均匀网格间距 OPL trapz(n3d, 3) * dz; % 沿第三维积分 OPD OPL - mean(OPL(:), all); end这里默认z_grid是均匀网格。如果CFD在边界层内做了加密z_grid不是等间距的需要把trapz改成trapz(z_grid(:), n3d, 3)的形式。trapz沿第三维的基本用法是对每个[x,y]像素做一维数值积分返回尺寸为Nx×Ny的OPL矩阵。用CFD数据时另一个容易被忽略的点是密度是当地静密度还是总密度导出的无量纲密度必须乘回自由流密度否则OPD整体会偏离正确量级。3.4 指标计算与可视化把OPD分布和直方图打印出来是判断程序是否跑通的第一道检查。可视化代码放在主程序末尾figure(Name, Aero-Optic OPD); subplot(1,2,1); imagesc(opd*1e6); axis image; colorbar; xlabel(x网格); ylabel(y网格); title(OPD分布 (um)); subplot(1,2,2); histogram(opd(:)*1e6, 50); xlabel(OPD (um)); ylabel(像素数); title(OPD直方图);OPD的均方根值可以直接从std得到直方图看分布形态是否接近高斯。气动光学的OPD统计在多数情况下近似高斯如果直方图明显偏斜或出现多峰通常说明合成场里混入了过大的低频成分或直流残留。程序输出的核心指标可以归纳为下表后续做参数扫描时统一按这套指标做回归。输出量符号单位读取方式光程差均方根OPD RMSμmstd(opd(:))*1e6Strehl比SR无量纲abs(mean(exp(1i*phi),all))^2相位结构函数D_φ(r)rad²第4.2节代码计算4. 气动光学仿真程序的参数验证与调试程序能跑只是第一步。仿真参数和真实流场对不上输出结果再漂亮也不能用于光学设计。这一章按参数表、结构函数验证、频域合成陷阱、确定性用例四个方向展开。4.1 关键参数表与调整原则气动光学程序中真正决定输出量级的参数并不多列成一张表便于快速定位问题。参数符号常见量级输出异常的典型表现折射率结构常数C_n²1e-141e-12OPD RMS整体偏大或偏小湍流层厚度δ550 mmOPD方差随δ近似线性变化计算域边长L50200 mm低频倾斜不足或伪周期条纹网格数N1281024高频细节不足或噪声过重工作波长λ0.810.6 μmStrehl比变化梯度明显调参的原则是先固定几何参数L和N只扫C_n²和δ把OPD RMS的曲线标定出来。这两个参数一个决定幅值一个决定积分长度与OPD标准差的关系近似为一次方。网格继续增大到N1024以上时计算量按O(N²)上涨但OPD RMS的增量往往已经进入噪声区不必要盲目增加网格。4.2 用波前结构函数验证统计特性验证程序是否正确光看RMS不够更严格的做法是让生成的OPD场满足Kolmogorov湍流的相位结构函数。结构函数定义是两点相位差的均方D_φ(r) ⟨[φ(xr) − φ(x)]²⟩对Kolmogorov湍流理论值等于6.88(r/r₀)^(5/3)其中r₀是大气相干长度。对均匀薄层由C_n²和δ可以求出r₀ [0.423·(2π/λ)²·C_n²·δ]^(−3/5)。验证代码可以直接复用上一阶段生成的phase_map% validate_aero_optic.m % 从OPD场计算相位结构函数并与理论曲线比较 r0 (0.423 * (2*pi/lambda)^2 * Cn2 * delta)^(-3/5); max_lag min(64, floor(N/4)); D_phi_sim zeros(1, max_lag); for lag 1:max_lag diff_field phase_map(1:end-lag, :) - phase_map(1lag:end, :); D_phi_sim(lag) mean(diff_field(:).^2); end r_lag (1:max_lag) * (L/N); D_phi_theory 6.88 * (r_lag / r0).^(5/3); loglog(r_lag, D_phi_sim, o-); hold on; loglog(r_lag, D_phi_theory, r--); xlabel(分离距离 r (m)); ylabel(相位结构函数 D_\phi (rad^2)); legend(仿真, 6.88(r/r0)^{5/3}, Location, northwest);验证的思路是结构函数只与两点间距有关不受平均活塞影响也不依赖绝对相位值。仿真曲线在与理论线重合的范围内程序统计上是可靠的。低频段出现偏差是因为计算域L截断了外尺度高频段是因为网格分辨率截断了内尺度这两处偏差本身就是程序中尺度参数设置的反映。4.3 频域合成时的三个常见坑第一个坑是计算域太小。L只有20 mm而实际口径是50 mm合成OPD场缺少足够大的低频起伏倾斜项不足Strehl比偏高。修法是让L至少覆盖口径的1.5到2倍否则波前斜率和整体像差都会失真。第二个坑是网格数N不够。N64时高频空间频率上限低OPD场太平滑结构函数在短间距处上不去。一般取N256起步做参数扫描时视情况降为128即可。第三个坑是零频保护不当。f(N/21,N/21)1e-9这行如果不写零频处幅度无穷大ifft之后会出现一大片常值偏移虽然减均值能扣掉活塞但欠采样的低频残余会污染结构函数中段。这个保护看起来不起眼却是每个从大气光学转做气动光学的人最容易漏的地方。4.4 先做确定性用例再做统计用例最容易上手的调试方法是先构造一个确定性OPD场而不是直接扔随机场。比如把整个孔径设置成常数0.05 μm的OPD程序计算出的Strehl必须等于1再构造一个倾斜面Strehl会随倾斜量衰减可以与解析的sinc调制结果对比。如果程序连常数和倾斜都不正确后面的统计结果没有意义。% 确定性用例: 常数OPD必须给出Strehl1 N 256; opd_flat 0.05e-6 * ones(N); % 常数OPD, 只贡献活塞 strehl_flat abs(mean(exp(1i*2*pi/lambda*opd_flat), all))^2; fprintf(Strehl for flat OPD %.6f\n, strehl_flat);这个用例能同时排查两件事一是相位换算公式是否写错二是exp和mean的维度处理是否一致。确定性用例通过后再进入随机OPD屏的统计验证问题定位会快很多。5. 进阶把静态气动光学程序扩展为动态时序仿真工程上最关心的往往不是单帧OPD而是时间序列。聚焦光束在湍流流场里的抖动、气动光学效应的时间频率谱这些都要靠动态仿真来评估。5.1 Taylor冻结假设生成时间序列常见做法是使用Taylor冻结涡假设把空间场平移成时间序列。假设流场以当地平均速度U平流时间间隔dt对应的空间位移是U·dt把上一帧OPD平移几个像素就能得到下一帧% 帧间平移生成时序OPD dx L / N; shift_px round(U_flow * dt / dx); opd_next circshift(opd, shift_px, 1);参数说明U_flow取边界层外缘速度或当地对流速度dt是采样间隔。把shift_px控制在12像素能保证帧间连续太大则帧间跳跃明显太小则帧间几乎不变化。这个近似只在对流马赫数不高、湍流演变时间远大于对流通过时间的情况下成立。速度梯度较大的边界层底部Taylor假设会有偏差更精细的做法是对每帧叠加一个独立的小随机场来模拟湍流演化。5.2 与光学仿真软件和实测数据的衔接时序OPD生成后保存为通用格式即可供光学设计软件读取% 保存OPD到mat和csv save(opd_sequence.mat, opd_seq, params); writematrix(opd, opd_frame.csv);把OPD帧导入Zemax OpticStudio或Code V时需要注意单位程序输出用米光学软件一般按微米或毫米读取。实测数据可以先做倾斜剔除再和仿真OPD的自相关时间对比检查时间尺度是否一致。动态气动光学程序的核心始终是帧间相关性和时间功率谱的衰减不能靠肉眼判断时间序列是否合理要把时间功率谱打出来与流场频谱特征对照统计范围内正确的仿真结果才真正可用。本文还有配套的精品资源点击获取 MATLAB气动光学仿真：从密度场到OPD与Strehl比的全流程实现
简介这份资源是一套基于 MATLAB 的气动光学效应仿真程序面向大气光学、激光传输与成像探测领域的科研人员和学生。气动光学关注大气湍流、温度梯度等非均匀条件对光传播的干扰本资源重点模拟高斯光束、涡旋光束及环形光束在湍流大气中的传输特性。涡旋光束具有螺旋相位结构环形光束中心为暗区它们在大气中的稳定性与畸变行为是当前研究热点。压缩包内共 3 个 m 脚本文件大小仅 2KB三个脚本分别承担光束传播计算、多截面光场系列分析和动态 GIF 可视化结构精简便于直接运行和二次修改。已有 299 人学习下载。运行这批脚本可直观观察不同光束在大气中的强度分布演化、相位畸变和光束扩展现象帮助理解大气对光场质量的影响机制为激光通信、遥感探测等系统设计提供初步参考适合作为课程设计或科研预研的入门模板。1. 气动光学程序在MATLAB里到底解决什么问题高超声速飞行器的光学窗口、红外导引头和激光定向器会遇到同一个问题高速流场在光窗表面形成边界层和剪切层这些区域的密度不均匀会让光束发生偏折波前产生畸变。这个现象的量化评估工程上离不开把CFD密度场或风洞实验数据转换成光程差OPD、波前参数和Strehl比这一步。气动光学程序就是搭起这条数据链路的工具。用MATLAB写气动光学效应仿真门槛不高数据处理、频谱分析和可视化都顺手修改参数再跑一遍的周期也短。这里按理论模型、程序实现、参数验证和动态扩展的顺序讲清楚整个体系的搭建路径适合正在做气动光学评估、光机热控或者想快速评估湍流对成像系统影响的工程人员。2. 气动光学效应的物理模型从密度起伏到波前畸变写程序之前先把物理模型的链条理清湍流造成密度起伏密度起伏通过Gladstone-Dale关系变换成折射率变化折射率沿光束路径积分得到光程差OPD再由波长换算成相位畸变。链条上的每一步都有对应的MATLAB表达式。哪个环节出错调试时最先暴露的就是OPD数值偏大偏小或者Strehl比不符合物理直觉。2.1 Gladstone-Dale公式是程序的数据红线从空气动力学角度看流场输出的是压力、速度和温度从光学角度看只有折射率直接作用于光传播。气体折射率与密度的关系由Gladstone-Dale公式给出n 1 K_GD · ρ对空气在0.42.0 μm波段Gladstone-Dale常数K_GD约为2.23×10⁻⁴ m³/kg随波长和气体组分略有变化。这个换算出现在程序的数据输入和输出两端把CFD密度场转成折射率场要乘K_GD把光程差反推到密度起伏时也要除K_GD。单位不统一是这条线上最常见的问题密度用g/cm³而K_GD用m³/kg光程差直接差三个数量级。% 密度场(kg/m^3)到折射率场, 逐元素运算 n ones(size(rho)) K_GD .* rho;这里的rho可以是二维或三维数组ones(size(rho))保证折射率场的维度不变。如果rho是风洞数据导出的归一化密度必须先乘回参考密度再代入公式这一步漏掉会直接污染后面所有统计量。2.2 光程差与相位畸变的换算光沿z方向穿过厚度δ的扰动区域时光程OPL的定义是沿路径的折射率积分。实际操作中扣掉孔径内的平均活塞项得到反映波前空间畸变的OPDOPL(x,y) ∫₀^δ n(x,y,z) dzOPD(x,y) OPL(x,y) − ⟨OPL⟩⟨OPL⟩是孔径上的空间平均。由OPD换相位畸变直接乘波数即可φ(x,y) (2π/λ) · OPD(x,y)同样的OPD在0.8 μm近红外下比在10.6 μm长波红外下的相位误差大一个量级这就是气动光学对短波长系统更敏感的原因。Strehl比是评价效应的核心指标小扰动近似下等于exp(−σ_φ²)更严谨的写法是复指数平均的模方% 从OPD场求相位畸变与Strehl比 phase_map 2*pi / lambda * opd; strehl abs( mean( exp(1i*phase_map), all) )^2;mean的all选项在MATLAB 2018b之后可用它会沿所有维度取平均。程序里用严格形式算Strehl不需要事先判断相位误差大小这行代码比近似式更稳。2.3 气动光学用Kolmogorov谱还是von Kármán谱湍流密度起伏的统计特性是合成随机流场模型的基础。对充分发展的湍流折射率起伏的三维功率谱常用Kolmogorov形式Φ_n(κ) 0.033·C_n²·κ^(−11/3)C_n²是折射率结构常数单位m^(−2/3)。气动光学问题里C_n²并不是常数而是随边界层动量厚度、当地密度和温度梯度变化的量通常由CFD换算或半经验公式给出。实际合成随机场时如果完全按这个纯幂律谱低频能量会发散所以程序中更稳的是von Kármán谱Φ_n(κ) 0.033·C_n²·(κ² κ₀²)^(−11/6)·exp(−κ²/κ_m²)其中κ₀2π/L₀L₀是外尺度与边界层厚度同一量级κ_m5.92/ℓ₀ℓ₀是内尺度对应Kolmogorov微尺度。外尺度控制OPD场的低频倾斜和活塞量内尺度压掉高频尾翼两者取错会让频谱形状整体变形。谱参数数学表达对OPD程序的影响外尺度L₀κ₀ 2π/L₀决定低频倾斜能量内尺度ℓ₀κ_m 5.92/ℓ₀决定高频截止位置结构常数C_n²谱幅度决定OPD总体幅值湍流层厚δ积分路径长度与OPD方差近线性2.4 从CFD密度场反推C_n²的程序化处理实际工程里C_n²并不是一个直接给定的常数。拿到CFD密度场之后需要先沿等z面做Gladstone-Dale换算得到折射率场再提取折射率结构函数D_n(r)用D_n(r) C_n²·r^(2/3)做log-log域直线拟合才能得到这个平面上的当地C_n²。这一步放到MATLAB里就是对结构函数曲线做一阶线性回归斜率应接近2/3由截距反解出C_n²。这也是气动光学程序区别于大气光学程序的地方大气湍流通常把C_n²当成路径上的慢变量气动光学则必须逐层处理边界层内的强梯度变化。3. 用MATLAB写气动光学仿真程序的基本流程程序的最小可运行版本按四段组织参数定义、OPD场生成、指标计算、可视化。下面给出一套可以直接落地的框架模块拆分保持独立后面替换CFD数据源或改成时序仿真时不至于重写主程序。3.1 主程序框架与参数定义把参数集中放在主程序顶部改动时不用去找散落各处的常数。% main_aero_optic.m % 气动光学效应最小可运行程序 clear; clc; close all; % ---- 光学参数 ---- lambda 1.064e-6; % 工作波长(m), 取Nd:YAG激光 K_GD 2.23e-4; % Gladstone-Dale常数(m^3/kg), 密度场输入时使用 % ---- 湍流层参数 ---- Cn2 5e-13; % 折射率结构常数(m^-2/3) delta 0.01; % 沿光路的湍流层厚度(m) % ---- 计算域与网格 ---- L 0.05; % 计算域边长(m), 对应50mm口径 N 256; % 每边网格数 % ---- 生成OPD场 ---- opd generate_opd_phase_screen(N, L, Cn2, delta); % ---- 波前指标 ---- phase_map 2*pi/lambda * opd; strehl abs(mean(exp(1i*phase_map), all))^2; fprintf(OPD RMS %.3f um\n, std(opd(:))*1e6); fprintf(Strehl %.4f\n, strehl);选型逻辑说明C_n²取5×10⁻¹³是临近空间高超声速边界层中常见的量级做配平和趋势研究可以先用它起步。δ取0.01 m对应薄边界层假设实际使用时应与CFD边界层厚度对齐。主程序里保留K_GD变量是给第3.3节接入密度场数据时预留的纯统计合成路径用不到它。3.2 谱反演法生成OPD屏OPD屏的生成采用频域滤波法也叫功率谱反演法。先按功率谱形状构造频域滤波函数再给每个频点配上复高斯随机数通过逆傅里叶变换回到空间域。function opd generate_opd_phase_screen(N, L, Cn2, delta) % 生成符合Kolmogorov谱的二维OPD屏 % 输入: N网格数, L计算域边长(m), Cn2结构常数(m^-2/3), delta层厚(m) % 输出: opd光程差分布(m) fx (-N/2 : N/2-1) / L; % 空间频率(cycles/m) [fx, fy] meshgrid(fx); f sqrt(fx.^2 fy.^2); f(N/21, N/21) 1e-9; % 对零频做保护 % 二维OPD功率谱密度, 薄层近似, 连续谱单位m^4 Phi_opd 0.033 * Cn2 * delta * (2*pi*f).^(-11/3); % 频域采样间隔与离散方差 df 1/L; random_amp (randn(N) 1i*randn(N)) .* sqrt(Phi_opd * df^2); % 逆FFT, 乘N^2还原傅里叶级数系数 opd real( ifft2( ifftshift(random_amp) ) ) * N^2; % 去除平均活塞 opd opd - mean(opd(:)); end这段代码有三个容易看漏的点。第一f(N/21,N/21)1e-9是保护零频用的不处理会因除零产生NaN或者出现一个很大的直流伪影。第二ifftshift把零频挪回矩阵左上角与fftshift配对才能让频率坐标和数组索引正确对应。第三sqrt(Phi_opd * df^2)把连续功率谱变成离散频率槽内的幅度ifft2自带的1/N²归一化需要用乘N²抵消。最后减均值去掉活塞项它对成像质量没有影响只改变整体光程常数。3.3 从CFD三维密度场积分得到OPD功率谱反演法的局限是没有真实空间结构。工程上拿到CFD密度场后更常用的路径是直接沿光束方向积分。这里建议单独封装一个函数后续无论换算例还是改网格都只需要动这一处。function [OPL, OPD] compute_opd_from_density(rho3d, z_grid, K_GD) % 从三维密度场计算光程和光程差 % rho3d: Nx×Ny×Nz密度场(kg/m^3) % z_grid: Nz×1, 光束方向坐标(m) % K_GD: Gladstone-Dale常数(m^3/kg) n3d 1 K_GD .* rho3d; % 折射率场 dz z_grid(2) - z_grid(1); % 均匀网格间距 OPL trapz(n3d, 3) * dz; % 沿第三维积分 OPD OPL - mean(OPL(:), all); end这里默认z_grid是均匀网格。如果CFD在边界层内做了加密z_grid不是等间距的需要把trapz改成trapz(z_grid(:), n3d, 3)的形式。trapz沿第三维的基本用法是对每个[x,y]像素做一维数值积分返回尺寸为Nx×Ny的OPL矩阵。用CFD数据时另一个容易被忽略的点是密度是当地静密度还是总密度导出的无量纲密度必须乘回自由流密度否则OPD整体会偏离正确量级。3.4 指标计算与可视化把OPD分布和直方图打印出来是判断程序是否跑通的第一道检查。可视化代码放在主程序末尾figure(Name, Aero-Optic OPD); subplot(1,2,1); imagesc(opd*1e6); axis image; colorbar; xlabel(x网格); ylabel(y网格); title(OPD分布 (um)); subplot(1,2,2); histogram(opd(:)*1e6, 50); xlabel(OPD (um)); ylabel(像素数); title(OPD直方图);OPD的均方根值可以直接从std得到直方图看分布形态是否接近高斯。气动光学的OPD统计在多数情况下近似高斯如果直方图明显偏斜或出现多峰通常说明合成场里混入了过大的低频成分或直流残留。程序输出的核心指标可以归纳为下表后续做参数扫描时统一按这套指标做回归。输出量符号单位读取方式光程差均方根OPD RMSμmstd(opd(:))*1e6Strehl比SR无量纲abs(mean(exp(1i*phi),all))^2相位结构函数D_φ(r)rad²第4.2节代码计算4. 气动光学仿真程序的参数验证与调试程序能跑只是第一步。仿真参数和真实流场对不上输出结果再漂亮也不能用于光学设计。这一章按参数表、结构函数验证、频域合成陷阱、确定性用例四个方向展开。4.1 关键参数表与调整原则气动光学程序中真正决定输出量级的参数并不多列成一张表便于快速定位问题。参数符号常见量级输出异常的典型表现折射率结构常数C_n²1e-141e-12OPD RMS整体偏大或偏小湍流层厚度δ550 mmOPD方差随δ近似线性变化计算域边长L50200 mm低频倾斜不足或伪周期条纹网格数N1281024高频细节不足或噪声过重工作波长λ0.810.6 μmStrehl比变化梯度明显调参的原则是先固定几何参数L和N只扫C_n²和δ把OPD RMS的曲线标定出来。这两个参数一个决定幅值一个决定积分长度与OPD标准差的关系近似为一次方。网格继续增大到N1024以上时计算量按O(N²)上涨但OPD RMS的增量往往已经进入噪声区不必要盲目增加网格。4.2 用波前结构函数验证统计特性验证程序是否正确光看RMS不够更严格的做法是让生成的OPD场满足Kolmogorov湍流的相位结构函数。结构函数定义是两点相位差的均方D_φ(r) ⟨[φ(xr) − φ(x)]²⟩对Kolmogorov湍流理论值等于6.88(r/r₀)^(5/3)其中r₀是大气相干长度。对均匀薄层由C_n²和δ可以求出r₀ [0.423·(2π/λ)²·C_n²·δ]^(−3/5)。验证代码可以直接复用上一阶段生成的phase_map% validate_aero_optic.m % 从OPD场计算相位结构函数并与理论曲线比较 r0 (0.423 * (2*pi/lambda)^2 * Cn2 * delta)^(-3/5); max_lag min(64, floor(N/4)); D_phi_sim zeros(1, max_lag); for lag 1:max_lag diff_field phase_map(1:end-lag, :) - phase_map(1lag:end, :); D_phi_sim(lag) mean(diff_field(:).^2); end r_lag (1:max_lag) * (L/N); D_phi_theory 6.88 * (r_lag / r0).^(5/3); loglog(r_lag, D_phi_sim, o-); hold on; loglog(r_lag, D_phi_theory, r--); xlabel(分离距离 r (m)); ylabel(相位结构函数 D_\phi (rad^2)); legend(仿真, 6.88(r/r0)^{5/3}, Location, northwest);验证的思路是结构函数只与两点间距有关不受平均活塞影响也不依赖绝对相位值。仿真曲线在与理论线重合的范围内程序统计上是可靠的。低频段出现偏差是因为计算域L截断了外尺度高频段是因为网格分辨率截断了内尺度这两处偏差本身就是程序中尺度参数设置的反映。4.3 频域合成时的三个常见坑第一个坑是计算域太小。L只有20 mm而实际口径是50 mm合成OPD场缺少足够大的低频起伏倾斜项不足Strehl比偏高。修法是让L至少覆盖口径的1.5到2倍否则波前斜率和整体像差都会失真。第二个坑是网格数N不够。N64时高频空间频率上限低OPD场太平滑结构函数在短间距处上不去。一般取N256起步做参数扫描时视情况降为128即可。第三个坑是零频保护不当。f(N/21,N/21)1e-9这行如果不写零频处幅度无穷大ifft之后会出现一大片常值偏移虽然减均值能扣掉活塞但欠采样的低频残余会污染结构函数中段。这个保护看起来不起眼却是每个从大气光学转做气动光学的人最容易漏的地方。4.4 先做确定性用例再做统计用例最容易上手的调试方法是先构造一个确定性OPD场而不是直接扔随机场。比如把整个孔径设置成常数0.05 μm的OPD程序计算出的Strehl必须等于1再构造一个倾斜面Strehl会随倾斜量衰减可以与解析的sinc调制结果对比。如果程序连常数和倾斜都不正确后面的统计结果没有意义。% 确定性用例: 常数OPD必须给出Strehl1 N 256; opd_flat 0.05e-6 * ones(N); % 常数OPD, 只贡献活塞 strehl_flat abs(mean(exp(1i*2*pi/lambda*opd_flat), all))^2; fprintf(Strehl for flat OPD %.6f\n, strehl_flat);这个用例能同时排查两件事一是相位换算公式是否写错二是exp和mean的维度处理是否一致。确定性用例通过后再进入随机OPD屏的统计验证问题定位会快很多。5. 进阶把静态气动光学程序扩展为动态时序仿真工程上最关心的往往不是单帧OPD而是时间序列。聚焦光束在湍流流场里的抖动、气动光学效应的时间频率谱这些都要靠动态仿真来评估。5.1 Taylor冻结假设生成时间序列常见做法是使用Taylor冻结涡假设把空间场平移成时间序列。假设流场以当地平均速度U平流时间间隔dt对应的空间位移是U·dt把上一帧OPD平移几个像素就能得到下一帧% 帧间平移生成时序OPD dx L / N; shift_px round(U_flow * dt / dx); opd_next circshift(opd, shift_px, 1);参数说明U_flow取边界层外缘速度或当地对流速度dt是采样间隔。把shift_px控制在12像素能保证帧间连续太大则帧间跳跃明显太小则帧间几乎不变化。这个近似只在对流马赫数不高、湍流演变时间远大于对流通过时间的情况下成立。速度梯度较大的边界层底部Taylor假设会有偏差更精细的做法是对每帧叠加一个独立的小随机场来模拟湍流演化。5.2 与光学仿真软件和实测数据的衔接时序OPD生成后保存为通用格式即可供光学设计软件读取% 保存OPD到mat和csv save(opd_sequence.mat, opd_seq, params); writematrix(opd, opd_frame.csv);把OPD帧导入Zemax OpticStudio或Code V时需要注意单位程序输出用米光学软件一般按微米或毫米读取。实测数据可以先做倾斜剔除再和仿真OPD的自相关时间对比检查时间尺度是否一致。动态气动光学程序的核心始终是帧间相关性和时间功率谱的衰减不能靠肉眼判断时间序列是否合理要把时间功率谱打出来与流场频谱特征对照统计范围内正确的仿真结果才真正可用。本文还有配套的精品资源点击获取