Особенности решения задачи рефракции на градиентных волоконных световодах различных типов Текст научной статьи по специальности «Физика»

ОСОБЕННОСТИ РЕШЕНИЯ ЗАДАЧИ РЕФРАКЦИИ НА ГРАДИЕНТНЫХ ВОЛОКОННЫХ СВЕТОВОДАХ РАЗЛИЧНЫХ ТИПОВ

© Андросик А.Б.*, Мировицкая С.Д.*

Московский государственный открытый университет, г. Москва

В работе рассмотрен новый подход к решению задачи моделирования основных геометрических и оптических характеристик градиентных волоконных световодов. Рассмотрена задача рефракции светового пучка на световодах с различными типами профилей показателя преломления.

1. Анализ различных типов градиентных волоконных световодов

Многообразие различных применений градиентных волоконных световодов (ВС) в приборостроении, волоконно-оптических линиях связи, в волоконно-оптических системах передачи (ВОСП), волоконных датчиках и др. привело к появлению большой группы профилей показателя преломления.

* Доцент кафедры Управления и информатики в технических системах, кандидат технических наук.

* Доцент кафедры Управления и информатики в технических системах, кандидат технических наук.

П

Рис. 2. Некруговое поперечное сечение световода

Рис. 3. Световод с двойной сердцевиной

Здесь обозначено п1, п2 - показатели преломления оболочки и поверхности сердцевины, г1, г2 - радиусы оболочки и сердцевины соответственно, р, g - показатели степени (рис. 1), ах, ау - размеры фигуры вдоль осей X, У -соответственно (рис. 2), Д - показатель профиля, ё - расстояние между осями световодов (рис. 3), Ь, с - константы.

Кривые зависимости:

= е = /М

Пп

для одиннадцати профилей показателя преломления ВС показаны на рис. 4.

Рис. 4. Основные типы профилей показателя преломления световодов

где

1 - обратно-квадратичный п = п€Е [1 — 2Д(1 —

2 - линейный п = n^-Jl — 2Д/х/.

3 - гиперболический-тангенсный п = псел/ 1 — 2&tk2(_x)\

4 - гауссов п = nCB [ 1 - 2А (l - е ~х*)];

5 - гиперболический секансный п = nCB5ch(x}:

6 - эллиптический п = пСЕ

7 - тангенсный п = псеЛ/! 1 — 2itgJ(ai);

8 - обобщенно-параболический п = пСЕ %/1 — 2Д.?"'; m = 2 — 2Д;

9 - неограниченно-степенной п — nCE^J 1 — 0 < Q < со;

10 - с ловушкой KARO п — пСЕ [1 — 2Д(—Ьх — сх2)]; 11-с ловушкой EXPO n = nCE[l- 2i(l - ebK+irar1)].

2. Исследование особенностей рефракции узкого пучка на градиентных ВС различных типов

Ниже для наиболее распространенных типов градиентных ВС определены границы изменения параметра Д и проанализированы зависимости угла рефракции от распределения показателя преломления по сечению ВС.

Эллиптический профиль показателя преломления. Обобщенно-эллиптический профиль показателя преломления:

где Д - относительная разность между n0 и показателем преломления оболочки световода; а- параметр, характеризующий профиль показателя преломления ВС.

Показатель преломления на оси ВС п0 = п(т^)/ (1 — 2 Д)1/2, тогда:

В процессе математического моделирования проанализированы параметры а и Д; при этом выбрано волокно с радиусом r0 = 100 мкм и показателем преломления на границе волокна n(r0) = 1,45. Угол рефракции определяется по формуле:

где rmi„ - решение уравнения rminn(rmiri~) = у (n.J.

Предположив, что Д и а = а - ЬД - ожидаемые параметры волокна, которым соответствует углы рефракции ф(у), можно установить отклонение реальных параметров от ожидаемых.

Допустимое отклонение параметра а от ожидаемого должно быть меньше 0,05, поэтому определяется отклонение а на значение 0,01. Вычисления показывают, что при изменении Д от 0 до 0,06, а а от 2 - 2Д - 0,05 до 2 - 2Д + 0,05, угол ф(у) будет изменяться так:

при у = 70, 80, 90 с увеличением а возрастает (фу), а с уменьшением а убывает ф(у);

при у = 70, 20, 30, 40 с увеличением а убывает ф(у), а с уменьшением а возрастает ф(у), кроме того, с увеличением Д всегда возрастает ф(у) (рис. 5, а -

а = 0,8, Ь - а = 0,9).

<55" 7,65 1,75 185 /г95 г/] 5

Рис. 5, а. Кривые зависимости угла рефракции от изменения параметра профиля; а = 0,8

Таким образом, чем больше у, тем более чувствительным к изменению а становится угол ф(у) (рис. 6). Поскольку при значениях у, близких к 100 мкм, возможно влияние краевых эффектов на результаты исследования, более подробно рассмотрена зависимость углов рефракции от параметров при

у = 80 мкм. При Д = 0,03 изменение а на 0,01 приводит к изменению (р(у) примерно на 0,05о, зависимость практически линейна при изменении а не более чем на 0,05 в ту или в другую сторону.

Рис. 5, б. Кривые зависимости угла рефракции от изменения параметра профиля; а = 0,9

При фиксированном значении Д = 0,03, у = 80 мкм и отклонении луча от расчетного пути не более, чем на 0,05о а отклоняется от ожидаемого значения не более чем на 0,01. Такого же порядка результаты получаются при других значениях Д.

Ниже рассмотрен случай, когда значение Д не совпадает с ожидаемым и отклонение параметра от ожидаемого не превосходит 10-3. Поскольку необходимо определить два параметра, выбирается два луча: один с координатой у = 80, а другой берется из интервала, где увеличение у приводит к уменьшению р(у), т.е. при у < 40. Расчеты показывают, что наиболее чувствительным к изменению а является луч с координатой у = 30. Отклонения параметров а, Д и угла р от ожидаемых значениях обозначены через уага, уагД и уагр. При |уага| < 0,05, |уагД| < 0,001, |уагр < 0,05 зависит от уага и уагД линейно, т.е. уагр(80) = к1, уага + Ь1, уагД; уагр(30) = -к2, уага + Ь2, уагД, где коэффициенты к1, к2, Ь1, Ь2 определяются методом неопределенных коэффициентов.

Рис. 6. Кривые зависимости <y) в функции а

Перед k2 во втором выражении системы стоит знак минус, так как при y = 30, vara< 0 и <р< 0.

Решив систему уравнений, определяются vara и varA через var<(80) и var<(30).

Пусть ожидаемое значение A = 0,03; а = 2 - 2; A = 1,94. Тогда при y = 80 мкм и изменении а на 0,01 < меняется примерно на 0,05о, а при изменении A на 10-3< меняется примерно на 0,11о. Если же y = 30, то при изменении а на 0,01 угол < меняется на 0,003о, а при изменении A на 10-3 угол < меняется примерно на 0,048о. В этом случае k = 0,5; k2 = 0,3; b = 110; b2 = 48, и получается система уравнений вида: var<(80) = 0,5, vara + 110, varA; var<(30) = -0,3, vara + 48, varA. Решая систему уравнений, определяется vara = [48 var<(80) -- 110var<(30)] / 57.

Обобщенный параболический профиль. Закон распределения показателя преломления: «(г) = п„(1— 2Д(г/г;3)Е) при 0 < г < г0. Расчеты проводились для п0 = 1,53; г9 = 100 мкм; а = 2 — (12/5}Д.

Гауссовский профиль. Закон распределения показателя преломления: >:;:■; = : ■ е 1 Ji :.:■/:; :"' Исходные данные для расчета выбирались аналогично предыдущему профилю. Результаты расчета представлены в табл. 1 с точностью не менее 10~4. Поскольку n(r0) > 1, Д< (lnn0)/2 = 0.228712-,

Обобщенный эллиптический профиль с провалом на оси. Закон распределения показателя преломления:

Вычисления проводились для реального профиля с п0 = 1,58; г0 = 100 мкм, w = 0,05г0, а = 2 - 2Д, пс = 1,50 с точностью не менее 10-4.

Графики зависимости ср(у) для обобщенно-эллиптического профиля (линия верхняя Д = 0,1) и профиля обобщенно-эллиптического с провалом на оси (линия нижняя Д = 0,1) показаны на рис. 7, а.

Рис. 7, а. Результаты расчета для обобщенно-эллиптического профиля

Штриховая линия соответствует Д = 0,0135 для обобщенно-эллиптического профиля. Результирующие зависимости п = п(г) для трех типов профилей: для обобщенных эллиптического (линия 1 соответствует Д = 0,1; 2 - Д = 0,25), параболического (линия 3 - Д = 0,1; линия 4 - Д = 0,15) и эллиптического с провалом на оси (линия 5 - Д = 0,1; линия 6 - Д = 0,25) также представлены на рис 7, б. При г > 15 графики для обобщенно-эллиптического профиля с провалом на оси при Д = 0,1 и Д = 0,25 совпадают с графиком для обобщенно-эллиптического профиля.

Исследования полученных зависимостей позволяют определить границы изменения параметра профиля Д:

1. для обобщенно-эллиптического профиля из условия /7(г0) > 1, т.е. VI - 2Д> 1 следует Д< (1 - 1 /и^)/2. При щ = 1,58, Д = 0,2997116;

2. для обобщенно-параболического профиля из условия п(г0) > 1, т.е. VI - 2Д> 1, следует, что при и0(1 - 2Д) > 1 и п0 = 1,58, А < (1 — 1 / п0) / / 2 = 1,1835;

3. для обобщенно-эллиптического профиля с провалом на оси при п) = 0,1. ¡'и = 5 и г > 20; ехр(—г:/м/г) < Ю-3 показатель преломления для случая обобщенно-эллиптического профиля, поэтому значение Д здесь те же, что и в первом случае.

Рис. 7, б. Результаты расчета для трех типов профилей

3. Алгоритм вычисления угла рефракции при зондировании градиентных ВС

Как показали исследования, при зондировании градиентных ВС угол отклонения луча записывается в виде интеграла Абеля, описывающего зависимость угла рефракции от показателя преломления исследуемой осе-симметричной неоднородности. Такой интеграл имеет особенность в одном из пределов интегрирования, и точность его вычисления зависит от способа разбиения интервала интегрирования.

В случае ВС с обобщенно-эллиптическим профилем показателя преломления угол рефракции описывается квадратурной формулой Абеля:

где у - координата точки входа зондирующего пучка в ВС;

пъ, п0 - значения показателя преломления на границе и на оси ВС,

соответственно; Ъ - радиус световода; т - регулируемый параметр;

Д - параметр, характеризующий относительную разность показателей преломления на оси и на поверхности световода; г0 - параметр, вычисляемый из уравнения:

У«Ь =г0п0[1-2&(г0/Ь)т]и2 Кроме того, введены следующие обозначения:

А== (п1 - п1)(2п20) = [ 1 - (пь/п0)г]/2 откуда получается:

Произведя подстановки в формуле (1), г = х/Ь, гд — гь/Ъ. д = у/Ъ, можно записать:

(2)

¡агссоя^г,

Для упрощения вводится следующее обозначение

Анализ этого интеграла показывает, что 2 зависит от отношения расстояния у от оптической оси до точки входа луча в неоднородность к радиусу Ъ оптического волокна. При сохранении отношения у / Ъ величина 2 не изменяется с уменьшением радиуса. Кроме того, как видно из интегральной формулы (1), значение 2 зависит от параметра Д, характеризующего разность показателей преломления на границе пъ и на оси п0 волокна.

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

Для вычисления Q интервал [г0, 1], разбивается на два подинтервала [¿о, О], [О, 1] в отношении 1 : W (рис. 8), т.е.:

Рис. 8. Пример деления интервала на подинтервалы и кривая поведения ошибки интегрирования при изменении числа разбиений

На подинтервале (Ц 1) задается число разбиений К и вычисляется частичный интеграл, являющийся первой составляющей частью искомого интеграла Q. Затем аналогично осуществляется разбиение подинтервала (10, О) в отношении 1 : Ж На правой его составляющей вычисляется следующий частичный интеграл I и добавляется к Q и т.д. Прерывание этого процесса осуществляется при выполнении условия по точности I / Q < Е.

Обобщенная формула Симпсона имеет вид:

где к = (Ь - а) / 2к, к - коэффициент Симпсона.

Выражение для остаточного члена показывает, что формула Симпсона точна, если Ах) - многочлен до третьей степени. Эта частная особенность формулы Симпсона объясняет ее преимущественное применение.

На фиксированном интервале интегрирования методом Симпсона ошибка ст интегрирования при изменении числа разбиений интервала интегрирования имеет вид кривой с минимумом (рис. 8), зависящим от вида подинте-гральной функции / Подобный вид кривой объясняется тем, что при малых значениях К не обеспечивается достаточная точность аппроксимации полиномом, а при больших К начинает сказываться ошибка, вызванная увеличением числа вычислительных операций.

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

(3)

Истинное значение этой величины известно и равно л/3; с ним сравнивалось значение вычисленного интеграла Q.

Точность вычислений определяется точностью вычисления первой составляющей интеграла правой части (выражения (2). Вычисляя значение угла рефракции по выражению (2) при а = 0, можно получить значение абсолютной ошибки выражений. Такой подход не дает объективной оценки точности вычислений и позволяет оценить только методическую ошибку и ее зависимость от коэффициента Симпсона к и коэффициента разбиения на подинтервалы Ж, а также требуемой точности частичного интегрирования Е, т.е. с изменением параметра Д форма кривой подинтегральной функции изменяется (в частности, изменяется ее крутизна).

Для оценки точности вычислений проводится интегрирование коэффициентов к и V (используя значение показателя профиля Д, не равное нулю).

Алгоритм вычисления угла рефракции в зависимости от изменения параметров ВС, обеспечивает ввод следующих параметров: E - требуемая точность частичного интегрирования, W - коэффициент разбиения на по-динтервалы, k - коэффициент Симпсона, r - радиус волокна, ymin, ymax, hy -диапазон изменения значения показателя профиля и шаг разбиения внутри диапазона.

После ввода указанных выше параметров и ответа на запрос о виде вывода, выполнятся подпрограмма вывода таблицы и графика (рис. 9).

Алгоритм подпрограммы вывода таблицы позволяет выводить таблицу значений функций (в данном случае функции urefrj) двух переменных. В циклах for осуществляется изменение аргументов функции от их минимального до максимального значения (в данном случае а и y) с заданным шагом и вывод значения функции.

ЛАГОтТМ ВЫЧИСЛЕН«!

Рис. 9. Алгоритмы вычисления и исследования точности

Блок-схема подпрограммы вывода графика позволяет ввести семейство графиков зависимости функции двух переменных (в данном случае urefq = =fa, y). Подпрограмма запрашивает цену деления по оси X, после ее ввода формируется массив значений функций при изменения аргумента, отложенного по оси X от min до max с шагом, равным введенной цене деления. При формирования массива определяется min и max значений функции, после чего запрашивается min, max значений и шаг по оцифровке оси Y. После нормализации оси Y осуществляется инициализация графического режима, определяются границы поля графика и обеспечивается автоматическое по-

строение осей координат и выбор начала координат в зависимости от знаков аргумента и функции. График функции строится по сформированному массиву значений функции. По завершению построения автоматически в конце каждой кривой печатается значение второго аргумента.

Алгоритм подпрограммы вычисления угла рефракции (рис. 10). Сначала вычисляется переменная V, затем в цикле а0 реализуется метод простых итераций, позволяющий вычислить значение г0. Далее обнуляется значение искомого интеграла g, вычисляется интервал интегрирования ^ в цикле «о интервал t разбивается на подинтервалы в соотношении 1 : Ж, определяются пределы интегрирования частичного интеграла а и Ь, в подпрограмме вычисляется частичный интеграл I и добавляется к g. При получении из подпрограммы значение ]г = 1, программа выходит из цикла. В этом случае, а также при нормальном окончании цикла функции ыгв/г возвращается значение угла рефракции, вычисленное согласно формуле (2).

-2агсс<и|ул1

ZZI_

2 ACOSty/E) ( »ЕТМШ)

Рис. 10. Алгоритм вычисления угла рефракции

Алгоритм подпрограммы вычисления частичного интеграла составлен на основании выражения (3) и реализует интегрирование методом Симпсо-на (рис. 11). Интервал интегрирования разбивается на 2k частей, обнуляется значение частичного интеграла i и в цикле for обеспечивается приращение X с шагом h, вычисление в подпрограмме значения подинтегральной функции и в зависимости от шага цикла добавления к искомому частичному интегралу i значения подинтегральной функции (при x = a и x = b), значения подинтегральной функции, умноженное на 2 при четных шагах умноженное на 4 при нечетных. Для определения четности шага применяется операция

«остаток». При получении из подпрограммы значения jr = 1 происходит выход их цикла. В этом случае и при нормальном окончании цикла функции intpr присваивается значение i *h/3.

Алгоритм вычисления подинтегральной функции. Переменной f присваивается значение подкоренного выражения (2), если это значение отрицательное или равно нулю, переменной jr присваивается значение 1, в противном случае вычисляется значение подинтегральной функции и присваивается к функции fun.

Рис. 11. Алгоритм вычисления интеграла

Алгоритм исследования точности вычисления (рис. 9) обеспечивает получение зависимости модуля абсолютной ошибки вычисления от параметров, определяющих точность: Е, Ж, Я, т.е обеспечивается ввод параметров Е, г, у, Ктш, Ктах, п, Я - диапазона изменения коэффициента Симпсона и шага по разбиению внутри поддиапазона; Жтш, Жтах, к, w - диапазона изменения коэффициента разбиения на подинтервалы и шага по разбиению внутри диапазона; значение показателя профиля

При вводе параметра & = 0 осуществляется переход к подпрорамме вывода графика, при этом выводится зависимость модуля абсолютной ошибки от параметров Е, Ж, к. При вводе Ло осуществляется переход к подпрограмме вывода таблицы и распечатываются значения углов рефракции для вариации параметров Ж и к, оценка точности при этом осуществляется по числу первых повторяющих цифр.

Анализ кривых (рис. 12) показывает, что с увеличением коэффициента К ошибка вычислений ст сначала уменьшается, при определенном значении достигает минимума и начинает возрастать, причем минимум имеет ярко выраженный характер. С увеличением коэффициента Ж минимум ошибки сдвигается в сторону увеличения коэффициента к (по оси У на графиках отложены значения порядка модуля абсолютной ошибки, т.е. 1§|ст|).

Результаты моделирования представлены на рис. 12.

Е«1Е-08 Я»ЮО <!•() У»50

20 22 24 2К 28 30 32 34 36 К

Рис. 12. График зависимости модуля абсолютной ошибки от моделируемых параметров

Анализ кривых (рис. 12) показывает, что с увеличением коэффициента К ошибка вычислений ст сначала уменьшается, при определенном значении достигает минимума и начинает возрастать, причем минимум имеет ярко выраженный характер. С увеличением коэффициента Ж минимум ошибки сдвигается в сторону увеличения коэффициента к (по оси У на графиках отложены значения порядка модуля абсолютной ошибки, т.е. 1§|ст|).

Список литературы:

1. Лазарев Л.П., Мировицкая С.Д. Контроль геометрических и оптических параметров волокон. - М.: Радио и связь, 1988. - 280 с.

2. Андросик А.Б., Воробьев С.А., Мировицкая С.Д. Анализ основных типов оптических волокон // Приволжский научный вестник. - 2011. - № 4. -С. 7-19.

3. Андросик А.Б., Мировицкая С.Д. Анализ особенностей метода расе-яния в переднюю полусферу на двухслойных световодах // Актуальные вопросы современной науки. Вып. 21, ч. 2. - Новосибирск. - С. 29-44.