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

• Хохлов Николай Игоревич
• доктор наук
• 2022

Оглавление диссертации доктор наук Хохлов Николай Игоревич

Введение

Глава 1. Определяющие уравнения

1.1 Введение

1.2 Система уравнений линейной теории упругости

1.3 Система уравнений линейной теории упругости: тензорная запись

1.4 Система уравнений линейной теории акустики

1.5 Двухбереговая модель трещины

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

1.7 Запись систем линейной теории упругости в матричном виде

1.8 Запись систем линейной теории акустики в матричном виде

1.9 Выбор системы координат

1.10 Запись систем линейной теории упругости и акустики в выбранной системе координат

1.11 Основной алгоритм

1.12 Диагонализация матриц уравнения упругости

1.13 Диагонализация матриц уравнения акустики

1.13.1 Граничные и контактные условия

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

2.1 Введение

2.2 Уравнение переноса с постоянными коэффициентами

2.3 Численный схемы для уравнения переноса

2.4 Компактные сеточно-характеристические схемы

2.4.1 Продолженная система уравнений

2.4.2 Интерполяционные полиномы

2.4.3 Исследование поведения интерполяционных полиномов

2.4.4 Гибридные схемы

2.4.5 Сеточно-характеристическая монотонизация

2.4.6 Гибридная схема третьего порядка точности

2.5 Тестирование схем

2.5.1 Импульс формы прямоугольного треугольника

2.5.2 Сеточная сходимость на равномерной сетке

2.5.3 Сеточная сходимость на сетке с неравномерным шагом

2.5.4 Сеточная сходимость в случае разрывного решения

2.5.5 Численное исследование консервативности полученных

схем

2.5.6 Численное исследование монотонности схем

2.5.7 Выводы

2.6 Сеточно-характеристические схемы для систем уравнений с переменными коэффициентами

2.7 Системы с переменными коэфициентами

2.8 Разностные схемы

2.9 Тестирование метода

2.9.1 Двуслойная среда с постоянным импедансом

2.9.2 Двуслойная среда с переменным импедансом

2.9.3 Трехслояная среда с переменным импедансом

2.9.4 Многослойная среда с переменным импедансом

2.10 Сеточная сходимость

2.11 Линейное уравнение переноса с кусочно-постоянной скоростью

2.12 Построение сеточно-характеристической схемы повышенного

порядка для кусочно-постоянной скорости переноса

2.12.1 Тестирование

2.13 Схема для уравнение акустики с переменными коэффициентами

2.14 Тестирование нового алгоритма для уравнения акустики

2.14.1 Выводы

2.15 Моделирование трещиноватых неоднородностей на структурных сетках

2.16 Подход 1: расположение трещины в ячейках расчетной сетки

2.16.1 Двумерный случай

2.16.2 Трехмерный случай

2.17 Подход 2: дублирование узлов в области трещины

2.18 Подход 3: использование наложенных или Химерных сеток

2.19 Использование наложенных сеток для выделения топографии

2.20 Верификация численных методов и алгоритмов

Глава 3. Программный комплекс для численного

моделирования волновых возмущений

3.1 Оптимизация работы сеточно-характеристического метода на CPU

3.1.1 Введение

3.1.2 Применение тайлинга к сеточно-характеристическому методу

3.1.3 Использование потоковых инструкций процессора и распараллеливание в общей памяти

3.2 Распараллеливание в системах с распределенной памятью

3.2.1 Введение

3.2.2 Декомпозиция расчетных сеток

3.2.3 Тестирование на одноблочной сетке

3.2.4 Тестирование на многоблочных сетках

3.3 Распараллеливание сеточно-характеристического метода на GPU

3.3.1 Введение

3.3.2 Особенности реализации сеточно-характеристического метода

3.3.3 Распараллеливание на нескольких графических процессорах

3.4 Выводы

Глава 4. Численное решение задач рапространения динамических возмущений в гетерогенных

геологических средах

4.1 Численное моделирование поверхностных сейсмических волн

4.2 Численный расчет сейсмического отклика от кластера трещин

4.3 Численное исследование анизотропии отклика сейсмического сигнала от трещиноватого слоя

4.4 Численное исследование влияние трещиноватых включений на примере модели Overthrust

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

4.6 Численное решение задач рапространения динамических возмущений с явным выделением криволинейных границ

Глава 5. Численное решение ряда задач распространение

волновых возмущений

5.1 Моделирования распространения волновых возмущений при землетрясениях

5.2 Решение задач глобальной сейсмики

5.3 Численное решение задач динамических разрушений

5.4 Другие приложения программного комплекса

Заключение

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

Введение диссертации (часть автореферата) на тему «Моделирование распространения динамических волновых возмущений в гетерогенных средах на высокопроизводительных вычислительных системах»

Введение

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

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

Математическое моделирование проводится в различных геологических средах, в том числе, в слоистых средах и в средах с наличием различных неодно-родностей (например, трещин или каверн). Задачи такого рода представляются очень ресурсоемкими с точки зрения вычислительных ресурсов. Область вычисления, как правило, представляет собой сейсмический куб с длинной ребра от 1 км до 10 км. В тоже время, неоднородности могут быть размером в несколько метров.

При моделировании задач сейсмостойкости также приходится сталкиваться с множеством неоднородностей. В современном мире все больше внимания уделяется оценке и обеспечению безопасности населения. Активно развиваются методы имитационного моделирования, позволяющие провести всесторонний анализ рисков и выделить наиболее значимые из них, на снижении которых необходимо сконцентрироваться. Очевидно, что любое моделирование основывается на некотором приближении (модели) реальных физических процессов. Если ранее модели, в основном, были простыми, обеспечивающими решение определяющих уравнений аналитическими методами, то с развитием высокопроизводительных вычислительных систем произошло их качественное изменение. Все больше усложняется постановка задач, включающая описание движения отдельных атомов и молекул, с целью моделирования таких макроскопических явлений, как течение воды в реках, сход оползневых тел и т.д. Одним из процессов, притягивающих внимание научного сообщества, является естественная сейсмичность, выраженная в происхождении землетрясений. Только лишь за последние два года во всем мире произошло более 10 крупных землетрясений (Китай, Турция, Япония), причинивших огромные материальные убытки и приведших к гибели тысяч человек. Таким образом, актуальность исследования процесса сейсмической активности с целью ее прогноза и быстрой оценки возможных разрушений не вызывает сомнений. Для этого по всей планете расположено множество специальных сооружений (сейсмостанций), в которых ведется непрерывная регистрация движений поверхности Земли. Наиболее плотные сети измерений расположены в Калифорнии и Японии. Как правило, измеряется зависимость вектора скорости или ускорения (или части их проекций) от времени в точке расположения сейсмического датчика (сейсмометра). При инициации землетрясения из его очага, расположенного обычно на глубинах 5-30 км, начинают распространяться объемные упругие волны, которые, достигая поверхности, приводят к разрушению строений. На настоящий

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

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

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

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

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

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

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

тивных сред и при использовании детального описания геосреды отличаются, не только количественно, но и качественно.

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

Так, в работах [6; 7] рассматривается распространение упругих волн в гетерогенных средах. Распространение волн Рэлея в гетерогенных средах численно моделируется в [8], однако в работе рассмотрен только двумерный случай. В работах [7; 9] проводится численное моделирование процессов распространения волн с использованием метода конечных элементов, в том числе с повышенным порядком точности.

Наибольшее распространение при решении таких задач получили метод конечных элементов (МКЭ), характеристические и сеточно-характеристические методы и конечно-разностные методы.

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

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

Одними из наиболее широко используемых вычислительных методов для решения систем уравнений гиперболического типа механики сплошных сред являются сеточно-характеристические методы, подробное описание и обзор которых можно найти в [10]. Значительно возросшая за последние десятилетия производительность вычислительных систем позволили увеличить размерность

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

Этот метод, строго говоря, не является прямым методом характеристик, а использует условия совместности в фиксированной форме для фиксированной расчетной сетки. Анализ характеристических свойств уравнений механики сплошных сред используется для формулировки обратного метода характеристик на основе известных понятий и терминов теории разностных схем. Эти идеи были опубликованы в работе [10], в которых было показано, что используя определенные комбинации условий совместности для многомерной системы уравнений в частных производных гиперболического типа общего вида, можно построить численные схемы обратного метода характеристик первого и второго порядков точности. Причем они являются естественным обобщением метода [12] и обладают высокой эффективностью при решении конкретных динамических задач, поскольку позволяют корректно строить вычислительные алгоритмы на границах области интегрирования и контактных границах, в определенной степени учитывать физику задачи (распространение возмущений вдоль характеристик). Предложенные схемы полностью описываются в терминах метода сеток и получили название сеточно-характеристических методов, что позволило исследовать соответствующие схемы на устойчивость, аппроксимацию, монотонность. Анализ метода для простейшего уравнения переноса показал, что схема первого порядка точности (сеточно-характери-стического метода) имеет наименьшую аппроксимационную вязкость среди явных разностных схем для простейшего шаблона [13]. В этой же работе был разработан аппарат исследования разностных схем — их анализ в пространстве неопределенных коэффициентов. Это позволило существенно расширить возможности сеточно-характеристического метода и рассмотреть разностные схемы для произвольного наперед заданного сеточного шаблона, в том числе гибридные методы высокого порядка точности [5; 14]. При помощи сеточно-

характеристического метода и его гибридных модернизаций с использованием прямоугольных и треугольных расчетных сеток было получено численное решение широкого класса задач динамики деформируемого твердого тела: [5; 15—18]. Для численного решения задач разведывательной сейсмологии се-точно-характеристические методы, с использованием треугольных расчетных сеток, разработаны в [4; 19; 20]. Особенностью этих работ является корректное выделение многочисленных контактных или свободных границ (для флюидона-сыщенных и пустых трещин и каверн) в геологических средах для исследования их структуры, что позволяет наиболее корректно строить расчетные алгоритмы в этих областях используя сеточно-характеристические методы. Эти подходы по этой же причине оказываются эффективными и для решения задач скважин-ной сейсмологии, в том числе, и при наличии трещин вблизи скважины. В этих работах рассматривались динамические процессы, как в кристаллических, так и в осадочных породах. Были получены фронты волн-откликов на сейсмосиг-нал, фронты волн в самих кластерах с большим количеством неоднородностей, дифракционные картины при отражении упругих волн, как от кластеров, так и от изолированных флюидонасыщенных и пустых неоднородностей. Показаны различия в физических картинах при использовании эффективных моделей с выделением неоднородностей. Успешное решение этого класса задач в немалой степени обусловлено не только возможностью построения корректных алгоритмов расчета функций в граничных точках и в точках, принадлежащих контактным границам, но и важнейшим свойством монотонности используемых гибридных сеточно-характеристических схем.

Характерной чертой уравнений гиперболического типа является возможность появления разрывных решений, в нелинейном случае — даже при гладких краевых условиях. Важным классом для таких уравнений являются монотонные разностные схемы (схемы с положительной аппроксимацией или мажорантные схемы), переводящие монотонный профиль численного решения на нижнем временном слое в монотонный профиль на верхнем временном слое. Наиболее заметно различия в численных решениях при расчетах по монотонности и немонотонным схемам проявляются на грубых пространственных сетках и при больших временах расчета. Существуют различные критерии монотонности разностных схем и способы регуляризации (монотонизации) разрывных численных решений, начало которым положено введением искусственной вязкости [21; 22], операторов "сглаживания" (фильтрацией высоких

частот) высокочастотных осцилляций ([22; 23] и др.). Для уравнений и систем уравнений гиперболического типа понятие о таких схемах было введено Фри-дрихсом [24]. К монотонным по Фридрихсу схемам принадлежат известные методы ([25; 26] и др.).

В работах [27] для монотонизации вводится понятие TVD-схем (схем с изменением полной вариации) и вводится понятие "ограничителя"или "лимит-ра"для переключения расчетных схем. В дальнейшем этим же автором были предложены ENO-схемы (essentially non-oscillatory schemes) и UNO (uniformly non-oscillatory schemes).

При моделировании сейсмических явлений важным аспектом является выделение границ раздела контактных поверхностей. Данная задача может рассматриваться как задача погруженного интерфейла (immersed interface) и является проблемой, возникающей в большом количестве областей. Для решения такого рода задач уже предложен ряд численных методов. Постановки можно раздлелить на два типа, когда граница раздела сред согласована с сеткой (fitted meshes), т.е. скачок среды приходится на ячейку или узел, и может быть не согласована с сеткой (non-fitted meshes). Следует заметить, что метод погруженного интерфейса достаточно полно представлен для задач параболического и эллиптического типа [28—32] и др. Для задач гиперболического типа он представлен не так широко. Для случая, когда граница раздела сред согласована с конечно-разностной сеткой существует ряд подходов для повышения порядка в области скачка параметров. Так в работе [33] рассматривается модификация конечно-объемного метода типа Годунова на примере уравнений акустики и линейной теории упругости. В работе приводится точное решение задачи распада-разрыва на границе интерфейса и метод, учитывающий это. Рассматриваются TVD методы, обладающие монотонным поведением. Однако, как описано в самой работе, порядок схемы в области скачка коэффициентов падает. Чем больше скачок параметров, тем больше ошибка. В работе [34] рассматривается подход повышения порядка точности в области скачка параметров для конечно-разностных методов на сдвинутых сетках типа FDTD [35]. Рассматривается подход усреднения (гомогенизации) коэффициентов конечно-разностной схемы для получения второго порядка на границе раздела сред. В простейшей реализации методы данного класса дают первый порядок точности в области скачка параметров. Также для решения такого рода задач на неструктурных сетках широко применяется метод разрывного Галеркина (discontinuous

Galerkin method) [36]. В работе [36] рассматривается данный метод на двумерных неструктурных сетках и показан подход, позволяющий строить методы произвольного порядка. В частности данный метод в том числе позволяет осуществлять расчет границы раздела сред и между различными физическими средами, типа акустическая-упругая среда [37].

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

В работе [41] авторами предложен так называемый явный упрощенный метод погруженного интерфейса (explicit simplified interface method - ESIM). Его отличие от классического IIM метода заключается в том, что на всей сетки используется один и тот же метод повышенного порядка точности и только в точках вблизи интерфейса применяются некоторые корректировки, которые не зависят от используемого метода. В работе [41] авторы рассматривают одномерный случай. Далее авторы рассмотрели двумерный случай [42], где рассматривают расширение данного подхода на двумерный случай. Метод позволяет производить расчет границы произвольной формы на структурной прямоугольной сетке. В дальнейшем авторы применяли данный подход и для других систем гиперболических уравнений в частных производных [43]. Метод выглядит достаточно универсальным, однако требует препроцессинга для вычисления матриц коэффициентов на границе и хранения дополнительных данных для границы. Для задач, где параметры среды меняются плавно или достаточно много границ раздела сред, применение такого рода методов может быть затруднительно.

В работе [44] предложен так называемый Correction Function Method (CFM). Авторы привели метод четвертого порядка точности для двумерного уравнения акустики. В подходе, предложенным авторами, к узлам вблизи границы раздела сред прибавляется компенсирующая добавка после каждого шага, чтобы учесть скачок параметров. Потенциально метод может быть расширен

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

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

Идея предложенного метода построена на сеточно-характеристическом методе или grid-characteristic method (GCM) [45] и на кусочно-полиномиальной интерполяции [46]. Сеточно-характеристический метод применяется для достаточно широкого круга задач. В частности он применяется для расчета различных геофизических задач [Stognii2021; Golubev2020], в том числе с наличием трещиноватых неоднородностей [49], задач движения поездов [50], задач сейсмики и сейсмостойкости [51] и др.

Одним из направлений для построения численных методов повышенного порядка точности является использование продолженных систем уравнений [52]. Схемы, построенные для таких систем, используют дифференциальные следствия исходных уравнений [53—56], что позволяет использовать для численной аппроксимации исходного уравнения переноса шаблон с минимальным набором точек, одновременно сохраняя повышенный порядок точности. Схемы такого типа принято называть компактными [57]. При построении компактных разностных схем для систем гиперболических уравнений в основном используются трех точечные [56; 57] и двухточечные шаблоны [53; 58] по координате. Максимальный порядок схем для систем гиперболических уравнений на двухточечном шаблоне, из представленных в литературе, - четвертый [53; 58]. В работе [58] используются интегральные следствия исходных дифференциальных уравнений, что позволило авторам предложить монотонную схему первого

порядка по времени и монотонизированную схему третьего порядка по времени, причем, обе схемы имеют четвертый порядок по координате. В обзорной работе [53] четвертый порядок схемы также достигается за счет использования искомой системы и ее дифференциального следствия. Разностная схема для исходной системы использует значение из дифференциального следствия на двух временных слоях. Схема монотонизированна и использует сеточно-характери-стический критерий монотонности [59].

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

Список литературы диссертационного исследования доктор наук Хохлов Николай Игоревич, 2022 год

— - и

р • ^ Л (пг • и) I + ц (пг ® и + и ® щ)

для систем уравнений линейной теории упругости и в виде

(1.55)

А,

— и

-жц

Л (пг • и)

(1.56)

для системы уравнений акустики.

1.11 Основной алгоритм

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

Рассмотрим независимо системы вида

щ + Агиг — 0, (1.57)

где г € 1, 2, 3 в зависимости от размерности системы (1.29). Рассмотрим в начале случай системы с постоянными коэффициентами. Для применения сеточно-характеристического метода рассмотрим диагонализацию матрицы А — А^ (в дальнейшем, где понятно из контекста, индекс г будем опускать)

А — П-1ЛП, (1.58)

здесь матрица Л - диагональная матрица размером N х N, состоящая из собственных чисел Лj матрицы А, матрица ^ - матрица, состоящая из собственных векторов матрицы А. В работе рассматриваются только системы уравнений гиперболического типа, в этом случае матрица А имеет полный набор действительных собственных чисел и N линейно независимых собственных векторов (соответственно матрица ^ имеет обратную).

В классических работах [77], описывающих сеточно-характеристический

метод, также вводится матрица X вида

х = , (1.59)

здесь г) - ]-й столбец матрицы П, а г^ - ]-я строка матрицы П-1. Тогда, из-за характеристических свойств гиперболических уравнений в частных производных (1.57) значение вектора неизвестных *(£, х) в момент времени £ + т можно выразить как

*(£ + т, х) = Е Xju(t ,х — Л^-т), (1.60)

здесь т - шаг интегрирования по времени. Данная формула (1.60) может быть использована для поиска решения *(£ + т, х) в момент времени £ + т исходя из начальных условий *(£, х). Легко заметить, что

Е х = I. (1.61)

Предположим, что матрица А имеет 3+ положительных собственных чисел, 3— отрицательных собственных чисел и нулевых собственных чисел соответственно. Тогда сумма матриц Xj, отвечающих за нулевые собственные значения может быть получена как

Е х = I — Е х, — Е Xj. (1.62)

Исходя из (1.60) и (1.62) можно записать следующую форму сеточно-характеричтисеского метода [45]

*(£ + т, х) = *(£, х)+ ^^ Xj (*(£,х — Л^-т) — *(£, х)). (1.63)

Запись (1.63) выгоднее (1.60) тем, что не требует переноса компонент вектора и, соответствующих нулевым собственным значениям Л.

При практической реализации записи вида (1.63) или (1.60) редко когда используются. Поскольку матрицы X имеют плотную структуру и вычисление их действия на вектор требует относительно много вычислительных операций.

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

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

ш (г,х) — Ш(Ь,х). (1.64)

Переменные ш иногда называют инвариантами Римана или характеристическими переменными.

За вторе действие необходимо найти вектор ш(р + т,х)

+ Т, х) — х) + X — Л^т) — Ш^, х)) , (1.65)

что эквивалентно решению системы

ш г + Леи х — 0, (1.66)

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

шг,х + Лгшг,ж — 0. (1.67)

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

ф + т,х) — (Ь + т,ж). (1.68)

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

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

1.12 Диагонализация матриц уравнения упругости

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

[33] приведен вид матриц для двумерного уравнения упругости. Одно из наиболее полных описаний вывода данных матриц в тензорном виде приводится в работе [78]. Рассмотрим вид данных матриц. Для трехмерного случая и произвольных осей координат он довольно громоздок и удобно записывать его в тензорном виде и в виде действия оператора на вектор. Тогда для уравнения упругости действие матрицы ^ будет иметь вид

П

и а

Н • U - ¿Р ^00 : а П • и + ¿Р ^0 : а

• и - ^01 : а

cs p

П1 • V + ¿p N01 : а П2 • и - ^ ^02 : а

CS p

П2 • v + ^ N02 : а N12 : а

(N11 - N22) : а

N11 + N22 - ^N00) : а Здесь : - двойное скалярное произведение (свертка), a : b = ^a,ijbij. Обратная матрица к (1.69) будет

(1.69)

У

^-1Ш = -

(Ш1 + Ш2) п + (шэ + Ш4) Щ + (Шб + Шб) П2

Р (Ш2 - Ш1) ((ср - Сэ) ^00 + Сэ1) + 2csp (Ш4 - Шэ) ^01 +

+2csp (Шб - Шб) N02 + 4шт^12 + Шв (Nn - N22) + Ш9 (I - N00)

(1.70)

где сэ =

Л+2ц

Со.

Вид действия операторов (1.69) и (1.70) в тензорной записи одинаков для двумерного и трехмерного случаев. Также он не зависит от осей координат. Покомпонентная запись будет иметь вид [51]

ш = fi • и т (N00 : а),

cpp

Шэ,4 = П1 • и т — (N01 : а), csp

Шб,6 = П2 • и т — (N02 : а),

Csp

(1.71)

(1.72)

(1.73)

Л

ш7 = N12 : а, (1.74)

ш8 = (N11 - N22) : а, (1.75)

Шд = ^N11 + N22 - у3N0^ : а. (1.76) и обратная к ней

и = 1 ((Ш1 + Ш2) п + (шз + Ш4) П1 + (Ш5 + Шб) П2). (1.77)

1

(р (Ш2 - Ш1) ((ср - Сз) Дю + Сз^ + 2С5р (Ш4 - Шз) Д)1+2Сзр (Шб - Ш5) N02 +

(1.78)

Данная запись может быть удобнее при программной реализации. В любом случае стоит смотреть на язык программирования и его возможности, чтобы реализовать данные операторы наиболее эффективно. Так в языке программирования С++ есть перегрузка операторов и реализация в виде (1.69) и (1.70) может быть предпочтительнее и проще. Однако покомпонентная реализация имеет больше возможностей для оптимизации. В любом случае надо ориентироваться на конкретную реализацию. Более того, для случая когда система координат является ортогональной и расчетная сетка не криволинейная, выражения могут быть существенно упрощены. Тензоры вида N11 и др. будут достаточно разрежены и операции с ними можно существенно оптимизировать. Для задач, которые рассчитываются на больших некриволинейных сетках такая оптимизация может дать существенный выигрыш. Покомпонентная запись в таком случае будет иметь вид (для оси ж1)

Ш1,2 = и + — аш (1.79)

р р

Шз,4 = и + — а12, (1.80)

С,р 1

Ш5,6 = из + — а1з, (1.81)

С,р

Ш7 = а2з,

(1.82)

Ш8 = а22 - азз, (1.83)

2 з

Шд = а22 + азз--аш (1.84)

р

и обратное преобразование

и = 1 (Ш1 + Ш2), (1.85)

и = 1(шз + Ш4) , (1.86)

из = 2(Ш5 + Шб), (1.87)

ац = 2рСр (Ш2 - Ш1), (1.88)

а22 = 2 (рСз (Ш2 - Ш1) + Ш8 + Шд) , (1.89)

азз = 1 (рСз (Ш2 - Ш1) - Ш8 + Шд) , (1.90)

а12 = 2Сар (Ш4 - Шз), (1.91)

а1з = 2С*р (Шб - Ш5), (1.92)

а2з = Ш7. (1.93)

Для других осей (координат) запись будет выглядеть аналогично. Матрица собственных значений Л будет иметь вид

Л = (Над {Цср, -Цср, Цс8, -Цс8, Цс8, -Цс8, 0,0,0} . (1.94)

Здесь Ц определяется из (1.53). Для случая, когда уравнения рассматриваются не в криволинейной системе координат, выражение для Л упрощается, и данная матрица перестает зависеть от оси

Л = (Над {ср, -ср, с8, -св, с8, -св, 0,0,0} . (1.95)

1.13 Диагонализация матриц уравнения акустики

Аналогичным образом может быть диагонализирована система уравнений акустики (1.56). Уравнения акустики возникают, если тензор напряжений имеет гидростатическую форму (шаровую). Спектральное разложение для уравнений акустики легко получается из формул для уравнения упругости, если произвести подстановку а = —р1, Ц = 0 (и вычеркнуть строки, теряющие смысл) [78].

Так преобразование (1.64) будет

ш1 2 = п • и ± —, (1.96)

ср

ш3 = п1 • и, (1.97)

Ш4 = п2 • и, (1.98)

и обратное к нему (1.68)

и = 1 ((ш1 + Ш2) п + Ш3Я1 + Ш4Я2), (1.99)

2

р =2ср (Ш1 - Ш2). (1.100)

И для случая не криволинейной сетки (для одной из осей, на примере ж1)

Ш1, 2 = и1 ± —, (1.101)

ср

Шз = и2, (1.102)

И обратное

Ш4 = из. (1.103)

и = 1(Ш1 + Ш2), (1.104)

и = Шз, (1.105)

из = Ш4, (1.106)

Р=2ср (Ш1 - Ш2). (1.107) Матрицы собственных чисел будут

Л = ( ад {Цс, -Цс,0,0} , (1.108)

и

Л =( ад {с, - с,0,0} , (1.109) соответственно для произвольных осей и ортономированного базиса.

1.13.1 Граничные и контактные условия

Рассмотрим далее контактные и граничные условия, используемые в работе. Подробно они описаны в работах [45; 51], здесь приведем лишь основные формулы.

Условие свободной границы для уравнения упругости можно записать в

виде

а •т = 0. (1.110)

Здесь т - единичный вектор нормали к поверхности. Контактное условие полного слипания между границами раздела поверхностей:

иь = -ид, (1.111)

аь • т = ад •т. (1.112)

Здесь и далее индексы Ь и Я отвечают за материал слева и справа от границы раздела поверхностей (контактной границы), т отвечает за внешнюю нормаль к левой поверхности. Условие на корректировку контактного условия полного слипания используется через граничное условие заданной скорости. Скорость может и может быть найдена по формуле

и = L 1 rR { pL (cLp - cLs) (m • vlin ) m + pLcLsvLIN+

p LcL + pR CR

^R { ^R ^RR\ f\ i ^R^tR i /_.R

+ p R (cR - cR) (m • URn) mm + pRcRuRv + (ctRn - арж) • m -

Р - о£ + Ш - • < + ^^ + (а™ - а1ш) • ™)М>.

р Ср + р Ср

/

Затем корректировка для левой стороный модет быть выполнена

= и, (1.114)

ар+1 = а^ + pL {((и - -uLn) • ») ((cpL ~ 2ср - cp) (m ® m) + cLl) +

+cp (m ® (v - vpN) + (и - upN) ®m)} ,(1.115)

и для правой

uR+1 = U, (1.116)

аR+l = аRRN + pR {((u^n - и) • m) ((cR - 2cR - cR) (m ®m) + cRl) + +cR (m ® (uRN - и) + (uRN - и) ® m) + cR (m ® (uRN - u) + (uRN - и) ® m¡Ц.117)

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

/J f f RIGHT

п+\

0)i// \ \чО)2 \

1 /Юз,5/ Ш4Д 0)3.5 \°4.б\

Ю7,8.9 X j

а а

Рисунок 1.3 — Входящие, нулевые и выходящие характеристики для уравнения

упругости.

Аналогично для уравнений акустики, условие свободной границы

р = 0. (1.118)

Контактные условия можно записать в виде

и1 •т = ид •т, (1.119)

/ = рк. (1.120)

Рассмотрим пример корректора контактного условия полного слипания. Значение Р на границе находится по формуле

р = СЬрЬр?м + СЕрКрЬш + сУсУ1' ((Ц^ - ЦДу) • 77/)

СЬрЬ +СДрД . '

Затем производится корректировка в левом материале

и^+1 = Ц« + Щ-^т, (1.122)

Й+1 = Р, (1.123)

и в правом

= ^ + (1.124)

= Р. (1.125)

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

а • т + рт = 0, (1.126)

иА • т = Vе • т (1.127)

В формулах и далее по тексту индексы А и Е относятся к акустической и упругой среде соответственно, т внешняя нормаль к упругой среде. Значение Р можно найти как

= Ср рЕсрА (и^ - -иАУ) • т + Ср рЕр1Н - срА (т • аш • т) (1128) Р = ср рЕ + срА . ( 8)

Затем корректировка для акустической среды проииходит по формуле

иА+1 = НА« + Р-АЕл, (1.129)

р

Рп+1 = Р

(1.130)

и для упругой

ъ = аш • т + Рт,

(1.131)

иЕ+1 =

^ъ + те(^ - (ъ ^)т,

рЕс: рЕ \с: Ср

(1.132)

ап+1 = аш - (ъ 0 т + т 0 ъ) + (ъ • т) ( — I - (1 + — ) (т 0 т)) .(1.133)

р р

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

2.1 Введение

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

2.2 Уравнение переноса с постоянными коэффициентами

Как уже упоминалось ранее, для решения одномерных уравнений переноса, которым удовлетворяют характеристические переменные, можно использовать различные методы. Рассмотрим простейшее однородное линейное уравнение переноса для функции и = и(х,Ь)

U + Лих = 0.

(2.1)

Будем полагать, что Л = const > 0, поскольку для отрицательных Л все построения аналогичны и могут быть осуществлены заменой Л на —Л и сеточного шаблона на симметричный по х относительно рассчитываемой точки (tn,xm).

m-2

m-1

m

i i

m+1

-át

n+1

n

xm-2 xm-1 xm xm+1

Рисунок 2.1

Базовый шаблон, показывающий принцип построения сеточно-характеристических схем

На рисунке 2.1 приведен базовый шаблон, показывающий принцип построения сеточно-характеристических схем. В основе сеточно-характеристического метода лежит принцип переноса искомых функций вдоль характеристик. Так, чтобы решить уравнение переноса (2.1), из узла т, в котором требуется решить уравнение на временном шаге п + 1 опускается характеристика на временной слой п [19]. Из точки пересечения опущенной характеристики с временным слоем п переносится соответствующее значение искомой функции и в точку ип+1:

ип+1(хт) = ип(хт — Лт). (2.2)

Здесь т - шаг интегрирования по времени. Также введем параметры дискретизации, а именно пространственный шаг по координате И, и число Куранта а = Лт/И.

Не ограничивая общности можно считать хт = 0, тогда можно записать более просто ип+1 = ип(—Лт).

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

Далее рассмотрим ряд схем, применявшихся в данной работе.

2.3 Численный схемы для уравнения переноса

В работе [53] приводится достаточно подробное описание принципа построения сеточно-характеристических схем. В ней приводится алгоритм построения монотонных схем для одномерного уравнения переноса, систем гиперболических уравнений. Рассматриваются вопросы консервативности.

В данной работе рассматриваются только двухслойные схемы по времени. Простейшей и классической схемой для уравнения переноса можно считать схему первого порядка точности типа Куранта-Изаксона-Риса (КИР, СЖ, данные обозначения будут использоваться далее по тексту) [52]

и ит \{ит ит-1).

(2.3)

Тест на качественное поведение разностных схем проводился для начальных данных в виде импульса сложной формы [79—81]:

ехр(- 1п2(ж + 0.7)2/0.0009), -0.8 < ж < -0.6,

м(0, х) = <

1 -|10ж - 1|,

(1 - 100(ж - 0.5)2)1/2,

0,

-0.4 < ж < -0.2, 0.0 < ж < 0.2, 0.4 < ж < 0.6,

в противном случае.

(2.4)

Данное начальное условие состоит из разрывного прямоугольного импульса, треугольного импульса, импульса формы Гауссово распределения и половины эллипса. В дальнейших тестах используются следующие параметры расчета, если не оговорено иное. А именно, число узлов сетки 400, скорость переноса Л = 1, число Куранта а = 0.2, время, до которого ведется расчет £ = 2. Периодические граничные условия.

Рисунок 2.2 — Решение, полученное схемой ОШ (2.55) для импульса (2.4).

Также рассмотрим классические схемы второго порядка аппроксимации, такие как Лакса-Вендорффа

и,

п+1

а

= ит 2((^т+1 ит-1) а(ит+1 2ит + ит-1)),

(2.5)

и явная схема Бима-Уорминга, являющаяся измененной схемой Мак-Кор-

мака

1

+1 ^ ^

ит = ит — 2(3ит — 4ит-1 + ит-2) + у(Ит — 2ит-1 + ит-2))' (2'6)

Еще одна достаточно интересная схема из классических схем - схема

Фромма, она может быть получена путем полусуммы схем (2.6) и (2.5)

2

+1 ^ ^

ит = ит — 4(ит+1 + 3ит — 5ит-1 +4ит-2) + (ит+1 — ит-ит-1 + ит-2))' (2.7)

2.4 Компактные сеточно-характеристические схемы

В данной главе рассматриваются компактные схемы второго-третьего порядка точности построенные с использованием интерполяционных полиномов. Решение системы аппроксимируется с помощью набора интерполяцонных полиномов различного порядка, затем, для обеспечения монотонного поведения решения, схема гибридизируется [82; 83], путем выбора того или иного полинома, в зависимости от характера решения. Используется гибридизация на основе сеточно-характеристического критерия [59] и критерия, основанного на поиске экстремума интерполяционного полинома на выбранном интервале, который предложен в данной работе. Все схемы построены на минимальном двухточечном, по координате, шаблоне, поэтому их можно отнести к классу бикомпактных [58] разностных схем. В результате подобного построения разностной схемы удалось добиться монотонного поведения разностной схемы и слабого "размытия"фронта разрывного решения со временем. В данной работе предложены новые сеточно-характеристические схемы, работы, где они рассмотрены [84; 85].

2.4.1 Продолженная система уравнений

Помимо уравнения (2.1) также рассмотрим его дифференциальное следствие. Введем в дополнение к и(Ь, х) новую искомую функцию у(Ь, х) = их(Ь, х) и, продифференцировав уравнение (2.1) по х, получим для у(Ь, х) аналогичное

Рисунок 2.3 — Решение, полученное схемами (2.5), (2.6) и (2.7) для импульса

(2.4).

уравнение переноса:

Уг + ЛУх = 0.

(2.8)

2.4.2 Интерполяционные полиномы

т — 1 т

Рисунок 2.4 — Шаблон разностной схемы

Построение продолженных разностных схем будем рассматривать на двухточечном шаблоне (рис. 2.4) с пространственным шагом Н и временным шагом т:

{Ь ,хт-\)) ,жто),{Ь ,жто). (2.9)

Будем использовать систему координат, в которой точка (£п,жто) имеет координаты {0,0), а точка {1п,хт—\) - {0, —Н) соответственно.

Рассмотрим на интервале {—Н,0) интерполяционные полиномы /{ж), которые аппроксимирует функцию и{х) на заданном интервале. В дальнейшем индекс по времени будем опускать, где это возможно. Решение уравнения (2.1) на временном шаге п + 1 может быть найдено как

<+ = / {—Лт), а решение продолженной системы будет:

= Л—Лт).

Полином максимальной степени, который можно построить на данном шаблоне - это полином третьей степени. Представим его в виде:

= а3х3 + Ь3х2 + с3х + ¿3 (2.10)

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

Г3 {х) = 3азх2 + 2Ьзх + сз. (2.11)

Коэффициенты данного полинома могут быть найдены из условий:

Д-^) = <_!, /(0) = ипт,

(2.12)

f(-h)= УПт_i, /'(0) = УПт.

Исходя из (2.12) находим:

п _ vm+vm — 1 __и^— 1)

«3 = р К3 ,

°3 =-h---h2-,

Сз = Ут,

(1з ит.

Разностная схема, использующая полином (2.10), по всей видимости, впервые была предложена в работе [86], в ней данная схема описана как для линейных, так и для квазилинейных систем гиперболических уравнений. В оригинальной статье данная схема названа CIP (Cubic-Interpolated Pseudo-particle); в дальнейшем мы будем придерживаться этого названия. В последующих работах схема CIP была распространена на многомерный случай, а также предложены некоторые усовершенствования схемы [87]. Во всех работах схема CIP используется либо без ограничителей, либо с ограничителем, построенным на расширенном шаблоне, использующем также точку (tn,xm+\) [86]. В данной работе предложен способ добиться от схемы монотонного поведения, путем понижения порядка интерполяционного полинома (2.10) в области разрывных решений, используя при этом компактный шаблон (2.9).

Результат работы схемы CIP для импульса сложной формы представлен на рис. 2.5б.

Наряду с полиномом (2.10), рассмотрим полиномы второго порядка точности, обозначив их, как F2 (х) и F2r (х):

F21 (х) = a,2ix2 + b2ix + С21, (2.13)

F2r (х) = а,2Гх2 + Ь2Гх + С2г. (2.14)

Полином F2i (х) построен с использованием точек ит—\,ит, vm—1, а F2r (х) -с использованием точек ит—1,ит, vm. Коэффициенты полиномов (2.13) и (2.14) можно найти из условий:

F21 (0) = ит, F21 (—h) = um—i,F2i(—h) = vm—1, F2r(0) = ит, F2r(—h) = ит—1, F!lr(0) = vm.

Они имеют следующий вид:

~ _ ^т ^ш— 1 ^т— 1

а.21 — —р---,

7 2(мт-цт_1)

02/ — -1--^то-1,

~ _ Ущ _ ^т ^т — 1

«2г — к ,

С2/ — С2г — ^то.

Для дифференциального следствия полиномы можно получить дифференцированием (2.13) и (2.14):

¥2Ж {х) — 2а,21 х + Ъ21, (2.15)

(Ж) — 2й2г Ж + &2г. (2.16)

Схемы, построенные на данных полиномах, имеют второй порядок точности, однако обладают заметной дисперсией. В дальнейшем будем обозначать схемы, использующих полиномы (2.13) и (2.14), как С1Р2Ь и С1Р2К соответственно.

Наряду с этими полиномами, рассмотрим также полином первого порядка точности на рассматриваемом шаблоне

^1(ж) — а1х + Ь1. (2.17)

Исходя из условий на границах интервала

^1(0) — Что, — Что-1

находим его коэффициенты:

а 1 —-к—,

&1 — Что.

Аппроксимация с помощью данного полинома представляет собой схему первого порядка типа КИР (Куранта-Изаксона-Риса). Дифференциальное следствие (2.8), в этом случае, будет аппроксимироваться по формуле:

^ (ж) — а1.

Обозначим данную схему как СШ.

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

2.4.3 Исследование поведения интерполяционных полиномов

Аппроксимируем первую производную функции и(£, К) на интервале (-К, 0), используя формулу первого порядка точности:

V* = = (ито - ит-\)/к.

Далее рассмотрим поведение различных полиномов, в зависимости от значений искомых величин на выбранном шаблоне.

Знаки V*, уто и Уто-\ совпадают, что означает выполнение условий:

УтоУто-1 ^ ° /о 1 сЛ

V*Vто ^ 0.

С целью повышения порядка точности, целесообразнее использовать кубический полином Р3(х), однако при выполнении условия (2.18), у функции Рз(х) возможны экстремумы на интервале (-К, 0). Условие наличия экстремума равносильно выполнению условия изменения знака производной на рассматриваемом интервале. Данное условие можно записать, например, в виде (значение в вершине параболы по знаку не совпадает со значением на границах интервала, и вершина параболы лежит внутри интервала):

—К ^ х0 ^ 0, ,

0 ' (2.19)

Р3(хо)у* < 0,

где х0 = — - координата вершины параболы.

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

тгп(ут, Vто—1) < V* ^ тах(у

то, Уто—1 ), (2.20)

то полиномы (2.13) и (2.14) лежат по разные стороны от прямой (2.17) на всем интервале (—К, 0). Выпуклая комбинация данных полиномов также будет полиномом второго порядка:

Р21г(х) = аР2(х) + (1 - а)^2Г(х), (2.21)

для всех 0 ^ а ^ 1. Поскольку полиномы лежат по разные стороны от прямой (2.17), то можно выбрать такое а, что в интервале (—Н, 0) будет выполнено Р21г(ж) — ^1(ж), и используя монотонный полином первой степени, можно аппроксимировать решение вторым порядком. Это получается потому, что коэффициент перед х2 зануляется в полиноме (2.21).

Если условие (2.20) не выполнено, то полиномы (2.13) и (2.14) лежат по одну сторону от прямой (2.17) на интервале (—Н, 0). В данной ситуации можно выбрать тот из полиномов, который не имеет экстремума на интервале (—Н, 0); в случае, если оба полинома имеют экстремум, можно воспользоваться схемой первого порядка (2.17). Однако, как показали тестовые расчеты, по крайней мере для линейных уравнений, можно выбрать тот из полиномов, который построен на производной, лежащей ближе к у*; данное условие можно записать как

(ж), Дто_ 1 ^ Дто 2( ^ то 1 " то (2.22)

^2Г(ж), Д то-1 > Д то,

где Д { — |Уi — V*!. Использование такого подхода позволяет обеспечить второй порядок точности, однако не гарантирует отсутствие экстремума на отрезке.

Рассмотрим следующий случай. Знаки уто и уто- 1 различны, т. е.:

УтоУто— 1 < (2.23)

В данном случае, при использовании интерполяции выше первого порядка, на рассматриваемом интервале всегда будет экстремум. Однако, поскольку производные меняют знак, то наличие экстремума не является нефизичной осцилляцией и, как показывают тестовые расчеты, применение полиномов порядков выше первого, не вызывают появления нефизичных осцилляций в решении. Полиномы второго порядка при таких условиях лежат по одну сторону от прямой, заданной полиномом (2.17) на отрезке (—Н, 0); по этой причине невозможно воспользоваться их выпуклой комбинацией. Однако, можно использовать полином третьей степени (2.10), либо выбрать полином второй степени, исходя из условия (2.22). Для получения монотонного решения, необходимо использовать интерполяцию первого порядка (2.17).

Наконец, пусть знаки уто, уто_ 1 одинаковые, но отличны от знака у*, т. е.:

утоуто-1 ^ 0, (2 24)

У*Уто < 0.

Тогда поскольку полином третьего порядка (2.10) имеет два экстремума в интервале (-К, 0), то он не может быть использован. Полиномы второго порядка (2.13) и (2.14) также имеют экстремумы на рассматриваемом интервале, однако лежат по разные стороны от прямой (2.17). Используя их выпуклую комбинацию (2.21), всегда можно занулить коэффициент перед квадратичным членом в полученном полиноме и использовать монотонный линейный полином (2.17), дающий в данном случае второй порядок аппроксимации.

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

2.4.4 Гибридные схемы

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

2.4.5 Сеточно-характеристическая монотонизация

Сеточно-характеристические схемы опираются на использовании характеристического критерия монотонности [19; 53]:

тгп(ит,ит-1) ^ ^ тах(ито,ито-1). (2.25)

Подробно построение разностных схем, в том числе для продолженных систем уравнений, основанных на использовании характеристического критерия монотонности (2.25), описано в [53]. Там же приводятся конечно-разностные схемы третьего и четвертого порядка точности для продолженных систем уравнений (2.1) и (2.8) на компактном шаблоне (2.9). В данной работе, в отличии от [53], используется интерполяционный полином, для решения системы (2.1), и его производная, для решения дифференциального следствия (2.8). В работе [53]

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

Гибридизация, основанная на сеточно-характеристическом критерии (2.25), для схем описанных в данной работе, осуществляется следующим образом. Вычисляется значения и^1, с использованием некоего интерполяционного полинома; для полученного этого значения проверяется выполнение данного критерия монотонности. В случае если критерий не выполнен, порядок полинома понижается до его выполнения. Данная операция производится последовательно над полиномами ^3(ж), (ж), ^2г(ж) и ^1(ж). После выбора полинома, для интерполяции вычисляется значение уто+1, используя его производную. Для данной схемы будем в дальнейшем использовать обозначение БИС1 (бикомпактная интерполяционная схема).

2.4.6 Гибридная схема третьего порядка точности

В работе также использовалась гибридизация схемы третьего порядка точности, исходя из поведения полиномов на рассматриваемом интервале. Схема строится следующим образом. Если выполнены условия (2.18) и (2.19), используется схема второго порядка исходя из (2.22). Если выполнено (2.18), но не выполнено (2.19), используется схема третьего порядка (2.10). Если выполнено условие (2.20), используется полином первого порядка (2.17), но, как показано выше, он будет давать второй порядок. В случае выполнения условия (2.23), используются полиномы второго порядка, исходя из (2.22). Наконец, при выполнении условия (2.24) используется полином первого порядка (2.17), также дающий второй порядок. Обозначим данную схему как БИС2. Схема тестировалась на начальных данных в виде импульса сложной формы (2.4), при различных значениях числа Куранта, во всех случаях схема показывала монотонное поведение.

2.5 Тестирование схем

Тест на качественное поведение разностных схем проводился для начальных данных в виде импульса сложной формы (2.4). Использовалась сетка, состоящая из 200 узлов и с периодичными граничными условиями. Решение представлено в момент времени £ = 4.0, шаг по времени выбирался исходя из числа Куранта 0.4. Скорость переноса во всех тестах Л = 1. Результаты для всех схем представлены на рис. 2.5.

Ь)

• 4 > 1

-

1 1 1 1 1

111111111

: : : : А

• •

>

.....: 1 -1 -:..... ...... 1.....-:1--1:.....I •• -I.....•

: т Т: • * 1 1 4 4 ■ * ? т Ч • :

1 1 1 1 1 1 1 1 1

Рисунок 2.5 — Численное (точки) и точное (линия) решение линейного уравнения переноса (2.1) для начальных данных (2.4) через 1000 шагов разностной схемы. а) - СШ; б) - С1Р; в) - С1Р2Ь; г) - С1Р2Я; д) - БИС1; е) - БИС2.

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

Таблица 1 — Нормы ошибок для различных схем

^2 Аэо

СШ 0.562 0.092 0.810

С1Р 0.055 0.019 0.389

С1Р2Ь 0.262 0.050 0.635

С1Р2Ь 0.271 0.051 0.644

БИС1 0.068 0.024 0.417

БИС2 0.054 0.020 0.429

имеет наименьшие осцилляции и, в ряде случаев, может быть применена без монотонизатора. Гибридные схемы имеют монотонное поведение, схема БИС1 монотонна по сеточно-характеристическому критерию, схема БИС2 также показывает монотонное поведение; к тому же, на коротких импульсах, показывает более лучшее поведение, чем схема БИС1. При других значениях числа Куранта схемы показывают похожий результат, осцилляции для схем БИС1 и БИС2 также отсутствуют.

Нормы ошибок, для данного теста, приведены в таблице 1. Использовались следующие нормы: Ь1 — ^ Ь2 — )1/2, ^^ — тах^^,

г г

Х{ — щ — ^еог, и{ - решение полученное численно, - точное решение,

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

Как видно из расчетов, наименьшие нормы дает схема С1Р, однако из не осциллирующих схем, меньшие ошибки имеют расчеты, полученные по схеме БИС2.

2.5.1 Импульс формы прямоугольного треугольника

Дополнительно проводился тест на импульсе в форме прямоугольного треугольника (2.26).

,(х + 0.4)/0.8, —0.4 < 0.4, и(0,х) — {к " (2.26)

0, в противном случае.

Остальные параметры аналогичны предыдущему тесту. Результаты для схем СШ, С1Р, БИС1 и БИС2 представлены на рис. 2.6.

a) b) c) d)

- - 1-----