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

1 3 МА'Л

На правах рукописи

ТЕН Аркадий Алексеевич

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ СДВИГОВЫХ ДЕФОРМАЦИЙ НЕОДНОРОДНЫХ ГОРНЫХ ПОРОД

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

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

Новосибирск -1997

Работа выполнена и Институте минералогии и петрографии СО РАН

Научные руководители:

Официальные оппоненты:

Ведущая организация:

доктор физико-математических наук В.И. Кузин

кандидат геолого-минералогнческих наук В.В. Хлестов

доктор физико-математических наук Доровский В.Н.

кандидат физико-математических наук Х.Х. Имомназаров

Институт вычислительных технологий СО РАН (г. Новосибирск)

Зашита диссертации состоится «20» мая 1997 года в 1600 на заседании диссертационного совета Д.002.10.02 при Вычислительном центре СО РАН по адресу: 630090, г. Новосибирск, пр. Академика Лаврентьева, 6

С диссертацией можно ознакомиться в библиотеке Вычислительного центра СО РАН.

Автореферат разослан «18» апреля 1997 года Ученый секретарь

диссертационного совета, Г.И.Забиняко

ст*

ОБЩАЯ ХАРАКТЕРИСТИКА РЛВОТЫ

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

Цель работы. Создать алгоритм для численного двухмерного моделирования нестационарного ползущего течения среды (горной породы) учитывающий:

• Степенную реологию среды;

• Неоднородность среды с высокими градиентами вязкости, позволяющую моделирование полифазных реологических потоков;

Выявить основные закономерности поведения физических полей и структурной эволюции включений в деформируемой нелинейной неоднородной среде.

Фактически!V материал и методы исс.кдоттт. Основой для физической постановки задачи являются реологические экспериментальные исследования на используемых в технике и представленных в природе материалах при малых скоростях деформаций. Реализация полученной математической модели осуществляется с помощью конечно-разностных методов второго порядка по пространству к четвертого по времени. С целью уменьшения эффектов численной диссипации при интегрировании по времени в областях высоких граднентов опробовались методы расщепления с применением уплотненных сеток для высокоградиентных полей и обратного хода по характеристикам с использованием сплаиновой- и вейвлет-интерполяцпп. Защищаемые положении.

■ Создана, оттестирована и опробована численная схема для расчета деформаций в неоднородных по вязкости нелинейных средах с высокими градиентами в стоксовом приближении.

■ Деформации неоднородных сред приводят к аномальной локализации давления и напряжения, величина которых и расположение зон локализации зависит от расположения неоднородностей

■ Структурная эволюция включений в существенной мере зависит от реологии среды и. в частности, в нелинейной жидкости носит осциллирующий характер.

Научная нотпна. Исходная нелинейная задача решается с применением метода установления, который простыми средствами позволяет регуляризовать численную задачу. Применение метода переменных направлений, с другой стороны, позволяет свести разностную задачу к решению серии систем линейных уравнений с трехдиагональной матрицей.

С помощью метода расщепления реализован эффективный алгоритм второго порядка по времени для решение уравнения переноса, требующий решения только двух разностных уравнений.

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

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

Применяя программный интерфейс РУМ и диалекты Фортрана для параллельных вычислений, полученный алгоритм был адаптирован для использования на параллельных вычислительных комплексах с различной архитектурой.

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

1 Панин В.Е.и др. Эффект локализации деформации у границ зерен при ползучести поликристаллов//ДАН - 1990,- т.ЗЮ. - №1. - с.78-83.

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

Публикации. Основные результаты диссертации опубликованы в 16 работах, список которых приведен в конце автореферата.

Апробации работы. Основные результаты были доложены и обсуждены на конференции Европейского геофизического союза (Эдинбург, Англия, 1992), конференции Европейского союза наук о Земле (Страсбург, Франция, 1993), конференции Европейской исследовательской группы высоких давлений (Белфаст, Англия, 1993), тектоническом совещании: «Тектоника и стрессметаморфизм» (Москва, 1994), 16 ежегодном совещании Международной минералогической ассоциации (Пиза, Италия, 1994), международной конференции: «Структура и свойства высокодеформированных зон в горных породах» (Вербания, Италия, 1996), международной конференции Американского геофизического союза (Сан-Франциско, США, 1996), а также в ряде семинаров Объединенного института геологии, геофизики и минералогии СО РАН и Университета Миннесоты (США), Института систем информатики СО РАН, ВЦ СО РАН, Института теплофизики СО РАН.

/. Структура и Объем работы. Диссертация состоит из введения, трех глав, заключения и списка литературы, состоящем из 127 наименований. Текст работы содержит 127 страниц, 27 рисунков.

СОДЕРЖАНИЕ РАБОТЫ

Введение,. Изложено состояние вопроса на текущий момент. Сформулирована цель диссертационной работы, показана актуальность поставленных задач. Описана структура диссертации.

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

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

s

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

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

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

Е

¿„»А ехр(-— )aj¡. {1)

где - тензор скоростей деформаций, си- тензор напряжений, п - степень реологического закона, Е -энергия активации крипа, R - универсальная газовая постоянная, Т - абсолютная температура, А - материальная константа, которая, вообще говоря, в предположении неоднородности среды является функцией зависящей от пространственных координат. Данное соотношение принято за базовое в данной работе. $2. Общая формулировка задачи. В параграфе производится условная лианеризация реологического закона и формулируется начально - краевая задача. Реологическое соотношение (1) является нелинейным, однако его можно преобразовать в терминах нелинейной вязкости. Выразим (I) относительно напряжения о:

ст„=А "ехр(--)ёу (2)

оч = А"»ехр(—)ё«~б„. („

'пЯТ

Предполагая что среда является изотропной, выражение (2) можно преобразовать к виду:

_Е_

Тогда предполагая, что эффективная вязкость может быть выражена как:

-1 Е 1-. П = А пехр(^)ё" ,

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

сту=т1ёи- (5)

Такая, условно линейная формулировка позволяет применить хорошо разработанные математический и численные методы гидродинамики для построения моделей неньютоновских сред. Учитывая вышеизложенное, полная безразмерная система может быть записана в стоксовом приближении:

0 = -ёгас1(р)+Цл)У, (б)

где Цп) в плоском случае определяется как:

( д д д д д д

2-г—(Л——) (Л^—) (ЛТ-)

Эх, Эх, Эх2 дк2 Эх2 Эх,

8 д д д д д ЪГ^кг} аГ^ь? +2 э^7(11 эЙ

Шу(У) = 0. (8)

где V - вектор скорости, р - давление, х - пространственные координаты, р -плотность. В предположении, что свойства лагранжевой точки не меняются во времени, можно установить взаимнооднозначное соответствие между плотностью и материальной константой. Это означает, что если мы будем

■V)

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

р(х,0) = ро(х) А(х) = ф(р(х)) п(х) = ф(р(х)) . (10) Е(х) = п(р(х))

Система уравнешш(4),(б)-(10) является замкнутой. Для ее решения необходимо добавить начальные условия доя плотности и вязкости:

р(Х,0) = Ро(х) А(х,0) = А0(х) п(х,0) = п0(х) , (и) Е(х,0) = Е0(х)

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

невозмущенного простого сдвига с линейным профилем скорости:

= ["'

^mcuuw.rpeaei» ""|0[ ^^

V я V =

левля.гралпи лр»м_гр«|>ша

$3. Оиенки генерируемого давления и рамки применимости модели. Параграф посвящен оценке генерируемого давления и условиям применимости модели. Модель построена в предположении континуальности среды, поэтому

2y-i О

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

-' Е ¿о !.,

где с!о • характерный масштаб системы, хо - температуропроводность. Для того, чтобы не происходило существенной аккумуляции дисснпатнвного тепла, необходимо потребовать М<1. Тогда, используя параметры характерные, например, для кислых пород, получаем ограничение на мощность сдвигаемого слоя:

с!0 <3.8*Ю10ст.

что позволяет применять модель без ограничений при среднетектонических скоростях деформаций. Исходя из уравнения движения (6), можно получить оценку на давление:

-I Е ' Р«А пехр(-—)ёп пК.1

Глава 2. Численные методы решения.

$4. Аппроксимация уравнений. Рассмотрена аппроксимация исходных уравнений (б)-(9), которая производилась двумя способами. Первый -наиболее простой и наглядный. Для аппроксимации используются две несмещенные сетки разной плотности

П„:х:х = ¡Ь;у = ,)'Ь;0 ^ 1 £ N¡0 j 5 М

П;„:х:х = (21 - 1)Ь;у = 1)Ь;05! ^ 2М;0 5 ] * 2М

Сетка П), совпадает в узлах с нечетными узлами сетки Лзь- Компоненты скорости и давления заданы в узлах Пь, а кинетические параметры (вязкость, плотность и т.д.) определяются на П2Ь. Для получения конечно разностных

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

узлы для производных вида —. И, хотя использование интерполяции

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

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

Как видно, различие погрешности аппроксимации по некоторым членам довольно существенно, что для высокоградиентных сред приводит к значительному ухудшению сходимости, поскольку в остатке есть члены, пропорциональные производным по вязкости, и, в случае наиболее высоких градиентов, ошибка может достигать значительных величин. Таким образом, использование более плотной сетки для вязкости не меняет порядок аппроксимации, но улучшает ее качество. Члены высоких порядков производных в (13) существенно ниже, что дает лучшую аппроксимацию в области перегибов поля вязкости. Сравнение скорости сходимости показывает, что применение более плотной сетки для вязкости приводит к 10-50% - ому увеличению скорости сходимости в зависимости от геометрии градиентных зон и величин градиентов в них. Второй метод аппроксимации известен как метод ячеек. Компоненты скорости, давление и вязкость задается на разнесенных сетках. Практика показывает, что С-аппрокснмация, при всех равных прочих условиях, дает более высокую скорость сходимости нежели аппроксимация, рассмотренная в первом

д ди

(13)

(14)

методе. В литературе это хорошо известно для случая уравнений с конвективными членами, однако в силу того, что давление н дивергенция скорости в методе, примененном в данной работе, связаны, и поле давления имеет высокие градиенты, повышение качества аппроксимации уравнения неразрывности также приводит к увеличению скорости сходимости. Необходимо отметить, что уравнения сохранения импульса также аппроксимируются на плошади в два раза меньшей, нежели в первом методе. Как и в предыдущем случае, трудно дать точную количественную оценку, поскольку увеличение скорости сходимости зависит о расположения и величин генерируемых давлений и колеблется в пределах 7-35%. Одним из недостатков метода является его меньшая прозрачность, что приводит к большим временным затратам прп сопровождении программы. Кроме того, смещенные сетки требуют аппроксимации граничных условий. В целом это выражается в том, что алгоритм становится менее гибким, что приводит к трудностям при модификации алгоритма на многопроцессорных вычислительных системах и его модернизации.

£5. Численный метод решения уравнения движения. В параграфе рассматривается итерационный метод решения разностных уравнений движения. Система решается методом установления с введением искусственной сжимаемости. Метод показал хорошие результаты, является гибким и успешно реализуется для параллельных вычислений. Для этого в исходные квазистацнонарные уравнения движения и несжимаемости вводится фиктивная производная по итерационному времени т:

ЭУ

~ = -ёгас1(р)+Цг1)У. (15)

-1 Е 1-1 Л = Апехр(—)ё» _

ЭР

— = а сКу(У)

(16)

(17)

Векторное уравнение (15) решается методом переменных направлений. Для этого вся система расписывается на трех слоях итерационного времени. На одном слое записываются производные компонент скорости в направлении

х ма верхнем временном слое, в другом - в направлении у. Это позволяет привести исходную разностную задачу к серии систем линейных уравнений с диагональной матрицей. При достижении стационарного решения (15)-(17) вырождаются в исходные уравнения (6)-(8).

Система (15)-(17) формирует итерационный процесс для нахождения скорости и давления. Полный итерационная схема строится следующим образом:

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

2. Используя поле скорости, находится распределение эффективной вязкости в системе.

3. Подставив вычисленное приближенное значение скорости в уравнения (17), находится новое приближенное поле давления.

4. Если условие сходимости соблюдено, выполняется шаг 5

Иначе выполняется шаг 1.

5. Используя вновь вычисленные значения поля скорости, находится распределение вязкости, плотности и других необходимых параметров на следующем временном слое. Расчетное время увеличивается на шаг по времени.

6. Если необходимое расчетное время достигнуто, то расчеты останавливаются.

Иначе выполняется шаг 1 ^'б. Численные методы решения уравнении переноса. Описаны и предложены различные методы решения уравнение переноса, показаны их преимущества и недостатки в применении к неоднородным сильноградиентным средам. Уравнение переноса (9) является одним из важнейших в данной системе. Именно его посредством реализуется нестационарность системы в целом, поэтому к нему предъявляются повышенные требования по точности решения, поскольку, как известно, ошибка при продвижении по времени накапливается и не компенсируется по ходу движения. Кроме того, как уже упоминалось ранее, переносимые поля обладают высокими градиентами, а следовательно, метод решения должен быть устойчив к численной диффузии.

Для решения был опробован метод растепления с экстраполяцией скорости на верхний временной слои. Уравнение записывается на трех временных слоях к, к+1/2, к+1. Экстраполяция скорости требуется на слон к+1/2 и может быть проведена с помощью линейных методов. В отличии от использованной Н.И.Мартыновым схемы2, требуется решение только двух уравнений вместо четырех. Порядок аппроксимации 0(Дх!,Д1:). Полученная неявная схема безусловно устойчива. К недостаткам метода можно отнести то, что в вычислениях используются высокоградиентные поля, что приводит к существенной численной диффузии.

Методами, которые позволяют избежать прямого влияния градиентов поля плотности, являются методы характеристик. В данной работе был использован метод обратного хода по характеристикам. Суть метода состоит в нахождении координаты узловой точки текущего временного слоя к на предыдущем к-1 временном слое. Для этого решается уравнение:

^ -V

Л • (|9)

назад по времени. После нахождения координаты точки на предыдущем временном слое хы с помощью интерполяции находится значение плотности ры или других кинетических параметров в ней и значения их переносятся в узлы текущего временного слоя (рк=ры). Для решения уравнения (19) использовался метод Рунге-Кутта четвертого порядка по времени. После нахождения координаты хы, которая, вообще говоря, не совпадает с узлами сетки, необходимо найти значение плотности в ней. Поскольку поле плотности может содержать большие градиенты, точность интерполяции фактически определяет точность решения уравнения (18) в целом. Линейная н билинейная интерполяции приводят к быстрому размыванию градиентных зон и демонстрируют точность хуже, нежели метод расщепления с аппроксимацией скорости. Интерполяция с использованием бикубических сплайнов более устойчива в градиентных

' Мартынов Н.И. Эволюция поверхности раздела при тейлоровской неустойчивости двумерного ползущего движения: Днс.канд. фнз.-мат. Наук. -Алма-Ата, 1985,- 107 с.

~r

0.00

4,00

~Г~ 8 00

I

12.00

Рисунок 1. Базовая функция и функция порожденная смещение и сжатием базовой (пунктирная линия)

полях3. Точность,

обеспечиваемая сплайнами, позволяет хорошо

сохранять градиенты исходного поля (Ten at al..l996 b,d). Тем не менее, существование базиса аппроксимации накладывает принципиальное ограничение на величину градиента, который может быть аппроксимирован.

Разрешить данную проблему можно введением адаптивных сеток или применением вейвлет-интерполяции. Одним из важнейших, в данном случае, свойств функционального базиса, построенного ка вейвлетах (wavelet - название принято из зарубежной литературы), является отсутствие локального носителя фиксированного масштаба. Таким образом, функциональное пространство построенное на вейвлетных функциях (Рисунок 1), не имеет характерного размера, а следовательно, теоретически может отражать детали с любой степенью разрешения. В основе вейвлстной интерполяция'* лежит фрактальный принцип - принцип вложенности (Ten et al., 1996с, Ten et al. 1997), что делает их похожими, в некотором смысле, на Фурье - преобразование. В тоже время вейлет-интерполяция обладает свойством локальности, которое не присуще приближению рядами Фурье. Таким образом, мультнмасшабностъ и локальность создает функциональное пространство, которое обладает лучшими свойствами рядов Фурье и сплайновой аппроксимации. Функциональный базис вейвлет - интерполяции порождается сжатием и сдвигом базовой функции (сплошная линия Рисунок I). Таким образом,

1 Наймарк Б.М., Малевскин A.B. Экономичный метод бикубической спланн-интсрполяции/ДПогичсские и вычислительные методы в сейсмологии. - М:

Наука, 1984. С.141-149

' Grossmann Л., Morlet J. //SIAM J. Math. Anal. - 1984. - v. 15 - p.723

потенциально таком базис имеет бесконечный спектр и определен, и данном случае, на всей действительной осн.

■У 7. Испо.уоояиние многопроцессорных вычислительны, у комюексан. Сформулированы основные концепты организации параллельных вычислений и их применение к задаче диссертации. Прозрачность схемы и достаточная простота аппроксимации делают алгоритм достаточно удобным для эффективного использования на многопроцессорных вычислительных комплексах. Алгоритм был модифицирован для использования как на комплексах с распределенной памятью (distributed memory model), так и на комплексах с разделяемой оперативной памятью (shared memory model). Для организации параллельных вычислений на многопроцессорных комплексах SGI (Power Challenge) использовался многопроцессорный диалект Фортрана, позволяющий параллельные вычисления в циклах. В рамках распределенной вычислительной среды реализация параллельных вычислений проводилась с помощью программного интерфейса PVM (parallel virtual machine), который позволяет использовать последние модификации Cray T3D , Cray ТЗЕ, кластерные комплексы на базе IBM RS6000/SP2 и, вообще говоря, любой набор компьютеров, объединенных тем или иным образом в компьютерную сеть (в том числе распределенная вычислительная сеть может быть создана на базе IBM PC совместимых компьютеров).

»'8. Тесты модели. Тестирование модели производилось на аналитических решениях. В качестве первого теста использовалось аналитическое решение обтекания цилиндра в бесконечной ньютоновской среде. В модели задавалось цилиндрическое включение с вязкостью, в 20 раз превышавшей вязкость вмещающей среды. За пределами градиентной зоны ( переход от вязкости среды к вязкости включения) различие составляло порядка квадрата пространственного шага. Это соотношение соблюдалось при дроблении сетки. В нелинейной области тестирование проводилось на аналитическом решении движения степенной жидкости между двумя пластинами под действием градиента давления. Сравнение численного и аналитических решений показало их совпадение также с точностью порядка квадрата шага.

Рисунок 2. Поле давления с наложенным векторным полем скорости вокруг одиночного включения при рачлнчноН реологии среды, а) п=| 6) п-3 в) п=5 г) п=7

имсщаюшеЛ среды к вязкости включения осуществлялся в узкой

градиентной зоне, составлявшей менее 3% от радиуса включения.

Еоиинчное жесткое включение посвящен исследованию поведения физических полей и эволюции одиночного жесткого включения. Под жестким включением подразумевается включение, вязкость которого выше вязкости вмещающей среды. Сдвиговый поток формирует неоднородное поле давления вблизи включения (Рисунок 2). Зоны максимального отклонения давления локализованы близ неоднородности. С ростом степени реологического закона происходит рост эффективного контраста за счет падения сдвиговых напряжении во включении, что, в свою очередь, приводит к тому, что деформация включения приближается к вращению тердого тела. Песмофя на ю. ч то эффективная вязкость вблизи включения

| гиЧ

— М м, джл«им

• Ммсммшмы* направим • матриц« - «. - МиМШ/ЬиММЛриММЯаОМЛе'МИИИ

~г-600

Рисунок 3. Поведение давления и сдвиговых напряжений в зависимости от реологии среды.

налает. амплитуда

генерируемых давлений остается относительно стабильной, так же как н касательные напряжения (Рисунок 3).

В У/0. Единичное мягкое включение рассматривается мягкое включение, вязкость которого ниже

эффективной вязкости вмещающей среды.

Генерация давления и локализация напряжений

происходит и в том случае, если среда содержит мягкое включение. Одной из важных особенностей формируемого течения является наличие существенного переноса вещества в поперечном к сдвигу направлении, которое происходит на мягком включении. В целом, структурно, поле давления аналогично жесткому включению (Рисунок 2), но является инверсным, то есть, там, где для жесткого включения амплитуды положительны, для мягкого - отрицательны. При вязком контрасте, равном 0.5, амплитуды генерируемого давления несколько ниже, чем в случае с жестким включением и вязком контрасте, равном 2.

$11. Пара шаиндрических включений. В параграфе рассмотрены особенности взаимодействия включений друг с другом. Присутствие в системе второго включения изменяет структуру потока и, соответственно, формируемые поля. Взаимодействие включений увеличивает генерируемое давление. Зона максимального повышенного давления располагается между включениями, и давление достигает почти 5 безразмерных единиц. Отмечается, что амплитуда давления зависит от структурных особенностей системы, в частности, от ориентации днпольной системы включений, а также от относительного расстояния между включениями. Поведение дипольной системы приближается с уменьшением расстояния между включениями к поведению одиночного включения (Рисунок 4 на следующей странице).

Рисунок- 4. Функция тока прн различном расстоянии включений друг от друга. Результаты моделирования позволят предположить, что значимое взаимодействие включений прекращается примерно на расстоянии М = 3.5, где 2 - расстояние мезду центрами включений и с! - их диаметры. Это, в свою очередь, приводит к заключению, что при фазовой плотности включений менее 0.125 система может рассматриваться как разбавленная (невзаимодействующие включения).

(12. Гипотеза о генерации высокого давления В параграфе рассматривается один из возможных аспектов применения результатов моделирования. Рассмотрены современные гипотезы, объясняющие происхождение метаморфических комплексов сверхвысоких давлений и низких температур (ВДНТ). Численное моделирование, свидетельствующее о том, что прн деформации неоднородной среды могут возникать динамические, локализованные во времени и пространстве, зоны повышенного давления, позволяет сформулировать гипотезу возникновения ВДНТ за счет деформаций в нижней части земной коры.

Заключение. В заключении сформулированы основные результаты диссертации:

1. Предложен и оттестирован алгоритм для расчета деформации неоднородной нелинейной среды с произвольным расположением и

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

2. Решение уравнения переноса методом обратного хода по характеристикам с интерполяцией на базе бикубических сплайнов является наиболее эффективным в смысле отношения количество операций/точность решения, в то время как вейвлет-ннтерполяцпя более эффективно подавляет численную диффузию.

3. Алгоритм адаптирован с применением параллельного диалекта Фортрана и программного интерфейса PVM к использованию па многопроцессорных вычислительных машинах и в распределенной вычислительной среде.

4. Деформация неоднородной среды приводит к динамической генерации давления в локализованных областях. Величина избыточного давления не испытывает существенных изменений при изменении реологии среды благодаря тому, что происходит компенсация падения эффективной вязкости с ростом степени реологического закона за счет локализации потока и роста эффективного вязкого контраста.

5. Показан осциллирующий характер эволюции включения в сдвиговом нелинейном потоке.

6. Показано, что S - образные структуры формируются только в случае нелинейной среды с включением, эффективная вязкость которого выше эффективной вязкости вмещающей среды.

7. Результаты моделирования объясняют происхождение структур, известных под названием «тени давления», за счет формирования областей разрежения вблизи включения.

8. Предлагается гипотеза генерации избыточного давления в земной коре на неоднородностях при ее деформации.

Основные результаты диссертации опубликованы в работах:

1. Теп A. Computer simulation of the rock shearing// Annales Gcophysicac -1992,-vlO-Supl l.-Pt. 1. - p. C98.

2. Тен А. Динамическая модель генерации высоких давлений при сдвиговых деформациях горных пород (результаты численного эксперимента)// ДАН - 1993.-т.328.-№3-с. 322-324

3. Ten Л. The computer model of shear zone in the earth's crust// Terra Nova, Abstract supplement. - 1993. - v.5. I. - p. 210

4. Ten A. A method of determing of rheologic charactristics// Annales Geophysicae - 1993 v. 14 - Supl 1 - Pi. I - p C434

5. Ten A. A new melhod of pressure generation// EHPRG XXXlst meeting, Belfast,-1993,-scc D.

6. Ten А., Лепетюха В.В. Стресс - метаморфическая модель образования высокобарических минеральных ассоциаций в земной коре// Тектоническое совещание тектоника и метаморфизм, тезисы докладов/ГФ МГУ. 1994. с. 22-25

7. Lepetyukha V.V., Reverdatto V.V., Ten A High pressure and ultrahigh-pressure metamorphic rocks in northen Kazakhstan// International mineralogical association, 16"1 General meeting Abstracts - 1994, p.239

8. Ten A A new dynamic model of high pressure genration in the earth's crust// International mineralogical association, 16,h General meeting Abstracts - 1994 -p.406-407

9. Лепетюха В.В., Ревердатто В.В., Тен А., Хлестов В.В Полиметаморфизм в северовосточной части Кокчетавского массива (Северный Казахстан)// ДАН - 1995 - т.340 - №5. - с. 653-658

10. Ten A, Yuen D., Podladchikov Yu.Yu., Larsen T.B., Malevski A.V. The evolution of material surfaces in convection with variable viscosity as monitored by a characteristics-based method// University of Minnesota Supercomputer Institute Research Report - 1996. UMSI 96/10 - p. 1-8.

11. Ten A, Yuen D., Podladchikov Yu.Yu., Larsen T.B., Malevski A.V. Chaotic mixing of non-newtonian and newtonian convection in the earths's mantle// University of Minnesota Supercomputer Institute Research Report - 1996. -UMSI 96/16 - p. 1-9

12. Ten A, Yuen D„ Podladchikov Yu.Yu., Larsen T.B., Pachepsky E„ Malevski A.V. Fractal features in mixing of non-newtonian and newtonian mantle convection// University of Minnesota Supercomputer Institute Research Report - 1996ю - UMSI 96/166 p. 1-19.

13. Ten A., Yuen D.A., Larsen T.B., Malevski A.V The Evolution of Material Surfaces in Convection With Variable Viscosity As Monitored By A

Characteristics - based method// Jour. Geophys. lies. Lett - 1996. - vol 23 - №16 ■p. 2001-2004

Ten A., Podladchikov Yu.Yu Stress distribution and structure evolution luring the deformation of inhomogeneous rocks// Uiversita degli studi di Viilano, An international conference on sturcture and properties of high strain :ones in rocks (Abstract volume). - 1996 - sup. 107 - plOO Y.Y. Podladchikov, A. Ten, D A Yuen, E Pachepsky, T B Larsen, A V dalevsky, Fractal Featutes in Mixing of non-Newtonian and Newtonian Mantle Convection //AGU Fall Meeting - 1996. - T42A-4 - p.F750 A. Ten, D.A. Yuen, Yu.Yu. Podladchikov, T.B. Larsen, E. Pachepsky, A. V. /lalevsky Fractal features in mixing of non-newtonian and newtonian mantle onvection// Earth and Planetary Science Letters - 1997 - vol 146 - issue 3-4 -401-414

Подписано к печати 16.04.97. эормат 60x84/16. Бумага офсет N 1. Гарнитура Тайме. Офсетная печать. Печ. л. 1,4. Тираж 70. Заказ148.

Новосибирск, 90, Университетский просп., 3 НИЦОИГГМСОРАН