CC BY
12
Оригинальная статья
УДК 502/504:626.82:691.11
DOI: 10.26897/1997-6011-2022-3-71-76
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ТЕПЛОВЫХ ПРОЦЕССОВ В СИСТЕМЕ «БЕТОН-ГРУНТ»
ЖАРНИЦКИЙ ВАЛЕРИЙ ЯКОВЛЕВИЧд-р техн. наук, профессор
zharnitskiy@mail.ru
КОРНИЕНКО ПАВЕЛ АЛЕКСАНДРОВИЧ, заведующий лабораторией
kornienko.p.a@mail.ru
СМИРНОВ АЛЕКСАНДР ПЕТРОВИЧ, канд. техн. наук, доцент
sxodnyasmirnov@mail.ru
Российский государственный аграрный университет — МСХА имени К.А. Тимирязева; 127434, г. Москва, Тимирязевская ул., 49, Россия
На расчет режимов выдерживания бетона и разработку математических моделей большое влияние оказывает температурное поле в бетоне облицовок каналов. Задача сводится к решению нелинейного уравнения теплопроводности, учитывающего экзотермические тепловыделения бетона, фазовое превращение влаги, переменные по времени граничные условия, позволяющие учесть воздействие на бетон внешней среды в процессе укладки и ухода за ним. Переменные теплофизические коэффициенты позволяют учесть неоднородность среды (в случае укладки бетона на грунт) и изменение агрегатного состояния вещества при достижении температуры фазового превращения. Поскольку получить аналитическое решение в общем виде невозможно, используется численный метод решения, основанный на сочетании конечно-разностного решения с методом расчета тепловыделений и прочности бетона по соответствующим полям изотермических кривых, полученных экспериментальным путем. При построении разностной схемы используется интегро-интерполяционный метод (метод баланса), основанный на законе сохранения тепла. Для протяженного тела достаточно больших размеров. Процесс переноса тепла в нем принимается линейным, а система координат с центром - на оси тела. Представленная математическая модель тепловых процессов в системе «Бетон-грунт» позволяет прогнозировать режимы выдержки монолитного бетона для достижения необходимых технологических требований, а также применять наиболее экономичные режимы.
Ключевые слова: бетон, облицовка каналов, температурное поле, математические модели, уравнение теплопроводности, тепловыделения бетона, прочность бетона, «Бетон-грунт», фазовое превращение (переход)
Формат цитирования: Жарницкий В.Я., Корниенко П.А., Смирнов П.А Математическое моделирование тепловых процессов в системе «Бетон-грунт» // Природообустройство. - 2022. -№ 3. - С. 71-76. DOI: 10.26897/1997-6011-2022-3-71-76.
© Жарницкий В.Я., Корниенко П.А., Смирнов П.А., 2022
Original article
MATHEMATICAL MODELING OF THERMAL PROCESSES IN THE «CONCRETE-SOIL» SYSTEM
ZHARNITSKIY VALERIJ YAKOVLEVICHH, doctor of technical sciences, professor
zharnitskiy@mail.ru
KORNIENKO PAVEL ALEKSANDROVICH, head of the laboratory
kornienko.p.a@mail.ru
SMIRNOV ALEKSANDR PETROVIICH, candidate of technical sciences, associate professor
sxodnyasmirnov@mail.ru
Russian state agrarian university — MAA named after C.A. Timiryazev; 127434, Moscow, Timiryazevskaya, 49, Russia
The calculation of concrete curing modes and the development of mathematical models are greatly influenced by the temperature field in the concrete of channel linings. The problem is reduced to solving a nonlinear heat equation that takes into account the exothermic heat release of concrete and the phase
transformation of moisture, and time-varying boundary conditions that allow taking into account the impact on concrete of the external environment during laying and maintenance. Variable thermophysical coefficients make it possible to take into account the inhomogeneity of the medium (in the case of laying concrete on the ground) and the change in the aggregate state of the substance when the phase transformation temperature is reached. Since it is impossible to obtain an analytical solution in a general form, a numerical solution method is used, based on a combination of a finite-difference solution with a method for calculating heat release and concrete strength from the corresponding fields of isothermal curves obtained experimentally. When constructing a difference scheme, an integro-interpolation method (balance method) is used, based on the law of conservation of heat. For an extended body of sufficiently large dimensions, the process of heat transfer in it is assumed to be linear, and the coordinate system with the center is taken to be on the axis of the body. The presented mathematical model of thermal processes in the "concrete-soil" system makes it possible to predict the modes of holding monolithic concrete to achieve the necessary technological requirements, as well as to apply the most economical modes.
Keywords: concrete, channel lining, temperature field, mathematical models, heat conduction
equation, concrete heat release, concrete strength, «concrete-soil», phase transformation (transition)
Format of citation: Zharnitskiy V.Ya., Kornienko P.A., Smirnov A.P. Mathematical modeling of thermal processes in the «concrete-soil» system // Prirodoobustrojstvo. - 2022. - № 3. - S. 71-76. DOI: 10.26897/1997-6011-2022-3-71-76.
Введение. Температурное поле в бетоне облицовок каналов объясняется как теплообменом бетона с окружающей средой, так и тепловым взаимодействием бетона с грунтом основания, который может быть промороженным или не промороженным [1-4].
Данное обстоятельство в значительной степени усложняет задачу по расчету режимов выдерживания бетона, и разработка математических моделей этих процессов является важнейшим этапом в оценке температурного поля в бетоне облицовок каналов.
Материалы и методы исследований. Задача определения температурного поля, твердеющего бетона (или двухкомпонентной среды «Бетон-грунт») протяженных конструкций сводится к решению нелинейного уравнения теплопроводности с источниками [5-10], обусловленными экзотермическими тепловыделениями бетона и фазовым превращением влаги, и переменными по времени граничными условиями, позволяющими учесть воздействие на бетон внешней среды в процессе укладки и ухода за ним, то есть
дт
А
дх
A(x,t,T
dt (х,т)
дх
+ ц Q (x,t,T) | д9фП
дт
дт
= (Т
a(Ln ,t,T
дх
Qt (ln ,t) q (т.
-дХ-= QrpN (т;
t (х, O) = tHa4 (х), Q (х, t, O) = Q„H (х, t).
(1)
(2)
(3)
(4)
Переменные теплофизические коэффициенты в уравнениях (1-4) позволяют учесть неоднородность среды (в случае укладки бетона на грунт) и изменение агрегатного состояния вещества при достижении температуры фазового превращения.
Поскольку получить аналитическое решение уравнений (1-4) в общем виде невозможно, используется численный метод решения, основанный на сочетании конечно-разностного решения уравнений (1-4) с методом расчета тепловыделений и прочности бетона по соответствующим полям изотермических кривых, полученных экспериментальным путем. При построении разностной схемы решения уравнений (1-4) используется интегро-интерпо-ляционный метод (метод баланса), основанный на законе сохранения тепла.
Для протяженного тела достаточно больших размеров процесс переноса тепла в нем принимается линейным, а система координат с центром - на оси тела. Для реализации задачи введены следующие обозначения: т - время; Ьо£х£Ьк -линейные размеры тела; g(x, I, т) - удельный вес материала тела; Ц - расход цемента на единицу объема; Ц^(х, t, т) - влажность; Яудфг1 - удельная теплота фазового перехода; ^х, т) - температура; С (х, t, т) - теплоемкость; 2.(х, t, т) - теплопроводность; р(г, t, т) - относительная прочность бетона (% от й28, где й28 - прочность в 28-дневном возрасте); Яэ(у, t, т) - количество тепла, выделенного за счет экзотермии; {Оэ(^, тк)}, е = 0,1,..., пЯ1, к = 0,1,., - совокупность экспериментальных зависимостей тепловыделений от времени при постоянных температурах; {Р(^ тк)}, е = 0,1,., пЯ1, к = 0,1,., тЯ1 - то же для прочности (% от й28);
Ьн(х) - начальное распределение температуры (К); Яэн(х, Ь) - начальное экзотермическое тепловыделение; Рж(х, Ь) - начальное распределение прочности (% от К28).
Расчетное тело разбивается на N элементарных объемов {к; х 1 х 1} с линейными размерами по х{к;}, I = 0,1,..., N-1 и центральными расчетными точками {х}, ; = 0,1,., N-1. Граничные условия на поверхности Ьк(К = 0, N представляются в общем виде как поток тепла через Ьк в единицу времени:
Огк (г) = ^пов.к + ЬГК (Г)" Ь (Хк >Г)
+ Г (0"
Гк (т)-
1 ,^изк W
К
(5)
«к (т) 4зк (т) 2Л(ОкЛ ,т)
где 5изк(т) - толщина изоляционного слоя на поверхности ЬК; Лшк(т) - теплопроводность изоляции; дпое к (т) - поток тепла извне через ЬК; ЬГК (г) внешняя температура возле поверхности ЬК; ак(т) - коэффициент конвективного теплообмена на поверхности ЬК; ГК1 (г), ГК2(г) вспомогательные коэффициенты граничных условий, позволяющие одной формулой представлять граничные условия разного рода.
Пусть Ат - расчетный шаг по т. На сетке
х0 = к + к0/2; х = хг _ 1+(Яг _ 1 + к)/2; I = 1,..., N - 1; т0 = 0; Т = т + Ат; } = 1, 2,.
Принимаются сеточные аналоги функций
К, = Ф)> = Оэн (Х Л; ) М = Ь (Х; Г ),
Я; = о (х ),
Л/ =Я(Х; Л! ),
При обозначении через АЯэ и АР изменения за Ат величин Яз и Р представляется возможным определить ДЯ^.спользуя принцип равных тепловыделений по значению О}-1 на-
Я (1е, г) д^ (¿е+1, г)
ходим —э—---э——-из эксперименталь-
дг дг
ных зависимостей {Яэ(Ье, тк)} таким образом,
чтобы Ье < Ь} < Ье. Поскольку изотермические кривые Яэ(Ье, тк) не пересекаются, предполагаем их пропорциональность значениям Ье. Тогда
AQi -Ат
8Q3 (,т)_
8т
Q (tg+i, т) 8Q3 (te ,т)
8Q (te,т) , 8т
8т
8т
te+1 te
<((' "te)
(6)
•Ат.
Аналогично и для APi
(7)
при достижении температуры фазового перехода 0°С определяется влажностью среды и равно
ДО, .п.; = +^Чд.ф.пЛ (8)
Предполагается, что вся влага в элементарном расчетном объеме, поскольку время мгновенно, превращается в лед (или наоборот) по достижении температуры фазового превращения. Это допущение возможно, поскольку время промерзания (оттаивания) незначительно в сравнении с общей длительностью процесса, а размеры элементарного объема невелики.
Записав начальные и граничные условия, а также уравнения баланса для каждого элементарного объема, получим следующую разностную схему для решения задачи определения температурного поля и полей экзотермических тепловыделений и прочности бетона протяженных конструкций:
^0 = Ьн;, ОО = Я:н;, РО = Р;, I = 0,1.....N-1; (9)
е^оэ+доэ, р;+=Р;+ДР;, ;=о,1.....N-1; (10)
t-j+i" t;+i n L0
h0 + 2Л0 2Л(
-Qi
^а
Aт+AQjЦ • h +AQ; ,ni =
(11)
= Cjrjhj (tj+1" tj ),i - 0;
о' о о V о о /' '
j1 - j t J - j
h+1 , К hi-1 , к
2Л+1 2Л- 2^1 2Л-
Ат + AQj • Ц • V
(12)
+AQ^ ,n.i = Cj (tj1 "tj ),
tj+1 "t
i-2
h
'N-2 _
¿An-2
h
-+Q/
rN-1
W-1
2*N-1
Aт+AQ!^N-1 -Ц-hw-1-
(13)
+AQ(J).n.N-1 - Civ-1^iN-1hN-1 (tij+-11 tNN-1)' i - N -1.
Количество тепла, необходимое для превращения воды в лёд (или наоборот)
Уравнения системы (12-13) определяют неявную двухслойную разностную схему задачи теплопроводности. Она абсолютно устойчива и равномерно сходится со скоростью (Ат + к2тсх), где ктах = max{кj}0£i£N-1. Система (9-13) линейна по индексу и решается методом прогонки. В процессе расчета определяются температурное поле, поле тепловыделений и поле твердеющего бетона, а также влажностное состояние бетона (и грунта в случае укладки бетона на грунт). Это позволяет прогнозировать режимы выдержки монолитного бетона для достижения в заданный срок необходимых
технологических требований, а также выбирать наиболее экономичные режимы. Математическое решение задачи тепловых процессов в системе «Бетон-грунт» положено в основу специализированного программного комплекса.
Результаты и их обсуждение. Расчеты режимов выдерживания бетона облицовок каналов выполнены при начальной температуре бетонной смеси на момент укладки 343, 323 и 298 К и температуре наружного воздуха 253, 263, 273, 293 и 313 К. В целях унификации методов зимнего и летнего бетонирования расчеты для летних условий, как и для зимних, проведены для схемы защиты свежеуложенно-го бетона тепловой изоляцией.
Расчеты выполнялись применительно к механизированным способам укладки. Анализ существующих методов ускоренного твердения бетона показывает, что наиболее предпочтителен форсированный разогрев свежеуложенного бетона и его последующего выдерживания под пассивным или активным термопокрытием, которое перемещается вслед за бетоноукладочным комплексом.
В расчетах принималась длина покрытия 100 м, что обеспечивает в зимних условиях работу комплекса со скоростью передвижения в среднем 10 м / ч. Результаты расчетов для облицовок толщиной 0,15 м представлены на рисунках 1-5. Сравнение данных расчета и результатов физического моделирования показывает, что расхождение по температуре составляет 2...2,5°С, а расхождение во времени достижения нулевой температуры не превышает 7 ч.
На рисунке 1 представлено изменение температуры бетона и грунта при различной глубине промерзания. Анализ полученных результатов показывает, что наиболее благоприятными являются условия укладки бетона при глубине промерзания 0,05 м. При глубине промерзания 0,2 и 0,5 м температура бетона к моменту снятия защитного покрытия на 10.12 К ниже, чем при глубине промерзания 0,05 м. Прочность бетона при промерзании грунта на 0,05, 0,2 и 0,5 м равна соответственно 53, 46 и 44% от й28. Кроме того, прослеживается, что при глубине промерзания грунта 0,2 и 0,5 м температура бетона отличается на 1.2 К. Такое различие в температуре связано с тем, что при глубине промерзания 0,2 м теплосодержание бетона достаточно для осуществления фазового перехода в промерзшем грунте, а при глубине промерзания 0,5 м тепла для фазового перехода по всей толщине промерзшего грунта не хватает, и не оттаявший грунт выполняет
роль термосопротивления. Это подтверждается данными об оттаивании грунта, представленными на рисунке 2, из которого следует, что при промерзании на 0,5 м грунт оттаивает только на глубину 0,3 м. Полученный результат свидетельствует о том, что в случае укладки бетона в зимнее время на основание, которое подготовлено в летних условиях, необходимо знать возможную глубину промерзания грунта для конкретных условий строительства.
Рис. 1. Влияние глубины промерзания грунта (□ - 0,5 м; •о - 0,2 м; ▲ - 0,05 м) на температуру бетона:
1, 4 — температура на поверхности бетона;
2, 5 — температура на контакте с грунтом;
3, 6 - температура грунта Fig. 1. Influence of soil freezing depths (□ - 0.5 m; «о - 0.2 m; ▲ - 0.05 m) on concrete temperature: 1.4 - temperature on the surface of concrete; 2, 5 - temperature in contact with the soil; 3, 6 - soil temperature
Глубина оттаивания, м 0,4]-
0,3
0,2
0,1
у' .--- *
i/a // // It
р н fn
10
20
30B 40
Время, час.
Рис. 2. Характер оттаивания грунта при различной глубине промерзания:
1 - глубина промерзания 0,2 м;
2 - глубина промерзания 0,5 м Fig. 2. The nature of soil thawing
at different freezing depths:
1 - freezing depth of 0.2 m; 2 - freezing depth 0.5 m
Кроме того, расчеты были проведены для условий выдерживания бетона под активным покрытием в предположении того, что после
окончания бетонирования бетон выдерживается в течение суток.
Результаты расчета (рис. 3) показывают, что при небольшой мощности 50.100 ккал / м2 ■ ч скорость остывания бетона в два раза меньше, чем при укрытии бетона обычной изоляцией. К моменту окончания выдерживания температура бетона составляла 309, 300 и 288 К соответственно при мощности 100,5 ккал / м2 ■ ч и защите изоляцией толщцной 0,05 м. В аналогичных условиях выдерживания прочность бетона к моменту его замерзания составляла 61, 56 и 45% от И,
'28'
Рис. 3. Влияние способов ухода за бетоном, нормирование температурного поля бетона:
2, 2* - тепловой поток 70 ккал/м2 ■
ч;
3, 3* — тепловой поток 50 ккал/м ■ ч Fig. 3. Influence of methods of concrete care, normalization of the temperature field of concrete:
1, 1 * — heat flow of 100 kcal / m2 ■ h;
2, 2* — heat flow of 70 kcal/m2 ■ h;
3, 3 * — heat flow 50 kcal/m2 ■ h
После отключения энергии от активного покрытия остывание бетона происходит аналогично остыванию бетона, выдерживаемого под пассивной изоляцией.
Расчеты для разных начальных температур разогрева бетона (323 и 343 К) свидетельствуют о том, что к моменту окончания выдерживания прочность бетона отличается на 16% от R28 (рис. 4). Поэтому при низких температурах окружающей среды бетон необходимо разогревать до максимально возможной температуры (рис. 5).
Результаты, представленные на рисунках 1-5, свидетельствуют о том, что при отрицательных температурах в случае использования составов бетона, применяемых в летний период строительства, прочность бетона к моменту его замерзания практически не превышает 50% от R28.
Поскольку расчеты выполнялись при максимально возможных параметрах
выдерживания бетона, выявляется необходимость при более высоких требованиях к бетону по прочностным показателям проектировать такой состав, при котором бетон в обеспечиваемом температурном режиме наберет требуемую к моменту замерзания прочность.
Температура, К
333
313
XL - - - «1»
Vv V * ч. / / J -< J- - - - _ ^i — H
г / * |"ЧГ 11 1 . . x 1 г
9 У,л г
Время, час.
Рис. 4. Влияние температуры предварительного разогрева на прочность и температуру бетона:
1, 2, 3 — температура бетона; 1*, 2*, 3* — прочность бетона, % от R28 Fig. 4. Influence of preheating temperature on strength and temperature of concrete:
1, 2, 3 — temperature of concrete; 1*, 2*, 3* — concrete strength in % of R28
Рис. 5. Влияние температуры наружного воздуха на формирование температурного поля бетона:
1, 2, 3 — температура бетона;
1*, 2*, 3* — прочность бетона, % от R28 Fig. 5. The influence of outdoor temperature on the formation of the temperature field of concrete:
1, 2, 3 — the temperature of concrete;
1*, 2*, 3* — concrete strength in % of R28
Выводы
Представленная математическая модель тепловых процессов в системе «Бетон-грунт» позволяет прогнозировать режимы выдержки монолитного бетона для достижения необходимых технологических требований, а также применять наиболее экономичные режимы.
|7б|
1, 1* — тепловой поток 100 ккал/м ■ ч
Библиографический список
1. Грозав В.И., Кулиев К.А. Интенсификация бетонирования монолитных облицовок каналов // Материалы научно-технической конференции «Природо-обустройство сельскохозяйственных территорий». - М.: МГУП, 2001. - 85 с.
2. Жарницкий В.Я., Корниенко П.А. Обоснование линейной математической модели сдвигового течения бетонной смеси под действием силы тяжести по наклонной поверхности откоса // Природообустрой-ство. - 2020. - № 1. - С. 88-93.
3. Жуков В.Д., Гаврилина НА., Никитин И.А. Разработка рецептур и режимов тепловой обработки теплоизоляционных плитных изделий из перлита // Сб .трудов № 25. Совершенствование режимов тепловой обработки строительных материалов. - М.: РосНИИНСа, 1962.
4. Интенсификация структурообразования бетонов автоклавного твердения / И.В. Удачкин, Ю.Д. Ка-циевский, П.И. Захарченко и др. // Строительные материалы. - 1984. - № 1.
5. Борделяну Г.В. Экспериментально-статистические исследования деформаций ползучести заводского бетона с построением математических моделей второго порядка для их вычисления и прогнозирования: Авто-реф. дис. ... канд. техн. наук. - Кишинев, 1974. - 23 с.
6. Вайнберг А.И. Решение динамических задач оползневого сдвига // Гидротехническое строительство. - 2002. - Вып. 5. - С. 24-28.
7. Калашников В.И. Как превратить бетоны старого поколения в высокоэффективные бетоны нового поколения // Бетон и железобетон. - 2012. - № 1. - С. 82.
8. Краснов Л.М., Федосов С.В., Акулова М.В. Влияние высокого наполнения мелкозернистого бетона на структурную прочность // Строительные материалы. - 2009. - № 1. - С. 48-50.
9. Жарницкий В.Я., Корниенко П.А. Экспериментальные исследования безвибрационного метода укладки бетонной смеси в облицовках каналов // При-родообустройство. - 2021. - № 3. - С. 88-94.
10. Жарницкий В.Я., Корниенко П.А. Кинетика гидратации цемента, пластическая прочность бетона облицовки канала и ее термонапряженное состояние // Природообустройство. - 2021. - № 5. - С. 85-90.
Критерии авторства
Жарницкий В.Я., Корниенко П.А., Смирнов А.П. выполнили теоретические исследования, на основании которых провели обобщение и написали рукопись, имеют на статью авторское право и несут ответственность за плагиат. Конфликт интересов
Авторы заявляют об отсутствии конфликтов интересов Статья поступила в редакцию 14.04.2022 Одобрена после рецензирования 11.05.2022 Принята к публикации 24.05.2022
References
1. Grozav V.I., Kuliev К.А. Intensifikatsiya betoni-rovaniya monolitnyh oblitsovok kanalov: Materialy nauch-no-tehnicheskoj konferentsii. - M.: MGUP, 2001. - 85 s.
2. Zharnitsky V.Ya., Коrnienko P.Ä. Obosnovanie linejnoj matematicheskoj modeli sdvigovogo techeniya betonnoj smesi pod dejstviem sily tyazhesti po naklon-noj poverhnosti otkosa // Prirodoobustrojstvo. - 2020. -№ 1. - S. 88-93.
3. Zhukov V.D., Gavrilina NA., Nikitin IA Rarabotka retseptur i rezhimov teplovoj obrabotki tep-loizolyatsionnyh plitnyh izdelij iz perlita. Sb. trudov № 25. - M.: RosNIINSa, 1962.
4. Intensifikatsiya strukturoobrazovaniya beto-nov avtoklavnogo tverdeniya / Udachkin I.V., Katsiev-skyj Yu.D., Zakharchenko P.I. i dr. // Stroitelnye materialy. - 1984. - № 1.
5. Bordelyanu G.V. Experimentalno-statisticheskie issledovaniya deformatsij polzuchesti zavodskogo betona s postroeniem matematicheskih modelej vtorogo poryad-ka dlya ih vychisleniya i prognozirovaniya: avtoreferat dis. kand. tehn. nauk. - Kishinev: 1974. - 23 s.
6. Vainberg А.I. Reshenie dinamicheskih za-dach opolznevogo sdviga // Gidrotehniheskoe stroi-telstvo. - 2002. - Vyp. 5. - S. 24-28.
7. Каlashnikov V.I. Kak prevratit betony starogo pokoleniya v vysokoeffektivnye betony novogo pokole-niya // Beton i zhelezobeton. - 2012. - № 1.
8. ^asnov L^., Fedosov S.V., Äkulova М.V. Vliya-nie vysokogo napolneniya melkozernistogo betona na struk-turnuyu prochnost // Stroitelnye materialy. - 2009. -№ 1. - S. 48-50.
9. Zharnitsky V.Ya., ^mien^ PA. Experi-mentalnye issledovaniya bezvibratsionnogo metoda ukladki betonnoj smesi v oblitsovkah knalov // Prirodoobustrojstvo. - 2021. - № 3. - S. 88-94.
10. Zharnitsky V.Ya., ^mien^ PA. Kinetika gidratatsii tsementa, plasticheskaya prochnost betona oblitsovki kanala i ee termonapryazhennoe sostoyanie // Prirodoobustrojstvo. - 2021. - № 5. - S. 85-90.
Criteria of Authorship
Zharnitskiy V.Ya., Kornienko P.A., Smirnov A.P. carried out theoretical studies, on the basis of which they generalized and wrote the manuscript. Zharnitskiy V.Ya., Kornienko P.A., Smirnov A.P. have a copyright on the article and are responsible for plagiarism. Conflict of interests
The authors declare that there are no conflicts of interests The article was submitted to the editorial office 14.04.2022 Approved after reviewing 11.05.2022 Accepted for publication 24.05.2022
