авторефераты диссертаций БЕСПЛАТНАЯ  БИБЛИОТЕКА

АВТОРЕФЕРАТЫ КАНДИДАТСКИХ, ДОКТОРСКИХ ДИССЕРТАЦИЙ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ

Математическое моделирование динамики процесса подземного выщелачивания в неоднородном рудоносном слое

МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ им. М.В. ЛОМОНОСОВА ФАКУЛЬТЕТ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И КИБЕРНЕТИКИ Кафедра вычислительных методов

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

Канцель Антон Алексеевич МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ДИНАМИКИ ПРОЦЕССА ПОДЗЕМНОГО ВЫЩЕЛАЧИВАНИЯ в НЕОДНОРОДНОМ РУДОНОСНОМ СЛОЕ Специальность 05.13.18 – Математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ

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

Москва – 2010

Работа выполнена на кафедре вычислительных методов факультета Вычислительной математики и кибернетики Московского государственного университета им. М.В.Ломоносова и в Обществе с ограниченной ответственностью «Интегра Груп.Ру» (г. Москва).

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

Научный консультант:

Куркина Елена Сергеевна доктор физико-математических наук, профессор

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

Елизарова Татьяна Геннадьевна, Институт математического моделирования РАН доктор технических наук, профессор Медведев Александр Сергеевич, ФГОУВПО «Национальный исследовательский технологический университет - Московский институт стали и сплавов» Московский физико-технический институт

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

(государственный университет)

Защита состоится 31 марта 2010 г. в 15 часов 30 минут на заседании диссертационного совета Д501.001.43 при Московском государственном университете им. М.В. Ломоносова по адресу: 119991, Москва, ГСП-1, Ленинские горы, Московский Государственный Университет им. М. В.

Ломоносова, второй учебный корпус, факультет ВМиК, ауд.685.

С диссертацией можно ознакомиться в научной библиотеке факультета ВМиК Московского государственного университета им. М.В. Ломоносова.

Автореферат разослан «» _2010 г.

Ученый секретарь диссертационного совета Д 501.001. доктор физико-математических наук, профессор Захаров Е. В.

-2

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

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

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

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

с ~15% в 90-х годах до 27,7% в 2007 г. (см. [1],[2]) Главной особенностью технологии добычи полезных ископаемых способом скважинного подземного выщелачивания является то, что добыча ведется in situ, без выемки руды из недр. Для осуществления добычи с поверхности в рыхлый обводненный рудоносный пласт, зажатый между водоупорами, пробуриваются технологические скважины. Через группу скважин, называемых закачными, в пласт под давлением подается реагентный раствор, который фильтруется, как правило – в горизонтальном направлении.

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

не происходит нарушения целостности почвы и исключается образование отходов переработки руды. Кроме того, при реализации способа ПВ из производственного цикла исключаются наиболее капитало- и энергоемкие операции, такие как выемка руды из недр, её транспортировка и механическое разрушение (см. [3]).

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

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



Имеет место и обратное явление – образование т.н. «промоин» – устойчивых русел, по которым раствор достигает откачных скважин, не проработав нужного количества руды (см. [5]). И в том и в другом случае снижаются показатели извлечения металла, возрастают дополнительные затраты.

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

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

(либо, Характеристики химико-кинетического взаимодействия раствора реагента и вещества в твердой фазе описываются, как правило, стационарными уравнениями.

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

численное решение уравнений модели гидродинамики без каких либо значимых упрощений и ограничений;

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

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

возможность учета переотложения ранее растворенного металла;

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

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

Научная новизна работы Научная новизна работы состоит в том, что в ней предложена оригинальная модель кинетики химического взаимодействия раствора реагента с минералами урана, находящимися в твердой фазе. Эта модель основана на том, что подземное выщелачивание есть гетерогенный процесс взаимодействия раствора реагента с твердым веществом, протекающий на границе раздела жидкой и твердой фаз. Таким образом, скорость взаимодействия зависит от количества реагирующего вещества и определяется общей площадью контакта жидкой и твердой сред, которая соответствует площади поверхности всех минеральных зерен, обтекаемых раствором реагента в каждый момент времени (см. [5]).

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

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

Предложенная методика является новшеством и не имеет аналогов в исследованных опубликованных материалах.

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

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

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

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

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

Кроме того, разработанная модель и полученные результаты являются фундаментом для решения «обратной» задачи прогнозирования зон локализации вторичных (регенерированных) рудных залежей, образующихся после отработки основных запасов. Гипотеза о существовании вторичного, техногенно природного уранового рудообразования интенсивно исследуется в настоящее время во многих странах (см. например [6],[7]).





-6 Апробация работы и публикации Результаты исследований доложены на:

XI международной конференции Компьютер.

"Математика.

Образование" (Дубна, 2004 г.);

Научно-исследовательском семинаре лаборатории математического моделирования в физике факультета ВМиК Московского государственного университета им. М.В.Ломоносова (2005 г.);

Заседании кафедры вычислительных методов ВМиК МГУ (2009 г.);

Научно-техническом совещании ООО «Интегра Груп.Ру» (2009 г.);

XVII международной конференции Компьютер.

"Математика.

Образование" (Дубна, 2010 г.).

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

Структура диссертации Диссертация состоит из перечня основных сокращений и обозначений, введения, трех глав, перечня основных результатов и списка цитируемых источников. Работа изложена на 127 страницах, включает 32 рисунка. В список цитируемых источников входит 57 наименований.

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

Благодарности Автор выражает сердечную признательность своему отцу и наставнику, доктору геолого-минералогических наук Алексею Викторовичу Канцелю. Автор также благодарит своего научного руководителя д. ф.-м. н. Елену Сергеевну Куркину за постоянное внимание, всестороннюю помощь и поддержку в данной работе, проф. Владимира Ивановича Дмитриева за ценные комментарии и критические замечания и д. ф.-м. н. Николая Владимировича Пескова за консультации в области численных алгоритмов и программной реализации.

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

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

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

Колчин, В. П. Коптелов, Л. П. Лаверов, Л. А. Линцер, Д. П. Лобанов, В. С.

Ломовский, Ю. В. Нестеров, В. В. Новосельцев, С. Н. Пыхарев, Е. А. Толстов, М. И. Фазлуллин, В. Я. Фарбер, P. M. Bommer, R. S. Schechter, L. W. Lake, R. D.

Schmidt, S. E. Follin, M. I. Kabir и многим другим.

Глава I. Структура и основные уравнения математической модели подземного выщелачивания §1 данной главы посвящен обсуждению свойств процесса ПВ как объекта математического моделирования, описанию структуры математической модели динамики процесса ПВ и принятых в ней ограничений и допущений.

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

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

модель гидродинамики процесса;

-8 модель массопереноса и кинетики химического взаимодействия.

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

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

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

Далее, область D, внутри которой проводятся все расчеты, выбирается достаточно большой, чтобы поток раствора через границу был равен (или близок) нулю. Соответственно режим закачки – откачки раствора полагается равнодебетным, т.е. Qi + Q j = 0, где Qi дебеты закачных, а Q j - откачных скважин. Сам раствор мы считаем идеальной жидкостью (т.е. внутреннее трение отсутствует). Течение раствора в пласте считаем напорным и ламинарным.

Известно, что процесс фильтрации раствора подчиняется уравнению неразрывности и закону фильтрации Дарси.

В §2 напоминаются основные сведения из теории фильтрации, используемые в исследовании (см. [8]). Рассматривается фильтрация в однородном пласте при наличии истоков и стоков, затем обсуждается фильтрация в неоднородном пласте.

Т.к. для несжимаемой идеальной жидкости справедливо уравнение неразрывности и с учетом Закона Дарси имеет место:

div( k grad h ) = 0 (1.1) где k коэффициент фильтрации, h функция напора. Выражение (1.1) суть основное уравнение, описывающее установившееся движение жидкости в пористой среде в отсутствии источников и стоков. Компоненты вектора скорости фильтрации потока являются касательными к линиям тока и ортогональны к линиям уровня функции напора h:

h h h V ( x, y,z ) = k,,.

x y z В случае однородного проницаемого пласта процесс фильтрации раствора подчиняется неоднородному двумерному уравнению Пуассона, в правой части которого стоит функция плотности источников и стоков:

n div ( k grad h ) = Q ( x x, y y ), ( x, y ) D i i i 2T i -9 где ( xi, yi ) обобщенная дельта-функция Дирака, взятая в точке забоя i ой скважины, T – мощность пласта, Q1,…, Qn - дебеты скважин, причем дебет закачной скважины берется со знаком «–», а откачной – со знаком «+».

Это уравнение имеет известное аналитическое решение (см. [9],[10]):

n Q ln ( x x ) + ( y yi ) h( x, y ) = (1.2) i i 4 kT i = В этом случае компоненты скорости фильтрации в точке с координатами ( x, y ) определяются по закону Дарси и имеют вид:

n n Q (x x ) Qi ( y yi ) 1 i R 2 i, Vy ( x, y ) = 2 T Vx ( x, y ) = (1.3) Ri 2 T i =1 i = i где Ri расстояние от точки ( x, y ) до i ой скважины.

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

dy V y ( x, y ) Уравнение для линий тока имеет вид: =.

dx Vx ( x, y ) В случае неоднородного пласта аналитическое решение не существует, и расчет поля напора необходимо выполнять численно. Для этого необходимо решить (численно) задачу Неймана для уравнения неразрывности потенциального течения жидкости в пористой среде, которая в общем случае имеет вид:

n ( k ( x, y ) h( x, y ) ) = Q ( x x, y y ) i i i 2T i = (1.4) h ( x, y ) = dn ( x, y ) Г где n внешняя нормаль к границе Г области D 2. В качестве области D, в которой рассматривается система (1.4), выбирается достаточно большая область, размеры которой заметно превышают расстояние между наиболее удаленными скважинами, чтобы скорость течения через границу, вызванная действием скважин, была пренебрежимо мала.

Найдя значение функции напора, легко рассчитать поле скоростей фильтрации: V ( x, y ) = k ( x, y )h ( x, y ).

Необходимо подчеркнуть следующий существенный момент. Как упоминалось ранее, под влиянием процессов растворения и переотложения полезного компонента значение коэффициента фильтрации k ( x, y ) может изменяться со временем, и система (1.4), вообще говоря, должна учитывать эту нестационарность. Однако известно, что значение k ( x, y ) изменяется -10 относительно медленно. Поэтому мы полагаем модель гидродинамики (1.4) квазистационарной. При численных расчетах «обновление» поля скоростей фильтрации производится не на каждом шаге времени, а при существенном (превосходящим заданную величину) изменении коэффициента фильтрации в какой-либо точке области D.

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

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

Такие процессы называются метасоматическими. Огромный вклад в развитие теории метасоматической зональности внес академик Д. С. Коржинский в середине прошлого века (см. [11]). Различают две крайние формы явления метасоматоза: диффузионный и инфильтрационный. При чисто диффузном метасоматозе перенос вещества осуществляется посредством диффузии через застойные поровые растворы. При инфильтрационном метасоматозе компоненты переносятся течением водных растворов, просачивающихся через поры и трещины горных пород. В природе инфильтрационные процессы всегда сочетаются с диффузионными.

Соответственно, считаем, что в процессе фильтрации растворяется один минерал и одновременно выпадает в осадок другой, так что пористость среды не изменяется. (Впоследствии мы учтем медленные изменения коэффициента и коэффициента фильтрации). Обозначим C p концентрация урана в твердой фазе, С – концентрация этого же вещества в жидкой фазе. Уменьшение C p дает увеличение C и наоборот. Обозначим также Cs концентрация урана в насыщенном растворе, Cbk «фоновая» концентрация, имеющая место до начала растворения.

Сначала рассматривается одномерная задача, в которой чистый раствор подается в точке x = 0 и фильтруется с постоянной скоростью в положительном направлении оси x. Из баланса массы легко получить следующее уравнение:

C p ( x, t ) C C C, x ( 0, l (t ) ) = D +V t x x t x (1.5) C ( x,0) = C, C (0, t ) = s где D коэффициент диффузии, V скорость фильтрации раствора, l (t ) движущаяся граница зоны растворения. При x l ( t ) имеем насыщенный раствор, то есть C = Cs. При ПВ, как правило, скорость переноса V много больше скорости диффузии D.

-11 Учитывая гетерогенную природу взаимодействия раствора реагента с веществом в твердой фазе, система (1.5) имеет вид:

C C = (Cs C ), x (0, ), t (0, ) +V t x (1.6) C (0, t ) = Cbk ;

C ( x, t ) Cs, при x Для моделирования процесса массопереноса и химического взаимодействия необходимо найти выражение константы через параметры среды и процесса, а также изучить её зависимость от времени. Для этого введем следующее предположение о строении урановой минерализации в рудоносном слое. Будем считать, что уран находится в виде небольших зерен (минералов), распределенных по всему объему рудоносного пласта с некоторой известной пространственно неоднородной плотностью ( x, y ). Будем также считать, что эти минералы имеют форму шара, и растворение интересующего нас металла в потоке реагента, омывающего минеральное зерно, идет с поверхности с постепенным уменьшением радиуса этого шара. В разделе 3.2 показано, что в этом случае имеет место зависимость:

Cs d = (1 c ), = (1.7) dt r0C p где 0 ( x, t ) 1 - безразмерный радиус минерального зерна ( = r ro, где r0 первоначальный радиус зерна), 0 c( x, t ) 1 - безразмерная концентрация урана в растворе, стационарным образом зависит от состава раствора и других характеристик. Тогда система (1.6), переписанная относительно нормированной концентрации полезного компонента, принимает вид:

c( x, t ) c( x, t ) = (1 c) +V t x, 0 x Vt, t (1.8) ( x, t ) = (1 c ), t 3 Cs где =, а соответствующие граничные и начальные условия – r c(0, t ) = 0, при t c( x, t ) = 1, ( x, t ) = 1, при x Vt (1.9) c( x,0) = 1, ( x,0) = 1, при x Система (1.8) – (1.9) легко обобщается на двумерный случай, как показано в разделе 3.3, и имеет вид:

-12 c( x, y, t ) + ( V, c( x, y, t ) ) = 2 (1 c( x, y, t ) ) t,t tf (1.10) ( x, y, t ) = (1 c( x, y, t ) ) t где t f ( x, y ) время, за которое поток раствора, фильтрующийся вдоль соответствующей линии тока, доходит до точки с координатами ( x, y ).

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

c( x, y, t ) = i i c( x, y, t ) = cbk, t t f ( x, y ) ( x, y, t ) = 1, t t f ( x, y ) где ( xi, yi ) координаты забоя i-ой скважины.

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

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

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

C p ( x, t ) C C C, x ( 0, l (t ) ) = D +V (2.1) t x x t x C ( x,0) = Cs, C (0, t ) = 0 (2.2) где D коэффициент диффузии, V скорость фильтрации раствора, l (t ) движущаяся граница зоны растворения.

-13 Вначале приводятся уравнения диффузионного и инфильтрационного метасоматоза, имеющие известные автомодельные решения (см. [9],[10]). В разделе 1.1 показано, что в случае чисто диффузионного метасоматоза (т.е. при D V ), при скорости растворения Vd = и D = D0 = const исходная задача имеет вид:

C C = D0, x (0, l (t )), t (0, ) t x x (2.3) C (0, t ) = 0, C ( x,0) = Cs Задача (2.3) имеет известное автомодельное решение:

C Cbk 2 4 D x e d + Cbk ;

= t C ( x, t ) = s (2.4) D где Cs концентрация насыщенного раствора, Cbk фоновая концентрация полезного компонента. Процесс растворения идет на границе x = l (t ), которая представляет собой плоскость и движется по известному закону:

l (t ) = 4 Dt (2.5) В случае чисто инфильтрационного метасоматоза ( D V ), рассмотренного в разделе 1.2, исходная задача имеет вид:

C C = 0, x ( 0, l (t ) ), t (0, ) +V t x, (2.6) C (0, t ) = 0;

C ( x,0) = Cs автомодельное решение C ( x, t ) = C ( ), = x Vt, а граница раздела l (t ) = Vt.

Из практики известно (см. [5]), что при закачке раствора реагента и распространении его в рудоносном пласте распределение концентрации урана в растворе имеет некоторые характерные особенности. Имеется зона полного выщелачивания (1), зона активного выщелачивания (2), зона фильтрации насыщенного раствора и Рис. 1 Особенности распределения переотложения (3) и зона фронта концентрации вдоль линии тока фильтрации раствора (см. Рис. 1).

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

Исследование движения переднего фронта приводит к смешанной диффузионно-инфильтрационной задаче, рассмотренной в разделах 1.3 – 1.5. С учетом конечной скорости растворения Vd эта задача имеет вид:

C 2C C D 2 = (Cs C ), Vt L x Vt, +V t x x (2.7) C ( x = Vt, t ) = C ;

C ( x = Vt L, t ) = C, bk s где L ширина зоны растворения, которая определяется скоростью растворения Vd, а константа, определяемая свойствами среды. Эта задача, как показано в разделе 1.5, имеет автомодельное решение вида:

Cs, x Vt L sh / D ( L Vt + x ), Vt L x Vt C ( x, t ) = Cs (Cs Cbk ) (2.8) sh / D L Vt x 0, Найденное решение дает наглядное представление о характере нарастания концентрации полезного компонента на переднем фронте фильтрации, который движется со скоростью фильтрации V.

Исследование движения заднего фронта (зоны активного выщелачивания) позволяет установить связь между скоростью фильтрации и скоростью движения заднего фронта:

1M Vp = V 1 + (2.9) Cs где M константа, характеризующая запасы полезного компонента. Эта формула отражает важное практическое наблюдение, состоящее в том, что скорость движения заднего фронта выщелачивания всегда меньше скорости фильтрации и линейным образом от нее зависит (см. [5]).

Далее, в разделе 1.6 обсуждается исследование динамики распределения концентрации вдоль потока. Для этого исследуется система (2.10) с начальными и граничными условиями (2.11):

c( x, t ) c( x, t ) = ( x, t ) (1 c( x, t ) ) +V t x, 0 x l (t ), t (2.10) ( x, t ) = (1 c( x, t ) ), t -15 c(0, t ) = 0 при t 0;

c( x, t ) = 1, ( x, t ) = 1 при x l (t ) (2.11) c( x,0) = 0, ( x,0) = 1 при x 0.

где, напомним, c( x, t ) - нормированная концентрация урана в растворе, ( x, t ) нормированный радиус минерального зерна, = ( x, t ) пористость 3Cs Cs среды, = константы.

, = r0 r0C p Показано, что автомодельное решение задачи (2.10)-(2.11), описывающее установившееся движение зоны выщелачивания, задается формулами (2.12) и (2.13). Последняя в неявном виде определяет зависимость от автомодельной переменной = x V pt :

c( ) = 3 ( ) (2.12) 1 (1 3 ) d 2 + 1 3 = 6 ln (1 )3 3 arctg 3 = Vp + const (2.13) Эти полученные в диссертационной работе соотношения позволяют сопоставить свойства модели с еще одним важным экспериментальным наблюдением, состоящим в том, что эффективная ширина зоны выщелачивания, отнесенная к скорости фильтрации p V не зависит от скорости фильтрации и обратно пропорциональна коэффициенту растворения (см. [5]):

p x V p 1 M r = ln10 = ln10 (2.14) ( C s + M ) V V V Найденные автомодельные решения и аналитические зависимости позволили провести качественное сравнение и удостовериться в соответствии основных свойств предложенной модели известным экспериментальным данным. Соотношение (2.8) и исследование системы (2.10)-(2.11) позволяют также составить представление о характере изменения концентрации полезного компонента в растворе на переднем (см. Рис. 2(б)) и на заднем (Рис. 2(а)) фронте фильтрации, а также проследить изменение распределения концентрации вдоль потока со временем (см. Рис. 2(в) и (г), на рис. (г) изображены профили концентрации в разные моменты времени). Легко видеть, что демонстрируемые моделью распределения концентрации также согласуется с данными эксперимента о существовании четырех зон концентрации: зоны полного выщелачивания, активного выщелачивания, переотложения и зоны фронта фильтрации раствора.

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

-16 В §3 исследована динамика извлечения урана в продуктивные растворы в среде с неоднородным распределением предельной концентрации урана в насыщенном растворе (насыщенная концентрация). Показано, что в зонах с повышенной насыщенной концентрацией происходит растворение соединений урана, а в зонах с пониженной насыщенной концентрацией – их переотложение и аккумуляция.

а б в г Рис. 2 Характер изменения концентрации полезного компонента в растворе Глава III. Численное исследование динамики ПВ в двумерном случае.

Комплекс программ для ЭВМ В данной главе рассматривается алгоритм численного решения системы уравнений математической модели ПВ и вопросы его реализации в виде комплекса программ для ЭВМ. Обсуждаются результаты расчетов, полученные для некоторых практически важных частных случаев – фильтрация в пористой среде при наличии зоны пониженной проницаемости и наличии зоны со сниженной насыщенной концентрацией. Также рассматривается механизм образования т.н. «тромбов» и «промоин».

-17 §1 посвящен напоминанию основных уравнений модели и описанию алгоритма расчета. Как обсуждалось в Главе I, модель гидродинамики процесса ПВ описывается системой (3.1):

n ( k ( x, y ) h ( x, y ) ) = Q ( x x, y y ) i i i 2T i, ( x, y ) D (3.1) h ( p ) = dn p D Данная система решается численно на основе метода конечных элементов. Предварительно выполняется построение треугольной сетки, которое выполняется методом Делоне.

Решив систему (3.1) численно и найдя значения функции напора h( x, y ) и используя закон Дарси, рассчитаем поле скоростей:

V ( x, y, t ) = k ( x, y, t ) h ( x, y, t ) (3.2) и линии тока:

dy V y ( x, y ) = (3.3) dx Vx ( x, y ) Далее, для каждой точки M ( x, y ) D (точнее, для центра масс каждого треугольника построенной сетки) рассчитаем время t f ( x, y ), за которое раствор, фильтрующийся вдоль линии тока, достигает этой точки.

Математическая модель массопереноса и кинетики химического взаимодействия описывается системой (3.4) относительно нормированной концентрации полезного компонента в растворе и нормированного радиуса ураноносного минерального зерна, с начальными и граничными условиями (3.5):

c( x, y, t ) + ( V, c( x, y, t ) ) = 2 (1 c( x, y, t ) ) t ( x, y, t ) = (1 c( x, y, t ) ) (3.4) t t t f ( x, y ) c( xi, yi, t ) = c( x, y, t ) = cbk, t t f ( x, y ) (3.5) ( x, y, t ) = 1, t t ( x, y ) f Уравнения системы (3.4) справедливы для тех точек области D, куда в текущий момент времени дошел раствор от закачных скважин.

-18 Система (3.4) методом прямых сводится к системе обыкновенных дифференциальных уравнений относительно значений неизвестных функций в узлах сетки. Для интегрирования полученной системы ОДУ используется неявный метод Эйлера либо метод Розенброка четвертого порядка с контролем точности и с переменным шагом по времени. Для предварительных расчетов можно выбрать метод Эйлера, а более точные расчеты выполняются методом Розенброка.

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

Алгоритм расчета Диаграмма 1. Блок-схема расчета динамики процесса ПВ динамки процесса ПВ схематично представлен на Диаграмме 1. Для реализации этого алгоритма разработан комплекс программ для ЭВМ, объединенных интерактивной графической оболочкой. Программы написаны на языке С++ и работают в среде Windows®.

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

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

-19 а- Концентрация U в растворе, t = 2000 б- Концентрация U в твердой фазе, t = в- Кривые извлечения U из откачных скважин (s4…s7) Рисунок 3. Влияние зоны плохой проницаемости на динамику процесса ПВ На Рис. 3(а) оттенками серого цвета показано распределение концентрации урана в растворе спустя 2000 единиц времени. Черный цвет на рисунке соответствует насыщенному раствору, светло-серый – чистому раствору. Хорошо видны очертания более крутого переднего фронта продвижения растворов и пологого заднего фронта – границы зоны полного выщелачивания. На Рис. 3(б) показано распределение урана в твердой фазе.

Светло-серый цвет соответствует низкому содержанию, то есть его полному растворению, а черный цвет– высокому (т.е. где растворение еще не началось).

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

На кривых извлечения полезного компонента (Рис. 3 (в)) мы видим, что к моменту времени t=4000 из двух скважин практически прекратилось извлечение растворенного урана, а из двух других заметное извлечение продолжается. На эти две (центральные) скважины влияет зона плохой проницаемости, процесс -20 извлечения оказывается неравномерным, затянутым;

кривые извлечения широкие, с невысокими максимумами и медленно спадающими «хвостами».

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

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

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

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

n Cs0, ( x, y ) D j, Csn Csj Cs Cs = (3.6) j = C j, ( x, y ) D, j = 1, 2,..., n s j Рассмотрим, как будет происходить процесс фильтрации растворов в среде с переменной насыщенной концентрацией. Вблизи закачных скважин, из которых вытекает чистый раствор, процесс растворения идет быстро, и концентрация в движущемся растворе быстро достигает насыщенной, равной Cs0. Как только концентрация достигла насыщения, процесс растворения прекращается, и фильтрация растворов по линиям тока идет без изменения концентрации. Однако когда раствор достигает зоны с пониженной насыщенной концентрацией урана Cs, которая меньше чем текущая концентрация металла в растворе, равная Cs0, запустится процесс переотложения. Оно будет продолжаться до тех пор, пока концентрация в растворе не уменьшится до значения Cs. По выходе растворов из области D1 опять начинается растворение, и концентрация урана в растворе снова станет увеличиваться.

-21 Уравнения модели (3.4), описывающие процессы метасоматоза в фильтрационном потоке жидкости в среде с переменной насыщенной концентрацией урана (3.6) принимают вид:

c( x, y, t ) + ( V, c( x, y, t ) ) = ± 2 ( Cs c( x, y, t ) ) (3.7) t ( x, y, t ) = ± ( Cs c ( x, y, t ) ) (3.8) t +, если c Cs где, ± константы скоростей ± = (1 / max ), если c Cs растворения и осаждения урана в разных областях, = 3C p константа, характеризующая запасы полезного компонента.

Приведенные выше рассуждения иллюстрируются на Рис. 4.

а- Концентрация U в растворе б- Концентрация U в твердой фазе в- Кривые извлечения U из откачных скважин Рисунок 4. Влияние зоны пониженной насыщенной концентрации полезного компонента в растворе, t = В рассматриваемом случае растворы, пройдя сквозь область D1, где переотложилась часть растворенного металла, быстро дошли до откачных скважин, и концентрация растворенного урана в них не успела восстановиться достигнуть предельного значения. По линиям тока, обходящим зону D1, -22 движение растворов с концентрацией Cs0 продолжается вплоть до откачных скважин. Такая динамика имела место пока задний фронт движения растворов, соответствующий границе полного растворения, где концентрация близка к нулю, не подошла к зоне D1. Только тогда и в этой области началось растворение.

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

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

Обратим внимание, что на основе анализа кривых извлечения можно выделить группу откачных скважин, находящихся в «тени» области D1. Для них кривые извлечения отличаются наличием длинного «хвоста». Расстояние особой зоны до откачной скважины влияет на величину пиковой концентрации на графике извлечения;

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

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

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

Простейшая модель, учитывающая изменение коэффициента фильтрации имеет вид:

-23 k ( x, y, t ) = ( x, y ) (1 c( x, y, t ) ) (3.9) t где ( x, y ) - некоторый коэффициент, характеризующий скорость изменения коэффициента фильтрации. Введение такого коэффициента оправдано, т.к. известно, что различные геотехнологические типы руд, встречающихся при ПВ, могут быть существенно неоднородными по данному показателю. Предполагается, что скорость изменения коэффициента фильтрации также пропорциональна скорости растворения урана (с противоположным знаком) (см. второе уравнение системы (3.4)).

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

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

2. Разработан комплекс программ для расчетов динамики ПВ по этой модели в двумерном случае;

3. В одномерном приближении найдены аналитически автомодельные решения, описывающие движение переднего фронта фильтрации раствора, зоны активного растворения и заднего фронта – границы зоны полного растворения;

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

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

-24 Список цитируемой литературы [1] Forty Years of Uranium Resources, Production and Demand in Perspective, OECD NEA and IAEA, 2006.

[2] Uranium 2007: Resources, Production and Demand, OECD NEA and IAEA, 2008.

[3] В. А. Мамилов, Р. П. Петров, и Г. Р. Шушания. Добыча урана методом подземного выщелачивания. Москва: Атомиздат, 1980.

[1] Перспективы Ядерной Энергии 2008. Резюме для Руководства. Основные Положения, OECD NEA, 2008.

[5] В. И. Белецкий, Л. К. Богатков, и Н. И. Волков. Справочник по геотехнологии урана. Москва: Энергоатомиздат, 1997.

[6] Х. К. Каримов и В. П. Купченко. О возможности природно-техногенного рудоформирования на пластово-инфильтрационных месторождениях (урана) // Узбекский геологический журнал, 1996.

[7] Х. К. Каримов и В. П. Купченко. Природно-техногенный рудогенез на ранее добытых залежах урановорудных месторождений учкудукского (песчаникового) типа // г. Ташкент: 1996.

[8] П. Я. Кочина и Н. Н. Кочина. Гидродинамика подземных вод и вопросы орошения. Москва: Физматлит, 1994.

[9] А. Н. Тихонов и А. А. Самарский. Уравнения математической физики.

Москва: Издательство Московского Университета, 1999.

[10] В. С. Владимиров. Уравнения математической физики. Москва: Наука, 1981.

[11] Д. С. Коржинский. Теория метасоматической зональности. Москва: Наука, 1969.

-25 Публикации автора по теме диссертации [1] A. A. Kantsel’ and E. S. Kurkina. Modeling dissolution in a filtration flow// Computational Mathematics and Modeling, Vol. 17, No. 3, 2006. pp. 211-225.

Канцель А. А., Куркина Е. С. «Моделирование процессов растворения в фильтрационном потоке жидкости»// Прикладная математика и информатика № 21, М.: Изд-во факультета ВМиК МГУ, 2005, С.30-47.

[2] V. I. Dmitriev, A. A. Kantsel’ and E. S. Kurkina. Mathematical modeling of infiltration metasomatism in a porous medium// Computational Mathematics and Modeling, Vol. 19, No. 3, 2008. pp. 239-247. Дмитриев В. И., Канцель А. А., моделирование процессов Куркина Е. С. «Математическое инфильтрационного метасоматоза в пористой среде»// Прикладная математика и информатика № 26, М.: Изд-во факультета ВМиК МГУ, 2007, С. 5-17.

[3] V. I. Dmitriev, A. A. Kantsel’, E. S. Kurkina and N. V. Peskov. Two-dimensional filtration problem for solutions in a nonhomogeneous porous medium// Computational Mathematics and Modeling, Vol. 20, No. 2, 2009. pp. 122-143.

Дмитриев В. И., Канцель А. А.,Куркина Е. С., Песков Н. В. «Двумерная задача о фильтрации растворов в неоднородной пористой среде»// Прикладная математика и информатика № 29, М.: Изд-во факультета ВМиК МГУ, 2008, С. 39-48.

[4] Дмитриев В. И., Канцель А. А., Куркина Е. С. «Математическое моделирование процессов растворения и отложения при фильтрации растворов в пористой среде»// Вестник МГУ, сер. XV, №1, 2009, С. 31-42.

[5] Канцель А. А. «Математическое моделирование процесса подземного выщелачивания»// Тезисы 11-й Международной конференции «Математика.

Компьютер. Образование». г. Дубна, 26-31 января, 2004, с. 54.

[6] Канцель А. А. «Математическое моделирование динамики процесса подземного выщелачивания в неоднородном рудоносном слое»// Тезисы 17-й Международной конференции «Математика. Компьютер. Образование». г.

Дубна, 25-30 января, 2010, с. 125.

-26

 

Похожие работы:





 
2013 www.netess.ru - «Бесплатная библиотека авторефератов кандидатских и докторских диссертаций»

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