Авторефераты по всем темам >>
Авторефераты по техническим специальностям
На правах рукописи
Гурьянов Вадим Владимирович
МАТЕМАТИЧЕСКИЕ МОДЕЛИ РАСПРОСТРАНЕНИЯ ПЛОСКИХ СЕЙСМИЧЕСКИХ ВОЛН В НЕЛИНЕЙНЫХ УПРУГИХ И ФЛЮИДО-НАСЫЩЕННЫХ СПЛОШНЫХ СРЕДАХ
Специальность 05.13.18 - Математическое моделирование, численные методы и комплексы программ
Автореферат диссертации на соискание ученой степени доктора физико-математических наук
Саратов 2007
Работа выполнена в ГОУ ВПО Саратовский государственный технический университет
Научный консультант: доктор технических наук, профессор Крысько Вадим Анатольевич
Официальные оппоненты: доктор физико-математических наук, член-корреспондент РАН Николаев Алексей Всеволодович доктор физико-математических наук, профессор Губатенко Валерий Петрович доктор физико-математических наук, профессор Ерофеев Владимир Иванович
Ведущая организация: Институт теоретической и прикладной механики СО РАН (г. Новосибирск)
Защита диссертации состоится 31 октября 2007 г. в 14:00 на заседании диссертационного совета Д 212.242.08 при ГОУ ВПО Саратовский государственный технический университет по адресу: 410054, г.Саратов, ул. Политехническая, 77, Саратовский государственный технический университет, корп. 1, ауд. 419.
С диссертацией можно ознакомиться в научно-технической библиотеке ГОУ ВПО Саратовский государственный технический университет.
Автореферат разослан "_____"____________ 2007 г.
Ученый секретарь диссертационного совета Терентьев А.А.
ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ
Актуальность работы. Сейсмические волны - колебания горных пород в Земле, возникающие в результате естественных (землетрясения) или искусственных процессов их возбуждения. Изучение сейсмических волн, порождаемых землетрясениями, необходимо для понимания природы землетрясений и их предсказания. С другой стороны, сейсмические волны являются главным источником информации о глубинном строении земных недр.
Длительное время при математическом моделировании сейсмических волновых процессов использовалась линейно-упругая модель. Об этом свидетельствуют научные школы акад. Г.А. Гамбурцева, акад. А.С.
Алексеева, акад. Н.Н. Пузырева, акад. С.В. Гольдина, чл-корр. РАН Б.Г.
Михайленко, проф. И.С. Берзон, проф. Г.И. Петрашень и др. Современная теория распространения сейсмических волн, основанная на линейноупругой модели, изложена в классической монографии К.Аки и П. Ричардса.
В конце прошлого столетия накопилось много экспериментальных фактов, показавших недостаточность линейно-упругой модели описания волновых процессов в горных породах. Эти факты были обнаружены благодаря существенному повышению точности измерения регистрируемых волн различного происхождения: тектонических процессов в Земле и возбуждения колебаний техническими средствами с высокой управляемостью и возможностью контроля излучения.
Было обнаружено нелинейное упругое поведение горных пород при прохождении через них сейсмических волн: обращение волнового фронта, возникновение высокочастотных колебаний (сотни и тысячи герц) на расстояниях 100 - 300 км от очага землетрясения (сейсмическая эмиссия), регулярное изменение спектрального состава волн и другие. Экспериментально показано, что, начиная с некоторого рубежа, пренебрежение нелинейными эффектами приводит к существенным отклонениям решения от истинного явления, - писал член-корреспондент РАН А.В. Николаев - создатель и руководитель научной школы нелинейной сейсмики. В числе первых экспериментальных работ следует отметить работы А.А.Гвоздева, В.В. Кузнецова, Л.Н. Рыкунова, О.Б. Хаврошкина, В.В. Цыплакова.
Работы теоретического плана по нелинейным сейсмическим упругим волнам содержали методы линеаризации и возмущений для получения приближенных уравнений динамики волн (работы И.А. Береснева, Г.М.
Шалашова, Б.Я. Гуревича, А.Л. Литвина, И.Д. Цванкина, Т.З. Вербицкого).
О неупругом поведении горных пород свидетельствуют экспериментальные работы по изучению процессов, предшествующих землетрясениям. Было обнаружено уменьшение скорости распространения сейсмических волн, и этот феномен был связан с образованием зон дилатансии гор ных пород и опубликован акад. М.А. Садовским в 1979 г. Далее был выяснен процесс формирования флюидо-насыщенных резервуаров в результате разрушения горной породы, проникновения в зону разрушения флюидов и последующего геохимического преобразования минералов в этой зоне с образованием диссипативно-дисперсных флюидо-насыщенных резервуаров (В.Н. Николаевский, обзор Дж. Райс). Неупругое поведение этих резервуаров, особенно в зоне их контакта с упругой средой, представляет собой наиболее значительный интерес; оно вызывает дисперсию скоростей в виде связи поглощения с частотой колебаний, необъяснимое уменьшение скорости распространения волн и разрастание амплитуд колебаний (феномен ляркого пятна на сейсмограммах). По этим признакам можно осуществлять поиск и разведку резервуаров. Экспериментальные работы в этом направлении проводились в ОАО Центральная геофизическая экспедиция под руководством доктора физико-математических наук В.Б. Левянта.
Математическим моделям исследования волновых процессов в твердых телах с микроструктурой, представляющих собой диссипативнодисперсные сплошные среды, посвящены научные работы и монографии доктора физико-математических наук В.И. Ерофеева.
Перечисленные факты нелинейного упругого и неупругого поведения горных пород до сих пор являются малоизученными с теоретической точки зрения. Необходимость создания новых математических моделей распространения сейсмических волн в моногенных и гетерогенных горных породах, учитывающих нелинейность и флюидо-насыщенность, была подтверждена на Международной конференции Математические методы в геофизике, прошедшей в 2003 г. в г.Новосибирске в Сибирском отделении Российской академии наук.
Предметом исследований диссертационной работы является математическое моделирование процессов распространения сейсмических волн в нелинейных упругих средах, во флюидо-насыщенных резервуарах и в зоне контакта резервуара с упругой средой.
Цель диссертационной работы: построение математических моделей процессов распространения сейсмических волн в сплошных средах:
нелинейных упругих изотропных и анизотропных, неупругих флюидонасыщенных резервуарах; зонах контакта линейных упругих сред и неупругих флюидо-насыщенных резервуаров, изучение особенностей распространения в них сейсмических волн, объяснение наблюдаемых природных явлений, таких как обращение волнового фронта, происхождение сейсмической эмиссии, регулярное изменение спектрального состава волн, связь поглощения с частотой колебаний, феномен ляркого пятна на сейсмограмме, значительное уменьшение скоростей волн в флюидо-насыщенных резервуарах и появление запаздывающих волн. Сравнение модельных результатов с наблюдаемыми в природе явлениями с целью установления адекватности модели природным явлениям и дальнейшего изучения этих явлений с помощью выбранной модели.
Формулирование и обоснование идеи использования интегральных законов сохранения в прямых задачах количественной сейсмологии с целью упрощения решения этих задач с получением результатов в усредненном виде.
Направление исследований. Построение математических моделей плоских сейсмических волн в нелинейных упругих и флюидо-насыщенных средах. Изучение особенностей распространения таких волн. Сопоставление теоретических результатов с экспериментальными.
Использование интегральных законов сохранения для решения задач распространения волн.
Методы исследований. При решении поставленных задач в диссертации использовались методы математической физики, теории дифференциальных уравнений и методы компьютерного моделирования.
Научная новизна:
1. Введено понятие монотипных плоских упругих волн конечных деформаций: квазипродольных, квазипоперечных и поперечных - волн, распространяющихся независимо друг от друга. Определено условие их существования - направленность векторов внутренних, внешних сил и силы инерции по собственному вектору матрицы уравнений движения. Компоненты собственного вектора содержат нелинейно-упругие параметры среды. В случае линейно-упругой среды квазипродольная волна вырождается в продольную, квазипоперечная - в поперечную, поляризованную с другой поперечной волной. Получено аналитическое общее решение дифференциальных уравнений динамики монотипных волн на основе инвариантов Римана. Дифференциальные уравнения характеристик, описывающие кинематику монотипных волн, решаются численно модифицированным методом Рунге-Кутта. Таким образом, предложен метод решения задачи распространения монотипных волн.
Получены волны Римана как частный случай монотипных волн.
Волна Римана распространяется в одном направлении. Монотипная волна представляется как результат взаимодействия двух волн Римана одинакового типа, распространяющихся в противоположных друг другу направлениях. В нелинейных средах в процессе движения происходит локальное перераспределение энергии (что показано на примере волны Римана), приводящее к образованию ударной волны. Построена ударная волна в виде ее фронта, несущего динамические переменные, заданные инициирующей их известной волной Римана. Локальное перераспределение энергии в распространяющейся волне и образование ударной волны объясняют феномены обращения волнового фронта (возникновение обратной волны в среде без границ разрыва упругих параметров), сейсмической эмиссии, регулярного изменения спектрального состава волн.
Дано обобщение результатов на анизотропные упругие нелинейные среды.
2. Построена линейная математическая модель распространения волн во флюидо-насыщенных резервуарах и зонах контакта их с упругой средой с учетом диссипативно-дисперсных свойств резервуаров. Обосновано использование в модели плоских волн. Модель объясняет наблюдаемые в природе эффекты понижения скорости распространения волн в резервуаре, смещение амплитудного спектра волны в сторону низких частот, затухание волн в резервуарах, разрастание амплитуд колебаний в окрестности границы контакта упругой среды и резервуара (феномен ляркого пятна на сейсмограмме).
Полученные результаты по нелинейным волнам и волнам в флюидонасыщенных резервуарах обобщены в виде модели, которая содержит нелинейную, диссипативную и дисперсную части и является комбинацией уравнений нелинейной упругости, Бюргерса и Кортевега де Фриза. Уравнение Кортевега де Фриза в качестве решения дает солитоны, которые возникают при движении вместо ударных волн и тем сильнее обеспечивают появление очень яркого пятна на сейсмограмме.
3. Предложен принципиально новый подход к построению математической модели явления распространения сейсмических волн в геологической среде. Это явление изучается не в точках среды, а в их окрестностях в усредненном виде на основе использования интегральных законов сохранения механики сплошной среды. Это согласуется с реальной гетерогенностью среды и с тем, что измеряемые механические параметры среды имеют усредненные значения. Предложенный подход формализован в виде математической модели, которая представляется обыкновенными дифференциальными уравнениями вместо уравнений с частными производными и позволяет естественным образом учитывать неоднородности среды типа границ контакта различных горных пород. При этом форму и объем усреднения исследователь может выбирать, исходя из конкретного типа гетерогенности среды.
Разработан численный метод решения смешанных краевых задач для системы обыкновенных дифференциальных уравнений распространения сейсмических волн в рамках предложенной модели. Метод полностью ориентирован на компьютерную обработку с едиными алгоритмами расчетов для линейных и нелинейных волн, гетерогенных, упругих и вязкоупругих сред и т.п., причем условия контакта различных сред удовлетворяются автоматически.
Достоверность результатов работы. Достоверность и обоснованность научных положений и выводов обеспечивается строгим соблюдением законов сохранения механики сплошной среды и определяющих уравнений, а также применением аппарата математического анализа. Сопос тавление модельных результатов с наблюдаемыми в природе говорит о достоверности и обоснованности научных положений и выводов.
Научная и практическая значимость работы. Теоретическая значимость работы заключается в получении результатов, объясняющих наблюдаемые экспериментально явления нелинейного упругого и неупругого поведения горных пород. Это явление сейсмической эмиссии, происходящей за счет образования и распада ударных волн, что доказано на моделях и спектрах вибросейсморазведки, дисперсия скоростей в рыхлых породах, феномен ляркого пятна, объясняемый появлением солитонов и т.д.
Практическая значимость работы заключается в возможности перехода в сейсморазведке от косвенного метода поиска зон скопления углеводородов сейсмическим профилированием к прямому методу сейсмического зондирования, особенно при поиске неструктурных зон скопления углеводородов, что значительно дешевле.
Результаты использовались в совместных работах в ОАО Центральная геофизическая экспедиция (г. Москва) и Нижневолжском НИИ геологии и геофизики (г. Саратов) и используются в учебном процессе при чтении спецкурсов по геофизике, прикладной математике и механике в Саратовском государственном университете.
На защиту выносятся следующие положения:
1. Математическая модель изоэнтропических плоских монотипных нелинейных волн конечных деформаций (квазипродольных, квазипоперечных, поперечных) на основе инвариантов Римана позволяет решать задачи распространения нелинейных сейсмических волн с получением результатов, совпадающих с наблюдаемыми природными явлениями.
2. Математическая модель движения плоских волн в зоне контакта линейных упругой и общей флюидо-насыщенной сред при нормальном падении волн на границу раздела этих сред дает возможность для реализации на ее основе метода адаптивной вибросейсморазведки, при котором режим работы сейсмовибратора определяется по отклику среды.
3. Нелинейная модель движения плоских волн во флюидонасыщенной среде и ее низкочастотное приближение, учитывающая все изученные в рамках данной работы эффекты нелинейного и неупругого поведения горных пород, при прохождении в них сейсмических волн, в том числе и солитоны.
4. Идея использования интегральных законов сохранения в прямых задачах сейсмологии, которая предполагает изучение явления распространения сейсмических волн не в точках среды, а в их окрестностях в усредненном виде. В этом случае не делается переход к дифференциальным уравнениям в частных производных, а для окрестностей точек из интегральных законов сохранения формируется система обыкновенных дифференциальных уравнений по времени, описывающая динамику окрестностей, заключенных в ограниченной области.
5. Разработанный численный метод решения смешанных краевых задач для обыкновенных дифференциальных уравнений движения сейсмических волн в функции времени позволяет решать задачи распространения сейсмических волн единообразно как в однородных средах, так и в средах с границами разрыва параметров, описывающих их свойства.
Апробация результатов работы. Основные результаты исследований, выполненных в диссертации, докладывались на научных конференциях: XXV Гагаринские чтения (Москва, 1999); Проблемы колебаний (Москва, 2001); Petropatch (New Delhi, 2001); Поверхностные волны в анизотропных и сложных средах и обнаружение дефектов (Москва, 2002);
Математические методы в геофизике (Новосибирск, 2003); систематически докладывались на научных семинарах: Института физики Земли РАН, Института вычислительной математики и математической геофизики СО РАН, госуниверситета и технического университета г. Саратова; применялись на договорных началах в ОАО Центральная геофизическая экспедиция (г. Москва) и Нижневолжском НИИ геологии и геофизики (г. Саратов), что отражено в совместных работах; использовались при чтении спецкурсов для студентов в Саратовском госуниверситете.
Публикации. По теме диссертации опубликованы 31 научные работы, в том числе 9 работ в изданиях, рекомендованных ВАК РФ.
Структура и объем диссертации. Диссертация состоит из введения, 5 глав, заключения, списка использованной литературы и приложений.
Общее число страниц ____. Диссертация иллюстрирована ___ рисунками.
КРАТКОЕ СОДЕРЖАНИЕ РАБОТЫ
Во введении представлен обзор основных работ по нелинейному моделированию сейсмических волн в упругих и неупругих средах применительно к диссертационным исследованиям. Сформулированы цели исследования и положения, выносимые на защиту, дано краткое описание работы по главам.
В первой главе поставлена задача распространения плоских нелинейных сейсмических волн и предложен численный метод ее решения. Для постановки задачи сравниваются два классических способа получения дифференциальных уравнений движения: 1) на основе аналитической связи между тензорами напряжений и деформаций, 2) энергетический способ.
Сравнение способов производится на матричных уравнениях представления плоских нелинейных волн и показывается, что матрицы соответствующие различным связям между напряжениями и деформациями имеют совершенно одинаковый вид с точностью до обозначения элементов - функций от производной вектора смещения по пространственной пе ременной. Поэтому с точки зрения теоретического изучения общих свойств волновых процессов безразлично, какая связь между напряжениями и деформациями имеет место.
В п.1.1 диссертации приведены уравнения движения плоских волн на основе связи между тензорами напряжений и деформаций.
В п.1.2 принят энергетический подход как более прозрачный, компактный, обеспечивающий изучение волн конечных деформаций. Использована идея К.Фридрихса о переходе от дифференциальных уравнений второго порядка к дифференциальным уравнениям первого порядка.
Введены следующие обозначения:
ui u = (u1,u2,u3) - смещение, p = ( p1, p2, p3), pi =, x u1 u2 u1 2 p = p1,h = ( p2 + p3 ), q =,,, t t t W = W ( p, h) - функция энергии динамического деформирования, W = 0U ( p,h,S) |S =0, где 0 - начальная плотность среды, U - внутренняя энергия среды, S - энтропия.
Wpp p2Wph p3Wph 1 A(p) = p2Wph p2Whh +Wh p2 p3Whh . (1) 0 p3Wph p2 p3Whh p3Whh +Wh В принятых обозначениях уравнения движения записываются в матричной форме для вектора смещения в виде utt - A(ux)uxx = или для вектора градиента смещения p и вектора скорости смещения q в виде qt - A(p)px = 0, (2) qx - pt = 0.
Уравнения системы (2) являются основными уравнениями распространения изоэнтропических плоских волн в изотропной среде. Функция W = W ( p, h) энергии динамического деформирования является заданной функцией для каждой конкретной сплошной среды, x - координата Лагранжа, t - время, нижние индексы x, t у обозначений функций обозначают производные по этим переменным.
В математической формулировке задача распространения волн эквивалентна задаче Коши для системы дифференциальных уравнений гиперболического типа первого порядка (2) с данными Коши t = t0( ), x = x0( ), q(x0( ),t0( )) = q0( ), (3) где p0(), q0() Ч заданные вектор-функции.
По известным вектор-функциям p(x,t), q(x,t) можно найти интегрированием смещение u(x,t).
В силу того, что А Ч симметричная матрица третьего порядка, у нее существуют три различных действительных собственных значения vi2 и три собственных вектора li:
Wh 22 2 v1,2 = (Wpp + + 2hWph ),v3 =, (4) 0 где = (-Wpp + 2hWhh +Wh);
l1 = (1,mp2,mp3), l2 = (-2hm, p2, p3), l3 = (0, p3, p2), (5) Wph где m = m( p,h) =.
+ 2hWph - В п.1.3 система уравнений движения (2) известным методом Куранта приводится к нормальной форме, что позволяет воспользоваться известным численным методом характеристик с использованием интерполяции и метода Эйлера интегрирования уравнений движения по характеристикам (Г.И. Петровский). В п.1.4 этот метод обобщен на метод Рунге-Кутта.
П.1.5 содержит основные уравнения нелинейной динамической теории упругости по монографии Д.Бленда Нелинейная динамическая теория упругости, которые приводятся справочно.
Во второй главе вводится понятие монотипных волн. Это квазипродольные, квазипоперечные и поперечные волны, распространяющиеся независимо друг от друга. Получено точно общее решение дифференциальных уравнений динамики этих волн в виде инвариантов Римана. Это позволило точно решить задачу распространения монотипных волн с учетом и без учета действия внешних сил частного вида.
Получены волны Римана как частный случай монотипных волн.
Волна Римана распространяется в одном направлении. Монотипная волна представляется как результат взаимодействия двух волн Римана одинакового типа, распространяющихся в противоположных друг другу направлениях. В случае линейно-упругой среды монотипная квазипродольная волна вырождается в продольную, квазипоперечная - в поперечную, поляризованную с другой поперечной волной; эти линейные волны распадаются каждая на две не взаимодействующие друг с другом волны, распространяющиеся в противоположных направлениях.
В нелинейных средах в процессе движения происходит локальное перераспределение энергии (что показано на примере волны Римана), приводящее к образованию ударной волны - источника колебаний, который, как всегда, дает волны, распространяющиеся в разные стороны по линии их движения. Построена продольная ударная волна в виде ее фронта, несущего динамические переменные, заданные инициирующей их известной волной Римана. Локальное перераспределение энергии в распространяющейся волне и образование ударной волны объясняют феномены обращения волнового фронта (возникновение обратной волны в среде без границ разрыва упругих параметров), сейсмической эмиссии, регулярного изменения спектрального состава волн.
В п.2.1 излагаются основные положения математической модели монотипных волн.
Выберем векторы px и pt так, чтобы они были направлены по собственному вектору l(p) = (l1, l2, l3) матрицы A(p) (1) системы (2), соответствующему собственному значению v2(p), т.е. px = bx(x, t)l, pt = bt(x, t)l, b(x,t) Ч дважды непрерывно дифференцируемая функция. Тогда p = p(b) Ч решение обыкновенного дифференциального уравнения:
dp = l(p). (6) db Система (2) в этом случае примет вид qt - v2lbx = 0, (7) qx - lbt = 0.
В зависимости от направления вектора l система (6)-(7) определяет монотипные волны: поперечные, квазипоперечные или квазипродольные.
Пусть (x, t)=const, (x, t)=const Ч два семейства характеристических кривых системы (7), тогда, переходя к характеристическим переменным и , получим систему уравнений движения в этих переменных:
dk а) x + v(p(b))t = 0, в) (qk +k (b)) = 0, д) = v(p(b))lk (p(b)), db (8) б) x - v(p(b))t = 0, г) (qk -k (b)) = 0, (k =1,2,3, lk 0).
Проинтегрировав уравнения (8в) и (8г), получим инварианты Римана k(), k():
k (b(,)) + qk (,) = k ( ), (9) k (b(,)) - qk (,) = k () не изменяющие своих значений на характеристиках (x,t)=const и (x,t)=const.
Инварианты Римана имеют важнейшее значение в изучении монотипных волн, поскольку представляют собой общее решение динамических уравнений (8в) и (8г) с произвольными функциями k(), k(). Эти произвольные функции выбираются в зависимости от конкретной решаемой задачи.
Например, решение задачи распространения сейсмических волн с данными Коши ==, b(,)=b0(), q(,)=q0() на нехарактеристической кривой x = x0(), t = t0() имеет вид:
qk (,) = (qk ( ) + qk () +k (b0( )) -k (b0()), (10) -1 0 b(,) =k 1 (qk ( ) - qk () +k (b0( )) +k (b0()), -где k - функция, обратная k.
По функции b(,) определяются p=p(b(,)),i(b),qi(,) (i k, i = 1, 2, 3).
Переход к переменным x,t осуществляется в процессе решения уравнений (8а) и (8б). Из них для t=t(,) получается линейное уравнение 2 t + (ln v) t + (ln v) t = 0. (11) Однако эффективнее другой способ решения уравнений характеристик (8), удобный для численной реализации с использованием метода Рунге-Кутта.
Введем в рассмотрение функции 1 и 2:
w1(,) = x(,) -vt(,), (12) w2(,) = x(,) + vt(,).
Тогда уравнения характеристик (8а) и (8б) для этих функций запишутся в виде системы обыкновенных дифференциальных уравнений dw1 v = (w1 - w2) (13) d 2v вдоль характеристик =const, dw2 v = (w2 - w1) (14) d 2v вдоль характеристик =const.
Таким образом, численное решение уравнений (13) и (14) дает решение кинематической части задачи распространения монотипных волн, решение динамической части определяется инвариантами Римана.
Монотипные волны представляют собой результат взаимодействия волн Римана (простых волн), бегущих в противоположных направлениях по оси 0x. Теория этих волн хорошо известна. Наиболее компактное их представление описывается в п.2.2 инвариантами Римана (9):
для волны, бегущей в направлении оси 0x:
() =0, p=p0(), q=q0(), b=b0(), x=x0()+v(p0())(t - t0());
для волны, бегущей в противоположном направлении:
() =0, p=p0(), q=q0(), b=b0(), x=x0()-v(p0())(t - t0()).
x = x0(), t = t0() - кривая, на которой задаются данные Коши при ==.
В п.2.3 рассматриваются продольные изоэнтропические ударные волны.
Для построения ударной продольной волны система уравнений продольной монотипной волны записывается в дивергентной форме:
qt -v2( p) px = qt - ( p)x = 0, qx - pt = 0, ( p) = v2( p) и осуществляется переход к характеристическим переменным:
x + v( p)t = 0, x + v( p)t = ( p(,)) + q(,) =( ), ( p(,)) - q(,) = (), ( p) = v( p), ( p) = v2( p).
Принимается, что волна Римана бежит в направлении оси 0x и x = x0(), t=t0(), () =0, p=p0(), q=q0().
Характеристики - однопараметрическое семейство прямых x = + v()t. (15) Для краткости кодификации сложных функций имена функций не меняются, поэтому v( p()) = v() и т.п.
Рис. 1. Схема расположения ударной волны В общем случае прямые, описываемые соотношением (15), непараллельны. Есть зоны их сходимости и расходимости с ростом времени t (см.
рис. 2). Пусть для определенности прямые (15) сходятся для =- s (рис.1).
По теореме о существовании неявной функции =(x,t), определенной равенством (15), эта функция будет неоднозначной при условии 1+ v'()t = 0, (16) которое определяет огибающую С характеристик (15) с параметром =- в виде t = tc() = -1/ v'(), x = xc() = + v()(tc() - ts). (17) Таким образом, если исходить из картины на рис.1, то функции p(-) и q(-) будут однозначными всюду, кроме клиновидной области D, ограниченной огибающей С и характеристикой R, проходящей через точку s. В этой области возможно нарушение однозначности. Точка S (острие клина) на рис.1 - точка начала разрыва функций p() и q() вследствие неоднозначности определения =(x,t). Координаты точки S x = xc(), t = tc(), являются решением трансцендентного уравнения v''() = 0. Характеристика R имеет уравнение x =s + v(s)t.
Далее известным способом определяются локально соотношения для функций, определяющих динамику ударной волны в случае p(+) = q(+) 0.
Ударная волна располагается на предполагаемой кривой W:
x = xw( ), t = tw( ) и на ней известны соотношения p = p( ), V = ( p) / p, q =- pV. Их связь с кривой W следующая:
' ' V = dx / dt = xw /tw.
Для выхода за пределы локальности определим кривую W в неявной форме (x,t) = const. Для V = V ( ) будем иметь дифференциальное уравнение V ( ) + = 0, характеристики которого удобно взять такими:
x t x = +V ( )t. Они отличаются от характеристик (15) только скоростью V ( ). Связь между параметрами и определяется ts = 0 и равенством =- для p(+) = q(+) 0.
Известно, что огибающая однопараметрического семейства решений дифференциального уравнения в частных производных первого порядка также является решением этого уравнения, поэтому параметрические уравнения кривой W определяются так же как кривая С:
t = tw( ) = -1/V '( ), x = xw( ) = +V ( )tw( ), =-, разница только в скоростях V и v.
На ударной волне инварианты Римана принимают вид ( p) =( p) - pV ( p), ( p) =( p) + pV ( p), ( p) 0.
Поэтому за счет скачка смещения и образования скорости V < v волна q Римана порождает монотипную волну с 0 и 0, то есть за ударной волной образуется отраженная от нее волна, бегущая в противоположном оси 0х направлении (ее характеристики располагаются в области F на рис.1). Это явление согласуется с феноменом обращения фронта волны, наблюдаемым в натурных экспериментах. Кроме этого, в зонах сгущения характеристик происходит повышение плотности энергии, приводящее к генерации высоких частот с образованием кратных частот в зависимости от характера нелинейности среды (см. рис.4 и 5). Разрыв p и q порождает весь спектр частот, в дальнейшем эти разрывы сглаживаются в результате диссипативных процессов в реальной среде, в спектре начинают преобладать низкие частоты.
Образование ударной волны и ее движение описывают феномен сейсмической эмиссии.
Для того чтобы монотипная волна оставалась таковой при действии внешних сил (п.2.4), их вектор F должен быть коллинеарен собственному вектору: F = f (x,t)l. Инварианты Римана в этом случае такие:
k (b(,)) + qk (,) = f (x(,),t(,)) k ( p)td = k ( ), k (b(,)) - qk (,) = f (x(,),t(,)) k ( p)t d = k ().
Интегралы вычисляются численно совместно с решением дифференциальных уравнений характеристик.
В п.2.5 дано обобщение полученных результатов на анизотропные среды.
Пусть U - внутренняя энергия деформируемой идеально-упругой однородной сплошной среды при нулевом значении энтропии. Не уменьшая общности, можно считать, что U является полиномом в случае изотропной среды - инвариантов (I1, I2, I3), а в случае анизотропной среды Ч компонентов (ij) тензора конечных деформаций, записанного в декартовой лагранжевой системе координат (x, y, z). Коэффициенты полинома в упомянутых случаях являются упругими константами.
Функция энергии динамического деформирования для плоских волн W = 0U =W(p1,p2,p3)=W(p) в отличие от случая изотропного тела, когда W=W(p1,h(p2,p3)).
Система уравнений, описывающих явление распространения безударных плоских волн, записывается так же, как для изотропной среды:
qt - A(p)px = 0, (18) qx - pt = 0.
1 2W Только в этой системе A(p) = Ч квадратная симметрич 0 pipj ная матрица третьего порядка, у которой все ее собственные значения (k, k=1,2,3) вещественные и положительные, чем и обеспечивается гиперболичность системы. В общем случае волна u(x,t) состоит из одной квазипродольной и двух квазипоперечных волн. В изотропной среде одна волна обязательно поперечная.
Уравнения (18) ничем не отличаются по форме от уравнений для изотропной среды. Коренное отличие заключается в том, что при повороте системы координат (x, y, z) и переходе при этом к системе (xТ, yТ, zТ) коэффициенты полиномов W=W(p1, p2, p3) и W=W(pТ1, pТ2, pТ3), где pi =ui / xi, будут различны, что и приводит к различию скоростей v и vТ. В изотропной среде скорость v не зависит от направления распространения плоской волны.
Таким образом, при качественном изучении особенностей распространения упругих волн в анизотропных средах можно пользоваться методами, разработанными в данной главе для изотропных сред. При количественном изучении нужно иметь в виду, что упругие константы зависят от принятого направления оси 0х и это скажется на скоростях распространения волн vk и собственных векторах lk (k=1,2,3). В этом вся сложность изучения явления распространения волн в конкретных анизотропных средах.
Следует отметить, что для анизотропных сред дифференциальные уравнения всегда строго гиперболические в изоэнтропическом случае в отличие от изотропных сред, для которых дифференциальные уравнения перестают быть строго гиперболическими.
В третьей главе рассмотрены частные модели сред со степенной нелинейностью, которые заданием конкретного вида функции энергии динамического деформирования получаются из модели, предложенной в главе 2. Рассматриваются среды с квадратичной нелинейностью (модель Мурнагана) и с кубичной нелинейностью. Модель среды с кубичной нелинейностью построена, исходя из предположения, что сжатие элемента среды равно его растяжению, взятому с обратным знаком, т.е. диаграмма растяжения-сжатия имеет вид кубической параболы. Получены формулы для решения динамической части задачи распространения монотипных волн, формулы для волн Римана. Для получения решения кинематической части задачи распространения монотипных волн, т.е. решения уравнений характеристик, предложена модификация метода Рунге-Кутта и описан алгоритм решения. Приведены в графическом виде результаты компьютерных расчетов по этому алгоритму с помощью созданной программы - моделирование процесса распространения монотипной волны с формированием из нее волн Римана. На полученных картинках представлено изменение формы волны, зоны сгущения характеристик - зоны возникновения ударных волн. Получены спектры продольных волн Римана в средах с квадратичной и кубичной нелинейностью. Показано совпадение полученного теоретически спектра для среды с квадратичной нелинейностью со спектром сигнала, зарегистрированного в результате эксперимента в поле.
В п.3.1 рассмотрены разделенные волны.
Положим Wph = 0 в матрице A. Тогда W ( p,h) = 0 представляется в виде W ( p,h) =W1( p) +W2(h). (19) Волны, как монотипные, разделяются на продольную и две поперечные, распространяющиеся с различными скоростями v1, v2, v3.
2 В случае W2 = 0 имеем v2 = v3 = / 0, т.е. одинаковые скорости линейных поперечных волн.
В п.3.2 приведены примеры волн в конкретных средах.
Функция энергии динамического деформирования (W ) разлагается в обобщенный ряд Тейлора по инвариантам. Для плоских волн получается представление:
1 1 1 1 1 1 W = ( + )m2 + n2 + 0 ( Am3 + Bmn2 + Lm4 + Mm2n2 + Pn4 ), (20) 2 2 6 2 12 12 2 2 2 в котором m = p1 = ux1, n2 = p2 + p3 = ux + ux, А,B,L,M,P - реологические 2 константы.
В п.3.2.1 рассмотрен пример среды с квадратичной нелинейностью (Мурнагана).
Примем в (20) следующие обозначения:
2 + 2 =, 0 A = , 0B = , m = p1 = p, n2 = p2 + p3 = 2h.
В этих обозначениях частные производные функции энергии деформаций (W ) будут такими:
Wpp = + p, Wph = , Whh = 0, Wh = + p.
Выбор знака у модуля определяется из экспериментов.
При = 0 в соответствии с п.3.1 нелинейной будет только продольная волна (v2 = ( + p) / 0).
Принято p(b) = b и в соответствии с п. 2.2 / + p 2 ( + p)3/ 2 1 3(,) -, (b) = ( p) = dp =, p() = 1/ 0 3 20 q(,) = (q0() + q0() +( p0()) -( p0())), (21) (,) = (q0() - q0() +( p0()) +( p0())) при данных Коши = = , x = x0 ( ), t = t0 ( ), p = p0 ( ), q = q0 ( ).
Формулы (21) определяют решение динамической части задачи Коши.
Приведен конкретный вид волн Римана.
В п.3.2.2 рассмотрен пример среды с кубичной нелинейностью, для которой поставлено требование, чтобы диаграмма растяжения удовлетворяла условию: сжатие есть растяжение с обратным знаком. Вследствие этого условия функция энергии деформаций W ( p,h) должна быть четной по p.
Тогда Wpp = + p2 + h, Wph = p, Whh =, Wh = + p2 + h. Здесь введенные в п. 3.2 обозначения 1 + 2 =, 0L = , 0M = , 0P = , A = B = 0.
3 Уравнение, описывающее продольные волны, имеет вид ( +ux ) uxx -utt = 0. При данных Коши t = 0, x = , u = u0( ) для волны Римана, бегущей в направлении оси 0x получаются формулы:
u(,t) = u0( ) + (1 u0( )v - Q)t, x = + vt, arcsin( u0( )) < Q =.
2 0 Arsh( u0( )) > Для волны, бегущей в противоположном направлении, нужно сменить знак у скорости v.
В п.3.3 описан алгоритм решения задачи распространения монотипной волны (задачи Коши). Решение задачи разделяется на динамическую и кинематическую части. В п. 2.1 главы 2 получено решение (10) динамической части задачи Коши. Для получения кинематики (характеристик) предложено решать обыкновенные дифференциальные уравнения вдоль характеристик (13), (14). Для этого случая модифицирован обычный метод Рунге-Кутта и дано его описание.
В п.3.4 представлены результаты компьютерного моделирования распространения монотипных волн. На рис. 2 представлены характеристики квазипродольной волны в среде с кубичной нелинейностью. Видно, что характеристики близки к прямым в зоне взаимодействия волн и прямые в зоне их превращения в волны Римана. Это свойство характеристик обусловило численный метод, описанный ранее.
На рис. 3 представлена продольная составляющая квазипродольной волны в среде с кубичной нелинейностью. Виден переход от взаимодействия волн к их независимому движению в виде волн Римана.
x t Рис. 2. Характеристики квазипродольной волны в среде с кубичной нелинейностью Рис. 3. Продольная составляющая квазипродольной волны в среде с кубичной нелинейностью В п.3.5 исследуются спектры продольных волн Римана в средах со степенной нелинейностью. На рис. 4 приведен теоретический амплитудный спектр продольных волн Римана для среды с квадратичной нелинейностью. Спектр представлен в логарифмическом масштабе. Спектр волны без учета нелинейности представлен пунктирной линией.
Этот спектр получен на 10 лет ранее полевого спектра на рис. 5 в нормальном масштабе. Полное совпадение кратных частот 2, 3, 4 исходной частоты определяет выбор модели с квадратичной нелинейностью для описания данной геологической среды.
Рис. 4. Спектр сигнала для среды с квадратичной нелинейностью Рис. 5. Спектр сигнала по результатам полевых наблюдений На рис. 6 приведен амплитудный спектр для среды с кубичной нелинейностью. Видно возникновение кратных частот , 3, 5 и т.д.
Рис. 6. Спектр сигнала для среды с кубичной нелинейностью В четвертой главе, носящей прикладной характер, описаны математические модели распространения сейсмических волн во флюидонасыщенных резервуарах и зонах их контакта с упругими средами.
Различают два случая неупругого поведения гетерогенной двухфазной геологической среды:
а) взаимодействие деформационно-диффузионных процессов движения жидкости в горных породах, подверженных дилатансии и образованию зон деструкции за счет флюидо-геохимических преобразований минералов с движением флюидов. В результате взаимодействия образуется диссипативно-дисперсная среда. Этот процесс очень медленный и длится столетиями и более. Важно, что при математическом моделировании описанного взаимодействия используется шаровая часть тензоров напряжений и деформаций и закон фильтрации Дарси. В среде происходят неупругая деформация и релаксация давления;
б) кратковременное волновое воздействие на геологическую среду, применяемое в сейсмике, когда изучаются сейсмические волны во флюидо-насыщенной среде. В этом случае для математического моделирования волновых движений в такой среде используют девиаторы соосных тензоров напряжений и деформаций. В этом случае квадраты интенсивностей касательных напряжений и деформаций сдвига пропорциональны вторым инвариантам соответствующих девиаторов.
Вследствие существенно разных временных масштабов процесса б на фоне процесса а, последний можно считать стационарным и не учитывать.
Если релаксационные процессы, происходящие в сплошной среде, обусловливаются взаимодействием между упругой (твердой) и вязкой (жидкой) фазами квазиизотропной среды, то соответствующие соотношения между напряжениями и деформациями называются вязкоупругими. В зависимости от того, определяется ли при таком взаимодействии доминирование твердой или жидкой фазы, простые вязкоупругие среды называются упругозапаздывающими или релаксирующими.
Простейшая упругозапаздывающая среда - КельвинаЦФойгта.
Среда с линейной релаксацией - Максвелла.
Приведено общее линейное дифференциальное уравнение описания механических движений во флюидо-насыщенном резервуаре в диссипативно-дисперсном состоянии. Принято доминирование твердой фазы, и для изучения выбрано волновое уравнение для диссипативно-дисперсной среды. Поскольку поперечные волны очень быстро затухают (опыты Н.Н.
Пузырева), то рассматриваются только продольные волны.
Для изучения свойств решения использован метод разделения переменных. Показано, что все особенности, связанные с диссипативными и дисперсными процессами, описываются уравнением для функции времени в разделенных переменных. Уравнения для функции пространственных переменных одинаковы для любых сред. Поэтому изучаются только плоские волны.
Смоделирована природная ловушка для флюидов (резервуар) и рассмотрены случаи отражения волн от окрестности резервуара (от упругой среды) и от разных зон внутри резервуара. Показано, что диссипативнодисперсные свойства среды влияют на уменьшение скорости распространения волн, на сдвиг частот в сторону низких в амплитудном спектре волны, на увеличение амплитуды отраженных от резервуара волн (феномен яркого пятна на сейсмограмме) и на появление запаздывающих по фазе волн. Тем самым обеспечивается адаптивный метод вибросейсморазведки, т.е. подбор частоты вибратора в зависимости от свойств исследуемой среды по параметрам отраженной волны. Эти результаты были получены в рамках работы по заказу ОАО Центральная геофизическая экспедиция при участии ее представителя В.Б. Левянта.
Полученные ранее результаты по нелинейным волнам и волнам во флюидо-насыщенных резервуарах обобщены в виде модели, содержащей нелинейную, диссипативную и дисперсную части и являющейся комбинацией уравнений нелинейной упругости, Бюргерса и Кортевега де Фриза, которые хорошо изучены. Уравнение Кортевега де Фриза в качестве решения дает солитоны, которые возникают при движении вместо ударных волн и тем сильнее обеспечивают появление очень яркого пятна на сейсмограмме.
В п.4.1 описана линейная математическая модель флюидонасыщенных сплошных сред.
Общий случай упругосжимаемой линейной вязкоупругой среды описывается определяющими уравнениями 2 P = 2R + (( + )P - R)divuE - (3 + 2)PQE (22) 3 Здесь P и R Ч линейные дифференциальные операторы m n P = a0 + a1 +... + am m, R = b0 + b1 +... + bn n, t t t t ak (k = 0,1,..., m), bk (k = 0,1,...,n) - материальные (реологические) константы, являющиеся комбинациями времен релаксации и модулей сдвига, - тензор напряжений, - тензор малых деформаций, E - единичный тензор, u - вектор смещений, Q - абсолютная температура, , - константы Ламе, - коэффициент объемного температурного расширения, t - время.
Следует подчеркнуть, что уравнение (22) содержит только равенства девиаторов тензоров напряжений и деформаций с включением влияния объемного температурного расширения. Учитывая кратковременность волнового воздействия на вязкоупругую среду, принимается Q = const и в дальнейшем P =1, т.е. в среде доминирует твердая фаза в дисперсном состоянии.
В п. 4.2 рассмотрена задача распространения волн в зоне контакта резервуара с упругой средой.
Упругозапаздывающая диссипативная среда КельвинаЦФойгта опре деляется оператором R = 1+ ( Ч время запаздывания).
t При плотности массовых сил F = f + g уравнения движения в потенциалах и A такие:
(( + 2) + t ) - tt + f = 0, (23) (A +At) - Att + A = 0, divA =0. (24) Здесь Ч плотность, = Ч коэффициент вязкости.
Векторный потенциал описывает поперечные волны. Экспериментальными работами Н.Н. Пузырева было показано их быстрое затухание, и поэтому их изучение опускается.
К этим уравнениям можно применить метод разделения переменных.
Для скалярных волн (23) при f = 0 полагаем = (x, y, z)T (t).
В результате применения метода разделения переменных получаем уравнения + k2 = 0, (25) k2 T + (( + 2)T + T ) = 0.
При = 0 уравнения (25) превращаются в уравнения колебаний упругой среды, т.е. первое уравнение (25) одинаково как для упругой, так и для вязкоупругой сред. Свойства решений этого уравнения достаточно хорошо изучены. Второе уравнение (25) описывает специфику распространения волн в вязкоупругой среде и по форме совпадает с уравнением T + 2nT + a2T = 0 (n > 0, a > 0) (26) затухающих упругих колебаний классической теоретической механики.
Уравнение для смещения u = u(x,t), описывающее движение плоских волн в направлении оси 0x представляется в виде v2uxx + 2buxxt - utt = 0. (27) Здесь v2 = ( + 2) , b=(2) (3), = (1- kп)M, M - плотность среды до дилатансии (плотность монолитной горной породы), kп Ч коэффициент пористости, - вычисляемая плотность вязкоупругой среды.
Дифференциальное уравнение (27) при b 0 превращается в уравнение плоских продольных волн, при v 0 Ч в уравнение диффузии, т.е. оно соединяет в себе упругие и диффузионные свойства.
Метод разделения переменных дает:
u = X (x)T(t) (28) X + k X = 0 (29) 2 T + 2bk T + k v2T = 0 (30) Заметим, что уравнение (29) не зависит от коэффициентов v2 и 2b исходного уравнения (27), и его общее решение запишем так:
X (x) = Asin(kx + ), (31) в этом решении A, k, Ч пока произвольные константы.
Для исследования свойств решений уравнения (30) введем обозначения n = bk, a = kv и получим уравнение (26). Оно описывает свободные упругие колебания при наличии сопротивления, а при n 0 Ч свободные гармонические колебания. Его общее решение зависит от знака разности a - n и представляется так:
при a - n > 0 : T = Be-nt sin(t a2 - n2 + ) (32а) при a - n = 0 : T = e-nt (A1t + A2) (32б) при a - n < 0: T = A1exp(-(n - n2 - a2 )t) + A2 exp(-(n+ n2 - a2 )t) (32в) Здесь A1, A2, B, Ч произвольные константы.
Случаи (32б) и (32в) не являются колебаниями, а представляют собой процесс диффузионного рассеяния (диссипации) механической энергии. На временном разрезе регистрации результатов наблюдений в сейсморазведке случаи (32б) и (32в) дают нерегулярные помехи. Сосредоточим свое внимание на случае (32а), т.е. принимаем T = Be-nt sin(t a2 - n2 + ), (a - n > 0). (33) Общее решение исходного уравнения (27) представляется суперпозицией бегущих волн вида x x Xk (x)Tk (t) = e-nt C1 cos( (t )) + C2 sin( (t ). (34) VV Здесь V = k. Из (27), (33) и (34) следует, что угловая частота колебаний меньше, чем частота свободных гармонических колебаний (идеально упругая среда), амплитуда колебаний уменьшается за счет увеличения как параметра k, так и времени t :
= a2 - n2 = kV < kv. (35) При увеличении параметра n увеличивается затухание колебаний (поглощение), чем и описывается диссипативность среды.
Рассмотрим в окрестности границы ( x = 0 ) раздела упругой и вязкоупругой сред явление отражения и преломления при нормальном падении на эту границу продольной упругой волны. Считаем известными следующие константы, характеризующие среды: упругие - 1, 1, 1; п - принятый в сейсморазведке коэффициент поглощения; вязкоупругие - , , , .
Используя эти константы, вычисляем n 1 +21 +2 2 v1 =, v2 =, b =, n = пv1, k =, V = v2 - b2k2, = kV.
3 1 b Основные константы подобраны так, чтобы можно было согласовать отдельные затухающие гармоники в контактирующих средах с тем, чтобы легче было разобраться в физической сущности изучаемого явления.
Известную падающую (uп ) волну при условии, что uп 0 перед фронтом волны t - x v1 (т.е. среда находится в покое), искомые отраженную (uо ) и преломленную (uпр ) волны представим так:
x x e-n(t - x v1 ) sin((t - )) t, x : t - > v1 v uп = (36) 0 t, x : t - x v x uо = e-n(t + x v1 )((t + )), vx x uпр = e-nt C1 cos((t - )) + C2 sin((t - ).
VV Следует отметить, что падающая волна - не импульс, а затухающая синусоида, определенная на координатной полуоси, поэтому она имеет большую протяженность во времени.
Функция определяется из условия непрерывности смещений на границе x = uп + uо = uпр и имеет вид (t) = C1 cost + (C2 -1)sint.
Константы C1 и C2 определяются из условия непрерывности нормальных составляющих тензоров напряжений упругой и вязкоупругой сред:
(1 + 21)(uп, x + uо, x ) = ( + 2)uпр, x + uпр, xt.
Так как cos = -sin( - 2), отраженная волна в окрестности границы отраженияЦпреломления приобретает вид uо = e-nt C sint - C1 sin(t - ) (C1 > 0), т.е. состоит из двух волн: первой Ч основной и второй Ч запаздывающей по фазе на 2. Здесь уместно заметить, что если время запаздывания равно нулю, то и вязкость среды ( ), поглощение n и коэффициент Cбудут нулевыми, следовательно никакой диссипации энергии не будет.
Таким же образом рассматривается движение волн из вязкоупругой среды в упругую с плотностью M и скоростью распространения волн vM.
В п. 4.3 описана математическая модель адаптивного метода вибрационной сейсморазведки.
Получена связь между падающей нормально к границе вязкоупругого слоя плоской продольной волной в виде затухающей синусоиды и волной, отраженной от этого слоя. Вязкоупругий слой обладает не только диссипативными свойствами, но и дисперсными. Это приводит к уменьшению скорости распространения и поглощению волн в слое. Рассмотренная модель ориентирована на адаптивную форму вибросейсморазведки, поскольку можно определять по форме возбуждаемой волны форму регистрируемой отраженной от слоя волны. Это позволяет подобрать режим работы виброустановки.
Уравнение движения берем в форме ( + 2)uxx + (1uxxt +2uxxtt ) - utt = 0. (37) Здесь введен параметр 2, учитывающий дисперсность среды.
Применив к уравнению (37) метод разделения переменных (u(x,t) = X (x)T (t) ), получим, как обычно, два уравнения:
X + k X = 0, T + 2nT + a2T = 0 (n > 0, a > 0).
Здесь:
1 k2v2 2 + 2 2 n =, a2 =, = k2vs, v2 =, vs =, (38) 1+ 22 2 p 1+ 22 v2, vs - квадраты скоростей продольных и поперечных волн. Параметр k p разделения переменных - волновое число.
Решение u(x,t) уравнения (37) представляется в форме бегущих простых волн u = e-nt A1 cos(t kx)) + A2 sin(t kx)). (39) () Из формул (38) и (39) следует соотношение = a2 - n2 = kV < ka, и влияние реологической константы такое, что с ее ростом скорость V уменьшается, следовательно уменьшается частота колебаний . Это объясняет сдвиг частоты в сторону низких частот, увеличение амплитуды отраженных от резервуара волн (феномен ляркого пятна) (см. рис.7, 8, 9).
В окрестностях границ x = 0 и x = h вязкоупругого слоя толщиной h, лежащего между упругими средами, рассмотрено явление отражения и преломления при нормальном падении на эти границы продольной упругой волны.
В п.4.4 рассмотрен общий случай упругосжимаемой линейной вязкоупругой среды при постоянной температуре.
Уравнение движения имеет вид 2 ( + + R)uxx - utt = 0.
3 После разделения переменных u(x,t) = X (x)T (t) для T (t) при выборе констант в операторе R такими, чтобы b0 =1, b1 =1, b2 =2, а остальные константы bj ( j = 3,4,Е) удовлетворяли бы равенству s+k2 2 4 2 + + R T = 2n + a2 T = 0, T + + 3 3 t2 t получается решение s j T (t) = t exp(-nt it a2 - n2 ) (s = 0,1,2,...).
A j j=Здесь Aj - произвольные феноменологические константы, n и a2 определены формулами (38), i - мнимая единица. Феноменологические константы позволяют по форме изучаемых конкретных источников получить их влияние на форму регистрируемых колебаний.
Смещение uk (x,t), соответствующее константе k разделения переменных, имеет вид uk (x,t) = c1(t)cos(t kx) + c2 (t)sin(t kx).
s s Здесь c1(t) = tie-nt, c2 (t) = C C tie-nt, V = / k Ч скорость распро1,i 2,i i=0 i=странения волн в резервуаре, = a2 - n2 = kV - угловая частота, k - волновое число.
Получена формула = - uxtt, упрощающая разрешение условий конk такта упругой и вязкоупругой сред. Эти условия разрешены в общем виде движения из среды упругой в вязкоупругую и наоборот при нормальном падении волн на границу контакта.
В п. 4.5 приведены примеры расчетов волновых полей в зоне контакта флюидо-насыщенных резервуаров с упругой средой.
Для конкретных сред получены результаты в графическом виде для сопоставления форм падающей упругой волны (импульс Берлаге) на границу флюидо-насыщенного резервуара, отраженной волны от него и от монолитной породы в зависимости от характера содержащихся в резервуаре флюидов: вода, переходная зона вода-нефть, нефть. Схематично резервуар изображен на рис. 7. Здесь монолит - горная порода, в которой образовался резервуар, осадочный чехол - слои пород, покрывающие резервуар. Монолит и осадочный чехол обладают упругими свойствами. Скорость распространения волн в монолите больше, чем в чехле.
Осадочный чехол т т Нефть Вода Монолит Монолит 0 800 1600 2000 3000 3400 4200 50Резервуар Рис. 7. Схема резервуара Результаты представлены в виде отдельных кривых (трасс) и в виде временных разрезов, используемых в геофизике. В случае вязкого флюида (нефть) видны запаздывающая по фазе волна и уменьшение ее амплитуды.
Пример представлен на рис. 8. На рис. 9 представлен синтетический сейсмический разрез резервуара, на котором выделяются ляркие пятна - увеличенные амплитуды отраженных от резервуара волн.
Эти результаты были получены в рамках работы по заказу ОАО Центральная геофизическая экспедиция при участии ее представителя В.Б. Левянта. На рис. 10, 11 представлены экспериментальные данные, полученные В.Б.Левянтом: трасса временного разреза, полученного для резервуара с нефтью, и ее амплитудный спектр.
Вода-нефть Вода-нефть Вода Вода Рис.8. Отраженные волны от резервуара Рис. 9. Синтетический сейсмический разрез резервуара Рис. 10. Трасса временного разреза, полученного для резервуара с нефтью (отраженный сигнал) Рис. 11. Амплитудно-частотный спектр сигнала, отраженного от резервуара с нефтью В п.4.6 рассмотрена нелинейная модель волнового движения во флюидо-насыщенном резервуаре.
Нелинейное дифференциальное уравнение с совместным действием нелинейных, диссипативных и дисперсных членов в правой части имеет вид:
2 n utt - uxx = auxuxx + buxxt + c2uxxtt. (40) ( ) p Известным методом получается его низкочастотная аппроксимация в координатах = x - t, T = t :
p 2 qT + aqnq - b q + c2 q = 0. (41) p p p Связь между уравнениями (40) и (41) осуществляется в соответствии с заданием начального профиля волны, т.е. u(x,0) = f (x,0), f '(x,0) = q(x,0).
В зависимости от значений параметров n, c, b получаются дифференциальные уравнения: n =1, c = 0 - Бюргерса для диссипативных сред, c =1, b = 0 - Кортевега де Фриза порядка n для сред в дисперсном состоянии, а при n =1 - классическое уравнение Кортевега де Фриза.
Решением уравнения Кортевега де Фриза являются солитоны, появляющиеся в дисперсной среде вместо ударной волны, возникающей в монотипной среде. Солитоны обладают очень большими амплитудами.
Таким образом, в четвертой главе:
- построена и исследована математическая модель явления распространения сейсмических волн во флюидо-насыщенных резервуарах с доминированием части твердого тела. Все особенности волн заключены в функции T(t), не зависят от геометрии границы резервуара и точно описываются плоскими волнами;
- получены точное общее представление плоских волн, которое можно рассматривать как обобщение импульса Берлаге, и выражение тензора напряжений для случая плоских волн;
- определены границы основных измеряемых параметров модели, определяющие особенности волнового движения: диапазон изменения скорости распространения волн в резервуаре (V), влияние вязкости среды на затухание волн, влияние частоты и скорости на волновое число, важное для представления движения в форме бегущих или стоячих волн при прикладном использовании модели.
Скорость в вязкоупругой среде (V) уменьшается не только за счет увеличения вязкости среды, но и за счет времени дисперсии волн. Вариацией их значений можно получить любую скорость V (в допустимом диапазоне ее изменения), вплоть до близкой к нулю, т.е. получить отсутствие волнового движения в вязкоупругой среде или, иными словами, полное отражение от границы резервуара. Резкое понижение скорости V в случае линейной модели и возникновение солитонов в нелинейной объясняет феномен ляркого пятна.
В пятой главе предложен принципиально новый подход к построению математической модели явления распространения сейсмических волн в геологической среде. Это явление изучается не в точках среды, а в их окрестностях в усредненном виде на основе использования интегральных законов сохранения механики сплошной среды, что ближе к реальности и позволяет естественным образом учитывать неоднородности среды типа границ контакта различных горных пород. Математическая модель представляется обыкновенными дифференциальными уравнениями вместо уравнений с частными производными. Такой подход может быть использован и для изучения других эволюционных процессов в Земле (например, электромагнитных), в основе которых лежат интегральные законы сохранения.
Геологическая среда весьма неоднородна по своему составу и структуре, в сущности гетерогенна. Ее механические параметры имеют усредненные значения в силу способа их измерения. При математическом моделировании волновых процессов в геологических средах по традиции используют дифференциальные уравнения, невольно считая при этом среду сплошной (континуальной).
Между тем в основе механики деформируемого твердого тела лежат интегральные законы сохранения и определяющие уравнения, предписывающие изучать явления распространения волн не в точках, а в их окрестностях в усредненном виде. При этом форму и объем усреднения исследователь может выбирать, исходя из конкретного типа гетерогенности среды.
Следует добавить, что границы разрыва материальных параметров (границы контакта геологических сред с различными горными породами) учитываются при вычислении интегралов, входящих в законы сохранения, и должным образом влияют на процесс распространения волн. Заметим также, что граничные условия контакта, необходимые для построения математической модели в виде дифференциальных уравнений, получаются из законов сохранения. Таким образом, интегральные законы сохранения не только ближе к реальности, но и упрощают постановку задач распространения волн в слоистых средах.
В случае малых деформаций, симметрии тензора напряжений, пренебрежения термодинамикой деформаций и использования энергетически допустимых определяющих уравнений (что сделано для более прозрачного изложения сути предлагаемого численного метода) остается только закон сохранения импульса и определяющее уравнение:
d 1 1 vdg = f dg, n ds + G dt G G (42) G G G = (x, y, z, p1, p2, p3), в которых - плотность среды, v = (v1,v2,v3 ) = ut - скорость смещения u = u(x, y, z,t), pm = um = ( p1m, p2m, p3m ), - тензор напряжений, f - плотность объемных сил, x, y, z - декартовы координаты, G, G, | G | Ч область усреднения, ее граница, объем в пространстве {x, y, z}, n - единичная внешняя к G нормаль, t - время.
Согласование векторов v и pm (m = 1,2,3) возьмем также в интегральной форме, пригодной в случае непрерывности упругих параметров среды:
d pmdg = nvmds, (m =1,2,3). (43) dt G G G G Нормировка объемных интегралов в (42) и (43) произведена с целью показа, что эти интегралы можно рассматривать как математические ожидания функций случайных величин с постоянной плотностью их распределения или как средние функции, широко используемые в математическом анализе. Фактически (42)Ц(43) Ч уравнения для средних значений искомых величин. Соотношения (42)Ц(43) точные и являются основой построения предлагаемого численного метода.
Отнесем к введенной системе координат {x, y, z} открытую связную область D с кусочно-гладкой границей S. В этой области ищется распределение усредненных величин v и pk (k = 1,2,3). По ним, в случае необходимости, можно вычислить усредненные величины и u.
В качестве основы построения областей усреднения (G ) в (42)Ц(43) возьмем куб Cx, y,z с центром в точке (x, y, z) и ребром длиной 2H. Ребра куба направим по осям координат.
Введем регулярную сетку узлов {xi, y, zk } (i, j,k = 0,1,2,...) с поj стоянным шагом h = H / n (n = 1,2,3,...). Сетка такая, что i+1 = i + h, ( = x, y, z).
Ребро 2H определяет размер области окрестности точки (x, y, z), шаг h Ч точность аппроксимации интегралов численными формулами с условием, чтобы не нарушался баланс соотношений (42)Ц(43).
Назовем куб Cx, y,z основным элементом.
Если Cx, y,z S = и (x, y,z) D, то такой основной элемент будем называть внутренним и обозначать Dx, y,z.
Если Cx, y,z S и (x, y,z) D, то такой основной элемент будем называть граничным и обозначать Sx, y,z. Таким образом, Sx, y,z = Cx, y,z (S D).
Все регулярные узлы сетки {xi, y, zk }, расположенные внутри обj ласти D, будем считать в (42)Ц(43) центрами элементов Di, j,k = Dxi, y j,zk и Si, j,k = Sxi, y j,zk. Таким образом, в (42)Ц(43) область G возьмем в виде Gi, j,k = Dxi, y j,zk или Gi, j,k = Sxi,y j,zk. Покроем этими элементами замкнутую область D S. Покрытие обозначим H,n ; очевидно, оно многократно перекрывает все точки покрываемой области. Ясно, что покрытие области D S можно осуществить многими способами, например, ввести регулярную сетку узлов с переменным шагом, криволинейную систему координат, в качестве основных элементов выбрать равнореберные тетраэдры и т.п.
Для всех элементов покрытия H,n запишем уравнения (42)Ц(43).
Получим систему обыкновенных дифференциальных уравнений с независимой переменной t и искомыми усредненными величинами vi, j,k и pm,i, j,k ( m = 1,2,3), отнесенными условно ко всем точкам (xi, y, zk ) D.
j Проиллюстрируем метод на случае, когда определяющим уравнением является закон Гука = (x, y, z) uE + 2(x, y, z), а вектор плотности объемных сил f равен нулю. Здесь тензоры E, Ч единичный и малых деформаций.
Переменную плотность = (x, y, z) аппроксимируем в рамках областей Gi, j,k постоянной i, j,k = (xi, y, zk ) и запишем уравнения (42) - j (43) так:
d 1 n ds, v dg = dt Gi, j,k Gi, i, j,k Gi, j,k Gi, j,k j,k d mm dg = ( n v + 2 nmvm ) ds + dt Gi, j,k Gi Gi, j,k Gi, j,k, j,k 1 1 + n v + 21- nmvm d (m =1,2,3), 1- 2 2 Gi, j,k i , j,k d qr dg = (nqvr + nrvq ) ds + dt Gi, j,k Gi Gi, j,k Gi, j,k, j,k 1 + (1- )(nqvr + nrvq ) d (q = 1, 2; r = 2,3; q r).
Gi, j,k i, j,k Здесь поверхность i, j,k является пересечением поверхности ( ) разрыва упругих параметров (контакта различных горных пород) и элемента Gi, j,k.
Этот элемент разрезается границей на части G1,i, j,k и G2,i, j,k. Часть G1,i, j,k Ч та, в которой располагается центр элемента Gi, j,k ; она имеет внешнюю единичную нормаль n = (n1,n 2,n3) и для нее определены упругие параметры 1, 1. Часть G2,i, j,k имеет упругие параметры 2, 2.
Алгоритмы расчетов для этих частей те же, что и для граничных элементов покрытия. При необходимости по вычисленному тензору напряжений ( ), используя известные формулы, можно найти тензор деформаций ( ) и вектор смещений (u ).
Важно отметить, что неоднородность геологической среды обычно вызывает большие вычислительные сложности при решении сейсмических задач из-за того, что в дифференциальные уравнения движения входят не только функциональные материальные параметры, но и их производные;
на границах разрыва этих параметров требуется разрешение граничных условий; в случае ограниченных твердых тел со сложной геометрией границы или финитных данных Коши необходимо преобразование физической области в расчетную. Компьютерная реализация вычислительных алгоритмов сложна и требует большой ручной работы при решении конкретных задач. Предложенное покрытие с введением граничных элементов освобождает от ручной работы по преобразованию физической области в расчетную, как это делается обычно. В данном случае ручная подготовка вычислительного процесса минимальна, однако колоссально возрастает объем автоматических вычислений, но современные компьютеры могут справиться с этой работой.
В заключении содержатся результаты, полученные в диссертационной работе.
Из известной математической модели распространения нелинейно упругих волн впервые выделена математическая модель волн четко определенного типа - квазипродольных, квазипоперечных и поперечных, в совокупности названных монотипными. Общая нелинейно-упругая волна представляется как результат взаимодействия монотипных волн. Монотипные волны распространяются независимо одна от другой.
Построена теория монотипных волн на основе полученного общего решения нелинейных дифференциальных уравнений, описывающих их динамику, в форме инвариантов Римана, без учета действия внешних сил и с учетом их. Получены волны Римана как частный случай монотипных волн.
В фазовой плоскости получено уравнение фронта ударной волны и все динамические величины на ней в зависимости от формы волны Римана, породившей ударную волну.
Нелинейная упругая среда при прохождении через нее сейсмических волн вызывает локальное перераспределение энергии внутри волн в виде зон увеличения и уменьшения ее плотности. Локальное увеличение плотности энергии в волне вызывает появление сейсмической эмиссии волн и источников колебаний. Сейсмическая эмиссия проявляется в спектрах волн в виде появления кратных частот. Источники колебаний образуются в виде ударных волн с феноменом уменьшения скорости их распространения и появлением лотраженных волн от ударной.
Эти процессы проливают свет на геофизическое понятие лэнергонасыщенности реальной среды, наблюдаемое, но не объясненное.
Построена линейная и нелинейная математические модели распространения сейсмических волн в неупругих диссипативно-дисперсных средах - флюидо-насыщенных резервуарах. Эти модели параметрические и объясняют происхождение эффекта увеличения амплитуд отраженных от резервуара волн - феномен ляркого пятна на сейсмограмме, уменьшение скорости и частоты колебаний, теоретически вплоть до нуля. В случае наличия вязкости среды происходит затухание колебаний и появление запаздывающей волны по отношению к основной. По этим особенностям отраженных от резервуара волн возможно проводить их поиск.
Предложен близкий к реальности подход к построению математических моделей сейсмических процессов в Земле, описываемых интегральными законами сохранения. В основе лежит положение, что все явления изучаются не в точках среды, а в их окрестностях в усредненном виде с использованием интегральных законов сохранения механики сплошной среды, для чего рассматриваемая физическая область покрывается внутренними и граничными элементами с характерным размером (в качестве элементов могут быть выбраны, например, кубы, равнореберные тетраэдры и т.п.). Размер покрывающих областей определяет усреднение свойств среды и выбирается из условий конкретной задачи. Для этих элементов записываются интегральные законы сохранения. Получена система обыкновенных дифференциальных уравнений для элементов покрытия со временем в качестве независимой переменной.
Предложен один из вариантов такого покрытия и приведен пример численной реализации метода для задач распространения сейсмических волн. Метод полностью ориентирован на компьютерную обработку с едиными алгоритмами расчетов для линейных и нелинейных волн, гетерогенных, упругих и вязкоупругих сред и т.п., причем условия контакта различных сред удовлетворяются автоматически, а физическая и расчетная области совпадают.
Основные публикации по теме диссертации 1. Публикации в изданиях, рекомендованных ВАК РФ:
1. Гурьянов В.B. Математическая модель плоских сейсмических ударных волн / В.В. Гурьянов // Вестник Cаратовского государственного технического университета. 2007. №1. Вып. 2. С.7-14.
2. Гурьянов В.B. Особенности распространения сейсмических волн в коллекторах, влияющих на их выявление и дифференциацию. Ч.2 / В.В.
Гурьянов, В.М. Гурьянов, В.Б. Левянт // Геофизика: науч.-техн. журн.
ЕАГО. М., 2003. №4. С.6-10.
3. Гурьянов В.B. Интегральные законы сохранения в прямых задачах сейсмологии / В.В. Гурьянов // Доклады РАН. 2002. Т.385. №1. С.107-109.
4. Гурьянов В.B. Особенности распространения сейсмических волн в коллекторах, влияющих на их выявление и дифференциацию. Ч.1 / В.В.
Гурьянов, В.М. Гурьянов, В.Б. Левянт // Геофизика: науч.-техн. журн.
ЕАГО. М., 2001. №6. С.10-15.
5. Гурьянов В.B. Монотипные плоские нелинейные сейсмические волны / В.В. Гурьянов // Известия АН СССР. Сер. Физика Земли. 1992. №7.
С. 81-88.
6. Гурьянов В.B. Обобщение формулы Кирхгофа / В.В. Гурьянов, В.М. Гурьянов // Известия АН СССР. Сер. Физика Земли. 1992. №3.
С.93-96.
7. Гурьянов В.B. Монотипные плоские нелинейные сейсмические волны / В.В. Гурьянов // Доклады АН СССР. 1991. Т. 319. №1. С. 121Ц124.
8. Гурьянов В.B. Особенности распространения нелинейных сейсмических волн в изотропных средах / В.В. Гурьянов // Доклады АН СССР.
1991. Т. 316. №4. С. 875-879.
9. Гурьянов В.B. Взаимодействие плоских нелинейных сейсмических волн / В.В. Гурьянов // Известия АН СССР. Сер. Физика Земли. 1990. №11.
С. 57Ц70.
2. Другие публикации 10. Гурьянов В.B. Спектры нелинейных сейсмических волн и их особенности / В.В. Гурьянов, М.И. Буслаева // Математика, механика: сб. науч. тр. Вып.7. Саратов: Изд-во СГУ, 2005. С. 165-168.
11. Гурьянов В.B. Математическая модель сейсмических волн во флюидо-насыщенных резервуарах / В.В. Гурьянов, М.А. Антонова // Математика, механика: сб. науч. тр. Вып.6. Саратов: Изд-во СГУ, 2004. С.186189.
12. Гурьянов В.B. Особенности распространения сейсмических волн в коллекторах, влияющие на их выявление и дифференциацию / В.В. Гурьянов, В.М. Гурьянов, В.Б. Левянт // Математические методы в геофизике:
тр. Междунар. конф. Ч. 1. Новосибирск: Изд-во НВМиМГ СО РАН, 2003.
С.93-98.
13. Гурьянов В.B. Интегральные законы сохранения в прямых задачах сейсмологии / В.В. Гурьянов // Математические методы в геофизике:
тр. Междунар. конф. Ч. 1. Новосибирск: Изд-во НВМиМГ СО РАН, 2003.
С.64-68.
14. Гурьянов В.B. Совершенствование глубинных построений с использованием программ TL5 и RAZREZ / В.В. Гурьянов, О.И. Шкуратов, В.М. Гурьянов // Недра Поволжья и Прикаспия: регион. науч.-техн. журнал. Вып.34. Саратов: Изд-во НВНИИГГ, 2003. С.59-64.
15. Гурьянов В.B. Численный метод решения задач распространения сейсмических волн на основе интегральных законов сохранения / В.В.
Гурьянов // Недра Поволжья и Прикаспия: регион. науч.-техн. журнал.
Вып.34. Саратов: Изд-во НВНИИГГ, 2003. С.43-46.
16. Гурьянов В.B. Математическая модель эффективного метода вибрационной сейсморазведки для изучения особенностей распространения сейсмических волн в коллекторах / В.В. Гурьянов, С.И. Михеев, М.В. Живодрова // Недра Поволжья и Прикаспия: регион. науч.-техн. журнал.
Вып.31. Саратов: Изд-во НВНИИГГ, 2002. С.30-34.
17. Гурьянов В.B. Identification of cavernous fractured reservoirs in compact and consolidated rocks by seismic methods / В.В. Гурьянов, В.Б. Левянт // Petropatch: тез. докладов на Междунар. конф. геофизиков в НьюДели (Индия), январь 2001. New-Delhi, 2001.
18. Гурьянов В.B. Ударные волны обобщенного уравнения КарманаФальковича / В.В. Гурьянов, В.М. Гурьянов // Аэродинамика. Ударноволновые процессы: межвуз. сб. науч. тр. Вып.15(18). Саратов: Изд-во Сарат. ун-та, 2001. С.45-52.
19. Гурьянов В.B. Численный метод решения волнового уравнения на основе усреднения дифференциального оператора / В.В. Гурьянов, В.М. Гурьянов // Математика, механика: сб. науч. тр. Вып.2. Саратов: Издво СГУ, 2000. С.165-167.
20. Гурьянов В.B. Плоские упругие волны конечных деформаций в анизотропных средах / В.В. Гурьянов, В.М. Гурьянов // Математика, механика, математическая кибернетика: межвуз. науч. сб. Саратов: Изд-во СГУ, 1999. С.86-88.
21. Гурьянов В.B. Влияние термодинамики на механику плоских волн конечных деформаций / В.В. Гурьянов // Проблемы прочности материалов и конструкций, взаимодействующих с агрессивными средами: межвуз. науч. сб. Саратов: СГТУ, 1999. С.58-64.
22. Гурьянов В.B. Продольные плоские упругие волны конечных деформаций с учетом термодинамики / В.В. Гурьянов; СГУ. Саратов, 1998.
9 с. Деп. в ВИНИТИ 27.02.98, №596-В98.
23. Гурьянов В.B. Дифференциальные уравнения ударных плоских волн / В.В. Гурьянов // Математика, механика и их приложения: тр. науч.- практ. конф. Саратов: Изд-во Сарат. ун-та, 1998. С.63.
24. Гурьянов В.B. Плоские изэнтропические волны конечных деформаций: учеб. пособие для студентов механико-математического факультета. / В.В. Гурьянов, Ю.Д. Каплунов. Саратов: Изд-во Сарат. ун-та, 1997.
45 с.
25. Гурьянов В.B. Приведение системы дифференциальных уравнений нелинейных плоских волн к нормальной характеристической форме / В.В. Гурьянов, В.М. Гурьянов // Аэродинамика. Нелинейные проблемы:
межвуз. сб. науч. тр. Вып.14(17). Саратов. Изд-во Сарат. ун-та, 1997.
С.128-130.
26. Гурьянов В.B. Ударные волны уравнения Кармана-Фальковича / В.В. Гурьянов, В.М. Гурьянов // Аэродинамика. Нелинейные проблемы:
межвуз. сб. науч. тр. Вып.14(17). Саратов: Изд-во Сарат. ун-та, 1997.
С.34-39.
27. Гурьянов В.B. Ударные волны в упругой среде Мурнагана / В.В.
Гурьянов // Фундаментальные и прикладные проблемы механики деформируемых сред и конструкций: сб. науч. тр. Вып.2. Н. Новгород: Изд-во Нижегород. ун-та, 1995. С.79-86.
28. Гурьянов В.B. Монотипные плоские изэнтропические волны конечных деформаций / В.В. Гурьянов // Фундаментальные и прикладные проблемы механики деформируемых сред и конструкций: сб. науч. тр.
Вып.1. Н. Новгород: Изд-во Нижегород. ун-та, 1993. С.149-157.
29. Гурьянов В.B. Алгоритмы решения задач распространения безударных волн в средах Генки: отчет по НИР СГУ. Ч.1 / В.В. Гурьянов, Л.В.
Борисова, О.А. Шитова; СГУ. Саратов, 1994. 31 с. Деп. в ВНТИЦ, инв.
№02.930002392.
30. Гурьянов В.B. Распространение и взаимодействие нелинейных сейсмических волн / В.В. Гурьянов // Динамические задачи механики сплошной среды. Теоретические и прикладные вопросы вибрационного просвечивания Земли. Краснодар: Изд-во Кубанск. ун-та, 1990. Ч.1. С.71.
31. Гурьянов В.B. Распространение и взаимодействие нелинейных сейсмических волн в изотропных средах / В.В. Гурьянов; Институт физики Земли АН СССР. М., 1991. 37 с. Деп. в ВИНИТИ, 16.01.1991, №476-В91.
Авторефераты по всем темам >>
Авторефераты по техническим специальностям