Содержание к диссертации
Введение
ГЛАВА 1. Проблемы моделирования неизотермической фильтрации в областях со сложной геометрией - 9
1.1 Обзор математических моделей применяемых для процессов неизотермической фильтрации..- 9
1.2 Существующие численные методы для решения задач неизотермической фильтрации в областях со сложной геометрией - 15
1.3 Проблемы применения коммерческого программного обеспечения для задач пеизотермической фильтрации в областях со сложной геометрией - 21
1.4 Выводы -24 -
ГЛАВА 2. Моделирование однофазной неизотермической фильтрации жидкости или газа на основе мкэко - 26
2.1 Постановка задачи -26 -
2.2 Метод дискретизации расчетной области - 29 -
2.3 Особенности МКЭКО с равным порядком интерполяции скорости и давления - 35 -
2.4 Вывод дискретного аналога для уравнения неразрывности и движения - 37 -
2.5 Особенности дискретного аналога уравнения переноса и энергии - 45 -
2.6 Алгоритм численного решения - 46 -
2.7 Программная реализация метода - 47 -
2.8 Тестирование разработанного программного обеспечения - 49 -
2.9 Выводы -73 -
ГЛАВА 3. Моделирование двухфазной неизотермической фильтрации на основе МКЭКО - 74 -
3.1 Постановка задачи -74 -
3.2 Вывод дискретного аналога уравнения для определения давления - 76 -
33 Вывод дискретного аналога уравнения неразрывности нефтяной фазы - 81 -
3.4 Особенности дискретного аналога уравнения энергии - 83 -
3.5 Программная реализация метода - 87 -
3.6 Тестирование разработанного программного обеспечения - 90 -
3.7 Выводы - 100-
ГЛАВА 4. Моделирование неизотермической двухфазной фильтрации смешивающихся слабосжимаемых жидкостей -101 -
4.1 Постановка задачи и особенности математической модели - 101 -
4.2 Метод дискретизации расчетной области - 106 -
4.3 Гибридный МКЭКО с равным порядком интерполяции скорости и давления - 109 -
4.4 Вывод дискретного аналога уравнения неразрывности - 110-
4.5 Особенности дискретного аналога конвективно - диффузионного уравнения - 112 -
4.6 Алгоритм численного решения - 114 -
4.7 Тестирование метода - 115-
4.8 Выводы - 125-
ГЛАВА 5. Исследование влияния неизотермичности в прискважиннои области на процесс пробоотбора в изотропных и анизотропных пластах -126-
5.1 Постановка задачи - 126 -
5.2 Анализ результатов решения изотермической задачи пробоотбора для вертикальной скважины .-136-
5.3 Анализ результатов решения неизотермическои задачи пробоотбора для вертикальной скважины -138 -
5.4 Особенности постановки задачи для случая пробоотбора через горизонтальную скважину...- 141 -
5.5 Влияние толщины пласта на процесс пробоотбора - 144 -
5.6 Влияние неизотсрмичностн на процесс пробоотбора - 145 -
5.7 Влияние анизотропии плста на процесс пробоотбора - 147 -
5.8 Выводы - 149-
Заключение -150-
Литература -151 -
- Существующие численные методы для решения задач неизотермической фильтрации в областях со сложной геометрией
- Метод дискретизации расчетной области
- Вывод дискретного аналога уравнения неразрывности нефтяной фазы
- Гибридный МКЭКО с равным порядком интерполяции скорости и давления
Введение к работе
Актуальность проблемы.
Моделирование разработки нефтегазоносных пластов осложняется значительной неоднородностью и анизотропией свойств пород, слагающих продуктивный пласт коллектора, а также тем, что при тепловой обработке коллектора среды, содержащиеся в пластах, приобретают температуру, отличную от естественной пластовой температуры. Кроме того, значительное влияние на продуктивность пласта оказывает околоскважинная зона, свойства которой обычно значительно изменены в сравнении с остальной частью пласта. Таким образом, фильтрационные процессы определяются различными характерными размерами: размер скважин, околоскважинных зон, размер контура питания. Поэтому разработка специальных численных методов, позволяющих решать разномасштабные неизотермические задачи фильтрации, является актуальной.
Цель работы.
Цель работы заключается в исследовании процессов неизотермической одно-и двухфазной фильтрации с учетом геометрических особенностей пласта, околоскважинной зоны и при наличии анизотропии теплопроводности и проницаемости на основе разработки эффективного численного метода
Научная новизна работы:
Разработаны консервативные варианты конечно - элементного метода контрольного объёма для численного решения задач неизотермической однофазной фильтрации жидкости или газа и двухфазной фильтрации несмешивающихся слабосжимаемых жидкостей в двумерных областях и двухфазной смешивающейся фильтрации в трехмерных областях.
Разработан новый способ учета тензора проницаемости и эффективной теплопроводности для анизотропных резервуаров.
Разработан новый метод "неявное давление - неявная насыщенность" для задач смешивающийся двухфазной фильтрации.
Установлено, что неизотермичность оказывает значительное влияние на очистку пробы на начальной стадии процесса отбора пластовых флюидов и не влияет на степень загрязнения пробы на конечной стадии пробоотбора как для горизонтальной, так и для вертикальной скважин.
Научная и практическая значимость работы.
Научная и практическая значимость работы заключается в разработке консервативных вариантов конечно-элементного метода контрольного объёма, позволяющих моделировать разномасштабные задачи одно- и двухфазной фильтрации с учетом геометрических особенностей, при наличии неоднородности свойств пласта и анизотропии проницаемости и теплопроводности, используя грубые сетки с локальным неструктурированным сгущением. Предложенный для случая смешивающейся фильтрации слабосжимаемых жидкостей численный метод решения многофазных задач "неявное давление - неявная насыщенность" позволяет использовать неявную разностную схему, которая является безусловно устойчивой.
Практическая ценность заключается в создании программных комплексов на языке Fortran для численного моделирования одно- и двухфазной фильтрации. Решены прикладные задачи о процессе пробоотбора в вертикальной и горизонтальной скважинах.
Апробация работы.
Основные результаты диссертационной работы докладывались и обсуждались на российских и международных научно-технических конференциях и семинарах:
Семинар лаборатории радиационной газовой динамики Института проблем механики им. А.Ю. Ишлинского РАН (Москва, 2009)
Заседание кафедры геофизики физического факультета Башкирского государственного университета (Уфа, 2009)
Вторая всероссийская конференция молодых ученых и специалистов «Будущее машиностроение России» (Москва, 2009)
IV-я международная научно-практическая конференция STAR-2009: «Компьютерные технологии решения прикладных задач тепломассопереноса и прочности» (Нижний Новгород, 2009)
XV международная научно-техническая конференция студентов и аспирантов «Радиоэлектроника, электротехника и энергетика» (Москва, 2009)
Российская техническая нефтегазовая конференция и выставка SPE 2008 (Москва, 2008)
X международная конференция «Тепловое поле Земли и методы его изучения» (Москва, 2008)
Научно-практическая конференция «Математическое моделирование и компьютерные технологии в разработке месторождений» (Уфа, 2008)
9. Выставка научно-технического творчества студентов в рамках XII Всемирного Русского Народного Собора (Москва, 2008)
10.Вторая студенческая научно-инженерная выставка «Политехника» (Москва, 2007)
Публикации.
Основные результаты диссертации опубликованы в 6 печатных работах, список которых приводится в конце автореферата. По результатам работы подано 2 заявки на патент.
Структура и объем диссертации.
Диссертация состоит из введения, пяти глав, заключения и списка цитируемой литературы из 134 наименований. Общий объем работы составляет 165 страницы, включающих 83 рисунка и 15 таблиц.
Существующие численные методы для решения задач неизотермической фильтрации в областях со сложной геометрией
Процессы тепло- и массообмена в геологических пластах осложняются особенностями геометрии, такими как многопластовые коллекторы с трещинами и разрывами или прискважинные зоны. В этих случаях возникает необходимость формирования расчетных сеток, обладающих более гибкими и разнообразными свойствами с точки зрения возможности адаптации сетки к геометрии задачи и особенностям решения. Для этих целей применяют сетки с нерегулярной структурой. Среди нерегулярных сеток выделим два основных класса: 1. Структурированные сетки. Наиболее известным примером таких сеток являются нерегулярные четырехугольные сетки (Рис. 1.1), которые во многом наследуют свойства стандартных прямоугольных сеток. 2. Неструктурированные сетки. В качестве примера можно рассмотреть треугольные сетки. Для них шаблон разностной схемы не сохраняет структуру. В отличие от случая структурированной сетки, в этом случае для каждого узла расчетной области строится разностный шаблон с учетом расположения всех его соседей. Преимущества структурированных сеток связаны с сохранением канонической структуры сеточного шаблона для каждого узла сетки. Данный подход значительно упрощает построение дискретного аналога и программную реализацию. Среди структурированных сеток необходимо выделить важный класс ортогональных сеток. В этом случае очень хорошо проявляются преимущества структурированных сеток перед неструктурированными, такие как простота построения дискретных аналогов и быстрое решение системы линейных уравнений. Аппроксимация на ортогональных структурированных сетках может проводиться подобно аппроксимации на стандартных прямоугольных сетках. Для криволинейных неортогональных сеток, при применении новых независимых переменных [64] происходит преобразование координат. В этом случае нерегулярная в исходных координатах сетка преобразуется в регулярную в новых независимых координатах (Рис. 1.1а).
Треугольные сетки являются наиболее используемыми неструктурированными сетками. Использовать более сложные топологические типы [31] неструктурированных сеток чаще всего нет необходимости. При построении треугольных сеток (триангуляции) следует обращать внимание на два основных требования. Первое требование триангуляции состоит в том, чтобы полученные треугольники были близки к равносторонним. Вторым требованием является равномерность сетки. Этим требованиям соответствует триангуляция по Делоне, которая обладает рядом оптимальных свойств: при триангуляции по Делоне максимизируется минимальное значение внутренних углов треугольников [55]. Триангуляция по Делоне активно используется в вычислительной практике при построении схем конечных элементов. Имеется также много хорошо проработанных вычислительных методов генерации таких треугольных сеток [69, 93, 50]. Отметим некоторые общие подходы к построению разностных схем. на нерегулярных сетках: метод опорных операторов, метод конечных элементов, метод контрольного объёма. В случае метода опорных операторов [29] задача формулируется в терминах инвариантных операторов векторного анализа - дивергенция, градиент.
Далее один из них (который и называется опорным) аппроксимируется, а аппроксимации других согласуются с аппроксимацией опорного на основе интегральных соотношений векторного анализа. В идейном плане этот подход примыкает к смешанному методу конечных элементов [84, 78], когда исходное эллиптическое уравнение второго порядка записывается в виде системы уравнений первого порядка и в качестве неизвестных выступает само решение и его первые производные. Как и для общих неструктурированных сеток, так и для структурированных можно строить разностные схемы на основе конечно-элементных аппроксимаций [26]. В работах [75, 132, 104, 108, 127, 134, 125] для моделирования одно- и двухфазных течений в пористой среде авторы
Метод дискретизации расчетной области
Рассматриваемый вариант конечно-элементного метода контрольного объёма [65] требует разбиения расчетной области на множество трехузловых непересекающихся треугольных подобластей, называемых конечными элементами (Рис. 2.1). особенности: 1. Сложные области могут быть более эффективно разбиты на треугольные элементы. 2. Не требуют при вычислении интегралов использования изопараметрических преобразований, в отличие от четырехугольных элементов.
Формирование контрольных объемов около узлов i-J—k, находящихся в вершинах конечных элементов происходит следующим образом: 1) В конечном элементе определяется положение центра масс, который имеет следующие координаты где (xt, уі), (xj, yf), (xk, yk) - координаты вершин элемента в глобальной системе координат (х, у). Затем центр масс элемента, который соответствует точке 2 на Рис. 2.2а, соединяется с серединами всех сторон точками 1, 3, 4. Часть области, ограниченная ломаной і—1—2—3-і, является вкладом УЦ элемента i—j-k в контрольный объём /-го узла от 7-го элемента, часть области j—4—2—l—j является вкладом V:\ в у — го узла от 1-го элемента и часть к-4-2—3-к есть вклад Vu элемента в контрольный объем к — го узла от 7-го элемента. 2) Объединяем вклады от каждого элемента, который содержит данный узел, в одну область, которая будет называться контрольным объемом для данного узла. Контрольный объем V-, для узла і показан на Рис. 2.26. Многоугольные контрольные объёмы, аналогичные представленному на Рис. 2.26 обладают следующими полезными свойствами: они не перекрываются, полностью заполняя расчетную область; их границы не включают межэлементных границ, что позволяет использовать разные свойства (материалы) на элементах; -31 -они могут быть построены на сетке, включающей тупоугольные треугольники.
В декартовой системе координаты объемы вкладов конечного элемента в соответствующий контрольный объем следующие (Рис. 2.3) где А - площадь треугольника i-j-k. Размер расчетной области для двухмерного случая в направлении оси z составляет Az= 1м. Rgi Rgk R где Rg;- радиус центра масс і-1-2-3, R, - радиус узла /, R - радиус центра масс j-4-2-1, Ry - радиус узла , R - радиус центра масс к-3-2-4, Rk - радиус узла к, Аср = 1 радиан - угловой размер расчетной области. -На элементе i-j-k могут применяться следующие два вида интерполяции: линейная и ступенчатая. 1. Линейная интерполяция (Рис. 2.4) треугольного линейного элемента функции формы имеют следующий вид, представленный в работе [26]. 2. Ступенчатая интерполяция (Рис. 2.5) Если на элементе используется ступенчатая интерполяция, то считается, что значение функции Ф постоянно по всему контрольному объему для данного узла и равно среднеинтегральному значению функции в узле. В этом случае на элементе функция Ф определяется выражением -Разработанный метод может применяться при наличии анизотропии свойств пласта. Рассмотрим конечный элемент, главные оси проницаемости (х ,У) и теплопроводности (х", ") которого развернуты относительно глобальных осей координат (х,у) на углы х и # и сохраняют свое направление в любой точке элемента. Для удобства будем использовать локальные системы координат [X ,Y ), (X",Y") оси которых сонаправлены с главными осями проницаемости и теплопроводности элемента, а начало координат перенесено из центра масс элемента в начало глобальной системы координат (Рис. 2.6).
Вывод дискретного аналога уравнения неразрывности нефтяной фазы
Рассматриваемый вариант конечно-элементного метода контрольного объёма [65] требует разбиения расчетной области на множество трехузловых непересекающихся треугольных подобластей, называемых конечными элементами (Рис. 2.1). особенности: 1. Сложные области могут быть более эффективно разбиты на треугольные элементы. 2. Не требуют при вычислении интегралов использования изопараметрических преобразований, в отличие от четырехугольных элементов. Формирование контрольных объемов около узлов i-J—k, находящихся в вершинах конечных элементов происходит следующим образом: 1) В конечном элементе определяется положение центра масс, который имеет следующие координаты где (xt, уі), (xj, yf), (xk, yk) - координаты вершин элемента в глобальной системе координат (х, у). Затем центр масс элемента, который соответствует точке 2 на Рис. 2.2а, соединяется с серединами всех сторон точками 1, 3, 4. Часть области, ограниченная ломаной і—1—2—3-і, является вкладом УЦ элемента i—j-k в контрольный объём /-го узла от 7-го элемента, часть области j—4—2—l—j является вкладом V:\ в у — го узла от 1-го элемента и часть к-4-2—3-к есть вклад Vu элемента в контрольный объем к — го узла от 7-го элемента. 2) Объединяем вклады от каждого элемента, который содержит данный узел, в одну область, которая будет называться контрольным объемом для данного узла. Контрольный объем V-, для узла і показан на Рис. 2.26. Многоугольные контрольные объёмы, аналогичные представленному на Рис. 2.26 обладают следующими полезными свойствами: они не перекрываются, полностью заполняя расчетную область; их границы не включают межэлементных границ, что позволяет использовать разные свойства (материалы) на элементах; -31 -они могут быть построены на сетке, включающей тупоугольные треугольники.
В декартовой системе координаты объемы вкладов конечного элемента в соответствующий контрольный объем следующие (Рис. 2.3) где А - площадь треугольника i-j-k. Размер расчетной области для двухмерного случая в направлении оси z составляет Az= 1м. Rgi Rgk R где Rg;- радиус центра масс і-1-2-3, R, - радиус узла /, R - радиус центра масс j-4-2-1, Ry - радиус узла , R - радиус центра масс к-3-2-4, Rk - радиус узла к, Аср = 1 радиан - угловой размер расчетной области. -На элементе i-j-k могут применяться следующие два вида интерполяции: линейная и ступенчатая. 1. Линейная интерполяция (Рис. 2.4) треугольного линейного элемента функции формы имеют следующий вид, представленный в работе [26]. 2. Ступенчатая интерполяция (Рис. 2.5) Если на элементе используется ступенчатая интерполяция, то считается, что значение функции Ф постоянно по всему контрольному объему для данного узла и равно среднеинтегральному значению функции в узле. В этом случае на элементе функция Ф определяется выражением -Разработанный метод может применяться при наличии анизотропии свойств пласта. Рассмотрим конечный элемент, главные оси проницаемости (х ,У) и теплопроводности (х", ") которого развернуты относительно глобальных осей координат (х,у) на углы х и # и сохраняют свое направление в любой точке элемента. Для удобства будем использовать локальные системы координат [X ,Y ), (X",Y") оси которых сонаправлены с главными осями проницаемости и теплопроводности элемента, а начало координат перенесено из центра масс элемента в начало глобальной системы координат (Рис. 2.6). kx, Трудности метода равного порядка интерполяции связаны с тем, что если даже на равномерной сетке давление определяется с использованием линейной функции формы, то только разность давлений между чередующимися (каждым вторым) узлами включается в полную систему уравнений. Следовательно, разностные уравнения не могут обнаружить разницу между равномерным и шахматообразным полями давлений. В конечно-разностных формулировках МКО для решения уравнений Навье-Стокса предотвращение шахматного поля давлений осуществляется за счет использования смещенных сеток [44]. В конечно-элементных методах шахматное поле давления исключается посредством: применением разных сеток конечных элементов для определения скорости и давления [99] (Рис. 2.7); получения давления с использованием уравнения
Пуассона, как это сделано Шнейдером и др. [126]; использования некоторых специальных методик для фильтрации нереальных колебаний давления [123]; использования формулировок через штрафные функции [92]. В данной работе рассматривается модификация метода равного порядка интерполяции скорости и давления [118], позволяющего избежать шахматообразного поля давления для задач фильтрации. Основная идея метода заключается в том, что массовый расход, рассчитанный по скоростям w, определяемым в узлах, заменяется расходом, рассчитанным по скоростям w, ассоциированным с / —ым конечным элементом. Эта идея схожа с подходом использования смещенных сеток. Поле скорости w в локальной системе координат определяется на каждом элементе и имеет разрывы в его узлах (Рис. 2.8) где ех,,ег,- единичные орты системы координат yX ,Yj связанной с главными осями проницаемости на элементе. Компоненты поля скорости w на элементе / определяются следующим образом элементе i-j-k (1.20). Значения коэффициентов D" и DJ определяются в соответствии с (1.21). 2.4 вывод дискретного аналога для уравнения неразрывности и движения Перепишем уравнение неразрывности (1.11) в виде Применим локальный закон сохранения массы к каждому контрольному объему. Тогда для расчетной области, состоящей из М контрольных объемов, получается следующая система уравнений dt В то же время расчетная область Q состоит из N конечных элементов. Таким образом, можно записать следующее выражение для уравнения неразрывности где J 7 - вклад от / - элемента в контрольный объём / - го узла, Р\, - вклад от / -элемента в контрольный объём j - го узла, Vkl - вклад от / - элемента в контрольный объём к - го узла. Рассмотрим интеграл по вкладу V., в (1.22) Выражение (1-23) включает нестационарный ф др_ dt и дивергентныйdiv(p\v) члены. Суммируя их, получим дискретный аналог уравнения. Используя теорему Остроградского - Гаусса, получим
Гибридный МКЭКО с равным порядком интерполяции скорости и давления
В предложенном трехмерном гибридном МКЭКО, также как и в двумерном МКЭКО, используется метод равного порядка интерполяции скорости и давления, который заключается в том, что массовый расход, рассчитанный по скоростям w, определяемым в узлах, заменяется расходом, рассчитанным по скоростям w, определяемых на границах контрольного объёма w = we х, + veY, + wez, где ex,,eY,,ez- единичные орты системы координат \X\Y yZ\ связанной с главными осями проницаемости на элементе. Компоненты й и v поля скорости w определяются на элементе / аналогично двумерному случаю, разобранному в главе 2. дР где и - постоянные значения градиентов давления на 1-ом конечном Рассмотрим построение дискретного аналога уравнения неразрывности. Перепишем уравнение неразрывности (3.17) в виде Применим локальный закон сохранения массы к каждому контрольному объему. Тогда для расчетной области, состоящей из М контрольных объемов, получается следующая система уравнений вкладов от конечных элементов, находящихся в плоскости сечения и содержащих данный узел где Vn - вклад от / - го элемента в контрольный объём / - го узла, L - число конечных элементов, содержащих узел /. Применив теорему
Остроградского - Гаусса, перейдём от объёмного интеграла к поверхностному. При вычислении интеграла по поверхности Su (Рис. 4.4), разобьем его на три интеграла соответственно по боковым Su, лицевой »%/ и тыльной S3i поверхностям вклада от конечного элемента в контрольный объем /-го узла J(F,n)dS = J(F,n)dS + J(,n)dS + J(F,n)dS Slt Su S2/ s}/ При вычислении интегралов по лицевой и тыльной поверхностям дР применялась конечно - разностная аппроксимация производной [18], в результате чего 8Z 7—7 Р PD 1= /(F,n)dS = в./„; zz 1М dS. V/ dS, ffruyB = -uf.St ZPU ZP /=1 s}, [P&lu где при вычислении интегралов / pd0\- - dS, \ р" Su \ eff)d S3l zz UVJ dS плотности Pg,p" определяются в точках и и d с помощью линейной аппроксимации, а подвижности V їм vzz Jh )u при помощи гармонической аппроксимации, как показано в [44]. Тогда коэффициенты дискретного аналога определяются следующим образом lPD элементов, находящихся в плоскости сечения и содержащих данный узел где Vn - вклад от / - го элемента в контрольный объём / - го узла, L - число конечных элементов, содержащих узел /. Применив теорему Остроградского - Гаусса, перейдём от объёмного интеграла к поверхностному. При вычислении интеграла по поверхности Su (Рис. 4.4), разобьем его на три интеграла соответственно по боковым Su, лицевой »%/ и тыльной S3i поверхностям вклада от конечного элемента в контрольный объем /-го узла J(F,n)dS = J(F,n)dS + J(,n)dS + J(F,n)dS Slt Su S2/ s}/ При вычислении интегралов по лицевой и тыльной поверхностям дР применялась конечно - разностная аппроксимация производной [18], в результате чего 8Z 7—7 Р PD 1= /(F,n)dS = в./„; zz 1М dS. V/ dS, ffruyB = -uf.St ZPU ZP /=1 s}, [P&lu где при вычислении интегралов / pd0\- - dS, \ р" Su \ eff)d S3l zz UVJ dS плотности Pg,p" определяются в точках и и d с помощью линейной аппроксимации, а подвижности V їм vzz Jh )u при помощи гармонической аппроксимации, как показано в [44].
Тогда коэффициенты дискретного аналога определяются следующим образом lPD Для конвективно — диффузионного уравнения + div(pyt&) — af/v(rgrad ) Дискретный аналог для контрольного объёма вокруг узла Р получается аналогично уравнению неразрывности и будет иметь следующий вид 113 Az(KP - F) + АРФР - АриФри - АроФРО = О Рассмотрим Для конвективно — диффузионного уравнения + div(pyt&) — af/v(rgrad ) Дискретный аналог для контрольного объёма вокруг узла Р получается аналогично уравнению неразрывности и будет иметь следующий вид 113 Az(KP - F) + АРФР - АриФри - АроФРО = О Рассмотрим особенности дискретизации конвективно-диффузионного уравнения на примере уравнения энергии (3.19) при Ф=Г. По оси Z конвективный член уравнения энергии определяется следующим образом J(FT,n2)dS TiJ(F,n2)dS, J(FT,n,)dS TuJ(F,n3)dS. Для определения значения температур в точках и и d используется противопоточная аппроксимация rp,(F,n2) 0 TPU,(F,n2) 0 Td=\ rp,(F,n3) 0 r№,(F,n3) 0 Диффузионный член уравнения энергии вычисляется следующим образом J (-Agrad7 2)fi = f(—Agrad T,n2)dS = P PD = S2I r, T-±fM.dS, JPU Z...-Z.T, где {\zz)ddS , \{\zz)udS определяются с помощью гармонической аппроксимации, как показано в [44]. Введем следующий оператор [А,ВJ, аналогичный оператору тах(А,В)в алгоритмическом языке Fortran. Тогда коэффициенты дискретного аналога можно записать следующим образом
