CC BY
97
УДК 551.345: 536.421
Энтальпийный метод численного решения задач теплопроводности в промерзающих или протаивающих грунтах
Канд. техн. наук Н.А. БУЧ КО СПбГУНТиПТ
The problems of non-stationary heat conductivity with phase transitions, associated with soil freezing, are solved by enthalpy methods; the history of development of the method for the generalized problem of Stefan is presented. An algorithm of calculations of the processes of solidification-melting in moist soils, based on the explicit finite difference approximation of enthalpy equation, being the function of temperature and iciness is proposed. A graphical expression of solutions for different zones of a phase transition in large-dispersed and fine dispersed soils is presented. On the basis of the proposed algorithm a series of programs in different programming languages was developed, as well as the version of criterion form of solution , and by processing of a large number of numerical experiments, by the method of similarity “a generalized method of calculations of seasonal cooling device” has been created.
Техногенная деятельность людей, связанная с разработкой полезных ископаемых в северных регионах страны, нарушает энергетический баланс земной поверхности, приводя к нежелательным явлениям, которые называют общим термином «деградация вечной мерзлоты». Проблема усугубляется наблюдаемым в последние годы глобальным потеплением. Для предотвращения нежелательных последствий деградации приходится искусственно регулировать температурный режим фунтовых массивов, расположенных в основаниях различных сооружений, с помощью различных технических средств. Наиболее целесообразно использовать сезонные охлаждающие устройства (СОУ), работающие за счет естественного зимнего холода атмосферы. Но иногда этого оказывается недостаточно, в связи с чем все чаще рассматривается возможность комбинированного охлаждения: зимой — посредством СОУ, а летом — с помощью холодильных машин. В обоих случаях для выбора типа и конструкции охлаждающих устройств и определения их холодопроизводительно-сти приходится решать задачу определения температурных полей в грунтовых основаниях. В настоящей статье рассматривается математический метод, основанный на достаточно строгих и четких физических представлениях о процессах фазовых переходов во влажных грунтах, который позволяет разработать алгоритм и составить программу прогнозирования температурного режима в системах грунт - охлаждающее устройство — атмосфера.
Энтальпийный метод решения задач нестационарной теплопроводности с фазовыми переходами начал развиваться с середины прошлого столетия как способ, позволяющий распространить метод конечных разностей не только на классические задачи Фурье, но и на задачи нелинейные: с переменными термическими коэффициентами и подвижными границами фазового перехода.
Параллельно развиваются два направления: практическое — решение конкретных инженерных задач и теоретическое — доказательство существования и единственности решения, определение критериев устойчивости, сходимости и погрешности расчетного процесса, а также возможности распространить метод на более широкий класс задач, называемый обобщенной задачей Стефана.
Сейчас трудно указать, кто первым предложил энталь-пийную формулировку задач Стефана. В ее современной форме она складывалась постепенно. Так, Дюзинбер (1945) [19] предлагает следующий прием при решении методом элементарных балансов одномерной задачи с подвижной границей раздела фаз: начиная с момента времени, в который в данном элементарном объеме достигается температура, равная температуре фазового перехода 7^ необходимо суммировать количества тепла, поступающие в объем в каждый последующий интервал времени Ат, и сопоставлять суммарное количество с теплотой, потребной для фазового перехода в данном объеме:
0Га5= фрД^=ФДК где (р — теплота фазового перехода единицы массы;
Ф — теплота фазового перехода единицы объема;
V— объем.
Весь период, пока поступившее тепло меньше теплоты фазового перехода, в объеме поддерживается постоянная температура 7}. Разумеется, чтобы достаточно точно отследить начало и конец фазового перехода, интервал времени Дт должен быть достаточно малым. Разумеется также, что точность определения координаты расположения границы раздела фаз равна половине пространственного интервала Дх Алгоритм, построенный на таком физически очевидном представлении о процессе фазового перехода, уже, по существу, является энталь-
пийным методом расчета, хотя в работах Дюзинбера понятие «энтальпия» (теплосодержание) не упоминается.
Теоретическое обоснование энтальпийного метода, по нашему мнению, принадлежит С.Л.Каменомостской (1958) [8]. Именно она впервые ввела понятие обобщенной задачи Стефана и доказала существование и единственность решения данного класса задач. Ею также доказано, что классическая задача Стефана является частным случаем обобщенной.
По Каменомостской, дифференциальное уравнение теплопроводности в энтальпийной формулировке записывается следующим образом:
ЭЛ/Эг = сНу(£гас1(м)), (1)
гдеи(Т)= £ Ц/)<#; (2)
И(Т) = £ с(0р(/)Л при Г< 7>; (3)
И(Т) = Ф + с(/)р(/)<* ПРИ т^ Т/.
Из уравнения (3) следует, что при температуре 7} функция И( 7) имеет разрыв, однако в работе [8] указано, что функция И(Т) — не что иное, как теплосодержание, т.е. по физическому смыслу является монотонной функцией, зависящей от величины теплового потока к элементам с температурой Т= 7}в данный момент времени.
В термодинамике, как известно, термин «теплосодержание» является синонимом термина «энтальпия». Под ним понимают одну из функций состояния, т.е. это более емкое понятие, чем просто параметр состояния, как температура. При фазовых переходах энтальпия зависит от двух параметров состояния (например, при парообразовании - от Г и х, где х — степень сухости). Функцию же м( 7),определяемую выражением (2), в различных последующих работах называют либо потоком температур, либо функцией теплового потока. Очевидно, что уравнение (1) можно рассматривать и решать независимо от температуры, считая И = И(и), т.е. считая решением (1) поле энтальпий (или поле и), но по значениям И или и можно однозначно определить Т для любого момента времени.
Свои доказательства Каменомостская выполняет, анализируя конечно-разностную аппроксимацию уравнения
(1) по явной схеме:
(/#' -^)/Дт=«,.; + -2и!,)/Ах2 +
(4)
+М*у+1-2«*,)/Ау2-
Здесь уравнение (4) написано в форме, справедливой для двумерной задачи в прямоугольной системе координат, но переход к трехмерной задаче и другим системам координат не вызывает серьезных осложнений.
Из уравнений (1) и (4) также явствуют следующие преимущества энтальпийной формулировки задач теплопроводности: исключаются трудности решения, связанные с переменностью тепловых коэффициентов А.( 7), с( 7), р( 7); отпадает необходимость особо выделять границу раздела фаз, так как ее положение определяется значениями энтальпии в интервале И' и И". Если за начало координат принята температура 7}, то И' = 0, а И" = Ф (см. рисунок).
Характер зависимости функции И(Т,ш) — сплошная линия — и и(Т) — пунктир — от температуры: а — для крупнодисперсных грунтов (вся влага в свободном состоянии); б — для мелкодисперсных грунтов (связанная влага замерзает в интервале отрицательных температур)
В работе [8] показано, что условие устойчивости явной конечно-разностной аппроксимации энтальпийного уравнения(4) при решении обобщенных задач Стефана остается таким же, как и для задач без фазовых переходов. Данная работа не имела целью разработку расчетного метода, но в ней есть указание на то, что предложенная схема может быть использована как основа алгоритмов программ для ЭВМ. Однако в то время вычислительные машины имели сравнительно малый объем оперативной памяти и были недостаточно быстродейственны. На основании того, что метод Каменомостской не позволяет использовать неявную конечно-раз-ностную схему, в теоретических работах О.А. Олейник [ 111 и Л. Рубинштейна и др. его признали малоэффективным, и большинство авторов прикладных работ пошли иными путями.
Оригинальное решение для плит и цилиндров предложено Д.С.Бакстером (1962): решены одномерные задачи на аналоговой ВМ, причем все переменные, включая энтальпию и поток температуры, приведены к безразмерной (критериальной) форме. Знаковой является также работа А.А. Самарского и Б.Д. Моисеенко (1965) [16], в которой для решения многомерной задачи Стефана применен метод прогонки. Этот метод в разных вариациях, в сочетании с неявной схемой конечно-разностной аппроксимации и введением понятия эффективной теплоемкости применен и во многих последующих исследованиях (например, [1, 17, 20, 21]).
Под эффективной теплоемкостью сЭф(7) понимают фиктивную функцию, которая, будучи введена в уравнение теплопроводности, позволяет учесть теплоту фазового перехода не при постоянной температуре а на некотором произвольно выбираемом интервале А 7}. Этот прием иногда называют «размазыванием теплоты фазового перехода по температуре». Уравнение теплопроводности при этом едино для обеих фаз, но сугубо нелинейно из-за максимума («всплеска») на кривой с^Т). Этот всплеск тем выше, чем меньше интервал АТр т.е. чем больше приближается численное решение к модели фазового перехода с четкой границей раздела фаз (классическая задача Стефана). При этом под энтальпией пони-
мается функция одной переменной Т в соответствии с выражением
КО = Г [с(0р(0 + ФЬ(Т- Tf )]dt, (5)
V
где 5(х) — дельта-функция Дирака.
Производя замену переменных, вместо (1) получаем уравнение
[Сэф(и)Л(и)][ды/Эт) = div(grad(w)). (6)
Далее (6) аппроксимируется по неявной схеме:
-и* ;)/Дг=«, + utXj-2и‘:/)/Лх! +
Ци) (7)
и решается методом прогонки [16]. Расчеты по описанной модели называют методом сглаживания коэффициентов. Данная модель дает тем лучшие результаты, чем лучше параметры сглаживания соответствуют конкретным параметрам пространственно-временной сетки и другим условиям задачи.
Здесь хотелось бы отметить, что, хотя решения по обеим схемам конечно-разностной аппроксимации [уравнения (4) и (7)] называют энтальпийными, собственно энтальпия присутствует и, следовательно, вычисляется только по уравнению (4), т.е. по явной схеме, которая не требует введения ПОНЯТИЯ Сф.
Основным преимуществом неявных схем по сравнению явными является их абсолютная устойчивость, т.е. возможность выбирать величину шага по времени на порядок большей, чем это диктуется условием устойчивости явных схем. Это дает хорошие результаты при решении задач теплопроводности без подвижных границ раздела фаз: скорость расчетов увеличивается в 4 — 6 раз, что было особенно важно при малом быстродействии вычислительных машин.
Но уже тогда некоторые исследователи подвергали сомнению преимущество неявных схем над явными при решении задач Стефана [2, 12, 15]. Так, в работе [2], посвященной процессам замораживания грунтов, отмечаются следующие недостатки неявных методов, связанные с введением понятия эффективной теплоемкости:
V фазовый переход содержащейся в грунте свободной ВОДЫ происходит при ПОСТОЯННОЙ температуре Tj— const, поэтому размазывание теплоты фазового перехода по температуре тем больше искажает реальную картину поля температур, чем больше свободной воды содержится в грунте и чем шире произвольно выбранный интервал температур АТ/
Vфазовый переход связанной воды происходит при Т< Tf, тогда как всплеск сэф симметричен относительно Тр т.е. опять-таки искажение тем больше, чем шире АТ/ Vалгоритм программы должен включать подпрограмму итераций, устанавливающую соответствие между и текущей температурой на каждом временном шаге;
Vувеличение шага по времени Дт (что составляет преимущество неявных схем) может привести к тому, что приращение температуры элементарного объема за этот промежуток времени оказывается соизмеримым или даже превышает ATf. В этом случае расчетная схема мо-
жет «проскочить» всплеск С,ф или учесть его не полностью, что соответствует неполному учету теплоты фазового перехода. Между тем теоретически обоснованного метода оценки погрешности для такого случая не существует. В большинстве прикладных работ без каких-либо обоснований рекомендуется так выбирать параметры схемы, чтобы обеспечить одновременное попадание 4-8 узлов пространственной сетки в интервал Д 7^. В процессе счета необходимо постоянно проверять это условие и при его нарушении корректировать параметры пространственно-временной сетки.
Перечисленное настолько усложняет алгоритм и увеличивает число расчетных операций неявных схем, что обе схемы оказываются практически равноценными, а иногда явная схема даже предпочтительнее.
Из этих соображений алгоритм расчета в [2] базируется на явной конечно-разностной аппроксимации (4). Кроме того, под энтальпией понимается функция двух параметров. Поскольку нас интересуют процессы затвердевания-плавления во влажных грунтах, то в качестве этих параметров выбраны температура Т и степень льди-стости со, т.е. h = h{ Т,со) и Ah = (Э/г/Э Т)ш АТ + (ЭА/Эсо) т Дсо, (8)
а выражение (3) трансформируется в следующее трехстрочное выражение:
А(7\со) = f c{t)p(t)dt при Т< 7}, w = 1 = const;
А( 7',со) = А1 + Ф( 1 — со) при Т= Tf = const, 1 > со > 0; (9)
h(T,со) = A" + j7 c{t)p(t)dt при Т> Tf, со = 0 = const,
где принято со = 0 для талого грунта и со — 1 для твердомерзлого;
А1 — энтальпия при Т= Tj, со = 1;
А11 — энтальпия при Т— 7}, со = 0.
Переход от функций A = А( Т,со) и и( 7) к температурам и степеням льдистости (или обратно) для различных зон фазового перехода наглядно показан на рисунке.
При этом следует учитывать, что в зонах твердомерзлого и полностью талого грунта эти зависимости практически линейны, так как коэффициенты теплопроводности >.м и А,т, так же как и теплоемкости (рс)м и (рс)т, по опытным данным, практически не зависят от температуры (индексы «м» и «т» относятся к мерзлому и талому состояниям соответственно). В зоне фазового перехода свободной воды понятия с и X вырождаются. В зоне фазового перехода связанной воды (от Тм до Tj) экспериментальных данных пока нет (и их будет трудно получить), но с большой долей вероятности можно полагать, что коэффициенты X и (рс) изменяются в пределах указанных значений пропорционально степени льдистости со.
Более прост случай, когда вся влага в грунте является свободной (см. рис. а), но очевидно, что усложнения алгоритма, связанные с наличием в фунте зоны промерзания связанной воды, в которой одновременно изменяются Ти со (см. рис. б), принципиальных трудностей не вызывают. Для этого необходимо только иметь зависимость между 7’и со в указанной зоне. Чаще всего это экс-
периментальная кривая льдистости исследуемого грунта, из которой следует найти производную с1Т/с1<л как некоторую функцию температуры Р( 7). Учитывая, что на участке неизотермического фазового перехода изменение энтальпии происходит с изменением как Т, так и ю (т.е. ДИ = Дйт + Д/гш, где \ИТ = (рс)мД7', а Дйш = Ф Дсо), путем несложных преобразований получаем
До» = ДА/[(рс)м Р( Т) + Ф], Д7"= АИ Р{ 7)/[ (рс)м Д 7) + Ф|.
Создание алгоритма по явной конечно-разностной схеме (4) с учетом (8) и (9) оказалось возможным благодаря разработкам и рекомендациям известного математика, профессора СПбГУНТиПТ П.П.Юшкова [18].
На основе описанного алгоритма был разработан ряд программ на различных языках программирования: Алгол, различные версии Фортрана, в том числе 190. Точность расчетов проверялась сопоставлением результатов численных и аналитических решений (в частности классического решения Стефана).
Программы предназначались для выбора оптимальных технических решений при проектировании конкретных сооружений (мерзлотная завеса в Анадырской грунтовой плотине), для прогнозирования температурных полей в грунте при использовании технически различных способов замораживания: с помощью СОУ различных типов, а также жидким хладоносителем (рассолом), жидким азотом и др. [3 — 5]. Результаты расчетов по программам неоднократно сопоставлялись с результатами наблюдений на специальных полигонах и непосредственно на объектах строительства. Кроме того, был разработан вариант программы с безразмерными переменными (критериальная форма решения) и на этой основе путем обработки большого числа численных экспериментов методом подобия создан «Обобщенный метод расчета СОУ» [6|.
Энтальпийный метод численного решения задач теплопроводности с подвижными границами раздела фаз по явной схеме конечно-разностной аппроксимации продолжает развиваться и находит сторонников среди известных специалистов по замораживанию грунтов [17, 18]. Суть этого метода (с усовершенствованиями) очень точно описана и рекомендована к применению в учебнике по геокриологии МГУ [10]. Однако в некоторых работах можно найти и иную оценку, основанную на явном недоразумении, которое нуждается в разъяснении. Начало было положено А.А. Плотниковым [13], который решал задачу о промерзании грунта энтальпийным методом по неявной схеме (т.е. С применением ПОНЯТИЯ Сзф) [15], и, видимо, не уяснив четко различия между явной и неявной схемами, необоснованно критикует авторов работ [ 19] и
[2] за то, что они не вводят понятия с^, из-за чего якобы неточно учитывают теплоту фазового перехода свободной влаги в грунте. Но из изложенного выше очевидно, что решение по явной схеме в этом понятии не нуждается. К сожалению, Я.А. Кроник [9], ссылаясь на А.А. Плотникова, повторяет его нелепое утверждение.
Данная работа имеет целью, во-первых, прервать описанную порочную цепочку, а во-вторых, обратить внимание молодых исследователей на простой и надежный
метод решения тепловых задач по замораживанию-оттаиванию грунтов и других материалов, сравнительно легко реализуемый на современных компьютерах.
Список литературы
1. Будак Б.М., Соловьева Е.Н., Успенский А.Б. Разностный метод со сглаживанием коэффициентов для решения задач Стефана //ЖВМиМФ. Т.5. 1965. № 5.
2. Бучко Н.А. Алгоритм численного решения двухмерной задачи Стефана энтальпийным методом по трехслойной явной схеме // Холодильная и криогенная техника и технология. — М., 1975.
3. Бучко Н.А., Кузнецов А.Л., Гапеев С. И. Сезоннодействующие охлаждающие устройства и их использование в северном строительстве // Труды Третьей Международной конференции по мерзлотоведению. Т.1, - Канада, Эдмонтон, 1978.
4. Бучко Н.А. и др. Применение жидкого азота на строительстве Ленинградского метрополитена // Транспортное строительство. 1976. № 2.
5. Бучко Н.А., Турчина В.А. Искусственное замораживание фунтов. — М.: Информэнерго. 1978.
6. Бучко Н.А. Система критериев и обобщенные зависимости для расчета процессов замораживания грунтов с помощью сезонных охлаждающих устройств // Холодильная техника. 1978. № 1.
7. Вялое С.С. и др. Термосваи в строительстве на Севере. —Л.: Стройиздат, 1984.
8. Каменомостская С.Л. О задаче Стефана // Математический сборник. Т.53. 1961. № 4.
9. Кроник Я.А.,Демин И.И. Расчеты температурных полей и напряженно-деформационных состояний МКЭ. -М.: МИСИ, 1982.
10. Общая геокриология. Учебник под ред. Э.Д. Ершова. 4.5 Инженерная геокриология/ Г.П. Пустовойт, Л.Н. Хрусталев, главы 1,4. — М.: МГУ, 1999.
И. Олейник О.А. Об одном методе решения общей задачи Стефана //ДАН СССР Т. 135.1960. № 5.
12. Палькин Ю.С. К расчету на ЭВМ двухмерного нестационарного процесса теплопроводности в промерзающем фунте // Известия вузов. Сфоительство и архитектура. 1970. № 3.
13. Плотников А.А. Расчет температурного режима вечномерзлых оснований // Энергетическое строительство. 1978. № 8.
14. Регулирование температуры фунтов основания с помощью сезоннодействующих охлаждающих устройств,— Якутск: СО АН СССР, 1983.
15. Робертсон С.Р., Шенк Н. Исправление уравнения Лондона и Себана в случае отвердевания расплавленного металла // Transactions of ASME, Теплопередача. 1967. № 1.
16. Самарский А.А., Моисеенко Б.Д. Экономичная схема сквозного счета для многомерной задачи Стефана // ЖВМиМФ. Т.5. 1965. №5.
17. Шамсундар Н., Спэрроу Е.М. Применение метода энтальпий к анализу многомерной задачи теплопроводности при наличии фазового перехода // Transactions of ASME, Теплопередача. 1975. № 3.
18.Юшков П.П. О численном интефировании уравнения теплопроводности в случае, когда термические коэффициенты зависят от температуры // ИФЖ. Т.1. Минск. 1958. № 9.
19. Dusinberre G.M. Numerical Method for Transient Heat Flow // Trans.of ASME,v.67,1945.
20. Hashemi H. Т., Sliepcevich C.M. A Numerical Method for solving two-dimension Problems of Heat Conduction with Change of Phase // Chemical Engineering Progress, Symhosium Series, v.63,1967, № 79.
21. Volter V., Cross M. Accurate Solutions of Moving Boundary Problems Using the Enthalpy Method. -Int. J. Heat Mass Transfer, v.24, 1981.
