Сходимость разностной схемы для модели термодинамики морского льда Текст научной статьи по специальности «Математика»

Труды Карельского научного центра РАН № 3. 2010. С. 87-92

УДК517.954, 517.929.7, 519.б33.б

СХОДИМОСТЬ РАЗНОСТНОЙ СХЕМЫ

ДЛЯ МОДЕЛИ ТЕРМОДИНАМИКИ МОРСКОГО ЛЬДА

И. А. Чернов

Институт прикладных математических исследований Карельского научного центра РАН

В статье рассматривается локально-одномерная модель термодинамики морского льда, представляющая собой краевую задачу теплопроводности со свободной границей и нелинейными граничными условиями. Для этой задачи построена неявная разностная схема; доказана равномерная по отношению к шагам сетки ограниченность сеточного решения и его производных и вытекающая из этого сходимость сеточного решения к обобщенному решению задачи. Отдельной трудностью является граничное условие I рода на границе вода-лед: рассматривается последовательность задач с условиями III рода и соответствующий предельный переход.

К л ю ч е в ы е с л о в а : термодинамика морского льда, краевые задачи теплопроводности, задача Стефана, свободная граница, сходимость разностных схем.

I. A. Chernov. CONVERGENCE OF THE DIFFERENCE SCHEME FOR THE MODEL OF SEA ICE THERMAL DYNAMICS

We consider the locally one-dimensional model of thermal dynamics of sea ice. This is a boundary-value problem of heat conduction with a free boundary and nonlinear boundary condition. For this problem we construct the lattice numerical method and prove that the lattice solution and its lattice derivatives are bounded uniformly with respect to the lattice steps. This implies the convergence of the lattice solution to the weak solution to the boundary-value problem. A difficulty is the Dirichlet boundary condition of the water-ice boundary: we consider a sequence of problems with Neumann conditions with consequent passing to limit.

K e y w o r d s : thermal dynamics of sea ice, boundary-value problems of heat conduction, the Stefan problem, free boundary, convergence of lattice methods.

Введение

Термодинамика морского льда (процессы таяния и намерзания льда) обычно описывается локально-одномерными моделями [МаукШ;, Ийегйетег, 1971]. Это связано с крупномас-штабностью процессов в океане (море), в ре-

зультате чего поля температур и потоков теплоты и массы можно с высокой точностью приближать кусочно-постоянными полями. Мы рассматриваем локально-одномерную модель термодинамики льда при ряде упрощений. Она представляет собой краевую задачу для уравне-

ния теплопроводности со свободной границей. Свободной границей является поверхность раздела вода-лед, на которой ставятся граничное условие I рода и условие сохранения количества вещества - условие типа Стефана. Потоки теплоты на границах вода-лед и воздух-лед задаются и являются достаточно произвольными нелинейными монотонными функциями температуры льда на указанных поверхностях. Таким образом, имеем параболическую краевую задачу со свободной границей, нелинейными граничными условиями III рода и Стефана и граничным условием I рода. Свободную границу «выпрямляем» заменой переменных. Условие I рода заменяем последовательностью условий III рода с малым параметром при производной. Строим разностную схему и доказываем ее сходимость к обобщенному решению задачи [Ладыженская, 1973]. Предельный переход по малому параметру дает решение исходной задачи.

Постановка задачи

Пусть ось г направлена вниз, за ноль принимаем верхнюю поверхность льда. Толщина льда -Ь(), температура - Т (,г). Не учитываем солнечную радиацию и осадки, а также считаем температуру воздуха достаточно низкой, так что верхняя поверхность льда не меняется. Тогда распределение температуры льда описывается следующей краевой задачей со свободной границей [МаукШ:, Ийег^етег, 1971]:

дТ = ад]Т , ге (0,Ь()), ге (о,г), (1)

ад гТ (,0) = -Л (,Т(,0)), (2)

Т(,Ь()) = Т , (3)

Ь = ~/ъ (, Т((Ь)) + ад Т(, Ь(), Ь > Ь > 0, (4)

Т (0, г ) = То (г), Ь(0) = Ьо. (5)

Здесь а - коэффициент температуропроводности. Температура льда не может быть выше точки плавления Т , которая зависит от солености. Будем пока считать ее постоянной. На верхней поверхности задан тепловой поток (,Т). Считаем, что (, Т) (непрерывно дифференцируемый с ограниченным градиентом) монотонно убывает по Т, (,Т )< 0, (,-&)> 0. На

границе вода-лед температура равна точке плавления - граничное условие I рода. Поток тепла (положительный или отрицательный) /Ь (,Т) задается, его свойства аналогичны (,Т). Поло-

жительная разница между приходящим теплом и диффузионным оттоком от поверхности

вглубь льда идет на компенсацию скрытой теплоты таяния льда (удельную теплоту таяния единицы объема считаем единичной); аналогично отрицательная разница компенсирует тепловыделение при кристаллизации. Закон сохранения энергии приводит к условию типа Стефана (4). Таким образом, свободная граница может двигаться в обоих направлениях (утолщение либо вымывание льда), однако считаем, что при достижении минимально возможной толщины Ь лед исчезает; во всяком случае, его эволюция описывается уже другими моделями. Также лед не может превышать некоторую предельную толщину Ь . Итак, далее всюду Ь е [Ь,Ь ].

С целью устранить свободные границы сделаем замену переменной по формулам и(,х) = Т(,г), г = Ьх, Ьд2= дх ,

дТ = ди -дхихЬ-1Ь. Получим задачу на [0,1] с дополнительным слагаемым в правой части уравнения теплопроводности. Граничное условие первого рода осложняет исследование задачи, поэтому заменим «жесткие» условия I рода последовательностью условий III рода (8) с параметром аь ; кроме того, допустим зависимость коэффициента а от времени (требуя гладкости, ограниченности вместе с производной и отделенности от нуля):

Ь2ди = аОи + дхихЬЬ, хе (0,1), (6)

а( )д и (,0) = -Ь/3 (,и(,0))) (7)

а()дхи(,1)= -аъЬ((,1)-Т), (8)

Ь = -/„ (,и (,1))-аь (и (,1)-Т), (9)

и(0, х) = Т0 (х), Ь(0) = Ь) . (10)

Разностные аппроксимации

Введем в [0,1]х[0,г] равномерную сетку с шагами к и т, индекс / меняется от 0 до I, индекс п - от 0 до N. Введем обозначения и'п=и (, х{)

и аналогичные. Коэффициенты уравнения (6) и правую часть (9) аппроксимируем на предыдущем слое, остальные величины - на текущем. Исключая граничные узлы, получаем множество узлов. Границу обозначим дОы . Введем обозначения для сеточных производных

дип = (ип - ип-1)/т, дкип = ( - ип)/к,

д2кип= дкдкиг- . Представим выражение

- ((пт, (1)+ аЬ (-1 - Т ) = ) + Вп , где

в+ > о, Вп <0.

Заменяя производные сеточными аналогами, получим неявную схему

Ь2 д игп=ап д кК+Пкд иВ +

+^кд и:1 вп (11)

апди = -ЬП fs(пт,ип ) (12)

апди1-1 = -аьЬп (ип - Т), _ (13)

дтЬп+1 = -А(пт,и[_1)-аЬ (-1 -Т), (14)

и0=Т0 (Ь0/к). (15)

Эта разностная схема имеет первый порядок аппроксимации.

Принцип максимума для сеточной задачи

Докажем аналог принципа максимума для параболических задач и ряд следствий. Экстремумы далее могут быть нестрогими.

УТВЕРЖДЕНИЕ 1. Пусть сеточная функция и'п определена в Бы и удовлетворяет (11)-(15).

Тогда она ограничена в равномерно по к, т и аЬ .

ДОКАЗАТЕЛЬСТВО. Пусть и'п достигает минимального значения в узле (п,П )е . При-

меним (11). Левая часть неположительна, все слагаемые правой части неотрицательны. Противоречие не возникает только если и'п = и+ . Применяя рассуждение к узлу (п, п +1) и далее, получаем, что если достаточно малый минимум достигается в объеме, то он достигается и на границе также. В узлах вида (п,1) применим

граничное условие (13): при игп <Т имеем и[-1 < и1п , что противоречит минимальности. Аналогично рассуждение для узлов (п,0) (с применением (12)). Если же минимальное значение достигается при п = 0, то оно ограничено, причем оценка зависит только от начальных условий. Доказательство для максимума совершенно аналогично.

Пусть и'п удовлетворяет (11)-(15) в Бы .

Обозначим через ' с множество, на котором определена сеточная производная ди'п . Введем также Оы' аналогично Оы .

УТВЕРЖДЕНИЕ 2. Пусть и'п удовлетворяет (11)-(15) в при достаточно малых к и т. Тогда дкигп ограничена в Бы' равномерно по к и т.

ДОКАЗАТЕЛЬСТВО. Сначала убедимся, что значения дhU'n на границе DN' равномерно ограничены. Для узлов вида (n,0), (n,I -1) это сразу следует из граничных условий (12), (13) и утверждения (1). Для начальных узлов вида

(0,i) в силу (15) д hUn=boдT (0 (хг+ в)К где 0 < в < h , и потому ограничена. Произведем сеточное дифференцирование уравнения (11) по

h с учетом дhU'n+l = дhU'n + дhдhU'nh :

К д , д hU'n=an д h д hU'n+ (i + 1)hB+д h д hUn + + ihBn д h д hUn-1 +B+д hU:+B-n д hUn.

Введем следующую сеточную функцию ре-куррентно: Mn =Mn-1 (1 - rb-2B+ )-1 , M0 = 1. Эта функция - сеточный аналог экспоненты:

д т Mn = К bfMn . Пусть Vi = д hU'n / Mn и перепишем уравнение в виде

К Mjn-LдУ1 = anд2V‘ + (i + 1)hBn+ дhV‘ +

n ъ f т n n h n \ ' n h n

Mn

+ ihB- д ,V‘-1 +B-V‘.

nhnnn

Так как B+b-2 ограничено, то при достаточно малом т (обеспечивающем положительность знаменателя) функция Mn > 1. Поэтому минимальное и максимальное значения функции Vi не достигаются в DN': в противном случае левая и правая части уравнения имеют разные знаки. Возможен только нестрогий экстремум, достигаемый и на границе. Однако значения дhU'n,

а следовательно, и Vn , равномерно ограничены на границе, так что Vi равномерно ограничена в Dn . Отсюда следует и равномерная ограниченность функции дhU'n =My'n . В самом деле, при т < т0

Mn <(1 - 1т )-t/т < (1 - ИЧ) -'t/то,

где ц = max B +. Здесь использована монотонность функции (1 - x )-1/x при x е (0,1).

Пусть U'n удовлетворяет (11)-(15) в DN . Обозначим через DN с DN множество, на котором определена сеточная производная дtU'„ . Введем также D°N аналогично DN .

УТВЕРЖДЕНИЕ 3. Пусть U'n удовлетворяет (11)-(15) в Dn при достаточно малых h и т. Тогда дtU'n ограничена в DN равномерно по h и т.

ДОКАЗАТЕЛЬСТВО. Сеточное дифференцирование по т уравнения (11) приводит к соотношению, содержащему следующие слагаемые:

1. Слагаемое с производной д tд tU1n в левой

части, положительное, если дtU„ - максимальное значение;

2. Слагаемые с производными д2дtU'n,

дhдtU'n и дhдU-1 в правой части, отрицательные, если д U' - максимальное значение;

n

3. Слагаемые с производными дт B,+ , дт B'n в правой части и д тЬ„ в левой, которое перенесем в правую часть;

4. Слагаемое дтanд2Un в правой части.

Коэффициенты при указанных величинах

равномерно ограничены, что следует из доказанных выше утверждений. Слагаемое пункта 4 заменим, выразив д2Un из (11); получим слагаемое д^дтUn /an , где производная от a - в средней точке, а также несколько ограниченных слагаемых. Сеточные производные коэффициентов из пункта 3 преобразуем следующим образом:

дт B+ = B+ + B+ дт bn+BU дт Un-1 ,

аналогично дтB'n . Здесь величины B+ ,B+,BU

ограничены независимо от шагов сетки; в случае непрерывности производных коэффициентов это - производные в средних точках. Поэтому первые два слагаемые ограничены независимо от шагов сетки.

Предположим, что дzU„ - максимальное значение, и оценим слагаемые, содержащие сеточные производные функции Un по т. Они в совокупности оцениваются величиной (дta / an +Y)Un , где Y определяется оценками коэффициентов при сеточных производных функции Un по т. Без ограничения общности можно считать, что (д ta / an+Y )< 0 и отделена от нуля. В противном случае сделаем замену переменных в = Y-1exp(- 2f), y = exp(-2f)x/Y . Для простоты пусть в исходных переменных a = const. Оценка Y коэффициентов не изменится в новых переменных, а указанное неравенство будет выполнено. Тогда при достаточно большом положительном значении дzU„ сумма слагаемых, содержащая сеточные производные функции Un по т, будет отрицательной и дос-

таточно большой по абсолютной величине. Таким образом, правая часть отрицательна, если д и' - достаточно большое положительное

максимальное значение. Левая часть - неотрицательна, что приводит к противоречию. Поэтому д иП не может принимать в 0°ы слишком

больших значений, т. е. ограничена сверху. Отметим, что оценка не зависит от шагов сетки, а только от коэффициентов уравнения и интервала времени. Рассуждение для минимума дословно такое же.

Осталось проверить граничные узлы. Для п = 1 рассмотрим (11) с учетом

и 1 = и 0 +тдти 1:

Ъ\ д и 1 = ад и0 +1кд и0 В+ + ікд и0-1 В- +

+ т(д к д ( 1 +ікд „ д и 1В1+ +1кд „ д и г1 В[) Если д и 1 - максимальное (минимальное)

значение, то выражение в скобках справа неположительно (неотрицательно). Применяя теорему о среднем к сеточным производным от начального распределения, получаем с учетом ограниченности начального распределения в С2([0,Ъ0]), что ди 1 ограничено независимо от шагов.

Для і = 0 продифференцируем на сетке условие (12):

(ъ Л , ч

д к ди 0 =~дт — ^, и1 )-

V an У

b

^ (д ,f +д uf, ди)

С учетом ограниченности частных производных (они в средних точках) и коэффициентов, а также монотонности АЬ^,и) по второму аргументу, получаем, что д к ду1 > 0 при достаточно большом значении д и1 , что противоречит максимальности. Для минимума доказательство дословно повторяется. Случай п = I рассматривается аналогично (с помощью условия (13)).

В итоге получаем, что функция дти'п ограничена в Б0 равномерно относительно шагов

сетки к и т .

Для решения системы (11)-(15) применим стандартный метод прогонки. На каждом слое времени вычисляется Ьп из (14), затем прого-

ночные коэффициенты для линейных уравнений (11) и (13), и в (12) выражаем и\ через и° с помощью найденных коэффициентов. Получа-

a

n -1

ется одно нелинейное уравнение для и . Находим и0 и далее обратной прогонкой определяем последовательно и'п .

Сходимость аппроксимаций и обобщенное решение

Итак, сеточные функции и'п равномерно ограничены независимо от шагов сетки вместе со своими сеточными производными д ип и

дти'п . Построим в П = [О,Г]X[0,1] функцию и (?, х) сплайновой интерполяцией первого порядка сеточной функции и'п .

По построению и^, х)е Н1 (П) (соболевское пространство функций, обладающих первыми обобщенными производными из £2(П) [Ладыженская, 1973; Михайлов, 1976] при любых (достаточно малых) шагах к и т, причем семейство и ^, х) ограничено по норме Н1 (П). Выберем последовательность шагов (кк,тк), сходящуюся к нулю, и соответствующую последовательность ик (?, х). Она ограничена в гильбертовом пространстве Н1 (П), поэтому можно выделить слабо сходящуюся подпоследовательность, за которой сохраним прежнее обозначение ик (?, х); ее предел обозначим и(?; х)е Н1(П). При этом обобщенные производные от ик (^, х) сходятся к таковым от предела. Эта функция непрерывна, так как равномерно ограниченные непрерывные члены последовательности ик ^, х), обладая равномерно ограниченными кусочно-постоянными производными, равностепенно непрерывны, что влечет равномерную сходимость (возможно, после перехода к подпоследовательности - теорема Арцела).

Функции Ьк (^), полученные сплайновой интерполяцией первого порядка сеточных функций Ьп , равномерно ограничены и равностепенно непрерывны, поскольку обладают равномерно ограниченными (в силу ограниченности правой части (14)) кусочно-постоянными производными. Переходя к подпоследовательности (теорема Арцела), получаем непрерывную функцию Ь(0 .

Назовем пару функций и(^, х)е Н1(П),

Ь^)е С([0,Г]) обобщенным решением [Лады-

женская, 1973; Михайлов, 1976] краевой задачи (6)-(10), если и^, х) и Ь^) удовлетворяют интегральному тождеству

| V (0, х)Т0(Ь0 х)ёх +1 идуйхЛ

П

/ ид, (хV /ь (и('■ О)-ад У (1) лл = 0

ПЬ при любой непрерывной V(^, х)е Н1 (П), такой, что V (Г, х)= 0 и V (^1)= 0, и тождествам

Ь(0 = Ь0 -1(/ь(и(,1))-Ьди(^,1)^, (17)

I

| )д хи (,1)+аьь(и (,1)- Т )^ )Л = 0 (18)

0

для любой непрерывно дифференцируемой функции н( 0 .

Построенные функции и ( I, х) и Ь( 0 образуют обобщенное решение в смысле данного определения: на аппроксимациях тождества выполняются с погрешностью 0(к, т), равномерная сходимость аппроксимаций позволяет перейти к пределам Ьк ( 0^ Ь(0 , ик ( Г, х)^ и( г, х), а слабая сходимость - к пределам

д х и к ( *,0Н д и ( ф, д хик ( М)^ д и (М).

Поскольку в исходной задаче имелось краевое условие I рода (3), то следует перейти к пределу при аь ^ +^ . Интуитивно ясно, что это

повлечет и(г,0)= Т , покажем это.

Достаточно показать, что

I

| ад хи и,1^(г)& =

0

£ Г 1 Л

= Ц аЭхи(,0) +1ад2хи(,х)х w(t)dt =

0 V 0 )

£ I 1

= | ад хи (,0)>^(г )dt + Ц )д и (, х )dxdt =

0 0 0

г 1 Г

= |адхи(t,0)w(t)dt -||Эн(г)и(,х)dtdx 0 0 0

1

| (н(Г)и и?', х)- н(0)и (0, х ))х

+

+

ограничено; выкладку следует понимать на аппроксимациях с последующим переходом к пределу. Ограниченность левой части следует из

ограниченности всех слагаемых правой части: и (^ х) ограничено, а д хи (t,0) ограничено в силу ограниченности и (^ х) и (7). Тождество (18), таким образом, при аЬ ^+^ дает и(,1) = Т в

силу основной леммы вариационного исчисления.

Отметим, что построенная свободная граница Ь(0 является абсолютно непрерывной, поэтому обратная замена переменных является гладкой. Пару непрерывных функций Т), Ь(0 назовем обобщенным решением исходной задачи, если Т (,Ь(0) = Т и соответствующая пара и(^ х) и Ь(0 удовлетворяет тождествам (16) и (17) при любой непрерывной V(^ х)е Н, (П), такой, что V(Г, х)= 0 и V(^1)= 0. Полученная обратной заменой из построенной функции и (^ х) функция Т ) в паре с построенной функцией Ь(0 является обобщенным решением в смысле данного определения.

Выводы

Таким образом, мы построили неявную разностную схему для краевой задачи III рода термодинамики морского льда и доказали ее сходи-

СВЕДЕНИЯ ОБ АВТОРЕ:

Чернов Илья Александрович

старший научный сотрудник, к. ф.-м. н.

Институт прикладных математических исследований КарНЦ РАН

ул. Пушкинская, 11, Петрозаводск, Республика Карелия, Россия, 185910

эл. почта: [email protected] тел.: (8142) 766312

мость к обобщенному решению задачи. Тем самым доказано существование обобщенного решения этой краевой задачи со свободной границей, нелинейными граничными условиями III рода. Доказано также существование решения задачи с граничным условием I рода, представляющим особую сложность в силу своей «жесткости». Поскольку порядки величин тепловых потоков известны, предложенная схема может применяться для решения задачи при достаточно большом параметре аь (таком, чтобы отличия температуры на границе вода-лед от точки замерзания были малы сравнительно с точностью модели).

Литература

Maykut G. A., Untersteiner N. Some Results from a Time-Dependent Thermodynamic Model of Sea Ice // Journal of Geophysical Research. 1971. Vol. 76, N 6. P. 1550-1575.

Ладыженская О. А. Краевые задачи математической физики. М.: Наука, 1973. 408 с.

Михайлов В. П. Дифференциальные уравнения в частных производных. М.: Наука, 1976. 391 с.

Chernov, Ilya

Institute of Applied Mathematical Research, Karelian Research

Centre, Russian Academy of Science

11 Pushkinskaya St., 185910 Petrozavodsk, Karelia, Russia

e-mail: [email protected]

tel.: (8142) 766312