Алгоритм решения системы уравнений течения льда в трех-мерной математической модели Текст научной статьи по специальности «Физика»

УДК 910.1 551.8

алгоритм решения системы уравнении течения льда в трехмерной математической модели

© 2010 г. О.О. Рыбак,1 Е.А. Рыбак2

г

1Сочинский научно-исследовательский центр РАН, ул. Театральная, 8а, г. Сочи, Краснодарский край, 354000, snic@soch.ru

2Государственный южный научно-исследовательский полигон РАН, ул. Театральная, 8а, г. Сочи, Краснодарский край, 354000

1Scientific Research Centre RAS, Teatralnaya St., 8a, Sochi, Krasnodar Region, 354000, snic@soch.ru

2State Southern Scientific Research Polygon RAS, Teatralnaya St., 8a, Sochi, Krasnodar Region, 354000

Описание течения льда в горных ледниках и ледниковых щитах с большим пространственным разрешением требует максимального учета действующих напряжений. Это приводит к чрезмерному усложнению вычислительной задачи. Однако для широкого класса процессов задача может быть сведена к так называемой аппроксимации неполного второго порядка уравнений динамики и к решению системы двух нелинейных эллиптических уравнений. Предлагается алгоритм численного решения такой системы уравнений. В его основе лежит представление уравнений в виде суммы специальных операторов и расчет эффективной вязкости на смещенной сетке. Алгоритм позволяет построить эффективный программный код.

Ключевые слова: лед, ледник, ледниковый щит, течение льда, математическая модель.

Description of an ice flow in the glaciers and in the ice sheets with the high spatial resolution requires making use of as many stresses as possible. Because of this the numerical problem becomes overcomplicated. Nevertheless, to a wide class ofproblems one can apply so-called incomplete second order approximation of the dynamic equations, which leads to a system of two non-linear elliptical equations. Suggested in the paper, is the algorithm of a numerical solution of such a system. It is based on a conversion of the initial equations into a set of special operators and further calculation of the effective viscosity on a staggered grid. This algorithm enables to build an effective computer code.

Keywords: ice, glacier, ice sheet, ice flow, mathematical model.

Лед отличается от других природных твердых тел относительно высокими скоростями течения. Это обстоятельство было замечено очень давно, преимущественно на примере горных ледников. К настоящему времени сформировался общепринятый взгляд на процесс деформации природного льда. Под действием нагрузки (и если не превышен предел прочности) лед деформируется без разрушения и изменения объема. Деформация ледникового объекта может быть представлена в виде суммы упругой обратимой деформации и пластической остаточной деформации. Специфическое свойство льда (ползучесть) можно охарактеризовать как изменение пластической деформации во времени под действием напряжения (приложенной силы) или как процесс релаксации при постоянном напряжении [1]. Гипотеза, связывающая напряжения и деформации льда, в окончательном виде была сформулирована Дж. Гленом [2]. Она лежит в основе современного математического аппарата описания динамики природных ледниковых объектов. В общем случае математическая модель ледника или ледникового щита строится в предположении о том, что лед представляет собой вязкую, неньютоновскую, тепло-проводящую, несжимаемую жидкость [1, 3].

С.С. Григоряном и П.А. Шумским [4] для описания динамики трехмерного нестационарного ледника был предложен метод тонкого пограничного слоя, позволявший упростить описание течение льда в тех случаях, когда кривизной поверхности и ложа ледника можно было пренебречь. Позднее К. Хуттер [3] сформулировал аналогичную концепцию «аппроксимации мелкого льда» (Shallow Ice Approximation (SIA)). В рамках SIA считается, что деформации льда

обусловливаются лишь тангенциальными напряжениями в вертикальной плоскости, а остальные напряжения в балансе сил несущественны. Модели, в основе которых лежит SIA, вполне адекватно воспроизводят динамику ледниковых щитов в целом. Вместе с тем за рамки возможностей SIA выходит описание динамики отдельных областей щитов, таких как ледовые купола и ледоразделы, области вблизи линии налегания [5]. Ограничения в пространственном разрешении не позволяют учесть влияние сложной топографии подстилающей поверхности на поле скоростей течения льда. За рамками SIA остается динамика шельфовых ледников и переходной зоны между шельфовыми ледниками и ледниковым щитом, а также, разумеется, динамика горных ледников.

Учет всех напряжений в уравнениях динамики льда приводит к чрезмерному усложнению вычислительной задачи. В качестве компромисса в [6, 7] была предложена так называемая аппроксимация неполного второго порядка уравнений динамики льда. На ее основе была разработана трехмерная модель течения льда [5], которая была использована в различных прикладных исследованиях [8-10]. В настоящей работе рассматривается способ модификации модели [5] и предлагается новый алгоритм ее реализации.

1. Напряжения, деформации и закон течения льда. Деформации льда практически не зависят от гидростатического давления, поэтому его можно не рассматривать при описании течения льда под действием гравитации [11]. Пусть с - компоненты полного

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

Р = + ауу + а . Разложив тензор напряжений на

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

1 я т.. = а.. —а,,д..,

ij 3 J

(i)

где 8у - символ Кронекера. При соблюдении условия

изотропности и несжимаемости льда главные оси тензора скоростей деформации параллельны осям тензора напряжений, а компоненты обоих тензоров пропорциональны друг другу [1, 11]. Вторые инварианты тензоров напряжений и скоростей деформации, или эффективные напряжения и эффективные скорости деформации, определяются как

2е2 = е2 + е2 + е2 + 2(е2 + е2 + е2 ), (2)

е хх уу гг \ ху XI уг / т ^ '

2т2 = т2 +т2 +т2 + 2(г2 +т2 +т2 ). (3)

е хх уу гг \ ху хг уг] х '

Закон Дж. Глена [2] связывает тензор скоростей деформации с тензором напряжений

¿е = Ат^, (4)

где А зависит от температуры льда [11]. Значение п=3 соблюдается для широкой полосы напряжений в природных процессах [11].

2. Исходные уравнения модели. Сохранение импульса единичного объема льда можно записать как

с „

= + ' (5)

са

где V - трехмерный вектор скорости; t - время; а - симметричный тензор напряжений Коши с компонентами ст.; р - плотность льда [5]. Локальными и конвективными производными в уравнении (5) пренебрегают. Считая далее, что V- V = 0, даж!сх = 0 и = 0, получим да^/с! = р ,

можно переписать (5) в скалярной форме:

да„

— ду да,.

дт

dz = P g •

■ = 0,

да да,,

дс

ду

а

■ = 0,

(6)

— (2т

дс ■

XX + Т уу ) +

дтсу дт„

—у

—

d

= P g— : dX

д

(2т +т )+-

уу XX

дт дт

= P g-

d

(7)

смена режима течения, например, в непосредственной близости от ледоразделов [5].

Скорости деформации льда связаны с компонентами скорости течения. Считая, что по абсолютной величине градиенты вертикальной скорости много меньше градиентов горизонтальной скорости, т.е. д/ сдх «¿V/ д и д/ ¿у «¿V/ Сд , можно выразить

компоненты тензора скоростей деформации через градиенты скорости течения:

• (8)

Комбинация (2), (4), (7) и (8) дает систему нелинейных эллиптических уравнений для горизонтальных компонент скорости течения льда в приближении неполного второго порядка:

сх ^ сх) — ^ ¿у) ¿у ^ ¿у)

(■ ■ \

S е е

XX хУ xz

е е е

ух уу у2

е е е

V 21 2у zz

— 1 ( — —V ) 2 V— дс ) 1 —

дс 2 —

1 (— —V ) 2 V— дс ) — 1 д

ду 2 д—Г

1 — 1 —V д

2 — 2 — д

д ( — + — I Л—

—у V дс

д ( — + — I Л—

— V —

= P g

äs_ дс

(9)

ду V ду) ду ^ дсJ дс [ — J

—с I —С ) — I д) i —С

где

Л = 1 A(r ■ )-1

l/n

2 + (— + U

^ — — 1 ( — — —X —у 41 —у —X

4

4

(10)

(1-й)/ 2n

Исключим гидростатическое давление из тензора напряжений а и подставим выражение для девиатора (1) в первые два уравнения системы (6). Пренебрежем атмосферным давлением в левой части третьего уравнении из системы (6) и проинтегрируем его от подстилающей поверхности до поверхности щита, считая, что на поверхности напряжение а^ 1х,у) = 0. Подставим аг = -р^в(х, у) в первые два уравнения (6) и получим уравнения для девиаторов напряжений в приближении неполного второго порядка [6, 7]

¿у4 уу хх/ ¿у — ' ¿у

Физический смысл приближения (7) заключается в пренебрежении компонентами так называемого вариационного напряжения. Эти члены на два порядка меньше по величине, чем другие трансверсальные и нормальные члены даже в районах, где происходит

эффективная вязкость.

Более подробные выкладки приведены в [5]. В качестве краевого условия на нижней границе берется скорость базального скольжения, определяемая в соответствии с [12]. На верхней границе вводится условие равенства нулю суммы всех напряжений [5].

3. Введение безразмерной вертикальной координаты. Введем систему координат (х, у,С), где вертикальная координата нормирована на толщину льда, С = (в - г )/ Н .

Формулы, согласно которым преобразуются горизонтальные координаты, приведены ниже.

Преобразование вертикальной координаты в уравнениях течения льда к безразмерному виду приводится по [13]. Определим безразмерную вертикальную координату как С = (в -2)/Н, где 5 - абсолютная высота поверхности; Н - толщина льда. На поверхности щита С = 0, на подстилающей поверхности -С = 1. Операция преобразования координат переводит (х, у, г) ^ (х', у', С). Первая производная может быть записана как

2

+

X

2

2

2

1

1

+

+

+

+

+

С = + + (П)

ск ск' (Эх су' (Эх сС ск Аналогичные выражения получаются и для (/ ду. Допустим, что ск'/ск = 1 и д'( = 1, а ск'(, ду'/¿к, ск'/ск , су'/ск = 0. Эти соотношения верны в условиях, когда градиенты поверхности ледника и подстилающей поверхности малы. В этом случае (11) может быть записано как

¿f ¿f ¿f — = + ах

ск ск ' ¿Ç

(12) (13)

1 с

ск Н С ' Из (12) и (13) следует, что вторые производные можно записать как

с2С с2С , с 2 С „ с2С

7 - ■ + Ъ — + аI — + 2« 7

ск2

ск '

¿Ç к ¿Ç

ск ' dÇ

1 С2 f

cz2 H2 ¿Ç2

c2f

1

¿хс h

s

1 ¿Î с

H ск ' ¿Ç ск ' ¿Ç

С2 f С2 f> - a

¿Ç2

dLf дЧ cf дЧ дЧ дЧ

- + с„,— + а,,-— + а„-— + а„а„

дкду С ' f 1

fÇ у С ' fÇ к С ' fÇ

¿Ç2

¿S „ сН где a = — I--Ç-

âir ¿0x Ьк =—к + ак—-

к ¿к' к ¿Ç

J_

H ,

ск с

i ^2

¿0^ ск'

с = —— + а

ку^г к

¿ау

¿Ç

для с , ау, и Ъу.

Применим преобразования к первому из уравнений системы (9):

4

"с ' с Л + 4ах "с с | + 4ах " с ' с У

V ¿f )_ пс У V ¿f )_

¿с с сÇ

+ 4а2 —\п—I + —\п—I + а.

H2

с ( с с I с с I с

с ( ¿и —\ п—

с I

с I с

су [пс

+ ау —\ п— I + а2 —\ п — I + (14)

с I с

¿%\пс%

+ с п +пьу )=

¿s ~ ¿V

= Pg--3пску--2

¿к ху ¿Ç

с( ¿V У — 2а

_с 1 ¿У )_

у

— 2а,

_с_( с ÇV с

су V ¿Ç

—3акау

¿1 с

сÇ \п сÇ

-Цп*

¿ÇV с

ск [псÇ

с ( ¿V —\ п —

су ^ с

Аналогичным образом преобразуем выражение для эффективной вязкости (10):

п = 1 a{t * )-

(15)

¿1 ¿1 2

+ 4ак--+ \ 4ак +—^ +1

к ¿Ç ¿к \ к H2

с

+ 5

+ 4а,

¿V ¿V

¿с,'

+ \ 4а 2 +-

H2

■ +1

, âi dv âi dv âv âi âi âv + 4--+ 4a--+ 4a--+ 4a a--+

âx ây x âÇ ây y dÇ âx x y dÇ dÇ

âi âi âv âv ncU âv „ âi âv

+ 2a--+ 2a--+ 2--+ 2a--+

y âÇ ây x âÇ âx ây âx y âÇ âx

âv âi âi âv

+ 2ax---+ 2ayax--

x âÇ ây y x âÇ âÇ

4. Численное решение. Численное решение, предложенное в [5], предполагает раздельную пространственную дискретизацию компонент скорости и эффективной вязкости. Мы предлагаем воспользоваться тем, что (15) состоит в основном из операторов вида

п(М) МпЩ, yJ-L*-

¿kV ¿к) сх I сх

(16)

с2 s гс2Н „ сНЛ

--Ç--2а -

¿с'2 Ç ¿x'2 х ск

¿а

+ ау—^ . Аналогичные в^1ражения получаются

где (р = и, V, Л = х, у,^, / ф ] , и непосредственно перейти к дискретизации операторов (16). Произведя соответствующую замену, получим уравнение для u <Ю(и, х) + 4а х Т(и, х, С) + 4ах т(и, С, х) + (17)

+ ^ 4а2 + а2у + С) + П(и, у) + а у Т(и, у, С) +

+ а у Т(и, С, у )+(Ъу + 4ЪХ = ж с - х, у) -

- 2ауТ^, х, С) - 2ахТ^, С, у) - ЗахауQ(v, С) -

- Т^,у,х)- ахТ^,у,С)- ауТ^,С, х)- Зс^Л^

и аналогичное выражение для v

4^, у) + 4ау Т^, у, С) + 4ау Т(ч С, у) + (18)

+ 14а2 + ах2 ^ 2

+ |Q(v, С) + х) + ауТ(v, х, С) +

+ а у Т^, С, х) + (Ъх + 4Ъу Л с = РЕ с - 2Т(и, у, х) -

- 2ахТ(и, у, С) - 2ауТ(и, С, х) - Захау&-(и, С) -

- k, у) - ауТ(u, к, с) - ахТ(u, С, к) - '3скуЛ'~~ .

Система нелинейных уравнений (9), (10), приведенная к безразмерному виду и записанная в форме (17) и (18), решается методом последовательных подстановок аналогично описанному в [14]. Обозначим номера итераций через ^ ^+2, ... Каждая из итераций распадается на три шага:

- расчет эффективной вязкости л (15) с использованием u и v, найденных на предыдущей итерации. На первой итерации используются локальные ^1А) решения для u и v. Если в отличие от [5] рассчитывать эффективную вязкости на смещенной сетке (рис. 1), можно значительно сократить ошибку линейного приближения производных за счет сокращения пространственного шага вдвое;

X

2

2

5

+

X

2

2

1

+

+

+

1

+

— а

— а

х

у

- расчет и. При этом V и п в первом из уравнений (9) считаются независимыми от и. Первые и вторые производные рассчитываются как центральные разности внутри области и как односторонние разности на границе. Уравнения и граничные условия для и и аналогичное для V в системе координат (х, у, С), будучи переписанными в конечно-разностной форме, преобразуются в систему линейных алгебраических уравнений. В матричной форме их можно записать как Ах = Ь , где х - неизвестная переменная; А - квадратная матрица из элементов ау. В случае, если А Ф 0, единственным решением будет х = А-1Ь. Поскольку лишь незначительное элементов ау отлично от нуля, то матрица А является разреженной. Для нахождения х используется метод сопряженных градиентов, алгоритм которого подробно изложен в [15];

- расчет V проводится аналогично расчету и.

Q(m, х) = —-\Л—-дс V дс

1

2 Ас2

"i+U

Л 11 + Л

1 1 1 ч 1 1

i+—, J+-- i+--, J--

V 2 2 2 2

(

+ U

i-1,J

(19)

Л

Л 1.1+Л

\

1 1 1 1

i--J+— i--,J--

V 2 2 2 2

-Ui,j Л 1 1 +Л 1 1 +Л 1 1 +Л 1 1

—, J+— i+—, j— i—, J-+— i—, j—

V 2 2 2 2 2 2 2 2 )

t(u x, у)=! [Лд Ь

A^Ax (4-1/2,J+1/2 - Л+1/2, J+1/2 )+

+ Ui,J (^-1/2,J+1/2 + Л i+1[2,J-1j2 -7;-1/2,J-1/2 -7i+1/2,J+1/2 ) + + Ui-1, J (Лг-1/2, J-1/2 - Л>-1/2, J+1/2 )+ +Ui-1,J-tfi-1/2, J-1/2 + Ui+1,J+tfi+1/2,J+1/2 + + Ui+1,J Л+1/2,J+1/2 - r?i+1/2,J-1/2 )-

(20)

Ы1+1,7-1^+1/2,7-1/2 иг-1,7+1^-1/2,7+1/2 \

Схематически расчет (19) и (20) показан на рис. 2, 3. Заметим, что индексация по координате С опущена для простоты записи. Множители при им ., и; . и

и являются составляющими частями элементов матрицы А. Присвоим им номера, соответствующие индексам при и (таблица), тогда О.(и, х) = а(15)и 7+1 + + а(7)и

где a(15) =

(

',J+1 1 r-i,J \ / i,J-1

( \

(21)

Л 1

--, J +

V 2,J 2

1 +Л 1 1

— —, J--

2,J 2 )

a

1) = (7) =

Л 11 +Л

V i+1J+-2

11

i+-, J— 2 2

/

+ Л

\

11

i—, J+— 2 2

+ Л

11

2J--2 )

Л 1 1

i- 2,J+2

+ Л

11

i--, J--

2 2 )

индексу /+1, j соответст-

вует номер т = 5; индексу у - номер т = 11; индексу /-1, ] - номер т = 7. Аналогично можем записать т(и, х, у) = а(11)и. ^ + а(15)и. ^ + а(7)и, ^ + (22)

+ а(\9)и1+1,7 + а(3)и;-1,7 + а(21)и,.+1,7+1 + а(1)и,-1,7-1 +

+ a

(17)u

i+1,J-1

+ a

(5)ui-1!

j+1 •

Рис. 1. Смещенная сетка, примененная для решения системы уравнений (9), (10)

Операторы О и Т (16) определяются пространственными производными компонент скорости течения льда. Например, операторы О.(и, х) и Т(и, х, у) можно записать как

Рис. 2. Пространственная сетка для численной аппроксимации О.(и, х). • - узлы основной сетки, X - смещенной;

¿и , ч

] — рассчитывается в узлах смещенной сетки (слева),

¿х

¿и\ „

— I ] — I - в узлах основной сетки (справа)

ах I ¿х )

Рис. 3. Пространственная сетка для численной аппроксимации ^(u, x, у)

>

Порядковая нумерация в матрице коэффициентов

Номер Индекс Номер Индекс

1 i - 1, j - 1, k 12 i, j, k + 1

2 i - 1, j, k - 1 13 i, j, k + 2

3 i - 1, j, k 14 i, j + 1, k - 1

4 i - 1, j, k + 1 15 i, j + 1, k

5 i - 1, j + 1, k 16 i, j + 1, k + 1

6 i, j - 1, k - 1 17 i + 1, j - 1, k

7 i, j - 1, k 18 i + 1, j, k - 1

8 i, j - 1, k + 1 19 i + 1, j, k

9 i, j, k - 2 20 i + 1, j, k + 1

10 i, j, k - 1 21 i + 1, j + 1, k

11 i, j, k

Остальные операторы находятся аналогично. Уравнения (17) и (18) содержат по 14 операторов Q и Т каждое. Суммируя a(m) при каждом из um в левой части (17), получим элементы матрицы A. Суммируя a(m) при каждом из vm, получим элементы матрицы b. Последовательность операций легко реализуется в виде программного кода.

Литература

1. Шумский П.А. Динамическая гляциология. М., 1969. 171 с.

2. Glen J. W. Experiments on the deformation of ice // J. of Glaciology. 1952. Vol. 2. P. 111-114.

3. HutterK. Theoretical Glaciology: material science of ice and the mechanics of glaciers and ice sheets. Dordrecht, 1983. 510 p.

4. Григорян С.С., Шумский П.А. Простейшая математическая модель трехмерного нестационарного ледника // Научные труды / Ин-т механики МГУ. 1975. № 42. C. 43-53.

Поступила в редакцию

5. Pattyn F. A new three-dimensional higher-order thermo-mechanical ice sheet model: Basic sensitivity, ice stream development, and ice flow across subglacial lakes // J. of Geophysical Research. 2003. Vol. 108(B8). Р. 2382.

6. Blatter H. Velocity and stress fields in grounded glaciers: a simple algorithm for including deviatoric stress gradients // J. of Glaciology. 1995. Vol. 41(138). P. 333-344.

7. Baral D.R., Hutter K., Greve R Asymptotic theories of large-scale motion, temperature and moisture distribution in land-based polythermal ice sheets: A critical review and new developments // Applied Mechanics Review. 2001. Vol. 54(3). P. 215-256.

8. Региональная модель динамики льда. Часть 1: Описание модели, постановка численного эксперимента и современная динамика потока в окрестностях станции Конен / О. Рыбак [и др.] // Материалы гляциологических исследований. М., 2007. Вып. 102. С. 3-11.

9. Ice thinning, upstream advection, and non-climatic biases of the upper 89 % of the EDML ice core from a nested model of the Antarctic ice sheet / P. Huybrechts [et al.] // Climate of the Past. 2007. Vol. 3. P. 577-589.

10. Past and present accumulation rate reconstruction in the Eastern Dronning Maud Land, Antarctica / P. Huybrechts [et al.] // Annals of Glaciology. 2009. Vol. 51. P. 112-120.

11. Paterson W.S.B. The physics of glaciers, 3rd edition. Oxford; N.Y., 1994. 480 p.

12. Huybrechts P. The Antarctic ice sheet and environmental change // Berichte zur Polarforschung. 1992. Vol. 99. 241 p.

13. Lliboutry L. Very Slow Flow of Solids. Basics of Modeling in Geodynamics and Glaciology. Dordrecht, 1987. 510 p.

14. Hindmarsh R., Payne A. Time-step limits for stable solutions of the ice sheet equation // Annals of Glaciology. 1996. Vol. 23. P. 74-85.

15. Numerical Recipes / W.H. Press [et al.]. Cambridge, 1992. 963 p.

29 марта 2010 г