Введение

Волоконно-оптические преобразователи широко применяются в высокоточных измерительных системах благодаря универсальности и устойчивости к внешним воздействиям [1]. Управление оптическим потоком осуществляется модулирующими устройствами предельного и поглощающего типов. Среди них выделяются растровые решетки – периодические структуры с переменным пропусканием, служащие для преобразования линейных и угловых перемещений в цифровой код. Такие решетки обеспечивают линейность характеристик, высокую стабильность и адаптируемость к требуемому диапазону измерений, а также устойчивы к вибрации и температурным колебаниям. Применение растровых преобразователей устраняет влияние нестабильности питающего напряжения и температуры на выходной сигнал, а амплитудно-логический метод обработки гарантирует высокую точность в широком температурном диапазоне. При микронных размерах элементов открытого измерительного канала доминирующую роль начинают играть дифракционные эффекты на круглом отверстии диафрагмы. Некорректный учет этих эффектов ведет к ошибкам позиционирования решетки и снижению точности, что обусловливает актуальность разработки строгой математической модели дифракции в ближней зоне и создание на ее основе вычислительного инструментария для инженерного проектирования.

Цель исследования – построение математической модели распространения света в открытом канале волоконно-оптического растрового преобразователя на основе интеграла Рэлея – Зоммерфельда, численное моделирование и создание программного комплекса, позволяющего проводить многовариантные вычислительные эксперименты для определения рациональных геометрических параметров преобразователя.

Материалы и методы исследования

Постановка задачи математического моделирования

На рис. 1 представлен волоконно-оптический растровый преобразователь. Показан случай пропускания светового потока через пропускающий участок решетки. В случае с поглощающей частью решетки световой поток поглощается полностью, соответственно сигнал с отводящего оптического волокна не поступает.

Согласно исследованиям [2] выявлено, что для данного преобразователя целесообразно расположить растровую решетку и отводящее оптическое волокно в зоне наибольшей интенсивности излучения и резкого перепада сигнала, то есть на расстоянии до xкр, так как для дальнейшего преобразования необходимо только наличие либо отсутствие светового потока на отводящем торце. Для оптического волокна с диаметром сердцевины dc = 200 мкм и апертурным углом θNA = 12° расстояние xкр составит примерно 470 мкм. Соответственно, возникает одна из сложнейших задач расположения всех элементов преобразователя в пределах данного расстояния с учетом затухания светового потока и влияния дифракции [3].

Рис. 1. Распределение светового потока в открытом канале волоконно-оптического преобразователя: ПОВ – подводящее оптическое волокно, ООВ – отводящее оптическое волокно, dc – диаметр сердцевины оптоволокна, dо – диаметр отверстия диафрагмы, xкр – дистанция формирования светового пучка, sк – толщина клеевого слоя, sд – толщина диафрагмы, sр – толщина растровой решетки, xок – длина открытого оптического канала, x1 – расстояние между диафрагмой и растровой решеткой, x2 – расстояние между растровой решеткой и ООВ Примечание: составлен авторами по результатам данного исследования

В данном случае используется диафрагма из фольги с круглым отверстием диаметром dо = 8 мкм, ограничивающая расходимость пучка. Прошедшее через нее излучение попадает на растровую решетку с шириной полосы пропускания 8 мкм, где в зависимости от локального коэффициента пропускания поглощается либо поступает в отводящее оптическое волокно [4].

Математическая модель дифракции

Большое значение в оптических системах имеет явление дифракции. Понятие дифракции в оптике связано с нарушением прямолинейности распространения света в среде с неоднородностями, такими как отверстия, края экранов и т. д. Согласно скалярной теории дифракции основными методами, применяемыми при численных расчетах, являются приближения Френеля и Фраунгофера, а также метод углового спектра [5]. В ближнем поле более строгим считается интеграл Рэлея – Зоммерфельда первого типа, который представляет собой решение скалярного волнового уравнения в полупространстве с граничным условием Дирихле и тесно связан с методом углового спектра. Данный интеграл используется при малых расстояниях распространения волнового фронта порядка десятков длин волн [6].

Распределение электрического поля дифрагирующего света в точке (x, y, z) задается дифракционным уравнением Рэлея – Зоммерфельда первого типа [7, 8]:

, (1)

где E (0, y0, z0) – начальное поле в плоскости диафрагмы (x = 0),

– расстояние от точки в плоскости отверстия диафрагмы (0, y0, z0) до наблюдаемой точки (x, y, z),

i – мнимая единица,

λ – длина волны,

x – ось распространения световой волны,

k = 2π/λ – волновое число,

cosθ = x/r– направляющий косинус угла между x и r.

Рис. 2. Распределение светового потока при прохождении диафрагмы с круглым отверстием Примечание: составлен авторами по результатам данного исследования

Интегрирование ведется только по площади круглого отверстия S, так как поле вне отверстия равно нулю. Все геометрические размеры выражаются в микрометрах. Формула (1) соответствует интегралу Рэлея – Зоммерфельда первого рода в приближении kr≫1; по сравнению с точным осевым решением [9], вносимая этим погрешность осевой интенсивности не превышает 0,02 % при x ≥ 18 мкм и составляет около 0,5 % при x = 5 мкм.

На рис. 2 представлен случай распределения светового потока от излучающего торца ПОВ при прохождении через диафрагму с круглым отверстием, где θ – угол дифракционной расходимости, r – расстояние от точки в плоскости отверстия диафрагмы до наблюдаемой точки, dc – диаметр сердцевины ПОВ, dо – диаметр отверстия в диафрагме, x – ось распространения световой волны.

В модели принимается, что отверстие освещается пространственно когерентной плоской волной, падающей по нормали. Поскольку диаметр сердцевины ПОВ значительно превышает диаметр отверстия, поле в плоскости x = 0 при этом можно считать равномерным в пределах отверстия и равным нулю вне его (границы применимости допущения обсуждаются ниже):

(2)

Амплитуда E0 в модели полагается равной единице, так как в дальнейшем анализируется относительное распределение интенсивности I = |E|2.

Вне оптической оси интеграл Рэлея – Зоммерфельда в ближней зоне не выражается в замкнутом виде, поэтому необходимо построить численную схему [10], обеспечивающую компромисс между точностью и вычислительной эффективностью.

Реализация численной схемы

Область отверстия покрывается равномерной прямоугольной сеткой с шагом Δy0 = Δz0. Для адекватного воспроизведения фазовых соотношений шаг выбран из условия Δy0 ≤ λ/10; вопросы выбора шага дискретизации при численном интегрировании формулы Рэлея – Зоммерфельда рассмотрены в [6]. Шаг λ/2 соответствует пределу Найквиста для распространяющихся волн, пространственные частоты которых не превышают 1/λ. Однако при расстояниях порядка десятков длин волн фаза ядра eikr между соседними узлами меняется на величину до kΔy0sinθ, и для точного воспроизведения интерференционных членов в ядре интеграла требуется запас: при Δy0 = λ/10 этот набег не превышает 2 π/10 рад. Достаточность выбранного шага подтверждена исследованием сходимости (табл. 2).

Интеграл аппроксимируется двойной суммой по узлам сетки, попадающим в круглое отверстие:

(3)

где а индексы p, q нумеруют активные узлы.

Погрешность такой аппроксимации обусловлена заменой гладкой подынтегральной функции ее значениями в узлах сетки. Для функций с осциллирующим множителем eikr эта составляющая убывает с уменьшением шага по степенному закону; к ней добавляется погрешность ступенчатой аппроксимации границы круглого отверстия прямоугольной сеткой, которая и определяет фактический порядок сходимости; для уменьшения этой составляющей при криволинейной границе области применяют, например, треугольные сетки [11]. Погрешность оценивается методом повторного счета на последовательно сгущаемых сетках с экстраполяцией по правилу Рунге:

По данным табл. 2 разности решений на последовательных сетках уменьшаются примерно вдвое при уменьшении шага вдвое, то есть наблюдаемый порядок сходимости p ≈ 1, а не p = 2, характерный для квадратурной формулы при гладкой границе области), и именно это значение использовано в оценке. Для шага λ/10 погрешность дискретизации составляет около 0,4% по максимальной интенсивности и около 0,15% по диаметру пятна. Эта оценка справедлива в предположении асимптотического режима сходимости и сама по себе не является погрешностью относительно точного решения; близость к точному решению подтверждается сравнением с точным аналитическим решением на оси: при x = 30 мкм решение с шагом λ/10 отличается от него на 0,36 %, что согласуется с приведенной оценкой, а интегрирование в полярных координатах, не использующее ступенчатую аппроксимацию границы, воспроизводит его с расхождением менее 0,01 %.

Интенсивность вычисляется как

. (4)

Параметры численного эксперимента, использованные при моделировании, сведены в табл. 1.

Таблица 1

Параметры численного эксперимента

Примечание: составлена авторами на основе полученных данных в ходе исследования.

Оценка погрешности квадратурной формулы (3) дополнительно проведена сравнением с аналитическими решениями на оси и в дальней зоне и кросс-верификацией с независимым кодом; полученные расхождения не превышают 0,9 %.

Прямое вычисление двойной суммы для каждой точки наблюдения стандартными циклическими конструкциями приводит к неприемлемым временным затратам. Поэтому в основу программной реализации положена техника векторизации с использованием библиотеки NumPy языка Python [12].

Программный комплекс построен по модульному принципу и включает четыре функциональных блока: модуль исходных данных, задающий геометрические и волновые параметры; расчетное ядро, реализующее вычисление поля и интенсивности посредством векторизованных операций библиотеки NumPy, что заменяет циклы по узлам отверстия и по точкам строки наблюдения операциями над двумерными массивами и на порядки ускоряет расчет; модуль анализа, автоматически определяющий координаты максимума интенсивности, ширину пятна по уровню половинной интенсивности, угол дифракционной расходимости и число зон Френеля; модуль визуализации, формирующий двумерные карты интенсивности, радиальные сечения и картины распространения пучка вдоль оптической оси. Функциональная схема взаимодействия модулей комплекса представлена на рис. 3, фрагмент кода вычислительного ядра – на рис. 4.

Для подтверждения корректности разработанной численной схемы и ее программной реализации выполнена серия тестов, результаты которых сведены в табл. 3.

Рис. 3. Функциональная схема взаимодействия модулей программного комплекса Примечание: составлен авторами по результатам данного исследования

Рис. 4. Фрагмент кода расчетного ядра Примечание: составлен авторами по результатам данного исследования

Таблица 2

Исследование сходимости по шагу сетки (x = 30 мкм)

Примечание: составлена авторами на основе полученных данных в ходе исследования.

Таблица 3

Результаты верификации численной схемы

Примечание: составлена авторами на основе полученных данных в ходе исследования.

1. Сравнение с аналитическим решением: на больших расстояниях x ≫ d02 / λ распределение интенсивности должно соответствовать картине Эйри

(5)

где ρ – радиальная координата в плоскости наблюдения, J1 – функция Бесселя первого рода первого порядка. Поскольку диаметр центрального максимума при x = 200 мкм (около 52 мкм) не помещается в основное расчетное окно 40×40 мкм, для данного теста использовано отдельное окно 120×120 мкм с шагом 0,5 мкм в плоскости наблюдения; положение первого минимума уточнялось параболической интерполяцией по трем соседним узлам. Расхождение с теоретическим значением 2,44λx/d0 = 51,85 мкм не превысило 0,9% (табл. 3). Небольшое систематическое превышение объясняется тем, что при x = 200 мкм условие x ≫ d02 / λ ≈ 75 мкм выполняется лишь приближенно.

2. Кросс-верификация: осевая интенсивность и полуширина пятна в контрольных сечениях x = 30, 60 и 100 мкм, рассчитанные независимой реализацией схемы (3) в среде Matlab, совпали с результатами комплекса до четвертого значащего знака. Точность самой дискретизации дополнительно проверена интегрированием в полярных координатах (квадратура Гаусса – Лежандра, 300×720 узлов), не использующим ступенчатую аппроксимацию границы отверстия; расхождения не превышают 0,2 % (табл. 3).

3. Исследование сходимости: измельчение сетки от λ/2 до λ/40 (табл. 2) показало монотонную сходимость решения с достаточной точностью уже при шаге λ/10.

4. Сравнение с точным решением на оси: для ядра Рэлея – Зоммерфельда первого рода в точной форме осевое поле равномерно освещенного круглого отверстия выражается в замкнутом виде [9]:

Расхождение решения по схеме (3) с шагом λ/10 не превышает 0,53% (табл. 3), а положение осевого максимума (x = 18,8 мкм, n = 1) совпадает с результатом моделирования.

Видно, что при уменьшении шага с λ/10 до λ/20 максимальная интенсивность изменяется на 0,21 %, а диаметр пятна – на 0,07 %; при дальнейшем измельчении до λ/40 изменения уменьшаются примерно вдвое (0,11 и 0,04 %), что подтверждает сходимость схемы и достаточность шага λ/10 для практических расчетов. Интенсивность в таблице нормирована на значение, полученное с шагом λ/10. Сравнение результатов всех проведенных тестов приведено в табл. 3.

Для оценки устойчивости полученных результатов к неизбежным технологическим допускам проведен упрощенный анализ чувствительности. Варьировались основные параметры: диаметр отверстия d0 ±0,5 мкм, длина волны λ ±10 нм, боковое (поперечное) смещение растровой решетки и отводящего волокна ±2 мкм, а также толщина диафрагмы sд. Результаты представлены в табл. 4 (для x = 30 мкм; для ООВ – для x = 70 мкм, середины рекомендуемого диапазона его размещения).

Для бокового смещения растровой решетки в таблице приведена доля мощности пучка, проходящая через полосу пропускания шириной 8 мкм (при центрированной решетке – 0,85). Для бокового смещения ООВ приведена интенсивность в точке, соответствующей смещенному центру торца, относительно осевой в сечении x = 70 мкм; если диаметр сердцевины ООВ значительно превышает размер пятна, собираемый волокном поток при таком смещении практически не изменяется. В силу симметрии задачи результаты для смещений +2 и −2 мкм совпадают. Гауссов профиль освещения задан как

при той же амплитуде E0 в центре отверстия, поэтому снижение интенсивности отражает уменьшение потока, проходящего через отверстие.

Толщина диафрагмы (sд ≈ 5 мкм) учтена в предположении квазинормального падения излучения на отверстие введением дополнительного набега фазы в канале отверстия конечной длины:

δφ(ρ) = ksд

где θi(ρ) ≈ ρ/xкр – оценочный угол падения луча на точку отверстия на расстоянии ρ от оси. При этом допущении неоднородность набега фазы по апертуре при sд = 5 мкм не превышает 1,3·10−3 рад, а варьирование толщины в диапазоне sд = 4–6 мкм изменяет диаметр пятна и максимальную интенсивность менее чем на 0,1 % (табл. 4). Этот результат, как и остальные количественные результаты работы, относится к принятой модели пространственно когерентного освещения, близкого к нормальному. Излучение многомодового ПОВ содержит наклонные компоненты с углами до θNA и является частично когерентным: при заполненной апертуре волокна радиус когерентности в плоскости торца составляет порядка

≈ 2 мкм,

то есть меньше диаметра отверстия. В этом случае дифракционная картина формируется некогерентным наложением наклонных пучков и расширяется по сравнению с рассчитанной, а набег фазы в канале отверстия и частичное виньетирование стенками канала становятся заметными. Учет углового спектра и частичной когерентности излучения ПОВ вынесен в предмет дальнейших исследований.

С использованием созданного программного комплекса проведена серия вычислительных экспериментов для параметров λ = 850 нм, dо = 8 мкм, варьируемое расстояние x = 5–200 мкм. На рис. 5 показан результат моделирования (фрагмент x = 5–120 мкм) распределения интенсивности излучения при прохождении диафрагмы с круглым отверстием.

Таблица 4

Влияние вариаций параметров на диаметр дифракционного пятна и интенсивность

Примечание: составлена авторами на основе полученных данных в ходе исследования

Рис. 5. Распределение интенсивности излучения при прохождении диафрагмы с круглым отверстием Примечание: составлен авторами по результатам данного исследования

Световой поток после прохождения отверстия диафрагмы достаточно сколлимирован: полный угол дифракционной расходимости по уровню половинной интенсивности составил θ = 6,2° (угол полурасходимости θп = 3,1°), что согласуется с угловым размером центрального максимума картины Эйри по уровню 0,5: 1,03λ/d0 ≈ 6,3°. При x < 19 мкм (число зон Френеля n > 1) наблюдается чередование темных и светлых колец, включая темный центр при x ≈ 9 мкм (n ≈ 2) – эта зона непригодна для размещения решетки. Максимум осевой интенсивности достигается при x ≈ 18,8 мкм (n = 1). Далее формируется единое центральное пятно, интенсивность которого монотонно убывает с расстоянием. Растровую решетку целесообразно располагать в диапазоне x ≈ 20–30 мкм – сразу за областью n > 1, где интенсивность остается высокой, а кольцевая структура отсутствует. В качестве расчетного сечения далее принято x = 30 мкм.

Число зон Френеля рассчитывается как

, (6)

где dо – диаметр отверстия в диафрагме, λ – длина волны, x – расстояние от диафрагмы до растровой решетки.

При x = 30 мкм расчет показал, что количество зон Френеля составляет n = 0,627 (значению n = 1 соответствует x ≈ 18,8 мкм). При n < 1 отверстие не охватывает первую зону полностью, интерференция на его краях не дает существенного гашения или усиления, и наблюдается одно яркое пятно без кольцевой структуры [13].

а б

Рис. 6. Распределение интенсивности в дифракционном пятне после прохождения световым потоком диафрагмы с круглым отверстием: а – сечение на расстоянии x1 = 30 мкм, б – сечение на расстоянии x2 = 100 мкм Примечание: составлен авторами по результатам данного исследования

Для оценки размера дифракционного пятна и интенсивности излучения, попадающего на растровую решетку, были получены сечения на расстоянии x1 = 30 мкм и x2 = 100 мкм от диафрагмы, представленные на рис. 6.

В сечении а (x1 = 30 мкм) интенсивность на оси принята за единицу (табл. 1). Распределение неравномерно: диаметр центрального максимума по первым минимумам составляет около 8 мкм, при этом полуширина центральной области высокой интенсивности (уровень 0,5) составляет 3,35 мкм, за пределами которой интенсивность резко падает.

В сечении б (x₂ = 100 мкм) распределение интенсивности более однородно, резкие перепады отсутствуют, однако интенсивность примерно в 8 раз ниже, чем в сечении а, что затрудняет дальнейшее преобразование светового потока. Диаметр пятна по уровню половинной интенсивности составляет около 11 мкм (по первым минимумам – около 26 мкм).

Произведенный расчет, а также результаты моделирования показали, что диаметр дифракционного пятна по уровню 0,5 не превышает ширину зоны пропускания растровой решетки (8 мкм) вплоть до расстояния около 73 мкм; при этом с учетом падения интенсивности рационально располагать растровую решетку на расстоянии не более 60 мкм, так как при большем расстоянии снижается интенсивность излучения и увеличивается диаметр дифракционного пятна, что может привести к ложному срабатыванию счетчиков и внести неточность при дальнейшем расчете перемещения растровой решетки. Этот вывод подтверждается и энергетической оценкой: доля мощности пучка, попадающая в полосу шириной 8 мкм, уменьшается с 0,85 при x = 30 мкм до 0,77 при x = 60 мкм и 0,70 при x = 73 мкм, то есть при большем удалении решетки заметная часть потока попадает на соседние полосы. Для решеток с меньшей шириной полосы допустимое расстояние сокращается: например, при ширине полосы 4 мкм полуширина пятна достигает ее уже при x ≈ 36 мкм, а доля мощности в такой полосе при x = 30 мкм составляет 0,66. Учитывая, что интенсивность излучения в сечении, соответствующем расстоянию x2 = 100 мкм от диафрагмы, очень мала, а низкая интенсивность излучения наблюдается начиная с расстояния от диафрагмы x = 85 мкм, для выбранных параметров модели рациональным является диапазон размещения торца ООВ x = 60–80 мкм от диафрагмы [14]. При выборе положения торца ООВ следует также учитывать общие условия эффективного ввода излучения в волокно – согласование размера пятна и углового спектра излучения с сердцевиной и числовой апертурой волокна [15].

Заключение

Разработана и верифицирована математическая модель дифракции света на круглом отверстии в открытом канале волоконно-оптического растрового преобразователя, основанная на интеграле Рэлея – Зоммерфельда, и реализующий ее модульный программный комплекс на Python (численная схема прямого суммирования с векторизацией вычислений). Комплекс с контролируемой погрешностью рассчитывает распределения интенсивности в ближней зоне, углы расходимости пучка и рациональные зоны размещения элементов преобразователя. Показано, что при числе зон Френеля менее единицы (для принятых параметров – на расстояниях свыше 19 мкм) интерференционные кольца отсутствуют, однако дифракция по-прежнему определяет размер пятна и уровень сигнала. Учет дифракционных эффектов критически важен при проектировании высокоточных растровых преобразователей: он позволяет снизить вклад дифракционных искажений и обосновать выбор геометрии, что создает предпосылки для улучшения их метрологических характеристик. Разработанный инструментарий применим в системах автоматизированного проектирования волоконно-оптических измерительных устройств. Представленные результаты являются предварительными вычислительными рекомендациями; их экспериментальная проверка и учет дополнительных факторов (угловой спектр и частичная когерентность излучения многомодового оптического волокна, качество поверхностей, связь мод в волокне) составят предмет дальнейших исследований.