CC BY
7
УДК 532.685
К. Шаббир, О. Я. Извеков, Б. Вам,си
Московский физико-технический институт (национальный исследовательский университет)
Моделирование пропитки пористой среды с помощью двумерной сетевой модели
Реализована двумерная сетевая модель (network model) двухфазных течений в пеод-ноднородной пористой среде, состоящей из двух подсистем: низко проницаемого блока с тонкими капиллярами, окруженного областью высокопроницаемой среды с толстыми капиллярами. Рассматривается задача установления капиллярного равновесия в результате противоточной пропитки блока. Считается, что в начальный момент времени высокопроницаемая и низко проницаемая части пористой среды насыщены смачивающей и несмачиваюгцей несжимаемой жидкостью соответственно. В численных расчетах на основе сетевой модели исследуется зависимость от времени насыщенности подсистем смачивающей жидкостью и зависимость капиллярного давления от текущей насыщенности. Получено качественное соответствие известным экспериментальным и теоретическим результатам, что в дальнейшем позволит использовать модель для верификации осредненных моделей капиллярной неравновесности.
Ключевые слова: пропитка, капиллярное давление, многофазный поток, сетевые модели
K. Shabbir, O. Ya. Izvekov, B. Vamsi Moscow Institute of Physics and Technology
Modeling imbibition in porous media using a two-dimensional network model
A two-dimensional network model of two-phase flow in a heterogeneous porous media consisting of two regions is implemented: a region of low-permeability with thin capillaries, surrounded by a region of high-permeability with thick capillaries. The case of capillary-equilibrium being established as a result of countercurrent imbibition of a region is studied. Initially, the regions of high-permeability and low-permeability are saturated with wetting and non-wetting incompressible fluids, respectively. Numerical calculations based on this network model examining the dependence of saturation of wetting fluid with respect to time in the regions and the dependence of capillary pressure on the current saturation were performed. A qualitative agreement with known experimental and theoretical results has been obtained, which will further allow the model to be used to verify homogenized models of capillary nonequilibrium.
Key words: imbibition, capillary pressure, multiphase flow, network model
1. Введение
Моделирование двухфазного течения в пористых средах важно во множестве приложений в нефтедобыче, гидрологии, производстве электроэнергии и т.д. [1]. Пористая среда состоит из каркаса (обычно твердого) и пустот (также называемых порами). Пустоты соединены друг с другом капиллярами. В пустотах могут содержаться жидкости, например,
© Шаббир К., Извеков О. Я., Вамси В., 2024
(с) Федеральное государственное автономное образовательное учреждение высшего образования
«Московский физико-технический институт (национальный исследовательский университет)», 2024
вода, нефть или газ [2]. Насыщенность жидкостью Sk определяется как отношение объема ^занимаемого жидкостью F&, к общему объему пустотного пространства VVOid-
Sk = ^ - (1)
'void
Далее будем считать, что пустотное пространство заполнено только двумя жидкостями. Насыщенность более смачивающей жидкости (например, воды в гидрофильной породе) обозначим Sw, а менее смачивающей (например, нефть в гидрофильной породе) Snw. Пусть эти жидкости полностью заполняют пространство пор и капилляров, тогда Sw + Snw = 1.
S
нейной фильтрации флюидов в пористой среде является закон Дарси:
К
Q =--VP, (2)
№
здесь Q — скорость потока; К — проницаемость; № — коэффициент вязкости; a VP — градиент давления.
Проницаемость в случае многофазной фильтрации традиционно считается функцией насыщенности К = К (S) [3], [4].
Однако такое предположение справедливо в том случае, когда характерное время протекания процессов фильтрации намного превышает характерное время релаксации, то есть время перераспределения жидкости в пространстве пор и капилляров за счет сил поверхностного натяжения, и можно считать, что в каждый момент времени конфигурация флюидов соответствует минимуму поверхностной энергии. Таким образом, традиционные модели фильтрации флюидов в пористой среде можно назвать равновесными. Капиллярное равновесие может нарушаться, когда насыщенность меняется относительно быстро или время релаксации к равновесному состоянию достаточно велико, например, при континуальном описании фильтрации в трещиновато-пористой среде [5], [6]. В этих случаях предположение о том, что проницаемость является функцией насыщенности, является недостаточным. Существуют продвинутые модели [7], [8], в которых эффекты капиллярной неравновесно-
К
насыщенности:
К = К (S,^). (3)
Тем не менее процесс перераспределения флюидов в поровом пространстве может про-
S =
К = К(S, О- (4)
к равновесному значению при постоянной насыщенности:
f = ОД О- (5)
Для лучшего уяснения физического смысла неравновесных моделей и уточнения входящих в них параметров полезно рассмотреть движение флюидов в масштабе пор и капилляров. Распространенными методами моделирования движения жидкостей в масштабе пор являются решетчатый метод Больцмана [11], прямое решение уравнения Навье - Стокса, сетевые модели [12]. Прямое моделирование дает достаточно точные результаты по распределению скорости, давления и положениям интерфейсов [13], однако оно требует значительных вычислительных мощностей и времени. Сетевые модели с меньшими трудозатратами позволяют качественно воспроизводить наблюдаемые эффекты [14]. В настоящей работе представлена двумерная сетевая модель, позволяющая воспроизводить неравновесные явления при двухфазном течении в пористой среде [15].
2. Численная модель двухфазного течения несжимаемых жидкостей
Рис. 1. Сотовая модель: (а) расположения узлов и трубок, (б) фрагмент сети, состоящей из трубок различного радиуса. Более темным цветом показана несмачивающая жидкость
Рассматривается модель пористой среды, представляющая собой множество капиллярных трубок, соединенных между собой в узлах рехулярной сетки. Каждый узел соединен с четырьмя соседними узлами как показано на рис. 1а. Узел может быть соединен менее чем с четырьмя узлами, если они расположены в границах области (например, как на рис. 16 и 3). В дальнейшем двойной индекс у обозначает величины, относящиеся к капилляру т^, соединяющему узлы г и ] множества уз лов {Vк}, к = Капиллярные трубки считаются круговыми цилиндрами и в общем случае имеют различный радиус К^. Среда насыщена двумя несжимаемыми жидкостями с различной вязкостью и смачивающими свойствами. Считается, что вся жидкость сосредоточена в капиллярах. Распределение жидкостей в среде определяется положением менисков. В рамках используемой модели количество менисков в одной трубке не превышает двух. Поток жидкости (объемный расход) для трубки Ту есть
Qгj = Ац арЦ + Вц, (6)
где разность давлений на концах трубки АР^ = Рг — Р^,
А. . = г-> лг] =
8Мц I
в = ъК^ 2 в^а
г В^
(8)
здесь I - длина трубки; а - коэффициент поверхностного натяжения; М^ = (р\1\ + ^2h)гj/ - средняя вязкость жидкости в капилляре ггде ^ - коэффициент динамической вязкости жидкости Рк, (1к^ - длина капилляра т^, занимаемая жидкостью Рк, (к = 1, 2). Множитель учитывает количество и ориентацию менисков в капилляре г] при вычислении суммарного скачка капиллярного давления в трубке: = 0, если концы трубки заполнены одинаковой жидкостью, 8^ = 1, если конец трубки в узле г заполнен смачивающей жидкостью, а конец в узле ] - несмачивающей, 8 = —1 в противном случае.
Поскольку Му = а 8^ = — 8^ имеем
Л-- = А ■ ■ •
^г] — •
(9)
В гз — ■
(10)
Скорость потока в капилляре т^ есть
Яц (и)
4 ~ к К2-'
Сохранение полного потока в узлах для закрытой системы (при отсутствии источников/стоков в узлах сетки) приводит к системе N линейных алгебраических уравнений для определения давлений в узлах сетки:
^Яц =0; г = (1, ), (12)
3
здесь ^ — индекс узлов, подключенных к узлу г. Например, на рис. 1, г = 3 и j = (1, 2, 4, 5). С учетом (6) получим
^ ^(Рг - Р3) — - £ Вгз; г — (1, 2,..., N). (13)
з з
Определитель матрицы системы линейных алгебраических уравнений равен нулю (сумми-
0=0
выполняется добавлением постоянной к элементам произвольного столбца к матрицы системы, что соответствует дополнительному условию Рк = 0.
Координаты менисков при их наличии в трубке т^ удовлетворяют уравнению
ж = "«■ <">
Запишем уравнения в безразмерном виде, введя в рассмотрение линейный масштаб ]т, масштаб
ДЕВЛ6НИЯ Рш И коэффициента динамИЧбСКОЙ ВЯЗКОСТИ [Лт'. X — — — ь — tmt.} — ^т^г^ Р — РшР 1 -Мщ — .
, - ]тРт; =]т
ит — ; ь т — . V1"/
т Vт
Уравнения в безразмерном виде есть
(х~.
13 — Угз, (16)
8Мц1 г]
где N — безразмерный параметр, характеризующий соотношение вязких и капиллярных сил:
N _ ^т^т _ ]тРт (18)
а а
Замечание. Параметр N может быть выражен через капиллярное число течения в (наиболее тонком) капилляре радиуса Кге/ и длины 1ге$ при разности давлений на концах капилляра Рт:
*—А г*Й ■ <19)
Кгр /Рт Nr К2 /
N —-, а — — I . (20)
8' ге/а N ыге$Ьт
В задаче капиллярной пропитки выберем в качестве масштаба давления величину Рт — а/Кге$, а в качестве масштаба длины Ьт — Кге/, тогда N — 1 и уравнения (16) и (17) принимают размерный вид. Распределение смачивающей и несмачивающей жидкостей по исходящим потокам в каждом узле определяется из условия минимума энергии контакта несмачивающей жидкости со стенками капилляров. Краткая версия алгоритма:
1) Определение шага интегрирования по времени At.
2) Расчет давления в узлах сетки при текущем положении менисков из решения системы линейных алгебраических уравнений (13).
3) Интегрирование (11) на шаг по времени с использованием явной схемы.
4) Расчет распределения потоков в узлах.
5) Пересчет положения менисков (слияние капель) в случае, если на этапе 3 в трубке появляется третий мениск.
На этане 3 для каждого узла с учетом траектории менисков выполняется расчет объема смачивающей жидкости, поступающей в узел. Полученный объем распределяется по исходящим из узла потокам, в порядке увеличения диаметра трубок. Оставшийся исходящий ноток соответствует неемачивающей жидкости, что гарантируется выполнением общего баланса объемов (13). Эта процедура может приводить к появлению нового мениска. Если на этапе 4 в трубке появляется третий мениск, это означает присутствие в одном капилляре двух капель (участков капилляра), занимаемых каждой жидкостью. В этом случае выполняется перерасчет положения менисков с уменьшением их количества на единицу (слияние капель). Пусть 0 < х\ < Х2 <хз координаты менисков в порядке возрастания. Координаты менисков после слияния x'i < х'2 определяются из условия сохранения объемов жидкостей х'2 — x'i = Х3 — + Х\ и сохранения положения центров масс жидкостей в капилляре х'2 + х' = (ж| — х\ + ж2)/(жэ — х2 + х\). Такая процедура перерасчета сохраняет насыщенность и ноток каждой жидкости в капилляре.
Далее все величины безразмерны, поэтому знак «~» опускается.
0,5 0,4 0,3
>
0,2 0,1 0
2 3 4 5 R
Рис. 2. Распределение капилляров по радиусам во внутренней зоне. V - объемная доля капилляров
3. Эксперименты
На рисунке 16 показан пример расчета фильтрации в области с открытыми границами, на которых поддерживается постоянный перепад давления (градиент давления направлен снизу вверх) [16]. Более смачивающая жидкость показана более светлым тоном.
С помощью разработанного алгоритма была решена задача о пропитке (imbibition) включения внутри области с изолированными границами. Смачивающая жидкость располагалась во внешней области, а несмачивающая жидкость во внутренней области, как показано на рис. За и 4а. В процессе расчета вычислялась насыщенность смачивающей жидкости S (t) во внутренней области и среднее капиллярное давление Рс во внутренней области для различных предельных (равновесных) значений S.
||
В расчетах варьировалось начальное содержание смачивающей жидкости в полной системе.
Рис. 3. Расчет с низкой насыщенностью смачивающей жидкости в системе: (а) начальная конфигурация жидкостей, (б) равновесная конфигурация жидкостей
Рис. 4. Расчет с высокой насыщенностью смачивающей жидкости в системе: (а) начальная конфигурация жидкостей, (б) равновесная конфигурация жидкостей
Расчетная область состоит из сетки 30 х 30 узлов, как показано на рис. 3 и 4. Внешняя ПО ОТНОШеНИЮ К ВКЛЮЧеНИЮ Область СОДерЖИТ Трубки раДИуСОМ Router = 6. Внутренняя область (включение) состоит из более тонких трубок с радиусами в диапазоне Rinner = (2, 3, 4, 5).
Объем внутренней области приблизительно равен объему внешней области Vinner, Vinner ~ V>«ier, где Vsys = Vinner + Vouter- На рис. 2 показано распределение по радиусам
объемной доли трубок V = Vk/Vinner, где Vk - суммарный объем всех трубок одинакового радиуса Rk, к = (2,3,4, 5).
0,0 100,0 200,0, 300,0 400,0 500,0
(а) Насыщенность (б) Капиллярное давление
Рис. 5. (а) зависимость насыщенности смачивающей жидкости Б во внутренней области от времени Ь, (б) среднее капиллярное давление Рс в зависимости от предельной (равновесной) насыщенности смачивающей жидкости Бед во внутренней области
Смачивающая жидкость в основном находится во внешней области, а несмачивающая жидкость во внутренней. Когда насыщенность смачивающей жидкости во всей системе Бвув < УааЬег/Увув ~ 0, 5 (см. рис. За), внешняя область содержит несмачивающую жидкость. А когда Б8у8 > 0, 5, внешняя область содержит только смачивающую жидкость, а остальная часть смачивающей жидкости находится во внутренней области (см. рис. 4а).
К изначальным радиусам трубок Д0 добавляются случайные малые значения, то есть Ду — Д0+АД^. Таким образом, не допускается ситуация, когда к одному узлу примыкают совершенно одинаковые трубки. Значит, всегда существует уникальный способ распределения жидкостей при прохождении узла. В настоящей работе АД^/Д^ ~ ±е — 3. Значения констант модели, использованные в расчетах: — — 0.05, а — 1 и I — 20.
Было проведено 10 расчетов для разных значений Бзув так, чтобы предельная насыщенность смачивающей жидкости во внутренней области попадала в диапазон 0, 02 ^ Б ^ 0, 75. Начальные конфигурации жидкостей для различных Бзув показаны на рис. За и 4а. Во всех расчетах распределение радиусов трубок и константы модели были неизменны.
4. Результаты
Примеры начальной и предельной конфигурации жидкостей в системе трубок показано на рис. 3, 4. Видно, что смачивающая жидкость преимущественно проникает в трубки меньшего радиуса.
На рисунке 5а показана эволюция насыщенности смачивающей жидкости в системе капилляров с течением времени для случаев предельной насыщенности Бед — 0,17 0, 53, 0, 70
насыщенность стремится к равновесному значению.
Расчетные кривые могут быть достаточно хорошо аппроксимированы с помощью функ-
Б(Ь)— Сг +С2 (1 — е-Ся 1), (21)
где константы Сг,С2,Сз, соответствующие расчетным кривым, представленным на рис. 5а, сведены в табл. 1.
Таблица 1
Значения констант кривых тренда на рисунке 5а
Кривая С1 С2 Сз
А 0,00 0,12 0,70е-2
Б 0,16 0,36 1,23е-2
В 0,27 0,37 1,00е-2
На рисунке 56 видно, что среднее равновесное капиллярное давление во внутренней области увеличивается по мере уменьшения равновесной насыщенности смачивающей жидкости во внутренней области, что качественно согласуется с классическими капиллярными кривыми [17]. Из-за геометрических особенностей разработанной сетевой модели трубки с самыми тонкими радиусами занимают малую долю объема внутренней области. В результате капиллярное давление значительно повышается только тогда, когда мениски попадают в тонкие капилляры, то есть при достаточно низкой равновесной насыщенности смачивающей жидкости.
Замечание о накоплении численной ошибки. После решения системы линейных уравнений мы получаем давления в каждом узле. Эти значения не являются точными, поскольку при решении системы линейных уравнений накапливается численная ошибка. Это приводит к незначительному нарушению закона сохранения объема в узлах. Насыщенность конкретной жидкости во всей системе должна оставаться неизменной, поскольку у нас замкнутые границы Ssys = const. Однако в узлах, в которые поступает более одной жидкости, закон сохранения объема слегка нарушается из-за численной ошибки, что приводит к изменению насыщенности конкретной жидкости в системе dSsys/dt ~ 0. Использование типа данных float в С++ позволяет в 2-3 раза ускорить вычисления, но после 104 шагов мы заметили, что ASsys/Ssys ~ 10-2, что мы сочли значительным. Однако когда все действительные числа обрабатываются как double, ошибка уменьшается до ^^Sgyg /'Sgyg ^^ 10 .
5. Заключение
Разработана сетевая двумерная модель для моделирования двухфазного течения жидкостей в пористой среде с учетом капиллярных эффектов. Модель представляет собой регулярную сетку цилиндрических трубок (капилляров), пересекающихся в узлах. Радиус трубок можно варьировать произвольным образом. Считается, что жидкости полностью содержатся в трубках, объемом узлов пренебрегается. Предложен новый метод распределения различных жидкостей при преодолении узлов, который заключается в следующем. Для каждого узла с учетом движения менисков выполняется расчет объема смачивающей жидкости, поступающей в узел. Полученный объем распределяется по исходящим из узла потокам, в порядке увеличения диаметра трубок. Другими словами, в конце расчетного шага происходит выбор конфигурации жидкостей в окрестности узла, соответствующей минимуму поверхностной энергии. Оставшийся объем заполняется несмачивающей жидкостью в соответствии с общим балансом объемов. С помощью разработанной модели численно решена задача о пропитке включения с иными по сравнению с вмещающей средой характеристиками порового пространства (меньший средний радиус капилляров, иное распределение капилляров по радиусам) в условиях, когда внешняя граница вмещающей среды гидродинамически изолирована. Показано стремление насыщенности смачивающей жидкости во внутренней области (включении) к равновесному значению. Построена равновесная капиллярная кривая для включения (зависимость капиллярного давления от равновесной насыщенности).
Максимальное количество соединений, которое может иметь узел в разработанной двумерной сетевой модели, равно 4. Метод распределения жидкостей в узле можно распространить на произвольное количество соединений, что будет использовано при разработке пространственной сетевой модели.
Проект поддержан Российским научным фондом, грант № 23-21-00175.
Список литературы
1. Labed N., Bennamoun L., Fohr J. Experimental study of two-phase flow in porous media with measurement of relative permeability // Fluid Dvn. Mater. Process. Citeseer. 2012. V. 8, N 4. P. 423-436.
2. Su В., Sanchez C., Yang X. Insights into hierarchically structured porous materials: from nanoscience to catalysis, separation, optics, energy, and life science // Hierarchically Structured Porous Materials. Wiley Online Library. 2012. P. 1-27.
3. Coussy O. Poromechanics. John Wiley k, Sons, 2004. P. 315.
4. Kondaypoe В.И. Механика и термодинамика насыщенной пористой среды. Москва : МФТИ, 2007. С. 309.
5. Barenblatt G., Zheltov I., Kochina I. Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks // Journal of applied mathematics and mechanics. Pergamon. 1960. V. 24, N 5. P. 1286-1303.
6. Barenblatt G., Patzek Т., Silin D. The mathematical model of nonequilibrium effects in water-oil displacement 11 SPE journal. 2003. V. 8, N 4. P. 409-416.
7. Hassanizadeh S. Continuum description of thermodynamic processes in porous media: Fundamentals and applications // Modeling Coupled Phenomena in Saturated Porous Materials. 2004. P. 179-223.
8. Hassanizadeh S., Gray W. High velocity flow in porous media // Transport in porous media. Springer. 1987. V. 2. P. 521-531.
9. Kondaurov V.I. The thermodynamicallv consistent equations of a thermoelastic saturated porous medium // Journal of applied mathematics and mechanics. Elsevier. 2007. V. 71, N 4. P. 562-579.
10. Kondaurov V.I. A non-equilibrium model of a porous medium saturated with immiscible fluids // Journal of Applied Mathematics and Mechanics.Elsevier. 2009. V. 73, N 1. P. 88102.
11. Chen S., Doolen G. Lattice Boltzmann method for fluid flows // Annual review of fluid mechanics. 1998. V. 30, N 1. P. 329-364.
12. Blunt M. \et al.\. Pore-scale imaging and modelling // Advances in Water resources. Elsevier. 2013. V. 51. P. 197-216.
13. Liu J., Lin L., Song R., Zhao J. A pore scale modeling of fluid flow in porous medium based on navier-stokes equation // Disaster Advances. 2013. V. 6. P. 5-11.
14. Meakin P., Tartakovsky A. Modeling and simulation of pore-scale multiphase fluid flow and reactive transport in fractured and porous media // Reviews of Geophysics. WTilev Online Library. 2009. V. 47, N 3.
15. Ramstad Т., Berg C., Thompson K. Pore-scale simulations of single-and two-phase flow in porous media: approaches and applications // Transport in Porous Media. Springer. 2019. V. 130. P. 77-104.
16. Шаббир К. Моделирование двухфазного течения в пористых средах с использованием двумерной сетевой модели // Труды 65-й всероссийской научной конференции МФТИ. Москва : Физматкнига, 2023. Т. 78. С. 205-206.
17. Fatt I. The network model of porous media. 3. Dynamic properties of networks with tube radius distribution // Transactions of the American institute of mining and metallurgical engineers. 1956. V. 207, N 7. P. 164-181.
References
1. Labed N., Bennamoun L., Fohr J. Experimental study of two-phase flow in porous media with measurement of relative permeability. Fluid Dvn. Mater. Process. Citeseer. 2012. V. 8, N 4. P. 423-436.
2. Su В., Sanchez C., Yang X. Insights into hierarchically structured porous materials: from nanoscience to catalysis, separation, optics, energy, and life science. Hierarchically Structured Porous Materials. Wiley Online Library. 2012. P. 1-27.
3. Coussy O. Poromechanics. John Wiley k, Sons, 2004. P. 315.
4. Kondaurov V.I. Mechanics and thermodynamics of saturated porous media. Moscow : MIPT, 2007. P. 309. (in Russian).
5. Barenblatt G., Zheltov I., Kochina I. Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks. Journal of applied mathematics and mechanics. Pergamon. 1960. V. 24, N 5. P. 1286-1303.
6. Barenblatt G., Patzek Т., Silin D. The mathematical model of nonequilibrium effects in water-oil displacement. SPE journal. 2003. V. 8, N 4. P. 409-416.
7. Hassanizadeh S. Continuum description of thermodynamic processes in porous media: Fundamentals and applications. Modeling Coupled Phenomena in Saturated Porous Materials. 2004. P. 179-223.
8. Hassanizadeh S., Gray W. High velocity flow in porous media. Transport in porous media. Springer. 1987. V. 2. P. 521-531.
9. Kondaurov V.I. The thermodynamicallv consistent equations of a thermoelastic saturated porous medium. Journal of applied mathematics and mechanics. Elsevier. 2007. V. 71, N 4. P. 562-579.
10. Kondaurov V.I. A non-equilibrium model of a porous medium saturated with immiscible fluids. Journal of Applied Mathematics and Mechanics.Elsevier. 2009. V. 73, N 1. P. 88-102.
11. Chen S., Doolen G. Lattice Boltzmann method for fluid flows. Annual review of fluid mechanics. 1998. V. 30, N 1. P. 329-364.
12. Blunt M. \et al]. Pore-scale imaging and modelling. Advances in Water resources. Elsevier. 2013. V. 51. P. 197-216.
13. Liu J., Lin L., Song R., Zhao J. A pore scale modeling of fluid flow in porous medium based on navier-stokes equation. Disaster Advances. 2013. V. 6. P. 5-11.
14. Meakin P., Tartakovsky A. Modeling and simulation of pore-scale multiphase fluid flow and reactive transport in fractured and porous media. Reviews of Geophysics. WTilev Online Library. 2009. V. 47, N 3.
15. Ramstad Т., Berg C., Thompson K. Pore-scale simulations of single-and two-phase flow in porous media: approaches and applications. Transport in Porous Media. Springer. 2019. V. 130. P. 77-104.
16. Shabbir K. Simulation of Two-Phase Flow in Porous Media using a Two-Dimensional Network Model. Proceedings of the 65th all Russia Scientific Conference MIPT. Moscow : Fizmatkniga, 2023. V. 78. P. 205-206.
17. Fatt I. The network model of porous media. 3. Dynamic properties of networks with tube radius distribution. Transactions of the American institute of mining and metallurgical engineers. 1956. V. 207, N 7. P. 164-181.
Поступим в редакцию 18.04-2024
