CC BY
25
2018 ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА Т. 14. Вып. 3
ПРИКЛАДНАЯ МАТЕМАТИКА. ИНФОРМАТИКА. ПРОЦЕССЫ УПРАВЛЕНИЯ
ПРИКЛАДНАЯ МАТЕМАТИКА
УДК 532.517+532.542 МБС 80А22
О расчете оледенения поверхностей в морской воде
Г. И. Курбатова
Санкт-Петербургский государственный университет, Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7—9
Для цитирования: Курбатова Г. И. О расчете оледенения поверхностей в морской воде // Вестник Санкт-Петербургского университета. Прикладная математика. Информатика. Процессы управления. 2018. Т. 14. Вып. 3. С. 186-199. https://doi.org/10.21638/11702/ spbu10.2018.301
Работа посвящена выбору метода численного решения задачи Стефана, который позволял бы рассчитать толщину слоя намерзающего льда в морской воде с приемлемой точностью в течение длительного промежутка времени и допускал обобщение на многомерные задачи оледенения. Сравнение точности расчета толщины слоя льда разными вариантами метода сквозного счета проведено на примере классической двухфазной задачи Стефана, имеющей точное аналитическое решение. Приведены результаты численного решения этой задачи методом сквозного счета по неявной схеме при постоянных шагах по пространству, по времени и при постоянной величине интервала размазывания дельта-функции в обобщенном уравнении теплопроводности. Из сравнения численного и аналитического решений определена зависимость суммарной погрешности численного расчета толщины слоя льда от величины интервала размазывания. Приведен качественный анализ причины нарастания этой погрешности со временем при любом выборе постоянного интервала размазывания. Предложено условие, при выполнении которого можно обеспечить в численном решении методом сквозного счета приемлемую точность вычисления толщины слоя льда для всех моментов времени. Рассмотрен вариант метода сквозного счета с переменным интервалом размазывания дельта-функции, в котором такое условие выполняется. Продемонстрирована возможность рассчитать этим методом толщину слоя льда в морской воде в течение длительного промежутка времени с приемлемой точностью.
Ключевые слова: задача Стефана, оледенение, газопровод, варианты метода сквозного счета, численные решения, примеры расчета.
Введение. В северных морях возможно оледенение внешней поверхности газопровода. При проектировании морских газопроводов важна информация о тол© Санкт-Петербургский государственный университет, 2018
щине слоя намерзающего льда, зависящей от температуры газа в потоке, от теплофизических и конструктивных параметров слоев обшивки, от температуры и солености морской воды, а также от условий прокладки газопровода.
Для практических целей часто важнее знать толщину слоя льда, чем распределение в нем температуры. Известно, что с достаточной точностью рассчитать толщину слоя льда позволяют численные методы, явно выделяющие границу фазового перехода [1]. Один из таких методов использовался авторами, например, в [2]. Однако методы с выделением фронта трудно реализуемы в многомерных задачах. Если условия прокладки газопровода таковы, что важен учет отсутствия аксиальной симметрии процесса теплообмена потока с окружающей средой, то квазиодномерная модель транспортировки газа становится недостаточно точной. В этих случаях для расчета теплообмена необходимо решать многомерную задачу Стефана.
Целью настоящей работы является выбор метода решения задачи Стефана, допускающего сравнительно простое обобщение на многомерные варианты и эффективного для решения задач оледенения в морской воде. Как известно, переход к многомерным задачам проще всего реализуется в методах сквозного счета (энтальпийных методах), в которых поверхность раздела фаз явно не выделяется и задача сводится к решению обобщенного уравнения теплопроводности. Эти методы практически одновременно были предложены и обоснованы в 1965 г. А. А. Самарским и Б. Д. Мои-сеенко [3], Б. М. Будаком, Е. Н. Соловьевой и А. Б. Успенским [4]. В 60-х годах С. Л. Каменомостская ввела понятие обобщенного решения для многомерной задачи Стефана. Уравнения теплопроводности для различных фаз трактовались как одно уравнение для энтальпии, являющейся известной функцией температуры и терпящей разрыв первого рода при значении температуры, равном температуре фазового перехода. О. А. Олейник и С. Л. Каменомостской были доказаны теоремы существования и единственности обобщенного решения многомерной задачи Стефана. Для одномерных задач энтальпийный метод был рассмотрен ранее в работе [5], однако в ней отсутствовало доказательство того, что энтальпийная формулировка задачи Стефана эквивалентна классической. В настоящее время известно множество вариантов метода сквозного счета.
Постановка задачи. Для корректного сравнения точности расчета толщины слоя льда разными вариантами метода сквозного счета в настоящей работе выбрана одна из немногих двухфазных задач Стефана, имеющая точное аналитическое решение. На примере этой задачи в исследуемом диапазоне условий рассмотрены возможности разных вариантов метода сквозного счета по расчету толщины слоя намерзающего льда.
Введем следующие обозначения: Тж — температура морской воды (в К), Т* — температура фазового перехода морская вода—лед (в К), Т0 — температура на заданной поверхности (например, на внутренней поверхности газопровода) (в К), к — номер фазы (к = 1 соответствует льду, к = 2 — воде), ей — коэффициент удельной объемной теплоемкости к-й фазы (в Дж/(м3-К)), — коэффициент теплопроводности к-й фазы (в Вт/(м • К)), Q — удельная (на единицу объема льда) теплота фазового перехода (скрытая теплота плавления) (в Дж/(м3)).
В работе авторов [6] был предложен алгоритм выбора согласованного набора теплофизических параметров нарастающего морского льда. На основании этих данных в расчетах использовался следующий набор теплофизических параметров и условий (все величины указаны в вышеприведенных единицах):
Тж = 272.15, Т* = 271.236, Т0 = 265.236, С1 = 2150500, е2 = 3995450,
А1 = 2.0, Л2 = 0.6, Q = 744584912.5.
Для исследуемых задач оледенения в соленой морской воде характерны малые значения числа Стефана. Для набора параметров (*) число Стефана равно 0.0173.
Классическая одномерная двухфазная задача Стефана, на которой опробовались разные варианты метода сквозного счета, записывается так [7]:
с 1-
дТг
Ж
А1
д2Тх дх2
г> 0, х € (0, у(*));
г> 0, х = 0 : Т1 = Т0;
дТ2 Л д2Т2 , , ,
= ¿>0, хе(у(г), оо);
t = 0, х € (0, те) : Т2(х) = Тто;
0, х ^те : Т2 = Тто; г> 0, х = у(г): Т1 = Т2 = Т*;
А1
дТг
дх
—
у—0
дП
дх
Q
у+0
t = 0: у = 0,
(1) (2)
(3)
(4)
(5)
(6)
(7)
(8)
здесь Т^(х,€) — температура к-й фазы, у(€) — толщина слоя льда (координата границы фазового перехода), (1), (3) — одномерные уравнения теплопроводности, (2) — граничное условие для температуры льда на плоскости х = 0, (4) — начальное условие для температуры воды (предполагается, что в начальный момент времени фаза льда отсутствует и вся вода находится при температуре Тж), (5) — граничное условие для температуры воды на бесконечности, (6) — условие неизменности температуры на границе фазового перехода, (7) — уравнение баланса между потоками тепла на границе фазового перехода и количеством тепла, выделяющимся при образовании льда ((6), (7) — условия Стефана), (8) — начальное условие для толщины слоя льда.
Для рассматриваемых задач характерно выполнение условия
То <П< Т
(9)
при котором температура как льда, так и воды является монотонно возрастающей функцией по пространству и монотонно убывающей функцией по времени; у(Ь) — монотонно возрастающая функция времени.
Задача (1)-(8) имеет аналитическое решение [7]:
Т1(х,г) = Л1 + В1 ей
\2а1лД
Т2(х^) = Л2 + В2 ей
у(1) = а^Т. (10)
Постоянные а^ определяются по формуле а^ = (А&/с^)1/2, к = 1,2; постоянные Л1,Л2,В1,В2 известным образом [7] выражаются через параметры задачи (1)-(8) и величину а, которая рассчитывается по трансцендентному уравнению [7] и зависит от всех параметров задачи: Т0, с1,с2, А1, Л2, Q, Т*,ТЖ.
Рассмотрим, следуя работе [3], алгоритм численного решения задачи Стефана (1)-(8) методом сквозного счета. Для постоянных теплофизических коэффициентов
х
х
С1, е2, Л1, Л2 решение задачи в области х € (0, то) при граничных условиях (2), (5) и начальном условии (4) сводится к решению обобщенного уравнения теплопроводности
дТ д ( дТ
е(Т ) = Л(Т ) =
е1 при Т < Т*,
е2 при Т > Т*,
Л1 при Т <Т*,
Л2 при Т > Т*,
(11) (12) (13)
где 6(£) — дельта-функция Дирака.
Решение Т(х,Ь) уравнения (11) при условиях (12), (13), как показано в работе [3] для общего случая многомерных многофазных задач Стефана, удовлетворяет уравнениям теплопроводности (1), (3) вне границы фазового перехода. На границе фазового перехода поток тепла и энтальпия являются разрывными функциями, поэтому выполнение условия Стефана (7) доказывается в интегральном смысле. Доказательство основано на свойствах дельта-функции и на условии на подвижной границе х = у(Ь) фазового перехода. Это условие следует из равенства (6) и в одномерном варианте записывается таким образом:
= у(;)
¿Т = 0
дх <И дЬ
V дх )
у№
(К'
(14)
При переходе к численному решению уравнения (11) дельта-функция заменяется дельтаобразной (размазанной) дельта-функцией с величиной полуинтервала размазывания Д. Если в каждой из фаз коэффициенты С1 и С2 можно считать постоянными, простейшей интерполяцией скачка энтальпии при переходе фазовой границы является линейная, ей соответствует следующее представление эффективной теплоемкости {е(Т) + QS(T - Т*)} = е(Т, Д):
е(Т, Д)
е1 , С1 + С2
Т <Т* -Д,
Q
2
е2,
+ 2К' - А < Т < Т, + Л, Т>п + Д.
(15)
На этом же интервале температуры эффективный коэффициент теплопроводности полагается равным Л(Т) = Л(Т, Д),
Л(Т, Д) =
Л1, А1 +Л2 2 :
Л2,
Т <Т* -Д,
Т* -Д <Т < Т* + Д,
Т>п + Д.
(16)
Решение уравнения (11) заменяется решением уравнения со сглаженными коэффициентами
—>
—>
—>
дТ д дТ
с<т'А»Ж = а;(л<т'л>а?) <17>
с соответствующими начальными и граничными условиями для температуры Т(х,Ь, Д).
Вопрос о сходимости решения уравнения (17) к решению уравнения (11) при Д —> 0, как отмечается в работе [3], рассмотрен в статье [8].
Численное решение. Уравнение (17) для набора параметров (*) решалось методом сквозного счета по чисто неявной схеме, в которой сглаженные коэффициенты е(Т, Д) и Л(Т, Д) рассчитывались на предыдущем шаге по времени по формулам (15), (16). Метод сквозного счета, предложенный в работе [3], не дает ответа на вопрос о выборе величины интервала Д размазывания дельта-функции. Несмотря на множество работ по исследованию решения обобщенной задачи Стефана методом сквозного счета, выбор величины интервала Д до конца неясен.
В настоящей работе при постоянных значениях безразмерных шагов Н и т (по х и по ; соответственно) в результате вычислительного эксперимента находилась величина Д, при которой сумма всех отклонений численного значения толщины слоя льда от его аналитического значения является минимальной на исследуемом временном интервале.
Первый вывод, к которому привели расчеты, состоял в следующем: численное значение границы льда уп ближе всего к ее аналитической величине у(Ьп) для интервала времени (0,Ьиоп), равного 23.22 ч, если потребовать, чтобы интервал температуры от Т* - Д до Т* + Д для любого момента времени Ьп из интервала (0,Ьиоп) содержал не более двух узлов по х. Для выполнения этого требования при выбранных Н и т находилась область допустимых значений Д, для которых в каждый момент времени Ьп существует узел гп такой, что выполняются условия
тпп_1 < Т* -Д, Тп+1 >Т* + Д, Т"+1 - ТПП_1 > 2Д, у;п е (0,;коп), (18)
из которых следует, что
хгп-1 <уп <хгп + 1. (19)
Толщина слоя льда полагалась равной уп = .
Если не требовать выполнения неравенств (18) и допускать попадание в указанный интервал температуры трех и более узлов по х, положение границы льда уп определяется по найденному массиву температуры как координата х^ узла, в котором температура Т™ наиболее близка температуре фазового перехода Т*. Однако, как отмечалось выше, при таком подходе при прочих равных условиях погрешность в вычислении уп оказывается существенно больше, чем при расчете с использованием условий (18).
Среди всех допустимых величин интервала Д, удовлетворяющих (18) при заданных Н и т, выбиралась такая, при которой сумма отклонений уп от аналитического решения у(Ьп) минимальна на заданном интервале времени. Обсуждение результатов выбора Д, а также примеры расчетов приведены далее.
Задача решалась при следующих условиях: задавалась начальная толщина у0 слоя льда (кратная шагу по х), по ней рассчитывались начальный момент времени ¿0 из аналитического решения (10) и начальное распределение температуры Т(0)(х,Ьо) в слое льда. Кроме того, по аналитическому решению (10) находилось начальное распределение температуры Т^0)(х,Ьо) в слое воды длиной /(¿о). Так называемая длина
теплового влияния /(£), начиная с которой для температуры воды с заданной точностью выполняется условие на бесконечности (5), определяется как минимальная длина, отсчитываемая от х = 0, при которой в момент времени £ выполняются следующие условия:
Т2(М) - Та
< £1,
дТ
дх
< £2,
1,г
где £1, £2 — заданные малые величины.
При решении задачи на интервале (¿0,^коп) длина теплового влияния выбиралась максимальной /т, необходимой для выполнения этих неравенств в конечный момент времени Ьиоп. Величина /т находилась из аналитического решения задачи.
Сглаженное уравнение (17) при выбранных начальных условиях решалось численно по алгоритму, основанному на выполнении требований (18). Для разных значений интервала Д температуры, найденные в численном решении величины уп(Д) толщин слоя льда в моменты времени ¿п, сравнивались с точными значениями у(Ьп) = (10) при тех же начальных условиях. Примеры сравнений представ-
лены в табл. 1 для набора параметров (*) при у0 = 0.5 см, ¿0 = 13.2 мин, значениях безразмерных шагов Н = 0.1, т = 1 и при указанных величинах интервала Д температуры. Для набора параметров (*) величина а в аналитическом решении равнялась 0.0001777 м/с1/2. В табл. 1 для разных моментов времени, отсчитываемых от t = 0, приведены размерные значения у(Ьп) = и величин отклонений = у(£п) — уп(Дт), т = 1, 2, 3. На интервале (0,^оп) рассчитывалась суммарная погрешность Б(Д), определенная равенством
N
Б(Д) = Е
у(П — уп(Д)
где N — число шагов по времени. Приведем зависимость суммарной погрешности Б(Д) (в см) от Д (в К) на интервале времени (0,^оп) при ¿иоп = 23.22 ч:
Д • 103 5(Д)
6
420.5
5.7
396.8
5.6 392.4
5.5 390.5
5.4 389.7
5.3 391.9
5
415.9
4
741.2
Из этой зависимости следует, что для принятых параметров задачи суммарная погрешность Б(Д) в численном расчете толщины слоя льда уп(Д) минимальна при Д = 0.0054.
Чтобы продемонстрировать характер поведения во времени отклонений при разных Д, приведем точные значения толщин слоя льда у^п) = и величины отклонений 5т = у(£п) — уп(Дт), т = 1, 2, 3.
Проведенные расчеты, часть из которых представлена в табл. 1, привели к выводу о невозможности для всего интервала времени вычисления толщины слоя льда с приемлемой точностью по указанному варианту метода сквозного счета с постоянными Н, т и Д. Наибольшим недостатком для решения задач о намерзании слоя льда для больших промежутков времени явилась нарастающая со временем погрешность в вычислении уп(Д), которая имела место в расчетах при любом допустимом значении Д. Такое положение дел отмечалось также в ряде работ, например в [9].
Таблица 1. Значения у(Ьп) и отклонений 5т, т = 1, 2, 3
ч у(Г), см (Ах = 0.006) ¿1, см (Д2 = 0.0054) ¿2, см (Дз = 0.005) йз, см
0.55 0.793 0.193 0.193 0.193
0.89 1.004 0.304 0.304 0.304
1.22 1.177 0.377 0.377 0.377
2.52 1.692 0.392 0.492 0.492
4.82 2.340 0.440 0.440 0.540
7.12 2.844 0.344 0.444 0.444
9.42 3.272 0.172 0,372 0.372
11.72 3.649 0.049 0.249 0.349
14.02 3.991 -0.108 0.091 0.191
16.32 4.306 -0.193 0.006 0.106
18.62 4.602 -0.400 -0.200 0.000
20.92 4.876 -0.524 -0.324 -0.124
23.22 5.137 -0.663 -0.463 -0.263
Сделанный вывод остается в силе и при квадратичной аппроксимации дельта-функции в интервале Т* - Д < Т < Т* + Д, и при вычислении сглаженных коэффициентов с(Т, Д) и Л(Т, Д) на текущем шаге по времени с использованием итерационного метода решения нелинейной системы разностных уравнений. Проведенные ранее расчеты оледенения поверхностей в морской воде (см. [10]) свидетельствовали о том, что при применении метода сквозного счета с неизменным интервалом Д большее влияние на точность совпадения с известными решениями оказывает величина Д, чем выбор интерполяционного полинома. Конечно, можно добиться локального улучшения точности расчета за счет выбора полинома более высокого порядка при сглаживании дельта-функции и коэффициентов теплоемкости и теплопроводности, однако и в этом случае не удается избежать нарастания со временем погрешности в вычислении уп(Д) при любом допустимом постоянном значении Д.
Качественный анализ причины нарастающей погрешности в вычислении границы оледенения при постоянных Н, т, Д. Обсудим причину того, что для нестационарных процессов при постоянных шагах Н и т невозможно выбрать постоянное значение Д таким образом, чтобы обеспечить расчет положения границы оледенения методом сквозного счета с приемлемой точностью для всех моментов времени.
Оценим, с какой погрешностью выполняется условие Стефана (7) при численном решении сглаженного уравнения (17). Запишем интегральное уравнение баланса внутренней энергии £ (в Дж/кг) покоящейся среды при отсутствии источников (стоков) и внешних сил. Скорость изменения во времени внутренней энергии среды плотностью р в области О равна в этом случае суммарному потоку тепла через замкнутую поверхность Бо области О. Считая, что поток тепла ] определяется законом Фурье: 2 = -Л grad Т, приходим к интегральному равенству
I ^¡г^и = <!> \gradТ • пёв, (20)
а Sа
где п — внешняя нормаль к поверхности области О.
При сделанных предположениях давление р не зависит от времени, поэтому в уравнении (20) можно перейти от объемной плотности р£ внутренней энергии к объемной плотности ри энтальпии среды:
др£ дри
ре = ри-р - Ж = Ж- (21)
Обозначим объемную плотность ри энтальпии через Н (в Дж/м3). В одномерной задаче интегральное равенство (20) с учетом (21) для интервала [х^ — к, + к] записывается в виде
х^ +к
дН I дТ
дТ дх
(22)
хгп + к ,
При переходе через границу лед—вода энтальпия Н и поток тепла претерпевают разрыв первого рода, величина скачка энтальпии равна Q. При постоянных коэффициентах теплоемкости и постоянных плотностях льда и воды энтальпия Н зависит от температуры Т следующим образом:
Н (Т ) =
сх Т при Т <Т*, С1 Т* + Q + С2 (Т — Т*) при Т>Т*.
(23)
В работе [11] для трехмерного варианта двухфазной задачи Стефана доказано выполнение условия Стефана (7) в интегральном смысле при энтальпийной записи обобщенного уравнения (11).
Перейдем к сглаженным функциям энтальпии Н(Т, А) и температуры Т(х, Ь, А). Легко показать, что в случае линейной интерполяции скачка энтальпии (23) выражение для производной от сглаженной энтальпии записывается в виде
дН
~дГ
<1Я дТ дТ 1Й
с 1 + с2 0_ 2 2 А
дТ
Ж'
(24)
постоянная величина А в (24) выражается через сх, С2, Т*. В задачах оледенения поверхностей в морской воде характерное значение ^д минимум на три порядка больше, чем значения сх, С2 для любого А, удовлетворяющего неравенству А < 1, поэтому слагаемым (сх + С2)/2 в (24) можно пренебречь по сравнению со слагаемым Q/(2 А). Таким образом, равенство (22) для сглаженных функций Н(Т, А), Т(х,Ь, А) можно записать как приближенное равенство
+к
4" \
2Д У дЬ &
Хг„ -к ^
дТ
— Л2 т—
дх
(25)
Производная ^ от сглаженной функции Т(х, Ь, А) на интервале [х^ — к, + к]
дг
непрерывна, поэтому интеграл в левой части (25) равен
х¿„ +к
1Т дТ
2 к, х* е [х^ — к, хп + к].
дТ дТ
, -к
При к ^ 0 из условия (19) следует, что х* записать приближенное равенство
хг« + к
дТ дТ
, -к
2 к.
это позволяет
к
х
к
х
к
х
х
х
п
х
—-
У
У
х
Для сглаженной температуры должно выполняется условие (6): Т(уп,1п, Д) = т*
const Vt™ > 0. Пусть для производной ^
имеет место приближенное равенство
типа соотношения (14) на подвижной границе фазового перехода:
дТ
Производная по х от сглаженной температуры Т(х,Ь, Д) в точке х = уп непрерывна. С учетом (25) для момента времени Ьп приходим к приближенному аналогу условия Стефана (7):
~ (1у дТ
у„" dt tn дх yn
л дт
1 дх
yn-h
дТ дх
yn+h
dt
1 дТ ~К~дх
h
(26)
Рассмотрим правую часть (26). Выясним, какой пространственный интервал НА соответствует интервалу Д сглаживания температуры и энтальпии. Интервал НА определяется равенствами
т(уп,гп, Д) = т* угп > о, т(уп + нА,гп, Д) = т* + Д, т(уп - нА,гп, Д) = т* -Д угп > о.
Из разложения функции т(x,t, Д) в ряд Тейлора в окрестности точки х = уп вытекает
Т(yn + hA,tn, Л) — Т(yn,tn, Д) +
дх
дТ
дх h
hA
Т + Д-П +
ОТ
дх
hA
hA —
,dy
Л
дТ
дх
yn
Обозначим через К множитель при величине О —
dt
в правой части (26):
(27)
К=
I Л дх
h) hA
(28)
Как видно из (26), условие Стефана (7) при Н ^ 0 в численном решении сглаженного уравнения (17) выполняется для всех моментов времени, если выполнено условие
К hA
ytn > 0.
(29)
Проведенный анализ в рамках принятых допущений позволяет объяснить, почему в рассмотренной нестационарной задаче Стефана на больших временах ошибка в расчете уп указанным методом сквозного счета не может быть устранена при любом выборе постоянного интервала Д при постоянных Н и т. Дело в том, что величина
/гд, как следует из (27), непостоянна. При выполнении условия (9) производная ^
y
—-
—-
y
y
y
y
y
снижается со временем, соответственно при постоянных к, т, А значение кд растет, а величина К убывает. Уменьшение со временем множителя К в приближенном равенстве (26) свидетельствует об искажении условия Стефана при численном решении сглаженного уравнения (17) этим методом. За счет выбора А можно добиться выполнения условия К ~ 1 в какой-то момент времени, но с увеличением времени условие (29) будет неизбежно нарушено. Фактически это приводит к следующему: численный расчет уп происходит со все более искаженным условием Стефана, при уменьшении К это приводит в данной задаче к завышению рассчитанного численно значения уп. Действительно, можно показать, что при уменьшении множителя перед производной
—^ в условии Стефана (7) возрастают скорость оледенения и толщина слоя льда, например, для набора параметров (*) при понижении Q в 2 раза величина а (ас нею и у(Ь) = а\/1) увеличивается в 1.404 раза. Именно эту ситуацию демонстрируют данные табл. 1.
Заметим, что можно обеспечить выполнение условия (29) и при постоянном А за счет соответствующего выбора переменного шага по пространству.
Вариант метода сквозного счета, позволяющий с приемлемой точностью рассчитать границу фазового перехода на длительном интервале времени. Были рассмотрены различные варианты метода сквозного счета, допускающие обобщение на многомерные задачи Стефана. Наибольший интерес вызвал вариант, опубликованный в серии работ В. И. Васильева, А. М. Максимова, Е. Е. Петрова, Г. Г. Цыпкина (см., например, [12]). В рамках этого метода можно обеспечить выполнение условия (29) на всем интервале времени, и, как показано в работе [12], метод обобщается на многомерные задачи с фазовыми переходами.
Воспользуемся для решения задачи (11)—(13) вариантом метода сквозного счета, близким к предложенному в [12]. Конечно-разностный аналог задачи для внутренних узлов пространственно-временной сетки строится следующим образом. Разрывные коэффициенты теплоемкости, теплопроводности и дельта-функция Дирака сглаживаются по температуре Т(х, Ь) только в интервалах, содержащих границу фазового перехода, т. е. в интервалах, где выполняются неравенства Т™ ^ Т* < Т"+1. Для всех остальных интервалов дискретным аналогом уравнения (17) служит трехточечное разностное уравнение, аппроксимирующее обычное уравнение теплопроводности с постоянными коэффициентами в соответствующей фазе.
Аппроксимацию коэффициентов эффективной теплоемкости и эффективной теплопроводности запишем следующим образом:
с(ТП,Т:
впсл + (1 - вп)с2 + —-2Ц—,
?+ь (31)
(30)
С2, Тп > Т*, Т+1 >Т*,
< Тп
< 1 ¿+1,
(32)
Уп = (г + вп) к.
Сглаженное уравнение (17) с аппроксимациями (31), (32) решалось численно по неявной схеме. На каждом временном слое нелинейная система разностных уравнений решалась методом простых итераций, на каждой итерации система линейных уравнений решалась методом прогонки. Условием окончания итерационного процесса служило
выполнение неравенства тах
т/+1 — т^ < е* (е* — заданная точность вычисления
температуры, в — номер итерации). При е* = 0.001 К число итераций в расчетах не превышало 3. Остальные параметры, такие как начальные данные, набор (*), длина теплового влияния были такими же, как при решении этой задачи с постоянным Д. Нетрудно видеть, что введенный выше множитель К (28) при аппроксимации (31) остается близким к единице для всех моментов времени. Это дало основание ожидать, что отклонения 6 = у(Ьп) — уп величин уп, вычисленных по данному варианту метода сквозного счета, от величин у(Ьп), найденных из аналитического решения задачи, будут малыми для всех моментов времени. Проведенные расчеты подтвердили ожидание. Отметим, что такой метод требует предварительного анализа выбора шагов т и Н, а именно, необходимо, чтобы за один шаг по времени граница фазового перехода либо оставалась в прежнем интервале, либо переходила в соседний интервал. В табл. 2 представлены рассчитанные значения отклонений 6 = у(Ьп) — уп для разных моментов времени (безразмерные шаги Н и т были такими же, как при расчете данных табл. 1).
Таблица 2. Значения у(1п) и отклонений 5
ч у(гп), см <5, см
0.55 0.793 -0.031
0.89 1.004 -0.041
1.22 1.177 -0.047
2.52 1.692 -0.052
4.82 2.340 -0.053
7.12 2.844 -0.052
9.42 3.272 -0.049
11.72 3.649 -0.048
14.02 3.991 -0.046
16.32 4.306 -0.044
18.62 4.602 -0.0420
20.92 4.876 -0.040
23.22 5.137 -0.038
Сравнение табл. 1 и 2 демонстрирует преимущество метода, в котором обеспечено выполнение условия (29), для расчета положения границы оледенения.
При тех же параметрах задачи и тех же значениях Н и т, для того же интервала времени, равного 23.22 ч, что и при расчете методом с постоянным Д, суммарная погрешность Б в расчете методом, основанном на аппроксимациях (31), (32), составила величину
N
Б = ^ у(гп) — уп = 63.64,
п=1
которая почти на порядок меньше минимальной суммарной погрешности Б(Д = 0.0054) = 389.7 решения этой задачи методом сквозного счета с постоянным интервалом Д. Особенно важно при расчетах длительных процессов оледенения в морской воде то, что метод сквозного счета, основанный на аппроксимациях (31), (32), не приводит к увеличению со временем погрешности расчета границы оледенения.
Замечание 1. Если в аппроксимации (31) коэффициента эффективной теплоемкости заменить слагаемое QAT+k — ТП) на слагаемое (QOn)/(Tn+1 — ТП ), добавив множитель вп (30) (как это сделано для задач, рассмотренных в работе [12]), то расчет при тех же параметрах и шагах h и т приводит в исследуемой задаче к существенной погрешности в расчете границы оледенения yn. Возможная причина этого заключается в том, что при аппроксимации со слагаемым (Q6n)/(Tn+i — ТП) нарушается выполнение условия (29).
Замечание 2. Исследование сходимости описанного варианта метода сквозного счета с интервалом сглаживания А, пропорциональным изменению температуры, выходит за рамки настоящей работы. Отметим следующее: решение классической задачи Стефана (1)-(8), как известно, существует, единственно и совпадает с обобщенным решением. Как показано в настоящей работе, численное решение обобщенной задачи по предложенному варианту метода сквозного счета совпадает с классическим решением с точностью до десятых процента на всем интервале времени. Это свидетельствует о том, что, по крайней мере для решенной задачи Стефана, предложенный вариант метода сквозного счета сходится и притом сходится к классическому решению задачи.
Заключение. В работе рассмотрен выбор численного метода решения задачи Стефана, который допускал бы обобщение на многомерный случай и позволял рассчитать границу оледенения с приемлемой точностью на длительных интервалах времени. Приведены результаты численного решения классической задачи Стефана методом сквозного счета по неявной схеме при постоянных шагах по пространству, по времени и при постоянной величине интервала размазывания дельта-функции в обобщенном уравнении теплопроводности. Из сравнения численного и аналитического решений определена зависимость суммарной погрешности численного расчета толщины слоя льда от величины интервала размазывания. Представлен качественный анализ причины нарастания этой погрешности со временем при любом выборе постоянного интервала размазывания. Предложено условие, при выполнении которого можно обеспечить в численном решении методом сквозного счета приемлемую точность вычисления толщины слоя льда для всех моментов времени. Рассмотрен вариант метода сквозного счета с переменным интервалом размазывания дельта-функции, в котором предложенное условие выполняется. Продемонстрирована возможность рассчитать этим методом толщину слоя льда в морской воде в течение длительного промежутка времени с приемлемой точностью.
Автор благодарит профессора Санкт-Петербургского государственного университета Сергея Константиновича Матвеева за плодотворное обсуждение работы.
Литература
1. Vasilev F. P. On finite difference methods for the solution of Stefan's single-phase problem // U.S.S.R. Comput. Math. Math. Phys. 1963. Vol. 3(5). P. 1175-1191.
2. Kurbatova G. I., Ermolaeva N. N. Analysis of the cylinder glaciation models in seawater // Applied Mathematics and Information Sciences. 2017. Vol. 11(3). P. 925-930.
3. Samarskii A. A., Moiseyenko B. D. An economic continuous calculation scheme for the stefan multidimensional problem // U.S.S.R. Comput. Math. Math. Phys. 1965. Vol. 5(5). P. 43-58.
4. Budak B. M., Sobol'eva E. N., Uspenskii A. B. A difference method with coefficient smoothing for the solution of Stefan problems // U.S.S.R. Comput. Math. Math. Phys. 1965. Vol. 5(5). P. 59-76.
5. Eyres N. R., Hartree D. R., Ingham J. et al. The calculation of variable heat flow in solids // Philosophical Transactions of the Royal Society. Series A. 1946. Vol. 240. P. 1-57.
6. Ермолаева Н. Н., Курбатова Г. И. Нестационарная модель нарастания морского льда //
Вестн. С.-Петерб. ун-та технологии и дизайна. Сер. 1. Естественные и технические науки. 2017. № 1. С. 3-8.
7. Tichonov A. N., Samarskii A. A. Equations of mathematical physics. New York et al., USA: Pergamon Press Ltd., 1963. 781 p.
8. Олейник О. А. Об одном методе решения общей задачи Стефана // Докл. АН СССР. 1960. Т. 135, № 5. С. 1054-1057.
9. Будок Б. М., Васильев Ф. П., Успенский А. Б. Разностные методы решения некоторых краевых задач типа Стефана // Численные методы в газовой динамике. М.: ВЦ МГУ, 1965. С. 139-183.
10. Krivovichev G. V., Ermolaeva N. N., Kurbatova G. I., Miheev S. A. Mathematical modelling of the glaciation process // International Journal of Pure and Applied Mathematics. 2017. Vol. 113(5). P. 609-616.
11. Shamsundar N., Sparrow E. M. Analysis of multidimensional conduction phase change via the enthalpy model // J. Heat Transfer. 1975. Vol. 97(3). P. 333-340.
12. Васильев В. И., Максимов А. М., Петров Е. Е., Цыпкин Г. Г. Тепломассоперенос в промерзающих и протаивающих грунтах. М.: Наука, Физматлит, 1996. 224 с.
Статья поступила в редакцию 30 ноября 2017 г.; принята к печати 14 июня 2018 г.
Контактная информация:
Курбатова Галина Ибрагимовна — д-р физ.-мат. наук, проф.; [email protected]
On the calculation of surfaces glaciation in seawater
G. I. Kurbatova
St. Petersburg State University, 7—9, Universitetskaya nab., St. Petersburg, 199034, Russian Federation
For citation: Kurbatova G. I. On the calculation of surfaces glaciation in seawater. Vestnik of Saint Petersburg University. Applied Mathematics. Computer Science. Control Processes, 2018, vol. 14, iss. 3, pp. 186-199. https://doi.org/10.21638/11702/spbu10.2018.301
Selection of a method for the numerical solution to the Stefan problem is considered, which permits one to calculate the ice boundary with accepted accuracy during a long period of time, and can be extended to multidimensional glaciation problems. A comparison of the different versions of continuous schemes with an exact analytical solution to the classical Stefan problem is presented. The criterion that permits estimation, under given assumptions, of the error behavior of calculating the phase boundary in numerical solution to the Stefan problem using a continuous calculation scheme is suggested. The advantage of this continuous calculation scheme for phase boundary calculation in problems of surface glaciation in seawater is demonstrated.
Keywords : Stefan problem, glaciation, gas pipeline, versions of the coefficient smoothing method, numerical solution, computational examples.
References
1. Vasilev F. P. On finite difference methods for the solution of Stefan's single-phase problem. USSR. Comput. Math. Math. Phys., 1963, vol. 3(5), pp. 1175-1191.
2. Kurbatova G. I., Ermolaeva N. N. Analysis of the cylinder glaciation models in seawater. Applied Mathematics and Information Sciences, 2017, vol. 11 (3), pp. 925-930.
3. Samarskii A. A., Moiseyenko B. D. An economic continuous calculation scheme for the stefan multidimensional problem. USSR. Comput. Math. Math. Phys., 1965, vol. 5 (5), pp. 43-58.
4. Budak B. M., Sobol'eva E. N., Uspenskii A. B. A difference method with coefficient smoothing for the solution of Stefan probkems. USSR Comput. Math. Math. Phys., 1965, vol. 5 (5), pp. 59-76.
5. Eyres N. R., Hartree D. R., Ingham J. et al. The calculation of variable heat flow in solids. Philosophical Transactions of the Royal Society. Series A, 1946, vol. 240, pp. 1-57.
6. Ermolaeva N. N., Kurbatova G. I. Nestacionarnaja model narastania morskogo l'da [Nonstationary
model of the sea ice growing]. Vestnik of Saint Petersburg University of Technology and Design. Series Natural and Technical Sciences, 2017, iss. 1, pp. 3—8. (In Russian)
7. Tichonov A. N., Samarskii A. A. Equations of mathematical physics. New York et all., USA, Pergamon Press Ltd., 1963, 781 p.
8. Oleinik O. A. Ob odnom metode reshenia obshei zadachi Stefana [On one solution method of the general Stefan problem]. Paper of AN USSR, 1960, vol. 135, no. 5, pp. 1054-1057. (In Russian)
9. Budak B. M., Vasilev F. P., Uspenskii A. B. Raznostnie metodi reshenia nekotorih kraevih zadach tipa Stefana [Difference methods of the solution some boundary Stefan problems]. Chislennie metodi v gasovoj dinamiki [Numerical methods in gas-dynamics]. Moscow, Calculating Centre of Moscow State University Publ., 1965, pp. 139-183. (In Russian)
10. Krivovichev G. V., Ermolaeva N. N., Kurbatova G. I., Miheev S. A. Mathematical modelling of the glaciation process. International Journal of Pure and Applied Mathematics, 2017, vol. 113 (5), pp. 609-616.
11. Shamsundar N., Sparrow E. M. Analysis of multidimensional conduction phase change via the Enthalpy Model. J. Heat Transfer, 1975, vol. 97(3), pp. 333-340.
12. Vasilev V. I., Maksimov A. M., Petov E. E., Cipkin G. G. Teplomassoperenos v promerzajushih i protaivajushih gruntah [The heat and mass transfer in the freezing ground]. Moscow, Nauka, Physmatlit Publ., 1996, 224 p. (In Russian)
Author's information:
Galina I. Kurbatova — Dr. Sci. in Physics and Mathematics, Professor; [email protected]
