На правах рекламы:
ISSN 0236-235X (P)
ISSN 2311-2735 (E)

Авторитетность издания

ВАК - К1
RSCI, ядро РИНЦ

Добавить в закладки

Следующий номер на сайте

4
Ожидается:
09 Декабря 2024

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

Optimization problems solution based on superelement modeling of oil-field development
Дата подачи статьи: 30.03.2017
УДК: 622.276.1/4
Статья опубликована в выпуске журнала № 3 за 2017 год. [ на стр. 384-391 ]
Аннотация:Нефтяные месторождения России разрабатываются преимущественно с помощь заводнения. Большинство из них находятся на третьей или четвертой стадии разработки. Следовательно, обводненность продукции скважин составляет 80–90 % и более. В этих условиях с целью оптимизации разработки месторождений инженеры-нефтяники стараются уменьшить добычу и закачку воды при сохранении или увеличении добычи нефти. Для этого решаются задачи контроля и регулирования разработки месторождения с использованием различных математических моделей. В данной работе рассматривается суперэлементная математическая модель заводнения, основанная на модели двухфазной фильтрации слабосжимаемых несмешивающихся жидкостей (нефти и воды) в упругом пласте по закону Дарси. Система дифференциальных уравнений для давления и насыщенности аппроксимирована на сетке Вороного полностью явным образом. Размеры суперэлементов сопоставимы с расстоянием между скважинами. Это позволяет проводить вычисления без использования специального ПО. Для решения обратных задач (определения коэффициентов модели или оптимизации параметров разработки) в работе используются метод Ньютона и метод сопряженных градиентов. В классической постановке обратной задачи методы теории оптимизации должны применяться непосредственно к математической модели исследуемого процесса. Однако в случае решения задач разработки нефтяных месторождений количество параметров для оптимизации может быть очень велико, а сложность математической модели довольно высока, поэтому применение методов теории оптимизации непосредственно к математической модели может быть очень трудоемким. Для преодоления этого противоречия предлагается с помощью математической модели заводнения строить статистические зависимости показателей разработки от искомых параметров, а затем применять методы теории оптимизации уже не к математической модели, а к полученным статистическим зависимостям. Для иллюстрации такого подхода рассматривается решение задачи адаптации модели по абсолютной проницаемости. Установлено, что применение метода сопряженных градиентов непосредственно к модели заводнения дает ошибку в определении проницаемости 11,8 %. Применение того же метода к статистической зависимости ошибки адаптации модели (по накопленной добыче и закачке нефти и воды) от логарифма проницаемости дает ошибку в определении проницаемости лишь немногим больше – 15 %.
Abstract:Oil fields of Russia are mainly developed by waterflooding. Most of them are in the 3rd or 4th stage of development. Consequently, the water cutting of well production is 80–90% or more. In these conditions, in order to optimize the development of deposits, oil engineers try to reduce water production and injection while maintaining or increasing oil production. For this purpose, there are the tasks of field development control and regulation. These problems are solved using various mathematical models. This paper considers a superelement mathematical waterflooding model based on a two-phase filtration model for weakly compressible immiscible liquids (oil and water) in an elastic bed under the Darcy law. The system of differential equations for pressure and saturation is approximated on Voronoi diagram in an entirely explicit manner. The size of the superelements is comparable to the distance between the wells. This allows performing calculations without using special software. To solve inverse problems (determine model coefficients or optimize development parameters), the work uses Newton's method and the conjugate gradient method. In the classical setting of an inverse problem, the optimization theory methods should be applied directly to the mathematical model of the process under study. However, when solving oilfield development problems, the number of optimization parameters can be large, and the complexity of the mathematical model is quite high. Therefore, the application of optimization theory methods directly to a mathematical model can be very time-consuming. To overcome this contradiction, it is proposed to build statistical dependencies of the development indices on the required parameters using a mathematical model of waterflooding, and then to apply optimization theory methods no longer to the mathematical model, but to the statistical dependencies obtained. To illustrate this approach, we consider the solution of the problem of model adaptation to absolute permeability. It is established that the application of the conjugate gradient method directly to the waterflooding model gives an error in determining the permeability of 11,8 %. Applying the same method to a statistical dependence of a model adaptation error (on accumulated production and injection of oil and water) on the logarithm of permeability gives an error in determining permeability of only a little more, it is 15 %.
Авторы: Афанаскин И.В. (ivan@afanaskin.ru) - Федеральный научный центр Научно-исследовательский институт системных исследований РАН (ФНЦ НИИСИ РАН) (зав. группой), Москва, Россия, кандидат технических наук, Ялов П.В. (petryalov@gmail.com) - Федеральный научный центр Научно-исследовательский институт системных исследований РАН (ФНЦ НИИСИ РАН) (инженер), Москва, Россия, Гиацинтов А.М. (algts@inbox.ru) - Центр визуализации и спутниковых информационных технологий НИИСИ РАН (cтарший научный сотрудник), Москва, Россия, Родителев А.В. (avrod_94@mail.ru) - Федеральный научный центр Научно-исследовательский институт системных исследований РАН (ФНЦ НИИСИ РАН) (ведущий программист), Москва, Россия
Ключевые слова: обратные задачи, сетка вороного, суррогатная модель, суперэлементная модель, оперативное моделирование, экспресс-моделирование, заводнение
Keywords: inverse problems, voronoy grid, surrogate model, superelement model, prompt modeling, express modeling, waterflooding
Количество просмотров: 8572
Статья в формате PDF
Выпуск в формате PDF (21.91Мб)
Скачать обложку в формате PDF (0.59Мб)

Размер шрифта:       Шрифт:

Основным методом разработки нефтяных месторождений в России является заводнение. Большинство нефтяных месторождений находятся на третьей или четвертой стадии разработки. Следовательно, в продукции добывающих скважин присутствует большое количество воды. В этих условиях инженерам по разработке месторождений необходимо решать много задач по контролю и регулированию разработки. Для быстрого оценочного решения подобных задач может быть использована модель на базе концепции суперэлементов. Из-за большой неоднозначности исходных данных необходимо адаптировать модель к фактически замеренным показателям разработки (дебит нефти, газа и воды; забойное давление; пластовое давление и пр.). Задача адаптации модели сводится к решению задачи оптимизации. Кроме того, с помощью теории оптимизации и математического моделирования решают задачи оптимизации разработки нефтяных месторождений (оптимизации заводнения). Данная работа посвящена подходу, позволя- ющему с помощью простых методов (метод Нью- тона, метод сопряженных градиентов) решать задачи оптимизации при суперэлементном моделировании разработки нефтяных месторождений. Используя только Excel Microsoft, описываемым методом можно решать задачи оптимизации с количеством параметров до 16. Используя специальные программы, количество параметров можно увеличить.

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

Система уравнений, описывающая упругую двухфазную фильтрацию, состоит из двух уравнений сохранения количества основных компонентов (объемов нефти и воды в стандартных условиях) и обобщенного закона Дарси, а капиллярными и гравитационными силами пренебрегается [1–6]:

– уравнения сохранения объемов нефти и воды

,                             (1)

;                    (2)

– обобщенный закон Дарси

,                                        (3)

,                                       (4)

где m –  пористость; So и Sw, Bo и Bw,  и ,  и , kro и krw, mo и mw – насыщенность, объемный коэффициент, вектор скорости фильтрации, плотность источника (стока), моделирующего работу скважины, относительная фазовая проницаемость, вязкость нефти и воды соответственно;  – плотность источника (стока), моделирующего работу законтурной области; k – абсолютная проницаемость; P –  пластовое давление.

Систему уравнений (1)–(4) необходимо дополнить замыкающими соотношениями:

So + Sw =1,                                                           (5)

m = m0[1 + Cr(P – P0)],                                      (6)

Bo = Bo0[1 – Co(P – P0)],                                          (7)

Bw = Bw0[1 – Cw(P – P0)],                                   (8)

где m0 – пористость при начальном пластовом давлении; Bo0 и Bw0 – объемный коэффициент нефти и воды при начальном пластовом давлении; Cr, Co и Cw – сжимаемость пласта, нефти и воды; P0 – начальное пластовое давление.

Систему уравнений (1)–(4) с учетом соотношений (5)–(8) можно преобразовать к следующей системе дифференциальных уравнений для водонасыщенности S º Sw и давления P [7]:

                       (9)

                           (10)

Уравнения (9) и (10) дополняются начальными условиями P = P(x, y, z, t = 0), S = S(x, y, z, t = 0) и граничными условиями непротекания на внешних границах. Для моделирования законтурной водоносной области используются источниковые слагаемые  в уравнениях (9) и (10) [1].

Численная схема

Рассмотрим расчетную сетку, составленную из так называемых суперэлементов [8, 9] – ячеек с размерами, в плане сопоставимыми с расстоянием между скважинами (300–500 м). Тогда количество ячеек в модели будет примерно равно количеству скважин. Скважины используются как центры ячеек. В неразбуренной части объекта могут быть введены фиктивные скважины для построения сетки. Использование такой сетки позволяет сократить вычислительные затраты в тысячи раз [8].

Задача построения суперэлементной сетки в плане (2D-сетки) сводится к задаче построения для каждой скважины зоны дренирования, которая бы­ла введена А.П. Крыловым, а она, в свою очередь, к построению областей Вороного [7] (рис. 1):

-     нанесение на карту скважин и границ;

-     построение выпуклой оболочки по внешним скважинам;

-     построение триангуляции Делоне по скважинам;

-     построение областей Вороного по триангуляции.

Для преобразования суперэлементной сетки в 3D-сетку построенную в плане сетку копируют для разных слоев, изменяя при этом глубину залегания кровли ячеек и их толщину.

Большой размер ячеек позволяет использовать полностью явную схему [10]. Рассмотрим такую схему для уравнений (9) и (10) на сетке суперэлементов в плане (2D-плоская модель):

где n – номер шага по времени; i – номер ячейки; Dtn+1 – шаг по времени; Fi, hi и (m0)i – площадь, толщина и пористость при давлении P0 ячейки i;  – дебит фазы a = o, w скважины в ячейке i на шаге по времени n;  – приток воды из законтурной области в ячейку i;  – переток фазы a = o, w между ячейками i и j на шаге по време- ни n; Ni – количество соседних ячеек для ячейки i.

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

Переток между ячейками i и j (рис. 2) определяется как , a = o, w, где Aij – длина совместного ребра ячеек i и j; Lij – расстояние между центрами ячеек (скважинами) i и j; ;

 

; .

Моделирование скважин

Дебиты скважин по жидкости ql, нефти qo, воде qw, обводненность W и забойное давление Pw опре- деляются из следующих соотношений:

,

,

,

,

где  – коэффициент продуктивности скважины i по жидкости, который определяется как

,

где (rw)i – радиус скважины по долоту;  – скин-фактор скважины; (Rc)i – эквивалентный радиус блока, определяемый как .

Для добывающих скважин необходимо задать дебит жидкости или забойное давление, для нагнетательных – расход воды или забойное давление.

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

где (Sc)i – водонасыщенность на фронте вытес- нения; (Swcr)i – насыщенность связанной водой; (Sowcr)i – насыщенность остаточной нефтью; коэффициенты Ai, Bi, ai, bi определяются по результатам интерпретации исследований керна и могут зависеть от пористости или проницаемости.

Водонасыщенность на фронте вытеснения (Sc)i  определяется из соотношения

,

где f(Si) – функция Бакли–Леверетта:

.

Предложенная 2D-расчетная схема легко обобщается на случай 3D. При этом приток в скважину из каждого слоя определяется согласно его продуктивности при известном дебите скважины по жидкости или забойном давлении.

Решение обратных задач методом Ньютона

Пусть необходимо найти минимум функции многих переменных f(X), где X = (x1, x2, x3, …, xn). Эта задача эквивалентна задаче нахождения значений X, при которых градиент функции f(X) равен нулю:

grad(f(X)) = 0.                                                    (11)

Применим к (11) метод Ньютона:

,        (12)

где j = 1, 2, 3, …, m – номер итерации; H(X) – гессиан функции f(X).

Напомним, что гессиан функции – это симметричная квадратичная форма, описывающая поведение функции во втором порядке:

где , функция f(X) задана на n-мерном пространстве вещественных чисел.

В более удобном для вычислений виде формулу (12) можно представить в виде

.

Метод Ньютона прост в реализации, но нахождение матрицы Гессе сопряжено с большими вычислительными затратами из-за необходимости нахождения большого количества частных производных.

Метод сопряженных градиентов

Метод сопряженных градиентов является методом нахождения локального минимума (максимума) функции с использованием информации о ее значениях и градиенте.

Метод сопряженных градиентов математи- чески существенно сложнее метода Ньютона. Не углубляясь в особенности метода, рассмотрим итерационный алгоритм нахождения экстремума функции: , где k – номер направления поиска; j – номер итерации, условие перпендикулярности направлений поиска k+1 и k: , направление поиска: , где .

Метод сопряженных градиентов математически сложнее метода Ньютона, но требует меньше вычислений.

Метод Ньютона и метод сопряженных элементов реализованы в надстройках Excel Microsoft.

Программная реализация математической модели двухфазной фильтрации на основе концепции суперэлементов

Приведенная выше математическая модель, обобщенная на случай 3D, была описана в среде разработки Microsoft Visual Studio 15 на языке программирования C++. Для представления исходных данных, описывающих геолого-физические свойства пласта, был выбран формат в виде таблицы с использованием ключевых слов для обозначения блоков данных. Выбор представления обусловлен возможностью редактирования данных пользователем без использования специализированных программных средств (например в Microsoft Excel), а также экспорта данных из других симуляторов. При редактировании большого объема исходных данных пользователю необходимо хранить метаданные, относящиеся к ячейкам таблицы. Для реализации предложен алгоритм фильтрации, который удаляет метаданные, маркированные символами «//» из каждого блока при считывании файла. После успешного считывания исходных данных производятся проверки на целостность данных, размеры связанных блоков, размеры связанных интервалов перфораций и слоев, отсутствие пустых ячеек внутри таблицы. При успешном считывании и прохождении всех проверок пользователю выводится список выполненных операций. В случае ошибки пользователь также получает сообщение с указанием позиции в блоке данных и способах устранения ошибки.

На следующем этапе выполняется алгоритм моделирования заводнения нефтяных месторождений с применением многопоточности. Создаются отдельные потоки для расчета статических данных, а затем очередь обработки функций, рассчитывающих динамические данные. На данном этапе производится расчет изменяемого шага по времени; в целях сокращения используемой оперативной памяти хранятся данные только по трем шагам. Далее происходит формирование выходного файла, в котором присутствуют ключевые слова, описание размерности данных, диапазоны адресов ячеек данных для блоков, а также нумерация скважин и слоев. Предлагаемый формат выходных данных позволяет произвести их анализ и обработку (например в Microsoft Excel), в том числе построить диаграммы без дополнительной подготовки.

Созданная программа получила название Oil production calculator, или сокращенно OPC. Она была проверена с помощью сравнительного тестирования с симуляторами Dz10 Каца Р.М. и Волгина Е.Р., Rubis Kappa Engineering, Eclipse Schlumberger. Получено хорошее совмещение результатов расчетов с помощью ОРС и результатов, полученных с помощью вышеописанных программ.

Новый подход к решению задач оптимизации при моделировании разработки нефтяных месторождений

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

Рассмотрим простейшую секторную модель с 4 скважинами – 1 нагнетательной и 3 добывающими. Параметры модели следующие:

-     размеры модели в плане – 1 500´1 500 м;

-     абсолютная проницаемость – переменная (см. табл. 1);

-     пористость – переменная;

-     начальное пластовое давление – 339 атм.;

-     начальная водонасыщенность – 0,2 д.ед.;

-     объемный коэффициент воды – 1,02 м3/м3; вязкость воды – 0,36 сПз; сжимаемость воды – 4,7*10-5 1/атм.;

-     объемный коэффициент нефти – 1,55 м3/м3; вязкость нефти – 0,397 сПз; сжимаемость нефти – 1,01´10-4 1/атм.; растворимость газа в нефти – 203,5 м3/м3; давление насыщения нефти газом – 18,8 атм.;

-     сжимаемость породы – 4,7*10-5 1/атм.;

-     капиллярным давлением пренебрегается;

-     функции относительных фазовых проницаемостей (ОФП) приведены на рисунке 3;

-     дебит жидкости добывающих скважин – данные в таблице 1; минимальное забойное давление в добывающих скважинах – 18,8 атм.;

-     закачка воды в нагнетательную скважину без ограничений; забойное давление в нагнетательной скважине – данные в таблице 1;

-     скин-фактор для всех скважин – 0 ед.;

-     радиус скважин – 0,2 м;

-     срок разработки – 35 лет;

-     притока воды из-за контура нет.

Таблица 1

Параметры модели

Table 1

Model parameters

Параметр

Скважина (суперэлементная ячейка)

1

2

3

4

Площадь, тыс. м2

514

628

603

504

Толщина, м

25

22

19

15

Пористость, д. ед.

0,20

0,18

0,16

0,15

Абсолютная проницаемость, мД

70

53

44

38

Дебит, м3/сут.

-

236

168

114

Забойное давление, атм.

372

18,8

18,8

18,8

Примечание: 1 – нагнетательная скважина, 2, 3, 4 – добывающие.

Добывающие скважины управлялись дебитом, а нагнетательная – забойным давлением.

Расстановка скважин в плане и конфигурация суперэлементов показаны на рисунке 4.

Сведения о длине совместного ребра соседних ячеек для суперэлементной модели и расстоянии между центрами соседних ячеек (скважинами) приведены в таблице 2.

Таблица 2

Геометрические параметры суперэлементной модели

Table 2

Geometric parameters of a superelement model

Параметр

Пары ячеек

1–2

1–3

2–3

2–4

3–4

Длина совместного ребра ячеек, м

650

538

425

650

663

Расстояние между центрами соседних ячеек, м

638

889

925

1013

813

Итак, в данном случае задача оптимизации состоит в том, чтобы уточнить значения проницаемости, используемые в модели, зная историю разработки. Имея оценку проницаемости по данным геофизических и гидродинамических исследований, а также исследований керна, с помощью генератора случайных чисел рассчитаем n вариантов проницаемости по 4 штуки в каждой группе k1, k2, k3, k4, где цифры 1–4 соответствуют номерам ячеек. Примем n = 20. Рассчитаем значения показателей модели для каждой группы значений проницаемости. Запишем целевую функцию следующим образом:

где Qo – накопленная добыча нефти по скважине; Qw – накопленная добыча воды по скважине; Qiw – накопленная закачка воды; индекс c – расчетные показатели; индекс f – фактические показатели; индекс i – номер ячейки. Для изучения эффективности предлагаемого метода в качестве фактических показателей примем показатели, полученные путем моделирования при значениях проницаемостей k1 = 70, k2 = 53, k3 = 44, k4 = 38 мД.

Таким образом, для каждого набора проницаемостей ячеек 1–4 получено значение целевой функции F. Из опыта анализа разработки известно, что накопленная добыча и закачка прямо пропорциональны проницаемости либо ее логарифму. Составим следующие уравнения для целевой функции F:

,                      (13)

               (14)

         (15)

где a1–a9 – постоянные коэффициенты.

С помощью регрессионного анализа найдем коэффициенты соотношений (13)–(15), используя полученную таблицу F от k1, k2, k3, k4. Затем, применив метод сопряженных градиентов к (13)–(15), определим проницаемости k1, k2, k3, k4, приняв начальное приближение проницаемости 50 мД (для такого маленького числа параметров (4 штуки) не важно, какой метод использовать – Ньютона или сопряженных градиентов). Сравним результаты определения проницаемости, полученные с помощью применения методов Ньютона и сопряженных градиентов непосредственно к модели фильтрации, а также полученные путем применения метода сопряженных градиентов к регрессионным моделям (13)–(15) (см. табл. 3). Начальное приближение проницаемости – 50 мД. Для оценки достоверности определения проницаемости введем ошибку по проницаемости следующим образом:

,                                      (16)

где индекс c – расчетная проницаемость; индекс f – фактическая проницаемость; индекс i – номер ячейки. Фактически соотношение (16) представляет среднюю относительную ошибку определения проницаемости.

Из таблицы 3 видно, что наиболее точные значения проницаемости можно получить, применяя метод сопряженных градиентов к модели фильтрации (ошибка – 11,8 %). Однако в случае, когда это невозможно (например, из-за большого числа иско- мых параметров), используя регрессионную модель и метод сопряженных градиентов (либо метод Ньютона), можно получить значения проницаемости с ошибкой не более 15,0 %. При этом важно правильно выбрать вид регрессионной модели.

Заключение

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

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

Работа выполнена при поддержке Программы фундаментальных исследований РАН № I.33П, НИР № 0065-2015-0111.

Литература

1.     Азиз Х., Сеттари Э. Математическое моделирование пластовых систем. М.–Ижевск: Изд-во Ин-та комп. исследований, 2004. 416 с.

2.     Каневская Р.Д. Математическое моделирование гидродинамических процессов разработки месторождений углеводородов. М.–Ижевск: Изд-во Ин-та комп. исследований, 2002. 140 с.

3.     Кац Р.М., Волгин Е.Р., Афанаскин И.В. Численное моделирование двухфазной фильтрации нефти и воды // Тр. НИИСИ РАН, 2014. Т. 4. Вып. 2. С. 141–148.

4.     Coats K.H., Thomas L.K., Pierson R.G. Compositional and black oil reservoir simulation. SPE 29111, 13th SPE Sympos. on Reservoir Simulation, 1995, pp. 149–162.

5.     Crichlow G.B. Modern reservoir engineering – a simulation approach. New Jersey: Prentice Hall. Inc., 1979, p. 354.

6.     Peaceman D.W. Fundamentals of numerical reservoir simulation. Amsterdam: Elsevier Scientific Publishing Co., 1977, p. 191.

7.     Хисамутдинов Н.И., Хасанов М.М., Телин А.Г. [и др.]. Разработка нефтяных месторождений. В 4 т. Т. 1. Разработка нефтяных месторождений на поздней стадии. М.: Изд-во ВНИИОЭНГ, 1994. 240 с.

8.     Мазо А.Б., Булыгин Д.В. Суперэлементы. Новый подход к моделированию разработки нефтяных месторождений // Нефть. Газ. Новации. 2011. Вып. 11. С. 6–8.

9.     Булыгин Д.В., Мазо А.Б., Поташев К.А., Калинин Е.И. Геолого-технические аспекты суперэлементной фильтрационной модели нефтяных месторождений // Георесурсы, 2013. Вып. 3. С. 31–35.

10.   Афанаскин И.В., Егоров А.А., Колеватов А.А. Экспресс-моделирование заводнения нефтяных месторождений с помощью концепции суперэлементов // Вестн. кибернетики. 2016. Вып. 2. С. 153–163.


Постоянный адрес статьи:
http://swsys.ru/index.php?page=article&id=4306
Статья в формате PDF
Выпуск в формате PDF (21.91Мб)
Скачать обложку в формате PDF (0.59Мб)
Статья опубликована в выпуске журнала № 3 за 2017 год. [ на стр. 384-391 ]

Возможно, Вас заинтересуют следующие статьи схожих тематик: