CC BY
2
УДК 658.5
ТЕРМОДИНАМИЧЕСКАЯ МОДЕЛЬ ЛЕДОВОГО ВОРОТНИКА В ГИДРОТЕХНИЧЕСКОМ СТРОИТЕЛЬСТВЕ
© 2023 Д.А. Шарапов, Ю.С. Клочков
Санкт-Петербургский политехнический университет Петра Великого (национальный исследовательский университет), г. Санкт-Петербург, Россия
Статья поступила в редакцию 08.04.2023
Одной из приоритетных задач развития страны является освоение Арктики. Для этого необходимо как разрабатывать новую нормативную документацию по стандартизации, так и проводить совершенствование имеющихся норм и стандартов. Арктический регион является сложным с точки зрения внешних условий. Для Гидротехнических конструкций внешние воздействия в Арктике, в первую очередь, включают воздействия со стороны льда. Ледовые воздействия могут носить самый разный характер и поэтому недостаточно полно отражены в нормативных документах. Образование ледового воротника - это один из типов вмерзания, когда вокруг частично находящегося в воде гидротехнического сооружения, при отрицательных температурах воздуха, образуется утолщение льда, превышающее по своей толщине окружающий ровный лед. Данное утолщение образуется по причине теплового потока через тело сооружения и может быть рассчитано. Аналитический расчет ледового воротника сложно выполним, поэтому предполагается использовать термодинамическую модель для численного моделирования.
Ключевые слова: ледовый воротник, вмерзание в лед, математическая модель, численный расчет, уравнение баланса энергии, метод Стефана.
Б01: 10.37313/1990-5378-2023-25-3-107-113
ЕБ№ ЦЛКЬ
Исследование выполнено в рамках проекта FSEG-2020-0021
ВВЕДЕНИЕ
Возникновение ледовых воротников вокруг цилиндрических опор гидротехнических сооружений при постоянном уровне воды обусловлено изменением направления и интенсивности переноса тепловой энергии из воды в холодный воздух с отрицательной температурой ниже точки замерзания. При использовании материалов с большой теплопроводностью в составе указанных сооружений (сталь и бетон) через них происходит интенсивный ток тепловой энергии в атмосферу. Сооружение нагревается водой ниже поверхности ледового покрова и отдает тепловую энергию в атмосферу над поверхностью ледового покрова. При этом, за счет высокой теплопроводности и отсутствии теплоизоляции в виде ледового покрова интенсивность теплового потока выше, чем на окружающей акватории. В местах интенсивного теплового потока, при переходе че-
Шарапов Дмитрий Андреевич, кандидат наук, доцент Высшей школы гидротехнического и энергетического строительства. E-mail: sharapov@bk.ru Клочков Юрий Сергеевич, доктор технических наук, доцент, проректор по научно-организационной деятельности. E-mail: y.kloch@gmail.com
рез точку замерзания, происходит образование льда. В рассматриваемом случае лед образуется под поверхностью ледового покрова акватории и обволакивает опоры сооружения, образуя характерные «ледовые воротники».
Ледовые воротники опасны по ряду причин. Они могут затруднять проведение операций вблизи сооружений, на которых они образуются, а также являться причиной максимальных пиковых нагрузок на гидротехнические сооружения при подвижках льда.
В общем случае рассматриваются стальные сооружения цилиндрической формы, вертикального профиля. Внешние условия считаются известными. Под внешними условиями понимаются предпосылки для формирования ледяного покрова водоёма и ледяного воротника у сооружения и/или наличие уже образовавшихся ледовых полей и ледовых образований. Предполагается отсутствие приливов и волнения. Считается, что лёд на акватории не движется в момент образования ледяного воротника. Известны: соленость воды, изменение температуры воздуха за заданный период времени - характеризуемое морозо-днями (суммой средних за день отрицательных температур ниже температуры замерзания воды за рассматриваемый период времени).
МАТЕМАТИЧЕСКАЯ ПОСТАНОВКА
В системе вода-лёд-воздух можно выделить три принципиальных тепловых потока: первый со стороны «теплой воды» через лёд в сторону «холодного» воздуха, второй тепловой поток образуется в процессе образования льда и направлен в сторону «холодного» воздуха, третий поток является суммой первых двух и направлен со стороны поверхности лёд-воздух в сторону «холодного» воздуха.
В процессе образования льда выделяется тепловая энергия. Когда первый тепловой поток становится равным третьему, то тепловой энергии в процессе образования льда не выделяется, т.е. нового льда не образуется. Решение задачи базируется на уравнении энергетического баланса.
Основные расчетные зависимости
Уравнение баланса энергии, описывающее температурное распределение внутри льда, воды и движение границы смены фаз (вода -лёд) (Метпапоу, 1992 [1-41):
(1)
рассмотрены ниже); р - плотность материала; с - теплоемкость материала; к - коэффициент теплопроводности материала; причем ¡,м> - индексы показывающие принадлежность параметра ко льду или воде.
Уравнение (1) описывает теплоперенос внутри льда при Т < Т и в воде при Т > Т. Следовательно принимается, что (р, с, к) = (р; , с. , к. ) при Т < Тг и (р, с, к) = (рж , ск, кк ) при Т > Тг. При этом учитывается, что теплоёмкость морской воды и морского льда зависит от температуры и солености уравнения (2), (3), рис. 1, [5-81:
с = 4.1325 + 0.0974Г + 0.0049357"2 +
+ 0.00009!ТЭ(-^-) ,
1.6 -
7.43 565 к J
(2) (3)
""Эх' ду ' 8з
где t - время; х,у,г - горизонтальные и вертикальная координаты; Т - температура; Т - температура замерзания воды; I - теплота образования вещества; 8(Т - Т) - дельта функция; (дельта функция, и её параметры
1.22932 4 Г vHttf=-i7J
Уравнение (3) для определения теплоемкости морского льда выведена для солености 2 ppt. Турбулентное перемешивание морской воды учитывалось дополнительным коэффициентом турбулентной теплопроводности (kwt)■ Уравнение баланса энергии, описывающее температурное распределение внутри опоры записывается аналогично.
Граничные и начальные условия
Температура на границах лед(вода)-воздух и опора-воздух:
Т — Т
■
Температура на нижней границе выделенной области:
с , kJ/kg °С
с,, kJ/kg
4
3.5
3.6 3.4 Ъ2
а) -■
Т.* С
-30 -25 -20 -15 -1 0 5 0
20 15 10
Ь) 4рр«1
2 ppt Т, "С
-20 -1 5 <10 5 0
Рис. 1. Теплоемкость морского рассола (соленая вода) (а) и морского льда (Ь) от температуры; соленость морского льда 2 - 4 рр^ (Доронин и Хесин, 1975 [9])
1—1 bottom-
Тепловой поток внутри выделенной области через границы материалов:
-<ъ) = 0;<?i = ¿¡VT-.
Условия отсутствия теплового потока на внешних боковых границах всей выделенной области:
п ■ (-feVT1) = 0.
Начальное распределение температуры внутри выделенной области:
Г(*в> = ^initio!-где: Г - температура; Tair - температура воздуха; Tbtt - температура нижней границы выделенной области; п - вектор нормальный к границе; q - тепловой поток; к - коэффициент теплопроводности материала; £0 - начальный момент времени; Tiriitial - температура в начальный момент времени.
Численный метод
Для численного интегрирования уравнения (1) использовался персональный компьютер. Рассматривались как трехмерные, так и двухмерные аппроксимации- В случае цилиндрической опоры задача будет являться осесимме-тричной и решение может быть рассмотрено в плоскости г - z, проходящей через центр опоры.
Для математического описания процесса смены фазового состояния используется Дирак-дельта функция [10-14]. Она вводится в уравнение теплового баланса и принимает ненулевые значения на границе смены фаз. Дельта функция позволяет учитывать теплоту образования льда в заданном интервале температуры смены фазового состояния. Дирак-дельта функция аппроксимируется следующим образом (4):
5^-7}) =
(l-tmh(l(r-Tf - е)))(i+tenh (л{г-Tf+s) ))
3s
,(4)
где X и e подбираются в соответствии с условиями:
I
+ ОС
,
I
f(x)S1(x)dx = /(0).
(5)
(6)
График Дирак-дельта функции принимая Тг = Л = 100аС-1, £ = 0.01, пред-
ставлен на рис. 2 (А. МагсЬепко 2008 [3]).
Данная аппроксимация рассматривает интервал температур при фазовом переходе:
Т Е (—1.32ПС, — 1.73аС). (7)
Тэта функция (функция Хэвисайда), 9(х) -используется для математической записи изменения коэффициентов уравнения (1) при изменении фазового состояния. Она записывается в виде:
в(х) = 1 при х > 0; в(х) = 0 при х < 0. (8)
Теплоёмкость связана с энергией через плотность. При расчетах имеет смысл сразу рассматривать произведение плотности на теплоемкость - рс. Если учесть, что теплота образования может быть учтена путём увеличения теплоёмкости в заданном интервале (используя Дельта функцию), то произведение - может быть аппроксимировано следующим образом:
рс^ = рЛ0( 7>-Г) +
Коэффициент теплопроводности - k, может быть записан для области вода-лёд в следующем виде:
Таким образом, при температурах ниже точки замерзания тэта функция в (л) от отрицательного значения обращается в ноль и учитывается только коэффициент для льда к^ в случае же температуры выше точки замерзания - учитывается только коэффициент для воды к
Положение границы смены фаз определяется с точностью, зависящей от ширины (Д) простран-
Рис. 2. График функции ^(г - Tf ) (а) и рс(т, Т{ ) (Ь) принимая Tf = -l.BDC, А = 100DC_1, е = 0.01
ственной области, где Т Е (— 1.82ПС,— 1.78аС). Размер элементов сетки выбирался меньшим, чем Д. Если Д меньше размера элемента конечно-элементной сетки, то эффективная теплота образования уменьшается, вследствие численного распределения функции 31 по всему конечному элементу. При расчете использовались сетки, содержащие около 20000-500000 треугольных элементов в зависимости от величины расчетной области. Площадь элемента варьируется от 0.04 м2 до 7*10-5 м2.
Решение задачи осуществлялось как в трёхмерном случае, так и в двумерном. Трехмерный случай дает наглядное представление происходящего вокруг опоры. Для цилиндрических опор целесообразно использовать двумерное (осесимметричное), так как решение будет содержать все необходимые характеристики при меньших затратах времени и ресурсов.
При компьютерном моделировании решается уравнение (11). Его решение в двумерном случае связанно с переходом к цилиндрическим координатам.
Уравнения теплопроводности в декартовых и цилиндрических координатах записываются в виде:
где уравнение (11) записано в декартовых координатах, а уравнение (12) записано в цилиндрических координатах; р - плотность материала; с - теплоемкость материала; к - коэффициент теплопроводности материала; f - внутреннее теплообразование; Г - температура; t - время; х, у, г - горизонтальные и вертикальная декартовы координаты; г, <р,г- радиус, угол, вертикальная координата в цилиндрической системе координат.
Рассматриваемая постановка задачи обладая
= 0 .
ет осевой симметрией, следовательно:
д<р
Тогда уравнение (12) переписывается в виде:
где обозначения аналогичны приведенным ранее.
Раскрывая производную первого члена правой части уравнения (12) можно получить уравнение (13). Уравнение (13) отличается от двумерного
1 5 „ Т
случая уравнения (11) членом--(к-), кото-
г 5г 5г
рый быстро затухает при увеличении г. Поэтому для не слишком маленьких опор можно пользоваться двумерным решением уравнения (11):
где обозначения аналогичны приведенным ранее.
РЕЗУЛЬТАТЫ МОДЕЛИРОВАНИЯ
Для оценки достоверности результатов численного моделирования были проведены физические эксперименты. Результаты экспериментов использовались для оценки правильности компьютерной модели. На рис. 3 представлены некоторые из результатов физических экспериментов и математического моделирования: Опыт со стальным цилиндром длинной 30 см, диаметром 1 см, погруженным в воду на 20 см, вмораживаемым в течение 46 часов при температуре -10°С (рис. 3 слева). На рис. 3 также представлен пример результата 3Б моделирования и пример результата 2Б моделирования.
Формы ледовых воротников, полученных численным моделированием, соответствуют физическим экспериментам, что видно на рис. 3. Отмечены некоторые случаи, когда численное решение несколько завышает размеры воротников, однако, учитывая, что оно так же немного завышает толщину ровного окружающего льда (рис. 4) отношение этих величин в расчетах и на моделях близко. Опыты для более продолжительного периода показали такую же относительную ошибку результата - не более 5% в наибольшем случае.
Рис. 3. Результаты физического эксперимента (слева) (стальная труба) и примеры результатов математического моделирования (3Б - центр, 2Б - справа)
170 т-----------1' Формула Зубова
150 -130 -
П но -
•Формула Лебедева •Полученное решение (1D) Полученное решение (2D)
4J ш к
^q _. ♦ Экспериментальные данные
^ * ♦ i♦ : „
зо ----д/*-*^
500
1000
1500 МГД
2000
2500
Рис. 4. Сравнение толщин льда полученных численным расчетов с известными измеренными значениями при тех же исходных данных
3000
Рассмотренные случаи с положительной температурой внутри оборы показали возможность удаления ледовых воротников обогревом опоры [15, 16].
Для оценки достоверности абсолютных значений получаемых результатов и настройке расчетной схемы были проведены серии экспериментов по определению толщин термообразованного льда при известных параметрах; результаты сравнивались с известными значениями (рис. 4) [17]. Приведенные экспериментальные точки нанесены по результаты анализа данных из Главной Геофизической Обсерватории [18] и Государственного Гидрологического Института [19].
Как видно из приведённого графика - результаты, полученные по методу конечных элементов (Ряды: «Полученное решение (Ш)» и «Полученное решение (20)»), укладываются в рамки результатов, полученных расчётом по эмпирическим формулам и аналитическим решениям, а также экспериментальных данных.
Во избежание возможных погрешностей, связанных с различиями в масштабе времени между моделью и реальными процессами, размеры ледовых воротников находились в отношении к толщинам окружающего льда.
ВЫВОДЫ
Вмерзание сооружения в лед при постоянном уровне может сопровождаться образованием ледовых воротников. Ледовые воротники об-
разуются ниже уровня воды вокруг сооружений и опор, имеющих большую теплопроводность.
Для получения размеров и формы ледовых воротников была разработана термодинамическая модель и программа для ЭВМ[20]. Термодинамическая модель позволяет провести численную оценку размеров и формы ледового воротника для большого массива исходных данных.
Основными данными, влияющими на рост льда, являются продолжительность вмерзания и градиент температуры. Международным стандартом введено понятие морозоградусодней (freezing degree days), являющееся мерой энергии, затраченной(выделенной) при образовании льда и, тем самым, позволяющей упростить вид представления исходных данных для расчета обмерзания.
Для верификации термодинамической модели использованы физические эксперименты. Физические эксперименты проводились как в полевых условиях, так и в условиях холодной лаборатории. Экспериментальные и численные формы ледовых воротников удовлетворительно совпадают.
СПИСОК ЛИТЕРАТУРЫ
1. Meirmanov A. The Stefan problem // Walter de Gruyter, Berlin. - 1992. - T. 244.
2. Loset S., Marchenko A. Field studies and numerical simulations of ice bustles on vertical piles // Cold Regions Science and Technology. - 2009. - T. 58, № 1-2. - C. 15-28.
3. Marchenko A. Thermodynamic consolidation and melting of sea ice ridges // Cold Regions Science and Technology. - 2008. - T. 52, № 3. - C. 278-301.
4. Sharapov D., Shkhinek K., DelValls T. Â. Ice collars, development and effects // Ocean Engineering. -2016. - T. 115. - C. 189-195.
5. Thermal Conductivity of Sea Ice and Antarctic Permafrost. A thesis submitted for the degree of Doctor of Philosophy. / Daniel J. P.: Wellington, New Zealand, 2004. A thesis submitted for the degree of Doctor of Philosophy.
6. SchwerdtfegerP. The thermal properties of sea ice // J. Glac. - 1963. - T. 4. - C. 789-807.
7. Назинцев, Ю.Л. Некоторые данные к расчету тепловых свойств морского льда // Труды института Арктики и Антарктики / Ю.Л. Назинцев. - 1964. -T. 267. - C. 31-47.
8. Шарапов, Д.А. Гидротехническое строительство. Лед и его свойства для сооружений Арктики: учеб. пособие / Д. А. Шарапов, А.С. Большев - СПб: ПОЛИТЕХ-ПРЕСС, 2023. - 89 с.
9. Доронин, Ю.П. Морской лед. / Ю.П. Доронин, Д.Е. Хейсин: Л.: Гидрометеоиздат, 1975.
10. Ram P.K. The Dirac Delta Function and Delta Sequences // Mathematics in Science and Engineering / Ram P. K.Elsevier, 1983. - C. 1-19.
11. Salva N., Tarzia D. Explicit solution for a Stefan problem with variable latent heat and constant heat flux boundary conditions // Journal of Mathematical Analysis and Applications. - 2011. - T. 379, № 1. -C. 240-244.
12. Mccue S., Wu B., Hill J. Classical two-phase Stefan problem for spheres // Proceedings of the Royal Society a-Mathematical Physical and Engineering Sciences. - 2008. - T. 464, № 2096. - C. 2055-2076.
13. Asaithambi A. Numerical solution of Stefan problems using automatic differentiation // Applied Mathematics and Computation. - 2007. - T. 189, № 1. - C. 943-948.
14. Stefan J. Uber die theorie der eisbildung, insbesondere uber die eisbildung in polarmeere // Ann. Phys. -1981. - T. 42 (2). - C. 269-186.
15. ICEMELT, Свидетельство о регистрации программы для ЭВМ RU 2015617335, 07.07.2015. Заявка № 2015614290 от 19.05.2015. / Шарапов Д. А., Шхинек К. Н., 2015.
16. Шарапов Д. А. Численная оценка необходимой энергии для предотвращения образования воротников на морских причалах Арктики / Д.А. Шарапов, А.С. Большев // Научно-технический сборник Российского морского регистра судоходства. - 2016. - T. 44-45. - C. 20-25.
17. Sharapov D., Shkhinek K. Numerical calculation of the ice grow and empirical calculation results // Advanced Materials Research, Proceedings of 3rd International Conference on Materials and Products Manufacturing Technology (ICMPMT 2013), September 25-26, Changsha, China. - 2013.
18. Архивные данные по температуре // Справочник по климату СССР, Омская и Тюменская области за период с 1955 по 1976 год, Главная Геофизическая Обсерватория. - 1976.
19. Архивные данные по толщине льда // Гидрологический ежегодник станции Новый порт бассейна Карского моря за период с 1955 по 1976 год, Государственный Гидрологический Институт. - 1976.
20. Шарапов, Д.А Гидротехническое строительство - ледовые воротники, свидетельство о регистрации программы для ЭВМ RU 2023615858, 20.03.2023. Заявка № 2023613926 ОТ 22.02.2023. / Д. А. Шарапов, 2023.
THERMODYNAMIC MODEL OF ICE COLLAR FOR HYDROTECHNICAL CONSTRUCTION
© 2023 D.A. Sharapov, Yu.S. Klochkov
St. Petersburg Polytechnic University of Peter the Great (National Research University), St. Petersburg, Russia
Arctic development is one of the priorities for the country. For this, it is necessary both to develop new normative documentation on standardization and to improve existing norms and standards. The Arctic region is complex in terms of weather and environmental conditions. Environmental loads and actions on hydraulic structures are most important and primarily include impacts from ice. Ice impacts can be of a very different nature and therefore are not fully reflected in the normative documents. The formation of an ice collar is one of the types of freezing on the submerged hydraulic structures. Ice collar is a thicker ice around structure due to the high thermal conduction through the structure. This ice is formed due to the heat flow through the body of the structure and can be calculated. Analytical calculation of the ice collar is difficult to perform, so it is proposed to use a thermodynamic model for numerical simulation. Keywords: ice collar, freezing into ice, mathematical modelling, numerical calculation, energy balance equation, Stefan method. DOI: 10.37313/1990-5378-2023-25-3-107-113 EDN: LJAISL
REFERENCES
1. Meirmanov A. The Stefan problem // Walter de Gruyter, Berlin. - 1992. - T. 244.
2. Loset S., Marchenko A. Field studies and numerical simulations of ice bustles on vertical piles // Cold Regions Science and Technology. - 2009. - T. 58, № 1-2. - C. 15-28.
3. Marchenko A. Thermodynamic consolidation and melting of sea ice ridges // Cold Regions Science and Technology. - 2008. - T. 52, № 3. - C. 278-301.
4. Sharapov D., Shkhinek K., DelValls T. A. Ice collars, development and effects // Ocean Engineering. -2016. - T. 115. - C. 189-195.
5. Thermal Conductivity of Sea Ice and Antarctic Permafrost. A thesis submitted for the degree of Doctor of Philosophy. / Daniel J. P.: Wellington, New Zealand, 2004. A thesis submitted for the degree of Doctor of Philosophy.
6. SchwerdtfegerP. The thermal properties of sea ice // J. Glac. - 1963. - T. 4. - C. 789-807.
7. Nazincev, Yu.L. Nekotorye dannye k raschetu teplovyh svojstv morskogo l'da // Trudy instituta Arktiki i Antarktiki / Yu.L. Nazincev. - 1964. - T. 267. - C. 31-47.
8. Sharapov, D.A. Gidrotekhnicheskoe stroitel'stvo. Led i ego svojstva dlya sooruzhenij Arktiki: ucheb. posobie / D. A. Sharapov, A.S. Bol'shev - SPb: POLITEKH-PRESS, 2023. - 89 s.
9. Doronin, Yu.P. Morskoj led. / YU.P. Doronin, D.E. Hejsin: L.: Gidrometeoizdat, 1975.
10. Ram P.K. The Dirac Delta Function and Delta Sequences // Mathematics in Science and Engineering / Ram P. K.Elsevier, 1983. - C. 1-19.
11. Salva N., Tarzia D. Explicit solution for a Stefan problem with variable latent heat and constant heat flux boundary conditions // Journal of Mathematical Analysis and Applications. - 2011. - T. 379, № 1. - C. 240-244.
12. Mccue S., Wu B., Hill J. Classical two-phase Stefan
problem for spheres // Proceedings of the Royal Society a-Mathematical Physical and Engineering Sciences. - 2008. - T. 464, № 2096. - C. 2055-2076.
13. Asaithambi A. Numerical solution of Stefan problems using automatic differentiation // Applied Mathematics and Computation. - 2007. - T. 189, № 1.
- C. 943-948.
14. Stefan J. Uber die theorie der eisbildung, insbesondere uber die eisbildung in polarmeere // Ann. Phys. -1981. - T. 42 (2). - C. 269-186.
15. ICEMELT, Svidetel'stvo o registracii programmy dlya EVM RU 2015617335, 07.07.2015. Zayavka № 2015614290 ot 19.05.2015. / Sharapov D. A., Shkhinek K. N., 2015.
16. Sharapov D.A., Bol'shev A.S. chislennaya ocenka neobhodimoj energii dlya predotvrashcheniya obrazovaniya vorotnikov na morskih prichalah Arktiki / D.A. Sharapov, A.S. Bol'shev // Nauchno-tekhnicheskij sbornik Rossijskogo morskogo registra sudohodstva. - 2016. - T. 44-45. - C. 20-25.
17. Sharapov D., Shkhinek K. Numerical calculation of the ice grow and empirical calculation results // Advanced Materials Research, Proceedings of 3rd International Conference on Materials and Products Manufacturing Technology (ICMPMT 2013), September 25-26, Changsha, China. - 2013.
18. Arhivnye dannye po temperature // Spravochnik po klimatu SSSR, Omskaya i Tyumenskaya oblasti za period s 1955 po 1976 god, Glavnaya Geofizicheskaya Observatoriya. - 1976.
19. Arhivnye dannye po tolshchine l'da// Gidrologicheskij ezhegodnik stancii Novyj port bassejna Karskogo morya za period s 1955 po 1976 god, Gosudarstvennyj Gidrologicheskij Institut. - 1976.
20. Sharapov, D.A. Gidrotekhnicheskoe stroitel'stvo
- ledovye vorotniki, svidetel'stvo o registracii programmy dlya EVM RU 2023615858, 20.03.2023. Zayavka № 2023613926 OT 22.02.2023. / D.A. Sharapov, 2023.
Dmitry Sharapov, PhD, Associate Professor of the Higher School of Hydrotechnical and Energy Construction. E-mail: sharapov@bk.ru
Yury Klochkov, Doctor of Technical Sciences, Associate Professor, Vice-Rector for Scientific and Organizational Activities. E-mail: y.kloch@gmail.com
