Численное моделирование поверхностных волн над подвижным дном Текст научной статьи по специальности «Физика»

Вычислительные технологии Том 12, № 4, 2007

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ПОВЕРХНОСТНЫХ ВОЛН НАД ПОДВИЖНЫМ ДНОМ*

М. Г. Хлжоян Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: [email protected]

Results of numerical simulations of the surface waves generation by an underwater landslide moving alonga plane sloping bottom are presented. The numerical algorithm is based on a finite-difference scheme for the equations in a movingframe of reference, which describes potential ideal flows with a free boundary. Dependence of general behavior of the generated wave on the bottom inclination, width of a landslide and its acceleration is investigated.

Введение

Задача генерации волн движением подводного оползня на сегодняшний день весьма актуальна. В естественных условиях подводный оползень представляет собой движение некоторой массы грунта вдоль склона дна. Большие объемы движущейся массы вызывают на поверхности воды волны, близкие по своим характеристикам цунами. Специфика моделирования таких волн определяется большой длительностью перемещения оползня от берега в глубоководную часть и тем, что они зарождаются в прибрежной зоне с малой глубиной, сравнимой с вертикальным размером оползня.

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 06-05-64869) и Программы интеграционных фундаментальных исследований СО РАН (проект 2006.113).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2007.

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

1. Постановка задачи и вычислительный алгоритм

Рассмотрим плоскопараллельное течение жидкости с поверхностными волнами. Введем декартову систему координат Oxy так, чтобы ось Ox лежала на невозмущенной свободной поверхности, а ось Oy была направлена вертикально вверх. Снизу слой жидкости ограничен подвижным дном, заданным функцией y = -h(x,t), сверху жидкость ограничена свободной границей, описываемой однозначной функцией y = n(x,t), где t — время. Будем рассматривать слой жидкости, ограниченный слева и справа отрезками вертикальных прямых. Область, занятую жидкостью, обозначим через Q(t).

Математическая постановка задачи для потенциальных течений жидкости со свободной границей заключается в определении потенциала скорости p(x,y,t), удовлетворяющего в области Q(t) уравнению Лапласа, и функции n(x,t), описывающей свободную границу, на которой должны выполняться кинематическое и динамическое условия. При t = 0 задаются начальные условия для потенциала p0(x,y), профиля свободной границы no(x) и формы дна h0(x). На левой границе ставится условие непротекания

£ = 0 (1)

на правой x = L — неотражающее краевое условие, которое обеспечивает свободный проход волн через эту границу [3]. Дно подвижно, но непроницаемо для жидкости, поэтому нормальные составляющие векторов скорости жидкости и точек дна совпадают и граничное условие на дне принимает вид

U ■ П = Vbot ■ n = v (2)

или

dp dn

h , = v, (3)

где и = Ур — вектор скорости жидкости с компонентами и, V; Уьоь(х,1) — заданный вектор скорости точек дна с компонентами и^, др/ди = Ур • п, п — единичная внешняя нормаль к границе области П(£).

Поскольку скорость дна ^ь^(х,£) задана, то форму дна в любой момент времени £ > 0 можно найти, решив следующее уравнение относительно Л,(х,£):

Ы + иьо^х + ^ = 0. (4)

В частности, если дно движется как твердое тело, не деформируясь, так что вектор скорости Vbot(t) не зависит от переменной x, то решением уравнения (4) является функция

h(x,t) = h0 í x — J Ubot(r)dr j — J Vbot(r)dr. (5)

V о /0

На рис. 1 изображена область течения для рассматриваемой в настоящей работе задачи о генерации поверхностных волн движущимся оползнем. В прибрежной зоне дно является плоским откосом, образующим с горизонтом угол в и описываемым уравнением y = — x tan в. Слева область течения ограничена вертикальной непроницаемой стенкой, расположенной на расстоянии xa от начала координат. Глубина дна около этой стенки зависит от угла наклона откоса и равна xa tan в. В некоторой точке B плоский откос сопрягается с горизонтальным дном, глубина которого также зависит от в и равна H = xb tan в. Длина области AC и длина наклонного участка AB были во всех расчетах фиксированными и такими, что в процессе движения оползень всегда оставался на плоском откосе и останавливался, не доходя до точки пересечения откоса с горизонтальным дном.

В начальный момент времени вода покоилась, поэтому полагалось, что ^0(x,y) = 0, n0(x) = 0. Форма дна на участке AB при t = 0 задавалась формулой

h0(x) = x tan в — T-

1 + tan h ( ——+ b 1 0.5 cosв

1 — tan h ( ——— b 0.5 cosв

[1 +tan h(b)]2

(6)

Здесь Т — максимальная толщина оползня, соответствующая его вершине; х0 — абсцисса этой вершины; Ь — эффективная длина оползня (вдоль откоса).

При £ > 0 начинается движение оползня. При этом предполагается, что он движется по плоскому откосу вниз как твердое тело, все точки которого имеют один и тот же вектор скорости ?7ъ^(£) = (иъсл(¿)), параллельный поверхности откоса. Поэтому ^Ьсл (¿) = — Иъ<^(£^ап в и закон движения оползня однозначно определяется заданием первой компоненты скорости Иъ<^(£).

В настоящей работе рассматривается следующий закон движения оползня: равноускоренное движение до момента времени ^ с заданным в горизонтальном направлении ускорением а, далее равномерное движение до мгновенной его остановки ¿^р. Расчет

Рис. 1. Схема области течения и расчетная сетка

тем не менее продолжается и после остановки оползня, до момента времени ¿епа. Таким образом,

(а*, 0 < г < ¿1,

«¿1, ¿1 < t < г^, (7)

0 tstop < г < ¿епё-

Подставляя указанные выражения для Мьсл(г) и ■иьсл(г) в формулу (5) и учитывая (6), получаем следующее выражение для функции к(х,г) на участке АВ:

к(х,г) = х 1ап 9 — Т-

1 + 1ап МХ-^ + 6

0.5 сое 9

1 — ^ап М Х к ХС(л — 6

0.5 сое 9

[1 +1ап к(6)]2

где хС(г) — абсцисса вершины оползня:

Хс(г)

0 < г < г1,

аГ

Х0 + "2",

«г 1 / \

хо + ^ + «¿1(г — ¿1), ¿1 < г < ^

«г 1 / \

х0 +--2--+ аt1(tstop — tstop < t < ¿епа-

На левой границе рассматриваемой области ставится условие (1):

= 0,

дх

Х=ХА

на горизонтальном участке дна — условие непротекания

0,

у=-н

на подвижном участке дна — условие (3), при этом

(8)

(9)

кх

п

ubot

угТкХ' лД+кХ

лД+кХ

(—кх + 1ап 9).

Из последней формулы с учетом (8) следует, что функция V является знакопеременной на участке АВ: при х < хс(г) она положительна, а при х > хс (¿) — отрицательна. Это означает, что слева от вершины оползня дно движется в направлении внешней нормали, а справа — в противоположном к нормали направлении.

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

г = ¿, х = х(д ,г), у = ,г).

(10)

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

1

V

условия Дирихле для расчета потенциала внутри области, удовлетворяющего конечно-разностному аналогу уравнения Лапласа. С использованием полученных значений потенциала определялось новое положение свободной границы для данного временного слоя и строилась сетка для следующего временного слоя. Подробное описание алгоритма и вид конечно-разностных уравнений на подвижной криволинейной сетке приведены в работе [3]. Здесь же мы укажем способ вычисления координат расчетной сетки и вид в криволинейной системе координат условия (3) непротекания через подвижное дно:

к + к

дд1 ¿/д2

= -V/gПv, (11)

д2=0

где

, _ ^22 , _ , _ ^12 , _ ^11

к11 - - , к12 - к21 ---, к22 - - ,

11 у 12 21 у 22 7 ,

дав (а, в = 1, 2) — компоненты метрического тензора,

/ дх \2 + / ду \ 2 дх дх + ду ду / дх \ 2 + / ду

д11 \ дд1/ \ дд1/ ' д12 д21 дд1 дд2 дд1 дд2' д22 \ дд2/ \ дд2 7 — якобиан преобразования (10),

дх ду дх ду

7

дд1 дд2 дд2 дд1

Отметим, что условие (11) не аппроксимируется, тем самым не выписывается разностное краевое условие в узлах, лежащих на подвижном дне. Вместо этого в узлах выписывается разностное уравнение, аппроксимирующее уравнение Лапласа для потенциала. При получении этого разностного уравнения используется интегрально-интерполяционный метод [3], в котором интеграл по подвижной части дна от функции в левой части равенства (11) заменяется интегралом от известной функции в правой части этого равенства.

В расчетах использовались подвижные в вертикальном направлении криволинейные сетки, подстраивающиеся под границы области течения, при этом в области над наклонным дном шаг сетки в горизонтальном направлении был равномерным, а в области над горизонтальным дном возрастал по закону геометрической прогрессии. Построение сетки начинается с расстановки узлов х, одномерной сетки на отрезке АС. Пусть на отрезке АС находится N1 ячеек сетки, на ВС — Ыг ячеек. Тогда равномерная сетка на отрезке АВ имеет шаг Дх — АВ/(Ж1 — N). Первый шаг неравномерной сетки на ВС полагается равным Дх, а последующие шаги вычисляются по формуле

Д,-+1/2Х - Дх д3-^1-^, 3 - N1 — N + 1,... , N1 — 1. (12)

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

(1 - дМг)

ВС - Дх(-^.

1 — д

Найдя д, на основе формулы (12) определяем координаты узлов сетки на ВС:

х,-+1 - х, + Дхд3-(М1-Мг), 3 - N1 — N + 1,... , N1 — 1.

После того как будет построена одномерная сетка на отрезке АС, вычисляются координаты узлов хх™+/12 - (х3-13-2 ,уП1+2) двумерной сетки, при этом полагаем

+1 +1 пГ1 + Мх?1 ,Г+1)

П+1 — '31 у 31' ' ■

N2

хл,32 - х,1, у™1+1 - — Мх^^Н3-^-З2, 31 - 0,N1, 32 - 0,..., N2.

2

2. Результаты численного моделирования

Приведенные ниже результаты расчетов получены при следующих безразмерных значениях параметров:

xa = 1, AC = 35, AB = 30, xo = 3, b =1, ti = 10, ¿stop = 18.04, tend = 30,

при этом безразмерные значения линейных величин, в том числе амплитуд генерируемых волн и возвышения свободной границы, получены делением размерных величин на h0 — некоторую характерную глубину, а безразмерное время — умножением на g/h0, где g — ускорение свободного падения.

Проводились расчеты на последовательности сеток, показавшие что дальнейшее измельчение сетки с количеством ячеек N = 400, N2 = 40 практически не меняет результатов расчетов, поэтому именно на такой сетке была выполнена вся серия расчетов при варьировании трех параметров: максимальной толщины оползня T, его ускорения a и угла наклона плоского откоса в.

На рис. 2-4 показана форма свободной границы в конечный момент времени t = tend при изменении одного из варьируемых параметров и фиксированных значениях двух

0.06 0.04 0.02 0.00 -0.02 -0.04 -0.06

-0.08 ...................................

5 10 15 20 25 30 35 X

Рис. 2. Свободная поверхность в конечный момент времени: в = 10°, a = 0.06; 1 — T = 0.03, 2 — T = 0.09, 3 — T = 0.15

Рис. 3. Свободная поверхность в конечный момент времени: в = 10°, T = 0.09; 1 — a = 0.02, 2 — a = 0.06, 3 — a = 0.1

Рис. 4. Свободная поверхность в конечный момент времени: a = 0.06, T = 0.09; 1 — в = 5°, 2 — в = 10°, 3 — в = 15°

других. Из рис. 2 видно, что при в = 10° и а = 0.06 увеличение толщины оползня приводит к возрастанию амплитуд генерируемых волн. При этом скорость движения волн слабо зависит от толщины оползня. Для указанного угла наклона дна головная

А,

0.05

0.00

5 А,

4

0.10

3

2 0.05

1

0.00

А

0.05

0.00

0 0.03 0.06 0.09 0.12 0.15 Т 0 0.03 0.06 0.09 0.12 0.15 Т 0 0.03 0.06 0.09 0.12 0.15 Т

а б в

Рис. 5. Зависимость амплитуды первой волны от толщины оползня: а — в = 5°, б — в = 10° в — в = 15°; 1 — а = 0.02, 2 — а = 0.04, 3 — а = 0.06, 4 — а = 0.08, 5 — а = 0.1

0.10

0.10

Таблица 1. Зависимость амплитуды первых двух волн и глубины впадины между ними от толщины оползня и его ускорения при в = 5°

Т а ¿1 Н12 ¿2

X П X П X П

0.03 0.02 27.714 0.0013 17.624 -0.0020 14.315 -0.0003

0.04 27.962 0.0040 19.609 -0.0078 15.639 0.0024

0.06 28.706 0.0087 20.932 -0.0206 16.549 0.0222

0.08 30.113 0.0139 22.255 -0.0517 19.195 0.0521

0.1 32.124 0.0176 25.563 -0.0545 21.015 0.0272

0.06 0.02 27.631 0.0030 17.576 -0.0049 13.902 -0.0011

0.04 27.714 0.0096 19.03 -0.0181 15.06 0.0036

0.06 28.541 0.0214 20.518 -0.0442 16.218 0.0624

0.08 30.029 0.0332 22.007 -0.1097 19.03 0.0989

0.1 32.2124 0.0395 25.398 -0.109 20.766 0.0453

0.09 0.02 27.383 0.0052 17.127 -0.0088 12.744 -0.0019

0.04 27.383 0.0177 18.203 -0.0320 14.15 -0.0039

0.06 28.458 0.0387 18.037 -0.0723 15.97 0.1253

0.08 30.113 0.0578 21.842 -0.1717 18.782 0.1337

0.1 32.224 0.0655 25.315 -0.1634 20.436 0.0534

0.12 0.02 27.052 0.0082 16.549 -0.0136 11.504 -0.0012

0.04 27.135 0.0291 16.135 -0.0520 12.082 0.0569

0.06 28.541 0.0596 17.541 -0.1407 15.721 0.2176

0.08 30.282 0.0866 21.594 -0.2345 18.451 0.1496

0.1 32.427 0.0954 25.315 -0.2178 20.105 0.0506

0.15 0.02 26.473 0.0127 13.737 -0.0207 9.8494 0.0025

0.04 26.969 0.0438 12.744 -1.1487 11.752 0.1465

0.06 28.706 0.0831 17.127 -0.2245 15.473 0.3044

0.08 30.54 0.1184 21.345 -0.2925 18.12 0.1457

0.1 32.632 0.1292 25.233 -0.2714 19.774 0.0386

волна покинула расчетную область, поэтому в окрестности правой границы свободная поверхность имеет впадину, и первая после нее волна имеет амплитуду А1, меньшую чем амплитуда А2 второй волны, движущейся за первой. Такой же монотонный характер имеет зависимость амплитуды А1 первой волны от толщины оползня и при других значениях ускорения и угла наклона, что видно из рис. 5.

А,

0.10 -

0.05 -

А,

0.10 -

0.05 -

А,

0.10

0.05

-0.05 4-■-1-■-1-■-1-■-1-■-1-■ -0.05 4-■-1-■-1-■-1-■-1-■-1-■ -0.05 4-■-.-■-.-■-.-■-.-■-.-■

0.00 0.02 0.04 0.06 0.08 0.10 а 0.00 0.02 0.04 0.06 0.08 0.10 а 0.00 0.02 0.04 0.06 0.08 0.10 а

Рис. 6. Зависимость амплитуды первой волны от ускорения: а — в = 5°, б — в = 10°, в — в = 15°; 1 — Т = 0.03, 2 — Т = 0.06, 3 — Т = 0.09, 4 — Т = 0.12, 5 — Т = 0.15

0.00

0.00

0.00

а

в

Таблица 2. Зависимость амплитуды первых двух волн и глубины впадины между ними от толщины оползня и его ускорения при в = 10°

Т а ¿1 Н12 ¿2

X П X П X П

0.03 0.02 17.489 0.00027 14.506 -0.00029 12.186 0.00028

0.04 19.312 0.0012 15.086 -0.0016 12.766 0.0007

0.06 20.804 0.0034 16.163 -0.0044 14.009 0.0080

0.08 20.886 0.0104 17.738 -0.0150 15.418 0.0162

0.1 23.041 0.0210 19.229 -0.0189 16.412 0.0183

0.06 0.02 17.406 0.00061 14.506 -0.00057 12.103 0.00065

0.04 19.146 0.0026 14.921 -0.0034 12.518 0.0017

0.06 20.638 0.0070 15.915 -0.0097 13.926 0.0181

0.08 20.721 0.0225 17.738 -0.0313 15.335 0.0342

0.1 22.958 0.0437 19.229 -0.0383 16.329 0.0373

0.09 0.02 17.323 0.0010 14.423 -0.0009 12.02 0.0011

0.04 18.981 0.0043 14.838 -0.0055 12.435 0.0033

0.06 20.472 0.0108 15.749 -0.0165 13.843 0.0308

0.08 20.638 0.0364 17.655 -0.0492 15.335 0.0541

0.1 22.875 0.0680 19.146 -0.0580 16.246 0.0569

0.12 0.02 17.241 0.00149 14.34 -0.0012 11.855 0.00152

0.04 18.732 0.0062 14.672 -0.0079 12.186 0.0058

0.06 20.306 0.0145 15.583 -0.0253 13.76 0.0465

0.08 20.555 0.0526 17.572 -0.0686 15.252 0.0759

0.1 22.792 0.0942 19.064 -0.0779 16.246 0.0770

0.15 0.02 17.158 0.0021 14.258 -0.0015 11.689 0.0020

0.04 18.484 0.0085 14.4423 -0.0106 11.938 0.0093

0.06 20.141 0.0179 15.418 -0.0363 13.678 0.0656

0.08 20.389 0.0713 17.489 -0.0895 15.25 0.0991

0.1 22.709 0.1223 18.981 -0.0977 16.163 0.0971

На рис. 3 показана зависимость формы свободной поверхности от ускорения при 9 = 10° и Т = 0.09. При малом ускорении оползень проходит относительно небольшое расстояние, поэтому генерируемые волны имеют небольшую амплитуду. При увеличении ускорения генерируются волны значительно более высокой амплитуды. При этом увеличение ускорения приводит также к увеличению скорости движения волн. Кроме того, при большом значении ускорения амплитуда А1 может стать больше, чем А2. Как и при изменении толщины оползня, при увеличении ускорения наблюдается монотонный рост амплитуды первой волны при разных значениях 9 и Т, что демонстрирует рис. 6.

Из рис. 5 и 6 видно, что для малого значения угла наклона 9 = 5° при увеличении толщины оползня или его ускорения происходит значительный рост амплитуды первой волны, в то время как при относительно большом угле 9 = 15° монотонный характер роста А1 хотя и сохраняется, однако амплитуда первой волны при разных Т и а изменяется в небольшом диапазоне значений. Это происходит по причине того, что глубина расчетной области зависит от угла наклона области, и в конечный момент времени оползень находится на большей глубине при 9 = 15°, чем при 9 = 5°.

Зависимость свободной поверхности от угла наклона 9 при а = 0.06, Т = 0.09 показана на рис. 4. На этом рисунке видно, что если при 9 = 10° и 9 = 15° головная волна

Таблица 3. Зависимость амплитуды первых двух волн и глубины впадины между ними от толщины оползня и его ускорения при в = 15°

т а ¿1 Н12 ¿2

X П X П X П

0.03 0.02 22.05 0.00021 15.337 —0.000001 12.934 0.00023

0.04 22.298 0.0009 15.917 —0.0004 13.265 0.0007

0.06 23.044 0.0021 17.409 —0.0017 14.343 0.0035

0.08 23.044 0.0039 18.403 —0.0052 15.586 0.0076

0.1 24.205 0.0090 19.729 —0.0088 16.58 0.0104

0.06 0.02 18.569 0.00048 15.088 0.00005 12.934 0.00053

0.04 22.216 0.0019 15.917 —0.0009 13.265 0.0015

0.06 22.961 0.0043 17.326 —0.0036 14.26 0.0075

0.08 22.961 0.0082 18.321 —0.0110 15.586 0.0161

0.1 24.205 0.0186 19.729 —0.0181 16.58 0.0217

0.09 0.02 18.569 0.0008 15.088 0.0001 12.851 0.0009

0.04 22.133 0.0030 15.834 —0.0014 13.182 0.0024

0.06 22.961 0.0068 17.243 —0.0055 14.26 0.0122

0.08 22.879 0.0129 18.238 —0.0172 15.503 0.0254

0.1 24.122 0.0290 19.647 —0.0279 16.497 0.0338

0.12 0.02 18.652 0.0011 15.088 0.0001 12.851 0.0012

0.04 22.05 0.0042 15.834 —0.0020 13.1 0.0035

0.06 23.879 0.0095 17.077 —0.0079 14.177 0.0176

0.08 22.713 0.0181 18.238 —0.0240 15.503 0.0358

0.1 24.122 0.0402 19.647 —0.0382 16.497 0.0468

0.15 0.02 18.652 0.0015 15.088 0.0002 12.851 0.0016

0.04 21.967 0.0055 15.751 —0.0026 12.934 0.0048

0.06 22.796 0.0124 16.995 —0.0101 14.094 0.0238

0.08 22.63 0.0237 18.155 —0.0314 15.42 0.0473

0.1 24.039 0.0522 19.564 —0.0491 16.414 0.0608

покидает расчетную область и первой волной Ai считается волна, следующая за ней, то при в = 5° это не так: первая волна и головная при этом угле наклона совпадают. Кроме того, можно сделать вывод, что больший угол наклона приводит к большей скорости распространения волн.

Результаты численных экспериментов сведены в табл. 1-3, в которых указаны амплитуды первых двух волн (Ai и A2) и глубина впадины между ними (Hi2). Кроме того, для каждой из этих величин указано расстояние от начала координат. Эти таблицы могут быть полезны для количественного сравнения полученных нами результатов с данными экспериментов или расчетов по другим математическим моделям.

Список литературы

[1] Шокин Ю.И., Чубаров Л.Б. О подходах к численному моделированию оползневого механизма генерации волн цунами // Вычисл. технологии. 2006. Т. 11. Спец. выпуск, посвященный 85-летию со дня рожд. акад. Н.Н. Яненко, ч. 2. С. 100-111.

[2] Елецкий С.В., Майоров Ю.Б., Максимов В.В., Нуднер И.С., Федотова З.И., ХАжоян М.Г., ХАкимзянов Г.С., Чубаров Л.Б. Моделирование генерации поверхностных волн перемещением фрагмента дна по береговому склону // Мат. междунар. конф. "Вычислительные и информационные технологии в науке, технике и образовании". Ч. 2. Алматы, 2004. С. 194-206.

[3] ХАжоян М.Г., ХАкимзянов Г.С Численное моделирование взаимодействия поверхностных волн с подводными препятствиями // Вычисл. технологии. 2003. Т. 8, № 4. С. 108-123.

Поступила в редакцию 19 марта 2007 г.