Исследования долгопериодной динамики тепловыделяющих геофизических структур тема автореферата и диссертации по физике, 01.04.02 ВАК РФ

РГб од

РОССИЙСКАЯ АКАДЕМИЯ НАУК ИНСТИТУТ ТЕОРЕТИЧЕСКОЙ ФИЗИКИ им. Л.Д.ЛАНДАУ

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

БЯЛКО Алексей Владимирович

УДК 532.13

ИССЛЕДОВАНИЯ ДОЛГОПЕРИОДНОЙ ДИНАМИКИ ТЕПЛОВЫДЕЛЯЮЩИХ ГЕОФИЗИЧЕСКИХ СТРУКТУР

01.04.02 — теоретическая и математическая

фиоика

Диссертация на соискание ученой степени доктора физико-математических наук в форме научного доклада

Черноголовка — 1995

Работа выполнена в Институте теоретической физики им. Л.Д.Ландау РАН

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

член-корреспондент РАН, доктор физико-математических наук С.И.Анисимов

доктор физико-математических наук, профессор Л.А.Большов

доктор фиоико-математических наук, профессор Б.А.Трубников

Ведущее предприятие — Институт динамики геосфер РАН

Защита состоится 29 декабря 1995 г. в 12 ч. на заседании Специализированного совета Д 002.41.01 по защите диссертаций на соискание ученой степени доктора физико-математических наук при Институте теоретической физики им. Л.Д.Ландау РАН по адресу: 142432, Московская обл. Ногинский р-н, Черноголовка, ИТФ РАН.

С диссертацией можно ознакомиться в библиотеке Института теоретической физики им. Л.Д.Ландау РАН

Диссертация разослана " " 1995 г.

Ученый секретарь Специализированного совета

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

Л. А. Фальковский

ВВЕДЕНИЕ

Актуальность темы

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

Захоронение ядерных отходов — важная, но практически еще не решенная задача. Ее трудность состоит в необходимости предсказания поведения образующейся геологической структуры на временах вплоть до десяти миллионов лет. Поскольку периоды полураспада образующихся в реакторе нестабильных изотопов сильно различаются между собой, то зависимость от времени для их активности и тепловыделения оказывается не экспоненциальной функцией времени (как для каждого изотопа в отдельности), а близкой к степенному закону Проблема долговременной безопасности захоронений ядерных отходов связана с тем фактом, что показатель экспоненты а в этом усредненном степенном законе оказывается меньше единицы ( для отработанного ядерного топлива а равно 0.05, для ядерных отходов — 0.91). Суммарный потенциальный ущерб произведенных ядерных отходов есть интеграл но времени от их энерговыделения. Поскольку а < 1, то этот интеграл оказывается расходящимся. Конечно, физически нет необходимости интегрировать до бесконечности, так как существует уровень естественной радиоактивности. Конечным моментом для отходов можно считать тот момент, когда их тепловыделение сравняется с тепловыделением естественных месторождений урана и тория, а поверхностная радиация захоронений выйдет на уровень естественного фона. Это произойдет только через 10 млн. лет после момента захоронения, когда из всех изотопов радиоактивными останутся только протоактиний-231, торий-232 и два изотопа урана. Этот период относительно короток в геологической шкале времени, но велик для прямых экспериментов.

Рассмотрению всего комплекса вопросов, связанных с захоронениями ядерных отходов и динамикой образующихся геологических структур, посвящена монография автора [1].

Основные направления исследования

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

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

Представление горных пород очень вязкими жидкостями имеет тот недостаток, что заменяет непрерывными те перемещения, которые реально дискретны. Четвертый раздел диссертации, основанный на публикациях [1,6], посвящен теоретическому рассмотрению распределения дискретных смещений и выводу распределения распределения землетрясений по магнитудам (Гутенберга—Рихтера).

На основе результатов диссертации предложен ряд практических решений проблемы захоронения ядерных отходов, изложенных в монографии [1], трудах конференций [5,7,8], защищенных патентами [9-11] и доведенных до массового читателя [12,13]. Предметом данной диссертации являются теоретические разра-

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

Научная ценность и новизна

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

- установлена точность оценок — формул по порядку величины, не содержащих численных коэфициентов;

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

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

- решена задача о погружении тепловыделяющего шара в жидкости переменной вязкости — модификация задачи Стокса;

- найдена связь распределения дискретных перемещений по их величине с конкурентными распределениями. В предположении справедливости конкурентного распределения выведен закон Гутенберга—Рихтера для частотности землетрясений;

- получено решение задачи тепломассопереноса для режима потоковой загрузки последовательности тепловыделяющих тел через узкий вертикальный канал (скважину);

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

Апробация работы

Диссертация содержит результаты одной монографии, 7 научных работ и 3 патентов. На основе работ, представленных в диссертации, написаны три научно-популярные статьи. Работы опубликованны в России и за рубежом, а также докладывались на семинарах в ИБРАЭ РАН, ИФЗ РАН, ИД Г РАН, ВНИИНМ им. Бочвара, ЦФТИ МО, на кафедре минералогии МГУ, а также на международных конференциях по переработке ядерного топлива и обращению с отходами ИЕСОБ'Э! (Япония) и 11ЕС00'94

(Великобритания) и международных конференциях Ядерного общества (Москва, 1994, Обнинск 1994).

Основные результаты

настоящей работы получены в Институте теоретической физики им. Л.Д.Ландау РАН. Работы [7,9,10] выполнены в соавторстве с О.Б.Хаврошкиным (ИФЗ РАИ), соавтором работы [9] является также И.М.Халатников (ИТФ РАН).

Результаты работ представлены далее в следующих разделах:

1. Точность оценочных формул

2. Гидродинамика с переменными вязкостью и теплопроводностью

3. Гидродинамика погружения нагретого тела

4. Статистика дискретных смещений

5. Практические приложения

6. Выводы

1 Точность оценочных формул

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

1.1 Статистика физических коэффициентов

Если цифра — разновидность буквы, то числа — это слова. Статистика обычных слов по корпусам однородных текстов приводит к ряду закономерностей, из которых наиболее известен закон Ципфа: частотность слова (вероятность его появления) степенным образом зависит от его ранга — номера в последовательности, упорядоченной по частотности. Поэтому возникает естественный интерес к статистике чисел в физических и математических текстах [2,3].

Прежде всего надо заметить, что есть два вида чисел и их статистика принципиально различна. Есть числа, которые выражают измеренные величины в установленных системных еди-

ницах. Они, как правило, имеют размерность, а их ¡значение зависит от выбранной системы единиц. Вопрос о статистике всех размерных чисел представляется довольно абстрактным, но существует одно количественное следствие такой логарифмической широты. Это закон Ньюкомба—Бедфорда [14]: вероятность \У(п) того, что первая цифра размерного числа равна п, есть Ж(п) = 18((п + 1 )/п), где (п = 1.2...9).

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

Другой вид чисел — это коэффициенты в физических и математических формулах (например: 2,3,1/4, л/2, 7г, 1п2). Они возникают в результате чисто математических операций и поэтому принципиально безразмерны. Коэффициенты, как правило, не очень велики и не очень малы по абсолютной величине. Именно статистика коэффициентов и будет нас далее интересовать. Знак коэффициента будем считать несущественным. Отметим сразу, что закон Ньюкомба—Бедфорда окажется не применим к первым цифрам физических коэффициентов.

В качестве эталонного корпуса текстов был выбран "Курс теоретической физики" Ландау—Лифшица. Его 10 томов охватывают разнообразные разделы физики, он однороден по авторскому стилю и дает представительную статистику коэффициентов. Проводился учет коэффициентов во всех формулах, записанных в виде отдельных строк (вне текста), включая задачи и сноски. Технически это осуществлялось последовательной компьютерной записью коэффициентов без учета знака. Формулы со знаками "много больше", "много меньше" и "порядка величины" исключались из подсчета, поскольку они принципиально не содержат численных коэффициентов.

Коэффициенты, содержащие множитель тг и его степени, считались отдельно — каждый как единый коэффициент; так например, коэффициент 2тг не вносил вклад ни в статистику коэффициента 2, ни в статистику коэффициента тг. Аналогично этому, считались отдельно дробные коэффициенты с числителями, отлич-

ными от единицы — они не вносили вклада в статистику своих числителей и знаменателей. Производился раздельный подсчет взаимно обратных коэффициентов 2 и 1/2, 3 и 1/3, 4 и 1/4; все остальные коэффициенты подсчитывались вместе с им обратными. Тождественные формы записи одного коэффициента (например, \/8, 2\/2 и 23/2) не различались.

Частотности коэффициентов (вероятности их появления в корпусе текстов) оказались весьма упорядоченными, они с достаточной точностью повторялись от тома к тому. Анализ уже первых пробных подсчетов показал, что содержательную информацию представляет также общее число равенств, не содержащих коэффициентов. В дальнейшем вместе с обычными коэффициентами подсчитывался и не проявляющийся в алгебраической записи коэффициент 1 (единица). Суммарное (вместе с единицей) число всех коэффициентов в курсе Ландау—Лифшица равно N = 7.45 • 104.

При анализе экспериментальной статистики обращают на себя внимание несколько закономерностей.

1) Частотности тех обратных коэффициентов, которые подсчитывались, (1/2, 1/3, 1/4) оказываются примерно вдвое меньшими, чем частотность появления коэффициента вместе с ему обратным; ото и оправдывает совместный подсчет ддя коэффициентов к и fc-1.

2) Убывание частотности при последовательном удвоении коэффициента происходит по закону, близкому к степенному.

3) Частотности коэффициентов 3, 4, V2, я", 2ж, 4ж близки между собой по порядку величины.

4) При удвоении коэффициентов, кратных 3, (в последовательности (3, 6, 12, 24, 48) частотности падают по степенному закону с наклоном в двойных логарифмических координатах, близким к наклону степеней двойки (1, 2, 4, 8, 16, 32, 64). Так же ведет себя и последовательность (1, т/2, \/8, \/32)-

Эти эмпирические закономерности позволяют провести моделирование генезиса коэффициентов.

1.2 Модельное распределение коэффициентов

В первой работе по статистике коэффициентов [2] на основе этих закономерностей была сделана попытка осмысления результа-

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

3 = 2 + 1; 4 = 2*2; 3/2=1/2+1.

Их вероятности пропорциональны р2.

Рациональные коэффициенты третьего класса (с вероятностями, пропорциональными р3), образуются в семи процессах:

5 = 4 + 1; 6=2*3; 8=2*4; 4/3 = 2/3*2;

5/2 = 3/2+1 = 1/2 + 2; 5/4=1 + 1/4; 5/3 = 2/3 + 1.

Ясно, что любой рациональный коэффициент может быть образован сложением, умножением и обращением, причем множеством способов. На примере коэффициента 5/2 видно, что некоторые коэффициенты имеют несколько равновероятных путей образования уже в главном процессе. Кроме наиболее вероятных способов генерации коэффициента, существуют и более сложные пути (например, 2=3*2/3), они также вносят свой малый вклад, пропорциональный более высоким степеням элементарной вероятности р. В более высоких порядках каждый коэффициент имеют все возрастающее число способов образования. Поэтому для вычисления полной вероятности коэффициента надо взять сумму вероятностей по всевозможным путям, а затем полученные вероятности перенормировать. Сделаем следующий шаг: создадим компьютерную модель генерации коэффициентов, учитывающую эту вероятностную комбинаторику [3]. При этом

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

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

и>(к 1 + = ри>(к1)и)(к2)\ и>(&1&2) = гу(£-1)ш(&2)- (1)

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

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

Обратимся к корневым нррационалыюстям. Поскольку наша "экспериментальная" статистика показывает, что частотность коэффициента \/2 порядка р2, то исследуем справедливость простого общего принципа: вероятности корней уравнения равны вероятности образования самого уравнения. Начнем с вырожденного квадратного уравнения х2 — а = 0. Оно возникает в результате алгебраического сложения, поэтому его вероятность равна рги(а). Тогда в главном порядке по р получим:

ю(у/2) = 2рю(2) = 4 р2,

и>(\/3) = 2рш(3) = 8 р3, ю(\/8) = 2рш(8) = 32р4

и т.д.

Подсчитаем таким же образом вероятность полного квадратного уравнения х1 + ах + Ъ — 0, а следовательно, и его корней

ги(х2+ ах + Ь=0)=р2ю(а)и](Ь). (3)

Отметим, что из всех коэффициентов, возникающих из этого уравнения, наиболее вероятным оказывается "золотое сечение" — корни уравнения х"1 — х = 1.

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

Первым (и наиболее вероятным) уравнением с комплексными корнями является ж2 + 1 = 0. Формально вычисленная вероятность его корня х (мнимой единицы) оказывается равной р. Мнимая единица достаточно широко употребляется в промежуточных формулах теоретической физики, но производить ее статистический учет не имеет большого смысла, поскольку в физике и математике исторически сложилась тенденция по возможности избегать формул, содержащих мнимую единицу — достаточно вспомнить, что мнимая единица в неявном виде содержится в записи любой тригонометрической функции.

На комплексной плоскости процедура образования коэффициентов 7г, 2я- и Атг из г достаточно естественна, к ним приводят тождества Эйлера

— 1 = ехр(7гг); 1 = ехр(27гг); 1 = ехр(47и). (4)

Сравнение с "экспериментальной" статистикой приводит к заключению, что вероятность появления экспоненты также надо учитывать множителем р. Однако пока мы не стали вводить (или вычислять) вероятности элементарных функций (экспоненты,

логарифма, тригонометрических и т.д.), а формально задали вероятности первых коэффициентов, кратных тг, как

тг) = р2; т(2ж) = 2р2; w(4tt) = 2р2.

(5)

Это примерно соответствует экспериментальной статистике (рис.1), вероятности же остальных коэффициентов, содержащих 7г, вычислялись уже сформулированными правилами (сложения и умножения) с единственным ограничением. А именно: числа, кратные тг, нельзя складывать с рациональными числами.

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

I

ц-5

«V -\fl

mßfy' JL &

да. ^

-Л-—:-.-1-:-.

-5 V2 -4 -3 -2 -1 О

логарифм экспериментальной вероятности Ig

Рис.1. Сравнение экспериментальной частотности коэффициентов (ось ординат) с теоретической вероятностью (ось абсцисс) при значении параметра р = 0.06. Пунктирная прямая проведена методом наименьших квадратов.

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

w(kl'n) = w(k)w(n)p2, iü(ln(Jb)) = w(k)P2. (б)

Конечно, для таких классов коэффициентов гипотезы образо-

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

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

0.06. По оси ордипат отложен десятичный логарифм коэффициента, по оси абсцисс десятичный логарифм его теоретической вероятности. Каокдый вертикальный отрезок соответствует одному коэффициенту, всего в диапазоне к от 0.002 до 500 оказалось 1280 коэффициентов с вероятностями выше Ю-6, по на рисунке показана лишь часть распределения коэффициентов, меньших единицы, оно продолжается симметрично правой половине. На рисунке показана также гистограмма с шагом Ак = 0.11п 10, она находится в приближенном согласии с логнормальным распределением (штриховая кривая). Однако, высшие точки распределения спадают с возрастанием коэффициента не по логпормали, а близко к степенному закону.

1.3 Точность оценок

Модельное распределение коэффициентов (см. рис.2) выглядит довольно необычно и имеет характер фрактала: при увеличении масштаба проявляется все большее число коэффициентов с меньшими вероятностями, но общая картина выглядит примерно одинаково при любом увеличении. Однако это касается только коэффициентов с малыми вероятностями, а пики распределения (в основном степени двойки) ведут себя не самоподобно, а убывают с ростом коэффициента более регулярно, примерно степенным образом. Особой областью, также отклоняющейся от фрактального поведения, является центр распределения, окрестность единицы. Вблизи единицы распределение выглядит особым образом потому, что там оно дублирует в сжатом виде область очень малых коэффициентов (и очень больших): близкие к единице числа образуются в основном при сложении ¿2 = 1 + где ¿1 мало.

Посмотрим теперь на усредненную зависимость распределения от величины коэффициента. Обычный метод статистического анализа — построение гистограммы, т.е. суммирование всех отсчетов внутри равных интервалов. Очевидно, что в нашем случае интервалы должны быть логарифмическими равными. На графике показана такая гистограмма с логарифмическим шагом 0.11п 10. Выбор шага произволен, важно только, чтобы он не совпадал с логарифмами малых натуральных чисел. Однако при построении гистограммы нас подстерегает логическая опасность: сложение вероятностей, широко распределенных в логарифмической шкале, связано с существенной потерей информации (малые вероятности маскируются на фоне больших). Чем шире шаг гистограммы, тем хуже она отражает поведение большого массива малораспространенных коэффициентов. Чрезвычайно резкие скачки модельного распределения еще остаются заметны на гистограмме с малым выбранным шагом: с ростом коэффициента в среднем происходит убывание гистограммы, но оно не регулярно: отклонения от апрокснмнрующей ее параболы не нестолько малы, чтобы можно уверенно сказать, что распределение коэффициентов в среднем логнормально.

Посмотрим теперь, как выполняется для распределения коэффициентов закон Ципфа. Для этого расположим коэффициенты

в порядке убывания их вероятности, их номер в этой последовательности называется рангом. При этом происходит суммирование распределения для разных к. Результат этой операции неоднозначен: его можно интерпретировать и как стеленную зависимость ш ~ г1,8 (что соответствует закону Ципфа), и как параболу логнормального распределения.

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

<т2 =< (1п^)2 >= 0.70 (7)

(угловыми скобками обозначено усреднение). Это означает: пределах одного стандартного отклонения (с вероятностью 67%) неизвестный коэффициент находится в пределах от 0.43 до 2.3.

В этом результате (с учетом оговорки о его точности) не содержится много нового: традиционно точность оценок считалась равной полпорядку величины (от 0.32 = Ю-0'5 до 3.2 = Ю0 5). Однако важна принципиальная возможность вычисления этой величины, она придает оценкам по порядку величины уверенность, в ряде случаев достаточную для конкретных приложений.

2 Гидродинамика с переменными вязкостью и теплопроводностью

В этом разделе поведение горных пород на больших временах при источниках тепла рассматриваются как перемещение очень вязкой жидкости с произвольными или же эмпирическими зави-

снмостями от температуры вязкости

»](Т) = щ ехр(в/Т); 6 = 2 • 104 -Н 7 • 104 К. (8)

и теплопроводности

к(Т) = Коо(^г)"; " = 0.23Н-0.9; = 2.8-5-3.2 Дж/(м с К). (9) 2.1 Плоское сдвиговое течение

Будем решать совместно уравнение Навье—Стокса и уравнение непрерывности [15]

д1 Ук дхк дх{ дхк^\дхк дх{) дх{ \дхк ]

(10) с»)

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

Замыкает систему уравнение теплопроводности с тепловыделением, вызываемым вязким трением, а также внутренними источниками с объемной плотностью ю [Вт/м3]:

Мы получим точное решение системы (10-12) в простейшем случае, для одномерного течения. При отом, однако, вязкость и теплопроводность будут рассматриваться как произвольные функции температуры.

Пусть течение направлено по горизонтальной оси х. Тогда скорость есть ьх = V, ьу = 0, ьг = 0. Все функции считаются зависящими только от г, ось г направлена вниз. При этом уравнение непрерывности удовлетворяется автоматически, г-компонента уравнения (10) дает обычную зависимость давления от глубины ¿р/йг = др{г). После интегрирования ж-компоненты уравнения (10) градиент скорости выражается через вязкость:

•■¿Иг) - Е-ет <13>

Здесь постоянная интегрирования Е с размерностью кг/(м с2) имеет физический смысл сдвигового напряжения, сохраняющегося поперек текущего слоя Е = стхг = <тгх — const. С одной стороны очевидно, что напряжение S не должно превышать максимальный предел прочности сгт холодной породы. С другой стороны, наиболее интересные явления происходят именно при напряжениях порядка £ ~ сгт: при появлении единичной трещины, напряжение сбрасывается, пространство около трещины разогревается, вязкость падает и порода начинает течь, оаплавляя трещину.

Далее, из уравнения теплопроводности (12) получаем

При w — 0 ото уравнение не содержит явно z и разрешается в квадратурах. Перейдем к переменной Т, введя по определению функцию N(T) = dT/dz. Вместе с решением уравнения (13) для скорости, точное общее решение системы (10-12) уравнений гидродинамики — теплопроводности для стационарного плоского течения в поле тяжести имеет вид:

N {1) кКГ) J п(Т) ' р(г)=/ gjm-dT• (15)

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

2.2 Распределение температур

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

толще плиты. При этом граничное условие на поверхности фиксирует наблюдаемый температурный градиент ^(Го) = Nо при г = 0. Тогда решение примет вид:

= (И)

То

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

Решение (16), формально продолженное до конечной точки — температуры Т — Тт, дает связь этой новой постоянной с поверхностным градиентом:

2 Тт

/V2 - — / <ПлТ' т1)

~ к» У „сг'Г • (17)

Т0

Выражая в (16) постоянные Е, «о через Щ, получим выражение для обезразмеренного градиента температуры

2 / V N2(Г) п (т) -

2ггл „2 г _ А(т) ■

Л{тт)

Т

г = (18)

-¿о

к2СО

Здесь введена безразмерная функция А(т), по порядку величины близкая к единице. В общем случае она зависит от температурных зависимостей вязкости и теплопроводности, а при явном использовании приближенных выражений (8,9) зависит от параметров и и О = Тов. Она равна

р.)

В свою очередь, зависимость температуры от глубины получается обращением второго из уравнений (15), выписанного с учетом граничного условия Г|г=о = То

Отсюда следует, что характерная глубина, на которой существенно изменяется температура, порядка Но = То/Ио — 104 м = 10 км.

Для выяснения физических причин полученной глубинной зависимости температуры, рассмотрим, как изменяются с глубиной скорость и тепловой поток.

2.3 Тепловой поток и скорость

Для рассмотрения зависимостей от глубины г или (от температуры Т) скорости V и теплового потока q = квТ/йг = кИ удобнее от квадратур (15) вернуться к дифференциальным уравнениям (13,14). Их зависимости оказываются подобны, — обе обратно пропорциональны вязкости. Исключая вязкость, получим точно интегрируемую связь скорости и теплового потока

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

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

т

(20)

¿г

и имеет тот же порядок величины, что и относительные скорости плит, то есть и = 1-^10 см/год~ Ю-9 м/с.

Сдвиговое напряжение Е в решении (21) также поддается независимой оценке. Очевидно, что £ не должно превышать предела прочности горных пород (тт. В противном случае произойдет разрыв плиты. С одной стороны, разрывы плит реально происходят (Красное море, Байкал), с другой стороны, это случается относительно редко. Это и дает возможность считать, что Е ~ сгт.

Предел прочности даже изотропной и однородной горной породы не постоянен по глубине. Однако с увеличением глубины величина сгт подвержена противоположными тенденциям: она увеличивается с давлением, но падает с ростом температуры. Зависимость от давления грубо пропорциональная, а зависимость от температуры — более сильная, экспоненциальная, как ехр(0/Т), с той же эффективной энергией связи 0, что и вязкость. Поэтому своих максимальных значений предел прочности достигает, по-видимому, на тех же, порядка Но — То/No — 10 — 15 км глубинах, что и максимум температурного градиента. На этих глубинах прочность скальных пород по порядку величины равна ат = 100 МПа. Поэтому для характерного значения напряжения Е примем оценку на полпорядка меньше

Е ~ crm/3 ~ 3 • 107 Па.

Сравним теперь оценки уравнения (21):

v£ ~ u£ ~ Ю-9 х 3 х 107 = 3 х Ю-2 ~ g0 ~ Ю-1. (22)

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

3 Гидродинамика погружения нагре-

того тела

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

Классическая формула Стокса (^ = 6ж7]Яи) описывает силу сопротивления твердой сферы радиуса Я, движущейся со скоростью и в жидкости с постоянной вязкостью т]. Эта сила определяется гидродинамическим решением на расстояниях порядка радиуса шара. В разложении силы сопротивления по степеням числа Рейнольдса 11е большие расстояния дают поправку первого порядка, а малые — определяют решение в главном, нулевом приближении и дают поправку второго порядка. Если тепловая мощность шара достаточна для уменьшения вязкости вблизи его поверхности, то гидродинамика в этой области будет существенно отличаться от решения Стокса, что должно привести к модификации формулы сопротивления [1,4,5].

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

3.1 Тепловое поле

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

стационарное уравнение теплопроводности

на расстояниях, больших по сравнению с радиусом источника Л. На таких расстояниях скорость ьь и коэффициент теплопроводности к близки к своим значениям на бесконечности к = коо, Ыс = ик- В этом случае решение уравнения (23) в сферических координатах имеет вид:

Т(г,0)=Г«

. (1 — со вО)иг 1 + —-ехр(—--

г

" "(24)

Здесь Р есть тепловая мощность источника и х = коо/рСр-

Из решения (24) ясно, что на расстояниях, много мень-. ших тепловой длины (для г «С х/и) тепловое поле сферически-симметрично. Это также ясно из уравнения (23): для источников с центральной симметрией решение не будет зависеть от углов в случае малых чисел Пекле Ре= ги/х -С 1. Наличие сферической симметрии теплового поля будет необходимым условием последующего гидродинамического решения. Противоположный случай больших чисел Пекле, но постоянной вязкости расплава был рассмотрели Б-Етегшап и Б.ТигсоЦе [17].

Для горных пород число Прандтля велико: Рг = ц! РХ 1-Поэтому условие сферической симметрии теплового поля автоматически приводит к тому, что число Рейнольдса мало:

„ иЯр ^ и/£ „ /<ч_.

Ее = —- < Ре = — <1. (25)

П X

Известно [18], что в разложении сопротивления по степеням числа Рейнольдса малые расстояния дают вклад нулевого (и второго) порядка, а большие — первого. Это оправдывает использование сферического приближения для теплового поля. Однако напомним, что величина скорости на бесконечности и в (25) нам пока неизвестна. Она будет найдена после решения уравнений гидродинамики из равенства архимедовой силы силе сопротивления. По этой причине диапазон применимости решения в исходных физических параметрах задачи будет определен после получения решения.

Параметр /?1 в (24) пропорционален тепловой мощности Р. Для горячего шара отношение /¿1/Л будем считать, вообще говоря, порядка единицы, не исключая возможности Я. В этих случаях температурная зависимость коэффициента теплопроводности к становится существенной. Для реальных горных

пород в значительном температурном интервале она хорошо аппроксимируется степенным законом (9). При этом оказывается возможным получить аналитическое решение уравнения (23):

Т„ф\{г) = Гс<

Tv=i(r) =Тооехр(—).

г

(26)

В области Д < г < Л] < и/х решения (24) и (26) совпадают: Т—Тоо ~ г-1. Как мы увидим далее, для задачи вычисления силы сопротивления горячей сферы существенную роль играет логарифмическая производная температуры L = d InT/d lnr. Эта функция имеет максимум на расстоянии г = R\. Отметим, что решение (26) оказывается неверным, если одновременно v > 1 and R\ > Rf{y — 1). Это означает, что стационарное решение (23) несправедливо, если тепловая мощность слишком велика. Однако этот случай "теплового взрыва" не представляет большого интереса для гидродинамических потоков в геофизике.

3.2 Течение вокруг тепловыделяющего шара

При малых числах Рейнольдса уравнение Навье—Стокса линейно по скорости и,-. Входящая в него вязкость г] — функция температуры, — после решения уравнения теплопроводности (23) она становится заданной функцией координат. Кроме того, в уравнениях гидродинамики [15] необходимо оставить конвективный член, пропорциональный коэффициенту объемного расширения /?, и член, пропорциональный ускорению свободного падения <7,-:

яг)-«»™ <эт>

В системе отсчета, где шар неподвижен, скорость на его поверхности г = R равна нулю, на больших расстояниях (на бесконечности) она равна постоянной скорости обтекания «¿.Скорость щ пока неизвестна по величине (ее определение и есть наша цель), но известна по направлению — она параллельна уско-

рению свободного падения <7,-. Ищем решение (27) в виде ( f \

= «< + J ф(0 df - \гч ф(0; С = Чг/Ri), V (о )

(28)

которое автоматически удовлетворяет уравнению непрерывности dvk/dxk = 0 и требует выполнения граничных условий

оо

Ф(£о) = 0; Ф(оо) = 0; J Ф <Ц = 1-, = 1п(ВД). (29)

fo

Исключая из уравнения (27) градиент давления р операцией rot, получим неоднородное линейное уравнение третьего порядка

<Э3Ф гх32Ф ,т2 „г г дЬ^Ф

f

+(L2 - 5L - б + Щ)Ф = —2Grexp(— J(L - 1) d£). (30)

fo

Уравнение, совпадающее с этим, но без конвективного члена правой части, было получено J.Mathews [19], который решал задачу сопротивления для иона, движущегося в жидком гелии, где L(£) < 0. В уравнение (30) входит одна постоянная — число Грассхофа Gr и одна известная функция £(£). Число Грассхофа определено здесь, как

^_____ PpPR _ L<yn gpRiR

Функция L(£) есть логарифмическая производная вязкости:

a Inr а£

Известная зависимость температуры Т(г) (26) позволяют установить конкретный вид функции £(£), зависящий от двух

параметров 0 и V, приведем его для двух знамений показателя V коэффициента теплопроводности:

^ = (33)

Заметим, что выбор нуля £ в уравнениях (29, 30) произволен. Он выбран при г = для упрощения формул (33), при таком выборе максимум функции £(£) находится при £ = 0.

Давление р в жидкости выражается в общем виде через решение Ф(£) уравнения (30):

, ЩкПк

д2Ф /п т,дФ тж —+ (3 + Ь)- + 1Ф

(34)

3.3 Решение с модельной вязкостью

Уравнение (30) с граничными условиями (29) не удалось решить точно. Однако, наша цель — получить формулу для силы сопротивления. Она, как отмечалось, определяется течением в прилегающем к шару слое. Вдоль него логарифмическую производную вязкости можно считать постоянной, равной своему значению на поверхности шара Х(£) ~ ¿(£о) = Ьо- Иными словами, мы ищем точное решение уравнения (30) для модельной зависимости вязкости т](г) = 7](Л)(г/И)ь°. Такое модельное решение для Ф(£) даст возможность получить выражение для силы сопротивления, справедливое для больших Ьо " произвольных Ь(£) и одновременно даст правильный, стоксов предел при Ьо —* 0, хотя и не будет справедливо при произвольных £(£), если Ьо ~ 1. В случае Ь(£) = Ьо решение уравнения (30) легко получить:

з

Ф = £С,ехР(М£-&,))- (35)

¿ = 0

Здесь индексом 0 обозначено частное решение неоднородного уравнения:

к0 = -Ь0 + 1\ Со = Сг(£о + 4)-1. (36)

Остальные решения есть решения однородного (с правой стороной уравнения (30) равной нулю), они соответствуют трем

(в = 1,2,3) корням кубического уравнения

к3 + 2(Ь0 + 1)к2 + {Ь% + ЗЬ0 - Ь)к + {Ь20 - бЬо - 6) = 0. (37)

Одно из них, определяемое наибольшим корнем, к 1, возрастает с радиусом (или же, при Ь > 6, убывает недостаточно медленно), поэтому коэффициент при нем должен быть равен нулю С\ — 0. Два отрицательных корня позволяют получить решение, удовлетворяющее граничным условиям (29):

¿0 + 1 у/Ц- 2Ь0 + 25 2 ~ 2 2

<39>

Из (36) ясно, что для выполнения граничных условий должно выполняться одно из условий: или ¿0 > 1, или вг = 0.

Функция Ф(£) убывает экспоненциально. Это означает, что скорость V,- приближается к скорости на бесконечности щ как степенная функция радиуса. Если Ьд 1, то &о,2,з ~1" — ¿о и эффективная толщина слоя обтекания равна 8Л ~ И/Ьо- Это и оправдывает использование модельного решения (35) для вычислений силы сопротивления в физических ситуациях с реальной зависимостью £(£). Строго говоря, все корни в (37) должны быть разложены по степеням 1 <С 1. Эта операция должна выполняться с должной осторожностью, поскольку следует удерживать различное число членов при вычислении разных величин, Некоторые трудности при этом связаны с наличием ложной особенности при Ьо = 6, когда два корня совпадают =

При Ь ~ 1 модельное решение (35) имеет качественный характер, оно нуждается в поправках, зависящих от и и с1Ь/с1£(£о). Однако решение (35) имеет правильный предельный переход при Ь = 0 в решение Стокса (к2 = —3, к3 = —1). Заметим, что для этого перехода сначала положить равным нулю число Грассхофа, и только затем Ьо —» 0.

Заметим, что число Грассхофа (31) не мало при и —* 0. Известно, что решение Стокса нарушается на больших расстояниях и при больших скоростях. Однако кроме того, при любом

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

1.-3 1_»б

Рис. 3. Линии тока вокруг горячего шара при различных значениях Ьо, логарифмической производной вязкости на его поверхности. Случай Ьо = 0 соответствует решению Стокса.

Линии тока вокруг горячего шара иллюстрируются на рис. 3 при числе Грассхофа вг= 0 и трех значениях Ьо. Заметьте, что с ростом Ьо линии тока примыкают ближе к поверхности шара и область возмущений сокращается.

3.4 Сила сопротивления

Сила сопротивления сферы Fi вычисляется стандартным образом [15], как интеграл по поверхности шара от проекции на вертикаль суммы нормального давления и вязких напряжений:

Ди,- = щпк - ^ + ^ . (40)

Используя решение (35) и принимая во внимание (28, 34), по-

лучим:

2л-

^ = у

+ к\)

Напомним, что £2 дано в формуле (38), где Ьо следует вычислить по формулам (33) на поверхности шара, то есть при £ = £о — 1п(47гЛ/с00Т00/Р). Число Ьо пропорционально отношению ©/Г,» и в общем случае не мало.

Чтобы найти скорость движения шара в случае, когда его средняя плотность отличается от плотности жидкости на Ар, приравниваем силу сопротивления силе Архимеда А'кдАрЕ?/Ъ, вызываемой разностью плотностей шара и среды. Чтобы получить выражение, справедливое в физических приложениях, удерживаем только высшие члены разложения по Ьо. Помня, что число Грассхофа, (31) также содержит скорость и, получим выражение для Скорости погружения (при Ар > 0) или всплытия (при Ар < 0) тепловыделяющего шара:

и =

2дАрВ? т/(Я )Ц

1 - Ьо/ЗТа

АрК_

(42)

3.5 Случай чистой конвекции

Как следует из (42), возможен случай, когда скорость на бесконечности равна нулю. Это происходит, если

1 = Ьо/?Тсод^-^(То-Тоо)д^- (43)

Такое условие физически понятно: малая разность плотностей между жидкостью и шаром компенсируется убыванием плотности жидкости при ее тепловом расширении в слое ¿Я ~ Я/Ьо вокруг сферы (То — поверхностная температура).

Вообще говоря, в таком случае мы должны искать решение в форме, несколько отличающемся от (28, 29 ): и; должно быть пропорционально <7,- (и дкПкП;) а интеграл в (29) должен быть равен 0, а не 1. Однако можно просто взять предел и —► 0 в решении (35), приняв во внимание, что Со ~ в г ~ и-1, и ^ ~ щ.

Таким образом получается поле скоростей в случае чистой конвекции:

2А рВ?

ШЩ

9*

( * ^ + /

{о

2Зъпк

/

= ехр(А0« - 6)) + кЛеМЫ£ ~ Ы)

+

к0(к3 - к2) к3(к2 — к0)

к0(к3 - к2)

ехр(А3«-^о)).(44)

Здесь размерный коэффициент при скорости дан в физическом пределе Ьо 1, но функция Ф^) выписана в модельном представлении. Такое выражение более удобно, поскольку эта функция содержит ложную особенность при Ьо = 6:

Ф!(Х0 = 6) = ехр(—5(£ - £„)) - (1 + 14£/5)ехр(-7(£ - 6))- (45)

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

V

3.6 Условие применимости

Вернемся к критерию (23), подставив в него ответ для скорости обтекания (42), в предположении, что конвективный член мал:

Я <1,о

1/3 -Ьо №)1

. 9&Р . . »716 .

1/3

х 100 ш.

(46)

В последней численной оценке

X = Ю-6 м /с, Ар ~ 1 х 103 кг/м", Щб = Ю16 Па с. (47)

Таким образом, (41) может применяться в геофизике в весьма широком диапазоне параметров. В частности, оно используется для численных расчетов глубин погружения захоронения ядерных отходов.

4 Статистика дискретных смещений

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

Наблюдаемая энергия землетрясений Е измеряется, как правило, в логарифмической шкале магнитудой М, определяемой как

Энергия, выделяющаяся при разломе, в свою очередь выражается через его размеры и характеристики среды. Разлом происходит в точке, где один из компонентов недиагональных локальных напряжении а^к — <тц6ц превосходит прочность породы ат. Трещины имеют как правило структуру, близкую к плоской: их площадь 5 порядка квадрата линейных размеров Ь, а смещение вдоль разлома А1 всегда много меньше его длины Ь. По порядку величины энергия землетрясения (разлома) есть

В силу закона Гука длина разлома Ь и величина смещения Д/ связаны между собой приближенным соотношением Ы/Ь ~ <тт/К, где К — модуль всестороннего сжатия, который для большинства твердых тел имеет тот же порядок величины, что и модуль Юнга. Соотношения между линейными размерами землетрясений и энергией подтверждаются экспериментальной статистикой:

М = |(18Я[ Дж] — 4.8).

(48)

Е ~ <гтЬ2А1.

(49)

Для разнообразных логарнфмически-шнроких распределений, возникающих при делении целого на большое число разных по размерам частей установлен эмипирнческнй закон распределения (Б.А.Цзубников [21]), названный им законом распределения конкурентов: аддитивная величина ранга г (номера в убывающей последовательности) пропорциональна г-1. Нормируя, например, распределение N частей полной массы М, имеем приближенное выражение для г-той массы:

т" = г(С-НпАГ)- (51)

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

= = (52)

г = 1 Г атГ Г = 1

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

/д, N2 УуМ

(Л/г) = -Кг(С + ЬМУ (53)

Наконец, используя связь (50) между смещением и энергией, мы приходим к результату, совпадающему с законом Гутенберга—Рихтера для распределения землетрясений по маг-нитудам:

ЕГ^—(А1Г) - у гз(с + 1п^)з (54)

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

5 Практические приложения

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

Предложенная автором концепция захоронения возникла как продолжение давно известной [22], но еще продуктивной [23] идеи самозахоронения ядерных отходов. В ее развитие была предложена [7,9] идея захоронения высокоактивных отходов в глубоких скважинах, заполненных каким-либо водонерастворимым легкоплавким материалом. Новая концепция динамического захоронения ядерных отходов также открывает принципиальные возможности [10] для изучения глубинных слоев земной коры, недостижимых традиционным бурением. Затем было показано, что во многом оптимальным заполнителем скважин для захоронения может служить сера [8, 11, 12]. Использование серы в качестве заполнителя скважины дает ряд преимуществ по сравнению с традиционной схемой самозахоронения ядерных отходов [22]. В то же время существенно большая глубина и связанная с этим дополнительная безопасность захоронения дает принципиальные преимущества предлагаемому методу по сравнению с известными методами захоронения ядерных отходов в горных породах [24] при соизмеримой с ними емкости.

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

1. Диаметр скважины в скальных породах (гранитах) может не превышать 0.3---0.4 м в ее основном сечении и

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

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

3. Загрузка ядерных отходов начинается после заполнения скважины серой. Через интервалы времени порядка часа герметичные капсулы с отходами последовательно опускаются через устье скважины в серу и погружаются под действием собственного тепловыделения. Максимально возможный диаметр капсул должен быть примерно вдвое меньше диаметра скважины. Низкая температура плавления природной серы (113 — 120°С) и ее низкая теплопроводность (на порядок ниже теплопроводности гранита) обеспечивают возможность проплавления капсулами малых размеров. При этом температура на стенках скважины в течение всего процесса загрузки остается ниже температуры плавления серы. В случае, если временной интервал между последовательными капсулами не превышает суток, последовательное погружение капсул происходит по тепловому каналу, прогретому предыдущими капсулами, только время застывания серы в канале составляет несколько часов. Тепловыделение ядерных отходов с возрастом порядка 10 лет достаточно для последовательного проплавления серного канала при диаметре капсул около 0.2 м.

4. Приблизительно через год после начала загрузки отходов первые капсулы достигнут дна скважины и начнут там накапливаться. Это приведет к росту температуры, которая примерно через год возрастет до 300 -=-500°С. Очевидно, что вследствие замедления ядерного тепловыделения температура затем будет убывать. В условиях неразрушающейся скважины нетрудно рассчитать значение максимальной температуры, которая достигается через 5-10 лет после начала аккумуляции. Эта величина растет с увеличением диаметра скважины и падает с ростом начального возраста отходов.

5. Вследствие теплового расширения горных пород диаметр скважины в ее активной (нижней) части сперва начинает увеличиваться и в горных породах возникают дополнительные сдвиговые напряжения. При температурах выше примерно 300°С эти сдвиговые напряжения превысят предел прочности гранита. Более того, при этих же температурах начнутся химические реакции между серой и оксидами железа, содержащимися в гранитах, они также понизят прочность горных пород. Оба процесса приведут к тому, что начнется разрушение стенок скважины вблизи ее дна.

6. Тепловое расширение и эрозия стенок скважины приводят к потере стабильности и обрушению нижней части скважины. В этом случае вся накопленная тяжелая колонна отходов (несколько сот метров) начнет деформироваться и опускаться ко дну скважины, формируя ансамбль капсул, перемешанный с обломками горных пород, уже прореагировавшими со всей доступной серой. Основными продуктами этих реакций являются водонерастворимые пирит ГеБг и пирротит. Температура этого ансамбля будет постепенно повышаться вследствие того, что поверхность теплоотвода каплеобразной (близкой к сферической) структуры меньше по сравнению с начальной цилиндрической структурой. Поэтому возможно достижение условии, когда температура ансамбля превысит температуру размягчения гранитов (около 1600'С).

7. Поскольку плотности и пирита (5.5 г/см3), и ядерных отходов выше плотности гранита (2.7 г/см3), горячий ансамбль капсул вместе с расплавленным пиритом отделяется от основной скважины и начинает проплавлять горные породы в режиме са-мопогруження. Основная масса отходов с уменьшением тепловыделения и остыванием ансамбля постепенно замедляет свое погружение. После застывания водонерастворимого пирита при температуре около 1100°С он формирует вторичную матрицу, предохраняющую захоронение от дальнейшей утечки радионуклидов в окружающую среду.

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

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

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

5.1 Потоковый режим загруоки

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

Определим потоковый режим загрузки как такой, при котором расплав серы в стволе скважины не успевает существенно остыть за интервал времени между последовательными капсулами. Вследствие малой теплопроводности серы характерное время охлаждения серной колонны оказывается довольно большим. Так, для скважины радиуса Л = 0.2 м время охлаждения серы по порядку величины равно

т-5 ~ Я2/\'<? ~2х 105 с, (55)

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

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

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

Другим важным параметром является объемное тепловыделение отходов w, которое в оценках предполагается характерным для остеклованных ядерных отходов с временем после извлечения из реактора в диапазоне от 5 до 30 лет и стандартными концентрациями радиоактивных изотопов. Ее порядок величины составляет w ~ 104 Вт/м3.

После фиксации размера и тепловой мощности капсул остается еще один параметр потокового режима загрузки: временной интервал Дt между опусканием последовательных капсул в устье скважины. Условие тепловой "памяти" при потоковом режиме накладывает ограничение на величину этого интервала: At <С Ts ~ 105 с. Это означает, что характерный промежуток времени между капсулами должен не превышать нескольких часов.

Скорость погружения капсул как функция глубины z зависит от начальной температуры T(:0(z) горных пород на данной глубине (она же —■ температура на большом горизонтальном удалении от скважины). Как известно, температура недр приблизительно линейно возрастает с глубиной. Для пород кристаллического щита вдали от сейсмически активных областей геотермический градиент относительно невелик: dToo/dz = 20 30 К/км.

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

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

Решение тепловой задачи потоковой загрузки начинается с определения теплового поля в окружающих скважину горных породах. Флуктуации температуры в них, вызванные проходом отдельных капсул, относительно невелики, а тепловой поток от скважины в окружающие породы флуктуирует слабо. Введем величину У/ как удельную (на метр вертикального ствола) среднюю по быстрому времени тепловую мощность, выделяемую потоком капсул. Удельная мощность является функцией глубины 2 по той причине, что от нее зависит окружающая температура Тоо(г), а следовательно, и скорость погружения. Эта погонная тепловая мощность переносится через серу к стенкам скважины и затем отводится теплопроводностью горных пород. Пока неизвестная скорость погружения и, тепловая мощность отдельной капсулы Р и удельная мощность Ш связаны между собой очевидным соотношением У/ — Р/иМ.

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

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

Т(г > Е,

(Я < г < А), (56)

А2 = 4ехр(-С)хо(< - <0(г), С = 0.5772.

(57)

радиусом скважины. Поскольку в формулу (56) она входит под логарифмом, требования к точности ее вычисления не очень высоки. Заметим, что и ¡зависимость температуры в среде оказалась логарифмически медленной. Наконец, радиальный температурный градиент ¿Т/йг вблизи стенок скважины много больше, чем вертикальный, геотермический градиент йТ/йг. Эти результаты оправдывают предположение квазистационарности.

Момент ¿о(г) прохода данной глубины первой капсулой следующим образом выражается через пока неизвестную скорость погружения «о(^):

Поскольку он также входит в выражение (57) под логарифмом, к точности его вычисления не следует предъявлять завышенных требований. По этой причине для его вычисления мы можем использовать не истинную скорость первой капсулы ио(г), а скорость и(г,< —» 0), которая вскоре будет определена из нашего квазистационарного приближения.

Пз общей формулы (57) для температуры в среде в частности можно получить выражение для температуры Тц на стенках скважины, при г = Л: Эта зависимость от времени I и от глубины г послужит граничным условием для внутренней тепловой задачи. (Напомним, что пока граничное условие выражено через пока неизвестную удельную тепловую мощность \У(1,г)).

Чтобы найти ее и затем скорость погружения и, необходимо найти распределение температуры в серной колонне с заданным тепловыделением и граничным условием Т(Д) = Тц. Для описания теплового поля внутри серы можно использовать полученное ранее тепловое поле точечного погружающегося источника (24). Это решение справедливо для движения с большими числами Пекле Ре= игс/х 1 па расстояниях, больших радиуса капсулы г > гс. Перепишем его в цилиндрических координатах с началом отсчета в центре тепловыделяющей сферы для случая, когда расстояние между последовательными капсулами много больше радиуса скважины, т.е. иА1 Я > гс. При этом выражение для теплового поля предыдущей капсулы к моменту подхода последующей (<х ~ Д<) может быть упрощено разложением этой

2

(58)

о

формулы по малому параметру г/и1\:

= (59)

Абсолютная величина показателя в этом выражении при г ~ Я и <1 ~ <5£ оказывается велика по сравнению с единицей, что и оправдывает приближенную замену точного граничного условия на Т(Я,*1 ) = ТЛ.

Для нахождения суммарного теплового поля всех предшествовавших капсул необходимо просуммировать прирост температуры от N тех капсул, чье выделившееся тепло еще не перешло в окружающие горные породы, а нагревает серу. Его можно оценить из приближенного условия, что тепловая длина в сере в течение времени Л^Д^ равна радиусу скважины. Число N ~ 4/22/х5 оказывается достаточно большим, что позволяет заменить суммирование интегрированием и взять интеграл, используя равенство Эйлера. В результате получим внутреннее тепловое поле в квазистационарном приближении, которое на стенках скважины оно удовлетворяет решению внешней задачи:

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

и(1,г) =

4тг(Тт - Тоо)Д<

1 4ехР(-С)хо(< - и) , 2 Я

—1п-—--—1п—

Ко № «5 г с

(61)

Это выражение имеет несколько ограничений для пределов применимости, которые упоминались по мере вывода, но содержали в себе скорость погружения. Выпишем их теперь в явном виде. Во-первых, предполагалось, что число Пекле велико, а число Рейнольдса мало. Следовательно, игс хб и игс < ^б/йб-Совместимость этих условий возможна, поскольку у серы достаточно велико число Прандтля Рг(123°С) ~ 30. Наконец, вместе с условием потокового приближения (иД2 Я > гс) получим

совокупность ограничении:

< ^ < Д< « < —. (62)

ЧБ и ХБ Хэ

Численные оценки, введенные в интерактивную вычислительную программу 'Ь0АВ"\УА5Т' (дискета-приложение к монографии [1]), покалывают, что для скважины диаметром 2/2 = 0.4 м, заполняемой капсулами диаметра 2гс = 0.2 м выполнимость этих условий достигается при интервале времени 2 ч- 4 ч. между последовательными капсулами. Скорость погружения возрастает от 2 х Ю-5 м/с до 1.2 х Ю-4 м/с или 1.7-=-10 м/сут. Время достижения глубины 4 км вариируется в зависимости от мощности и диаметра скважины от 1 до 5 лет.

5.2 Термическое разрушение скважины

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

Изначально породы на глубине г считаются находящимися под давлением верхних слоев р^ ~ pgz, вызывающим однородные напряжения = —poo^ik- В процессе бурения скважины возникают дополнительные напряжения, вызванные тем, что заполнитель скважины (буровая жидкость) имеет плотность, меньшую чем горные породы. Будем считать давление внутри скважины р произвольным — как большим, так и меньшим Pca(z).

Начальная температура вокруг скважины равна Too- После заполнения (в момент f = 0) скважины ядерными отходами с известным тепловыделением, температура в окружающих горных породах начинает повышаться. Тепловое поле T(r > R) вокруг скважины радиуса R нетрудно вычислить, но пока будем считать его произвольной функцией радиуса.

Теория упругости [25] дает связь между вектором смещения щ, тензором деформации и тензором напряжений <7,1. В частности, смещения, вызванные термическими напряжениями, опре-

деляются уравнением:

и„ = divu = j + l/ J AT, AT = T(r) - Too • (63) J^l — v)

где /3 — коэффициент теплового расширения, ai/ — коэффициент Пуассона.

Обычные граничные условия для уравнений теории упругости выбираются таким образом, чтобы начальное состояние с нулевым давлением было недогруженным. Однако для земных недр на заданной глубине z удобнее считать нулевым смещение при начальной температуре Too(z) и начальном давлении Pco(z).

При цилиндрической симметрии тензоры деформации и на-пряжний имеют только две ненулевых компоненты (гг) и (фф), а вектор смещения направлен по радиусу и,- = u(r)r,-/r. Интегрирование уравнения (63) дает

г

= /АТ(Г')ГМГ' + -,

3(1 - v) г J г

(64)

а постоянная интегрирования В, в свою очередь определяется из решения уравнений Гука теории упругости с граничными условиями <гРГ(Я) = —Р + Роо- Таким образом возможно полностью решить задачу при произвольном тепловом поле цилиндрической симметрии:

(1 + |/)(Р-роо) (1 + и) /?ДТ(Я)Я2

3(1 - 2i/) г

(65)

«ТГ = —(р — Роо)—5" +

R2 Е(3

ДТ(г) - ДГ(Л)

0L

т-2

3(1-2//) .

Г

-W^-)HAnr'ydr'' (6б)

.Л2 EpAT(R)R2 EpAT(r)v (Тфф - (Р-Роо)~2- + -+

г2 3(1-21/) г2 3(1 - 2i/)(l - и)

г

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

(2оуг - (Тфф)/3, (2сгфф - <тгг)/3, (<7ГГ + (?фф)/3. (68)

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

г) ,дч_2(р-Рсо) , 2Д/?ЛГ(Д)

D°h{R)~-з-+ 9(1-2W(1-W ( }

Сравнение этой величины с прочностью горных пород — предельным напряжением сдвига сг„, позволяет сформулировать критерий термической устойчивости скважины. Небольшой (в обычных условиях) разностью давлений р — Роо, можно пренебречь по сравнению с ат. Кроме того, конкретизируем тепловое поле T(r,t) вокруг скважины. Для расстояний порядка радиуса скважины и моментов времени, когда тепловая длина Ат = у/хО достаточно велика (Ат Я) зависимость температуры от радиуса дается выражением:

wm ы

4тгк r¿

Здесь W(<) — удельная (на единицу длины) тепловая мощность отходов. Подставляя его в выражение (69) и сравнивая с предельным напряжением сдвига, получим критерий стабильности скважины, выраженный в двух вариантах — через температуру и через тепловую мощность:

E/3AT(R)

T(t,r,z) = T00(z) + ^-(\n^-C). (70)

3(1 - 2i/)(l — i/)

<<rm, (71)

WC) m

12™(l-2i/)(l-i/)v R2

Модуль Юнга горных пород имеет порядок величины Е ~ 6 х Ю10 Па, а коэффициент теплового расширения порядка /3 ~ 2х10-5 1/К. При и ~ 0.2 это дает предельный рост температуры равный АТ ~ 250° С. При рассматриваемых радиусах скважины и удельных тепловыделениях порядка 1 кВт/м требуется примерно год для такого разогрева. Более подробные вычисления при произвольных параметрах дает программа 'ЬОАБ\УА8Т' — дискета-приложение к монографии [1].

5.3 Погружение захоронения ядерных отходов

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

В режиме численного интерактивного моделирования при свободной вариации исходных параметров дана оценка глубин, которых асимптотически достигнет захоронение на больших временах. Оказывается, что оно экспоненциально зависит от исходной полной тепловой мощности. Так, при начальном среднем энерговыделении 1.5 • 103 Вт/м3 захоронение диаметра 4.2 м полной массы 150 т через 30 лет (период полураспада двух наиболее активных изотопов) достигает дополнинтелыюго углубления в 17 километров.

Выводы

1. Получена статистика коэффициентов физических формул по корпусу текстов "Курса теоретической физики" Ландау— Лифшица. Анализ статистики позволил установить ряд закономерностей для частотностей коэффициентов, в частности, возможность описания последовательного образования коэффициентов с единственным свободным параметром. Проведено численное моделирование генезиса коэффициентов, сравнение модельного и наблюдаемого распределений позволяет вычислить параметр генезиса. Распределение коэффициентов за исключением первых, наиболее вероятных членов, оказалось близким к логнормалыюму, что позволяет вычислить его логарифмическую дисперсию и тем самым установить точность оценочных формул (формул, не содержащих численных коэфициентов).

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

3. Найдено температурное поле в среде с движущимся точечным тепловым источником, задающее вязкость среды как функцию координат.

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

5. Показано, что распределение дискретных перемещений (разломов) твердых тел по их величине при фиксированном в каждый момент суммарном перемещении является частным случаем конкурентного распределения. На этой основе получена связь распределения Гутенберга—Рихтера для частотности землетрясений с законом конкурентного распределения.

6. Дано решение тепловой задачи потоковой загрузки последовательности тепловыделяющих капсул.

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

8. Численными методами рассмотрена динамика самопогружения захоронения ядерных отходов.

ЛИТЕРАТУРА

1. Byalko A.V. Nuclear waste disposals. Geophysical safety, CRC Press, 1995, 277 p.

2. Бялко А.В. Статистика коэффициентов в физических формулах// Теоретическая и математическая физика, 1991, т. 88, N1, 153-160.

3. Бялко А.В. Распределение коэффициентов// Природа,

1995, т\,$4Ж

4. Бялко А.В. Модификация формулы Стокса (ламинарное движение тепловыделяющего шара)// Письма в ЖЭТФ, 1992, т. 55, N3, 190-194.

5. Byalko A.V. Dynamic safety of IILW disposals// Proceedings of European safety and reliability conference '92, Copenhagen, Elsevier, 1992.

6. Бялко А.В. Конструктивность закона конкуренции// Природа, 1993, N11, 15-19.

7. Byalko A.V., Khavroshkin О.В. A method of disposal of high level nuclear waste by deep sinking// The Third International Conference on Nuclear Fuel Reprocessing and Waste Management, RECOD'91-Proceedings, Sendai, 1991, V.I, 486-490.

8. Byalko A.V. NEWDEEDS Concept for HLW Disposal// The Fourth International Conference on Nuclear Fuel Reprocessing and Waste Management, RECOD'94-Proceedings, London, 1994, V.II (no page numbering).

9. Бялко А.В., Хаврошкин О.Б., Халатников И.М. Способ захоронения высокоактивных нуклидов, патент SU N1725667, приоритет от 22.01.1990.

10. Бялко А.В., Хаврошкин О.Б. Способ получения информации о глубинных структурах литосферы Земли, патент SU N1787279, приоритет от 6.7.1990.

11. БялкоА.В. Способ глубинного захоронения высокоактивных ядерных отходов, заявка N4909190/25 от 20.1.1991, положительное решение от 3.6.91.

12. Бялко А.В. Состав земли не знает грязи (Захоронения ядерных отходов, новые идеи) // Наука и жизнь, 1993, N2, 60-64.

13. Бялко А.В. Лондон, апрель: конференция по ядерному топливу // Природа, 1994, N7, 13-14; Ложная посылка // Наука и жизнь, 1994, N8, 106-107.

14. Newcomb S. Note on the frequency of use of the different digits in natural numbers// Amer. J. Math. 1881, 4, 39; Bedford F. Proc. Amer. Philos. Soc. 1938, V.78, N39, 551.

15. Ландау Л.Д., Лифшиц E.M. Гидродинамика, M., 1988.

16. Бостанджиян С.А., Черняева С.М. О гидродинамическом тепловом "взрыве" неньютоновской жидкости// ДАН, 170, 301, 1966; Бостанджиян С.А., Мержанов А.Г., Худяев С.И. О гидродинамическом тепловом взрыве// ДАН, 163, 133, 1965.

17. Emerman S.II., Turcotte D.L. Stokes's problem with melting// International Journal of Heat and Mass Transfer, 1983, V.26, 1625.

18. Зельдович Я.Б. Избранные труды, т.1. Химическая физика и гидродинамика, М. 1984, 71; ЖЭТФ, 1937, т.7, 1466.

19. Mathews J. Drag force on a slow moving sphere in a medium with variable viscosity// Physics of Fluids, 1978, V.21, 876.

20. Журков C.H., Петров В.А. Физические основы температурной и временной зависимости прочности твердых тел// ДАН СССР, 1978, т. 239, N6, 1316.

21. Трубников Б.А. Закон распределения конкурентов//Природа, 1993, N11, 3.

22. Logan S.E. Deep self-burial of radioactive wastes by rock-melting capsules// Nuclear technology, 1973, V.21, 111.

23. Кащеев В.А., Никифоров А.С. Полуэктов П.П. Поляков А.С. Теория самозахоронения ВАО// Атомная энергия, 1992, т. 73, N 3, 215.

24. PAGIS (Performance Assessment of Geological Isolation Systems for Radioactive Waste), Summary, Commission of the European Communities, 1988, EUR 11775 EN.

25. Ландау Л.Д., Лифшиц E.M. Теория упругости, M., 1965.